Dear Mahesh,

I have a LAMMPS script that applies the Green-Kubo relation corresponding to the thermal conductivity (i.e. based on the self-correlation function of the internal heat). Note, however, that this approach underpredicts the value of the thermal conductivity in cases where electronic conduction is important (as in metals). I can pass along the script if you are interested.

Dear Jaime,

Thanks for replying. Objective of my study is to predict the value of thermal conductivity as close as experimental or literature value by modifying/changing the potentials to accommodate relevant interactions . It would be of great help if you could mail the input script. 

mail id: [email protected], [email protected]

Thanks in advance

Hi Mahesh,

I realize that I never sent you this. Sorry for the delay. It is for Ge so you have to adapt it to what you want. It is based on the Green-Kubo formulas for the time self-correlation function of the heat. Here it goes:

# Sample LAMMPS input script for thermal conductivity of solid Ge

units real

atom_style atomic

variable T equal 700

variable V equal vol

variable dt equal 4.0

variable p equal 200 # 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

dimension 3

boundary p p p

lattice diamond 5.654 &

orient x 1 0 0 &

orient y 0 1 0 &

orient z 0 0 1

region box block 0 2 0 2 0 2

create_box 1 box

create_atoms 1 box

read_dump amorph.dump 220000 x y z replace no add yes purge yes format native

###############

mass * 72.64

pair_style sw

pair_coeff * * Ge.sw Ge

timestep ${dt}

thermo $d

# equilibration and thermalization

velocity all create $T 329001 mom yes rot yes dist gaussian

fix NVT all nvt temp $T $T 10 drag 0.2

run 8000

# thermal conductivity calculation, switch to NVE if desired

#unfix NVT

#fix NVE all nve

reset_timestep 0

compute myKE all ke/atom

compute myPE all pe/atom

compute myStress all stress/atom 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

fix JJ all ave/correlate $s $p $d &

c_flux[1] c_flux[2] c_flux[3] type auto file J0Jt.dat ave running

variable scale equal ${convert}/${kB}/$T/$T/$V*$s*${dt}

variable k11 equal trap(f_JJ[3])*${scale}

variable k22 equal trap(f_JJ[4])*${scale}

variable k33 equal trap(f_JJ[5])*${scale}

thermo_style custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33

run 100000

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"

Thanks Jaime for script. It seems very useful for my work.

Hi Jaime,

cannot open file amorf.dump

Navid,

the lammps script was intended for calculating the thermal conductivity in amorphous Ge. "amorf.dump" is an input file with the coordinates of a block of amorphous Ge. You have to adapt it to your purposes. Change the Ge potential to whatever material you want and create your own input configuration.

Thanks  Jaime

I'm working on thermal conductivity of graphene. I read this paper entitled “Strain engineering of thermal conductivity in graphene sheets and nanoribbons: a demonstration of magic flexibility ". In figure 1-d of this paper the temperature of hot region value is larger than the other regions ، and using the airebo potential, But I use Tersoff potential for interactions. I have a question. Why the temperature of hot region is less than some other regions? Is it because of potential type and lattice, or Fixing mistakes or …?

The sample of my GNR is divided into 50 slabs along X. the temperature of hot region value is the smallest of Some parts in my logfile. Please let me know what is wrong in my calculation and refer to my code posted here:

units real

dimension 2

atom_style atomic

boundary p p p

#Atom definition

lattice custom 2.53 a1 1 0 0 a2 0 1.732051 0 a3 0 0 1.361789 &

basis 0 0 0 &

basis 0.5 0.5 0 &

basis 0.5 0.16666 0 &

basis 0 0.66666 0

neighbor 2.0 bin

neigh_modify every 3 delay 3

#region

region graphene block 0 100 0 50 0 3 units box

region cold1 block 0 1.5 0 50 0 3 units box

region r2 block 1.5 3.5 0 50 0 3 units box

region r3 block 3.5 5.5 0 50 0 3 units box

