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

src/COLLOID/fix_wall_colloid.cpp

@@ -33,7 +33,7 @@ using namespace LAMMPS_NS;
 FixWallColloid::FixWallColloid(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg)
 {
-  if (narg != 9) error->all("Illegal fix wall/colloid command");
+  if (narg != 8) error->all("Illegal fix wall/colloid command");
 
   scalar_flag = 1;
   vector_flag = 1;
@@ -63,11 +63,12 @@ FixWallColloid::FixWallColloid(LAMMPS *lmp, int narg, char **arg) :
   } else error->all("Illegal fix wall/colloid command");
 
   coord = atof(arg[4]);
+  //NOTE: this next variable should become Hamaker pre-factor
   epsilon = atof(arg[5]);
   sigma = atof(arg[6]);
-  diam = atof(arg[7]);
-  cutoff = atof(arg[8]);
+  cutoff = atof(arg[7]);
 
+  //NOTE: coeff2 uses diam, so will need to be computed below for each particle?
   coeff1 = -576.0/315.0 * epsilon * pow(sigma,6.0);
   coeff2 = -288.0/3.0 * 0.125*diam*diam*diam* epsilon;
   coeff3 = 144.0 * epsilon * pow(sigma,6.0)/7560.0;
@@ -105,6 +106,17 @@ int FixWallColloid::setmask()
 
 void FixWallColloid::init()
 {
+  if (!atom->shape)
+    error->all("Fix wall/colloid requires atom attribute shape");
+
+  // insure all particle shapes are spherical
+
+  for (int i = 1; i <= atom->ntypes; i++)
+    if ((atom->shape[i][0] != atom->shape[i][1]) || 
+	(atom->shape[i][0] != atom->shape[i][2]) ||
+	(atom->shape[i][1] != atom->shape[i][2]))
+      error->all("Fix wall/colloid requires spherical particles");
+
   if (strcmp(update->integrate_style,"respa") == 0)
     nlevels_respa = ((Respa *) update->integrate)->nlevels;
 }
@@ -136,12 +148,13 @@ void FixWallColloid::post_force(int vflag)
   double **x = atom->x;
   double **f = atom->f;
   int *mask = atom->mask;
+  int *type = atom->type;
   int nlocal = atom->nlocal;
 
   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;
+  double rad,rad2,rad4,rad8,diam;
   wall[0] = wall[1] = wall[2] = wall[3] = 0.0;
   wall_flag = 0;
 
@@ -151,7 +164,8 @@ void FixWallColloid::post_force(int vflag)
       else delta = coord - x[i][dim];
       if (delta <= 0.0) continue;
       if (delta > cutoff) continue;
-      rad = 0.5*diam;
+      rad = atom->shape[type[i]][0];
+      diam = 2.0*rad;
       rad2 = rad*rad;
       rad4 = rad2*rad2;
       rad8 = rad4*rad4;
@@ -162,8 +176,8 @@ void FixWallColloid::post_force(int vflag)
       r8inv = r4inv*r4inv;
       fwall = (coeff1*(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
-             - coeff2*r2inv) * side;
+                       + 21.0*rad2*rad*pow(delta,6.0))*r8inv - 
+	       coeff2*r2inv) * side;
       f[i][dim] -= fwall;
       r2 = 0.5*diam - delta;
       rinv2 = 1.0/r2;
@@ -176,9 +190,9 @@ void FixWallColloid::post_force(int vflag)
       r4inv3 = r2inv3*r2inv3;
       r6inv3 = r4inv3*r2inv3;
       wall[0] += coeff3*((-3.5*diam+delta)*r4inv2*r2inv2*rinv2
-			 + (3.5*diam+delta)*r4inv3*r2inv3*rinv3) 
-	- coeff4*((-diam*delta+r2*r3*(log(-r2)-log(r3)))*
-		  (-rinv2)*rinv3) - offset;
+			 + (3.5*diam+delta)*r4inv3*r2inv3*rinv3) -
+	coeff4*((-diam*delta+r2*r3*(log(-r2)-log(r3)))*
+		(-rinv2)*rinv3) - offset;
       wall[dim+1] += fwall;
     }
 }

src/COLLOID/fix_wall_colloid.h

@@ -34,7 +34,7 @@ class FixWallColloid : public Fix {
 
  private:
   int dim,side,thermo_flag,eflag_enable;
-  double coord,epsilon,sigma,diam,cutoff;
+  double coord,epsilon,sigma,cutoff;
   double coeff1,coeff2,coeff3,coeff4,offset;
   double wall[4],wall_all[4];
   int wall_flag;