Modifed ZBL to accomodate mixing

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

src/pair_zbl.cpp

@@ -199,11 +199,14 @@ void PairZBL::settings(int narg, char **arg)
 
 /* ----------------------------------------------------------------------
    set coeffs for one or more type pairs
+
 ------------------------------------------------------------------------- */
 
 void PairZBL::coeff(int narg, char **arg)
 {
-  if (narg != 3)
+  double z_one, z_two;
+
+  if (narg != 4)
     error->all(FLERR,"Incorrect args for pair coefficients");
   if (!allocated) allocate();
 
@@ -213,7 +216,8 @@ void PairZBL::coeff(int narg, char **arg)
   int jlo,jhi;
   force->bounds(arg[1],atom->ntypes,jlo,jhi);
 
-  double z_one = force->numeric(FLERR,arg[2]);
+  z_one = force->numeric(FLERR,arg[2]);
+  z_two = force->numeric(FLERR,arg[3]);
 
   // set flag for each i-j pair
   // set z-parameter only for i-i pairs
@@ -221,8 +225,13 @@ void PairZBL::coeff(int narg, char **arg)
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
     for (int j = MAX(jlo,i); j <= jhi; j++) {
-      if (i == j) z[i] = z_one;
+      if (i == j) {
+	if (z_one != z_two)
+	  error->all(FLERR,"Incorrect args for pair coefficients");
+	z[i] = z_one;
+      }
       setflag[i][j] = 1;
+      set_coeff(i, j, z_one, z_two);
       count++;
     }
   }
@@ -248,63 +257,8 @@ void PairZBL::init_style()
 
 double PairZBL::init_one(int i, int j)
 {
-
+  if (setflag[i][j] == 0)
-  double ainv = (pow(z[i],pzbl) + pow(z[j],pzbl))/(a0*force->angstrom);
+    set_coeff(i, j, z[i], z[j]);
-  d1a[i][j] = d1*ainv;
-  d2a[i][j] = d2*ainv;
-  d3a[i][j] = d3*ainv;
-  d4a[i][j] = d4*ainv;
-  zze[i][j] = z[i]*z[j]*force->qqr2e*force->qelectron*force->qelectron;
-
-  d1a[j][i] = d1a[i][j];
-  d2a[j][i] = d2a[i][j];
-  d3a[j][i] = d3a[i][j];
-  d4a[j][i] = d4a[i][j];
-  zze[j][i] = zze[i][j];
-
-  // e =  t^3 (sw3 + sw4*t) + sw5 
-  //   = A/3*t^3 + B/4*t^4 + C
-  // sw3 = A/3 
-  // sw4 = B/4 
-  // sw5 = C
-
-  // dedr = t^2 (sw1 + sw2*t) 
-  //      = A*t^2 + B*t^3
-  // sw1 = A 
-  // sw2 = B 
-
-  // de2dr2 = 2*A*t + 3*B*t^2
-
-  // Require that at t = tc:
-  // e = -Fc
-  // dedr = -Fc'
-  // d2edr2 = -Fc'' 
-
-  // Hence:
-  // A = (-3Fc' + tc*Fc'')/tc^2
-  // B = ( 2Fc' - tc*Fc'')/tc^3
-  // C = -Fc + tc/2*Fc' - tc^2/12*Fc''
-
-  double tc = cut_global - cut_inner;
-  double fc = e_zbl(cut_global, i, j);
-  double fcp = dzbldr(cut_global, i, j);
-  double fcpp = d2zbldr2(cut_global, i, j);
-
-  double swa = (-3.0*fcp + tc*fcpp)/(tc*tc);
-  double swb = ( 2.0*fcp - tc*fcpp)/(tc*tc*tc);
-  double swc = -fc + (tc/2.0)*fcp - (tc*tc/12.0)*fcpp; 
-  
-  sw1[i][j] = swa;
-  sw2[i][j] = swb;
-  sw3[i][j] = swa/3.0;
-  sw4[i][j] = swb/4.0;
-  sw5[i][j] = swc;
-
-  sw1[j][i] = sw1[i][j];
-  sw2[j][i] = sw2[i][j];
-  sw3[j][i] = sw3[i][j];
-  sw4[j][i] = sw4[i][j];
-  sw5[j][i] = sw5[i][j];
 
   return cut_global;
 }
@@ -434,3 +388,67 @@ double PairZBL::d2zbldr2(double r, int i, int j) {
   return result;
 }
 
+/* ----------------------------------------------------------------------
+   calculate the i,j entries in the various coeff arrays 
+------------------------------------------------------------------------- */
+
+void PairZBL::set_coeff(int i, int j, double zi, double zj)
+{
+  double ainv = (pow(zi,pzbl) + pow(zj,pzbl))/(a0*force->angstrom);
+  d1a[i][j] = d1*ainv;
+  d2a[i][j] = d2*ainv;
+  d3a[i][j] = d3*ainv;
+  d4a[i][j] = d4*ainv;
+  zze[i][j] = zi*zj*force->qqr2e*force->qelectron*force->qelectron;
+  
+  d1a[j][i] = d1a[i][j];
+  d2a[j][i] = d2a[i][j];
+  d3a[j][i] = d3a[i][j];
+  d4a[j][i] = d4a[i][j];
+  zze[j][i] = zze[i][j];
+  
+  // e =  t^3 (sw3 + sw4*t) + sw5 
+  //   = A/3*t^3 + B/4*t^4 + C
+  // sw3 = A/3 
+  // sw4 = B/4 
+  // sw5 = C
+  
+  // dedr = t^2 (sw1 + sw2*t) 
+  //      = A*t^2 + B*t^3
+  // sw1 = A 
+  // sw2 = B 
+  
+  // de2dr2 = 2*A*t + 3*B*t^2
+  
+  // Require that at t = tc:
+  // e = -Fc
+  // dedr = -Fc'
+  // d2edr2 = -Fc'' 
+  
+  // Hence:
+  // A = (-3Fc' + tc*Fc'')/tc^2
+  // B = ( 2Fc' - tc*Fc'')/tc^3
+  // C = -Fc + tc/2*Fc' - tc^2/12*Fc''
+  
+  double tc = cut_global - cut_inner;
+  double fc = e_zbl(cut_global, i, j);
+  double fcp = dzbldr(cut_global, i, j);
+  double fcpp = d2zbldr2(cut_global, i, j);
+  
+  double swa = (-3.0*fcp + tc*fcpp)/(tc*tc);
+  double swb = ( 2.0*fcp - tc*fcpp)/(tc*tc*tc);
+  double swc = -fc + (tc/2.0)*fcp - (tc*tc/12.0)*fcpp; 
+  
+  sw1[i][j] = swa;
+  sw2[i][j] = swb;
+  sw3[i][j] = swa/3.0;
+  sw4[i][j] = swb/4.0;
+  sw5[i][j] = swc;
+  
+  sw1[j][i] = sw1[i][j];
+  sw2[j][i] = sw2[i][j];
+  sw3[j][i] = sw3[i][j];
+  sw4[j][i] = sw4[i][j];
+  sw5[j][i] = sw5[i][j];
+}
+

src/pair_zbl.h

@@ -46,6 +46,7 @@ class PairZBL : public Pair {
   double e_zbl(double, int, int);
   double dzbldr(double, int, int);
   double d2zbldr2(double, int, int);
+  void set_coeff(int, int, double, double);
 };
 
 namespace PairZBLConstants {