small edits for reformatting

src/USER-MISC/fix_wall_reflect_stochastic.cpp

@@ -19,40 +19,40 @@
 #include "update.h"
 #include "modify.h"
 #include "domain.h"
+#include "force.h"
 #include "lattice.h"
 #include "input.h"
 #include "variable.h"
-#include "error.h"
-#include "force.h"
 #include "random_mars.h"
-
+#include "error.h"
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
 
 enum{NONE=0,DIFFUSIVE=1,MAXWELL=2,CERCIGNANILAMPIS=3};
 
-
 /* ---------------------------------------------------------------------- */
 
-FixWallReflectStochastic::FixWallReflectStochastic(LAMMPS *lmp, int narg, char **arg) :
+FixWallReflectStochastic::
+FixWallReflectStochastic(LAMMPS *lmp, int narg, char **arg) :
   FixWallReflect(lmp, narg, arg), random(NULL)
 {
-  int arginc,dir;
-
   if (narg < 8) error->all(FLERR,"Illegal fix wall/reflect/stochastic command");
   
-  if (domain->triclinic != 0) error->all(FLERR, "This fix wall/reflect/stochastic is not incompatible with triclinic simulation box");
+  if (domain->triclinic != 0) 
+    error->all(FLERR, "Fix wall/reflect/stochastic cannot be used with "
+               "triclinic simulation box");
 
   dynamic_group_allow = 1;
 
   // parse args
 
+  int arginc;
+
   nwall = 0;
   int scaleflag = 1;
   reflectionstyle = NONE;
   
-  
   if (strcmp(arg[3],"diffusive") == 0) {
     reflectionstyle = DIFFUSIVE;
     arginc = 6;
@@ -66,15 +66,15 @@ FixWallReflectStochastic::FixWallReflectStochastic(LAMMPS *lmp, int narg, char *
 
 
   seedfix = force->inumeric(FLERR,arg[4]);
-    if (seedfix <=0) error->all(FLERR,"Random seed must be a postive number");
+  if (seedfix <= 0) error->all(FLERR,"Random seed must be a postive number");
-    random = new RanMars(lmp,seedfix + comm->me);
 
   int iarg = 5;
   while (iarg < narg) {
     if ((strcmp(arg[iarg],"xlo") == 0) || (strcmp(arg[iarg],"xhi") == 0) ||
         (strcmp(arg[iarg],"ylo") == 0) || (strcmp(arg[iarg],"yhi") == 0) ||
         (strcmp(arg[iarg],"zlo") == 0) || (strcmp(arg[iarg],"zhi") == 0)) {
-        if (iarg+2 > narg) error->all(FLERR,"Illegal fix wall/reflect/stochastic command");
+      if (iarg+2 > narg) 
+        error->all(FLERR,"Illegal fix wall/reflect/stochastic command");
 
       int newwall;
       if (strcmp(arg[iarg],"xlo") == 0) newwall = XLO;
@@ -86,7 +86,8 @@ FixWallReflectStochastic::FixWallReflectStochastic(LAMMPS *lmp, int narg, char *
       
       for (int m = 0; (m < nwall) && (m < 6); m++)
         if (newwall == wallwhich[m])
-            error->all(FLERR,"Wall defined twice in fix wall/reflect/stochastic command");
+          error->all(FLERR,
+                     "Fix wall/reflect/stochastic command wall defined twice");
 
       wallwhich[nwall] = newwall;
       if (strcmp(arg[iarg+1],"EDGE") == 0) {
@@ -100,18 +101,19 @@ FixWallReflectStochastic::FixWallReflectStochastic(LAMMPS *lmp, int narg, char *
         coord0[nwall] = force->numeric(FLERR,arg[iarg+1]);
       }
 
-
       walltemp[nwall]= force->numeric(FLERR,arg[iarg+2]);
 
-        for (dir = 0; dir < 3; dir++) {
+      for (int dir = 0; dir < 3; dir++) {
         wallvel[nwall][dir]= force->numeric(FLERR,arg[iarg+dir+3]);
         int dim = wallwhich[nwall] / 2;
-          if ((wallvel[nwall][dir] !=0) & (dir == dim)) error->all(FLERR,"The wall velocity must be tangential");
+        if ((wallvel[nwall][dir] !=0) & (dir == dim)) 
+          error->all(FLERR,"The wall velocity must be tangential");
         
         if (reflectionstyle == CERCIGNANILAMPIS) {
           wallaccom[nwall][dir]= force->numeric(FLERR,arg[iarg+dir+6]);
         } else if (reflectionstyle == MAXWELL) {
-            wallaccom[nwall][dir]= force->numeric(FLERR,arg[iarg+6]);// One accommodation coefficient for all directions
+          // one accommodation coefficient for all directions
+          wallaccom[nwall][dir]= force->numeric(FLERR,arg[iarg+6]);
         }
       }
       
@@ -119,7 +121,8 @@ FixWallReflectStochastic::FixWallReflectStochastic(LAMMPS *lmp, int narg, char *
       iarg += arginc;
       
     } else if (strcmp(arg[iarg],"units") == 0) {
-        if (iarg+2 > narg) error->all(FLERR,"Illegal wall/reflect/stochastic command");
+      if (iarg+2 > narg) 
+        error->all(FLERR,"Illegal wall/reflect/stochastic command");
       if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
       else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
       else error->all(FLERR,"Illegal fix wall/reflect/stochastic command");
@@ -133,17 +136,21 @@ FixWallReflectStochastic::FixWallReflectStochastic(LAMMPS *lmp, int narg, char *
 
   for (int m = 0; m < nwall; m++) {
     if ((wallwhich[m] == XLO || wallwhich[m] == XHI) && domain->xperiodic)
-        error->all(FLERR,"Cannot use fix wall/reflect/stochastic in periodic dimension");
+      error->all(FLERR,"Cannot use fix wall/reflect/stochastic "
+                 "in periodic dimension");
     if ((wallwhich[m] == YLO || wallwhich[m] == YHI) && domain->yperiodic)
-        error->all(FLERR,"Cannot use fix wall/reflect/stochastic in periodic dimension");
+      error->all(FLERR,"Cannot use fix wall/reflect/stochastic "
+                 "in periodic dimension");
     if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->zperiodic)
-        error->all(FLERR,"Cannot use fix wall/reflect/stochastic in periodic dimension");
+      error->all(FLERR,"Cannot use fix wall/reflect/stochastic "
+                 "in periodic dimension");
   }
   
   for (int m = 0; m < nwall; m++)
     if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->dimension == 2)
       error->all(FLERR,
-                   "Cannot use fix wall/reflect/stochastic zlo/zhi for a 2d simulation");
+                 "Cannot use fix wall/reflect/stochastic zlo/zhi "
+                 "for a 2d simulation");
   
   // scale factors for CONSTANT walls, VARIABLE wall is not allowed
   
@@ -167,14 +174,16 @@ FixWallReflectStochastic::FixWallReflectStochastic(LAMMPS *lmp, int narg, char *
     }
   }
 
-  }
+  // random number generator
 
+  random = new RanMars(lmp,seedfix + comm->me);
+}
+
+/* ---------------------------------------------------------------------- */
 
 FixWallReflectStochastic::~FixWallReflectStochastic()
 {
-
+  delete random;
-delete random ;
-
 }
 
 /* ---------------------------------------------------------------------- */
@@ -188,7 +197,6 @@ void FixWallReflectStochastic::wall_particle(int m, int which, double coord)
   double *mass = atom->mass;
 
   // coord = current position of wall
-  // evaluate variable if necessary, wrap with clear/add
 
   double **x = atom->x;
   double **v = atom->v;
@@ -200,67 +208,67 @@ void FixWallReflectStochastic::wall_particle(int m, int which, double coord)
   side = which % 2;
 
   for (i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) {
+    if (!(mask[i] & groupbit)) continue;
       
     sign = 0;
+    if ((side == 0) & (x[i][dim] < coord)) sign = 1;
+    else if ((side == 1) & (x[i][dim] > coord)) sign = -1;
+    if (sign == 0) continue;
 
-      if ((side == 0) & (x[i][dim] < coord)){
+    // theta = kT/m
-        sign = 1;
-      } else if ((side == 1) & (x[i][dim] > coord)){
-        sign = -1;
-      }
 
-
+    if (rmass) theta = force->boltz*walltemp[m]/(rmass[i]*force->mvv2e);
-      if (sign != 0) {
-        if (rmass) theta = force->boltz*walltemp[m]/(rmass[i]*force->mvv2e);// theta = kT/m
     else theta = force->boltz*walltemp[m]/(mass[type[i]]*force->mvv2e);
     factor = sqrt(theta);
       
-        timecol= (x[i][dim] - coord)/v[i][dim]; // time travelling after collision
+    // time travelling after collision
 
-        if (reflectionstyle  == MAXWELL) {difftest = random->uniform(); }// for Maxwell model only
+    timecol = (x[i][dim]-coord) / v[i][dim]; 
+
+    // only needed for Maxwell model
+    
+    if (reflectionstyle == MAXWELL) difftest = random->uniform();
     
     for (dir = 0; dir < 3; dir++) {
 
-          x[i][dir] -= v[i][dir]*timecol; // pushing back atoms to the wall position and assign new reflected velocity
+      // pushing back atoms to wall position, assign new reflected velocity
       
-          // Diffusive reflection
+      x[i][dir] -= v[i][dir]*timecol; 
-          if (reflectionstyle  == DIFFUSIVE){
       
-            if (dir != dim) {
+      // diffusive reflection
+      
+      if (reflectionstyle  == DIFFUSIVE) {
+        if (dir != dim) 
           v[i][dir] = wallvel[m][dir] + random->gaussian(0,factor);
-            } else { v[i][dir] =  sign*random->rayleigh(factor) ; }
+        else v[i][dir] =  sign*random->rayleigh(factor);
-
-          }
         
       // Maxwell reflection
-          if (reflectionstyle  == MAXWELL){
           
+      } else if (reflectionstyle  == MAXWELL) {
         if (difftest < wallaccom[m][dir]) {
-              if (dir != dim) {
+          if (dir != dim) 
             v[i][dir] = wallvel[m][dir] + random->gaussian(0,factor);
-              } else { v[i][dir] =  sign*random->rayleigh(factor) ; }
+          else v[i][dir] =  sign*random->rayleigh(factor);
         } else {
-              if (dir == dim) {v[i][dir] = -v[i][dir];
+          if (dir == dim) v[i][dir] = -v[i][dir];
         }
-            }
-          }
-
 
       // Cercignani Lampis reflection
-          if (reflectionstyle  == CERCIGNANILAMPIS){
-            if (dir != dim) {
-              v[i][dir] = wallvel[m][dir] + random->gaussian((1-wallaccom[m][dir])*(v[i][dir]-wallvel[m][dir]),sqrt((2-wallaccom[m][dir])*wallaccom[m][dir]*theta)) ;
-            } else { v[i][dir] = random->besselexp(theta, wallaccom[m][dir], v[i][dir]) ; }
 
+      } else if (reflectionstyle  == CERCIGNANILAMPIS) {
+        if (dir != dim)
+          v[i][dir] = wallvel[m][dir] + 
+            random->gaussian((1-wallaccom[m][dir]) * 
+                             (v[i][dir] - wallvel[m][dir]),
+                             sqrt((2-wallaccom[m][dir]) * 
+                                  wallaccom[m][dir]*theta));
+        else v[i][dir] = random->besselexp(theta,wallaccom[m][dir],v[i][dir]);
       }
 
       // update new position due to the new velocity
-          if (dir != dim) {
+      
-            x[i][dir] += v[i][dir]*timecol; }
+      if (dir != dim) x[i][dir] += v[i][dir]*timecol;
-          else {x[i][dir] = coord + v[i][dir]*timecol;}
+      else x[i][dir] = coord + v[i][dir]*timecol;
-        }// for loop
+    }
-      }// if sign
-    }// if mask
   }
 }

src/USER-MISC/fix_wall_reflect_stochastic.h

@@ -20,7 +20,6 @@ FixStyle(wall/reflect/stochastic,FixWallReflectStochastic)
 #ifndef LMP_FIX_WALL_REFLECT_STOCHASTIC_H
 #define LMP_FIX_WALL_REFLECT_STOCHASTIC_H
 
-
 #include "random_mars.h"
 #include "fix_wall_reflect.h"
 
@@ -29,20 +28,16 @@ namespace LAMMPS_NS {
 class FixWallReflectStochastic : public FixWallReflect {
  public:
   FixWallReflectStochastic(class LAMMPS *, int, char **);
-  void wall_particle(int m,int which, double coord);
   virtual ~FixWallReflectStochastic();
 
-
+ private:
-
-
-
- protected:
-    class RanMars *random;
   int seedfix;
   double walltemp[6],wallvel[6][3],wallaccom[6][3];
   int reflectionstyle;
 
+  class RanMars *random;
 
+  void wall_particle(int m,int which, double coord);
 };
 
 }

src/fix_wall_reflect.cpp

@@ -18,30 +18,26 @@
 #include "update.h"
 #include "modify.h"
 #include "domain.h"
+#include "force.h"
 #include "lattice.h"
 #include "input.h"
 #include "variable.h"
 #include "error.h"
-#include "force.h"
-
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
 
-
-
 /* ---------------------------------------------------------------------- */
 
 FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg),
   nwall(0)
 { 
-  
-
   if (narg < 4) error->all(FLERR,"Illegal fix wall/reflect command");
 
-  if (strcmp(arg[2],"wall/reflect/stochastic") == 0) return;// child class can be stochastic
+  // let child class process all args
 
+  if (strcmp(arg[2],"wall/reflect/stochastic") == 0) return;
 
   dynamic_group_allow = 1;
 
@@ -89,13 +85,13 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
       nwall++;
       iarg += 2;
 
-
     } else if (strcmp(arg[iarg],"units") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal wall/reflect command");
       if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
       else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
       else error->all(FLERR,"Illegal fix wall/reflect command");
       iarg += 2;
+
     } else error->all(FLERR,"Illegal fix wall/reflect command");
   }
 
