Green-Kubo

Tony · July 12, 2021, 4:16pm

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

akohlmey · July 12, 2021, 4:21pm

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

Damir_Karimov · July 3, 2022, 7:18pm

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

Pragin · July 28, 2022, 5:45am

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.

Thanks

akohlmey · July 28, 2022, 10:27am

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.

Nur · February 13, 2023, 5:44pm

Why do you have no pair coefficient?

Nur

akohlmey · February 13, 2023, 8:59pm

They are included in the data file that is read with read_data.

Nur · February 13, 2023, 10:33pm

Thank you, Axel.

I have checked the LAMMPS example folder of KAPPA (the file in.heatflux.lmp). It uses lj units for LJ liquid. While his script above is real unit for water. According to Green-Kubo formula, I think he missed Boltzmann constant in the “variable scale”. Please correct me if I am wrong. I am also new to LAMMPS.

Nur

akohlmey · February 13, 2023, 10:36pm

I could not tell. The example was not written by me and I have not done this kind of calculation with LAMMPS. However, it should be trivial for you to confirm. Just set up a calculation for a simple system for which the result is known and see if you can reproduce the known result. If it is off by a factor, you could check if that corresponds to k_B.

Nur · February 13, 2023, 10:40pm

In LAMMPS manual, there is an input script example that calculates solid Ar thermal conductivity using real unit. In that, I can see the creator include kb in the “variable scale”. I run it and it gives me 0.26 W/mK, while the LAMMPS manual claims that it will give around 0.29 W/mK.

Nur.

akohlmey · February 13, 2023, 10:50pm

That is the result for running with a single MPI process. Try with 2 or 4.
Or run with OpenMP parallelization instead of MPI.

That should give you a sense for how little the result is statistically converged. The different parallelizations will lead to statistically equivalent, but different trajectories which then also produces results that may vary depending on how well they are converged.

Nur · February 13, 2023, 11:07pm

I used the LAMMPS MSMPI version of 3Nov2022 with 4 processors. Fyi, my processor has only 1 thread for each core, it doesn’t matter if I use the OpenMP or not, right? And I think the results are still acceptable.

Btw, I used the input script in the manual for my flexible SPCE water model suggested by O’Brien here [which is based on Zhang et al., Fluid Phase Equilibria, 262 (2007) 210-216]:

Working with water in LAMMPS - Christopher O’Brien, Ph.D. (google.com)

But for my case the water is in vapor phase which is at 452K. At this temperature, TC = 0.06 W/mK for experimental water. After 0.3ns production (now it is still running), my results gave me only 0.0047 W/mK, I can see it is still decreasing.

What do you think, Axel? Is Green-Kubo appropriate to calcualte TC for water vapor with flexible model?

Nur.

akohlmey · February 13, 2023, 11:11pm

For the sake of demonstration, you can always run with fewer processor cores than the maximum. Also, you can use OpenMP instead of MPI, e.g. run with 1 MPI task and 4 OpenMP threads, or 2 MPI tasks and 2 OpenMP threads. My point was just to demonstrate the volatility of the results (examples bundled with LAMMPS or provided in the manuals are not designed for high accuracy, but to run quickly and demonstrate the principles).

As I already stated, this is not my area of expertise. There are multiple settings and parameters in the input that have an impact on the convergence of the result. I just pointed out the differences by using different parallelizations. The boundary conditions and requirements will be different for any kind of system.

You have to either find an expert that is willing to tutor and advise you, or figure out ways to derive an assessment yourself.

Nur · February 13, 2023, 11:26pm

Thank you for your suggestion, Axel.

Nur.

Abdul_Al_Mahfuj · October 29, 2024, 1:14pm

Hello al I have a query .I want to find the thermal conductivity of Cu using Green-Kubo in LAMMPS. An example in LAMMPS manual for Ar is available. What kind of modification should i done for Cu.
Please help me. I am new to LAMMPS.