Calculate thermal conductivity by Green-Kubo method
If the thermal conductivity of a material is desired (at least some trends such as the temperature-dependency trend), ab initio calcultions are commonly the best choice due to its accuracy. However, for complex structures, it is no longer feasible, as the number of displacement patterns or/and the number of interatomic force constants can be too large to be calculated or fitted. In this case, to simply obtain the thermal conductivity, the Green-Kubo method is on the stage, which is an equilibrium method that uses only fluctuations of the heat current to calculate the thermal conductivity tensor:
\[ \kappa_{\mu\nu}(t) = \frac{1}{k_{\text{B}}T^2 V}\int_0^t \Braket{J_{\mu}(0)J_{\nu}(t')}dt' \]
where \(J_{\mu}\) is the heat current in the \(\mu\) direction. In this post, we will calculate the thermal conductivity of Silicon again, but this time by the EMD method.
LAMMPS provides a simple command to calculate the heat current \(J_{\mu}\), that is compute heat/flux
,
with the output of which the thermal conductivity can be easily
calculated. But to do so, basic knowledge of LAMMPS is required.
Input file: a script of commands
Different from other softwares such as Quantum ESPRESSO and VASP that
take the input file as a collections of parameters and values, LAMMPS
takes the input file as a script that contains a series of commands,
just like the shell script, and the commands are executed in the order
they appear in the script (Due to the similarity with bash scripts, I
personally prefer the text editors to highlight .lammps
files as bash scripts). For example, the following script is an example
of a LAMMPS input file that simulates a Lennard-Jones
fluid containing two kinds of atoms:
1 | # PART A - ENERGY MINIMIZATION |
The script is divided into two parts, the first part is the definition of the system and the energy minimization, and the second part is the molecular dynamics simulation. As you may notice from the script, a common LAMMPS calculation is composed of the following steps:
- Initialization: Set the units, dimensions, atom style, pair style, and boundary conditions. It contains the basic styles of the simulation.
- System definition: Define the simulation domain and atoms.
- Simulation settings: Define interatomic potentials, masses, and other parameters.
- Visualization: Set the output format and frequency of the desired output. It is not necessary, but is helpful for monitoring the simulation and depends on how you will analyze the simulation.
- Run: Run the simulation under certain conditions and for a certain number of steps.
Also, the script eventually outputs a file
dump.lammpstrj
that contains the trajectory of the
simulation, and the file can be visualized by some visualization
softwares such as VMD:
![](/images/2024/VMD.jpg)
Silicon crystal and lattice relaxation
Let's back to the topic. To simulate the thermal conductivity of Silicon, a crystal structure is required to be defined. The following script is used to (1) define styles and the system, and (2) thermalize the system to a certain temperature (300 K):
1 | units metal |
Here, I used the EDIP potential, as the others are not so accurate for the thermal conductivity. The lattice constant is a result of relaxation, which is a little bit tricky for MD, as the stress of a system fluctuactes too much. One possible stable way is to do a average calculation of the stress for different lattice constants:
1 | # Compute the average pressure |
and use a binary method to find the lattice constant that gives the average pressure closest to zero. For this case, I obtained a lattice constant of 5.433144532 Ã… at 300 K.
Heat current calculation
After the lattice relaxation, we may perform the following commands
to simulate the fluctuations of the heat current and record everything
in the log.lammps
file:
1 | reset_timestep 0 |
Here, the heat current is calculated by the
compute heat/flux
command, and corresponds to
Jx
, Jy
, and Jz
in the
log.lammps
file. The results are wrapped by the
====STARTING SIMULATION====
and
====ENDING SIMULATION====
lines, which helps us to locate
the results using either the grep
command or some glue
post-processing scripts. Another more direct way is to just calculate
everything inside LAMMPS, as shown by the documentation
of the compute heat/flux
command, but I prefer to print
everything and do the post-processing by Python to obtain a better
visualization.
Post-processing
We perform the calculation three times. As silicon is symmetric in the \(x\), \(y\), and \(z\) directions, we then have nine curves of the heat current autocorrelation. We change the upper limit \(t\) of the integral, and for each \(t'\), the integrand \(\Braket{J_{\mu}(0)J_{\nu}(t')}\) is calculated by taking the average within the whole 30 ns simulation time. The results are shown as follows:
![](/images/2024/k_EMD.png)
After 0.2 ns, the average of the nine curves becomes noisy and instable, but the results are at least very close to the saturated value, as seen from the figure. It is noteworthy that EMD actually takes a lot of time to obtain a stable and accurate result. For simple structures like silicon crystal, it is surely better to do an ab initio calculation, as EMD takes more time and is not accurate. But for complex structures, EMD is the only choice.