The ab inito calculation of thermal conductivity is very
complex and needs several softwares to work together. The basic idea is
to calculate the interatomic force constants (IFCs) for both the
harmonic and anharmonic terms, and then use the IFCs to calculate the
phonon dispersion relation and phonon scattering time using knowledges
that are summarized in the previous posts about anharmonicity and phonon dispersion.
Workflow
Overall, we use the following open-source softwares:
phonopy that helps
to generate the supercell we want.
Quantum Espresso
that carries out the ab initio calculation of electronic
structures.
alamode
which fits the IFCs from the ab initio calculation and
calculates the phonon dispersion relation and phonon scattering time and
finally thermal conductivity.
And the whole process would be like this:
graph TD
A["Construct unit cell"] -- Quantum Espresso --> B["Test convergence"]
A -- Quantum Espresso --> B2["Solve lattice constant for zero pressure"]
A -- phonopy --> C["Construct supercell"]
C -- alamode/alm/suggest --> D["Figure out displace patterns"]
D -- Quantum Espresso --> E["Calculate forces for each pattern"]
B -- ecutwfc, K_POINTS, conv_thr, etc. --> E
B2 -- a --> E
E -- alamode/alm/optimize --> F["Fit IFCs"]
F -- alamode/anphon/phonons --> H["Calculate phonon dispersion"]
F -- alamode/anphon/RTA --> G["Calculate thermal conductivity"]
Silicon as an example
In this post, we will use Silicon as an example to show how to
proceed the whole workflow, as Silicon crystal is simple and well
studied.
Monocrystalline silicon, with diamond structure, is a face centered
cubic (fcc) lattice with two atoms in the unit cell. The lattice
constant is \(a=5.43\,\mathrm{Ã…}\), or
\(10.26\,\mathrm{Bohr}\) (which is the
edge length of the cube shown in the following figure), and the two
atoms are located at \((0,0,0)\) and
\((0.25,0.25,0.25)\) in the unit
cell.
Construct unit cell
The very first step is to construct a unit cell of Silicon by hand.
Using the information above, we can create a Quantum Espresso input file
si.pw.in with the following content for calculating the
electronic structure as well as the forces:
Specifications of those parameters can be found in the pw.x
manual of Quantum Espresso. It is worth mentioning that we let
ibrav = 0 (free) rather than ibrav = 2 (fcc)
because we want to further generate supercells using
phonopy later, and phonopy only supports the
case of ibrav = 0. The pseudo-potential file
Si.pz-vbc.UPF can be downloaded from Quantum
Espresso website, and we take this one as it is recommended by a
senior researcher in our group. In addition, the 8 atoms in the unit
cell are inputed by hand, taking the nearest 7 neighbors of the \((0,0,0)\) atom.
Also to note that, the parameters tprnfor and
tstress are set to be .true. so that the
forces and stresses can be printed out in the output file, which is
crucial for the calculation of IFCs as well as the preliminary test of
convergence and lattice constant.
We can surely run a calculation by executing the following command,
generating the output file si.pw.out:
1 2 3 4 5
# Run on a single core pw.x < si.pw.in > si.pw.out
# Run on 16 cores mpiexec -n 16 pw.x < si.pw.in > si.pw.out
For users that use clusters with job scheduler, it is important and
recommended to evaluate the memory usage of the calculation, and it can
be done by simply running the above command for several seconds, and
then click Ctrl+C to terminate the calculation. In the file
si.pw.out, an estimated memory usage is indicated by some
lines like this:
si.pw.out
1
Estimated total dynamical RAM > 29.03 GB
Preliminary tests
using Quantum Espresso
Convergence test
The above baseline settings of Silicon conventional cell are actually
easy for Quantum Espresso to solve, and it takes about 1 minute to
finish the calculation on a single core. However, it is not always the
case, as we will solve very large supercell later. Therefore, now it is
the time to take advantage of the very fast speed of the baseline
calculation to test the convergence for the crucial parameters such
as:
ecutwfc: the plane wave cutoff kinetic energy (in Ry).
We change it from \(30\,\mathrm{Ry}\)
to \(150\,\mathrm{Ry}\).
K_POINTS: the number of \(k\) points in the Brillouin zone (in each
direction), we choose it from \(4\times
4\times 4\) to \(20\times 20\times
20\) (all with offsets so that it would be more accurate).
conv_thr: the convergence threshold for the
self-consistent field (SCF) calculation, we change it from \(1\times 10^{-7}\) to \(1\times 10^{-16}\).
Keep in mind that, we care more about the convergence of the forces
rather than the total energy, as the forces are used to calculate the
IFCs later. Therefore, to make the forces exist, we have to manually
perturb the system a bit, say \(0.001\)
crystal unit (which is about \(0.01\)
Bohr), for the first Silicon atom in the cell, by modifying the
ATOMIC_POSITIONS block in si.pw.in:
si.pw.in
1 2 3
ATOMIC_POSITIONS crystal + Si 0.001000000 0.000000000 0.000000000 - Si 0.000000000 0.000000000 0.000000000
To automatically perform the calculation, what I did is to write a
bash script to perform calculation from different configurations, so
that we can be free from the boring and repetitive work.
## conv_thr test cp si.pw.in si.pw.thr.in for ((thr = 7; thr <= 16; thr += 1)); do sed -i "s/conv_thr = 1.0d-[0-9]\+/conv_thr = 1.0d-$thr/" si.pw.thr.in echo"let conv_thr = 1.0d-$thr, calculating..." mpiexec -np 16 pw.x < si.pw.thr.in > ./out/si.pw.thr.$thr.out done
rm si.pw.thr.in
We perform the convergence test, and plot the results below. We
marked the zone in which the largest atom force is in a range within
\(10^{-6}\,\mathrm{Ry/Bohr}\) comparing
to the converged result (forces are marked by purple lines). Also,
variations of the total energy are marked by dark blue lines, and the
dashed zone indicates the range within \(10^{-3}\,\mathrm{Ry}\) comparing to the
converged result.
Convergence test
It can be noticed that:
ecutwfc=70 or larger is far enough for the
calculation.
K_POINTS setting as 8 8 8 1 1 1 is enough
for the calculation of a conventional cell. For supercells, we may
reduce it to 4 4 4 1 1 1 (\(2\times 2\times 2\) supercell) or
3 3 3 1 1 1 (\(3\times 3 \times
3\) supercell), to make sure that the density of \(k\) points remains the same or becomes even
larger.
conv_thr=1.0d-11 is way enough.
Lattice constant
The lattice constant \(a=10.26\,\mathrm{Bohr}\) that we set
previously, may not be applicable for our problem, as it is a rough
value, and it may not be consistent with the zero-pressure condition. In
fact, the author first set \(a=10.26\,\mathrm{Bohr}\) and carried out
all the remaining calculations, and then found that the thermal
conductivity at room temperature is about \(88\,\mathrm{W/mK}\), rather than \(\sim 140\,\mathrm{W/mK}\) that are reported
by various experiments. I then realized that the stress is about \(-14.08\,\mathrm{GPa}\), meaning that a very
strong force is applied to strentch the solid, leading to the result
that interatomic force is weeker, and thus the thermal conductivity is
smaller.
To obtain a more accurate lattice constant to make sure that the
zero-pressure condition is satisfied, we write a python script to
perform a simple binary search.
Remember, before running the script, you have to copy the
si.pw.in to si.stress.in, and make sure that
the ATOMIC_POSITIONS block is changed back to the original
value without any perturbation. Also, here we'd better consistent values
of ecutwfc, K_POINTS, and
conv_thr with the further calculations. Here we just take
70, 8 8 8 1 1 1, and 1.0d-11,
respectively.
Anyway, the script will yield the result that the lattice constant is
10.20753662109375 Bohr, and we just take it as the new
lattice constant.
Construct supercell
and displace patterns
Generate supercell using
phonopy
It is already quite boring to write the baseline input file for a
cell of 8 atoms, and it would be impoosible to imagine that we have to
write the input file for a cell of 64 atoms (\(2^3\times 8\)), or even 216 atoms (\(3^3\times 8\)). Fortunately,
phonopy provides a very convenient way to generate
supercells by one line of command. With it installed, simply run the
following command:
1
phonopy --qe -d --dim="2 2 2" -c si.pw.in
Once again, don't forget to check if you changed
si.pw.in to the unperturbed settings. Also, the
--dim option can be whatever you want. For silicon,
"2 2 2" (Si222) is pretty good, but Si333 will be better if
you have the resources and time to do the calculations. Copy and paste
the contents of generated files into a new file si222.pw.in
and modify prefix, nat,
CELL_PARAMETERS and K_POINTS, like the
following:
&CONTROL calculation = 'scf' prefix = 'si222' pseudo_dir = './pseudo/' outdir = './tmp/' verbosity = 'high' tprnfor = .true. tstress = .true. disk_io = 'none' / &SYSTEM ibrav = 0 nat = 64 ntyp = 1 ecutwfc = 70 / &ELECTRONS conv_thr = 1.0d-11 mixing_beta = 0.7 / CELL_PARAMETERS bohr 20.4150732421875 0.0000000000000 0.0000000000000 0.0000000000000 20.4150732421875 0.0000000000000 0.0000000000000 0.0000000000000 20.4150732421875 ATOMIC_SPECIES Si 28.0855 Si.pz-vbc.UPF ATOMIC_POSITIONS crystal Si 0.0000000000000000 0.0000000000000000 0.0000000000000000 Si 0.5000000000000000 0.0000000000000000 0.0000000000000000 Si 0.0000000000000000 0.5000000000000000 0.0000000000000000 Si 0.5000000000000000 0.5000000000000000 0.0000000000000000 Si 0.0000000000000000 0.0000000000000000 0.5000000000000000 Si 0.5000000000000000 0.0000000000000000 0.5000000000000000 Si 0.0000000000000000 0.5000000000000000 0.5000000000000000 Si 0.5000000000000000 0.5000000000000000 0.5000000000000000 Si 0.0000000000000000 0.2500000000000000 0.2500000000000000 Si 0.5000000000000000 0.2500000000000000 0.2500000000000000 Si 0.0000000000000000 0.7500000000000000 0.2500000000000000 Si 0.5000000000000000 0.7500000000000000 0.2500000000000000 Si 0.0000000000000000 0.2500000000000000 0.7500000000000000 Si 0.5000000000000000 0.2500000000000000 0.7500000000000000 Si 0.0000000000000000 0.7500000000000000 0.7500000000000000 Si 0.5000000000000000 0.7500000000000000 0.7500000000000000 Si 0.2500000000000000 0.0000000000000000 0.2500000000000000 Si 0.7500000000000000 0.0000000000000000 0.2500000000000000 Si 0.2500000000000000 0.5000000000000000 0.2500000000000000 Si 0.7500000000000000 0.5000000000000000 0.2500000000000000 Si 0.2500000000000000 0.0000000000000000 0.7500000000000000 Si 0.7500000000000000 0.0000000000000000 0.7500000000000000 Si 0.2500000000000000 0.5000000000000000 0.7500000000000000 Si 0.7500000000000000 0.5000000000000000 0.7500000000000000 Si 0.2500000000000000 0.2500000000000000 0.0000000000000000 Si 0.7500000000000000 0.2500000000000000 0.0000000000000000 Si 0.2500000000000000 0.7500000000000000 0.0000000000000000 Si 0.7500000000000000 0.7500000000000000 0.0000000000000000 Si 0.2500000000000000 0.2500000000000000 0.5000000000000000 Si 0.7500000000000000 0.2500000000000000 0.5000000000000000 Si 0.2500000000000000 0.7500000000000000 0.5000000000000000 Si 0.7500000000000000 0.7500000000000000 0.5000000000000000 Si 0.1250000000000000 0.1250000000000000 0.1250000000000000 Si 0.6250000000000000 0.1250000000000000 0.1250000000000000 Si 0.1250000000000000 0.6250000000000000 0.1250000000000000 Si 0.6250000000000000 0.6250000000000000 0.1250000000000000 Si 0.1250000000000000 0.1250000000000000 0.6250000000000000 Si 0.6250000000000000 0.1250000000000000 0.6250000000000000 Si 0.1250000000000000 0.6250000000000000 0.6250000000000000 Si 0.6250000000000000 0.6250000000000000 0.6250000000000000 Si 0.1250000000000000 0.3750000000000000 0.3750000000000000 Si 0.6250000000000000 0.3750000000000000 0.3750000000000000 Si 0.1250000000000000 0.8750000000000000 0.3750000000000000 Si 0.6250000000000000 0.8750000000000000 0.3750000000000000 Si 0.1250000000000000 0.3750000000000000 0.8750000000000000 Si 0.6250000000000000 0.3750000000000000 0.8750000000000000 Si 0.1250000000000000 0.8750000000000000 0.8750000000000000 Si 0.6250000000000000 0.8750000000000000 0.8750000000000000 Si 0.3750000000000000 0.1250000000000000 0.3750000000000000 Si 0.8750000000000000 0.1250000000000000 0.3750000000000000 Si 0.3750000000000000 0.6250000000000000 0.3750000000000000 Si 0.8750000000000000 0.6250000000000000 0.3750000000000000 Si 0.3750000000000000 0.1250000000000000 0.8750000000000000 Si 0.8750000000000000 0.1250000000000000 0.8750000000000000 Si 0.3750000000000000 0.6250000000000000 0.8750000000000000 Si 0.8750000000000000 0.6250000000000000 0.8750000000000000 Si 0.3750000000000000 0.3750000000000000 0.1250000000000000 Si 0.8750000000000000 0.3750000000000000 0.1250000000000000 Si 0.3750000000000000 0.8750000000000000 0.1250000000000000 Si 0.8750000000000000 0.8750000000000000 0.1250000000000000 Si 0.3750000000000000 0.3750000000000000 0.6250000000000000 Si 0.8750000000000000 0.3750000000000000 0.6250000000000000 Si 0.3750000000000000 0.8750000000000000 0.6250000000000000 Si 0.8750000000000000 0.8750000000000000 0.6250000000000000 K_POINTS automatic 4 4 4 1 1 1
Obtain displace patterns
The next step is to obtain the displace patterns for the supercell.
Create a new file si222.alm.in that will be used by
alamode to generate the patterns:
Here, the mode suggest means that we ask
alamode to generate patterns. And we set the cutoff
parameters empirically. After this, execute the following command:
1
alm si222.alm.in > si222.alm.log
The harmonic and anharmonic patterns will be saved as files
si222.pattern_HARMONIC and
si222.pattern_ANHARM3, respectively. With the anharmonic
patterns (16 patterns inside the file), we perturb again the input file
of Quantum Espresso, but this time automatically by using the
displace.py script provided by alamode in its
tools directory.
The --mag option specifies the magnitude, and the
-pf option specifies the pattern file. For Si222 (and also
Si333), 16 new input files named from disp01.pw.in to
disp16.pw.in will be generated.
Calculate forces using
pw.x
Just perform the calculations of all the input files generated by
displace.py! To be a bit lazy, I submit all the jobs using
a bash script that run through all the input files and submit them to
the queue system of the cluster in my university. It is something like
this:
_si.disp.submit.sh
1 2 3 4 5 6 7 8 9
for ((i=1;i<=16;i++)) do num=`echo$i | awk '{printf("%02d",$1)}'` cp _si.disp.sh _si.disp${num}.sh # match patterns like "input=disp.pw*.in", "output=disp.pw*.out", change to "input=disp.pw.01.in", "output=disp.pw.01.out" sed -i "s/input=disp.pw.in/input=disp${num}.pw.in/g" _si.disp${num}.sh sed -i "s/output=disp.pw.out/output=disp${num}.pw.out/g" _si.disp${num}.sh qsub _si.disp${num}.sh done
which manipulates a submission script _si.disp.sh:
The following script is just an example, and is only suitable for the
cluster of my university. You have to modify it according to your own
cluster. Or you can simply run the calculations on your own computer
manually...
#### submit_job.sh START #### #!/bin/bash #$ -cwd # error = Merged with joblog #$ -o joblog.$JOB_ID #$ -j y ## Edit the line below as needed: #$ -l h_rt=24:00:00,h_data=3G ## Modify the parallel environment ## and the number of cores as needed: #$ -pe shared 16
# echo job info on joblog: echo"Job $JOB_ID started on: " `hostname -s` echo"Job $JOB_ID started on: " `date ` echo" "
# load the job environment: . /u/local/Modules/default/init/modules.sh
## Edit the line below as needed: module load intel/2020.4 module load gcc/10.2.0 module load espresso/6.7.0
echo"Current Packages:"
module list
echo" " echo"Start running..."
## substitute the command to run your code ## in the two lines below: cd /PATH/TO/THE/WORKING/DIRECTORY/
input=disp.pw.in output=disp.pw.out
mpirun -np 16 pw.x < ${input} > ${output}
# try to determine if there is "JOB DONE" in the output file, if not, echo a warning message if grep -q "JOB DONE"${output}; then echo"PW finished successfully" echo"Job $JOB_ID for $input ended succesffully" >> __submit_job.log else echo"PW may not finished successfully" echo"Job $JOB_ID for $input FAILED!" >> __submit_job.log fi
# echo job info on joblog: echo"Job $JOB_ID ended on: " `hostname -s` echo"Job $JOB_ID ended on: " `date ` echo" " #### submit_job.sh STOP ####
When all the calculations are finished, we extract the forces from
the output files by using the extract.py script that is
also provided by alamode:
&kpoint 1 # KPMODE = 1: line mode G 0.0 0.0 0.0 X 0.5 0.5 0.0 51 X 0.5 0.5 1.0 G 0.0 0.0 0.0 51 G 0.0 0.0 0.0 L 0.5 0.5 0.5 51 /
It is worth noting that the kpoint field is selected to
be the line mode, meaning that anphon will calculate phonon dispersion
relation along the lines that are specified.
Perform the calculation by:
1
anphon si222.phonon.in > si222.phonon.log
Using another script plotband.py provided by
alamode, the phonon band structure can be quickly
plotted:
1
python plotband.py si222.bands
As for the density of states, just change the kpoint
field to:
si222.phonon.in
1 2 3 4 5 6 7 8
&kpoint - 1 # KPMODE = 1: line mode - G 0.0 0.0 0.0 X 0.5 0.5 0.0 51 - X 0.5 0.5 1.0 G 0.0 0.0 0.0 51 - G 0.0 0.0 0.0 L 0.5 0.5 0.5 51 + 2 + 20 20 20 /
perform the same command again, anphon will create
si222.dos and si.thermo that are the density
of states and the thermodynamic properties. To see the DoS, just use
plotdos.py, with the following command:
1
python plotdos.py --emax 550 --nokey si222.dos
Calculate thermal
conductivity
Now it is time to calculate the thermal conductivity. We still use
anphon, and let's create a config file
si222.RTA.in similar to the one of phonon dispersion
relation:
where mode=RTA means that we are going to use
relaxation-time approximation (RTA) to calculate the thermal
conductivity. The FCSXML option specifies the IFCs file
that we just obtained. The cell block, as you can see,
specifies a primitive cell of Silicon monocrystalline, and it is surely
different than the supercell that we used above. The kpoint
block specifies the number of wavevectors in each direction of the
reciprocal space, and it can be set to quite large comparing with that
of the supercell.
Perform the calculation by:
Although the official documentation of alamode claims
that anphon does support MPI, the parallel execution was
not working on both my computer (M1 Max Macbook Pro) and on the cluster
of my university (Intel x86). Luckily, it is not that slow to run on a
single core. Just be patient for a while.
1
anphon si222.RTA.in > si222.RTA.log
will give the thermal conductivity up to 1000 K, stored in
si222.kl, and all other information is stored in
si222.result. Here we plot the thermal conductivity, and
compare it with the data in literature, shown below:
If we redo all the stuff above for Si333, a nearly same result will
be obtained, and we plot the data from 250 K to 450 K, it is shown that
the difference between the calculated data for Si222 and Si333 is within
4.22 W/mK, and is about 2.5% of the value of Si222.
It is so amazing that the calculated thermal conductivity is just in
very good agreement with the experimental data. This is the power of
physics, that we can absoultely calculate anything in the ab
initio way, and they just do the job!