@@ -146,8 +142,6 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
     if (wallstyle[m] == VARIABLE) varflag = 1;
 }
 
-
-
 /* ---------------------------------------------------------------------- */
 
 FixWallReflect::~FixWallReflect()
@@ -194,14 +188,11 @@ void FixWallReflect::init()
 
 void FixWallReflect::post_integrate()
 {
-  //int m;
   double coord;
 
   // coord = current position of wall
   // evaluate variable if necessary, wrap with clear/add
 
-
-
   if (varflag) modify->clearstep_compute();
 
   for (int m = 0; m < nwall; m++) {
@@ -212,17 +203,19 @@ void FixWallReflect::post_integrate()
       else coord *= zscale;
     } else coord = coord0[m];
 
-
     wall_particle(m,wallwhich[m],coord);
-
+  }
 
   if (varflag) modify->addstep_compute(update->ntimestep + 1);
 }
-}
+
+/* ----------------------------------------------------------------------
+   this method may be overwritten by a child class
+------------------------------------------------------------------------- */
 
 void FixWallReflect::wall_particle(int m, int which, double coord)
 {
-  int i, dim, side,sign;
+  int i,dim,side;
 
   double **x = atom->x;
   double **v = atom->v;
@@ -232,9 +225,8 @@ void FixWallReflect::wall_particle(int m, int which, double coord)
   dim = which / 2;
   side = which % 2;
 
-    for (i = 0; i < nlocal; i++)
+  for (i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
-
       if (side == 0) {
         if (x[i][dim] < coord) {
           x[i][dim] = coord + (coord - x[i][dim]);
@@ -247,4 +239,5 @@ void FixWallReflect::wall_particle(int m, int which, double coord)
         }
       }
     }
+  }
 }

