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

src/COLLOID/fix_wall_colloid.cpp

@@ -16,14 +16,11 @@
 ------------------------------------------------------------------------- */
 
 #include "math.h"
-#include "stdlib.h"
 #include "string.h"
 #include "fix_wall_colloid.h"
 #include "atom.h"
 #include "atom_vec.h"
-#include "domain.h"
 #include "update.h"
-#include "output.h"
 #include "respa.h"
 #include "error.h"
 
@@ -32,74 +29,7 @@ using namespace LAMMPS_NS;
 /* ---------------------------------------------------------------------- */
 
 FixWallColloid::FixWallColloid(LAMMPS *lmp, int narg, char **arg) : 
-  Fix(lmp, narg, arg)
+  FixWall(lmp, narg, arg) {}
-{
-  if (narg != 8) error->all("Illegal fix wall/colloid command");
-
-  scalar_flag = 1;
-  vector_flag = 1;
-  size_vector = 3;
-  scalar_vector_freq = 1;
-  extscalar = 1;
-  extvector = 1;
-
-  if (strcmp(arg[3],"xlo") == 0) {
-    dim = 0;
-    side = -1;
-  } else if (strcmp(arg[3],"xhi") == 0) {
-    dim = 0;
-    side = 1;
-  } else if (strcmp(arg[3],"ylo") == 0) {
-    dim = 1;
-    side = -1;
-  } else if (strcmp(arg[3],"yhi") == 0) {
-    dim = 1;
-    side = 1;
-  } else if (strcmp(arg[3],"zlo") == 0) {
-    dim = 2;
-    side = -1;
-  } else if (strcmp(arg[3],"zhi") == 0) {
-    dim = 2;
-    side = 1;
-  } else error->all("Illegal fix wall/colloid command");
-
-  coord = atof(arg[4]);
-  hamaker = atof(arg[5]);
-  sigma = atof(arg[6]);
-  cutoff = atof(arg[7]);
-
-  coeff1 = -4.0/315.0 * hamaker * pow(sigma,6.0);
-  coeff2 = -2.0/3.0 * hamaker;
-  coeff3 = hamaker * pow(sigma,6.0)/7560.0;
-  coeff4 = hamaker/6.0;
-
-  double rinv = 1.0/cutoff;
-  double r2inv = rinv*rinv;
-  double r4inv = r2inv*r2inv;
-  offset = coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv;
-
-  if (dim == 0 && domain->xperiodic)
-    error->all("Cannot use wall in periodic dimension");
-  if (dim == 1 && domain->yperiodic)
-    error->all("Cannot use wall in periodic dimension");
-  if (dim == 2 && domain->zperiodic)
-    error->all("Cannot use wall in periodic dimension");
-
-  wall_flag = 0;
-  wall[0] = wall[1] = wall[2] = wall[3] = 0.0;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int FixWallColloid::setmask()
-{
-  int mask = 0;
-  mask |= POST_FORCE;
-  mask |= THERMO_ENERGY;
-  mask |= POST_FORCE_RESPA;
-  mask |= MIN_POST_FORCE;
-  return mask;
-}
 
 /* ---------------------------------------------------------------------- */
 
@@ -135,34 +65,33 @@ void FixWallColloid::init()
   MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
   if (flagall) error->all("Fix wall/colloid requires extended particles");
 
-  if (strcmp(update->integrate_style,"respa") == 0)
+  FixWall::init();
-    nlevels_respa = ((Respa *) update->integrate)->nlevels;
 }
 
 /* ---------------------------------------------------------------------- */
 
-void FixWallColloid::setup(int vflag)
+void FixWallColloid::precompute(int m)
 {
-  if (strcmp(update->integrate_style,"verlet") == 0)
+  coeff1[m] = -4.0/315.0 * epsilon[m] * pow(sigma[m],6.0);
-    post_force(vflag);
+  coeff2[m] = -2.0/3.0 * epsilon[m];
-  else {
+  coeff3[m] = epsilon[m] * pow(sigma[m],6.0)/7560.0;
-    ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
+  coeff4[m] = epsilon[m]/6.0;
-    post_force_respa(vflag,nlevels_respa-1,0);
+
-    ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+  double rinv = 1.0/cutoff[m];
-  }
+  double r2inv = rinv*rinv;
+  double r4inv = r2inv*r2inv;
+  offset[m] = coeff3[m]*r4inv*r4inv*rinv - coeff4[m]*r2inv*rinv;
 }
 
 /* ---------------------------------------------------------------------- */
 
