Calculate thermal conductivity by using NEMD

LAMMPS LAMMPS Beginners
1

Dear LAMMPS Users,

I am using NEMD method to calculate thermal conductivity. My script is modified based on the “in.langevin” example file given by the LAMMPS. However, at the end of the simulation, the cumulative energy in heat source and heat sink is not equal. I tried to change the damp parameter in fix langevin and increase the total simulation steps, but the result was not improved. I would like to ask what parameters can affect the cumulative energy in such a simulation? The script I use after the minimization and relaxation and the log file are attached here. Thank you!

echo            both
units           metal
dimension       3
processors      16 2 2
boundary        p p p
atom_style      atomic

read_restart    sicgraph.restart

neighbor        2.0 bin
neigh_modify    delay 0 every 1 check yes

group           SiAtom    type 1
group           CAtom     type 2 3

mass            1 28.08
mass            2*3 12.01

pair_style      hybrid tersoff lj/cut 8.5
pair_coeff      * * tersoff SiC.tersoff Si C C C C
# pair_coeff      * * airebo CH.airebo NULL NULL C C C
pair_coeff      3 4 lj/cut 0.0024 3.4
pair_coeff      3 5 lj/cut 0.0024 3.4
pair_coeff      4 5 lj/cut 0.0024 3.4

thermo_style    custom step temp pe ke etotal press vol
thermo          1000
timestep        0.0003

fix             5 all nvt temp 300 300 0.03
run             40000
unfix           5
fix             6 all npt temp 300 300 0.03 iso 0 0 0.3
run             100000
unfix           6

reset_timestep  0

# define heat source and heat sink

region          hot block 95 105 INF INF INF INF units box
region          cold block 295 305 INF INF INF INF units box
compute         Thot all temp/region hot
compute         Tcold all temp/region cold

variable        thi equal 310
variable        tlo equal 290

fix             1 all nve
fix             hot all langevin ${thi} ${thi} 0.03 59804 tally yes
fix             cold all langevin ${tlo} ${tlo} 0.03 59804 tally yes
fix_modify      hot temp Thot
fix_modify      cold temp Tcold

variable        tdiff equal c_Thot-c_Tcold
thermo_style    custom step temp c_Thot c_Tcold f_hot f_cold v_tdiff
thermo	        1000
run             200000

reset_timestep  0

compute         myKE all ke/atom
variable        atemp atom c_myKE/(1.5*8.621738*0.00001)

fix             hot all langevin ${thi} ${thi} 0.03 59804 tally yes
fix             cold all langevin ${tlo} ${tlo} 0.03 59804 tally yes
fix_modify      hot temp Thot
fix_modify      cold temp Tcold

fix             ave all ave/time 10 100 1000 v_tdiff ave running

compute         cc1 all chunk/atom bin/1d x lower 0.025 units reduced
fix             4 all ave/chunk 1 50000 50000 cc1 v_atemp file Temperature.txt
dump            1 all atom 50000 dump.nve
thermo_style    custom step temp c_Thot c_Tcold f_hot f_cold v_tdiff f_ave

run             800000
write_restart   sicgraph.restart

log.lammps (117.4 KB)

2

I’d say your force field is inconsistent and thus bogus. Why use a hybrid style when all elements are covered by the tersoff potential? Also, the override with lj/cut is only partial and creating inconsistent processing of the tersoff interactions.

3

Thank you for your help. Yes, in my case, all the elements can be covered by tersoff potential. My consideration is that if I use both tersoff and lj potentials, I need to use hybrid style, or maybe there is a better way? And you mentioned the override with lj can bring issues to the tersoff interactions, does that mean I cannot combine tersoff with lj? What if I use airebo potential instead of tersoff on element 3 to 5? Thank you.

4

Why do you want to use lj/cut when you can use Tersoff??
Besides what you do in your input is bogus. Full stop.

You could, but not like this. When using a group of atom types with a manybody potential, all internal interactions must be covered by the manybody potential and no inter-group interactions may be described by a manybody potential but must be described by a pair-wise additive potential.

The outcome would still be questionable. Are you trying to follow some published description?

The same “rules” that I mentioned (and that are also explained in the documentation) must be followed.

Please note that using a hybrid potential is almost always inferior to a consistent description with the same functional form.

5

Thank you for the reply. I use graphene and SiC in my system, and I want to use lj to describe the interactions between different graphene layers.
I read some papers, usually what they do is use tersoff for SiC, airebo for graphene (or tersoff), and lj for the interactions between SiC and graphene, and between different graphene layers. But here I want to show strong interactions between SiC and graphene, so I tried to follow this configuration but only changing the interactions between SiC and graphene from lj to tersoff. Is that okay with hybrid style?

6

The hybrid style will run, but it is wrong scientifically. You treat Tersoff like a pairwise additive potential but it is not. I already mentioned the “rules” and your setup violates it and thus is bogus.

If you want to follow the publications, you must either use Tersoff twice in the pair style or Tersoff once and AIREBO once and apply those for all atom types within those compounds. You may use LJ only for the mixed terms between the two compounds.

If you want to use Tersoff for the inter-compound interactions you must use it for everything. Everything else is leading to inconsistent interactions and thus bogus.

7

Thank you, I understand it now. Really appreciate your help!

8

Mind you. This is all documented.

9

Okay, I will read the manual more carefully. Thank you!