Thermal Conductivity of CdSe using Muller Plathe algorithm

I am trying to calculate the Thermal Conductivity (TC) of CdSe using Muller plathe algorithm.
The result I am getting is around 2.5~3.5 W/m-K, but from the literatures, the value should be around 9~11 W/m-K depending on the structure, I tried a lot of time and getting the same result, can you please help finding out the problem?

following is the lammps code I am using to determine the TC.

Muller-Plathe method via fix thermal_conductivity


variable x equal 1998.0894
variable y equal 43.9831
variable z equal 43.9831
variable t equal 300
variable first_run equal 100000

setup problem

units metal
dimension 3
boundary p p p
atom_style atomic


min_style cg
minimize 1e-3 1e-5 1000 10000

velocity all create $t 12354

pair_style sw
pair_coeff * *
comm_modify cutoff 14.70

1st equilibration run

fix 1 all nvt temp $t $t 0.1

reset_timestep 0
thermo 1000

run ${first_run}

unfix 1

2nd equilibration run

compute ke all ke/atom
variable temp atom c_ke/(1.5x8.621738x0.00001)

fix 1 all nve

thermo 10000

compute layers all chunk/atom bin/1d x lower 20 units box #reduced 0.025
fix 2 all ave/chunk 10 1000 10000 layers v_temp file
fix 3 all thermal/conductivity 500 x 100

variable tdiff equal f_2[51][3]-f_2[1][3]

variable hf equal f_3/(2x43.9831x43.9831)/(v_tdiff+0.0000000001)

thermo_style custom step temp epair etotal f_3 v_tdiff v_hf
thermo_modify colname f_3 E_delta colname v_tdiff dTemp_step colname v_hf hf_step

variable start_time equal time

fix print1 all print 10000 “${hf}” file hf.txt screen no

run 2000000

thermal conductivity calculation

fix ave all ave/time 1 1 10000 v_tdiff ave running

variable kappa equal (f_3/(({start_time}-{first_run}/1000)($y$z)2.0))(($x/2.0)/f_ave)*1602

thermo_style custom step temp epair etotal f_3 v_tdiff f_ave v_hf v_start_time v_kappa
thermo_modify colname f_3 E_delta colname v_tdiff dTemp_step colname f_ave dTemp colname v_hf hf_step

run 100000
print “Running average thermal conductivity: $(v_kappa:%.2f)”

which paper?

in this paper, they used same potential

But the abstract says 14-18 W/m.K, and yet your initial post said 9-11 W/m.K …

What I’m trying to get at is that it is impossible for us to know why a LAMMPS simulation isn’t fulfilling your expectations (or even whether your expectations are right or wrong) unless we know where your expectations come from. There are a lot of settings that go into an MD run and any of them, from force field to timestep to fixes to thermostat / barostat parameters, could be what makes the difference. Every year thousands of papers are published using water models which are pretty good at 30C, but which don’t freeze at 0C or boil at 100C (and that’s not a problem, since the simulations aren’t being run at that temperature), so if even plain water doesn’t behave the way it “should”, why should your simulation?

It is much easier to start investigating whether a simulation is fulfilling expectations if the expectations come from:

  1. Rigorous mathematical physics. If a system should conserve momentum or energy, but doesn’t, then that’s clearly wrong and will immediately arouse lots of interest. If a water model freezes when you heat it and boils when you pressurize it, that’s clearly wrong too. But if some number is X when it’s “supposed” to be 1.5*X, it’s exceedingly difficult to know why unless you happen to come across an expert in that particular subfield.

  2. A prior publication. But if you claim the result comes from that paper, you must first show that you have done your best to completely replicate the paper. Again, any change in parameter can lead to a change in observables, including box sizes and configuration geometry and so on.

So without knowing that you have tried to very closely replicate this paper, it’s very hard to know if your results are “as expected” or not.