TDEP workflow using VASP

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

Here we show an example with SrTiO3 (also as a test for the new workflow).

Stochastic sampling approach

Compute forces and positions using VASP

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 ensemble.
> canonical_configuration -n 24  -t 300  --quantum --debye_temperature 400

Note that debye_temperature is used to generate a fake force constant 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 forces on each atom.

In the case of VASP, we need to prepare INCAR, POTCAR, and KPOINTS. Submit jobs using a script like the following.

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

Extract force constants using TDEP

After the scf calculations are done, gather all the CONTAR_confxxxx files. To generate all the necessary input files for extract_forceconstants, include infile.forces, infile.positions, infile.meta, infile.stat, see tutorials on the TDEP GitHub page.

For polar materials, we need to prepare infile.lotosplitting for the long-range electrostatic corrections. 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. The way to obtain them from VASP outputs is also detailed in the TDEP parsing tutorials.

Extract the force constants (for non-polar materials, neglect --polar):

> extract_forceconstants -rc2 5 --polar

The resulting outfile.forceconstant contains all the information for the calculations of lattice dynamical properties, such as the phonon dispersion.

For polar materials, the short-ranged force constants extracted by fitting the displacement and subtracted forces and the quantities to compute the long-ranged (dipole-dipole) force constants, such as the polar correction type, dielectric tensor, and the Born effective charges. In other words, we subtract the long-range forces (due to dipole-dipole interactions) from the DFT forces, and only extract the effective short-ranged force constants using the least-squre fitting.

For example, in qe2pert-src/lattice_tdep.f90, adapted from TDEP src/libolle/lo_longrange_electrostatics_dynmat.f90, it contructs the dynamical matrix at a given q-points (for pc = 3), first obtaining the short-ranged dynamical matrix by Fourier transformation of the short-ranged force constants, and then adding back the long-ranged (dipole-dipole) contribution.

Test phonon dispersion relations

> 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 dispersions. (The path is generated automatically based on the unit cell information).

Iterate the sampling process

To refine and converge the force constants, we need to go throught a self-consistent sampling 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 file, compute the new phonon dispersion and check if it converged (by comparing with phonon dispersion in the previous step).

For more information on self-consistent sampling, please refer to TDEP tutorial page.

Ab-initio molecular dynamics approach

References:

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 requires 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 uses force constants 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 due to QE’s performace on large supercell (slow and easy to crash). But one can still use QE to compute forces and positions in the TDEP self-consistent iteration.

  • Since TDEP canonical_configuration does not support QE input formats, one can use -of 1 in running canonical_configuration to generate VASP CONTCAR files, which contains the supercell configuration. Then, one can use the script provided in the example, contcar_to_qe.py to convert the contcar_conf0001 file to qe_conf0001.

  • 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. The merging script can be found in the example provided in folder tdep/tools. Make sure to use the same unit cell structure for scf calculations in scf_ref.in, infile.ucposcar, and the scf.in in DFT and DFPT calculations from Quantum Espresso.

  • After running scf calculation using QE, one gets the output file scf.out with the force information. To extract this info, run collect.sh in the folder containing scf.out, which can concatenate forces and positions for multiple output files. It will generate infile.positions and infile.forces, needed in TDEP.

  • The parsing scripts and tutorials for preparing the input files for extract_forceconstant from QE ph.x outputs can be found here.

Technical issues

Referencing atomic positions

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

In the paper PRL 125, 045701(2020), titled anharmonicity and ultralow thermal conductivity in lead-free halide double perovskites, they 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 positions and displacements.

Born effective charge and Ewald parameter for polar materials

In both code, only the short-ranged force constants 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 force constants and dynamical matrices are computed on the fly.

Born effective charge

TDEP enforces the symmetry of the Born effective charges. At the bottom of the infile.forceconstant file for polar materials, there is a block of Born effective charges in the irreducible representation,

   3  # number of irreducible components in the Born charges
   2.554432935340        # irrep of Born charges, # of lines should agree with the number above
   7.282998859270
   -5.767556255820
     1.000000000000     0.000000000000     0.000000000000 # coefficients of Born charges, # of lines = 9 * # of atoms
     0.000000000000     0.000000000000     0.000000000000
     0.000000000000     0.000000000000     0.000000000000
     ...

followed by a second block of Born effective charges of each atom:

     2.554432935340     0.000000000000     0.000000000000 Born effective charge atom   1
    -0.000000000000     2.554432935340    -0.000000000000
     0.000000000000     0.000000000000     2.554432935340
     7.282998859270    -0.000000000000     0.000000000000 Born effective charge atom   2
    -0.000000000000     7.282998859270     0.000000000000
     0.000000000000    -0.000000000000     7.282998859270
     ...

These two representations are related. The coefficient matrix, coeff_Z, has the dimensions (9 * # of atoms, # of irrep of Born charges). Multiplying coeff_Z with the vector of irrep Born charges x_Z results in a vector, Z_flat, of 9 * # of atoms. Folding each set of 9 lines for each atom into a 3 * 3 matrix, we arrive at the Born effective charges of each atom in the unit cell.

One can set the flag, symmborn, in the input file of qe2pert.x to decide which block of information it uses.

In qe2pert.in, the default is symmborn = .true.. In this case, qe2pert reads in coeff_Z and x_Z, symmetrizes the inputs, and uses the resulting Born effective charges in subsequent calculations. If the user sets symmborn = .false., the second block of Born effective charges is read and used by qe2pert.x directly, without Hermiticity enforced.

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 (lo_longrange_electrostatics_dynmat.f90)

  inv4lambda2 = 1.0_r8/(4.0_r8*(ew%lambda**2))
 do ig = 1, ew%n_Gvector
     Gvec = ew%Gvec(:, ig)
     Kvec = (Gvec + q)*lo_twopi
     Kvec = lo_chop((Gvec + q)*lo_twopi, lo_sqtol)
     if (lo_sqnorm(Kvec) .lt. lo_sqtol) cycle
     ! K-norm with dielectric metric
     knorm = dot_product(Kvec, matmul(eps, Kvec))
     invknorm = 1.0_r8/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.