# sample LAMMPS input script for viscosity of liquid waterÕ³¶È¼ÆËã
# Green-Kubo method via fix ave/correlate

# settings

variable          den equal density
variable          v equal vol
variable          dt equal 1
variable          t equal 300 
variable          p equal 2000 # correlation length
variable          s equal 5 # sample interval
variable          d equal $p*$s # dump interval

# convert from LAMMPS real units to SI
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}

# problem setup

units		          real
boundary          p p p 
atom_style      	full

read_data         spce.data

pair_style        lj/charmm/coul/long 9 10 10
kspace_style      pppm 1.0e-4                   

pair_coeff        1 1 0.15535 3.166
pair_coeff        * 2 0.0000 0.0000

bond_style        harmonic
angle_style       harmonic
dihedral_style    none
improper_style    none

bond_coeff        1 1000.00 1.000
angle_coeff       1 100.0 109.47

#special_bonds     lj/coul 0.0 0.0 0.5
neighbor          2.0 bin 
neigh_modify      every 1 delay 10 check yes

fix               1 all shake 0.0001 10 1000 b 1 a 1

timestep          ${dt}
thermo            $d
thermo_style      custom step temp press v_den

# equilibration and thermalization

velocity          all create $t 102486 mom yes rot yes dist gaussian
fix               2 all nvt temp $t $t 100.0
run               100000

# viscosity calculation, switch to NVE if desired
#unfix NVT
#fix NVE all nve

reset_timestep    0
variable          pxy equal pxy
variable          pxz equal pxz
variable          pyz equal pyz
variable          tsp equal step
variable          te equal temp
fix               SS all ave/correlate $s $p $d &
                  v_pxy v_pxz v_pyz type auto file $s_S0S.txt 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}
variable          vis equal (v_v11+v_v22+v_v33)/3.0

fix               pr all print $d "${tsp} ${te} ${v11} ${v22} ${v33} ${vis} ${scale}" file vxx_$s.txt 
thermo_style      custom step temp press v_pxy v_pxz v_pyz v_v11 v_v22 v_v33
run               5000000
#write_restart     spce.rst 
#write_data        spce.data

variable          ndens equal count(all)/vol
print             "average viscosity: ${vis} Pa.s/ @ $t K, ${ndens} /A^3"