## 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 AIMDat 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.