TDEP workflow using VASP

We are developing a new and more user-friendly PERTURBO+TDEP workflow.

Note: since QE scales badly w.r.t to system size, e.g., QE is bad at supercell calculations, and TDEP really needs scf calculations a very large supercell. I will used VASP for force calculation on supercells. The combination of QE and VASP might cause some issues. But since we only needs the force-constant (gauge-invariant) from TDEP. As long as we make sure consistent parameter are used in both QE and VASP (exchange correction function, and lattice constant), and VASP generates forceconstants that are accurate, that should be fine.

I am using SrTiO3 as an example. (also as a test for the new workflow).

Stochastic sampling approach

We start with the crystal structure in cif format.

#generate unit cell file
> cif2cell srtio3-rt.cif -p vasp --vasp-format=5 -o infile.ucposcar

#generate supercell file
> generate_structure -d 4 4 4
> mv outfile.ssposcar  infile.ssposcar

#sample supercells emulating a canonical ensembels.
> canonical_configuration -n 24  -t 300  --quantum --debye_temperature 400

Note that debye_temperature is used to generate a fake forceconstant in the initial step, which is used to generate displacements. Starting from the second iteration, we use the extracted force constants to generate displacements, thus debye_temperature is no longer needed.

Once the displaced supercells are generated, we run scf calculations to compute the force on each atoms.

In the case of VASP, we need to prepare INCAR, POTCAR, and KPOINTS. Submit jobs using scripts.

bash
WDIR="rundir"
mkdir -p ${WDIR}

