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

src/PERI/Install.sh

@@ -4,11 +4,13 @@ if (test $1 = 1) then
 
   cp atom_vec_peri.cpp ..
   cp pair_peri_pmb.cpp ..
+  cp pair_peri_lps.cpp ..
   cp fix_peri_neigh.cpp ..
   cp compute_damage_atom.cpp ..
 
   cp atom_vec_peri.h ..
   cp pair_peri_pmb.h ..
+  cp pair_peri_lps.h ..
   cp fix_peri_neigh.h ..
   cp compute_damage_atom.h ..
 
@@ -16,11 +18,13 @@ elif (test $1 = 0) then
 
   rm -f ../atom_vec_peri.cpp
   rm -f ../pair_peri_pmb.cpp
+  rm -f ../pair_peri_lps.cpp
   rm -f ../fix_peri_neigh.cpp
   rm -f ../compute_damage_atom.cpp
 
   rm -f ../atom_vec_peri.h
   rm -f ../pair_peri_pmb.h
+  rm -f ../pair_peri_lps.h
   rm -f ../fix_peri_neigh.h
   rm -f ../compute_damage_atom.h
 

src/PERI/fix_peri_neigh.cpp

@@ -17,18 +17,21 @@
 
 #include "math.h"
 #include "fix_peri_neigh.h"
+#include "pair_peri_pmb.h"
+#include "pair_peri_lps.h"
 #include "atom.h"
+#include "domain.h"
 #include "force.h"
-#include "pair.h"
+#include "comm.h"
+#include "update.h"
 #include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
+#include "pair.h"
+#include "lattice.h"
 #include "memory.h"
 #include "error.h"
 
-#include "comm.h"
-#include "update.h"
-
 using namespace LAMMPS_NS;
 
 #define MIN(A,B) ((A) < (B)) ? (A) : (B)
@@ -39,6 +42,7 @@ using namespace LAMMPS_NS;
 FixPeriNeigh::FixPeriNeigh(LAMMPS *lmp,int narg, char **arg) : 
   Fix(lmp, narg, arg)
 {
+  restart_global = 1;
   restart_peratom = 1;
   first = 1;
 
@@ -51,6 +55,7 @@ FixPeriNeigh::FixPeriNeigh(LAMMPS *lmp,int narg, char **arg) :
   partner = NULL;
   r0 = NULL;
   vinter = NULL;
+  wvolume = NULL;
 
   grow_arrays(atom->nmax);
   atom->add_callback(0);
@@ -60,6 +65,10 @@ FixPeriNeigh::FixPeriNeigh(LAMMPS *lmp,int narg, char **arg) :
 
   int nlocal = atom->nlocal;
   for (int i = 0; i < nlocal; i++) npartner[i] = 0;
+
+  // set comm sizes needed by this fix
+
+  comm_forward = 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -78,6 +87,7 @@ FixPeriNeigh::~FixPeriNeigh()
   memory->destroy_2d_int_array(partner);
   memory->destroy_2d_double_array(r0);
   memory->sfree(vinter);
+  memory->sfree(wvolume);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -92,6 +102,8 @@ int FixPeriNeigh::setmask()
 
 void FixPeriNeigh::init()
 {
+  if (!first) return;
+
   // need a full neighbor list once
 
   int irequest = neighbor->request((void *) this);
@@ -127,7 +139,7 @@ void FixPeriNeigh::setup(int vflag)
   int *tag = atom->tag;
   int nlocal = atom->nlocal;
 
-  // only build list of bonds on very first run
+ // only build list of bonds on very first run
 
   if (!first) return;
   first = 0;
@@ -190,6 +202,7 @@ void FixPeriNeigh::setup(int vflag)
   for (i = 0; i < nlocal; i++) {
     npartner[i] = 0;
     vinter[i] = 0.0;
+    wvolume[i] = 0.0;
   }
 
   for (ii = 0; ii < inum; ii++) {
@@ -203,7 +216,7 @@ void FixPeriNeigh::setup(int vflag)
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
- 
+
       delx = xtmp - x[j][0];
       dely = ytmp - x[j][1];
       delz = ztmp - x[j][2];
@@ -211,14 +224,74 @@ void FixPeriNeigh::setup(int vflag)
       jtype = type[j];
 
       if (rsq <= cutsq[itype][jtype]) {
-	partner[i][npartner[i]] = tag[j];
+        partner[i][npartner[i]] = tag[j];
-	r0[i][npartner[i]] = sqrt(rsq);
+        r0[i][npartner[i]] = sqrt(rsq);
-	npartner[i]++;
+        npartner[i]++;
         vinter[i] += vfrac[j];
       }
     }
   }
 
+  // compute wvolume for each atom
+
+  double **x0 = atom->x0;
+  double half_lc = 0.5*(domain->lattice->xlattice);
+  double vfrac_scale;
+  PairPeriLPS *pairlps = dynamic_cast<PairPeriLPS*>(anypair);
+  PairPeriPMB *pairpmb = dynamic_cast<PairPeriPMB*>(anypair);
+  for (i = 0; i < nlocal; i++) {
+    double xtmp0 = x0[i][0];
+    double ytmp0 = x0[i][1];
+    double ztmp0 = x0[i][2];
+    jnum = npartner[i];
+    itype = type[i];
+
+    // loop over partners of particle i
+
+    for (jj = 0; jj < jnum; jj++) {
+
+      // if bond already broken, skip this partner
+
+      if (partner[i][jj] == 0) continue;
+
+      // lookup local index of partner particle
+
+      j = atom->map(partner[i][jj]);
+
+      // skip if particle is "lost"
+
+      if (j < 0) continue;
+
+      double delx0 = xtmp0 - x0[j][0];
+      double dely0 = ytmp0 - x0[j][1];
+      double delz0 = ztmp0 - x0[j][2];
+      double rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0;
+
+      jtype = type[j];
+      double delta = sqrt(cutsq[itype][jtype]);
+
+      // scale vfrac[j] if particle j near the horizon
+
+      if ((fabs(r0[i][jj] - delta)) <= half_lc)
+        vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) +
+          (1.0 + ((delta - half_lc)/(2*half_lc) ) );
+      else vfrac_scale = 1.0;
+
+      if (pairpmb != NULL) // define influence function to be 1.0
+        wvolume[i] += 1.0 * rsq0 * vfrac[j] * vfrac_scale;
+      else if (pairlps != NULL) // call the PairPeriLPS influence function
+        wvolume[i] += pairlps->influence_function(delx0,dely0,delz0) * 
+	  rsq0 * vfrac[j] * vfrac_scale;
+      else
+        error->all("Unknown peridynamic pair style in FixPeriNeigh.");
+
+    }
+  }
+
+  // communicate wvolume to ghosts
+
+  comm->forward_comm_fix(this);
+
   // bond statistics
 
   int n = 0;