region r4 block 5.5 7.5 0 50 0 3 units box

region r5 block 7.5 9.5 0 50 0 3 units box

region r6 block 9.5 11.5 0 50 0 3 units box

region r7 block 11.5 13.5 0 50 0 3 units box

region r8 block 13.5 15.5 0 50 0 3 units box

region r9 block 15.5 17.5 0 50 0 3 units box

region r10 block 17.5 19.5 0 50 0 3 units box

region r11 block 19.5 21.5 0 50 0 3 units box

region r12 block 21.5 23.5 0 50 0 3 units box

region r13 block 23.5 25.5 0 50 0 3 units box

region r14 block 25.5 27.5 0 50 0 3 units box

region r15 block 27.5 29.5 0 50 0 3 units box

region r16 block 29.5 31.5 0 50 0 3 units box

region r17 block 31.5 33.5 0 50 0 3 units box

region r18 block 33.5 35.5 0 50 0 3 units box

region r19 block 35.5 37.5 0 50 0 3 units box

region r20 block 37.5 39.5 0 50 0 3 units box

region r21 block 39.5 41.5 0 50 0 3 units box

region r22 block 41.5 43.5 0 50 0 3 units box

region r23 block 43.5 45.5 0 50 0 3 units box

region r24 block 45.5 47.5 0 50 0 3 units box

region r25 block 47.5 49.5 0 50 0 3 units box

region hot block 49.5 52.5 0 50 0 3 units box

region r27 block 52.5 54.5 0 50 0 3 units box

region r28 block 54.5 56.5 0 50 0 3 units box

region r29 block 56.5 58.5 0 50 0 3 units box

region r30 block 58.5 60.5 0 50 0 3 units box

region r31 block 60.5 62.5 0 50 0 3 units box

region r32 block 62.5 64.5 0 50 0 3 units box

region r33 block 64.5 66.5 0 50 0 3 units box

region r34 block 66.5 68.5 0 50 0 3 units box

region r35 block 68.5 70.5 0 50 0 3 units box

region r36 block 70.5 72.5 0 50 0 3 units box

region r37 block 72.5 74.5 0 50 0 3 units box

region r38 block 74.5 76.5 0 50 0 3 units box

region r39 block 76.5 78.5 0 50 0 3 units box

region r40 block 78.5 80.5 0 50 0 3 units box

region r41 block 80.5 82.5 0 50 0 3 units box

region r42 block 82.5 84.5 0 50 0 3 units box

region r43 block 84.5 86.5 0 50 0 3 units box

region r44 block 86.5 88.5 0 50 0 3 units box

region r45 block 88.5 90.5 0 50 0 3 units box

region r46 block 90.5 92.5 0 50 0 3 units box

region r47 block 92.5 94.5 0 50 0 3 units box

region r48 block 94.5 96.5 0 50 0 3 units box

region r49 block 96.5 98.5 0 50 0 3 units box

region cold2 block 98.5 100 0 50 0 3 units box

#create atoms&box

create_box 1 graphene

create_atoms 1 region graphene

#group

group cold1 region cold1

#group r2 region r2

group r3 region r3

#group r4 region r4

group r5 region r5

#group r6 region r6

group r7 region r7

#group r8 region r8

group r9 region r9

#group r10 region r10

group r11 region r11

#group r12 region r12

group r13 region r13

#group r14 region r14

group r15 region r15

#group r16 region r16

group r17 region r17

#group r18 region r18

group r19 region r19

#group r20 region r20

group r21 region r21

#group r22 region r22

group r23 region r23

group r24 region r24

group r25 region r25

group hot region hot

group r27 region r27

group r28 region r28

group r29 region r29

#group r30 region r30

#group r31 region r31

#group r32 region r32

group r33 region r33

#group r34 region r34

group r35 region r35

#group r36 region r36

group r37 region r37

#group r38 region r38

group r39 region r39

#group r40 region r40

group r41 region r41

#group r42 region r42

