Syntax:
-fix ID group-ID viscosity N vdim pdim Nbin +fix ID group-ID viscosity N vdim pdim Nbin keyword value ...-
swap value = Nswap = number of swaps to perform every N steps + vtarget value = V or INF = target velocity of swap partners (velocity units) ++
Examples:
-fix 1 all viscosity 100 x z 20 +fix 1 all viscosity 100 x z 20 +fix 1 all viscosity 50 x z 20 swap 2 vtarget 1.5Description:
@@ -42,27 +57,36 @@ system's response.The simulation box is divided into Nbin layers in the pdim -direction. Every N steps, two atoms are chosen in the following -manner. Only atoms in the fix group are considered. The atom in the -bottom layer with the most positive momentum component in the vdim -direction is the first atom. The atom in the middle later with the -most negative momentum component in the vdim direction is the second -atom. The vdim momenta components of these two atoms are swapped, -which resets their velocities, typically in opposite directions. Over +direction. Every N steps, Nswap pairs of atoms are chosen in the +following manner. Only atoms in the fix group are considered. Nswap +atoms in the bottom layer with positive velocity components in the +vdim direction closest to the target value V are selected. +Similarly, Nswap atoms in the middle later with negative velocity +components in the vdim direction closest to the negative of the +target value V are selected. The two sets of Nswap atoms are paired +up and their vdim momenta components are swapped within each pair. +This resets their velocities, typically in opposite directions. Over time, this induces a shear velocity profile in the system which can be measured using commands such as the following, which writes the profile to the file tmp.profile:
-compute c1 all attribute/atom vx -fix f1 all ave/spatial 100 10 1000 z lower 0.05 tmp.profile & - compute c1 units reduced +fix f1 all ave/spatial 100 10 1000 z lower 0.05 vx & + file tmp.profile units reduced+Note that by default, Nswap = 1 and vtarget = INF, though this can be +changed by the optional swap and vtarget keywords. When vtarget = +INF, one or more atoms with the most positive and negative velocity +components are selected. Setting these parameters appropriately, in +conjunction with the swap rate N, allows the momentum flux rate to be +adjusted across a wide range of values, and the momenta to be +exchanged in large chunks or more smoothly. +
As described below, the total momentum transferred by these velocity swaps is computed by the fix and can be output. Dividing this quantity by time and the cross-sectional area of the simulation box yields a momentum flux. The ratio of momentum flux to the slope of the shear velocity profile is the viscosity of the fluid, in -appropriate units. See the Muller-Plathe paper for +appopriate units. See the Muller-Plathe paper for details.
IMPORTANT NOTE: After equilibration, if the velocity profile you @@ -73,7 +97,7 @@
An alternative method for calculating a viscosity is to run a NEMD simulation, as described in this section of -the manual. NEMD simulations deform the simulation box via the fix +the manual. NEMD simulations deform the simmulation box via the fix deform command. Thus they cannot be run on a charged system using a PPPM solver since PPPM does not currently support non-orthogonal boxes. Using fix viscosity keeps the @@ -85,10 +109,10 @@ files. None of the fix_modify options are relevant to this fix.
-The cumulative momentum transferred between the bottom and middle of +
The cummulative momentum transferred between the bottom and middle of the simulation box (in the pdim direction) is stored as a scalar quantity by this fix. This quantity is zeroed when the fix is defined -and accumulates thereafter, once every N steps. The units of the +and accumlates thereafter, once every N steps. The units of the quantity are momentum = mass*velocity. This quantity can be accessed by various output commands, such as thermo_style custom. The scalar value calculated @@ -115,7 +139,7 @@ in a computation of alcohol molecule properties.
When running a simulation with large, massive particles or molecules -in a background solvent, you may want to only exchange momenta between +in a background solvent, you may want to only exchange momenta bewteen solvent particles.
Related commands: @@ -123,7 +147,9 @@
fix ave/spatial, fix nvt/sllod
-Default: none +
Default: +
+The option defaults are swap = 1 and vtarget = INF.
diff -Naur lammps-30Mar08/doc/fix_viscosity.txt lammps-31Mar08/doc/fix_viscosity.txt --- lammps-30Mar08/doc/fix_viscosity.txt 2008-02-29 18:13:20.000000000 -0700 +++ lammps-31Mar08/doc/fix_viscosity.txt 2008-03-28 09:45:51.000000000 -0600 @@ -10,18 +10,25 @@ [Syntax:] -fix ID group-ID viscosity N vdim pdim Nbin :pre +fix ID group-ID viscosity N vdim pdim Nbin keyword value ... :pre -ID, group-ID are documented in "fix"_fix.html command -viscosity = style name of this fix command -N = perform momentum exchange every N steps -vdim = {x} or {y} or {z} = which momentum component to exchange -pdim = {x} or {y} or {z} = direction of momentum transfer -Nbin = # of layers in pdim direction :ul +ID, group-ID are documented in "fix"_fix.html command :ulb,l +viscosity = style name of this fix command :l +N = perform momentum exchange every N steps :l +vdim = {x} or {y} or {z} = which momentum component to exchange :l +pdim = {x} or {y} or {z} = direction of momentum transfer :l +Nbin = # of layers in pdim direction :l + +zero or more keyword/value pairs may be appended :l +keyword = {swap} or {target} :l + {swap} value = Nswap = number of swaps to perform every N steps + {vtarget} value = V or INF = target velocity of swap partners (velocity units) :pre +:ule [Examples:] -fix 1 all viscosity 100 x z 20 :pre +fix 1 all viscosity 100 x z 20 +fix 1 all viscosity 50 x z 20 swap 2 vtarget 1.5 :pre [Description:] @@ -39,27 +46,36 @@ system's response. The simulation box is divided into {Nbin} layers in the {pdim} -direction. Every N steps, two atoms are chosen in the following -manner. Only atoms in the fix group are considered. The atom in the -bottom layer with the most positive momentum component in the {vdim} -direction is the first atom. The atom in the middle later with the -most negative momentum component in the {vdim} direction is the second -atom. The {vdim} momenta components of these two atoms are swapped, -which resets their velocities, typically in opposite directions. Over +direction. Every N steps, Nswap pairs of atoms are chosen in the +following manner. Only atoms in the fix group are considered. Nswap +atoms in the bottom layer with positive velocity components in the +{vdim} direction closest to the target value {V} are selected. +Similarly, Nswap atoms in the middle later with negative velocity +components in the {vdim} direction closest to the negative of the +target value {V} are selected. The two sets of Nswap atoms are paired +up and their {vdim} momenta components are swapped within each pair. +This resets their velocities, typically in opposite directions. Over time, this induces a shear velocity profile in the system which can be measured using commands such as the following, which writes the profile to the file tmp.profile: -compute c1 all attribute/atom vx -fix f1 all ave/spatial 100 10 1000 z lower 0.05 tmp.profile & - compute c1 units reduced :pre +fix f1 all ave/spatial 100 10 1000 z lower 0.05 vx & + file tmp.profile units reduced :pre + +Note that by default, Nswap = 1 and vtarget = INF, though this can be +changed by the optional {swap} and {vtarget} keywords. When vtarget = +INF, one or more atoms with the most positive and negative velocity +components are selected. Setting these parameters appropriately, in +conjunction with the swap rate N, allows the momentum flux rate to be +adjusted across a wide range of values, and the momenta to be +exchanged in large chunks or more smoothly. As described below, the total momentum transferred by these velocity swaps is computed by the fix and can be output. Dividing this quantity by time and the cross-sectional area of the simulation box yields a momentum flux. The ratio of momentum flux to the slope of the shear velocity profile is the viscosity of the fluid, in -appropriate units. See the "Muller-Plathe paper"_#Muller-Plathe for +appopriate units. See the "Muller-Plathe paper"_#Muller-Plathe for details. IMPORTANT NOTE: After equilibration, if the velocity profile you @@ -70,7 +86,7 @@ An alternative method for calculating a viscosity is to run a NEMD simulation, as described in "this section"_Section_howto.html#4_13 of -the manual. NEMD simulations deform the simulation box via the "fix +the manual. NEMD simulations deform the simmulation box via the "fix deform"_fix_deform.html command. Thus they cannot be run on a charged system using a "PPPM solver"_kspace_style.html since PPPM does not currently support non-orthogonal boxes. Using fix viscosity keeps the @@ -82,10 +98,10 @@ files"_restart.html. None of the "fix_modify"_fix_modify.html options are relevant to this fix. -The cumulative momentum transferred between the bottom and middle of +The cummulative momentum transferred between the bottom and middle of the simulation box (in the {pdim} direction) is stored as a scalar quantity by this fix. This quantity is zeroed when the fix is defined -and accumulates thereafter, once every N steps. The units of the +and accumlates thereafter, once every N steps. The units of the quantity are momentum = mass*velocity. This quantity can be accessed by various "output commands"_Section_howto.html#4_15, such as "thermo_style custom"_thermo_style.html. The scalar value calculated @@ -112,7 +128,7 @@ in a computation of alcohol molecule properties. When running a simulation with large, massive particles or molecules -in a background solvent, you may want to only exchange momenta between +in a background solvent, you may want to only exchange momenta bewteen solvent particles. [Related commands:] @@ -120,7 +136,9 @@ "fix ave/spatial"_fix_ave_spatial.html, "fix nvt/sllod"_fix_nvt_sllod.html -[Default:] none +[Default:] + +The option defaults are swap = 1 and vtarget = INF. :line diff -Naur lammps-30Mar08/src/fix_viscosity.cpp lammps-31Mar08/src/fix_viscosity.cpp --- lammps-30Mar08/src/fix_viscosity.cpp 2008-02-06 18:45:42.000000000 -0700 +++ lammps-31Mar08/src/fix_viscosity.cpp 2008-03-28 09:47:41.000000000 -0600 @@ -11,6 +11,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +#include "math.h" #include "mpi.h" #include "string.h" #include "stdlib.h" @@ -21,14 +22,19 @@ using namespace LAMMPS_NS; -#define BIG 1.0e20 +// needs to be big, but not so big that lose precision when subtract velocity + +#define BIG 1.0e10 + +#define MIN(A,B) ((A) < (B)) ? (A) : (B) +#define MAX(A,B) ((A) > (B)) ? (A) : (B) /* ---------------------------------------------------------------------- */ FixViscosity::FixViscosity(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { - if (narg != 7) error->all("Illegal fix viscosity command"); + if (narg < 7) error->all("Illegal fix viscosity command"); MPI_Comm_rank(world,&me); @@ -52,11 +58,50 @@ nbin = atoi(arg[6]); if (nbin < 3) error->all("Illegal fix viscosity command"); + // optional keywords + + nswap = 1; + vtarget = BIG; + + int iarg = 7; + while (iarg < narg) { + if (strcmp(arg[iarg],"swap") == 0) { + if (iarg+2 > narg) error->all("Illegal fix viscosity command"); + nswap = atoi(arg[iarg+1]); + if (nswap <= 0) error->all("Fix viscosity swap value must be positive"); + iarg += 2; + } else if (strcmp(arg[iarg],"vtarget") == 0) { + if (iarg+2 > narg) error->all("Illegal fix viscosity command"); + if (strcmp(arg[iarg+1],"INF") == 0) vtarget = BIG; + else vtarget = atof(arg[iarg+1]); + if (vtarget <= 0.0) + error->all("Fix viscosity vtarget value must be positive"); + iarg += 2; + } else error->all("Illegal fix viscosity command"); + } + + // initialize array sizes to nswap+1 so have space to shift values down + + pos_index = new int[nswap+1]; + neg_index = new int[nswap+1]; + pos_delta = new double[nswap+1]; + neg_delta = new double[nswap+1]; + flux = 0.0; } /* ---------------------------------------------------------------------- */ +FixViscosity::~FixViscosity() +{ + delete [] pos_index; + delete [] neg_index; + delete [] pos_delta; + delete [] neg_delta; +} + +/* ---------------------------------------------------------------------- */ + int FixViscosity::setmask() { int mask = 0; @@ -92,8 +137,8 @@ void FixViscosity::end_of_step() { - int i; - double p,coord; + int i,m,insert; + double p,coord,delta; MPI_Status status; struct { double value; @@ -113,8 +158,10 @@ slabhi_hi = boxlo + ((nbin-1)/2 + 1)*binsize; } - // ipos,ineg = my 2 atoms with most pos/neg momenta in bottom/top slabs - // map atom back into periodic box if necessary + // make list of nswap atoms with velocity closest to +/- vtarget + // only consider atoms in the bottom/middle slabs + // map atoms back into periodic box if necessary + // insert = location in list to insert new atom double **x = atom->x; double **v = atom->v; @@ -124,78 +171,113 @@ int *mask = atom->mask; int nlocal = atom->nlocal; - int ipos = -1; - int ineg = -1; - double pmin = BIG; - double pmax = -BIG; + npositive = nnegative = 0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if (mass) p = mass[type[i]] * v[i][vdim]; - else p = rmass[i] * v[i][vdim]; coord = x[i][pdim]; if (coord < boxlo && periodicity) coord += prd; else if (coord >= boxhi && periodicity) coord -= prd; + if (coord >= slablo_lo && coord < slablo_hi) { - if (p > pmax) { - pmax = p; - ipos = i; + if (v[i][vdim] < 0.0) continue; + delta = fabs(v[i][vdim] - vtarget); + if (npositive < nswap || delta < pos_delta[nswap-1]) { + for (insert = npositive-1; insert >= 0; insert--) + if (delta > pos_delta[insert]) break; + insert++; + for (m = npositive-1; m >= insert; m--) { + pos_delta[m+1] = pos_delta[m]; + pos_index[m+1] = pos_index[m]; + } + pos_delta[insert] = delta; + pos_index[insert] = i; + if (npositive < nswap) npositive++; } } + if (coord >= slabhi_lo && coord < slabhi_hi) { - if (p < pmin) { - pmin = p; - ineg = i; + if (v[i][vdim] > 0.0) continue; + delta = fabs(v[i][vdim] + vtarget); + if (nnegative < nswap || delta < neg_delta[nswap-1]) { + for (insert = nnegative-1; insert >= 0; insert--) + if (delta > neg_delta[insert]) break; + insert++; + for (m = nnegative-1; m >= insert; m--) { + neg_delta[m+1] = neg_delta[m]; + neg_index[m+1] = neg_index[m]; + } + neg_delta[insert] = delta; + neg_index[insert] = i; + if (nnegative < nswap) nnegative++; } } } - // find 2 global atoms with most pos/neg momenta in bottom/top slabs - // MAXLOC also communicates which procs own them - - mine[0].value = pmax; - mine[0].proc = me; - mine[1].value = -pmin; - mine[1].proc = me; - - MPI_Allreduce(mine,all,2,MPI_DOUBLE_INT,MPI_MAXLOC,world); - + // loop over nswap pairs + // find 2 global atoms with smallest delta in bottom/top slabs + // MINLOC also communicates which procs own them // exchange momenta between the 2 particles - // if I own both particles just swap, else point2point comm of mass,vel - + // if I own both particles just swap, else point2point comm of vel,mass + + int ipos,ineg; double sbuf[2],rbuf[2]; - if (me == all[0].proc && me == all[1].proc) { - sbuf[0] = v[ineg][vdim]; - if (mass) sbuf[1] = mass[type[ineg]]; - else sbuf[1] = rmass[ineg]; - rbuf[0] = v[ipos][vdim]; - if (mass) rbuf[1] = mass[type[ipos]]; - else rbuf[1] = rmass[ipos]; - v[ineg][vdim] = rbuf[0] * rbuf[1]/sbuf[1]; - v[ipos][vdim] = sbuf[0] * sbuf[1]/rbuf[1]; - - } else if (me == all[0].proc) { - sbuf[0] = v[ipos][vdim]; - if (mass) sbuf[1] = mass[type[ipos]]; - else sbuf[1] = rmass[ipos]; - MPI_Sendrecv(sbuf,2,MPI_DOUBLE,all[1].proc,0, - rbuf,2,MPI_DOUBLE,all[1].proc,0,world,&status); - v[ipos][vdim] = rbuf[0] * rbuf[1]/sbuf[1]; - - } else if (me == all[1].proc) { - sbuf[0] = v[ineg][vdim]; - if (mass) sbuf[1] = mass[type[ineg]]; - else sbuf[1] = rmass[ineg]; - MPI_Sendrecv(sbuf,2,MPI_DOUBLE,all[0].proc,0, - rbuf,2,MPI_DOUBLE,all[0].proc,0,world,&status); - v[ineg][vdim] = rbuf[0] * rbuf[1]/sbuf[1]; + double dflux = 0.0; + mine[0].proc = mine[1].proc = me; + int ipositive = 0; + int inegative = 0; + + int nactual = MIN(nnegative,npositive); + + for (m = 0; m < nactual; m++) { + if (ipositive < npositive) mine[0].value = pos_delta[ipositive]; + else mine[0].value = BIG; + if (inegative < nnegative) mine[1].value = neg_delta[inegative]; + else mine[1].value = BIG; + + MPI_Allreduce(mine,all,2,MPI_DOUBLE_INT,MPI_MINLOC,world); + + if (me == all[0].proc && me == all[1].proc) { + ipos = pos_index[ipositive++]; + ineg = neg_index[inegative++]; + rbuf[0] = v[ipos][vdim]; + if (mass) rbuf[1] = mass[type[ipos]]; + else rbuf[1] = rmass[ipos]; + sbuf[0] = v[ineg][vdim]; + if (mass) sbuf[1] = mass[type[ineg]]; + else sbuf[1] = rmass[ineg]; + v[ineg][vdim] = rbuf[0] * rbuf[1]/sbuf[1]; + v[ipos][vdim] = sbuf[0] * sbuf[1]/rbuf[1]; + dflux += rbuf[0]*rbuf[1] - sbuf[0]*sbuf[1]; + + } else if (me == all[0].proc) { + ipos = pos_index[ipositive++]; + sbuf[0] = v[ipos][vdim]; + if (mass) sbuf[1] = mass[type[ipos]]; + else sbuf[1] = rmass[ipos]; + MPI_Sendrecv(sbuf,2,MPI_DOUBLE,all[1].proc,0, + rbuf,2,MPI_DOUBLE,all[1].proc,0,world,&status); + v[ipos][vdim] = rbuf[0] * rbuf[1]/sbuf[1]; + dflux += sbuf[0]*sbuf[1]; + + } else if (me == all[1].proc) { + ineg = neg_index[inegative++]; + sbuf[0] = v[ineg][vdim]; + if (mass) sbuf[1] = mass[type[ineg]]; + else sbuf[1] = rmass[ineg]; + MPI_Sendrecv(sbuf,2,MPI_DOUBLE,all[0].proc,0, + rbuf,2,MPI_DOUBLE,all[0].proc,0,world,&status); + v[ineg][vdim] = rbuf[0] * rbuf[1]/sbuf[1]; + dflux -= sbuf[0]*sbuf[1]; + } } - // tally momentum flux - // sign of all[1].value was flipped for MPI_Allreduce + // tally momentum flux from all swaps - flux += all[0].value + all[1].value; + double dflux_all; + MPI_Allreduce(&dflux,&dflux_all,1,MPI_DOUBLE,MPI_SUM,world); + flux += dflux_all; } /* ---------------------------------------------------------------------- */ diff -Naur lammps-30Mar08/src/fix_viscosity.h lammps-31Mar08/src/fix_viscosity.h --- lammps-30Mar08/src/fix_viscosity.h 2007-10-16 12:08:31.000000000 -0600 +++ lammps-31Mar08/src/fix_viscosity.h 2008-03-28 09:47:41.000000000 -0600 @@ -21,7 +21,7 @@ class FixViscosity : public Fix { public: FixViscosity(class LAMMPS *, int, char **); - ~FixViscosity() {} + ~FixViscosity(); int setmask(); void init(); void end_of_step(); @@ -30,9 +30,15 @@ private: int me; int vdim,pdim,nbin,periodicity; + int nswap; + double vtarget; double prd,boxlo,boxhi; double slablo_lo,slablo_hi,slabhi_lo,slabhi_hi; double flux; + + int npositive,nnegative; + int *pos_index,*neg_index; + double *pos_delta,*neg_delta; }; }