Dear LAMMPS community,
I’ve followed the compute_heat_flux doc pages to generate the thermal conductivity of water using SPC/E water models. For some reasons, it estimates at 16 times more than the experimental value. I’ve come across some previous LAMMPS threads, e.g. http://lammps.sandia.gov/threads/msg15723.html , http://lammps.sandia.gov/threads/msg35433.html ,
http://lammps.sandia.gov/threads/msg30882.html , where similar problems occur.
AFAIK, GK formalism, in general, tends to yield a bit (<20%) over-estimation, but certainly not 16 times. There are papers that have roughly the same settings but produces much better results than mine. My question is: are my settings not rigorous enough, e.g time-step too large? SHAKE settings not rigorous enough? pppm not rigorous enough? correlation interval too short? I reduced the correlation length from 200 to 150 though.
At first the simulation cell is fixed with Langevin thermostat to equilibrate temperature, then the thermostat is unplugged and the simulation runs on NVE ensemble. Perhaps Langevin thermostat is not recommended and NVT should be use instead?
units real
atom_style full
pair_style hybrid lj/cut/coul/long 9.0
bond_style hybrid harmonic
angle_style hybrid harmonic
kspace_style pppm 0.0001
pair_modify mix geometric
read_data system.data
-------------------------------------------------------------------------------------
bond_coeff 1 harmonic 1000.0 1.0
angle_coeff 1 harmonic 1000.0 109.47
pair_coeff 1 1 lj/cut/coul/long 0.1553 3.16555529
pair_coeff 2 2 lj/cut/coul/long 0.0 2.058
group spce type 1 2
fix fShakeSPCE spce shake 0.0001 10 100 b 1 a 1
variable density equal mass(all)10/(6.0221417930vol) # density = [g/cm^3]
print “simulation density = ${density} g/cm^3” # experimental value ~ 998-1000 g/cm^3
print “experimental density = 1.000 g/cm^3”
----------------------------- fix nve + langevin -------------------------------
– minimization protocol –
print " equil via NVE + Langevin at 1K "
neighbor 2.0 bin
neigh_modify delay 0 every 1 check no cluster yes one 10000
fix 1 all nve
fix 2 all langevin 1.0 1.0 100.0 456483 # large damp factor -> small disspative friction
timestep 0.01
thermo 100
thermo_style custom step density vol temp press epair ebond eangle edihed pe etotal
thermo_modify norm no flush yes
dump 3 all custom 50 test.lammpstrj id mol type x y z ix iy iz
run 40000
write_restart system_after_nveLangevin.rst
----------------------------- refix langevin at 300K -------------------------------
print " equil via NVE + Langevin at 300K "
unfix 2
timestep 1.0
fix 2 all langevin 300.0 300.0 100.0 904124 # large damp factor -> small disspative friction
thermo 50
thermo_style custom step density vol temp press epair ebond eangle edihed pe etotal
thermo_modify norm no flush yes
run 50000
----------------------------- unfix langevin -------------------------------
print " equil via NVE "
params for Green-Kubo formalism
variable p equal 150 # correlation length
variable s equal 10 # sample interval
variable d equal {p}*{s} # dump interval
simulation settings
variable dt equal 1.0 # timestep
variable V equal vol
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}
variable T equal 300 # K
reset_timestep 0
timestep ${dt}
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 ID group-ID ave/correlate Nevery Nrepeat Nfreq
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
unfix 2
timestep 1.0
dump 2 all custom ${d} water_nve.lammpstrj id mol type x y z ix iy iz
thermo {d}
thermo_style custom step density vol temp press epair ebond eangle edihed pe etotal v_Jx v_Jy v_Jz v_k11 v_k22 v_k33
thermo_modify norm no flush yes
variable k equal (v_k11+v_k22+v_k33)/3.0
fix k all print {d} “${k}” file thermCond.txt screen no
run 6000000
variable ndens equal count(all)/vol
print “average conductivity: $k[W/mK] @ T K, {ndens} /A^3”
Many thanks,
Anh