group r43 region r43

#group r44 region r44

group r45 region r45

#group r46 region r46

group r47 region r47

group r48 region r48

group r49 region r49

group cold2 region cold2

group htot subtract all hot cold1 cold2

group graphene region graphene

#Settings and Simulation

pair_style tersoff

pair_coeff * * sic.tersoff C

mass * 12.011

###EQUILIBRATION 1###

fix 25 all nvt temp 300.0 300.0 100.0

thermo 2000

run 200000

###EQUILIBRATION 2###

fix 10 all nve

dump 1 all xyz 10 gf.xyz

run 10000

unfix 10

###EQUILIBRATION 3###

fix 1 htot nve

fix 3 cold1 heat 10 -0.05

fix 2 hot heat 10 0.1

fix 33 cold2 heat 10 -0.05

#run 1000

#output1

thermo 2000

thermo_style custom step temp ke etotal

timestep 0.5

run 20000

fix 28 all npt temp 300.0 300.0 100.0 iso 0.0 0.0 1000.0

thermo 2000

run 200000

#--compute--#

compute coldBath1 cold1 temp

compute coldBath2 cold2 temp

compute hotBath hot temp

compute T3 r3 temp

compute T5 r5 temp

compute T7 r7 temp

compute T9 r9 temp

compute T11 r11 temp

compute T13 r13 temp

compute T15 r15 temp

compute T17 r17 temp

compute T19 r19 temp

compute T21 r21 temp

compute T23 r23 temp

compute T24 r24 temp

compute T25 r25 temp

compute T27 r27 temp

compute T28 r28 temp

compute T29 r29 temp

compute T33 r33 temp

compute T35 r35 temp

compute T37 r37 temp

compute T39 r39 temp

compute T41 r41 temp

compute T43 r43 temp

compute T45 r45 temp

compute T47 r47 temp

compute T48 r48 temp

compute T49 r49 temp

#output2

thermo 2000

thermo_style custom step temp c_coldBath1 c_T3 c_T5 c_T7 c_T9 c_T11 c_T13 c_T15 c_T17 c_T19 c_T21 c_T23 c_T24 c_T25 c_hotBath c_T27 c_T28 c_T29 c_T33 c_T35 c_T37 c_T39 c_T41 c_T43 c_T45 c_T47 c_T48 c_T49 c_coldBath2

run 28000

This is my log file

Step Temp coldBath T3 T5 T7 T9 T11 T13 T15 T17 T19 T21 T23 T24 T25 hotBath T27 T28 T29 T33 T35 T37 T39 T41 T43 T45 T47 T48 T49 coldBath

458000 337.63682 165.08407 322.01872 305.54054 328.61349 215.70579 404.77457 414.60319 399.77107 400.23761 270.63301 360.45658 347.42522 437.33378 395.87029 244.82007 354.8689 358.9156 333.21225 341.92651 507.20771 262.9455 365.25513 334.81607 315.08318 303.57287 323.96405 325.82107 334.73165 144.66116

Thank you very much.

Best Regards,

navid

Dear Jamie, I have also used something similar to this. Do you have any published reference for this? Thanks

Alex

Hello friends

i saw all  of you are using Ge method to calculate thermal conductivity, while ,i want to use MP method, so did anyone try this before? i'd like to discuss with someone.

ps: i know there has one mp example of lammps, but i want to try to do one for graphene

Thank you  

Anyone can help me out to calculate thermal conductivity for liquid Aluminum using Green-Kubo method??

I tried with change unit from real to metal in input file given in manual and the value for conductivity is 100 times smaller.. Do anyone have idea about it..? Please help me out.

compute myStress all stress/atom virial is giving an error:- Could not find compute stress/atom temperature ID (../compute_stress_atom.cpp:61)

What can be the problem?

Dear Sandeep

Your answer can be found in lammps manual for compute stress. You are going to need a compute command for the temperature and call its ID right after "stress/atom" or just simply put "NULL" to make it work.