Hello everyone!
I try to use the compute heat/flux command to calculate the thermal conductivity of the water-saturated clay system in LAMMPS. I expect the thermal conductivity value should be around 1 W/mK (based on the literature data), however, the simulation results turn out to be around 1,000 W/mK, which is not very likely to be true. The LAMMPS version I used is 4May2022. I tried several ways to debug:
-
Intuitively, this might be just a unit problem. But I have checked my script many times and it seems not to be an issue (if I don’t miss anything). I have attached my script below, which is mainly based on the LAMMPS manual (compute heat/flux command — LAMMPS documentation).
-
Then, I start to think this issue might be caused by the fact that LAMMPS believes that the H atoms in the water molecules have a high potential energy (because the 'special_bonds lj/coul 0 1 1 ’ command tells it to calculate the H-H Coulomb interaction, but the 'fix shake’ command tells it to keep the H-O-H angle fixed). So I applied 'the neigh_modify exclude molecule/intra water’ command to ignore the interactions within each water molecule. However, the results are still 1,000 times larger than expected.
I am wondering if there is anything else I should be careful about for this computation. Any comments will be highly appreciated. Thanks a lot!
units metal
dimension 3
atom_style full
pair_style lj/cut/coul/long 12.0
pair_modify mix arithmetic
bond_style harmonic
angle_style harmonic
dihedral_style opls
read_restart NaMMT_01.save #this reads the previous equilibrated system
reset_timestep 0
#LJ pair coefficients for all atom types (the units are eV and Angstrom)
pair_coeff 1 1 5.76646e-8 4.27124 #Al
pair_coeff 2 2 7.98102e-8 3.30203 #Si
pair_coeff 3 3 6.73853e-3 3.16556 #O
#…
pair_coeff 19 19 4.33633e-3 2.58400 #Na
pair_coeff 20 20 6.73853e-3 3.16556 #Ow
pair_coeff 21 21 0.0 0.00000 #Hw
bond_coeff 1 24.0291 1.0
bond_coeff 2 24.0291 1.0
angle_coeff 1 1.30090 109.47
angle_coeff 2 1.98472 109.47
special_bonds lj/coul 0 1 1 #Ignore LJ and Coulomb interactions between atoms connected through a bond
kspace_style pppm 0.01
run_style verlet
group water type 20 21
group solutes type 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
variable T equal 314.8 #K
variable V equal vol #Angstrom^3
variable dt equal 0.001 #ps
variable p equal 200 # correlation length
variable s equal 10 # sample interval
variable d equal $p*$s # dump interval
variable kB equal 1.3806504e-23 # Boltzmann constant, J/K
variable eV2J equal 1.60218e-19 #convert eV to J
variable A2m equal 1.0e-10 #convert A to m
variable ps2s equal 1.0e-12 #convert ps to s
variable convert equal {eV2J}*{eV2J}/{ps2s}/{A2m}
thermo $d
timestep ${dt} #in ps
compute myKE all ke/atom
compute myPE all pe/atom
#compute myStress all stress/atom NULL virial
compute myStress all centroid/stress/atom NULL pair bond angle
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 k11 equal trap(f_JJ[3])/${kB} /$T /$T /$V*s*{dt}*${convert}
variable k22 equal trap(f_JJ[4])/${kB}/$T/$T/$V*s*{dt}*${convert}
variable k33 equal trap(f_JJ[5])/${kB}/$T/$T/$V*s*{dt}*${convert}
thermo_style custom step temp press vol ke pe v_Jx v_Jy v_Jz v_k11 v_k22 v_k33
fix 1 water shake 0.0001 20 0 b 2 a 2 #water molecules are kept rigid using the SHAKE algorithm
neigh_modify exclude molecule/intra water #ignore the interactions within each water molecule
fix 5 solutes npt temp $T $T 1 z 1 1 10
run 500000
unfix 5
variable k equal (v_k11+v_k22+v_k33)/3.0
print “average conductivity: $k [W/mK] @ $T K”
write_restart NaMMT_02.save #writes a restart file at the end of the simulation
LAMMPS runs calculations, it does dispense truths, nor does it have beliefs. You are using a fairly complicated model for Al/Si/O/Na/O. It is probably good for certain types of calculations, but it may not be good for thermal conductivity in water-saturated clay. I suggest you take a couple of steps back and try to use LAMMPS to calculate the thermal conductivity for a system where previously published results are available. Once you can match those results, you will have a much better chance of executing a correct calculation for the above system.