src/fix_wall_reflect.h

@@ -26,14 +26,14 @@ namespace LAMMPS_NS {
 
 class FixWallReflect : public Fix {
  public:
+  enum{XLO=0,XHI=1,YLO=2,YHI=3,ZLO=4,ZHI=5};
+  enum{NONE=0,EDGE,CONSTANT,VARIABLE};
+
   FixWallReflect(class LAMMPS *, int, char **);
   virtual ~FixWallReflect();
-  virtual void wall_particle(int m, int which, double coord);
   int setmask();
   void init();
   void post_integrate();
-  enum{XLO=0,XHI=1,YLO=2,YHI=3,ZLO=4,ZHI=5};
-  enum{NONE=0,EDGE,CONSTANT,VARIABLE};
   
  protected:
   int nwall;
@@ -43,6 +43,8 @@ class FixWallReflect : public Fix {
   int varindex[6];
   int varflag;
   double xscale,yscale,zscale;
+
+  virtual void wall_particle(int m, int which, double coord);
 };
 
 }

src/random_mars.cpp

@@ -90,7 +90,6 @@ double RanMars::uniform()
   return uni;
 }
 
-
 /* ----------------------------------------------------------------------
    gaussian RN
 ------------------------------------------------------------------------- */
@@ -127,7 +126,6 @@ double RanMars::gaussian(double mu, double sigma)
   return v1;
 }
 
