LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] Green-Kubo method for calculating viscosity
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] Green-Kubo method for calculating viscosity

From: "Smith, Micholas D." <smithmd@...79...>
Date: Mon, 19 Jun 2017 12:06:06 +0000

Viscosity estimates are known to have difficulty converging. You may need to run multiple short runs once the system is at equilibrium and compute the average viscosity from these instead of a single correlation measurement. 


 Yong Zhang, Akihito Otani, and Edward J. Maginn. "Reliable Viscosity Calculation from Equilibrium Molecular Dynamics Simulations: A Time Decomposition Method" JCTC

for suggested best-practicies

Micholas Dean Smith, PhD.
Post-doctoral Research Associate
University of Tennessee/Oak Ridge National Laboratory
Center for Molecular Biophysics

From: Majid Rezaei <majid-rezaei@...3054...>
Sent: Monday, June 19, 2017 6:45 AM
Subject: [lammps-users] Green-Kubo method for calculating viscosity
Dear all,

I use the following Green-Kubo formulation to calculate viscosity of water and NaCl. It seems that it works good. However, when I run the job twice, I obtain different viscosities (the difference is in order of 0.00007 Pa.s which is not negligible). I checked my results and found that this conflict comes from the difference appeared in variable V11 (the other variables are almost the same). 
Could you please help me to find the problem?


units       real
variable    T equal 300
variable    V equal vol
variable    dt equal 4.0
variable    p equal 400     # correlation length
variable    s equal 5       # sample interval
variable    d equal $p*$s   # dump interval

variable    kB equal 1.3806504e-23    # [J/K/** Boltzmann
variable    atm2Pa equal 101325.0
variable    A2m equal 1.0e-10
variable    fs2s equal 1.0e-15
variable    convert equal ${atm2Pa}*${atm2Pa}*${fs2s}*${A2m}*${A2m}*${A2m}

dimension    3
boundary     p p p
atom_style      full
bond_style      harmonic
angle_style     harmonic
pair_style      lj/cut/coul/long 11.0
kspace_style    pppm 1e-4

read_data  data.1

pair_coeff 1 1  0.0 0.0
pair_coeff 2 2  0.155 3.17 
pair_coeff 3 3  0.106 4.45
pair_coeff 4 4  0.0148 2.58   
pair_modify     mix geometric

bond_coeff   1  1000.0 1.0
angle_coeff 1  100.0 109.47
special_bonds   lj/coul 0.0 0.0 0.5

min_style cg
minimize 0 1.0e-6 100 1000    

timestep     ${dt}
thermo       $d

velocity     all create $T 102486 mom yes rot yes dist gaussian

fix          NVT all nvt temp $T $T 10 drag 0.2

fix        66 all shake 0.0001 20 0 b 1 a 1

run          1000000

reset_timestep 0
variable     pxy equal pxy
variable     pxz equal pxz
variable     pyz equal pyz
fix          SS all ave/correlate $s $p $d &
             v_pxy v_pxz v_pyz type auto file S0St.dat ave running
variable     scale equal ${convert}/(${kB}*$T)*$V*$s*${dt}
variable     v11 equal trap(f_SS[3])*${scale}
variable     v22 equal trap(f_SS[4])*${scale}
variable     v33 equal trap(f_SS[5])*${scale}
thermo_style custom step temp press v_pxy v_pxz v_pyz v_v11 v_v22 v_v33

run          5000000

variable     v equal (v_v11+v_v22+v_v33)/3.0
variable     ndens equal count(all)/vol
print        "average viscosity: $v [Pa.s] @ $T K, ${ndens} /A^3"