diff --git a/src/PERI/pair_peri_eps.cpp b/src/PERI/pair_peri_eps.cpp
index b5807c0e3c19c2c148e6b8f43033bb36a145722d..670e1d6937aaafce5962ca6f7a24bc2fb945f6fc 100644
--- a/src/PERI/pair_peri_eps.cpp
+++ b/src/PERI/pair_peri_eps.cpp
@@ -46,7 +46,7 @@ PairPeriEPS::PairPeriEPS(LAMMPS *lmp) : Pair(lmp)
ifix_peri = -1;
- nmax = 0;
+ nmax = -1;
s0_new = NULL;
theta = NULL;
@@ -209,7 +209,7 @@ void PairPeriEPS::compute(int eflag, int vflag)
for (i = 0; i < nlocal; i++) maxpartner = MAX(maxpartner,npartner[i]);
- if (atom->nmax > nmax) {
+ if (nlocal > nmax) {
memory->destroy(s0_new);
memory->destroy(theta);
nmax = atom->nmax;
@@ -220,13 +220,13 @@ void PairPeriEPS::compute(int eflag, int vflag)
// ******** temp array to store Plastic extension *********** ///
// create on heap to reduce stack use and to allow for faster zeroing
- double **deviatorPlasticExtTemp;
- memory->create(deviatorPlasticExtTemp,nlocal,maxpartner,"pair:plastext");
- memset(&(deviatorPlasticExtTemp[0][0]),0,sizeof(double)*nlocal*maxpartner);
+ double **deviatorPlasticExtTemp = NULL;
+ if (nlocal*maxpartner > 0) {
+ memory->create(deviatorPlasticExtTemp,nlocal,maxpartner,"pair:plastext");
+ memset(&(deviatorPlasticExtTemp[0][0]),0,sizeof(double)*nlocal*maxpartner);
+ }
// ******** temp array to store Plastic extension *********** ///
-
-
// compute the dilatation on each particle
compute_dilatation();
@@ -280,12 +280,10 @@ void PairPeriEPS::compute(int eflag, int vflag)
double fsurf = (tdnorm * tdnorm)/2 - pointwiseYieldvalue;
bool elastic = true;
- double alphavalue = (15 * shearmodulus[itype][itype]) /wvolume[i];
-
-
- if (fsurf>0) {
+ if (fsurf > 0) {
elastic = false;
- deltalambda = ((tdnorm /sqrt(2.0 * pointwiseYieldvalue)) - 1.0) / alphavalue;
+ deltalambda = ((tdnorm /sqrt(2.0 * pointwiseYieldvalue)) - 1.0) * wvolume[i]
+ / (15 * shearmodulus[itype][itype]);
double templambda = lambdaValue[i];
lambdaValue[i] = templambda + deltalambda;
}
@@ -348,10 +346,9 @@ void PairPeriEPS::compute(int eflag, int vflag)
( (omega_plus / wvolume[i]) + (omega_minus / wvolume[j]) ) *
(deviatoric_extension - edpNp1);
- if(elastic) {
+ if (elastic) {
rkNew = tdtrialValue;
- }
- else {
+ } else {
rkNew = (sqrt(2.0*pointwiseYieldvalue) * tdtrialValue) / tdnorm;
deviatorPlasticExtTemp[i][jj] = edpNp1 + rkNew * deltalambda;
}
@@ -402,10 +399,12 @@ void PairPeriEPS::compute(int eflag, int vflag)
memcpy(s0,s0_new,sizeof(double)*nlocal);
- memcpy(&(deviatorPlasticextension[0][0]),
- &(deviatorPlasticExtTemp[0][0]),
- sizeof(double)*nlocal*maxpartner);
- memory->destroy(deviatorPlasticExtTemp);
+ if (nlocal*maxpartner > 0) {
+ memcpy(&(deviatorPlasticextension[0][0]),
+ &(deviatorPlasticExtTemp[0][0]),
+ sizeof(double)*nlocal*maxpartner);
+ memory->destroy(deviatorPlasticExtTemp);
+ }
}
/* ----------------------------------------------------------------------
diff --git a/src/PERI/pair_peri_lps.cpp b/src/PERI/pair_peri_lps.cpp
index cd88b41825d3b15125116f9ed0d85640c5c76e86..0fe8f29f387fee4513cb9dd1ad60305561b7f012 100644
--- a/src/PERI/pair_peri_lps.cpp
+++ b/src/PERI/pair_peri_lps.cpp
@@ -33,6 +33,7 @@
#include "memory.h"
#include "error.h"
#include "update.h"
+#include "math_const.h"
using namespace LAMMPS_NS;
@@ -174,7 +175,7 @@ void PairPeriLPS::compute(int eflag, int vflag)
// 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]);
+ (MathConst::MY_PI * cutsq[itype][jtype] * cutsq[itype][jtype]);
rk = (kshort * vfrac[j]) * (dr / cut[itype][jtype]);
if (r > 0.0) fpair = -(rk/r);
@@ -285,13 +286,14 @@ void PairPeriLPS::compute(int eflag, int vflag)
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 ((wvolume[i] > 0.0) && (wvolume[j] > 0.0)) {
+ 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;
+ } else rk = 0.0;
if (r > 0.0) fbond = -(rk/r);
else fbond = 0.0;
@@ -305,9 +307,11 @@ void PairPeriLPS::compute(int eflag, int vflag)
double deviatoric_extension = dr - (theta[i]* r0[i][jj] / 3.0);
- if (eflag) evdwl = 0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) *
+ if (eflag && (wvolume[i] > 0.0))
+ evdwl = 0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) *
omega_plus*(deviatoric_extension * deviatoric_extension) *
vfrac[j] * vfrac_scale;
+ else evdwl = 0.0;
if (evflag) ev_tally(i,i,nlocal,0,0.5*evdwl,0.0,
0.5*fbond*vfrac[i],delx,dely,delz);
diff --git a/src/USER-OMP/pair_peri_lps_omp.cpp b/src/USER-OMP/pair_peri_lps_omp.cpp
index 4876e6b15f0f67c7fde3228c06745b626f5ac6c5..a471b47750ee7416ee494bee23fefba69554db20 100644
--- a/src/USER-OMP/pair_peri_lps_omp.cpp
+++ b/src/USER-OMP/pair_peri_lps_omp.cpp
@@ -310,12 +310,14 @@ void PairPeriLPSOMP::eval(int iifrom, int iito, ThrData * const thr)
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 ((wvolume[i] > 0.0) && (wvolume[j] > 0.0)) {
+ 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;
+ } else rk = 0.0;
if (r > 0.0) fbond = -(rk/r);
else fbond = 0.0;
@@ -327,9 +329,11 @@ void PairPeriLPSOMP::eval(int iifrom, int iito, ThrData * const thr)
// 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]) *
+ if (EFLAG && (wvolume[i] > 0.0))
+ evdwl = 0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) *
omega_plus*(deviatoric_extension * deviatoric_extension) *
vfrac[j] * vfrac_scale;
+ else evdwl = 0.0;
if (EVFLAG) ev_tally_thr(this,i,i,nlocal,0,0.5*evdwl,0.0,
0.5*fbond*vfrac[i],delx,dely,delz,thr);