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

src/PERI/compute_damage_atom.cpp

@@ -100,7 +100,6 @@ void ComputeDamageAtom::compute_peratom()
     if (mask[i] & groupbit) {
       jnum = npartner[i];
       damage_temp = 0.0;
-
       for (jj = 0; jj < jnum; jj++) {
         if (partner[i][jj] == 0) continue;
         
@@ -112,11 +111,11 @@ void ComputeDamageAtom::compute_peratom()
         
         damage_temp += vfrac[j];
       }
-    }
-    else damage_temp = vinter[i];
 
       if (vinter[i] != 0.0) damage[i] = 1.0 - damage_temp/vinter[i];
       else damage[i] = 0.0;
+
+    } else damage[i] = 0.0;
   }
 }
 

src/PERI/compute_dilatation_atom.cpp

@@ -0,0 +1,127 @@
+/* ----------------------------------------------------------------------
+   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: Rezwanur Rahman, John Foster (UTSA)
+------------------------------------------------------------------------- */
+
+#include "string.h"
+#include "compute_dilatation_atom.h"
+#include "atom.h"
+#include "update.h"
+#include "modify.h"
+#include "comm.h"
+#include "force.h"
+#include "pair.h"
+#include "pair_peri_lps.h"
+#include "pair_peri_pmb.h"
+#include "pair_peri_ves.h"
+#include "pair_peri_eps.h"
+#include "fix_peri_neigh.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeDilatationAtom::
+ComputeDilatationAtom(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg)
+{
+  if (narg != 3) error->all(FLERR,"Illegal compute Dilatation/atom command");
+
+  peratom_flag = 1;
+  size_peratom_cols = 0;
+
+  nmax = 0;
+  dilatation = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeDilatationAtom::~ComputeDilatationAtom()
+{
+  memory->destroy(dilatation);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeDilatationAtom::init()
+{
+  int count = 0;
+  for (int i = 0; i < modify->ncompute; i++)
+    if (strcmp(modify->compute[i]->style,"dilatation/peri") == 0) count++;
+  if (count > 1 && comm->me == 0)
+    error->warning(FLERR,"More than one compute dilatation/atom");
+
+  // check PD pair style
+
+  isPMB = isLPS = isVES = isEPS = 0;
+  if (force->pair_match("peri/pmb",1)) isPMB = 1;
+  if (force->pair_match("peri/lps",1)) isLPS = 1;
+  if (force->pair_match("peri/ves",1)) isVES = 1;
+  if (force->pair_match("peri/eps",1)) isEPS = 1;
+  
+  if (isPMB) 
+    error->all(FLERR,"Compute dilatation/atom cannot be used "
+               "with this pair style");
+
+  // find associated PERI_NEIGH fix that must exist
+
+  int ifix_peri = -1;
+  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(FLERR,"Compute dilatation/atom requires Peridynamic pair style");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeDilatationAtom::compute_peratom()
+{
+  invoked_peratom = update->ntimestep;
+
+  // grow dilatation array if necessary
+
+  if (atom->nlocal > nmax) {
+    memory->destroy(dilatation);
+    nmax = atom->nmax;
+    memory->create(dilatation,nmax,"dilatation/atom:dilatation");
+    vector_atom = dilatation;
+  }
+
+  // extract dilatation for each atom in group
+
+  double *theta;
+  Pair *anypair = force->pair_match("peri",0);
+  if (isLPS) theta = ((PairPeriLPS *) anypair)->theta;
+  if (isVES) theta = ((PairPeriVES *) anypair)->theta;
+  if (isEPS) theta = ((PairPeriEPS *) anypair)->theta;
+
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) dilatation[i] = theta[i];
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double ComputeDilatationAtom::memory_usage()
+{
+  double bytes = nmax * sizeof(double);
+  return bytes;
+}

src/PERI/compute_dilatation_atom.h

@@ -0,0 +1,48 @@
+/* ----------------------------------------------------------------------
+   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 COMPUTE_CLASS
+
+ComputeStyle(dilatation/atom,ComputeDilatationAtom)
+
+#else
+
+#ifndef LMP_COMPUTE_DILATATION_ATOM_H
+#define LMP_COMPUTE_DILATATION_ATOM_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeDilatationAtom : public Compute {
+  friend class PairPeriPMB;
+  friend class PairPeriLPS;
+  friend class PairPeriVES;   
+  friend class PairPeriEPS;   
+ public:
+  ComputeDilatationAtom(class LAMMPS *, int, char **);
+  ~ComputeDilatationAtom();
+  void init();
+  void compute_peratom();
+  double memory_usage();
+
+ private:
+  int nmax;
+  double *dilatation;
+  int isPMB,isLPS,isVES,isEPS;
+};
+
+}
+
+#endif
+#endif

src/PERI/compute_plasticity_atom.cpp

@@ -0,0 +1,111 @@
+/* ----------------------------------------------------------------------
+   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: Rezwanur Rahman, John Foster (UTSA)
+------------------------------------------------------------------------- */
+
+#include "string.h"
+#include "compute_plasticity_atom.h"
+#include "atom.h"
+#include "update.h"
+#include "modify.h"
+#include "comm.h"
+#include "force.h"
+#include "pair_peri_pmb.h"
+#include "fix_peri_neigh.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputePlasticityAtom::
+ComputePlasticityAtom(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg)
+{
+  if (narg != 3) error->all(FLERR,"Illegal compute plasticity/atom command");
+  
+  if (!force->pair_match("peri/eps",1))
+    error->all(FLERR,"Compute plasticity/atom cannot be used "
+               "with this pair style");
+
+  peratom_flag = 1;
+  size_peratom_cols = 0;
+
+  nmax = 0;
+  plasticity = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputePlasticityAtom::~ComputePlasticityAtom()
+{
+  memory->destroy(plasticity);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePlasticityAtom::init()
+{
+  int count = 0;
+  for (int i = 0; i < modify->ncompute; i++)
+    if (strcmp(modify->compute[i]->style,"plasticity/peri") == 0) count++;
+  if (count > 1 && comm->me == 0)
+    error->warning(FLERR,"More than one compute plasticity/atom");
+
+  // find associated PERI_NEIGH fix that must exist
+
+  ifix_peri = -1;
+  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(FLERR,"Compute plasticity/atom requires Peridynamic pair style");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePlasticityAtom::compute_peratom()
+{
+  invoked_peratom = update->ntimestep;
+
+  // grow damage array if necessary
+
+  if (atom->nlocal > nmax) {
+    memory->destroy(plasticity);
+    nmax = atom->nmax;
+    memory->create(plasticity,nmax,"plasticity/atom:plasticity");
+    vector_atom = plasticity;
+  }
+
+  // extract plasticity for each atom in group
+
+  double *lambdaValue = ((FixPeriNeigh *) modify->fix[ifix_peri])->lambdaValue;
+
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) plasticity[i] = lambdaValue[i];
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double ComputePlasticityAtom::memory_usage()
+{
+  double bytes = nmax * sizeof(double);
+  return bytes;
+}

src/PERI/compute_plasticity_atom.h

@@ -0,0 +1,44 @@
+/* ----------------------------------------------------------------------
+   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 COMPUTE_CLASS
+
+ComputeStyle(plasticity/atom,ComputePlasticityAtom)
+
+#else
+
+#ifndef LMP_COMPUTE_PLASTICITY_ATOM_H
+#define LMP_COMPUTE_PLASTICITY_ATOM_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputePlasticityAtom : public Compute {
+ public:
+  ComputePlasticityAtom(class LAMMPS *, int, char **);
+  ~ComputePlasticityAtom();
+  void init();
+  void compute_peratom();
+  double memory_usage();
+
+ private:
+  int nmax;
+  double *plasticity;
+  int ifix_peri;
+};
+
+}
+
+#endif
+#endif

src/PERI/fix_peri_neigh.cpp

@@ -20,6 +20,7 @@
 #include "pair_peri_pmb.h"
 #include "pair_peri_lps.h"
 #include "pair_peri_ves.h"
+#include "pair_peri_eps.h"
 #include "atom.h"
 #include "domain.h"
 #include "force.h"
@@ -41,10 +42,11 @@ using namespace FixConst;
 FixPeriNeigh::FixPeriNeigh(LAMMPS *lmp,int narg, char **arg) :
   Fix(lmp, narg, arg)
 {
-  isPMB = isLPS = isVES = 0;
+  isPMB = isLPS = isVES = isEPS = 0;
   if (force->pair_match("peri/pmb",1)) isPMB = 1;
   if (force->pair_match("peri/lps",1)) isLPS = 1;
   if (force->pair_match("peri/ves",1)) isVES = 1;
+  if (force->pair_match("peri/eps",1)) isEPS = 1;
   
   restart_global = 1;
   restart_peratom = 1;
@@ -59,6 +61,8 @@ FixPeriNeigh::FixPeriNeigh(LAMMPS *lmp,int narg, char **arg) :
   partner = NULL;
   deviatorextention = NULL;
   deviatorBackextention = NULL;
+  deviatorPlasticextension = NULL;
+  lambdaValue = NULL;
   r0 = NULL;
   vinter = NULL;
   wvolume = NULL;
@@ -92,6 +96,8 @@ FixPeriNeigh::~FixPeriNeigh()
   memory->destroy(partner);
   memory->destroy(deviatorextention);
   memory->destroy(deviatorBackextention);
+  memory->destroy(deviatorPlasticextension);
+  memory->destroy(lambdaValue);
   memory->destroy(r0);
   memory->destroy(vinter);
   memory->destroy(wvolume);
@@ -218,6 +224,8 @@ void FixPeriNeigh::setup(int vflag)
   memory->destroy(partner);
   memory->destroy(deviatorextention);
   memory->destroy(deviatorBackextention);
+  memory->destroy(deviatorPlasticextension);
+  memory->destroy(lambdaValue);  
   memory->destroy(r0);
   memory->destroy(npartner);
 
@@ -225,6 +233,8 @@ void FixPeriNeigh::setup(int vflag)
   partner = NULL;
   deviatorextention = NULL;
   deviatorBackextention = NULL;
+  deviatorPlasticextension = NULL;
+  lambdaValue = NULL;
   r0 = NULL;   
   grow_arrays(atom->nmax);
 
@@ -235,6 +245,7 @@ void FixPeriNeigh::setup(int vflag)
     npartner[i] = 0;
     vinter[i] = 0.0;
     wvolume[i] = 0.0;
+    if (isEPS) lambdaValue[i] = 0.0;
   }
 
   for (ii = 0; ii < inum; ii++) {
@@ -261,6 +272,8 @@ void FixPeriNeigh::setup(int vflag)
         if (isVES)
           deviatorextention[i][npartner[i]] = 
             deviatorBackextention[i][npartner[i]] = 0.0;
+        if (isEPS)
+           deviatorPlasticextension[i][npartner[i]] = 0.0;
         r0[i][npartner[i]] = sqrt(rsq);   
         npartner[i]++;
         vinter[i] += vfrac[j];
@@ -292,6 +305,7 @@ void FixPeriNeigh::setup(int vflag)
   PairPeriLPS *pairlps = static_cast<PairPeriLPS*>(anypair);
   PairPeriPMB *pairpmb = static_cast<PairPeriPMB*>(anypair);
   PairPeriVES *pairves = static_cast<PairPeriVES*>(anypair);
+  PairPeriEPS *paireps = static_cast<PairPeriEPS*>(anypair);
 
   for (i = 0; i < nlocal; i++) {
     double xtmp0 = x0[i][0];
@@ -319,6 +333,7 @@ void FixPeriNeigh::setup(int vflag)
       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];
@@ -332,7 +347,6 @@ void FixPeriNeigh::setup(int vflag)
       else vfrac_scale = 1.0;
 
       // for PMB, influence = 1.0, otherwise invoke influence function
-
       if (isPMB) 
         wvolume[i] += 1.0 * rsq0 * vfrac[j] * vfrac_scale; 
       else if (isLPS)
@@ -341,6 +355,9 @@ void FixPeriNeigh::setup(int vflag)
       else if (isVES)
         wvolume[i] += pairves->influence_function(delx0,dely0,delz0) *
           rsq0 * vfrac[j] * vfrac_scale;
+      else if (isEPS)
+        wvolume[i] += paireps->influence_function(delx0,dely0,delz0) *
+          rsq0 * vfrac[j] * vfrac_scale;    
     }
   }
 
@@ -377,12 +394,16 @@ double FixPeriNeigh::memory_usage()
 { 
   int nmax = atom->nmax;
   int bytes = nmax * sizeof(int);
-  bytes += nmax*maxpartner * sizeof(int);
+  bytes += nmax*maxpartner * sizeof(tagint);
   bytes += nmax*maxpartner * sizeof(double);
   if (isVES) {
     bytes += nmax*maxpartner * sizeof(double);
     bytes += nmax*maxpartner * sizeof(double);
   }  
+  if (isEPS) {
+    bytes += nmax*maxpartner * sizeof(double);
+    bytes += nmax * sizeof(double);
+  }  
   bytes += nmax * sizeof(double);
   bytes += nmax * sizeof(double);
   return bytes; 
@@ -402,7 +423,10 @@ void FixPeriNeigh::grow_arrays(int nmax)
      memory->grow(deviatorBackextention,nmax,maxpartner,
                   "peri_neigh:deviatorBackextention");
    }
+   if (isEPS) memory->grow(deviatorPlasticextension,nmax,maxpartner,
+                           "peri_neigh:deviatorPlasticextension");
    memory->grow(r0,nmax,maxpartner,"peri_neigh:r0");
+   if (isEPS) memory->grow(lambdaValue,nmax,"peri_neigh:lambdaValue");   
    memory->grow(vinter,nmax,"peri_neigh:vinter");
    memory->grow(wvolume,nmax,"peri_neigh:wvolume");
 }
@@ -420,8 +444,11 @@ void FixPeriNeigh::copy_arrays(int i, int j, int delflag)
       deviatorextention[j][m] = deviatorextention[i][m];
       deviatorBackextention[j][m] = deviatorBackextention[i][m];
     }  
+    if (isEPS)
+      deviatorPlasticextension[j][m] = deviatorPlasticextension[i][m];
     r0[j][m] = r0[i][m];
   }
+  if (isEPS) lambdaValue[j] = lambdaValue[i];
   vinter[j] = vinter[i];
   wvolume[j] = wvolume[i];
 }
@@ -443,10 +470,13 @@ int FixPeriNeigh::pack_exchange(int i, double *buf)
       buf[m++] = deviatorextention[i][n];
       buf[m++] = deviatorBackextention[i][n];
     } 
+    if (isEPS) buf[m++] = deviatorPlasticextension[i][n];
     buf[m++] = r0[i][n];
   }
   if (isVES) buf[0] = m/4;
+  else if (isEPS) buf[0] = m/3;
   else buf[0] = m/2;
+  if (isEPS) buf[m++] = lambdaValue[i]; 
   buf[m++] = vinter[i];
   buf[m++] = wvolume[i];
   return m;
@@ -466,8 +496,10 @@ int FixPeriNeigh::unpack_exchange(int nlocal, double *buf)
       deviatorextention[nlocal][n] = buf[m++];
       deviatorBackextention[nlocal][n] = buf[m++];
     }
+    if (isEPS) deviatorPlasticextension[nlocal][n] = buf[m++];
     r0[nlocal][n] = buf[m++];     
   }
+  if (isEPS) lambdaValue[nlocal] = buf[m++];
   vinter[nlocal] = buf[m++];
   wvolume[nlocal] = buf[m++];
   return m;
@@ -543,8 +575,8 @@ int FixPeriNeigh::pack_restart(int i, double *buf)
 {
   int m = 0;
   if (isVES) buf[m++] = 4*npartner[i] + 4;
+  else if (isEPS) buf[m++] = 3*npartner[i] + 5;
   else buf[m++] = 2*npartner[i] + 4;
-
   buf[m++] = npartner[i];
   for (int n = 0; n < npartner[i]; n++) {
     buf[m++] = partner[i][n];
@@ -552,8 +584,10 @@ int FixPeriNeigh::pack_restart(int i, double *buf)
       buf[m++] = deviatorextention[i][n];
       buf[m++] = deviatorBackextention[i][n];
     }  
+    if (isEPS) buf[m++] = deviatorPlasticextension[i][n];
     buf[m++] = r0[i][n];
   }
+  if (isEPS) buf[m++] = lambdaValue[i];
   buf[m++] = vinter[i];
   buf[m++] = wvolume[i];
   return m;  
@@ -581,8 +615,10 @@ void FixPeriNeigh::unpack_restart(int nlocal, int nth)
       deviatorextention[nlocal][n] = extra[nlocal][m++];
       deviatorBackextention[nlocal][n] = extra[nlocal][m++];
     }  
+    if (isEPS) deviatorPlasticextension[nlocal][n] = extra[nlocal][m++];
     r0[nlocal][n] = extra[nlocal][m++];
   }
+  if (isEPS) lambdaValue[nlocal] = extra[nlocal][m++];
   vinter[nlocal] = extra[nlocal][m++];
   wvolume[nlocal] = extra[nlocal][m++];  
 }
@@ -594,6 +630,7 @@ void FixPeriNeigh::unpack_restart(int nlocal, int nth)
 int FixPeriNeigh::maxsize_restart()
 {
   if (isVES) return 4*maxpartner + 4;
+  if (isEPS) return 3*maxpartner + 5;
   return 2*maxpartner + 4;  
 }
 
@@ -604,5 +641,6 @@ int FixPeriNeigh::maxsize_restart()
 int FixPeriNeigh::size_restart(int nlocal)
 {
   if (isVES) return 4*npartner[nlocal] + 4;
+  if (isEPS) return 3*npartner[nlocal] + 5;
   return 2*npartner[nlocal] + 4; 
 }

src/PERI/fix_peri_neigh.h

@@ -28,9 +28,11 @@ class FixPeriNeigh : public Fix {
   friend class PairPeriPMB;
   friend class PairPeriPMBOMP;
   friend class PairPeriLPS;
-  friend class PairPeriVES;   //NEW
+  friend class PairPeriVES;
+  friend class PairPeriEPS;
   friend class PairPeriLPSOMP;
   friend class ComputeDamageAtom;
+  friend class ComputePlasticityAtom;
 
  public:
   FixPeriNeigh(class LAMMPS *,int, char **);
@@ -63,14 +65,14 @@ class FixPeriNeigh : public Fix {
   tagint **partner;          // neighs for each atom, stored as global IDs
   double **deviatorextention; // Deviatoric extention     
   double **deviatorBackextention; // Deviatoric back extention 
+  double **deviatorPlasticextension; // Deviatoric plastic extension 
+  double *lambdaValue;
   double **r0;               // initial distance to partners
-  double **r1;               // Instanteneous distance to partners *** NEW ***
+  double **r1;               // instanteneous distance to partners
-  double *thetaOld;          // Dilatation Old one
+  double *thetaValue;        // dilatation 
   double *vinter;            // sum of vfrac for bonded neighbors
   double *wvolume;           // weighted volume of particle
-  int isPMB;
+  int isPMB,isLPS,isVES,isEPS;  // which flavor of PD
-  int isLPS;
-  int isVES;
 
   class NeighList *list;
 };

src/PERI/pair_peri_eps.cpp

@@ -0,0 +1,843 @@
+/* ----------------------------------------------------------------------
+   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: Rezwanur Rahman, John Foster (UTSA)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "stdlib.h"
+#include "string.h"
+#include "pair_peri_eps.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;
+
+/* ---------------------------------------------------------------------- */
+
+PairPeriEPS::PairPeriEPS(LAMMPS *lmp) : Pair(lmp)
+{
+  for (int i = 0; i < 6; i++) virial[i] = 0.0;
+  no_virial_fdotr_compute = 1;
+  single_enable = 0;
+
+  ifix_peri = -1;
+
+  nmax = 0;
+  s0_new = NULL;
+  theta = NULL;
+
+  bulkmodulus = NULL;
+  shearmodulus = NULL;
+  s00 = alpha = NULL;
+  cut = NULL;
+  m_yieldstress = NULL;
+
+  // set comm size needed by this Pair
+  // comm_reverse not needed
+
+  comm_forward = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairPeriEPS::~PairPeriEPS()
+{
+  if (ifix_peri >= 0) modify->delete_fix("PERI_NEIGH");
+
+  if (allocated) {
+    memory->destroy(setflag);
+    memory->destroy(cutsq);
+    memory->destroy(bulkmodulus);
+    memory->destroy(shearmodulus);
+    memory->destroy(s00);
+    memory->destroy(alpha);
+    memory->destroy(cut);
+    memory->destroy(m_yieldstress);
+    memory->destroy(theta);
+    memory->destroy(s0_new);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairPeriEPS::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,dr1,rk,rkNew,evdwl,fpair,fbond;
+  double ed,fbondElastoPlastic,fbondFinal;
+  double deltalambda,edpNp1;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  double d_ij,delta,stretch;
+
+  evdwl = 0.0;
+  if (eflag || vflag) ev_setup(eflag,vflag);
+  else evflag = vflag_fdotr = eflag_global = eflag_atom = 0;
+
+  double **f = atom->f;
+  double **x = atom->x;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+
+  double timestepsize = update->dt;
+  double *vfrac = atom->vfrac;
+  double *s0 = atom->s0;
+  double **x0 = atom->x0;
+  double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0;
+  double **deviatorPlasticextension = 
+    ((FixPeriNeigh *) modify->fix[ifix_peri])->deviatorPlasticextension;
+  tagint **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner;
+  int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner;
+  double *wvolume = ((FixPeriNeigh *) modify->fix[ifix_peri])->wvolume;
+  double *lambdaValue = ((FixPeriNeigh *) modify->fix[ifix_peri])->lambdaValue;
+
+  // 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 newton_pair = force->newton_pair;
+  int periodic = domain->xperiodic || domain->yperiodic || domain->zperiodic;
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  // 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 &= NEIGHMASK;
+
+      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 (periodic) 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 / cut[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*vfrac[i],delx,dely,delz);
+      }
+    }
+  }
+
+  // grow bond forces array if necessary
+  
+  int  maxpartner = 0;
+  for (i = 0; i < nlocal; i++) maxpartner = MAX(maxpartner,npartner[i]);
+  
+
+  if (atom->nmax > nmax) {
+    memory->destroy(s0_new);
+    memory->destroy(theta);
+    nmax = atom->nmax;
+    memory->create(s0_new,nmax,"pair:s0_new");
+    memory->create(theta,nmax,"pair:theta");
+    
+  }
+  
+  // ******** temp array to store Plastic extension *********** ///
+  double deviatorPlasticExtTemp[nlocal][maxpartner];
+  for (int ii = 0; ii < nlocal; ii++) {
+      for (int kk = 0; kk < maxpartner; kk++) {
+           deviatorPlasticExtTemp[ii][kk] = 0.0;
+      }
+  }     
+  // ******** temp array to store Plastic extension *********** ///
+                 
+   
+
+  // compute the dilatation on each particle
+
+  compute_dilatation();
+
+  // communicate dilatation (theta) of each particle
+  comm->forward_comm_pair(this);
+
+  // communicate weighted volume (wvolume) upon every reneighbor
+
+  if (neighbor->ago == 0)
+    comm->forward_comm_fix(modify->fix[ifix_peri]);
+
+  // volume-dependent part of the energy
+
+  if (eflag) {
+    for (i = 0; i < nlocal; i++) {
+      itype = type[i];
+      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, omega;
+
+  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;
+        
+
+    double yieldStress = m_yieldstress[itype][itype];
+    double horizon = cut[itype][itype];
+    double tdnorm = compute_DeviatoricForceStateNorm(i);
+    double pointwiseYieldvalue = 25.0 * yieldStress * 
+                            yieldStress / 8 / M_PI / pow(horizon,5);
+                                              
+        
+    double fsurf = (tdnorm * tdnorm)/2 - pointwiseYieldvalue;
+    bool elastic = true;
+
+    double alphavalue = (15 * shearmodulus[itype][itype]) /wvolume[i];
+       
+      
+    if (fsurf>0) {
+      elastic = false;
+      deltalambda = ((tdnorm /sqrt(2.0 * pointwiseYieldvalue)) - 1.0) / alphavalue;
+      double templambda = lambdaValue[i];
+      lambdaValue[i] = templambda + deltalambda;
+    } 
+
+    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 (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 (periodic) domain->minimum_image(delx0,dely0,delz0);
+      jtype = type[j];
+      delta = cut[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);
+             
+      //Elastic Part            
+      rk = ((3.0 * bulkmodulus[itype][itype]) * ( (omega_plus * theta[i] / wvolume[i]) +
+         ( omega_minus * theta[j] / wvolume[j] ) ) ) * r0[i][jj];
+        
+      if (r > 0.0) fbond = -((rk/r) * vfrac[j] * vfrac_scale);
+      else fbond = 0.0;
+      
+      //Plastic part
+              
+      double deviatoric_extension = dr - (theta[i]* r0[i][jj] / 3.0);
+      edpNp1 = deviatorPlasticextension[i][jj];
+
+  
+      double tdtrialValue = ( 15 * shearmodulus[itype][itype]) *
+        ( (omega_plus / wvolume[i]) + (omega_minus / wvolume[j]) ) * 
+           (deviatoric_extension - edpNp1);    
+         
+      if(elastic) {
+        rkNew = tdtrialValue;
+      }
+      else {
+        rkNew = (sqrt(2.0*pointwiseYieldvalue) * tdtrialValue) / tdnorm;
+        deviatorPlasticExtTemp[i][jj] = edpNp1 + rkNew * deltalambda;
+      }  
+      
+                       
+      if (r > 0.0) fbondElastoPlastic = -((rkNew/r) * vfrac[j] * vfrac_scale);
+      else fbondElastoPlastic = 0.0;
+
+
+      // total Force state: elastic +  plastic
+      fbondFinal=fbond+fbondElastoPlastic;
+      fbond=fbondFinal;
+      
+   
+      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
+
+      if (eflag) evdwl =  (0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) *
+                    omega_plus * deviatoric_extension * deviatoric_extension) + 
+                    (0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) *
+                       omega_plus * (deviatoric_extension - edpNp1) * 
+                      (deviatoric_extension-edpNp1)) * vfrac[j] * vfrac_scale;
+      if (evflag) ev_tally(i,i,nlocal,0,0.5*evdwl,0.0,
+                           0.5*fbond*vfrac[i],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;
+    }
+  }
+
+  // store new s0
+
+  for (i = 0; i < nlocal; i++) s0[i] = s0_new[i];
+  
+  for (i = 0; i < nlocal; i++) {
+      jnum = npartner[i];
+      for (jj = 0; jj < jnum; jj++) {
+         double temp_data = deviatorPlasticExtTemp[i][jj];
+         deviatorPlasticextension[i][jj] = temp_data;
+      }
+  }    
+         
+}
+
+/* ----------------------------------------------------------------------
+   allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairPeriEPS::allocate()
+{
+  allocated = 1;
+  int n = atom->ntypes;
+
+  memory->create(setflag,n+1,n+1,"pair:setflag");
+  for (int i = 1; i <= n; i++)
+    for (int j = i; j <= n; j++)
+      setflag[i][j] = 0;
+
+  memory->create(cutsq,n+1,n+1,"pair:cutsq");
+  memory->create(bulkmodulus,n+1,n+1,"pair:bulkmodulus");
+  memory->create(shearmodulus,n+1,n+1,"pair:shearmodulus");
+  memory->create(s00,n+1,n+1,"pair:s00");
+  memory->create(alpha,n+1,n+1,"pair:alpha");
+  memory->create(cut,n+1,n+1,"pair:cut");
+  memory->create(m_yieldstress,n+1,n+1,"pair:m_yieldstress");
+}
+
+/* ----------------------------------------------------------------------
+   global settings
+------------------------------------------------------------------------- */
+
+void PairPeriEPS::settings(int narg, char **arg)
+{
+  if (narg) error->all(FLERR,"Illegal pair_style command");
+}
+
+/* ----------------------------------------------------------------------
+   set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairPeriEPS::coeff(int narg, char **arg)
+{
+  if (narg != 8) error->all(FLERR,"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]);
+  double myieldstress_one = atof(arg[7]);
+
+  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;
+      m_yieldstress[i][j] = myieldstress_one;
+      setflag[i][j] = 1;
+      count++;
+    }
+  }
+
+  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+   init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairPeriEPS::init_one(int i, int j)
+{
+  if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
+
+  bulkmodulus[j][i] = bulkmodulus[i][j];
+  shearmodulus[j][i] = shearmodulus[i][j];
+  s00[j][i] = s00[i][j];
+  alpha[j][i] = alpha[i][j];
+  cut[j][i] = cut[i][j];
+  m_yieldstress[j][i] = m_yieldstress[i][j];
+  return cut[i][j];
+}
+
+/* ----------------------------------------------------------------------
+   init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairPeriEPS::init_style()
+{
+  // error checks
+
+  if (!atom->peri_flag) 
+    error->all(FLERR,"Pair style peri requires atom style peri");
+  if (atom->map_style == 0)
+    error->all(FLERR,"Pair peri requires an atom map, see atom_modify");
+
+  if (domain->lattice == NULL)
+    error->all(FLERR,"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(FLERR,"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(FLERR,"Fix peri neigh does not exist");
+
+  neighbor->request(this);
+}
+
+/* ----------------------------------------------------------------------
+  proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairPeriEPS::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);
+        fwrite(&m_yieldstress[i][j],sizeof(double),1,fp);
+      }
+    }
+}
+
+/* ----------------------------------------------------------------------
+  proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairPeriEPS::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);
+          fread(&m_yieldstress[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);
+        MPI_Bcast(&m_yieldstress[i][j],1,MPI_DOUBLE,0,world);
+      }
+    }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based arrays
+------------------------------------------------------------------------- */
+
+double PairPeriEPS::memory_usage()
+{
+  double bytes = 2 * nmax * sizeof(double);
+  return bytes;
+}
+
+/* ----------------------------------------------------------------------
+   influence function definition
+------------------------------------------------------------------------- */
+
+double PairPeriEPS::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(FLERR,"Divide by 0 in influence function");
+  omega = 1.0/r;
+  return omega;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairPeriEPS::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;
+  tagint **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner;
+  int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner;
+  double *wvolume = ((FixPeriNeigh *) modify->fix[ifix_peri])->wvolume;
+
+  int periodic = domain->xperiodic || domain->yperiodic || domain->zperiodic;
+
+  // compute the dilatation theta
+
+  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];
+
+    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 (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 (periodic) 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 = cut[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;
+
+    }
+
+    // if wvolume[i] is zero, then particle i has no bonds
+    // therefore, the dilatation is set to
+
+    if (wvolume[i] != 0.0) theta[i] = (3.0/wvolume[i]) * theta[i];
+    else theta[i] = 0;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairPeriEPS::compute_DeviatoricForceStateNorm(int i)
+{
+  int 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 tdtrial;
+  double norm = 0.0;
+
+  double **x = atom->x;
+  int *type = atom->type;
+  double **x0 = atom->x0;
+  double *s0 = atom->s0;
+  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;
+  tagint **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner;
+  int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner;
+  double *wvolume = ((FixPeriNeigh *) modify->fix[ifix_peri])->wvolume;
+  double **deviatorPlasticextension = 
+    ((FixPeriNeigh *) modify->fix[ifix_peri])->deviatorPlasticextension;
+
+  int periodic = domain->xperiodic || domain->yperiodic || domain->zperiodic;
+  
+  // compute the dilatation theta
+
+    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];
+    itype = type[i];
+
+    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;
+      }     
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      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 (periodic) domain->minimum_image(delx0,dely0,delz0);
+      r = sqrt(rsq);
+      dr = r - r0[i][jj];
+      if (fabs(dr) < 2.2204e-016) dr = 0.0;
+      
+      // scale vfrac[j] if particle j near the horizon
+      double vfrac_scale;
+      
+      jtype = type[j];
+      double delta = cut[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;
+      
+      double ed = dr - (theta[i] * r0[i][jj])/3;
+      double edPNP1 = deviatorPlasticextension[i][jj];
+
+      jtype = type[j];
+      delta = cut[itype][jtype];
+      
+      double omega_plus  = influence_function(-1.0*delx0,-1.0*dely0,-1.0*delz0);
+      double omega_minus = influence_function(delx0,dely0,delz0);
+      
+      double stretch = dr / r0[i][jj];
+      
+      tdtrial = ( 15 * shearmodulus[itype][itype]) *
+           ((omega_plus * theta[i] / wvolume[i]) +
+             ( omega_minus * theta[j] / wvolume[j] ) ) * (ed - edPNP1);
+          
+      norm += tdtrial * tdtrial * vfrac[j] * vfrac_scale;   
+    }
+  return sqrt(norm);
+}
+
+
+/* ----------------------------------------------------------------------
+   communication routines
+---------------------------------------------------------------------- */
+
+int PairPeriEPS::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 PairPeriEPS::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_eps.h

@@ -0,0 +1,111 @@
+/* ----------------------------------------------------------------------
+   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/eps,PairPeriEPS)
+
+#else
+
+#ifndef LMP_PAIR_PERI_EPS_H
+#define LMP_PAIR_PERI_EPS_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairPeriEPS : public Pair {
+ public:
+  double *theta;
+
+  PairPeriEPS(class LAMMPS *);
+  virtual ~PairPeriEPS();
+  int pack_comm(int, int *, double *, int, int *);
+  void unpack_comm(int, int, double *);
+
+  virtual 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 memory_usage();
+  double influence_function(double, double, double);
+  void compute_dilatation();
+  double compute_DeviatoricForceStateNorm(int);
+
+ protected:
+  int ifix_peri;
+  double **bulkmodulus;
+  double **shearmodulus;
+  double **s00, **alpha;
+  double **cut, **m_yieldstress;   //NEW: **m_yieldstress
+
+  double *s0_new;
+  int nmax;
+
+  void allocate();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+E: All pair coeffs are not set
+
+All pair coefficients must be set in the data file or by the
+pair_coeff command before running a simulation.
+
+E: Pair style peri requires atom style peri
+
+Self-explanatory.
+
+E: Pair peri requires an atom map, see atom_modify
+
+Even for atomic systems, an atom map is required to find Peridynamic
+bonds.  Use the atom_modify command to define one.
+
+E: Pair peri requires a lattice be defined
+
+Use the lattice command for this purpose.
+
+E: Pair peri lattice is not identical in x, y, and z
+
+The lattice defined by the lattice command must be cubic.
+
+E: Fix peri neigh does not exist
+
+Somehow a fix that the pair style defines has been deleted.
+
+E: Divide by 0 in influence function of pair peri/lps
+
+This should not normally occur.  It is likely a problem with your
+model.
+
+*/

src/PERI/pair_peri_lps.cpp

@@ -625,7 +625,6 @@ void PairPeriLPS::compute_dilatation()
 int PairPeriLPS::pack_comm(int n, int *list, double *buf,
                            int pbc_flag, int *pbc)
 {
-
   int i,j,m;
 
   m = 0;

src/PERI/pair_peri_lps.h

@@ -26,6 +26,8 @@ namespace LAMMPS_NS {
 
 class PairPeriLPS : public Pair {
  public:
+  double *theta;
+
   PairPeriLPS(class LAMMPS *);
   virtual ~PairPeriLPS();
   int pack_comm(int, int *, double *, int, int *);
@@ -52,7 +54,6 @@ class PairPeriLPS : public Pair {
   double **cut;
 
   double *s0_new;
-  double *theta;
   int nmax;
 
   void allocate();

src/PERI/pair_peri_ves.h

@@ -26,6 +26,8 @@ namespace LAMMPS_NS {
 
 class PairPeriVES : public Pair {
  public:
+  double *theta;
+
   PairPeriVES(class LAMMPS *);
   virtual ~PairPeriVES();
   int pack_comm(int, int *, double *, int, int *);
@@ -49,10 +51,11 @@ class PairPeriVES : public Pair {
   double **bulkmodulus;
   double **shearmodulus;
   double **s00, **alpha;
-  double **cut, **m_lambdai, **m_taubi;   //NEW: **m_lambdai,**m_taubi
+  double **cut;
+  double **m_lambdai;
+  double **m_taubi;
 
   double *s0_new;
-  double *theta;
   int nmax;
 
   void allocate();