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

src/CLASS2/bond_class2.cpp

@@ -197,12 +197,15 @@ void BondClass2::read_restart(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-double BondClass2::single(int type, double rsq, int i, int j)
+double BondClass2::single(int type, double rsq, int i, int j, double &fforce)
 {
   double r = sqrt(rsq);
   double dr = r - r0[type];
   double dr2 = dr*dr;
   double dr3 = dr2*dr;
   double dr4 = dr3*dr;
+  double de_bond = 2.0*k2[type]*dr + 3.0*k3[type]*dr2 + 4.0*k4[type]*dr3;
+  if (r > 0.0) fforce = -de_bond/r;
+  else fforce = 0.0;
   return (k2[type]*dr2 + k3[type]*dr3 + k4[type]*dr4);
 }

src/CLASS2/bond_class2.h

@@ -34,7 +34,7 @@ class BondClass2 : public Bond {
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
-  double single(int, double, int, int);
+  double single(int, double, int, int, double &);
 
  protected:
   double *r0,*k2,*k3,*k4;

src/MC/fix_bond_swap.cpp

@@ -650,8 +650,9 @@ double FixBondSwap::pair_eng(int i, int j)
 
 double FixBondSwap::bond_eng(int btype, int i, int j)
 {
+  double tmp;
   double rsq = dist_rsq(i,j);
-  return force->bond->single(btype,rsq,i,j);
+  return force->bond->single(btype,rsq,i,j,tmp);
 }
 
 /* ---------------------------------------------------------------------- */

src/compute_bond_local.cpp

@@ -27,6 +27,8 @@ using namespace LAMMPS_NS;
 
 #define DELTA 10000
 
+enum{DIST,ENG,FORCE};
+
 /* ---------------------------------------------------------------------- */
 
 ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
@@ -42,17 +44,22 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
   if (nvalues == 1) size_local_cols = 0;
   else size_local_cols = nvalues;
 
-  dflag = eflag = -1;
+  bstyle = new int[nvalues];
-  nvalues = 0;
 
-  int i;
+  nvalues = 0;
   for (int iarg = 3; iarg < narg; iarg++) {
-    i = iarg-3;
+    if (strcmp(arg[iarg],"dist") == 0) bstyle[nvalues++] = DIST;
-    if (strcmp(arg[iarg],"dist") == 0) dflag = nvalues++;
+    else if (strcmp(arg[iarg],"eng") == 0) bstyle[nvalues++] = ENG;
-    else if (strcmp(arg[iarg],"eng") == 0) eflag = nvalues++;
+    else if (strcmp(arg[iarg],"force") == 0) bstyle[nvalues++] = FORCE;
     else error->all(FLERR,"Invalid keyword in compute bond/local command");
   }
 
+  // set singleflag if need to call bond->single()
+
+  singleflag = 0;
+  for (int i = 0; i < nvalues; i++)
+    if (bstyle[i] != DIST) singleflag = 1;
+
   nmax = 0;
   vector = NULL;
   array = NULL;
@@ -64,6 +71,7 @@ ComputeBondLocal::~ComputeBondLocal()
 {
   memory->destroy(vector);
   memory->destroy(array);
+  delete [] bstyle;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -109,7 +117,8 @@ int ComputeBondLocal::compute_bonds(int flag)
 {
   int i,m,n,atom1,atom2;
   double delx,dely,delz,rsq;
-  double *dbuf,*ebuf;
+  double *dbuf,*ebuf,*fbuf;
+  double *ptr;
 
   double **x = atom->x;
   int *num_bond = atom->num_bond;
@@ -120,22 +129,10 @@ int ComputeBondLocal::compute_bonds(int flag)
   int nlocal = atom->nlocal;
   int newton_bond = force->newton_bond;
 
-  if (flag) {
-    if (nvalues == 1) {
-      if (dflag >= 0) dbuf = vector;
-      if (eflag >= 0) ebuf = vector;
-    } else {
-      if (dflag >= 0 && array) dbuf = &array[0][dflag];
-      else dbuf = NULL;
-      if (eflag >= 0 && array) ebuf = &array[0][eflag];
-      else ebuf = NULL;
-    }
-  }
-
   Bond *bond = force->bond;
+  double eng,fbond;
 
   m = n = 0;
-  double fforce; // unused
   for (atom1 = 0; atom1 < nlocal; atom1++) {
     if (!(mask[atom1] & groupbit)) continue;
     for (i = 0; i < num_bond[atom1]; i++) {
@@ -150,13 +147,29 @@ int ComputeBondLocal::compute_bonds(int flag)
         delz = x[atom1][2] - x[atom2][2];
         domain->minimum_image(delx,dely,delz);
         rsq = delx*delx + dely*dely + delz*delz;
-        if (dflag >= 0) dbuf[n] = sqrt(rsq);
+
-        if (eflag >= 0) {
+        if (singleflag) {
           if (bond_type[atom1][i] > 0)
-            ebuf[n] = bond->single(bond_type[atom1][i],rsq,atom1,atom2,fforce);
+            eng = bond->single(bond_type[atom1][i],rsq,atom1,atom2,fbond);
-          else ebuf[n] = 0.0;
+          else eng = fbond = 0.0;
+        }
+
+        if (nvalues == 1) ptr = &vector[m];
+        else ptr = array[m];
+
+        for (n = 0; n < nvalues; n++) {
+          switch (bstyle[n]) {
+          case DIST:
+            ptr[n] = sqrt(rsq);
+            break;
+          case ENG:
+            ptr[n] = eng;
+            break;
+          case FORCE:
+            ptr[n] = sqrt(rsq)*fbond;
+            break;
+          }
         }
-        n += nvalues;
       }
 
       m++;

src/compute_bond_local.h

@@ -35,6 +35,8 @@ class ComputeBondLocal : public Compute {
  private:
   int nvalues,dflag,eflag;
   int ncount;
+  int *bstyle;
+  int singleflag;
 
   int nmax;
   double *vector;

src/fix.h

@@ -122,6 +122,7 @@ class Fix : protected Pointers {
   virtual void final_integrate_respa(int, int) {}
 
   virtual void min_setup_pre_exchange() {}
+  virtual void min_setup_pre_neighbor() {}
   virtual void min_setup_pre_force(int) {}
   virtual void min_pre_exchange() {}
   virtual void min_pre_force(int) {}
@@ -189,11 +190,12 @@ namespace FixConst {
   static const int POST_FORCE_RESPA =        1<<12;
   static const int FINAL_INTEGRATE_RESPA =   1<<13;
   static const int MIN_PRE_EXCHANGE =        1<<14;
-  static const int MIN_PRE_FORCE =           1<<15;
+  static const int MIN_PRE_NEIGHBOR =        1<<15;
-  static const int MIN_POST_FORCE =          1<<16;
+  static const int MIN_PRE_FORCE =           1<<16;
-  static const int MIN_ENERGY =              1<<17;
+  static const int MIN_POST_FORCE =          1<<17;
-  static const int POST_RUN =                1<<18;
+  static const int MIN_ENERGY =              1<<18;
-  static const int FIX_CONST_LAST =          1<<19;
+  static const int POST_RUN =                1<<19;
+  static const int FIX_CONST_LAST =          1<<20;
 }
 
 }

src/modify.cpp

@@ -56,8 +56,9 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp)
   list_initial_integrate_respa = list_post_integrate_respa = NULL;
   list_pre_force_respa = list_post_force_respa = NULL;
   list_final_integrate_respa = NULL;
-  list_min_pre_exchange = list_min_pre_force =
+  list_min_pre_exchange = list_pre_neighbor = NULL;
-  list_min_post_force = list_min_energy = NULL;
+  list_min_pre_force = list_min_post_force = NULL;
+  list_min_energy = NULL;
 
   end_of_step_every = NULL;
 
@@ -104,6 +105,7 @@ Modify::~Modify()
   delete [] list_post_force_respa;
   delete [] list_final_integrate_respa;
   delete [] list_min_pre_exchange;
+  delete [] list_min_pre_neighbor;
   delete [] list_min_pre_force;
   delete [] list_min_post_force;
   delete [] list_min_energy;
@@ -150,6 +152,7 @@ void Modify::init()
             n_final_integrate_respa,list_final_integrate_respa);
 
   list_init(MIN_PRE_EXCHANGE,n_min_pre_exchange,list_min_pre_exchange);
+  list_init(MIN_PRE_NEIGHBOR,n_min_pre_neighbor,list_min_pre_neighbor);
   list_init(MIN_PRE_FORCE,n_min_pre_force,list_min_pre_force);
   list_init(MIN_POST_FORCE,n_min_post_force,list_min_post_force);
   list_init(MIN_ENERGY,n_min_energy,list_min_energy);
@@ -255,6 +258,9 @@ void Modify::setup_pre_neighbor()
   if (update->whichflag == 1)
     for (int i = 0; i < n_pre_neighbor; i++)
       fix[list_pre_neighbor[i]]->setup_pre_neighbor();
+  else if (update->whichflag == 2)
+    for (int i = 0; i < n_pre_neighbor; i++)
+      fix[list_min_pre_neighbor[i]]->min_setup_pre_neighbor();
 }
 
 /* ----------------------------------------------------------------------

src/modify.h

@@ -27,7 +27,8 @@ class Modify : protected Pointers {
   int n_final_integrate,n_end_of_step,n_thermo_energy;
   int n_initial_integrate_respa,n_post_integrate_respa;
   int n_pre_force_respa,n_post_force_respa,n_final_integrate_respa;
-  int n_min_pre_exchange,n_min_pre_force,n_min_post_force,n_min_energy;
+  int n_min_pre_exchange,n_min_pre_neighbor;
+  int n_min_pre_force,n_min_post_force,n_min_energy;
 
   int restart_pbc_any;       // 1 if any fix sets restart_pbc
   int nfix_restart_global;   // stored fix global info from restart file
@@ -110,8 +111,9 @@ class Modify : protected Pointers {
   int *list_initial_integrate_respa,*list_post_integrate_respa;
   int *list_pre_force_respa,*list_post_force_respa;
   int *list_final_integrate_respa;
-  int *list_min_pre_exchange,*list_min_pre_force;
+  int *list_min_pre_exchange,*list_min_pre_neighbor;
-  int *list_min_post_force,*list_min_energy;
+  int *list_min_pre_force,*list_min_post_force;
+  int *list_min_energy;
 
   int *end_of_step_every;