Hi all,
I am trying to calculate thermal conductivity using Green-Kubo based on LAMMPS input script for thermal conductivity of solid Ar. I have obtained the J0Jt.dat file. How can I obtain the curve for HCACF (heat current autocorrelation function)?

Hello all I have a query on a similar problem. I am trying to calculate the thermal conductivity of water using Green-Kubo method. The thermal conductivity value that I get is too small (order of 10^-8). I am giving the script that I have prepared by modifying the script that was available for calculating ‘kappa’ of Lennard Jones liquid.

#LAMMPS input file for calculating thermal conductivity of water
units real
atom_style full
boundary p p p
pair_style lj/cut/coul/long 10.0
pair_modify mix arithmetic
pair_modify tail yes
kspace_style pppm 1.0e-4
dielectric 1.0
special_bonds amber
bond_style harmonic
angle_style harmonic
dihedral_style none
improper_style none
read_data data.cos.1000SPCE
variable T equal 300
#variable P equal 1.0
velocity all create ${T} 12345 mom yes rot yes dist gaussian
timestep 0.1
variable dt equal 0.1
variable p equal 200 # correlation length
variable s equal 10 # sample interval
variable d equal $p*$s # dump interval
# Constraint ##################################
fix com all momentum 100 linear 1 1 1
fix rigid all shake 1e-4 20 0 b 1 a 1
# 1st equilibration run
fix 1 all nvt temp $T $T 0.5
thermo 100
run 1000
velocity all scale $T
unfix 1
# thermal conductivity calculation
reset_timestep 0
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 1 all nve
fix JJ all ave/correlate $s $p $d &
c_flux[1] c_flux[2] c_flux[3] type auto &
file profile.heatflux ave running
variable scale equal $s*dt/$T/$T/vol
variable k11 equal trap(f_JJ[3])*${scale}
variable k22 equal trap(f_JJ[4])*${scale}
variable k33 equal trap(f_JJ[5])*${scale}
variable kappa equal (v_k11+v_k22+v_k33)/3.0
thermo $d
thermo_style custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33 v_kappa
thermo_modify colname v_Jx Jx colname v_Jy Jy colname v_Jz Jz &
colname v_k11 kappa_11 colname v_k22 kappa_22 &
colname v_k33 kappa_33 colname v_kappa kappa
run 100000
print "Running average thermal conductivity: $(v_kappa:%.2f)"

I understand that the calculation of kinetic energy is done per atom as per this script and in case of water I might have to calculate it per molecule. As a complete beginner in LAMMPS I am unaware of any alternate ways in which this can be done. Please suggest any modifications that can be done to the input script.

There are a few obvious points where you have made rather questionable choices:

your timestep is 0.1fs but when a rigid water molecule using fix shake this is wasteful since you could be using 2fs and thus cover much more ground with the same computational effort.

your equilibration time is crazy short (1000 steps of 0.1 fs is only 0.1ps) and also your production simulation time is rather short. Unlike for a Lennard-Jonesium, water with its complex structure and much stronger interactions will take a much longer time to equilibrate and properly sample the relevant phase space

your thermostat constant is too short. please see the fix nvt documentation.

you should not be using fix momentum during production

Beyond that I can only give some general advice:

you need to make certain that your results are statistically converged, i.e. that you are neither suffering form significant finite size effects (due to a too small system) or that your accumulated data is noisy (due to too short an accumulation time and insufficient averaging)

since there are multiple methods to compute thermal conductivity, you can try multiple of those and see how well they agree with each other

you should check the published literature about limitations of the available methods, e.g. some cannot be applied to molecules with bonded interactions, other have problems when used with fix shake. This could apply here. I vaguely recall that compute centroid/stress/atom was implemented as alternative, but I am not certain about applying it to water

make sure you compare your results to those obtained for the same water potential. It is not likely that you will get an exact match to experimental data. Water is particularly difficult to model with LJ+Coulomb and 3 sites only.