Interface to TDEP and SSCHA | Perturbo

For materials with strong lattice anharmonicity, it is critical to incorporate anharmonic phonons in electron-phonon calculations, take SrTiO₃ (see this paper) and La₂CuO₄ (see this paper) for example. The temperature-dependent effective potential (TDEP) and self-consistent harmonic approximation (and its variants, for example SSCHA) are two frequently used approaches for simulating anharmonic lattice dynamics. PERTURBO provides interface to TDEP package and SSCHA package.

TDEP

Directory: example07-sto-tdep-sscha link

Currently, PERTURBO only reads 2nd-order force constants computed from TDEP.

Compilation

PERTURBO is now compiled independent of TDEP. Namely, there are no additional steps for the PERTURBO compilation.

Usage

To incorporate TDEP interatomic 2nd-order force constants in PERTURBO, one needs to prepare two key files from TDEP, namely infile.ucposcar and infile.forceconstant, and set tdep = .true. in qe2pert.in (the main input file for qe2pert.x).

Note: The required TDEP files must be in the following formats for proper reading and processing by PERTURBO. One can rewrite the outputs generated from earlier (or later) versions of TDEP to enforce compatibility. The TDEP version we tested and recommend is commit a37432a, committed on Dec 5, 2024.

Here is an example of infile.ucposcar, which should be the same as the input used in TDEP to obtain the force constants:

SrTiO3
90003602 
000000000000000   0.000000000000000   0.000000000000000 
000000000000000   1.000000000000000   0.000000000000000 
000000000000000   0.000000000000000   1.000000000000000
   Sr  Ti  O                                                
 1   3                                                
Direct                                                      
000000000000000   0.000000000000000   0.000000000000000 
500000000000000   0.500000000000000   0.500000000000000 
000000000000000   0.500000000000000   0.500000000000000
500000000000000   0.000000000000000   0.500000000000000 
500000000000000   0.500000000000000   0.000000000000000 

The first line is a comment line. The next line is the global scaling factor of lattice vectors (and of atomic positions if they are in Cartesian coordinates). The following three lines are lattice vectors. The next two lines indicate the names of atomic species and numbers of these atoms in the unit cell.

The next line starts with either Direct or Cartesian, indicating either the following atomic positions are in direct (fractional) or Cartesian coordinates(qe2pert.x reads the first alphabet c/d/C/D).

Here is an example of the beginning lines of infile.forceconstant, which should be the same as TDEP output, outfile.forceconstant:

          5               How many atoms per unit cell
    7.021350089479943     Realspace cutoff (A)
        119               How many neighbours does atom   1 have
          1               In the unit cell, what is the index of neighbour 1 of atom 1
  0.000000000000000E+000  0.000000000000000E+000  0.000000000000000E+000
   2.44565967557536       0.000000000000000E+000  0.000000000000000E+000
  0.000000000000000E+000   2.44565967557536       0.000000000000000E+000
  0.000000000000000E+000  0.000000000000000E+000   2.44565967557536     
          5               In the unit cell, what is the index of neighbour 2 of atom 1
  0.000000000000000E+000  -1.00000000000000       0.000000000000000E+000
 -0.278127281008224       0.323848538919098       0.000000000000000E+000
  0.323848538919098      -0.278127281008224       0.000000000000000E+000
  0.000000000000000E+000  0.000000000000000E+000  0.492528494328962     
  ...

Here is an example of the input file for qe2pert.x with the TDEP force constants:

&qe2pert
 prefix = 'sto-nosoc'
 outdir = './tmp'
 phdir = '../pw-ph-wann/phonon/save'
 nk1 = 4
 nk2 = 4
 nk3 = 4
 num_wann = 3
 dft_band_min = 21
 dft_band_max = 23
 lwannier = .true.
 tdep = .true.  ! flag to turn on tdep
/

In the PERTURBO-TDEP interface, qe2pert.x reads data from infile.ucposcar and infile.forceconstant, maps the TDEP force contants to the Wigner Seitz Supercell (WS) generated by PERTURBO, and then writes them to ‘prefix’_epr.h5. If any pair of force constants in infile.forceconstant cannot be mapped to WS, qe2pert.x will issue a runtime error. To solve this issue, one needs to either increase the coarse \(\mathbf{q}\)-grid or decrease the rc2 parameter (the cutoff for the 2nd force constants) in TDEP, since WS needs to be larger than rc2.

Once the TDEP force constants are employed in ‘prefix’_epr.h5, perturbo.x will use them automatically, there is no need to modify the input file for perturbo.x. If tdep = .false., qe2pert.x will generate ‘prefix’_epr.h5 with 2nd-order force constants from PHonon QE calculation.

Technical details

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

In both TDEP and PERTURBO, only the short-ranged force constants are stored in the form of interatomic force constants (IFC) for polar materials, while the long-ranged dipole-dipole contributions to the dynamical matrices are computed on the fly, using the dielectric tensor, Born effective charges, and Ewald parameters. For more details of the dipole-dipole corrections for polar materials, we refer the users to this paper and its Supplemental Material.

Here is an example to run TDEP:

extract_forceconstants -rc2 5 --polar

Here, --polar means the dipole-dipole correction is on for polar materials. If not specifying --polar, the material is assumed to be nonpolar.

