From 00f3ccf3b017765d024e594e828de6091ac8326d Mon Sep 17 00:00:00 2001 From: athomps Date: Mon, 18 Apr 2016 23:59:16 +0000 Subject: [PATCH] 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 | 23 ++++++--- examples/VISCOSITY/in.einstein.2d | 84 +++++++++++++++++++++++++++++++ 2 files changed, 101 insertions(+), 6 deletions(-) create mode 100644 examples/VISCOSITY/in.einstein.2d diff --git a/examples/VISCOSITY/README b/examples/VISCOSITY/README index 78ad48bfc1..88b874cd65 100644 --- a/examples/VISCOSITY/README +++ b/examples/VISCOSITY/README @@ -1,8 +1,8 @@ -This directory has 4 scripts that compute the viscosity (eta) of a -Lennard-Jones fluid using 4 different methods. See the discussion in +This directory has 5 scripts that compute the viscosity (eta) of a +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 -non-equilibrium methods; the last is an equilibrium method. +These are the 5 methods for computing viscosity. The first 3 are +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 diff --git a/examples/VISCOSITY/in.einstein.2d b/examples/VISCOSITY/in.einstein.2d new file mode 100644 index 0000000000..72ccf5d7a0 --- /dev/null +++ b/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}"