About thermal conductivity calculation of Graphene

LAMMPS LAMMPS Mailing List Mirror
1

Dear All,

I am trying to calculate thermal conductivity of single layer graphene using Muller Plathe method. Attached is my input script as well as data file for the same.

I have an error while execution: “ERROR: Divide by 0 in variable formula (…/variable.cpp:1426)”

This is in the command line:
“variable thermal_conductivity equal f_heat_swap1.60217e-19/(20.00051000001.0e-125051.77f_delta_t1.0e-10)*51.03”

Detailed input as per attached.

If I put “v_delta_t” instead of “f_delta_t”, it works fine but value of thermal conductivity is negative in some cases.

Help in resolving the same would be great. Also, inputs on any error in the input file can be suggested.

Best Regards,
Sarang

graphene_1008atoms.txt (280 KB)

graphene_thermal_conductivity.in.in (2.22 KB)

2

I’m guessing your f_delta_t is 0.0 at the point at which
the variable is evaluated, maybe for the first time on step 0?

If so, the denominator of your formula is 0.0 and LAMMPS

complains.

You have to write a valid formula. E.g. maybe you
should use (f_delta_t + epsilon).

Steve

3

Thanks a lot for the suggestion Steve. I changed the line of variable by adding a constant, so as to exclude the first step where value is zero.

Another major issue is when I am removing NVT thermostat and shifting to NVE for performing thermal conductivity calculations, the temperature of the system shoots up to thousands of kelvin. I could not understand a way to stabilize it. Heat flux is quite linear.

My input script and data file are as attached for the single layer of graphene structure.

Please let me know a way to do the same.

Best Regards,
Sarang

graphene_1008atoms.txt (280 KB)

graphene_thermal_conductivity.txt (2.1 KB)

4

Try reducing your heat flux.

5

Thanks a lot. I actually tried changing the thermostat and my revised file is as per attached “graphene_thermal_conductivity.in” and also mentioned below. However, thermal conductivity value of graphene is found to be very small as compared to literature.

Temperature profile along Z direction is also very linear as per attached plot “Temp (K) Vs. Z(Ang)”. Also, attached is kinetic energy graph versus time (ps) as per file “Kinetic Energy Swapped versus time (ps)”, which also seems fine. In order to calculate thermal conductivity, I calculated slope of temperature profile as dT/dJ = 0.6293 K / Ang, while slope of Kinetic energy exchanged versus time is 11.11 eV / ps. Simulation box is orthogonal with dimensions = 50 x 50 x 100 (Ang x Ang x Ang) with 1008 C atoms. Hence, cross sectional area for heat flow along z direction is 50 x 50 = 2500 Ang^2.

Thus, using formula: Thermal conductivity k = (dJ/dt) / (Cross Section Area)(dT/dJ) = 11.11 eV/ps / (50*50) * 0.6293 K/Ang ~ 11.3 W/mK.

Please let me know what is error in my calculation.

Input Script:

units metal

boundary p p p

atom_style full

#Atom definition

read_data graphene_1008atoms

#Creating Regions

region up block INF INF INF INF -5 5 units box

region down block INF INF INF INF -50 -40 units box

#Potentials

pair_style tersoff

pair_coeff * * SiC.tersoff C

mass 1 12

neighbor 2.0 bin

neigh_modify delay 5

dump snap all cfg 10000 dump.config.*.cfg id type xs ys zs vx vy vz

equilibrium

velocity all create 300 458127641 mom yes rot yes dist gaussian

fix NPT all npt temp 300 300 1 iso 0 0 0.1 drag 1

thermo_style custom step temp etotal vol lx ly lz press pxx pyy pzz

thermo_modify lost warn

thermo 1000

timestep 0.001

run 50000

reset_timestep 0

unfix NPT

fix temp all temp/berendsen 300 300 0.1

fix nve all nve

variable kb equal 1.3806504e-23

variable ev2J equal 1.60217e-19

compute ke all ke/atom

variable temp atom c_ke*{ev2J}/(1.5*{kb})

fix temp_profile_whole all ave/spatial 1 100000 100000 z -50 50 v_temp file temp.profile_whole units box # for whole structure

#calculate temperature difference between hottest and coldest regions

compute hottest_temp all temp/region up

compute coldest_temp all temp/region down

variable delta_temp equal c_hottest_temp-c_coldest_temp

fix delta_out all ave/time 1 100000 100000 v_delta_temp file delta_temp_graphene.dat

fix hottesttemp_out all ave/time 1 100000 100000 c_hottest_temp file hottesttemp_graphene.dat

fix coldesttemp_out all ave/time 1 100000 100000 c_coldest_temp file coldesttemp_graphene.dat

run 100001

unfix temp

fix heat_swap all thermal/conductivity 10 y 80

fix e_exchange all ave/time 10 1000 100000 f_heat_swap file e_exchange_graphene.dat

run 2000000

Best Regards,
Sarang

delta_temp_graphene.dat (384 Bytes)

Temp(K) Vs. Z(Ang).jpg

temp.profile_whole.doc (1.73 KB)

graphene_thermal_conductivity.in (1.88 KB)

Kinetic Energy Swapped Vs. Time.jpg

6

Think again. The area should be 50×3.4

That will make the conductivity of 10 nm graphene to be ~ 150 W/mK which is very close to what many reports have published.

Cheers

7

Thanks a lot for the help and advice.

Best Regards,
Sarang