-void FixWallColloid::min_setup(int vflag)
+void FixWallColloid::wall_particle(int m, double coord)
 {
-  post_force(vflag);
+  double delta,delta2,rinv,r2inv,r4inv,r8inv,fwall;
-}
+  double r2,rinv2,r2inv2,r4inv2,r6inv2;
+  double r3,rinv3,r2inv3,r4inv3,r6inv3;
+  double rad,rad2,rad4,rad8,diam,new_coeff2;
 
-/* ---------------------------------------------------------------------- */
-
-void FixWallColloid::post_force(int vflag)
-{
   double **x = atom->x;
   double **f = atom->f;
   double **shape = atom->shape;
@@ -170,21 +99,17 @@ void FixWallColloid::post_force(int vflag)
   int *type = atom->type;
   int nlocal = atom->nlocal;
 
-  double delta,delta2,rinv,r2inv,r4inv,r8inv,fwall;
+  int dim = m/2;
-  double r2,rinv2,r2inv2,r4inv2,r6inv2;
+  int side = m % 2;
-  double r3,rinv3,r2inv3,r4inv3,r6inv3;
-  double rad,rad2,rad4,rad8,diam,new_coeff2;
-  wall[0] = wall[1] = wall[2] = wall[3] = 0.0;
-  wall_flag = 0;
 
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
-      if (side == -1) delta = x[i][dim] - coord;
+      if (side == 0) delta = x[i][dim] - coord;
       else delta = coord - x[i][dim];
       if (delta <= 0.0) continue;
-      if (delta > cutoff) continue;
+      if (delta > cutoff[m]) continue;
       rad = shape[type[i]][0];
-      new_coeff2 = coeff2*rad*rad*rad;
+      new_coeff2 = coeff2[m]*rad*rad*rad;
       diam = 2.0*rad;
       rad2 = rad*rad;
       rad4 = rad2*rad2;
@@ -194,7 +119,7 @@ void FixWallColloid::post_force(int vflag)
       r2inv = rinv*rinv;
       r4inv = r2inv*r2inv;
       r8inv = r4inv*r4inv;
-      fwall = (coeff1*(rad8*rad + 27.0*rad4*rad2*rad*pow(delta,2.0)
+      fwall = (coeff1[m]*(rad8*rad + 27.0*rad4*rad2*rad*pow(delta,2.0)
 			  + 63.0*rad4*rad*pow(delta,4.0)
 			  + 21.0*rad2*rad*pow(delta,6.0))*r8inv - 
 	       new_coeff2*r2inv) * side;
@@ -204,59 +129,15 @@ void FixWallColloid::post_force(int vflag)
       r2inv2 = rinv2*rinv2;
       r4inv2 = r2inv2*r2inv2;
       r6inv2 = r4inv2*r2inv2;
-      r3 = delta+0.5*diam;
+      r3 = delta + 0.5*diam;
       rinv3 = 1.0/r3;
       r2inv3 = rinv3*rinv3;
       r4inv3 = r2inv3*r2inv3;
       r6inv3 = r4inv3*r2inv3;
-      wall[0] += coeff3*((-3.5*diam+delta)*r4inv2*r2inv2*rinv2
+      ewall[0] += coeff3[m]*((-3.5*diam+delta)*r4inv2*r2inv2*rinv2
 			     + (3.5*diam+delta)*r4inv3*r2inv3*rinv3) -
-	coeff4*((-diam*delta+r2*r3*(log(-r2)-log(r3)))*
+	coeff4[m]*((-diam*delta+r2*r3*(log(-r2)-log(r3)))*
-		(-rinv2)*rinv3) - offset;
+		   (-rinv2)*rinv3) - offset[m];
-      wall[dim+1] += fwall;
+      ewall[m+1] += fwall;
     }
 }
-
-/* ---------------------------------------------------------------------- */
-
-void FixWallColloid::post_force_respa(int vflag, int ilevel, int iloop)
-{
-  if (ilevel == nlevels_respa-1) post_force(vflag);
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixWallColloid::min_post_force(int vflag)
-{
-  post_force(vflag);
-}
-
-/* ----------------------------------------------------------------------
-   energy of wall interaction
-------------------------------------------------------------------------- */
-
-double FixWallColloid::compute_scalar()
-{
-  // only sum across procs one time
-
-  if (wall_flag == 0) {
-    MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world);
-    wall_flag = 1;
-  }
-  return wall_all[0];
-}
-
-/* ----------------------------------------------------------------------
-   components of force on wall
-------------------------------------------------------------------------- */
-
-double FixWallColloid::compute_vector(int n)
-{
-  // only sum across procs one time
-
-  if (wall_flag == 0) {
-    MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world);
-    wall_flag = 1;
-  }
-  return wall_all[n+1];
-}

src/COLLOID/fix_wall_colloid.h

@@ -14,31 +14,19 @@
 #ifndef FIX_WALL_COLLOID_H
 #define FIX_WALL_COLLOID_H
 
-#include "fix.h"
+#include "fix_wall.h"
 
 namespace LAMMPS_NS {
 
-class FixWallColloid : public Fix {
+class FixWallColloid : public FixWall {
  public:
   FixWallColloid(class LAMMPS *, int, char **);
-  ~FixWallColloid() {}
-  int setmask();
   void init();
-  void setup(int);
+  void precompute(int);
-  void min_setup(int);
+  void wall_particle(int, double);
-  void post_force(int);
-  void post_force_respa(int, int, int);
-  void min_post_force(int);
-  double compute_scalar();
-  double compute_vector(int);
 
  private:
-  int dim,side,thermo_flag,eflag_enable;
+  double coeff1[6],coeff2[6],coeff3[6],coeff4[6],offset[6];
-  double coord,hamaker,sigma,cutoff;
-  double coeff1,coeff2,coeff3,coeff4,offset;
-  double wall[4],wall_all[4];
-  int wall_flag;
-  int nlevels_respa;
 };
 
 }

src/fix_wall.cpp

@@ -0,0 +1,222 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "stdlib.h"
+#include "string.h"
+#include "fix_wall.h"
+#include "atom.h"
+#include "domain.h"
+#include "update.h"
+#include "output.h"
+#include "respa.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+enum{XLO,XHI,YLO,YHI,ZLO,ZHI};
+
+/* ---------------------------------------------------------------------- */
+
+FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  scalar_flag = 1;
+  vector_flag = 1;
+  size_vector = 6;
+  scalar_vector_freq = 1;
+  extscalar = 1;
+  extvector = 1;
+
+  // parse args
+
+  for (int m = 0; m < 6; m++) wflag[m] = 0;
+  vflag = 0;
+  aflag = 0;
+
+  int iarg = 3;
+  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+5 > narg) error->all("Illegal fix wall command");
+    int m;
+    if (strcmp(arg[iarg],"xlo") == 0) m = XLO;
+    else if (strcmp(arg[iarg],"xhi") == 0) m = XHI;
+    else if (strcmp(arg[iarg],"ylo") == 0) m = YLO;
+    else if (strcmp(arg[iarg],"yhi") == 0) m = YHI;
+    else if (strcmp(arg[iarg],"zlo") == 0) m = ZLO;
+    else if (strcmp(arg[iarg],"zhi") == 0) m = ZHI;
+    wflag[m] = 1;
+    coord0[m] = atof(arg[iarg+1]);
+    epsilon[m] = atof(arg[iarg+2]);
+    sigma[m] = atof(arg[iarg+3]);
+    cutoff[m] = atof(arg[iarg+4]);
+    iarg += 5;
+  } else if (strcmp(arg[iarg],"vel") == 0) {
+    if (iarg+2 > narg) error->all("Illegal fix wall command");
+    vflag = 1;
+    vel = atof(arg[iarg+1]);
+    iarg += 2;
+  } else if (strcmp(arg[iarg],"wiggle") == 0) {
+    if (iarg+3 > narg) error->all("Illegal fix wall command");
+    aflag = 1;
+    amplitude = atof(arg[iarg+1]);
+    period = atof(arg[iarg+2]);
+    iarg += 3;
+  } else error->all("Illegal fix wall command");
+
+  // error check
+
+  int flag = 0;
+  for (int m = 0; m < 6; m++) if (wflag[m]) flag;
+  if (!flag) error->all("Illegal fix wall command");
+
+  if (vflag && aflag) 
+    error->all("Cannot set both vel and wiggle in fix wall command");
+
+  if ((wflag[XLO] || wflag[XHI]) && domain->xperiodic)
+    error->all("Cannot use fix wall in periodic dimension");
+  if ((wflag[YLO] || wflag[YHI]) && domain->yperiodic)
+    error->all("Cannot use fix wall in periodic dimension");
+  if ((wflag[ZLO] || wflag[ZHI]) && domain->xperiodic)
+    error->all("Cannot use fix wall in periodic zimension");
+
+  // setup oscillations
+
+  if (aflag) {
+    double PI = 4.0 * atan(1.0);
+    omega = 2.0*PI / period;
+  }
+
+  eflag = 0;
+  for (int m = 0; m < 7; m++) ewall[m] = 0.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixWall::setmask()
+{
+  int mask = 0;
+  mask |= POST_FORCE;
+  mask |= THERMO_ENERGY;
+  mask |= POST_FORCE_RESPA;
+  mask |= MIN_POST_FORCE;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWall::init()
+{
+  dt = update->dt;
+
+  // setup coefficients
+
+  for (int m = 0; m < 6; m++)
+    if (wflag[m]) precompute(m);
+
+  if (strcmp(update->integrate_style,"respa") == 0)
+    nlevels_respa = ((Respa *) update->integrate)->nlevels;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWall::setup(int vflag)
+{
+  if (strcmp(update->integrate_style,"verlet") == 0)
+    post_force(vflag);
+  else {
+    ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
+    post_force_respa(vflag,nlevels_respa-1,0);
+    ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWall::min_setup(int vflag)
+{
+  post_force(vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWall::post_force(int vflag)
+{
+  eflag = 0;
+  for (int m = 0; m < 7; m++) ewall[m] = 0.0;
+
+  // coord0 = initial position of wall
+  // coord = current position of wall
+
+  double delta = (update->ntimestep - update->beginstep) * dt;
+  double coord;
+
+  for (int m = 0; m < 6; m++) {
+    if (wflag[m] == 0) continue;
+
+    if (vflag) {
+      if (m/2 == 0) coord = coord0[m] + delta*vel;
+      else coord = coord0[m] - delta*vel;
+    } else if (aflag) {
+      if (m/2 == 0) coord = coord0[m] + amplitude*sin(omega*delta);
+      else coord = coord0[m] - amplitude*sin(omega*delta);
+    } else coord = coord0[m];
+
+    wall_particle(m,coord);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWall::post_force_respa(int vflag, int ilevel, int iloop)
+{
+  if (ilevel == nlevels_respa-1) post_force(vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWall::min_post_force(int vflag)
+{
+  post_force(vflag);
+}
+
+/* ----------------------------------------------------------------------
+   energy of wall interaction
+------------------------------------------------------------------------- */
+
+double FixWall::compute_scalar()
+{
+  // only sum across procs one time
+
+  if (eflag == 0) {
+    MPI_Allreduce(ewall,ewall_all,7,MPI_DOUBLE,MPI_SUM,world);
+    eflag = 1;
+  }
+  return ewall_all[0];
+}
+
+/* ----------------------------------------------------------------------
+   components of force on wall
+------------------------------------------------------------------------- */
+
+double FixWall::compute_vector(int n)
+{
+  // only sum across procs one time
+
+  if (eflag == 0) {
+    MPI_Allreduce(ewall,ewall_all,7,MPI_DOUBLE,MPI_SUM,world);
+    eflag = 1;
+  }
+  return ewall_all[n+1];
+}

src/fix_wall.h

@@ -0,0 +1,51 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef FIX_WALL_H
+#define FIX_WALL_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixWall : public Fix {
+ public:
+  FixWall(class LAMMPS *, int, char **);
+  virtual ~FixWall() {}
+  int setmask();
+  virtual void init();
+  void setup(int);
+  void min_setup(int);
+  void post_force(int);
+  void post_force_respa(int, int, int);
+  void min_post_force(int);
+  double compute_scalar();
+  double compute_vector(int);
+
+  virtual void precompute(int) = 0;
+  virtual void wall_particle(int, double) = 0;
+
+ protected:
+  int wflag[6];
+  double coord0[6],epsilon[6],sigma[6],cutoff[6];
+  int vflag,aflag;
+  double vel,amplitude,period,omega;
+  int eflag;
+  double ewall[7],ewall_all[7];
+  int nlevels_respa;
+  double dt;
+};
+
+}
+
+#endif

src/fix_wall_lj126.cpp

@@ -11,212 +11,57 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
-/* ----------------------------------------------------------------------
-   Contributing author: Mark Stevens (SNL)
-------------------------------------------------------------------------- */
-
 #include "math.h"
-#include "stdlib.h"
-#include "string.h"
 #include "fix_wall_lj126.h"
 #include "atom.h"
-#include "domain.h"
-#include "update.h"
-#include "output.h"
-#include "respa.h"
-#include "error.h"
 
 using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
 FixWallLJ126::FixWallLJ126(LAMMPS *lmp, int narg, char **arg) : 
-  Fix(lmp, narg, arg)
+  FixWall(lmp, narg, arg) {}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallLJ126::precompute(int m)
 {
-  if (narg < 8) error->all("Illegal fix wall/lj126 command");
+  coeff1[m] = 48.0 * epsilon[m] * pow(sigma[m],12.0);
+  coeff2[m] = 24.0 * epsilon[m] * pow(sigma[m],6.0);
+  coeff3[m] = 4.0 * epsilon[m] * pow(sigma[m],12.0);
+  coeff4[m] = 4.0 * epsilon[m] * pow(sigma[m],6.0);
 
-  scalar_flag = 1;
+  double r2inv = 1.0/(cutoff[m]*cutoff[m]);
-  vector_flag = 1;
-  size_vector = 3;
-  scalar_vector_freq = 1;
-  extscalar = 1;
-  extvector = 1;
-
-  // set defaults
-
-  if (strcmp(arg[3],"xlo") == 0) {
-    dim = 0;
-    side = -1;
-  } else if (strcmp(arg[3],"xhi") == 0) {
-    dim = 0;
-    side = 1;
-  } else if (strcmp(arg[3],"ylo") == 0) {
-    dim = 1;
-    side = -1;
-  } else if (strcmp(arg[3],"yhi") == 0) {
-    dim = 1;
-    side = 1;
-  } else if (strcmp(arg[3],"zlo") == 0) {
-    dim = 2;
-    side = -1;
-  } else if (strcmp(arg[3],"zhi") == 0) {
-    dim = 2;
-    side = 1;
-  } else error->all("Illegal fix wall/lj126 command");
-
-  coord0 = atof(arg[4]);
-  epsilon = atof(arg[5]);
-  sigma = atof(arg[6]);
-  cutoff = atof(arg[7]);
-
-  // read options
-
-  vel = 0.0;
-
-  int iarg = 8;
-  while (iarg < narg) {
-    if (strcmp(arg[iarg],"vel") == 0) {
-      if (iarg+2 > narg) error->all("Illegal fix wall/lj126 command");
-      vel = atof(arg[iarg+1]);
-      iarg += 2;
-    } else error->all("Illegal fix wall/lj126 command");
-  }
-
-  coeff1 = 48.0 * epsilon * pow(sigma,12.0);
-  coeff2 = 24.0 * epsilon * pow(sigma,6.0);
-  coeff3 = 4.0 * epsilon * pow(sigma,12.0);
-  coeff4 = 4.0 * epsilon * pow(sigma,6.0);
-
-  double r2inv = 1.0/(cutoff*cutoff);
   double r6inv = r2inv*r2inv*r2inv;
-  offset = r6inv*(coeff3*r6inv - coeff4);
+  offset[m] = r6inv*(coeff3[m]*r6inv - coeff4[m]);
-
-  if (dim == 0 && domain->xperiodic)
-    error->all("Cannot use wall in periodic dimension");
-  if (dim == 1 && domain->yperiodic)
-    error->all("Cannot use wall in periodic dimension");
-  if (dim == 2 && domain->zperiodic)
-    error->all("Cannot use wall in periodic dimension");
-
-  wall_flag = 0;
-  wall[0] = wall[1] = wall[2] = wall[3] = 0.0;
 }
 
 /* ---------------------------------------------------------------------- */
 
-int FixWallLJ126::setmask()
+void FixWallLJ126::wall_particle(int m, double coord)
 {
-  int mask = 0;
+  double delta,rinv,r2inv,r6inv,fwall;
-  mask |= POST_FORCE;
-  mask |= THERMO_ENERGY;
-  mask |= POST_FORCE_RESPA;
-  mask |= MIN_POST_FORCE;
-  return mask;
-}
 
-/* ---------------------------------------------------------------------- */
-
-void FixWallLJ126::init()
-{
-  if (strcmp(update->integrate_style,"respa") == 0)
-    nlevels_respa = ((Respa *) update->integrate)->nlevels;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixWallLJ126::setup(int vflag)
-{
-  if (strcmp(update->integrate_style,"verlet") == 0)
-    post_force(vflag);
-  else {
-    ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
-    post_force_respa(vflag,nlevels_respa-1,0);
-    ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixWallLJ126::min_setup(int vflag)
-{
-  post_force(vflag);
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixWallLJ126::post_force(int vflag)
-{
   double **x = atom->x;
   double **f = atom->f;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  double delta,rinv,r2inv,r6inv,fwall;
+  int dim = m/2;
-  wall[0] = wall[1] = wall[2] = wall[3] = 0.0;
+  int side = m % 2;
-  wall_flag = 0;
-
-  // coord = current position of wall
-  // coord0 = initial position of wall
-  
-  double delt = (update->ntimestep - update->beginstep) * update->dt;
-  double coord = coord0 + delt*vel;
 
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
-      if (side == -1) delta = x[i][dim] - coord;
+      if (side == 0) delta = x[i][dim] - coord;
       else delta = coord - x[i][dim];
       if (delta <= 0.0) continue;
-      if (delta > cutoff) continue;
+      if (delta > cutoff[m]) continue;
       rinv = 1.0/delta;
       r2inv = rinv*rinv;
       r6inv = r2inv*r2inv*r2inv;
-      fwall = side * r6inv*(coeff1*r6inv - coeff2) * rinv;
+      fwall = side * r6inv*(coeff1[m]*r6inv - coeff2[m]) * rinv;
       f[i][dim] -= fwall;
-      wall[0] += r6inv*(coeff3*r6inv - coeff4) - offset;
+      ewall[0] += r6inv*(coeff3[m]*r6inv - coeff4[m]) - offset[m];
-      wall[dim+1] += fwall;
+      ewall[m+1] += fwall;
     }
 }
-
-/* ---------------------------------------------------------------------- */
-
-void FixWallLJ126::post_force_respa(int vflag, int ilevel, int iloop)
-{
-  if (ilevel == nlevels_respa-1) post_force(vflag);
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixWallLJ126::min_post_force(int vflag)
-{
-  post_force(vflag);
-}
-
-/* ----------------------------------------------------------------------
-   energy of wall interaction
-------------------------------------------------------------------------- */
-
-double FixWallLJ126::compute_scalar()
-{
-  // only sum across procs one time
-
-  if (wall_flag == 0) {
-    MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world);
-    wall_flag = 1;
-  }
-  return wall_all[0];
-}
-
-/* ----------------------------------------------------------------------
-   components of force on wall
-------------------------------------------------------------------------- */
-
-double FixWallLJ126::compute_vector(int n)
-{
-  // only sum across procs one time
-
-  if (wall_flag == 0) {
-    MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world);
-    wall_flag = 1;
-  }
-  return wall_all[n+1];
-}

src/fix_wall_lj126.h

@@ -14,34 +14,19 @@
 #ifndef FIX_WALL_LJ126_H
 #define FIX_WALL_LJ126_H
 
-#include "fix.h"
+#include "fix_wall.h"
 
 namespace LAMMPS_NS {
 
-class FixWallLJ126 : public Fix {
+class FixWallLJ126 : public FixWall {
  public:
   FixWallLJ126(class LAMMPS *, int, char **);
-  int setmask();
+  void precompute(int);
-  void init();
+  void wall_particle(int, double);
-  void setup(int);
-  void min_setup(int);
-  void post_force(int);
-  void post_force_respa(int, int, int);
-  void min_post_force(int);
-  double compute_scalar();
-  double compute_vector(int);
 
  private:
-  int dim,side,thermo_flag,eflag_enable;
+  double coeff1[6],coeff2[6],coeff3[6],coeff4[6],offset[6];
-  double coord,epsilon,sigma,cutoff;
-  double coeff1,coeff2,coeff3,coeff4,offset;
-  double wall[4],wall_all[4];
-  double vel,coord0;
-  int wall_flag;
-  int nlevels_respa;
 };
 
 #endif
 }
-
-

src/fix_wall_lj93.cpp

@@ -12,209 +12,59 @@
 ------------------------------------------------------------------------- */
 
 #include "math.h"
-#include "stdlib.h"
-#include "string.h"
 #include "fix_wall_lj93.h"
 #include "atom.h"
-#include "domain.h"
-#include "update.h"
-#include "output.h"
-#include "respa.h"
-#include "error.h"
 
 using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
 FixWallLJ93::FixWallLJ93(LAMMPS *lmp, int narg, char **arg) : 
-  Fix(lmp, narg, arg)
+  FixWall(lmp, narg, arg) {}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallLJ93::precompute(int m)
 {
-  if (narg < 8) error->all("Illegal fix wall/lj93 command");
+  coeff1[m] = 6.0/5.0 * epsilon[m] * pow(sigma[m],9.0);
+  coeff2[m] = 3.0 * epsilon[m] * pow(sigma[m],3.0);
+  coeff3[m] = 2.0/15.0 * epsilon[m] * pow(sigma[m],9.0);
+  coeff4[m] = epsilon[m] * pow(sigma[m],3.0);
 
-  scalar_flag = 1;
+  double rinv = 1.0/cutoff[m];
-  vector_flag = 1;
-  size_vector = 3;
-  scalar_vector_freq = 1;
-  extscalar = 1;
-  extvector = 1;
-
-  // set defaults
-
-  if (strcmp(arg[3],"xlo") == 0) {
-    dim = 0;
-    side = -1;
-  } else if (strcmp(arg[3],"xhi") == 0) {
-    dim = 0;
-    side = 1;
-  } else if (strcmp(arg[3],"ylo") == 0) {
-    dim = 1;
-    side = -1;
-  } else if (strcmp(arg[3],"yhi") == 0) {
-    dim = 1;
-    side = 1;
-  } else if (strcmp(arg[3],"zlo") == 0) {
-    dim = 2;
-    side = -1;
-  } else if (strcmp(arg[3],"zhi") == 0) {
-    dim = 2;
-    side = 1;
-  } else error->all("Illegal fix wall/lj93 command");
-
-  coord0 = atof(arg[4]);
-  epsilon = atof(arg[5]);
-  sigma = atof(arg[6]);
-  cutoff = atof(arg[7]);
-
-  // read options 
-
-  vel = 0.0;
-
-  int iarg = 8;
-  while (iarg < narg) {
-    if (strcmp(arg[iarg],"vel") == 0) {
-      if (iarg+2 > narg) error->all("Illegal fix wall/lj93 command");
-      vel = atof(arg[iarg+1]);
-      iarg += 2;
-    } else error->all("Illegal fix wall/lj93 command");
-  }
-
-  coeff1 = 6.0/5.0 * epsilon * pow(sigma,9.0);
-  coeff2 = 3.0 * epsilon * pow(sigma,3.0);
-  coeff3 = 2.0/15.0 * epsilon * pow(sigma,9.0);
-  coeff4 = epsilon * pow(sigma,3.0);
-
-  double rinv = 1.0/cutoff;
   double r2inv = rinv*rinv;
   double r4inv = r2inv*r2inv;
-  offset = coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv;
+  offset[m] = coeff3[6]*r4inv*r4inv*rinv - coeff4[6]*r2inv*rinv;
-
-  if (dim == 0 && domain->xperiodic)
-    error->all("Cannot use wall in periodic dimension");
-  if (dim == 1 && domain->yperiodic)
-    error->all("Cannot use wall in periodic dimension");
-  if (dim == 2 && domain->zperiodic)
-    error->all("Cannot use wall in periodic dimension");
-
-  wall_flag = 0;
-  wall[0] = wall[1] = wall[2] = wall[3] = 0.0;
 }
 
 /* ---------------------------------------------------------------------- */
 
-int FixWallLJ93::setmask()
+void FixWallLJ93::wall_particle(int m, double coord)
 {
-  int mask = 0;
+  double delta,rinv,r2inv,r4inv,r10inv,fwall;
-  mask |= POST_FORCE;
-  mask |= THERMO_ENERGY;
-  mask |= POST_FORCE_RESPA;
-  mask |= MIN_POST_FORCE;
-  return mask;
-}
 
-/* ---------------------------------------------------------------------- */
-
-void FixWallLJ93::init()
-{
-  if (strcmp(update->integrate_style,"respa") == 0)
-    nlevels_respa = ((Respa *) update->integrate)->nlevels;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixWallLJ93::setup(int vflag)
-{
-  if (strcmp(update->integrate_style,"verlet") == 0)
-    post_force(vflag);
-  else {
-    ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
-    post_force_respa(vflag,nlevels_respa-1,0);
-    ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixWallLJ93::min_setup(int vflag)
-{
-  post_force(vflag);
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixWallLJ93::post_force(int vflag)
-{
   double **x = atom->x;
   double **f = atom->f;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  double delta,rinv,r2inv,r4inv,r10inv,fwall;
+  int dim = m/2;
-  wall[0] = wall[1] = wall[2] = wall[3] = 0.0;
+  int side = m % 2;
-  wall_flag = 0;
-
-  // coord = current position of wall
-  // coord0 = initial position of wall
-
-  double delt = (update->ntimestep - update->beginstep) * update->dt;
-  double coord = coord0 + delt*vel;
 
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
-      if (side == -1) delta = x[i][dim] - coord;
+      if (side == 0) delta = x[i][dim] - coord;
       else delta = coord - x[i][dim];
       if (delta <= 0.0) continue;
-      if (delta > cutoff) continue;
+      if (delta > cutoff[m]) continue;
       rinv = 1.0/delta;
       r2inv = rinv*rinv;
       r4inv = r2inv*r2inv;
       r10inv = r4inv*r4inv*r2inv;
-      fwall = (coeff1*r10inv - coeff2*r4inv) * side;
+      fwall = (coeff1[m]*r10inv - coeff2[m]*r4inv) * side;
       f[i][dim] -= fwall;
-      wall[0] += coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv - offset;
+      ewall[0] += coeff3[m]*r4inv*r4inv*rinv - 
-      wall[dim+1] += fwall;
+	coeff4[m]*r2inv*rinv - offset[m];
+      ewall[m+1] += fwall;
     }
 }
-
-/* ---------------------------------------------------------------------- */
-
-void FixWallLJ93::post_force_respa(int vflag, int ilevel, int iloop)
-{
-  if (ilevel == nlevels_respa-1) post_force(vflag);
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixWallLJ93::min_post_force(int vflag)
-{
-  post_force(vflag);
-}
-
-/* ----------------------------------------------------------------------
-   energy of wall interaction
-------------------------------------------------------------------------- */
-
-double FixWallLJ93::compute_scalar()
-{
-  // only sum across procs one time
-
-  if (wall_flag == 0) {
-    MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world);
-    wall_flag = 1;
-  }
-  return wall_all[0];
-}
-
-/* ----------------------------------------------------------------------
-   components of force on wall
-------------------------------------------------------------------------- */
-
-double FixWallLJ93::compute_vector(int n)
-{
-  // only sum across procs one time
-
-  if (wall_flag == 0) {
-    MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world);
-    wall_flag = 1;
-  }
-  return wall_all[n+1];
-}

src/fix_wall_lj93.h

@@ -14,32 +14,18 @@
 #ifndef FIX_WALL_LJ93_H
 #define FIX_WALL_LJ93_H
 
-#include "fix.h"
+#include "fix_wall.h"
 
 namespace LAMMPS_NS {
 
-class FixWallLJ93 : public Fix {
+class FixWallLJ93 : public FixWall {
  public:
   FixWallLJ93(class LAMMPS *, int, char **);
-  ~FixWallLJ93() {}
+  void precompute(int);
-  int setmask();
+  void wall_particle(int, double);
-  void init();
-  void setup(int);
-  void min_setup(int);
-  void post_force(int);
-  void post_force_respa(int, int, int);
-  void min_post_force(int);
-  double compute_scalar();
-  double compute_vector(int);
 
  private:
-  int dim,side,thermo_flag,eflag_enable;
+  double coeff1[6],coeff2[6],coeff3[6],coeff4[6],offset[6];
-  double coord,epsilon,sigma,cutoff;
-  double coeff1,coeff2,coeff3,coeff4,offset;
-  double wall[4],wall_all[4];
-  double vel,coord0;
-  int wall_flag;
-  int nlevels_respa;
 };
 
 }