for i in {1..12}
do
   if [ "$i" -lt 10 ]; then
      ND="000${i}"
   else
      ND="00${i}"
   fi

   echo "running conf${ND} ..."

   #prepare input file
   rm -f ${WDIR}/*
   cp INCAR   ${WDIR}
   cp POTCAR  ${WDIR}
   cp KPOINTS ${WDIR}
   cp contcar_conf${ND} ${WDIR}/POSCAR

   #run calculations
   cd ${WDIR}
   mpirun -np $Nproc $VASP
   cd ..

   #save results
   cp ${WDIR}/OUTCAR  OUTCAR_conf${ND}
done

After the scf calculations are done, gather all the CONTAR_confxxxx file. Use process_contar_5.3.py to collect all the displacements and forces.

> process_outcar_5.3.py OUTCAR_conf*

This scripts will generate all the necessary input files for extract_forceconstant, include infile.forces, infile.positions, infile.meta, infile.stat.

For polar materials, we need to prepare infile.lotosplitting for the long-range electrostatic corrections. See an example for Bi2Te3. The first three lines is the dielectric tensor, after that three lines per atom in infile.ucposcar with the Born effective charges. The dielectric tensor is unitless, and the Born effective charges are in electron charges.

Extract the forceconstant (for non-polar materials, neglect --polar and -pc 3):

> extract_forceconstants -rc2 5 --polar -pc 3

here, --polar means add dipole-dipole corrections for polar materials, and -pc 3 indicates using the polar correction scheme introduced in my 2018 PRL paper. Namely, we subtract the long-range force (due to dipole-dipole interactions) from the DFT forces, and only extract the effective short-ranged force constant using the least-squre fitting.

The resulting outfile.forceconstant contains the short-ranged forceconstants extracted by fitting the displacement and subtracted forces and the quantities to compute the long-ranged (dipole-dipole) forceconstant, such as the polar correction type, dielectric tensor, and the Born effective charges. For example,

           3 # This is a forceconstant for a polar material.
       6.255018086086       0.000000000000       0.000000000000 Dielectric tensor xx xy xz
       0.000000000000       6.255018086086       0.000000000000 Dielectric tensor yx yy yz
       0.000000000000       0.000000000000       6.255018086086 Dielectric tensor zx zy zz
  0.99826965332031237       # Coupling parameter in Ewald summation
           4  # number of irreducible components in the Born charges

The number 3 in the first line indicates the polar correction type. The fifth line shows the Ewald summation parameters. (The born effective charge for each atoms are followed, but not shown here).

The outfile.forceconstant file includes all the informations for the calculations of lattice dynamical properties, such as the phonon dispersion.

> mv outfile.forceconstant infile.forceconstant
> phonon_dispersion_relations

phonon_dispersion_relations reads information from infile.forceconstant, and figures out the correct way to compute phonon dispersion. (the path is generated automatically based on the unit cell information). For example, in libolle/type_forceconstant_secondorder_dynamicalmatrix.f90, it contruct the dynamical matrix at a given q-points using different approach according the polar correction type, for pc = 3, it first obtained the short-ranged dynamical matrix by Fourier transformation of the short-ranged forceconstant, and then add back the long-ranged (dipole-dipole) contribution.

 case(3)
    ! Incredibly smart correction. First the shortranged-part:
    do a1=1,fc%na
    do i=1,fc%atom(a1)%n
        a2=fc%atom(a1)%pair(i)%i2
        Rv=fc%atom(a1)%pair(i)%lv2
        qdotr=dot_product( Rv,qv )
        exp_iqr=cmplx( cos(qdotr), sin(qdotr) ,flyt)
        ! Normal dynamical matrix
        Dt(:,:,a1,a2)=Dt(:,:,a1,a2)+exp_iqr*fc%atom(a1)%pair(i)%m
    enddo
    enddo
    ! Then the polar part
    call longrange_dynamical_matrix(fc,PDt,p,qpoint%w,reconly=.true.)

    ! And the final non-analytical part
    if ( nonanalytical ) then
        call dynmat_at_gamma(fc,p,qdir,GDt)
    endif

    Dt=Dt+PDt+GDt
end select

To refine and converge the forceconstants, we need to go throught an iterative process. We generate a new set of supercell configurations from the newly obtained infile.forceconstant using canonical_configuration, repeat the above process to obtained a new forceconstant files, compute the new phonon dispersion and check if it converged (by comparing with phonon dispersion in the previous step).

Ab initio molecular dynamic approach

Reference:

How to run molecular dynamic using VASP, using Cs2AgBiBr6 as example, at T=300K:

First, prepare infile.ucposcar file, same approach as above. To run molecular dynamics, you need a larger simulation cell.

> generate_structure -na 320
> mv outfile.ssposcar infile.ssposcar

which generates a supercell with 320 atoms.

Starting molecular dynamics from ideal lattice positions usually require a long equilibration time that has to be discarded. This time can be minimized by using a better seed, which can be generated using the tool canonical_configuration.

canonical_configuration: Use forceconstants or a Debye temperature to generate uncorrelated supercell configurations emulating a canonical ensemble. These configurations can be used to either start ensemble runs of molecular dynamics with negligible equilibration time, or be used to directly sample phase space.

> canonical_configuration -n 1 -t 300 --debye_temperature 400

Which will generate a configuration contcar_conf0001, and we will start from this structure to run molecular dynamic simulation using VASP.

TDEP workflow with Quantum Espresso

TDEP does not support QE officially, partially before QE’s performace on large supercell is terrible (slow and easy to crash). But, in my case, sometimes to keep TDEP+PERTURBO calculations consistent, I need to use QE for the force calculations as well. For this reason, I added support for QE in TDEP.

Adjustment comparing to the VASP workflow:

  • generate displacement configurations, -of 8 specifies the output format of Quantum Espresso.
> canonical_configuration -n 24  -t 300  --quantum --debye_temperature 400 -of 8
  • To perform DFT calculation on the generated supercell, we need to generate a proper input file for QE. To do that, we put all the computational setting in scf_ref.in, and merge the generate structure data and the setting parameters.
> process_pwscf -m input --prefix scf qe_conf0001 scf_ref.in

The first argument qe_conf0001 is the supercell structure generated by canonical_configuration, the second file contains the computational setting. It will generate a new input file scf.in which has the structure data and setting parameters.

  • After running scf calculation using QE, one gets the output file scf.out with the force information. To extract this info, run
> process_pwscf -m output --prefix infile scf.in. scf.out

It will generate infile.positions and infile.forces, needed in TDEP. (To append the position and forces to existing files, add the option: -a)

The issue with Born effective charge

TDEP enforce the symmetry of the Born effective charge. In the case of SrTiO3, it enforce the Born effective charge to a wrong value. To bypass this issue, I modified the TDEP code to skip the symmetrization of Born effective charge.

When I run extract_forceconstant, I need to add a new tag --nosymmborn

> extract_forceconstants -rc2 5 --polar -pc 3 --nosymmborn

The generated forceconstant file will become

           3 T # This is a forceconstant for a polar material.
           ...

The T after 3 means that the forceconstant was generated using un-symmetrized Born effective change, which is written in this forconstant file. the fc%readfromfile procedure will read this tag, if T then we will initialize the fc%atom(a1)%loto%Z with the Born effective change stored in this file, if F, then it will generate the Born effective change from the irreducible coefficent (as TDEP normally does). Since phonon_dispersion, canonical_configuration, PERTURBO interface are all initialize Z via calling fc%readfromfile, so they will all work properly.

The issue with the reference atomic positions in TDEP

When extracting the displacements, we need a reference equilibrium atomic position. Usually we a static structural relaxed position as a reference atomic positions, but this sometimes does not work very well, especially for low-symmetry structure, such as tetragonal phase of SrTiO3.

In the paper PRL 125, 045701(2020), titled anharmonicity and ultralow thermal conductivity in lead-free halide double perovskites, the used a different approach, as discribed in the supplemental materaials.

The fitting is constrained to impose the full symmetry of the underlying double perovskite structure on the IFCs. The reference positions of the Cs, Ag, and Bi ions are thus completely fixed by symmetry, while a single degree of freedom (DOF) remains for the reference positions of the Br ions. This DOF specifies the position of the Br ion along the Ag-Br-Bi bond. We use the average value of this DOF extracted from AIMD at the temperature of interest as reference positions in Eq. (1).

Conventionally, this DOF could be determined by a static structural relaxation, potentially at a thermally expanded volume. For Cs2AgBiBr6 we find, however, that using such statically relaxed positions of the Br ions produces a large, artificial, softening of a high frequency mode around the Γ-point. This is a further indicator of the strong anharmonicty of Cs2AgBiBr6.

We have to use the AIMD approach, instead of the random sampling, to generate the equilibrium position and displacement.

Interface to PERTURBO

Due to some unknown bugs introduced into tdep-devel via Commit on Sep 29, 2018 (commit #6628fcf), TDEP compiled using gfortranon tempo crashes (Segmentation fault).

Besides, some bugs related to intel compiler are introduced via commits on Jun 21, 22, 2018 (commit #91d9363, #1f1cdb5, #e91af27). TDEP compiled using intel compiler on tempo crashes ( SOLVING FOR FORCECONSTANTS, segmentation fault occurred).

Thus, PERTURBO is interfaced to the version commited on Jun 14, 2018, (commit #ef3d150). Will be upgraded to the last commit once the bugs are fixed in TDEP.

Treatment of long-range dipole-dipole correction in TDEP and PERTURBO

In both code, only the short-ranged forceconstants are stored in IFC. Since the long-ranged dipole-dipole interactions are determined by dielectric tensor, Born effective charges, and Ewald parameter, only these quantities are stored, and the long-ranged forceconstants and dynamical matrices are computed on the fly.

In the PERTURBO-TDEP interface, qe2pert.x will read data from infile.forceconstant, map it to the Wigner Seitz Supercell (WS) generated by PERTURBO, and then write them to prefix_ephwan.h5 file. Note that the WS should be larger than the rc2 in TDEP, if anyone of the pairs in infile.forceconstant cannot be mapped to WS, qe2pert.x will throw an error.

qe2pert.x will check the polar correction type. if PC is 0, then polar correction is turned off. If correction type is 3, polar correction is turned on, otherwise, throw an error. if correction type is 3, write dielectric tensor, Born effective charnges, and Ewald parameter to prefix_ephwan.h5.

No change is needed for perturbo.x, since perturbo.x will go throught the same process no matter it is TDEP or not. But one needs to make sure perturbo.x is doing the equivalent operations as TDEP.

PERTURBO:

polar_correction.f90
subroutine dyn_mat_longrange_3d

TDEP:

type_forceconstant_secondorder_dynamicalmatrix.f90
subroutine dynamicalmatrix
! which calls
call longrange_dynamical_matrix(fc,PDt,p,qpoint%w,reconly=.true.)

subroutine longrange_dynamical_matrix is defined in type_forceconstant_secondorder_loto_gonze.f90 (moved to type_forceconstant_secondorder_loto.f90 in later version).

Design idea

  • Follow the same stragy as PERTURBO-fortran.
  • only infile.ucposcar and infile.forceconstant are required.
  • only qe2pert.x will work with TDEP. and it only work with pc = 3. perturbo.x is not affected.
  • add a new parameter tdep
  • qe2pert.x will read data from infile.forceconstant, and replace the lattice dynamical info from phonon
  • The Ewald parameter for phonon dispersion is also from infile.forceconstant
  • add __TDEP precompile option, if __TDEP, then libolle.a is required.
  • map TDEP ifc to PERTURBO/wsz, so that WSZ in PERTURBO should be larger enought to contain all the TDEP-ifcs, if not, throw an error.
  • in TDEP mode, use dielectric tensor and Born effective charge from TDEP, and use it for both ph and eph.
  • use different Ewald parameters for ph and eph, in TDEP mode, ph-ewald is set to TDEP-ewarld, otherwise use 1.0, so that it is consistent with that used in rgd-blk. polar_alpha and loto_alpha
  • write tdep flag in epr file, just a signal that this file is generated using TDEP, the tdep flag will not be used in perturbo.x.
  • for now, TDEP does not work with 2D polar correction. add a warning if when using TDEP for 2D polar materials.

Check List

  • map TDEP IFC to PERTURBO WS
  • read polar parameter from TDEP, and save it to PERTURBO/ph
  • polar correction type
  • number of atoms
  • lattice vectors
  • atomic positions
  • check if alpha in PERTURBO has the same meaning as lambda in TDEP
  • check the order of (i,j) in fc_m (TDEP) and in ifc (PERTURBO)
  • The polar part in PERTURBO breaks symmetry (degeneracy is lifted by a tiny bit), while that of TDEP does not, matdyn also has good symmetry. Fixed.
  • Make sure longrange_dynamical_matrix in TDEP is equivalent to dyn_mat_longrange_3d in perturbo-src.
  • Check PERTURBO and tdep/phonon-dispersion generate the exact same dynamical matrix (the two approaches should be equivalent), both the longe range and short range.

Testing

Test on STO, the e-ph matrix elements are consistent with previous results published on PRL.

Ewald parameter

PERTURBO (polar_correction.f90)

 !alpha is the Ewald summation parameter (in unit of (tpiba)**2 )
 real(dp) :: alpha

 falph = 4.0E0_dp * alpha

 !xq is crystal coordinate of q-point. xqr is (q+G)
 xqr(1:3) = xq(1:3) + (/m1, m2, m3/)
 ! transform to cartesian coordinate (in unit of tpiba)
 call cryst_to_cart(1, xqr(1:3), bg, 1)
 ! (q+G)*epsilon*(q+G)
 qeq = dot_product(xqr(1:3), matmul(epsil(1:3,1:3), xqr(1:3)))
 qfac = exp( -qeq/falph ) / qeq

TDEP (type_forceconstant_secondorder_loto.f90)

 write(u,*) fc%loto%lambda,' # Coupling parameter in Ewald summation'

 inv4lambda2=1.0_flyt/(4.0_flyt*(fc%loto%lambda**2))

 Gvec=fc%loto%qvec(:,qi)
 Kvec=(Gvec+q)*lo_twopi
 ! K-norm with dielectric metric
 knorm=dot_product(Kvec,matmul(fc%loto%eps,Kvec))
 invknorm=1.0_flyt/knorm
 expLambdaKnorm=exp(-knorm*inv4lambda2)
 ! Half of the Chi-function
 partialChi=expLambdaKnorm*invknorm

So alpha * (tpiba**2) in PERTURBO is equal to lambda**2 in TDEP.