[lammps-users] A question about thermal conductivity calculation using compute heat/flux command

Dear All,

I recently wanted to calculate the thermal conductivity of eutectic of
Li2CO3 and K2CO3. I used the sample script in the compute heat/flux
command page. What I found from the result is the value of the thermal
conductivity was negative. It was -0.33 W/mK.
(The literature value is 0.5 W/mK)

Anyone has an idea why I got the negative value?

Does it mean my simulation went totally wrong?

Please give me some advice.

The following is my input script for your reference.

Thank you for your kind reply in advance.

Dong,

## MINIMIZATION

units real
dimension 3
boundary p p p
atom_style full
pair_style buck/coul/long 15.0
pair_modify shift yes
kspace_style ewald 1e-4
bond_style harmonic
angle_style harmonic
improper_style cvff

read_data CAR_3600.lammps05
pair_coeff 1 1 31602.558709 0.289855072 0
pair_coeff 1 2 31119.783663 0.289855072 0
pair_coeff 1 3 9108.290566 0.289855072 0
pair_coeff 1 4 55624.468654 0.289855072 0
pair_coeff 2 2 21684.477780 0.289855072 0
pair_coeff 2 3 6556.870769 0.289855072 0
pair_coeff 2 4 45853.882966 0.289855072 0
pair_coeff 3 3 1975.676210 0.289855072 0
pair_coeff 3 4 13629.974296 0.289855072 0
pair_coeff 4 4 89024.097848 0.289855072 0

reset_timestep 0
timestep 0.5
neigh_modify delay 5 every 1 one 10000
dump dump1 all atom 1000 dump.min.Car
velocity all create 0.0 2705163 dist uniform
thermo 1000
min_style sd
min_modify dmax 0.1 line quadratic
minimize 0.0 0.0 200000 1000000
## Finish
write_restart restart.min.Car
print "Finished with minimization!"
clear

## NVE

units real
dimension 3
boundary p p p
atom_style full
pair_style buck/coul/long 15.0
pair_modify shift no tail yes
kspace_style ewald 1e-4
bond_style harmonic
angle_style harmonic
improper_style cvff

read_restart restart.min.Car
pair_coeff 1 1 31602.558709 0.289855072 0
pair_coeff 1 2 31119.783663 0.289855072 0
pair_coeff 1 3 9108.290566 0.289855072 0
pair_coeff 1 4 55624.468654 0.289855072 0
pair_coeff 2 2 21684.477780 0.289855072 0
pair_coeff 2 3 6556.870769 0.289855072 0
pair_coeff 2 4 45853.882966 0.289855072 0
pair_coeff 3 3 1975.676210 0.289855072 0
pair_coeff 3 4 13629.974296 0.289855072 0
pair_coeff 4 4 89024.097848 0.289855072 0

reset_timestep 0
timestep 0.5
neigh_modify delay 5 every 1 one 10000
dump dump2 all atom 1000 dump.NVE.Car
velocity all create 800.0 2705163 dist uniform
thermo 1000
fix 1 all nve
run 150000
write_restart restart.NVE.Car
print "Finished with NVE!"
clear

## NPT

units real
dimension 3
boundary p p p
atom_style full
pair_style buck/coul/long 15.0
pair_modify shift no tail yes
kspace_style ewald 1e-4
bond_style harmonic
angle_style harmonic
improper_style cvff

read_restart restart.NVE.Car
pair_coeff 1 1 31602.558709 0.289855072 0
pair_coeff 1 2 31119.783663 0.289855072 0
pair_coeff 1 3 9108.290566 0.289855072 0
pair_coeff 1 4 55624.468654 0.289855072 0
pair_coeff 2 2 21684.477780 0.289855072 0
pair_coeff 2 3 6556.870769 0.289855072 0
pair_coeff 2 4 45853.882966 0.289855072 0
pair_coeff 3 3 1975.676210 0.289855072 0
pair_coeff 3 4 13629.974296 0.289855072 0
pair_coeff 4 4 89024.097848 0.289855072 0

reset_timestep 0
timestep 0.5
neigh_modify delay 5 every 1 one 10000
dump dump3 all atom 1000 dump.NPT.Car
velocity all create 800.0 2705163 dist uniform
thermo 10

variable DENSITY equal 1.66053886*mass(all)/vol
fix f1 all ave/time 1 999 1000 c_thermo_temp v_DENSITY
file DENavg.Car
fix f2 all ave/time 1 1 10 c_thermo_temp v_DENSITY
file DENint.Car
fix f3 all npt temp 800.0 800.0 100.0 iso 1.0 1.0 200.0
compute c1 all pe
log thermo_pe_enthalpy.log
thermo_style custom step temp enthalpy c_c1 vol
run 150000
write_restart restart.NPT.Car
print "Finished with NPT!"
clear

## Thermal Conductivity

units real
dimension 3
boundary p p p
atom_style full
pair_style buck/coul/long 15.0
pair_modify shift no tail yes
kspace_style ewald 1e-4
bond_style harmonic
angle_style harmonic
improper_style cvff

read_restart restart.NPT.Car
pair_coeff 1 1 31602.558709 0.289855072 0
pair_coeff 1 2 31119.783663 0.289855072 0
pair_coeff 1 3 9108.290566 0.289855072 0
pair_coeff 1 4 55624.468654 0.289855072 0
pair_coeff 2 2 21684.477780 0.289855072 0
pair_coeff 2 3 6556.870769 0.289855072 0
pair_coeff 2 4 45853.882966 0.289855072 0
pair_coeff 3 3 1975.676210 0.289855072 0
pair_coeff 3 4 13629.974296 0.289855072 0
pair_coeff 4 4 89024.097848 0.289855072 0

variable kB equal 1.3806504e-23 # [J/K] Boltzmann
variable kCal2J equal 4186.0/6.02214e23
variable T equal 800
variable V equal vol
variable dt equal 0.5
variable p equal 200 # correlation length
variable s equal 10 # sample interval
variable d equal $p*s \# dump interval \# \-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\- 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 50000
# -------------- flux calculation ---------------
reset_timestep 0
compute myKE all ke/atom
compute myPE all pe/atom
compute myStress all stress/atom 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 {kCal2J}*\{kCal2J\}/{kB}/$T/$T/$V*s\*{dt}*1.0e25
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
run 150000
variable k equal (v_k11+v_k22+v_k33)/3.0
write_restart restart.TC.Car
print "average conductivity: $k [W/mK]"

I don't know. Your script is long and complicated - no one
else is likely to debug it for you. You do reallize that
the long-range Coulombics are not tallied in the per-atom
energy and stress, as the doc pages for those commands
indicate? So they are not part of your heat flux.

Steve