The infile.forceconstant file from TDEP contains the short-ranged force constants extracted by fitting the displacement and subtracted forces and the quantities to compute the long-ranged (dipole-dipole) force constant, such as the polar correction type, the dielectric tensor, the Born effective charges, and so on.

This is an example of the end of the infile.forceconstant file:

           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
   1.00143890380859       # Coupling parameter in Ewald summation
           3  # number of irreducible components in the Born charges
   2.554432935340        # irrep of Born charges
   7.282998859270
   -5.767556255820
     1.000000000000     0.000000000000     0.000000000000 
     0.000000000000     0.000000000000     0.000000000000
     0.000000000000     0.000000000000     0.000000000000
     ...
     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
     ...

The number 3 in the first line indicates the polar correction type. The fifth line shows the Ewald summation parameter. The next block shows the information on Born effective charges. The first block contains coefficients of Born charges, with number of lines equal to 9 * number of atoms. qe2pert.x checks the polar correction type. If it is 0, then the polar correction is turned off. If the correction type is 3, polar correction is turned on, otherwise, qe2pert.x will issue a runtime error. If the correction type is 3, the dielectric tensor, the Born effective charges, and the Ewald parameter will be written to the ‘prefix’_epr.h5 file.

For more details, including running TDEP with VASP and Quantum Espresso, the issue with Born effective charge, etc., please refer to this document.

SSCHA

Directory: example07-sto-tdep-sscha link

Installation

Please refer to SSCHA website for details about download and installation. To avoid introduing more complexity, we reformulate the output of SSCHA in the format of QE output files, which is made possible by a patch to SSCHA source code.

Copy the modified file Phonons.py provided in the example folder into SSCHA,

cp example07/sscha/patch/Phonons.py SSCHA_PATH/CellConstructor/cellconstructor

Or you can find the patch code block in the provided Phonons.py file labelled by the keyword PERTURBO and then directly copy/paste into Phonons.py in “SSCHA/CellConstructor/cellconstructor”. (in /site-packages/ if installed with pip)

Usage

To incorporate SSCHA interatomic 2nd-order force constants in PERTURBO, one could first run DFPT calculations using Quantum Espresso, and then replace prefix_dyn.xml files with those generated from SSCHA using the same q-grid. No additional input parameters for qe2pert.x is needed. The example folder includes example calculations using SSCHA.

Next, we will break down the entire simulation into five steps. For each step, there are corresponding scripts provided in the example folder.

Step 1: Sample NVT ensembles starting from DFPT calculations

After you performed the DFPT calculation with the output file in the conventional format (not xml), we are able to make the NVT ensemble.

We start with reading the DFPT results,

 # Specify the number of irreducible q points
 dyn = CC.Phonons.Phonons("./phonon-not-xml/sto-nosoc.dyn", nqirr=10)

Then, SSCHA allows you to impose symmetry and eliminate the soft modes,

 dyn.Symmetrize()
 dyn.ForcePositiveDefinite()

Finally, we are able to create the NVT ensembles with a given temperature and number of configurations.

 # Set the temperature
 ensemble = sscha.Ensemble.Ensemble(dyn, T0=200, supercell=dyn.GetSupercell())
 
 # Set the number of configs to generate
 ensemble.generate(N=100)
 
 # Set the saving directory for the configs
 ensemble.save("data_ensemble_manual", population = 1)

You can execute the script to finish all the steps,

python 1-sscha-generate.py

Step 2: Simulate the configurations in the NVT ensembles

The script we provided for this step is merely a demo and highly dependent on the cluster where you submit the jobs. After proper modification, submit the job, like the command below for SLURM,

sbatch 2-config-submit.slurm

Note that scf.in is already provided in the example folder example07-sto-tdep-sscha/sscha, please notice the differences between this scf.in file and the regular file used for QE calculation without a supercell.

Step 3: Extract the force, stress and energy

This step is just to extract the information of force, stress and energy from the simulations of NVT ensembles, which will be used the for minimization process in next step.

python 3-sscha-collect.py

Step 4: Minimization

The minimization process of the self-consistent harmonic approximation method is done in this step. Please refer to the script 4-sscha-minimize.py for details while we will just emphasize several important parts.

The information of force et. al. will be loaded by

ensemble.load("data_ensemble_manual", population=1, N=100)

Please initialize the minimizer,

minimizer = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)

Please adjust the following two parameters to achieve better convergence,

minimizer.min_step_dyn = 1 
 
minimizer.kong_liu_ratio = 0.5 

The values specified above should serve as a good starting point.

Finally, let’s save the dynamical matrices in the way QE/PERTURBO wants

minimizer.dyn.save_qe("dyn_pop1_")
minimizer.dyn.save_qe_xml("../pw-ph-wann/phonon/save/sto-nosoc.dyn", "dyn_pop3_")

The first parameter of the last line serves as a template for the script to create output files.

Though this step looks a bit complicated, we provide a script as well,

python 4-sscha-minimize.py

Step 5: Self-consistent iteration

Generally, we need to repeat the entire process multiple times to achieve convergence. You can execute the fifth script to move essential files to a new folder from the previous running.

bash 5-new-phonon-folder.sh

And then you are going to repeat Step 1-4, until the final convergence of the phonon dispersions.

After you get convergent results, you can specify the phdir in qe2pert.in as that “save” folder obtained by SSCHA. Everything else is the same as the regular workflow of PERTURBO.