# Lammps thermal conductivity

Hello people!
I am trying to calculate the thermal conductivity of phosphorene. I used the Green kubo methode with a tersoff potential. I have a problem with the simulation code because I am not able to achieve the desired value.
if you have some remarks to improve my simulation I ll be thankful.
#################################################################################
units real
variable T equal 300
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
#region box block 0 24 0 16 0 3
#create_box 1 box
#create_atoms 1 box
pair_style tersoff
pair_coeff * * P.tersoff P
timestep \${dt}
thermo \$d

# equilibration and thermalization

velocity all create \$T 102486 mom yes rot yes dist gaussian
fix NVT all nvt temp \$T \$T 10 drag 0.2
#run 8000
run 1000

# 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 NULL 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}
dump dump_1 all custom 100 yt.dump id type x y z ix iy iz vx vy vz
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 k equal v_k11
variable o equal v_k22
variable ndens equal count(all)/vol
print â€śAverage thermal conductivity armchair: \$k[W/mK] @ T K, {ndens} /A^3â€ť
print â€śAverage thermal conductivity ZigZag: \$o[W/mK] @ T K, {ndens} /A^3â€ť

1 Like

`````````
Like this,
```
``````

so that the special characters do not get misformatted.

Also, we do not know what you mean when you say your results are wrong. What do you expect to see? How do you know you should expect that? And what do you see instead?

I meant that I have found results that differ from the values presented in the lectures.
It means that I got wrong result.
I you have remarks on my input, because I think something is wrong in there. just I want to know that the code simulation is ok and then I will look if I can modify the tersoff parameters. Thanks!

``````units       real
variable    T equal 300
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
#region       box block 0 24 0 16 0 3
#create_box   1 box
#create_atoms 1 box
pair_style   tersoff
pair_coeff   * * P.tersoff P
timestep     \${dt}
thermo       \$d

# equilibration and thermalization

velocity     all create \$T 102486 mom yes rot yes dist gaussian
fix          NVT all nvt temp \$T \$T 10 drag 0.2
#run         8000
run          1000

# 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 NULL 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}
dump dump_1 all custom 100 yt.dump id type x y z ix iy iz vx vy vz
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     k equal v_k11
variable     o equal v_k22
variable     ndens equal count(all)/vol
print        "Average thermal conductivity armchair: \$k[W/mK] @ \$T K, \${ndens} /A\^3"
print        "Average thermal conductivity ZigZag: \$o[W/mK] @ \$T K, \${ndens} /A\^3"
``````

Did the lecture you referenced mentions how it gets the thermal conductivity values? There are a handful of factors that can affect the results of thermal conductivity. To list a few I can think of:

1. experimental factors (defects, impure samples, â€¦)
2. the theoretical method (Green-Kubo vs. NEMD vs. BTE â€¦)
3. the force field
4. the length of simulation for sampling heat flux, as well as the correlation length for ave/correlate (both should be long enough for convergence)
5. statistical errors (you may run the script with different random seed for initial velocities to check that)
6. NVE vs. NVT (remove comment for that â€śswitch to NVEâ€ť part to check)
7. the method of integration (sometimes the correlation function fluctuates heavily and itâ€™s tricky to do the integration correctly)
8. the size of system simulated (you may need a sufficiently large supercell to get the result converged)
9. method of calculating atomic stress (`stress/atom` vs. `centroid/stress/atom`; you probably have no choice here, but it can be a source of error since you have manybody forces)
1 Like

You should plot the correlation that you are integrating, and make sure that it actually decays to zero. Check the examples here: lammps-transport/thermal-conductivity/lj at main Â· rohskopf/lammps-transport Â· GitHub

Plot the heat flux autocorrelation - does yours look like this?

You can also verify that your average heat flux is zero.

1 Like

Thank you!

I took a (5-5) phosphorene sheet as a starting example. It is clear in the code that I used Green-kubo method with a tersoff potential (if someone has a reference of a useful tersoff parameters) . I got very small value of the thermal conductivity

almost this is how it looks!

Something I notice from the picture you posted:

1. the authors used NVE for calculating heat current, while yours use NVT;
2. the authors used much more steps for both equilibration and NVE runs than you;
3. the smallest cell in the figure is 20x20 UCs(=unit cell?).

btw: I believe youâ€™re aware of this already, but just in case, make sure that youâ€™re calculating in the same directions as the cited paper did. Sometimes external tools mess up the definition of axes when exporting lammps data files (Iâ€™ve personally encountered this).

the authors run both NVT and then NVE to calculate the heat current. I did the same I fixed to NVT for 10000 for equilibration and thermalization and then I fixed to NVE for 1000000 to calculate the heat current. But I still got a very weak value (~0.01â€¦W/mk).

Your input card looks good for me.
Did you try with a bigger sheet? I reported the simulations of graphene sheets with a similar protocol than this that you are using in the SI of
Can graphene improve the thermal conductivity of copper nanofluids? - Physical Chemistry Chemical Physics (RSC Publishing).
and I remember the TC changed when the sheet increased until a certain size.
Maybe it could be helpful to you.
OlguĂ­n-Orellana et al. - 2023 - Can graphene improve the thermal conductivity of c.pdf (2.3 MB)

1 Like

You worked with Lennard-jones potential, in my case I am trying to use tersoff potential to find the thermal conductivity of phosphorene. If the input file is fine then it should be something wrong with defining tersoff parameters for phosphorene?

How large, exactly, is your simulation cell? (Alternatively post the first few â€śdata blocksâ€ť of your LAMMPS data file.)

This is my data file, it contains (10Ă—10) UCs (400 atoms).

With (1) such a small z-direction (2) periodicity in the z-direction, it is very likely that your sheet is interacting with z-periodic copies of itself. I would recommend trying a z-height of 60 angstroms.

We increased the size of the sheet to (20-20) UCs (1600 atoms) also as you recommend we adjusted the z-height to 60 angstroms. The final result is still not good (0.20 W/mk, 0.23 W/mk for armmchair and zigzag resperctively)

Hello brother.