@@ -251,6 +324,7 @@ double FixPeriNeigh::memory_usage()
   bytes += nmax*maxpartner * sizeof(int);
   bytes += nmax*maxpartner * sizeof(double);
   bytes += nmax * sizeof(double);
+  bytes += nmax * sizeof(double);
   return bytes;
 }
 
@@ -267,6 +341,8 @@ void FixPeriNeigh::grow_arrays(int nmax)
   r0 = memory->grow_2d_double_array(r0,nmax,maxpartner,"peri_neigh:r0");
   vinter = (double *) memory->srealloc(vinter,nmax*sizeof(double),
 				       "peri_neigh:vinter");
+  wvolume = (double *) memory->srealloc(wvolume,nmax*sizeof(double),
+                                       "peri_neigh:wvolume");
 }
 
 /* ----------------------------------------------------------------------
@@ -281,6 +357,7 @@ void FixPeriNeigh::copy_arrays(int i, int j)
     r0[j][m] = r0[i][m];
   }
   vinter[j] = vinter[i];
+  wvolume[j] = wvolume[i];
 }
 
 /* ----------------------------------------------------------------------
@@ -301,6 +378,7 @@ int FixPeriNeigh::pack_exchange(int i, double *buf)
 
   buf[0] = m/2;
   buf[m++] = vinter[i];
+  buf[m++] = wvolume[i];
   return m;
 }
 
@@ -317,9 +395,73 @@ int FixPeriNeigh::unpack_exchange(int nlocal, double *buf)
     r0[nlocal][n] = buf[m++];
   }
   vinter[nlocal] = buf[m++];
+  wvolume[nlocal] = buf[m++];
   return m;
 }
 
+/* ---------------------------------------------------------------------- */
+
+int FixPeriNeigh::pack_comm(int n, int *list, double *buf,
+                             int pbc_flag, int *pbc)
+{
+  int i,j,m;
+
+  m = 0;
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    buf[m++] = wvolume[j];
+  }
+  return 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPeriNeigh::unpack_comm(int n, int first, double *buf)
+{
+  int i,m,last;
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) {
+    wvolume[i] = buf[m++];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   pack entire state of Fix into one write
+------------------------------------------------------------------------- */
+
+void FixPeriNeigh::write_restart(FILE *fp)
+{
+  int n = 0;
+  double list[2];
+  list[n++] = first;
+  list[n++] = maxpartner;
+
+  if (comm->me == 0) {
+    int size = n * sizeof(double);
+    fwrite(&size,sizeof(int),1,fp);
+    fwrite(list,sizeof(double),n,fp);
+  }
+}
+
+/* ----------------------------------------------------------------------
+   use state info from restart file to restart the Fix
+------------------------------------------------------------------------- */
+
+void FixPeriNeigh::restart(char *buf)
+{
+  int n = 0;
+  double *list = (double *) buf;
+
+  first = static_cast<int> (list[n++]);
+  maxpartner = static_cast<int> (list[n++]);
+
+  // grow 2D arrays now, cannot change size of 2nd array index later
+
+  grow_arrays(atom->nmax);
+}
+
 /* ----------------------------------------------------------------------
    pack values in local atom-based arrays for restart file
 ------------------------------------------------------------------------- */
@@ -327,13 +469,14 @@ int FixPeriNeigh::unpack_exchange(int nlocal, double *buf)
 int FixPeriNeigh::pack_restart(int i, double *buf)
 {
   int m = 0;
-  buf[m++] = 2*npartner[i] + 2;
+  buf[m++] = 2*npartner[i] + 4;
   buf[m++] = npartner[i];
   for (int n = 0; n < npartner[i]; n++) {
     buf[m++] = partner[i][n];
     buf[m++] = r0[i][n];
   }
   buf[m++] = vinter[i];
+  buf[m++] = wvolume[i];
   return m;
 }
 
@@ -357,6 +500,7 @@ void FixPeriNeigh::unpack_restart(int nlocal, int nth)
     r0[nlocal][n] = extra[nlocal][m++];
   }
   vinter[nlocal] = extra[nlocal][m++];
