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

src/PERI/fix_peri_neigh.cpp

@@ -120,6 +120,15 @@ void FixPeriNeigh::init_list(int id, NeighList *ptr)
   list = ptr;
 }
 
+/* ---------------------------------------------------------------------- 
+   For minimization: setup as with dynamics 
+------------------------------------------------------------------------- */
+
+void FixPeriNeigh::min_setup(int vflag)
+{
+  setup(vflag); 
+}
+
 /* ----------------------------------------------------------------------
    create initial list of neighbor partners via call to neighbor->build()
    must be done in setup (not init) since fix init comes before neigh init
@@ -231,6 +240,20 @@ void FixPeriNeigh::setup(int vflag)
     }
   }
 
+  // Sanity check: Does any atom appear twice in any neigborlist?
+  // Should only be possible if using pbc and domain not at least of width 2 \delta 
+  if (domain->xperiodic || domain->yperiodic || domain->zperiodic) {
+    for (i = 0; i < nlocal; i++) {
+      jnum = npartner[i];
+      for (jj = 0; jj < jnum; jj++) {
+        for (int kk = jj+1; kk < jnum; kk++) {
+          if (partner[i][jj] == partner[i][kk]) error->one("Duplicate particle in bond family. Check that box is greater than size 2*delta in periodic dimensions.");
+        }
+      }
+    }
+
+  }
+
   // compute wvolume for each atom
 
   double **x0 = atom->x0;

src/PERI/fix_peri_neigh.h

@@ -36,6 +36,7 @@ class FixPeriNeigh : public Fix {
   void init();
   void init_list(int, class NeighList *);
   void setup(int);
+  void min_setup(int);
 
   double memory_usage();
   void grow_arrays(int);

src/PERI/pair_peri_lps.cpp

@@ -44,6 +44,7 @@ using namespace LAMMPS_NS;
 PairPeriLPS::PairPeriLPS(LAMMPS *lmp) : Pair(lmp)			
 {
   for (int i = 0; i < 6; i++) virial[i] = 0.0;
+  no_virial_compute=1;
 
   ifix_peri = -1;
 
@@ -120,7 +121,7 @@ void PairPeriLPS::compute(int eflag, int vflag)
 
   int nall = atom->nlocal + atom->nghost;
   int newton_pair = force->newton_pair;
-  int nonperiodic = domain->nonperiodic;
+  int periodic = domain->xperiodic || domain->yperiodic || domain->zperiodic;
 
   inum = list->inum;
   ilist = list->ilist;
@@ -156,7 +157,7 @@ void PairPeriLPS::compute(int eflag, int vflag)
       delx0 = xtmp0 - x0[j][0];
       dely0 = ytmp0 - x0[j][1];
       delz0 = ztmp0 - x0[j][2];
-      if (nonperiodic == 0) domain->minimum_image(delx0,dely0,delz0);
+      if (periodic) domain->minimum_image(delx0,dely0,delz0);
       rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0;
       jtype = type[j];
  
@@ -194,8 +195,7 @@ void PairPeriLPS::compute(int eflag, int vflag)
         }
 
         if (eflag) evdwl = 0.5*rk*dr;
-	if (evflag) ev_tally(i,j,nlocal,newton_pair,
+	if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair*vfrac[i],delx,dely,delz);
-			     evdwl,0.0,fpair,delx,dely,delz);
       }
     }
   }
@@ -265,12 +265,12 @@ void PairPeriLPS::compute(int eflag, int vflag)
       delx = xtmp - x[j][0];
       dely = ytmp - x[j][1];
       delz = ztmp - x[j][2];
-      if (nonperiodic == 0) domain->minimum_image(delx,dely,delz);
+      if (periodic) domain->minimum_image(delx,dely,delz);
       rsq = delx*delx + dely*dely + delz*delz;
       delx0 = xtmp0 - x0[j][0];						
       dely0 = ytmp0 - x0[j][1];						
       delz0 = ztmp0 - x0[j][2];						
-      if (nonperiodic == 0) domain->minimum_image(delx0,dely0,delz0);   
+      if (periodic) domain->minimum_image(delx0,dely0,delz0);   
       jtype = type[j];
       delta = sqrt(cutsq[itype][jtype]);
       r = sqrt(rsq);
@@ -309,8 +309,7 @@ void PairPeriLPS::compute(int eflag, int vflag)
       if (eflag) evdwl = 0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) * 
 		   omega_plus*(deviatoric_extension * deviatoric_extension) *
 		   vfrac[j] * vfrac_scale;
-      if (evflag) ev_tally(i,i,nlocal,0,
+      if (evflag) ev_tally(i,i,nlocal,0,0.5*evdwl,0.0,0.5*fbond*vfrac[i],delx,dely,delz);
-			   0.5*evdwl,0.0,0.5*fbond,delx,dely,delz);
 
       // find stretch in bond I-J and break if necessary
       // use s0 from previous timestep
@@ -330,8 +329,6 @@ void PairPeriLPS::compute(int eflag, int vflag)
     }
   }
 
-  if (vflag_fdotr) virial_compute();
-
   // store new s0
   for (i = 0; i < nlocal; i++) s0[i] = s0_new[i];
 
@@ -542,7 +539,8 @@ double PairPeriLPS::single(int i, int j, int itype, int jtype,
   delx0 = x0[i][0] - x0[j][0];
   dely0 = x0[i][1] - x0[j][1];
   delz0 = x0[i][2] - x0[j][2];
-  if (domain->nonperiodic == 0) domain->minimum_image(delx0,dely0,delz0);
+  int periodic = domain->xperiodic || domain->yperiodic || domain->zperiodic;
+  if (periodic) domain->minimum_image(delx0,dely0,delz0);
   rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0;
 
   d_ij = MIN(0.9*sqrt(rsq0),1.35*lc);
@@ -663,7 +661,7 @@ void PairPeriLPS::compute_dilatation()
   int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner;
   double *wvolume = ((FixPeriNeigh *) modify->fix[ifix_peri])->wvolume;
 
-  int nonperiodic = domain->nonperiodic;
+  int periodic = domain->xperiodic || domain->yperiodic || domain->zperiodic;
 
   // compute the dilatation theta
 
@@ -693,12 +691,12 @@ void PairPeriLPS::compute_dilatation()
       delx = xtmp - x[j][0];
       dely = ytmp - x[j][1];
       delz = ztmp - x[j][2];
-      if (nonperiodic == 0) domain->minimum_image(delx,dely,delz);
+      if (periodic) domain->minimum_image(delx,dely,delz);
       rsq = delx*delx + dely*dely + delz*delz;
       delx0 = xtmp0 - x0[j][0];
       dely0 = ytmp0 - x0[j][1];
       delz0 = ztmp0 - x0[j][2];
-      if (nonperiodic == 0) domain->minimum_image(delx0,dely0,delz0);
+      if (periodic) domain->minimum_image(delx0,dely0,delz0);
 
       r = sqrt(rsq);
       dr = r - r0[i][jj];

src/PERI/pair_peri_pmb.cpp

@@ -45,6 +45,7 @@ using namespace LAMMPS_NS;
 PairPeriPMB::PairPeriPMB(LAMMPS *lmp) : Pair(lmp)
 {
   for (int i = 0; i < 6; i++) virial[i] = 0.0;
+  no_virial_compute=1;
 
   ifix_peri = -1;
 
@@ -112,7 +113,7 @@ void PairPeriPMB::compute(int eflag, int vflag)
 
   int nall = atom->nlocal + atom->nghost;
   int newton_pair = force->newton_pair;
-  int nonperiodic = domain->nonperiodic;
+  int periodic = (domain->xperiodic || domain->yperiodic || domain->zperiodic);
 
   inum = list->inum;
   ilist = list->ilist;
@@ -145,7 +146,7 @@ void PairPeriPMB::compute(int eflag, int vflag)
       delx0 = xtmp0 - x0[j][0];
       dely0 = ytmp0 - x0[j][1];
       delz0 = ztmp0 - x0[j][2];
-      if (nonperiodic == 0) domain->minimum_image(delx0,dely0,delz0);
+      if (periodic) domain->minimum_image(delx0,dely0,delz0);
       rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0;
       jtype = type[j];
  
@@ -178,8 +179,7 @@ void PairPeriPMB::compute(int eflag, int vflag)
         }
 
         if (eflag) evdwl = 0.5*rk*dr;
-	if (evflag) ev_tally(i,j,nlocal,newton_pair,
+	if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair*vfrac[i],delx,dely,delz);
-			     evdwl,0.0,fpair,delx,dely,delz);
       }
     }
   }
@@ -224,7 +224,7 @@ void PairPeriPMB::compute(int eflag, int vflag)
       delx = xtmp - x[j][0];
       dely = ytmp - x[j][1];
       delz = ztmp - x[j][2];
-      if (nonperiodic == 0) domain->minimum_image(delx,dely,delz);
+      if (periodic) domain->minimum_image(delx,dely,delz);
       rsq = delx*delx + dely*dely + delz*delz;
       jtype = type[j];
       delta = sqrt(cutsq[itype][jtype]);
@@ -254,8 +254,7 @@ void PairPeriPMB::compute(int eflag, int vflag)
       // since I-J is double counted, set newton off & use 1/2 factor and I,I 
 
       if (eflag) evdwl = 0.5*rk*dr;
-      if (evflag) ev_tally(i,i,nlocal,0,
+      if (evflag) ev_tally(i,i,nlocal,0,0.5*evdwl,0.0,0.5*fbond*vfrac[i],delx,dely,delz);
-			   0.5*evdwl,0.0,0.5*fbond,delx,dely,delz);
 
       // find stretch in bond I-J and break if necessary
       // use s0 from previous timestep
@@ -267,16 +266,12 @@ void PairPeriPMB::compute(int eflag, int vflag)
       if (first)
          s0_new[i] = s00[itype][jtype] - (alpha[itype][jtype] * stretch);
       else
-         s0_new[i] = MAX(s0_new[i],
+         s0_new[i] = MAX(s0_new[i],s00[itype][jtype] - (alpha[itype][jtype] * stretch));
-			 s00[itype][jtype] - (alpha[itype][jtype] * stretch));
       first = false;
     }
   }
 
-  if (vflag_fdotr) virial_compute();
-
   // store new s0
-
   for (i = 0; i < nlocal; i++) s0[i] = s0_new[i]; 
 }
 
@@ -474,7 +469,8 @@ double PairPeriPMB::single(int i, int j, int itype, int jtype, double rsq,
   delx0 = x0[i][0] - x0[j][0];
   dely0 = x0[i][1] - x0[j][1];
   delz0 = x0[i][2] - x0[j][2];
-  if (domain->nonperiodic == 0) domain->minimum_image(delx0,dely0,delz0);
+  int periodic = domain->xperiodic || domain->yperiodic || domain->zperiodic;
+  if (periodic) domain->minimum_image(delx0,dely0,delz0);
   rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0;
 
   d_ij = MIN(0.9*sqrt(rsq0),1.35*lc);