Tutorial: Thermal conductivity from EMD

Jump to navigation Jump to search


In this example, we use the EMD (Green-Kubo) method to calculate the lattice thermal conductivity of graphene at 300 K and zero pressure.

Preparing the Inputs

The xyz.in file

  • Note that the thickness of the graphene sheet is set to 3.35 A according to the convention in the literature. This thickness is needed to calculate an effective 3D thermal conductivity for a 2D material.

The run.in file

  • The run.in input file is given below:
   potential    potentials/tersoff/Graphene_Lindsay_2010_modified.txt
   velocity     300
   # equilibration
   ensemble     npt_ber 300 300 0.01 0 0 0 0.0005
   time_step    1
   dump_thermo  1000
   run          1000000
   # production
   ensemble     nve
   compute_hac  20 50000 10
   run          10000000
  • The second line of the command with the velocity keyword tells that the velocities will be initialized with a temperature of 300 K.
  • There are two runs.
    • The first run serves as the equilibration stage, where the NPT ensemble (the Berendsen barostat) is used. The temperature is 300 K and the pressures are zero in all the directions. The coupling constants are 0.01 (dimensionless) and 0.0005 (in the natural unit system adopted by GPUMD) for the thermostat and the barostat, respectively. The time step for integration is 1 fs. There are [math]10^6[/math] steps (1 ns) for this run and the thermodynamic quantities will be output every 1000 steps.
    • The second run is for production. Here, the NVE ensemble is used. The line with the compute_hac keyword means that heat currents will be recorded every 20 steps (20 fs), 50000 HAC data (the maximum correlation time is then about 1 ns) will be calculated, and the HAC are averaged for every 10 data points before written out. The production time is 10 ns ([math]10^7[/math] steps), which is 10 times as long as the maximum correlation time. This is a reasonable choice.

Results and discussion


Thermal conductivity results for pristine graphene at 300 K from EMD simulations. (a) Normalized HAC as a function of correlation time for the in-plane and out-of-plane components. (b) Individual (thin lines) and averaged (thick line) RTC as a function of correlation time for the in-plane component. (c) Individual (thin lines) and averaged (thick line) RTC as a function of correlation time for the out-of-plane component. (d) Averaged RTC as a function of correlation time for various components.
  • The above figure shows the results from three independent simulations, which took about two hours in total using a Tesla K40 card.
  • As the system is essentially isotropic in the planar directions, we can average over the two in-plane directions.
  • From (a), we can see that the in-plane component and the out-of-plane component of the HAC have different time scales. The latter decays much more slowly.
  • Panel (b) shows the individual and averaged RTCs for the in-plane component [math]\kappa^{\rm in}(t)[/math]. The averaged RTC converges to about 1000 W/mK at around 100 ps.
  • Panel (c) shows the individual and averaged RTCs for the out-of-plane component [math]\kappa^{\rm out}(t)[/math], and the convergence property is not very clear here. This is because the out-of-plane component converges very slowly and three independent simulations (each with 10 ns) are not enough to give accurate results.
  • Summing up [math]\kappa^{\rm in}(t)[/math] and [math]\kappa^{\rm out}(t)[/math], we get [math]\kappa^{\rm tot}(t)[/math], as shown in panel (d).


  • Accurately calculating thermal conductivity of graphene using the EMD method can be a very time consuming task. The results we presented here are from three independent simulations with a total production time of 30 ns. It can been seen that the HAC data already become very noisy when the correlation time is 100 ps. To obtain accurate results, one needs to do many independent simulations. Much more accurate data were presented in Fig. 2 of this paper. Here are the simulation parameters used in the reference paper which differ from those used in this example:
    • The simulation cell size used in the reference paper is larger, which is about 25 nm x 25 nm (24000 atoms), instead of 15 nm x 15 nm (8640 atoms) here.
    • The maximum correlation time used in the reference paper is larger, which is 10 ns, instead of 1 ns here.
    • The production time used in the reference paper for one independent simulation is larger, which is 50 ns, instead of 10 ns here.
    • There were 100 independent simulations in the reference paper, instead of three here.
    • Each independent simulation in the reference paper took about 10 GPU hours (using Tesla K40) and about 1000 GPU hours were used to obtain the results shown in Fig. 2 of the reference paper.

Related pages