Single layer graphene thermal conductivity - small exchanged energy

Dear Lammps users,

Hi, I am trying to calculate thermal conductivity of single layer graphene (100 x 11 x 0.34 nm) using M-P NEMD method. Please find the copied input/output scripts below (highlighted numbers are ‘exchanged energy’, ‘delta T’, and ‘thermal conductivity’, respectively). I am getting exceptionally small exchanged energy between the coldest and hottest slabs, which results in very low thermal conductivity, around 40 W/m-K. Please let me know if you find anything suspicious. Thanks.

Best Regards,

-------input script----------------------------------------------------------------

variable T equal 300
variable V equal vol
variable dt equal 0.001
#variable p equal 1000 # correlation length
#variable s equal 5 # 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 eV2J equal 1.602e-19
variable A2m equal 1.0e-10
variable ps2s equal 1.0e-12
variable convert equal {eV2J}*{eV2J}/{ps2s}/{A2m}
variable alat equal 2.46

setup problem

units metal
boundary p p p
atom_style atomic

lattice custom ${alat} a1 0.86603 0.5 0 a2 0.86603 -0.5 0 a3 0 0 10 basis 0.33333 0.33333 0 basis 0.66667 0.66667 0
region box block 0 1000 0 110 -50 50 units box
create_box 1 box
region tf block INF INF INF INF -5 5 units box
create_atoms 1 region tf
mass 1 12

pair_style tersoff
pair_coeff * * /home/kpark39/programs/lammps-1Feb14/potentials/SiC_1994.tersoff C
neighbor 2.0 bin
neigh_modify delay 5
timestep ${dt}

dump 1 all cfg 10000 dump.config.*.cfg mass type xs ys zs vx vy vz x y z

1st equilibration run

velocity all create $T 4928459 mom yes rot yes dist gaussian units box
fix 1 all npt temp $T $T 1 iso 0 0 0.1 drag 1
thermo 1000
run 50000
write_restart restart.NPT20k.equil
unfix 1
#velocity all scale $t

2nd equilibration run

compute ke all ke/atom

units of Kb = J/K, units of ke = eV or J, so temp becomes K

variable temp atom c_ke*{eV2J}/(1.5*{kB})
fix 1 all nve
fix 2 all ave/spatial 10 100 1000 x lower 5 v_temp file profile.mp units box
fix 3 all thermal/conductivity 10 x 190 swap 2
fix 4 all temp/berendsen 300 300 100
variable tdiff equal f_2[96][3]-f_2[1][3]

thermo_style custom step temp epair etotal f_3 v_tdiff
thermo 1000
dump 2 all cfg 10000 dump.config.*.cfg mass type xs ys zs vx vy vz x y z v_temp
run 50000
unfix 4

thermal conductivity calculation

fix 3 all thermal/conductivity 10 x 190 swap 2
fix ave all ave/time 1 1 1000 v_tdiff ave running file ave.dat
variable K equal f_3*{eV2J}*1000/2/(110*3.4)/{A2m}/{dt}/100000/{ps2s}/f_ave
thermo_style custom step temp epair etotal f_3 v_tdiff f_ave v_K
thermo 1000
run 100000

-------output script----------------------------------------------------------------

LAMMPS (1 Feb 2014)
Lattice spacing in x,y,z = 4.26087 2.46 24.6
Created orthogonal box = (0 0 -50) to (1000 110 50)
16 by 2 by 2 MPI processor grid
Created 42255 atoms
Setting up run …
Memory usage per processor = 2.04764 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 300 -182631.38 0 -180992.85 28749.175 11000000
1000 1674.1317 -223215.2 0 -214071.49 -43.559921 11573511
2000 1443.4808 -225168.02 0 -217284.07 3.7486431 11231479
3000 1270.692 -226736.2 0 -219795.98 0.71739626 10904871

.
.
.
.

