# Green-Kubo

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)?

Thank you,
Tony

The format of files generated by fix ave/correlate is described on its doc page in the manual: fix ave/correlate command — LAMMPS documentation

Tony, can please share with the answer to your question, if you find it out?
Thanks

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

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/vol
variable        Jy equal c_flux/vol
variable        Jz equal c_flux/vol

fix             1 all nve
fix             JJ all ave/correlate \$s \$p \$d &
c_flux c_flux c_flux type auto &
file profile.heatflux ave running

variable        scale equal \$s*dt/\$T/\$T/vol
variable        k11 equal trap(f_JJ)*\${scale}
variable        k22 equal trap(f_JJ)*\${scale}
variable        k33 equal trap(f_JJ)*\${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.

Thanks

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.
1 Like