+  wvolume[nlocal] = extra[nlocal][m++];
 }
 
 /* ----------------------------------------------------------------------
@@ -365,7 +509,7 @@ void FixPeriNeigh::unpack_restart(int nlocal, int nth)
 
 int FixPeriNeigh::maxsize_restart()
 {
-  return 2*maxpartner + 3;
+  return 2*maxpartner + 4;
 }
 
 /* ----------------------------------------------------------------------
@@ -374,5 +518,5 @@ int FixPeriNeigh::maxsize_restart()
 
 int FixPeriNeigh::size_restart(int nlocal)
 {
-  return 2*npartner[nlocal] + 3;
+  return 2*npartner[nlocal] + 4;
 }

src/PERI/fix_peri_neigh.h

@@ -26,6 +26,7 @@ namespace LAMMPS_NS {
 
 class FixPeriNeigh : public Fix {
   friend class PairPeriPMB;
+  friend class PairPeriLPS;
   friend class ComputeDamageAtom;
 
  public:
@@ -41,10 +42,15 @@ class FixPeriNeigh : public Fix {
   void copy_arrays(int, int);
   int pack_exchange(int, double *);
   int unpack_exchange(int, double *);
+  void write_restart(FILE *);
+  void restart(char *);
   int pack_restart(int, double *);
   void unpack_restart(int, int);
   int size_restart(int);
   int maxsize_restart();
+  int pack_comm(int, int *, double *, int, int *);
+  void unpack_comm(int, int, double *);
+
 
  private:
   int first;                 // flag for first time initialization
@@ -53,6 +59,7 @@ class FixPeriNeigh : public Fix {
   int **partner;             // neighs for each atom, stored as global IDs
   double **r0;               // initial distance to partners
   double *vinter;            // sum of vfrac for bonded neighbors
+  double *wvolume;           // weighted volume of particle
 
   class NeighList *list;
 };

src/PERI/pair_peri_lps.cpp

@@ -0,0 +1,745 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Mike Parks (SNL)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "values.h"
+#include "stdlib.h"
+#include "string.h"
+#include "pair_peri_lps.h"  			
+#include "atom.h"
+#include "domain.h"
+#include "lattice.h"
+#include "force.h"
+#include "update.h"
+#include "modify.h"
+#include "fix.h"
+#include "fix_peri_neigh.h"
+#include "comm.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "memory.h"
+#include "error.h"
+#include "update.h"                       
+
+using namespace LAMMPS_NS;
+
+#define MIN(a,b) ((a) < (b) ? (a) : (b))
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
+
+/* ---------------------------------------------------------------------- */
+
+PairPeriLPS::PairPeriLPS(LAMMPS *lmp) : Pair(lmp)			
+{
+  for (int i = 0; i < 6; i++) virial[i] = 0.0;
+
+  ifix_peri = -1;
+
+  nmax = 0;
+  s0_new = NULL;
+  theta = NULL;				
+
+  bulkmodulus = NULL;			
+  shearmodulus = NULL;			
+  s00 = alpha = NULL;                    
+  cut = NULL;
+
+  // set comm size needed by this Pair
+  // comm_reverse not needed
+  comm_forward = 1;  // For passing dilatation (theta)
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairPeriLPS::~PairPeriLPS()						
+{
+  if (ifix_peri >= 0) modify->delete_fix("PERI_NEIGH");
+
+  if (allocated) {
+    memory->destroy_2d_int_array(setflag);
+    memory->destroy_2d_double_array(cutsq);
+    memory->destroy_2d_double_array(bulkmodulus);		
+    memory->destroy_2d_double_array(shearmodulus); 		
+    memory->destroy_2d_double_array(s00); 		
+    memory->destroy_2d_double_array(alpha); 		
+    memory->destroy_2d_double_array(cut);
+    memory->sfree(theta);					
+    memory->sfree(s0_new);					
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairPeriLPS::compute(int eflag, int vflag)				
+{
+  int i,j,ii,jj,inum,jnum,itype,jtype;
+  double xtmp,ytmp,ztmp,delx,dely,delz;
+  double xtmp0,ytmp0,ztmp0,delx0,dely0,delz0,rsq0;
+  double rsq,r,dr,rk,evdwl,fpair,fbond;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  double d_ij,delta,stretch;                            
+
+  evdwl = 0.0;
+  if (eflag || vflag) ev_setup(eflag,vflag);
+  else evflag = vflag_fdotr = 0;
+
+  double **f = atom->f;
+  double **x = atom->x;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+
+  double *vfrac = atom->vfrac;
+  double *s0 = atom->s0;
+  double **x0 = atom->x0;
+  double **r0   = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0;
+  int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner;
+  int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner;
+  double *wvolume = ((FixPeriNeigh *) modify->fix[ifix_peri])->wvolume;
+
+  // lc = lattice constant
+  // init_style guarantees it's the same in x, y, and z
+
+  double lc = domain->lattice->xlattice;
+  double half_lc = 0.5*lc;
+  double vfrac_scale = 1.0;
+
+  // short-range forces
+
+  int nall = atom->nlocal + atom->nghost;
+  int newton_pair = force->newton_pair;
+  int nonperiodic = domain->nonperiodic;
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  double dt = update->dt;                   
+
+  // loop over neighbors of my atoms
+  // need minimg() for x0 difference since not ghosted
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    xtmp0 = x0[i][0];
+    ytmp0 = x0[i][1];
+    ztmp0 = x0[i][2];
+    itype = type[i];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      j %= nall;
+ 
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+
+      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);
+      rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0;
+      jtype = type[j];
+ 
+      r = sqrt(rsq);
+
+      // short-range interaction distance based on initial particle position
+      // 0.9 and 1.35 are constants
+
+      d_ij = MIN(0.9*sqrt(rsq0),1.35*lc);
+
+      // short-range contact forces
+      // 15 is constant taken from the EMU Theory Manual
+      // Silling, 12 May 2005, p 18
+
+      if (r < d_ij) {
+        dr = r - d_ij;
+
+        // kshort based upon short-range force constant of the bond-based theory used in PMB model
+        double kshort = ( 15.0 * 18.0 * bulkmodulus[itype][itype] ) /                                  
+                 ( 3.141592653589793 * cutsq[itype][jtype] * cutsq[itype][jtype] );             
+        rk = ( kshort * vfrac[j]) * (dr / sqrt(cutsq[itype][jtype]));				
+
+        if (r > 0.0) fpair = -(rk/r);
+        else fpair = 0.0;
+
+        f[i][0] += delx*fpair;
+        f[i][1] += dely*fpair;
+        f[i][2] += delz*fpair;
+        if (newton_pair || j < nlocal) {
+          f[j][0] -= delx*fpair;
+          f[j][1] -= dely*fpair;
+          f[j][2] -= delz*fpair;
+        }
+
+        if (eflag) evdwl = 0.5*rk*dr;
+	if (evflag) ev_tally(i,j,nlocal,newton_pair,
+			     evdwl,0.0,fpair,delx,dely,delz);
+      }
+    }
+  }
+
+  // grow bond forces array if necessary
+
+  if (atom->nmax > nmax) {
+    memory->sfree(s0_new);				
+    memory->sfree(theta);				
+    nmax = atom->nmax;
+    s0_new = (double *) memory->smalloc(nmax*sizeof(double),"pair:s0_new");		
+    theta = (double *) memory->smalloc(nmax*sizeof(double),"pair:theta");		
+  }
+
+  // Compute the dilatation on each particle				
+  compute_dilatation();						
+
+  // communicate dilatation (theta) of each particle	
+  comm->forward_comm_pair(this);
+  // communicate wighted volume (wvolume) upon every reneighbor
+  if (neighbor->ago == 0)
+    comm->forward_comm_fix(modify->fix[ifix_peri]);
+
+  // Volume-dependent part of the energy
+  for (i = 0; i < nlocal; i++) {   
+    itype = type[i];
+    if (eflag) {
+      if (eflag_global) eng_vdwl +=  0.5 * bulkmodulus[itype][itype] * (theta[i] * theta[i]);  
+      if (eflag_atom) eatom[i] += 0.5 * bulkmodulus[itype][itype] * (theta[i] * theta[i]);     
+    }
+  }
+
+  // loop over my particles and their partners
+  // partner list contains all bond partners, so I-J appears twice
+  // if bond already broken, skip this partner
+  // first = true if this is first neighbor of particle i
+
+  bool first;
+  double omega_minus, omega_plus;
+
+  for (i = 0; i < nlocal; i++) {
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    xtmp0 = x0[i][0];			
+    ytmp0 = x0[i][1];		
+    ztmp0 = x0[i][2];			
+    itype = type[i];
+    jnum = npartner[i];
+    first = true;
+
+    for (jj = 0; jj < jnum; jj++) {
+      if (partner[i][jj] == 0) continue;
+      j = atom->map(partner[i][jj]);
+
+      // check if lost a partner without first breaking bond
+
+      if (j < 0) {
+        partner[i][jj] = 0;
+        continue;
+      }
+
+      // compute force density, add to PD equation of motion
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      if (nonperiodic == 0) 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);   
+      jtype = type[j];
+      delta = sqrt(cutsq[itype][jtype]);
+      r = sqrt(rsq);
+      dr = r - r0[i][jj];
+
+      // avoid roundoff errors
+
+      if (fabs(dr) < 2.2204e-016) dr = 0.0;
+
+      // scale vfrac[j] if particle j near the horizon
+
+      if ((fabs(r0[i][jj] - delta)) <= half_lc)
+        vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + 
+	  (1.0 + ((delta - half_lc)/(2*half_lc) ) );
+      else vfrac_scale = 1.0;
+
+      omega_plus  = influence_function(-1.0*delx0,-1.0*dely0,-1.0*delz0);				
+      omega_minus = influence_function(delx0,dely0,delz0);		
+      rk = ( (3.0 * bulkmodulus[itype][itype]) - (5.0 * shearmodulus[itype][itype]) ) * vfrac[j] * vfrac_scale * ( (omega_plus * theta[i] / wvolume[i]) + ( omega_minus * theta[j] / wvolume[j] ) ) * r0[i][jj] ; 
+      rk +=  15.0 * ( shearmodulus[itype][itype] * vfrac[j] * vfrac_scale )  *  ( (omega_plus / wvolume[i]) + (omega_minus / wvolume[j]) ) * dr; 
+
+      if (r > 0.0) fbond = -(rk/r);
+      else fbond = 0.0;
+
+      f[i][0] += delx*fbond;
+      f[i][1] += dely*fbond;
+      f[i][2] += delz*fbond;
+
+      // since I-J is double counted, set newton off & use 1/2 factor and I,I 
+
+      double deviatoric_extension = dr - (theta[i]* r0[i][jj] / 3.0);
+      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,
+			   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
+
+      stretch = dr / r0[i][jj];
+      if (stretch > MIN(s0[i],s0[j])) partner[i][jj] = 0;
+
+      // update s0 for next timestep
+
+      if (first)
+         s0_new[i] = s00[itype][jtype] - (alpha[itype][jtype] * stretch);
+      else
+         s0_new[i] = MAX(s0_new[i],
+                         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];
+
+}
+
+/* ----------------------------------------------------------------------
+   allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairPeriLPS::allocate()					
+{
+  allocated = 1;
+  int n = atom->ntypes;
+
+  setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag");
+  for (int i = 1; i <= n; i++)
+    for (int j = i; j <= n; j++)
+      setflag[i][j] = 0;
+
+  cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq");
+  bulkmodulus = memory->create_2d_double_array(n+1,n+1,"pair:bulkmodulus");		
+  shearmodulus = memory->create_2d_double_array(n+1,n+1,"pair:shearmodulus");		
+  s00 = memory->create_2d_double_array(n+1,n+1,"pair:s00");		
+  alpha = memory->create_2d_double_array(n+1,n+1,"pair:alpha");		
+  cut = memory->create_2d_double_array(n+1,n+1,"pair:cut");
+
+}
+
+/* ----------------------------------------------------------------------
+   global settings
+------------------------------------------------------------------------- */
+
+void PairPeriLPS::settings(int narg, char **arg)		
+{
+  if (narg) error->all("Illegal pair_style command");
+}
+
+/* ----------------------------------------------------------------------
+   set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairPeriLPS::coeff(int narg, char **arg)			
+{
+  if (narg != 7) error->all("Incorrect args for pair coefficients"); 	
+  if (!allocated) allocate();
+
+  int ilo,ihi,jlo,jhi;
+  force->bounds(arg[0],atom->ntypes,ilo,ihi);
+  force->bounds(arg[1],atom->ntypes,jlo,jhi);
+
+  double bulkmodulus_one = atof(arg[2]);	
+  double shearmodulus_one = atof(arg[3]);		
+  double cut_one = atof(arg[4]);			
+  double s00_one = atof(arg[5]);			
+  double alpha_one = atof(arg[6]);			
+
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++) {
+    for (int j = MAX(jlo,i); j <= jhi; j++) {
+      bulkmodulus[i][j] = bulkmodulus_one;		
+      shearmodulus[i][j] = shearmodulus_one;		
+      cut[i][j] = cut_one;
+      s00[i][j] = s00_one;
+      alpha[i][j] = alpha_one;
+      setflag[i][j] = 1;
+      count++;
+    }
+  }
+
+  if (count == 0) error->all("Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+   init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairPeriLPS::init_one(int i, int j)			
+{
+  if (setflag[i][j] == 0) error->all("All pair coeffs are not set");
+
+  cutsq[i][j] = cut[i][j] * cut[i][j];
+  cutsq[j][i] = cutsq[i][j];
+
+  // set other j,i parameters
+
+  bulkmodulus[j][i] = bulkmodulus[i][j];		
+  shearmodulus[j][i] = shearmodulus[i][j];		
+  s00[j][i] = s00[i][j];		
+  alpha[j][i] = alpha[i][j];		
+
+  return cut[i][j];
+}
+
+/* ----------------------------------------------------------------------
+   init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairPeriLPS::init_style() 				
+{
+  // error checks
+
+  if (atom->map_style == 0) 
+    error->all("Pair peri requires an atom map, see atom_modify");
+
+  if (atom->style_match("peri") == 0)
+    error->all("Pair style peri_lps requires atom style peri");		
+
+  if (domain->lattice == NULL)
+    error->all("Pair peri requires a lattice be defined");
+  if (domain->lattice->xlattice != domain->lattice->ylattice || 
+      domain->lattice->xlattice != domain->lattice->zlattice || 
+      domain->lattice->ylattice != domain->lattice->zlattice)
+    error->all("Pair peri lattice is not identical in x, y, and z");
+
+  // if first init, create Fix needed for storing fixed neighbors
+
+  if (ifix_peri == -1) {
+    char **fixarg = new char*[3];
+    fixarg[0] = (char *) "PERI_NEIGH";
+    fixarg[1] = (char *) "all";
+    fixarg[2] = (char *) "PERI_NEIGH";
+    modify->add_fix(3,fixarg);
+    delete [] fixarg;
+  }
+
+  // find associated PERI_NEIGH fix that must exist
+  // could have changed locations in fix list since created
+
+  for (int i = 0; i < modify->nfix; i++)
+    if (strcmp(modify->fix[i]->style,"PERI_NEIGH") == 0) ifix_peri = i;
+  if (ifix_peri == -1) error->all("Fix peri neigh does not exist");
+
+  int irequest = neighbor->request(this);
+}
+
+/* ----------------------------------------------------------------------
+  proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairPeriLPS::write_restart(FILE *fp)		
+{
+  int i,j;
+  for (i = 1; i <= atom->ntypes; i++)
+    for (j = i; j <= atom->ntypes; j++) {
+      fwrite(&setflag[i][j],sizeof(int),1,fp);
+      if (setflag[i][j]) {
+	fwrite(&bulkmodulus[i][j],sizeof(double),1,fp); 	
+	fwrite(&shearmodulus[i][j],sizeof(double),1,fp);	
+	fwrite(&s00[i][j],sizeof(double),1,fp);	
+	fwrite(&alpha[i][j],sizeof(double),1,fp);	
+	fwrite(&cut[i][j],sizeof(double),1,fp);
+      }
+    }
+}
+
+/* ----------------------------------------------------------------------
+  proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairPeriLPS::read_restart(FILE *fp)		
+{
+  allocate();
+
+  int i,j;
+  int me = comm->me;
+  for (i = 1; i <= atom->ntypes; i++)
+    for (j = i; j <= atom->ntypes; j++) {
+      if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
+      MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
+      if (setflag[i][j]) {
+	if (me == 0) {
+	  fread(&bulkmodulus[i][j],sizeof(double),1,fp);		
+	  fread(&shearmodulus[i][j],sizeof(double),1,fp);		
+	  fread(&s00[i][j],sizeof(double),1,fp);                        
+	  fread(&alpha[i][j],sizeof(double),1,fp);                        
+	  fread(&cut[i][j],sizeof(double),1,fp);
+	}
+	MPI_Bcast(&bulkmodulus[i][j],1,MPI_DOUBLE,0,world);		
+	MPI_Bcast(&shearmodulus[i][j],1,MPI_DOUBLE,0,world);		
+	MPI_Bcast(&s00[i][j],1,MPI_DOUBLE,0,world);                      
+	MPI_Bcast(&alpha[i][j],1,MPI_DOUBLE,0,world);                      
+	MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
+      }
+    }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairPeriLPS::single(int i, int j, int itype, int jtype, double rsq,  		
+			   double factor_coul, double factor_lj,
+			   double &fforce)
+{
+  double delx0,dely0,delz0,rsq0;
+  double d_ij,r,dr,rk,vfrac_scale;
+
+  double *vfrac = atom->vfrac;
+  double **x0 = atom->x0;
+  double **r0   = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0;
+  int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner;
+  int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner;
+  double *wvolume = ((FixPeriNeigh *) modify->fix[ifix_peri])->wvolume;
+
+  double lc = domain->lattice->xlattice;
+  double half_lc = 0.5*lc;
+
+  double kshort;				
+ 
+  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);
+  rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0;
+
+  d_ij = MIN(0.9*sqrt(rsq0),1.35*lc);
+  r = sqrt(rsq);
+
+  double energy = 0.0;
+  fforce = 0.0;
+
+  if (r < d_ij) {
+    dr = r - d_ij;
+    // kshort resembles the short-range force constant of the bond-based theory in 3-D !
+    kshort = (15.0 * 18.0 * bulkmodulus[itype][itype]) / 				 	
+             ( 3.141592653589793 * cutsq[itype][jtype] * cutsq[itype][jtype] ); 	 
+    rk = ( kshort * vfrac[j]) * (dr / sqrt(cutsq[itype][jtype]));		 	 
+    if (r > 0.0) fforce += -(rk/r);
+    energy += 0.5*rk*dr;	
+  }
+ 
+   if (atom->nmax > nmax) {
+    memory->sfree(theta);
+    nmax = atom->nmax;
+    theta = (double *) memory->smalloc(nmax*sizeof(double),"pair:theta");
+  }
+
+  // Compute the dilatation on each particle			
+  compute_dilatation();					
+
+  // communicate dilatation (theta) of each particle
+  comm->forward_comm_pair(this);
+  // communicate wighted volume (wvolume) upon every reneighbor
+  if (neighbor->ago == 0)
+    comm->forward_comm_fix(modify->fix[ifix_peri]);
+
+  double omega_plus, omega_minus;               
+
+  int jnum = npartner[i];
+  for (int jj = 0; jj < jnum; jj++) {
+    if (partner[i][jj] == 0) continue;
+    if (j < 0) continue;
+    if (j == atom->map(partner[i][jj])) {
+      dr = r - r0[i][jj];
+      if (fabs(dr) < 2.2204e-016) dr = 0.0;
+
+      // scale vfrac[j] if particle j near the horizon
+
+      if ( (fabs(r0[i][jj] - sqrt(cutsq[itype][jtype]))) <= half_lc)
+	vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + 
+	  (1.0 + ((sqrt(cutsq[itype][jtype]) - half_lc)/(2*half_lc)));
+      else vfrac_scale = 1.0;
+
+      omega_plus  = influence_function(-1.0*delx0,-1.0*dely0,-1.0*delz0);			
+      omega_minus = influence_function(delx0,dely0,delz0);	
+      rk = (3.0* bulkmodulus[itype][itype] -5.0 *  shearmodulus[itype][itype] )* vfrac[j] * vfrac_scale  * ( (omega_plus * theta[i] / wvolume[i]) + ( omega_minus * theta[j] / wvolume[j] ) ) * r0[i][jj] ;
+      rk +=  15.0 * ( shearmodulus[itype][itype] * vfrac[j] * vfrac_scale )  *   ( (omega_plus / wvolume[i]) + (omega_minus / wvolume[j]) ) * dr; 
+
+      if (r > 0.0) fforce += -(rk/r);
+        energy += 0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) * omega_plus * (  dr - theta[i]* r0[i][jj] / 3.0 ) * (  dr - theta[i]* r0[i][jj] / 3.0 ) * vfrac[j] * vfrac_scale; 
+
+     }
+  }
+
+  return energy;
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based arrays 
+------------------------------------------------------------------------- */
+
+double PairPeriLPS::memory_usage()
+{
+  double bytes = 2 * nmax * sizeof(double);
+  return bytes;
+}
+
+/* ----------------------------------------------------------------------
+   Influence function definition
+------------------------------------------------------------------------- */
+
+double PairPeriLPS::influence_function(double xi_x, double xi_y, double xi_z)
+{
+  double r = sqrt(xi_x*xi_x + xi_y*xi_y + xi_z*xi_z);
+  double omega;
+
+  if (fabs(r) < 2.2204e-016) error->one("pair_peri_lps - influence_function: Divide by 0");
+
+  omega = 1.0/r;
+
+  return omega;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairPeriLPS::compute_dilatation()
+{
+  int i,j,jj,jnum,itype,jtype;
+  double xtmp,ytmp,ztmp,delx,dely,delz;
+  double xtmp0,ytmp0,ztmp0,delx0,dely0,delz0;
+  double rsq,r,dr;
+  double delta;
+
+  double **x = atom->x;
+  int *type = atom->type;
+  double **x0 = atom->x0;
+  int nlocal = atom->nlocal;
+  double *vfrac = atom->vfrac;
+  double vfrac_scale = 1.0;
+  
+  double lc = domain->lattice->xlattice;
+  double half_lc = 0.5*lc;
+
+  double **r0   = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0;
+  int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner;
+  int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner;
+  double *wvolume = ((FixPeriNeigh *) modify->fix[ifix_peri])->wvolume;
+
+  int nonperiodic = domain->nonperiodic;
+
+  // compute the dilatation theta
+  // loop over my particles
+  for (i = 0; i < nlocal; i++) {
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    xtmp0 = x0[i][0];
+    ytmp0 = x0[i][1];
+    ztmp0 = x0[i][2];
+    jnum = npartner[i];
+    theta[i] = 0.0;
+    itype = type[i];
+
+    // loop over partners of particle i
+    for (jj = 0; jj < jnum; jj++) {
+
+      // if bond already broken, skip this partner
+      if (partner[i][jj] == 0) continue;
+
+      // Look up local index of this partner particle
+      j = atom->map(partner[i][jj]);
+
+      // Skip if particle is "lost"
+      if (j < 0) continue;
+
+      // Compute force density and add to PD equation of motion
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      if (nonperiodic == 0) 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);
+
+      r = sqrt(rsq);
+      dr = r - r0[i][jj];
+      if (fabs(dr) < 2.2204e-016) dr = 0.0;
+
+      jtype = type[j];
+      delta = sqrt(cutsq[itype][jtype]);
+
+      // scale vfrac[j] if particle j near the horizon
+
+      if ((fabs(r0[i][jj] - delta)) <= half_lc)
+        vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) +
+          (1.0 + ((delta - half_lc)/(2*half_lc) ) );
+      else vfrac_scale = 1.0;
+
+      theta[i] += influence_function(delx0, dely0, delz0) * r0[i][jj] * dr * vfrac[j] * vfrac_scale;
+
+    } // end loop over all neighbors jj
+
+    // If wvolume[i] is zero, then particle i has no bonds. Therefore, the dilatation is set to 0. 
+    if (wvolume[i] != 0.0) theta[i] = (3.0/wvolume[i]) * theta[i];
+    else theta[i] = 0;
+
+  } // end loop over all particles i
+
+}
+
+/* ----------------------------------------------------------------------
+   communication routines
+   ---------------------------------------------------------------------- */
+ 
+int PairPeriLPS::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
+{
+
+  int i,j,m;
+ 
+  m = 0;
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    buf[m++] = theta[j];
+  }
+  return 1;
+}
+ 
+/* ---------------------------------------------------------------------- */
+ 
+void PairPeriLPS::unpack_comm(int n, int first, double *buf)
+{
+  int i,m,last;
+ 
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) {
+    theta[i] = buf[m++];
+  }
+}
+
+