-
 /* ----------------------------------------------------------------------
    Rayleigh RN
 ------------------------------------------------------------------------- */
@@ -136,16 +134,13 @@ double RanMars::rayleigh(double sigma)
 {
   double first,v1;
  
-  if (sigma <= 0) {
+  if (sigma <= 0) error->all(FLERR,"Invalid Rayleigh parameter");
-    error->all(FLERR,"Invalid Rayleigh parameter");
+
-  } else {
   v1 = uniform();
   first = sigma*sqrt(-2.0*log(v1));
   return first;
-  }
 }
 
-
 /* ----------------------------------------------------------------------
    Bessel exponential RN
 ------------------------------------------------------------------------- */
@@ -153,20 +148,21 @@ double RanMars::rayleigh(double sigma)
 double RanMars::besselexp(double theta, double alpha, double cp)
 {
   double first,v1,v2;
-  if ((theta < 0) || (alpha < 0) || (alpha >1)) {
+
+  if (theta < 0.0 || alpha < 0.0 || alpha < 1.0)
     error->all(FLERR,"Invalid Bessel exponential distribution parameters");
-  } else {
+  
   v1 = uniform();
   v2 = uniform();
-    if (cp < 0) {
+
-      first = sqrt((1.0-alpha)*cp*cp - 2.0*alpha*theta*log(v1)
+  if (cp < 0.0)
-                   + 2.0*sqrt(-2.0*theta*(1.0-alpha)*alpha*log(v1))
+    first = sqrt((1.0-alpha)*cp*cp - 2.0*alpha*theta*log(v1) + 
-                   * cos(2.0*MathConst::MY_PI*v2)*cp);
+                 2.0*sqrt(-2.0*theta*(1.0-alpha)*alpha*log(v1)) * 
-    } else {
+                 cos(2.0*MathConst::MY_PI*v2)*cp);
-      first = -sqrt((1.0-alpha)*cp*cp - 2.0*alpha*theta*log(v1)
+  else
-                    - 2.0*sqrt(-2.0*theta*(1.0-alpha)*alpha*log(v1))
+    first = -sqrt((1.0-alpha)*cp*cp - 2.0*alpha*theta*log(v1) - 
-                    * cos(2.0*MathConst::MY_PI*v2)*cp);
+                  2.0*sqrt(-2.0*theta*(1.0-alpha)*alpha*log(v1)) * 
-    }
+                  cos(2.0*MathConst::MY_PI*v2)*cp);
+
   return first;
-  }
 }