Total # of neighbors = 806268
Ave neighs/atom = 19.081
Neighbor list builds = 1035
Dangerous builds = 22
System init for write_restart …
Setting up run …
Memory usage per processor = 4.21604 Mbytes
Step Temp E_pair TotEng 3 tdiff
50000 299.66342 -234816.49 -233179.8 0 0
51000 298.48589 -234809.4 -233179.14 26.991906 132.34869
52000 300.24937 -234818.29 -233178.4 48.134262 216.23663
53000 300.35156 -234817.94 -233177.49 67.511275 253.24118
54000 300.65795 -234818.65 -233176.53 86.98012 273.35123
55000 300.92432 -234819.46 -233175.88 104.72854 278.42473
56000 299.14405 -234809.44 -233175.59 122.06699 284.53212
57000 299.67826 -234811.71 -233174.94 139.3028 300.54634

.
.
.
.

Total # of neighbors = 806524
Ave neighs/atom = 19.0871
Neighbor list builds = 550
Dangerous builds = 0
Setting up run …
Memory usage per processor = 4.33457 Mbytes
Step Temp E_pair TotEng 3 tdiff ave K
100000 306.93395 -234827.59 -233151.19 0 460.32424 460.32424 0
101000 307.6306 -234831.05 -233150.85 9.3210409 459.66693 459.99558 0.43398195
102000 307.84754 -234831.72 -233150.33 18.620655 441.17254 453.72124 0.87895524
103000 307.32652 -234828.09 -233149.55 28.463156 446.35555 451.87982 1.3490282
104000 308.20542 -234831.94 -233148.59 38.011678 446.43801 450.79145 1.8059357
105000 308.46388 -234832.33 -233147.57 47.386046 467.22107 453.52972 2.2377194
106000 308.23792 -234830.44 -233146.92 56.758293 460.10769 454.46943 2.6747646

.
.
.
.
.

193000 323.04911 -234840.79 -233076.37 835.50758 507.25976 474.51753 37.710218
194000 322.08415 -234834.39 -233075.24 844.18327 502.84934 474.81576 38.077859
195000 322.99552 -234838.66 -233074.54 853.05529 470.6863 474.77275 38.481529
196000 322.15434 -234833.62 -233074.09 862.24385 482.1616 474.84892 38.889787
197000 323.03706 -234837.66 -233073.31 870.78038 479.84933 474.89995 39.270591
198000 322.63209 -234834.51 -233072.37 879.81093 488.43577 475.03667 39.666431
199000 322.23855 -234831.55 -233071.56 888.75904 489.39559 475.18026 40.05775
200000 321.5094 -234826.35 -233070.34 897.6081 479.46705 475.2227 40.452978
Loop time of 513.501 on 64 procs for 100000 steps with 42255 atoms

Pair time () = 233.329 (45.4388) Neigh time () = 0.946099 (0.184245)
Comm time () = 175.881 (34.2514) Outpt time () = 0.464756 (0.0905072)
Other time (%) = 102.88 (20.035)

Nlocal: 660.234 ave 1333 max 0 min
Histogram: 12 6 5 5 4 3 7 6 5 11
Nghost: 646.625 ave 1287 max 22 min
Histogram: 2 1 1 10 24 11 8 5 1 1
Neighs: 0 ave 0 max 0 min
Histogram: 64 0 0 0 0 0 0 0 0 0
FullNghs: 12612.3 ave 25399 max 0 min
Histogram: 12 7 4 5 4 2 8 5 7 10

Total # of neighbors = 807186
Ave neighs/atom = 19.1027
Neighbor list builds = 1273
Dangerous builds = 0

Looks like you’ve used the KAPPA/in.mp script as a starting point. When performing the sampling (4th run of your calculation), your system is clearly not at equilibrium. You can see that since the total energy drifts in the NVE ensemble. Also, your MP run is not long enough. Have you checked if the temperature gradient and the heat flux density are converged? You need to have solid convergence criteria to look at, both on the temperature gradient (eg difference with step n-1), and on the heat flux (exchanged energy convergence, you can also look at the slope of the accumulated exchanged energy vs time). I would expect much longer simulation time, several hundreds of ps for the equilibration, and� a couple of ns for the actual NEMD run.

There is also something weird with your tdiff variable, it’s pretty high (~450K), when the mean temperature of your system is around 300K. Usually, in NEMD simulations the temperature difference is about 60K - which is by the way already extremely high for the small samples we study with molecular dynamics, when compared to experiments. There is probably something wrong with the way you’ve set up tdiff - don’t forget that KAPPA/in.mp uses LJ units.

Arthur