src/PERI/pair_peri_lps.h

@@ -0,0 +1,65 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS                    
+
+PairStyle(peri/lps,PairPeriLPS)      
+
+#else                                 
+
+#ifndef LMP_PAIR_PERI_LPS_H          
+#define LMP_PAIR_PERI_LPS_H          
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairPeriLPS : public Pair {    
+ public:
+  PairPeriLPS(class LAMMPS *);       
+  ~PairPeriLPS();                                   
+  int pack_comm(int, int *, double *, int, int *);  
+  void unpack_comm(int, int, double *);             
+
+  void compute(int, int);
+  void settings(int, char **);
+  void coeff(int, char **);
+  double init_one(int, int);
+  void init_style();
+  void write_restart(FILE *);
+  void read_restart(FILE *);
+  void write_restart_settings(FILE *) {}
+  void read_restart_settings(FILE *) {}
+  double single(int, int, int, int, double, double, double, double &);
+  double memory_usage();
+  double influence_function(double, double, double);
+  void compute_dilatation();   
+
+ protected:
+  int ifix_peri;
+  double **bulkmodulus;                
+  double **shearmodulus;               
+  double **s00, **alpha;                  
+  double **cut;
+ 
+  double *s0_new;                  
+  double *theta;                      
+  int nmax;
+
+  void allocate();
+};
+
+}
+
+#endif
+#endif

src/PERI/pair_peri_pmb.cpp

@@ -31,6 +31,7 @@
 #include "comm.h"
 #include "neighbor.h"
 #include "neigh_list.h"
+#include "neigh_request.h"
 #include "memory.h"
 #include "error.h"
 
@@ -399,6 +400,7 @@ void PairPeriPMB::init_style()
 
   for (int i = 0; i < modify->nfix; i++)
     if (strcmp(modify->fix[i]->style,"PERI_NEIGH") == 0) ifix_peri = i;
+  if (ifix_peri == -1) error->all("Fix peri neigh does not exist");
 
   int irequest = neighbor->request(this);
 }