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 runningcanonical_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.