Commit c57c14d9 authored by sjplimp's avatar sjplimp
Browse files

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5092 f3b2605a-c512-4ea7-a41b-209d697bcdaa
parent e0302dd5
Loading
Loading
Loading
Loading
+27 −11
Original line number Diff line number Diff line
@@ -72,15 +72,10 @@ void FixWallColloid::init()

void FixWallColloid::precompute(int m)
{
  coeff1[m] = -4.0/315.0 * epsilon[m] * pow(sigma[m],6.0);
  coeff2[m] = -2.0/3.0 * epsilon[m];
  coeff1[m] = 4.0/315.0 * epsilon[m] * pow(sigma[m],6.0);
  coeff2[m] = 2.0/3.0 * epsilon[m];
  coeff3[m] = epsilon[m] * pow(sigma[m],6.0)/7560.0;
  coeff4[m] = epsilon[m]/6.0;

  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;
}

/* ----------------------------------------------------------------------
@@ -96,6 +91,7 @@ void FixWallColloid::wall_particle(int m, int which, double coord)
  double r2,rinv2,r2inv2,r4inv2,r6inv2;
  double r3,rinv3,r2inv3,r4inv3,r6inv3;
  double rad,rad2,rad4,rad8,diam,new_coeff2;
  double eoffset;

  double **x = atom->x;
  double **f = atom->f;
@@ -136,20 +132,40 @@ void FixWallColloid::wall_particle(int m, int which, double coord)
				 + 21.0*rad2*rad*pow(delta,6.0))*r8inv - 
		      new_coeff2*r2inv);
      f[i][dim] -= fwall;
      r2 = 0.5*diam - delta;

      r2 = rad - delta;
      rinv2 = 1.0/r2;
      r2inv2 = rinv2*rinv2;
      r4inv2 = r2inv2*r2inv2;
      r6inv2 = r4inv2*r2inv2;
      r3 = delta + 0.5*diam;
      r3 = delta + rad;
      rinv3 = 1.0/r3;
      r2inv3 = rinv3*rinv3;
      r4inv3 = r2inv3*r2inv3;
      r6inv3 = r4inv3*r2inv3;
      ewall[0] += coeff3[m]*((-3.5*diam+delta)*r4inv2*r2inv2*rinv2
			     + (3.5*diam+delta)*r4inv3*r2inv3*rinv3) -
	coeff4[m]*((-diam*delta+r2*r3*(log(-r2)-log(r3)))*
		   (-rinv2)*rinv3) - offset[m];
	coeff4[m]*((diam*delta-r2*r3*(log(-r2)-log(r3)))*
		   (-rinv2)*rinv3);

      // offset depends on particle size

      r2 = rad - cutoff[m];
      rinv2 = 1.0/r2;
      r2inv2 = rinv2*rinv2;
      r4inv2 = r2inv2*r2inv2;
      r6inv2 = r4inv2*r2inv2;
      r3 = cutoff[m] + rad;
      rinv3 = 1.0/r3;
      r2inv3 = rinv3*rinv3;
      r4inv3 = r2inv3*r2inv3;
      r6inv3 = r4inv3*r2inv3;
      eoffset = coeff3[m]*((-3.5*diam+cutoff[m])*r4inv2*r2inv2*rinv2
			   + (3.5*diam+cutoff[m])*r4inv3*r2inv3*rinv3) -
	coeff4[m]*((diam*cutoff[m]-r2*r3*(log(-r2)-log(r3)))*
		   (-rinv2)*rinv3);
      ewall[0] -= eoffset;

      ewall[m+1] += fwall;
    }