Calculate thermal conductivity by Green-Kubo method

Posted on In , Word count in article: 1.2k

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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
# PART A - ENERGY MINIMIZATION
# 1) Initialization
units lj
dimension 3
atom_style atomic
pair_style lj/cut 2.5
boundary p p p

# 2) System definition
region simulation_box block -20 20 -20 20 -20 20
create_box 2 simulation_box
create_atoms 1 random 1500 341341 simulation_box
create_atoms 2 random 100 127569 simulation_box

# 3) Simulation settings
mass 1 1
mass 2 1
pair_coeff 1 1 1.0 1.0
pair_coeff 2 2 0.5 3.0

# 4) Visualization
thermo 10
thermo_style custom step temp pe ke etotal press

# 5) Run
minimize 1.0e-4 1.0e-6 1000 10000

# PART B - MOLECULAR DYNAMICS
# 4) Visualization
thermo 50
dump mydmp all atom 100 dump.lammpstrj

# 5) Run
fix mynve all nve
fix mylgv all langevin 1.0 1.0 0.1 1530917
timestep 0.005
run 10000

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:

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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
units       metal
variable    T equal 300
variable    dt equal 0.001
variable    s equal 10

dimension    3
boundary     p p p
lattice      diamond 5.433144532
region       box block 0 10 0 10 0 10
create_box   1 box
create_atoms 1 box
mass         1 28.0855
pair_style   edip
pair_coeff   * * ../Si.edip Si
timestep     ${dt}

velocity     all create $T 12345 mom yes rot no dist gaussian
fix          NVT all nvt temp $T $T $(100.0*dt) drag 0.2
thermo       100
thermo_style custom step temp press
run          1000

reset_timestep 0

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
2
3
4
5
6
# Compute the average pressure

fix          ave all ave/time 10 200 2000 c_thermo_temp c_thermo_press ave running
thermo_style custom step f_ave[1] f_ave[2]
thermo       2000
run          100000

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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
reset_timestep 0
compute      myKE all ke/atom
compute      myPE all pe/atom
compute      myStress all stress/atom NULL virial
compute      flux all heat/flux myKE myPE myStress
variable     Jx equal c_flux[1]/vol
variable     Jy equal c_flux[2]/vol
variable     Jz equal c_flux[3]/vol

thermo_style custom step temp v_Jx v_Jy v_Jz press

thermo       $s

print        "====STARTING SIMULATION===="
run          30000000   # 30ns
print        "====ENDING SIMULATION===="

variable     ndens equal count(all)/vol
print        "Volume: ${V} [A^3]"
print        "Number density: ${ndens} /A^3"

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:

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.