Added Einstein version of Green-Kubo

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14846 f3b2605a-c512-4ea7-a41b-209d697bcdaa

examples/VISCOSITY/README

@@ -1,8 +1,8 @@
-This directory has 4 scripts that compute the viscosity (eta) of a
+This directory has 5 scripts that compute the viscosity (eta) of a
-Lennard-Jones fluid using 4 different methods.  See the discussion in
+Lennard-Jones fluid using 5 different methods.  See the discussion in
 Section 6.21 of the manual for an overview of the methods and pointers
 to doc pages for the commands which implement them.  Citations for the
-various methods can also be found in the manaul.
+various methods can also be found in the manual.
 
 These scripts are provided for illustration purposes.  No guarantee is
 made that the systems are fully equilibrated or that the runs are long
@@ -10,16 +10,17 @@ enough to generate good statistics and highly accurate results.
 
 -------------
 
-These are the 4 methods for compute viscosity.  The first 3 are
+These are the 5 methods for computing viscosity.  The first 3 are
-non-equilibrium methods; the last is an equilibrium method.
+non-equilibrium methods; the last 2 are equilibrium methods.
 
 in.wall = move a wall to shear the fluid between two walls
 in.nemd = use fix deform and fix nvt/sllod to perform a NEMD shear simulation
 in.mp = use fix viscosity and the Muller-Plathe method
 in.gk = use the Green-Kubo method
+in.einstein = use the Einstein version of Green-Kubo method
 
 All the systems have around 800 atoms.  The NEMD methods run for short
-times; the G-K system needs to run longer to generate good statistics.
+times; the G-K and Einstein systems need to run longer to generate good statistics.
 
 The scripts were all run on a single processor.  They all run in a
 minute or so and produce the accompanying log files and profile files
@@ -78,3 +79,13 @@ heat/flux doc page - the resulting value prints at the end of the run
 and is in the log file
 
 eta = 1.07
+
+(5) in.einstein.2d
+
+eta is computed directly within the script, by performing a time
+integration of the formula discussed in Section 6.21 of the manual,
+analagous to the formula for thermal conductivity given on the compute
+heat/flux doc page - the resulting value prints at the end of the run
+and is in the log file
+
+eta = 1.07

examples/VISCOSITY/in.einstein.2d

@@ -0,0 +1,84 @@
+# sample LAMMPS input script for viscosity of 2d LJ liquid
+# Einstein form of Green-Kubo
+
+# settings
+
+variable	x equal 20
+variable	y equal 20
+
+variable	rho equal 0.6
+variable        t equal 1.0
+variable	rc equal 2.5
+
+variable    p equal 400     # correlation length
+variable    s equal 5       # sample interval
+variable    d equal $p*$s   # dump interval 
+
+# problem setup
+
+units		lj
+dimension	2
+atom_style	atomic
+neigh_modify	delay 0 every 1
+
+lattice         sq2 ${rho}
+region          simbox block 0 $x 0 $y -0.1 0.1
+create_box      1 simbox
+create_atoms    1 box
+
+pair_style      lj/cut ${rc}
+pair_coeff      * * 1 1
+
+mass            * 1.0
+velocity        all create $t 97287
+
+# equilibration run
+
+fix             1 all nve
+fix	        2 all langevin $t $t 0.1 498094
+fix	        3 all enforce2d
+
+thermo          $d
+run	        10000
+
+velocity	all scale $t
+
+unfix		2
+
+# Einstein viscosity calculation
+
+reset_timestep  0
+
+# Define distinct components of symmetric traceless stress tensor
+
+variable         pxy equal pxy
+variable         pxx equal pxx-press
+
+fix              avstress all ave/time $s $p $d v_pxy v_pxx ave one file einstein.dat
+
+# Diagonal components of SS are larger by factor 2-2/d,
+# which is 4/3 for d=3, but 1 for d=2.
+# See Daivis and Evans, J.Chem.Phys, 100, 541-547 (1994)
+
+variable         scale equal vol/(2.0*$t*dt*$d)
+variable         diagfac equal 2-2/2
+variable         deltasqxy equal (f_avstress[1]*$d*dt)^2
+variable         deltasqxx equal (f_avstress[2]*$d*dt)^2/${diagfac}
+
+# compute mean square displacements as running averages
+
+fix              avdeltasq all ave/time $d 1 $d v_deltasqxy v_deltasqxx ave running
+
+# convert to viscosities
+
+variable         vxy equal f_avdeltasq[1]*${scale}
+variable         vxx equal f_avdeltasq[2]*${scale}
+
+thermo_style     custom step temp pe press pxy v_vxy v_vxx
+
+run              500000
+
+variable         etaxy equal v_vxy
+variable         etaxx equal v_vxx
+variable         eta equal 0.5*(${etaxy}+${etaxx})
+print            "running average viscosity: ${eta}"