# Thermal conductivity of CNT

Dear lammps users,

Hello! I’m working on thermal conductivity prediction of Single Wall Carbon Nanotube by Charmm forcefiled. I equilibrated the CNT systems at first, and then did the thermal conductivity prediction by Green-Kubo method. However, I got a really small thermal conductivity of CNT. The thermal conductivity predicted result is only 0.6[W/mK]. I pasted my scripted at the end . In the thermal conductivity calculation. I tried different correlation length. However , it gave similar result about 0.6~10[W/mK]. I also tried to used Tersoff forcefield to do thermal conductivity calculation. But the result is still similar. I’m trying to find the problem in my scripts. Can anyone give me some comments?

So much thanks for the help.

Kerwin

Thermal conductivity calculation script in charmm forcefield:

variable T equal 273
variable dt equal 0.5
variable V equal vol
variable p equal 1000 # correlation length
variable s equal 10 # sample interval
variable d equal \$p*\$s # dump interval

#convert from Lammps real units to SI

variable kB equal 1.3806504e-23 # [J/K] Boltzmann
variable kCal2J equal 4186.0/6.02214e23
variable A2m equal 1.0e-10
variable fs2s equal 1.0e-15
variable convert equal {kCal2J}*{kCal2J}/{fs2s}/{A2m}

# setup problem

units real
dimension 3
atom_style full

neighbor 2.0 bin
neigh_modify every 2 delay 0 check yes page 100000

pair_style lj/charmm/coul/long 10.00 10.10
kspace_style pppm 1.0e-6
bond_style harmonic
angle_style charmm
dihedral_style charmm
improper_style harmonic
kspace_modify gewald 0.225
atom_modify sort 0 0.0

velocity all scale \$T

# thermal conductivity calculation

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/vol
variable Jy equal c_flux/vol
variable Jz equal c_flux/vol

fix 1 all nve
fix JJ all ave/correlate \$s \$p d & c_flux c_flux c_flux type auto & file profile.heatflux-run02 ave running dump 2 all xyz 1000 xyz-2000.file-run02 dump 3 all custom 1000 dump_vmd-thermal-run02.lammpstrj id mol type x y z ix iy iz variable scale equal {convert}/\${kB}/\$T/\$T/\$V*s*{dt}
variable k11 equal trap(f_JJ){scale} variable k22 equal trap(f_JJ)*{scale}
variable k33 equal trap(f_JJ)
\${scale}

thermo \$d
thermo_style custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33 density press etotal

timestep 0.5
restart \$d heatflux-1-run02 heatflux-2-run02
run 5000000

variable k equal (v_k11+v_k22+v_k33)/3.0
variable ndens equal count(all)/vol
print “average conductivity: \$k[W/mK] @ T K, {ndens} /A^3”

reset_timestep 0

Thermal conductivity calculation script by tersoff forcefield:

# Variables

variable input index lammps_data01

# INITIALIZATION

units metal
boundary p p p
atom_style atomic

variable T equal 300.0
variable dt equal 0.0002 # ps
variable p equal 1000 # correlation length
variable s equal 10 # sample interval
variable d equal \$p*\$s # dump interval

thermo 100000
timestep \${dt} # ps

# FORCE FIELD

pair_style tersoff
pair_coeff * * /home/tntech.edu/czhan42/lammps-16Mar18/potentials/SiC.tersoff C

# SETTINGS

neighbor 2.0 bin
neigh_modify delay 1

# Calculating Volume

variable diameter equal (sqrt(3)1.42sqrt(6^2+6*6+6^2))/PI #(6,6) armchair CNT; assuming C-C bond distance of 1.42
variable radius equal {diameter}/2 variable height equal "lz" #Height is equal to z-dimension length variable V equal 2*PI*{radius}3.4\${height}

# Convert from metal units to SI

variable kB equal 1.3806504e-23 # [J/K] Boltzmann
variable eV2J equal 1.602176565e-19
variable A2m equal 1.0e-10
variable ps2s equal 1.0e-12
variable convert equal {eV2J}*{eV2J}/{ps2s}/{A2m}

group Carbon type 1

# ------------------------------------------------------------------------------

variable temp equal temp

# output for VMD

dump 1 all custom 10000 dump.uhmwpe-coord id type x y z
dump 2 all xyz 10000 dump.uhmwpe-xyz-format
dump_modify 2 element C

## — Relaxation —

fix 51 all nvt temp \$T \$T 0.1
run 100000

## — Thermal Conductivity Calculation —

unfix 51
fix 51 all nve
thermo d
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/\$V
variable Jy equal c_flux/\$V
variable Jz equal c_flux/\$V
fix JJ all ave/correlate \$s \$p d & c_flux c_flux c_flux type auto file profile.dat ave running variable scale equal {convert}/\${kB}/\$T/\$T/\$V*s*{dt}

variable k11 equal trap(f_JJ){scale} variable k22 equal trap(f_JJ)*{scale}
variable k33 equal trap(f_JJ)
\${scale}
thermo_style custom step temp cpu v_Jx v_Jy v_Jz v_k11 v_k22 v_k33 density press etotal
run 5000000

variable k equal (v_k11+v_k22+v_k33)/3.0
print “average conductivity: \$k[W/mK] @ \$T K” Hello,

I have some comments on some of your input statements (the Tersoff version) below.

variable dt equal 0.0002 # ps

A time step of 0.5 fs (or even 1 fs) is ok for Tersoff carbon.

variable p equal 1000 # correlation length
variable s equal 10 # sample interval
variable d equal \$p*\$s # dump interval

So your correlation time is only 1010000.2fs = 2 ps. For Tersoff CNT, it needs a correlation time of about 1 ns (1000 ps) to achieve a converged thermal conductivity. Also, it is better to make \$d larger because it corresponds to the production time for each correlation function. Good parameters for Tersoff carbon:
timestep=0.001 # time step = 1 fs
s=10 # sampling interval = 10 fs
p=10^5 # correlation time = 1 ns
d=10^7 # production time = 10 ns
Also you need to do at least 10 independent runs. Otherwise, you cannot expect to get statistically meaningful results. Thermal conductivity of CNT is very difficult to calculate. It takes a total production of about 1000 ns to obtain an accurate estimation of the thermal conductivity using Green-Kubo.

pair_coeff * * /home/tntech.edu/czhan42/lammps-16Mar18/potentials/SiC.tersoff C

As far as I know, the version by Lindsay and Broido is currently the best potential for thermal transport: https://doi.org/10.1103/PhysRevB.81.205441

compute flux all heat/flux myKE myPE myStress
variable Jx equal c_flux/\$V
variable Jy equal c_flux/\$V
variable Jz equal c_flux/\$V
*fix JJ all ave/correlate \$s \$p d &** **c_flux c_flux c_flux type auto file profile.dat ave running** __variable scale equal {convert}/\${kB}/\$T/\$T/\$Vs*{dt}__

Here, you have put the volume V in the wrong position. It should be in the numerator instead of the denominator, because you have already divided J by V. If you have not divided J by V, it is correct to put V in the denominator.

variable k equal (v_k11+v_k22+v_k33)/3.0
You cannot use this formula for quasi-one-dimensional CNT. This is only for 3D isotropic systems.

print "average conductivity: \$k[W/mK] @ \$T K"

Did you report your thermal conductivity value based on this printed number? You need to analyse the heat current autocorrelation function and the corresponding running thermal conductivity carefully. This single number is not very meaningful.

Best,
Zheyong