This tutorial aims to introduce the process of performing density functional theory calculations using ABACUS and calculating lattice thermal conductivity using the ShengBTE software. During the entire calculation process, the following external software are also used: 1) Phonopy for calculating second-order force constants, 2) ASE for converting atomic structures, 3) ShengBTE’s thirdorder program for calculating third-order force constants, and 4) finally using ShengBTE to calculate the material’s lattice thermal conductivity.
Here is the announcement of ShengBTE with ABACUS: ShengBTE - The ABACUS DFT software can now be used with ShengBTE
Some external packages that need to be combined are mentioned above, and here it is recommended to read the relevant documentation and instructions of these packages:
The ABACUS software package provides an example of using ABACUS+ShengBTE to calculate the lattice thermal conductivity in the
examples/interface_ShengBTE folder. The example includes two folders:
LCAO (Linear Combination of Atomic Orbitals) which uses numerical atomic orbitals and
PW (Plane wave) which uses plane wave basis vectors. Each folder contains three subfolders:
shengbte, which respectively store the relevant files for calculating second-order force constants using Phonopy (
2nd), third-order force constants using the thirdorder program (
3rd), and lattice thermal conductivity using ShengBTE (
How to use
LCAO folder as an example, we provide the test case of a diamond structure Si structure containing 2 atoms with the norm-conserving pseudopotential
Si_ONCV_PBE-1.0.upf and the atomic orbital file
Si_gga_7au_100Ry_2s2p1d.orb (GGA functional, 7 a.u. cut-off radius, 100 Ry energy cut-off, and DZP orbitals containing 2s2p1d).
1. Calculating the second-order force constants
It would be best to combine Phonopy and ASE with ABACUS to calculate the second-order force constants. First, enter the
1.1 Structure optimization
Before performing lattice thermal conductivity calculations, it is necessary to optimize the atomic configuration of the simulated material system. The following is the atomic configuration file
STRU obtained after structure optimization (relax) using ABACUS. In this example, for simplicity, a 2*2*2 Brillouin zone k-point sampling was used in the structure optimization process, with an energy cutoff value of 100 Ry for plane waves (the plane wave basis vector is also used in LCAO). Note that the actual calculation should use more convergent k-point sampling.
ATOMIC_SPECIES Si 28.0855 Si_ONCV_PBE-1.0.upf NUMERICAL_ORBITAL Si_gga_7au_100Ry_2s2p1d.orb LATTICE_CONSTANT 1.88972612546 LATTICE_VECTORS 0 2.81594778072 2.81594778072 #latvec1 2.81594778072 0 2.81594778072 #latvec2 2.81594778072 2.81594778072 0 #latvec3 ATOMIC_POSITIONS Direct # direct coordinate Si #label 0 #magnetism 2 #number of atoms 0.875 0.875 0.875 m 0 0 0 0.125 0.125 0.125 m 0 0 0
1.2 Calculating the second-order force constants
The Phonopy software is called to generate multiple atomic configurations of the supercell and corresponding perturbations needed for calculation with the following command:
phonopy setting.conf --abacus -d
setting.conf file reads:
DIM = 2 2 2 ATOM_NAME = Si
In this Si example, we only need to generate one perturbed atomic configuration,
STRU-001. Perform SCF calculations (SCF stands for Self-Consistent Field and represents the iterative self-consistent calculation of density functional theory) on all perturbed configurations (in this case, there is only one for Si) to obtain the forces on the atoms. Afterward, use the following command to generate the
phonopy -f OUT.DIA-50/running_scf.log
Tip: In the input file
INPUT of ABACUS, you can set the variable
stru_file, which corresponds to the atomic configuration file
STRU-001, and ABACUS will read the structure file directly.
Next, set the
band.conf file to calculate the phonon spectrum and the second-order force constants:
phonopy -p band.conf --abacus
band.conf file mentioned here contains the following contents (you can refer to the Phonopy documentation for specific parameter meanings):
ATOM_NAME = Si DIM = 2 2 2 MESH = 8 8 8 PRIMITIVE_AXES = 1 0 0 0 1 0 0 0 1 BAND = 0.0 0.0 0.0 0.5 0.0 0.5 0.625 0.25 0.625, 0.375 0.375 0.75 00 0.0 0.0 0.5 0.5 0.5 BAND_POINTS = 101 BAND_CONNECTION = .TRUE. FORCE_CONSTANTS = WRITE FULL_FORCE_CONSTANTS = .TRUE.
After this step, the Phonopy software will generate
band.yaml (for plotting the phonon spectrum) and the
FORCE_CONSTANTS file. The data contained in the
FORCE_CONSTANTS file is the second-order force constants. It is important to set
FULL_FORCE_CONSTANTS = .TRUE., which outputs all the second-order force constants. Otherwise, there may be errors when ShengBTE reads the data.
In addition, you can use the following command to output the gnuplot format of the phonon spectrum for plotting:
phonopy-bandplot --gnuplot > pho.dat
Note that ShengBTE software requires the unit of the data in the
FORCE_CONSTANTS_2ND file to be eV/Å^2, but the unit of the
FORCE_CONSTANTS calculated by ABACUS combined with Phonopy is eV/(Å*au), where au is the atomic unit system and 1 au = 0.52918 Å. You can use the provided
au2si.py script in the
2nd directory to convert the units and generate the
FORCE_CONSTANTS_2ND file. The command is as follows:
FORCE_CONSTANTS_2ND file is provided in the shengbte folder for reference to the calculation results.
2. Calculating the third-order force constants
To calculate the third-order force constants, you need to combine with the thirdorder program and output the third-order force constant file
FORCE_CONSTANTS_3RD. However, thirdorder currently only supports reading input and output files from VASP and QE. Therefore, we are using thirdorder by converting ABACUS’s structure and output force files to
vasprun.xml, respectively. Please enter the
3rd folder first, and the specific steps will be described below.
2.1 Obtaining perturbed configurations
First, convert the optimized
STRU file from ABACUS software to
POSCAR (the converted
POSCAR file is already provided in the directory, or you can do this conversion by yourself).
Then, run the
thirdorder_vasp.py program to generate a series of atomic configuration files
3RD.POSCAR.* after perturbation. For example, in this example, a total of 40 configurations were generated:
thirdorder_vasp.py sow 2 2 2 -2
pos2stru.py to convert the above
STRU file. Note that this script calls functions from the ASE software package (ASE needs to be installed in advance):
Note: The dpdata software cannot be called here to perform the conversion. This is because the dpdata software forces the lattice to change into a lower triangular matrix, which is equivalent to rotating the lattice and leads to a corresponding rotation in the direction of the interatomic forces, which will cause errors.
2.2 Calculation of atomic forces for perturbation configurations
You can refer to the
run_stru.sh script provided in the directory to batch generate
SCF-* folders and submit calculations. Here, ABACUS needs to perform SCF calculations on 40 atomic configurations, which may take some time. It is recommended to run each SCF separately in the
SCF-* folder. The
scf_thr parameter in INPUT file should be set to at least 1e-8 to obtain converged results.
After the calculations are complete, run
aba2vasp.py to package the atomic forces calculated by ABACUS into the
vasprun.xml format and place them in each
SCF-\* folder with the following command:
vasprun.xml format is illustrated as follows:
<modeling> <calculation> <varray name="forces"> <v>1.865e-05 -0.04644196 -0.00153852</v> <v>-1.77e-05 -0.00037715 -0.00149635</v> <v>1.973e-05 0.002213 -0.00149461</v> <v>-1.976e-05 0.00065303 -0.0014804</v> <v>8.31e-06 -0.0003306 -0.00024288</v> <v>-8.25e-06 -0.00038306 -0.00025385</v> <v>1.071e-05 0.00060621 -0.00025797</v> <v>-1.05e-05 -0.00014553 -0.00027532</v> <v>0.00668053 0.00645634 -0.04642593</v> <v>-0.00668085 0.00645595 -0.00040122</v> <v>-0.00650454 0.00628877 -0.00025123</v> <v>0.00650504 0.00628892 -0.00028948</v> <v>-0.00039591 2.479e-05 0.00223371</v> <v>0.00039608 2.426e-05 0.0006732</v> <v>0.0003264 3.122e-05 0.00052874</v> <v>-0.00032589 3.415e-05 -0.00023577</v> <v>-2.908e-05 -0.00832477 0.00635709</v> <v>3.737e-05 -0.00125057 -7.444e-05</v> <v>-2.582e-05 0.00656076 0.00636285</v> <v>2.566e-05 -0.00049974 -6.661e-05</v> <v>-5.431e-05 0.00502637 0.00639077</v> <v>4.553e-05 -0.00180978 0.0001325</v> <v>-3.609e-05 -0.00676473 0.00638092</v> <v>3.806e-05 5.503e-05 0.00012759</v> <v>-0.00670704 0.00646596 0.01310437</v> <v>0.00670119 3.673e-05 0.00602948</v> <v>0.00036366 0.00627899 -0.00657272</v> <v>-0.00036508 2.288e-05 0.00026009</v> <v>0.00648649 0.0064463 -0.00036521</v> <v>-0.00648098 1.594e-05 0.00671469</v> <v>-0.00034493 0.00630074 0.00662932</v> <v>0.00034331 4.157e-05 -0.0002028</v> </varray> </calculation> </modeling>
Finally, execute the following command:
find SCF-* -name vasprun.xml|sort -n|thirdorder_vasp.py reap 2 2 2 -2
Then, the third-order force constant file
FORCE_CONSTANTS_3RD can be obtained by running the above command. The
FORCE_CONSTANTS_3RD file is provided in the shengbte folder for reference in calculating the results.
3. Run ShengBTE to obtain lattice thermal conductivity
shengbte folder, in which the three files
CONTROL (parameter file of ShengBTE),
FORCE_CONSTANTS_2ND (second-order force constant file), and
FORCE_CONSTANTS_3RD (third-order force constant file) have been prepared. Next, run ShengBTE with the following command to obtain the lattice thermal conductivity, where the calculation results are given in the Ref folder for reference:
mpirun -n 10 ShengBTE
For using plane wave (PW) in ABACUS to perform ShengBTE calculations, similar procedures should be followed. However, the
scf_thr parameter in the
INPUT file for calculating the third-order force constant needs to be set to at least 1e-12. The experimental lattice thermal conductivity of Si at 300 K is around 150 W/(m K), while the calculated thermal conductivity of Si at 300 K is around 100 W/(m K) by using the provided example. This is because, as a demo, a 2*2*2 expanded cell and a 2*2*2 K-point are used in the example, but the results are not converged yet with respect to the given system size and k-points. In actual research, the size of the supercell and the sampling scheme of K-points need to be tested to obtain converged results.