# ShengBTE

## Introduction

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:

ShengBTE：https://bitbucket.org/sousaw/shengbte/src/master/

phonopy：http://abacus.deepmodeling.com/en/latest/advanced/interface/phonopy.html

ASE：http://abacus.deepmodeling.com/en/latest/advanced/interface/ase.html

thirdorder: https://bitbucket.org/sousaw/thirdorder/src/master/

## Prepare

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: `2nd`

, `3rd`

, and `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 (`shengbte`

).

## How to use

Taking the `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 `2nd`

folder.

### 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
```

where the `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 `FORCE_SET`

file:

```
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
```

The `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
```

### 1.3 Post-processing

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:

```
python au2si.py
```

The `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 `POSCAR`

and `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
```

Run `pos2stru.py`

to convert the above `POSCAR`

to `STRU`

file. Note that this script calls functions from the ASE software package (ASE needs to be installed in advance):

```
python pos2stru.py
```

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:

```
python aba2vasp.py
```

The `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

Enter the `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
```

## Conclusion

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.