From dd84b5c97596e15f1abf95a9ed2eacc6b0843a38 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Fri, 30 Nov 2018 07:50:57 -0500
Subject: [PATCH] make compute stress/mop and stress/mop/profile compatible
with per-atom masses
---
src/USER-MISC/compute_stress_mop.cpp | 35 ++++++++-----
src/USER-MISC/compute_stress_mop_profile.cpp | 52 +++++++++++++-------
2 files changed, 59 insertions(+), 28 deletions(-)
diff --git a/src/USER-MISC/compute_stress_mop.cpp b/src/USER-MISC/compute_stress_mop.cpp
index e496fd4dfa..cb2a0b9713 100644
--- a/src/USER-MISC/compute_stress_mop.cpp
+++ b/src/USER-MISC/compute_stress_mop.cpp
@@ -70,7 +70,7 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[4],"center")==0) {
pos = 0.5*(domain->boxlo[dir]+domain->boxhi[dir]);
} else pos = force->numeric(FLERR,arg[4]);
-
+
if ( pos < (domain->boxlo[dir]+domain->prd_half[dir]) ) {
pos1 = pos + domain->prd[dir];
} else {
@@ -118,7 +118,7 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) :
if (pos >domain->boxhi[dir] || pos <domain->boxlo[dir])
error->all(FLERR, "Plane for compute stress/mop is out of bounds");
-
+
// Initialize some variables
values_local = values_global = vector = NULL;
@@ -168,13 +168,13 @@ void ComputeStressMop::init()
// Timestep Value
dt = update->dt;
-
+
// Error check
-
+
// Compute stress/mop requires fixed simulation box
if (domain->box_change_size || domain->box_change_shape || domain->deform_flag)
error->all(FLERR, "Compute stress/mop requires a fixed simulation box");
-
+
// This compute requires a pair style with pair_single method implemented
if (force->pair == NULL)
@@ -252,6 +252,7 @@ void ComputeStressMop::compute_pairs()
int *ilist,*jlist,*numneigh,**firstneigh;
double *mass = atom->mass;
+ double *rmass = atom->rmass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@@ -388,9 +389,15 @@ void ComputeStressMop::compute_pairs()
fi[2] = atom->f[i][2];
//coordinates at t-dt (based on Velocity-Verlet alg.)
- xj[0] = xi[0]-vi[0]*dt+fi[0]/2/mass[itype]*dt*dt*ftm2v;
- xj[1] = xi[1]-vi[1]*dt+fi[1]/2/mass[itype]*dt*dt*ftm2v;
- xj[2] = xi[2]-vi[2]*dt+fi[2]/2/mass[itype]*dt*dt*ftm2v;
+ if (rmass) {
+ xj[0] = xi[0]-vi[0]*dt+fi[0]/2/rmass[i]*dt*dt*ftm2v;
+ xj[1] = xi[1]-vi[1]*dt+fi[1]/2/rmass[i]*dt*dt*ftm2v;
+ xj[2] = xi[2]-vi[2]*dt+fi[2]/2/rmass[i]*dt*dt*ftm2v;
+ } else {
+ xj[0] = xi[0]-vi[0]*dt+fi[0]/2/mass[itype]*dt*dt*ftm2v;
+ xj[1] = xi[1]-vi[1]*dt+fi[1]/2/mass[itype]*dt*dt*ftm2v;
+ xj[2] = xi[2]-vi[2]*dt+fi[2]/2/mass[itype]*dt*dt*ftm2v;
+ }
// because LAMMPS does not put atoms back in the box
// at each timestep, must check atoms going through the
@@ -406,9 +413,15 @@ void ComputeStressMop::compute_pairs()
//approximate crossing velocity by v(t-dt/2) (based on Velocity-Verlet alg.)
double vcross[3];
- vcross[0] = vi[0]-fi[0]/mass[itype]/2*ftm2v*dt;
- vcross[1] = vi[1]-fi[1]/mass[itype]/2*ftm2v*dt;
- vcross[2] = vi[2]-fi[2]/mass[itype]/2*ftm2v*dt;
+ if (rmass) {
+ vcross[0] = vi[0]-fi[0]/rmass[i]/2*ftm2v*dt;
+ vcross[1] = vi[1]-fi[1]/rmass[i]/2*ftm2v*dt;
+ vcross[2] = vi[2]-fi[2]/rmass[i]/2*ftm2v*dt;
+ } else {
+ vcross[0] = vi[0]-fi[0]/mass[itype]/2*ftm2v*dt;
+ vcross[1] = vi[1]-fi[1]/mass[itype]/2*ftm2v*dt;
+ vcross[2] = vi[2]-fi[2]/mass[itype]/2*ftm2v*dt;
+ }
values_local[m] += mass[itype]*vcross[0]*sgn/dt/area*nktv2p/ftm2v;
values_local[m+1] += mass[itype]*vcross[1]*sgn/dt/area*nktv2p/ftm2v;
diff --git a/src/USER-MISC/compute_stress_mop_profile.cpp b/src/USER-MISC/compute_stress_mop_profile.cpp
index 76bf487630..5499a37a46 100644
--- a/src/USER-MISC/compute_stress_mop_profile.cpp
+++ b/src/USER-MISC/compute_stress_mop_profile.cpp
@@ -168,20 +168,20 @@ void ComputeStressMopProfile::init()
// timestep Value
dt = update->dt;
-
+
// Error check
// Compute stress/mop/profile requires fixed simulation box
-
+
if (domain->box_change_size || domain->box_change_shape || domain->deform_flag)
error->all(FLERR, "Compute stress/mop/profile requires a fixed simulation box");
-
+
//This compute requires a pair style with pair_single method implemented
if (force->pair == NULL)
error->all(FLERR,"No pair style is defined for compute stress/mop/profile");
if (force->pair->single_enable == 0)
error->all(FLERR,"Pair style does not support compute stress/mop/profile");
-
+
// Warnings
if (me==0){
@@ -232,7 +232,7 @@ void ComputeStressMopProfile::compute_array()
MPI_Allreduce(&values_local[0][0],&values_global[0][0],nbins*nvalues,
MPI_DOUBLE,MPI_SUM,world);
- int ibin,m,mo;
+ int ibin,m,mo;
for (ibin=0; ibin<nbins; ibin++) {
array[ibin][0] = coord[ibin][0];
mo=1;
@@ -257,9 +257,10 @@ void ComputeStressMopProfile::compute_pairs()
double delx,dely,delz;
double rsq,fpair,factor_coul,factor_lj;
int *ilist,*jlist,*numneigh,**firstneigh;
- double pos,pos1;
+ double pos,pos1;
double *mass = atom->mass;
+ double *rmass = atom->rmass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@@ -411,9 +412,15 @@ void ComputeStressMopProfile::compute_pairs()
fi[2] = atom->f[i][2];
//coordinates at t-dt (based on Velocity-Verlet alg.)
- xj[0] = xi[0]-vi[0]*dt+fi[0]/2/mass[itype]*dt*dt*ftm2v;
- xj[1] = xi[1]-vi[1]*dt+fi[1]/2/mass[itype]*dt*dt*ftm2v;
- xj[2] = xi[2]-vi[2]*dt+fi[2]/2/mass[itype]*dt*dt*ftm2v;
+ if (rmass) {
+ xj[0] = xi[0]-vi[0]*dt+fi[0]/2/rmass[i]*dt*dt*ftm2v;
+ xj[1] = xi[1]-vi[1]*dt+fi[1]/2/rmass[i]*dt*dt*ftm2v;
+ xj[2] = xi[2]-vi[2]*dt+fi[2]/2/rmass[i]*dt*dt*ftm2v;
+ } else {
+ xj[0] = xi[0]-vi[0]*dt+fi[0]/2/mass[itype]*dt*dt*ftm2v;
+ xj[1] = xi[1]-vi[1]*dt+fi[1]/2/mass[itype]*dt*dt*ftm2v;
+ xj[2] = xi[2]-vi[2]*dt+fi[2]/2/mass[itype]*dt*dt*ftm2v;
+ }
for (ibin=0;ibin<nbins;ibin++) {
pos = coord[ibin][0];
@@ -424,14 +431,25 @@ void ComputeStressMopProfile::compute_pairs()
sgn = copysign(1.0,vi[dir]);
//approximate crossing velocity by v(t-dt/2) (based on Velocity-Verlet alg.)
- double vcross[3];
- vcross[0] = vi[0]-fi[0]/mass[itype]/2*ftm2v*dt;
- vcross[1] = vi[1]-fi[1]/mass[itype]/2*ftm2v*dt;
- vcross[2] = vi[2]-fi[2]/mass[itype]/2*ftm2v*dt;
-
- values_local[ibin][m] += mass[itype]*vcross[0]*sgn/dt/area*nktv2p/ftm2v;
- values_local[ibin][m+1] += mass[itype]*vcross[1]*sgn/dt/area*nktv2p/ftm2v;
- values_local[ibin][m+2] += mass[itype]*vcross[2]*sgn/dt/area*nktv2p/ftm2v;
+ if (rmass) {
+ double vcross[3];
+ vcross[0] = vi[0]-fi[0]/rmass[i]/2*ftm2v*dt;
+ vcross[1] = vi[1]-fi[1]/rmass[i]/2*ftm2v*dt;
+ vcross[2] = vi[2]-fi[2]/rmass[i]/2*ftm2v*dt;
+
+ values_local[ibin][m] += rmass[i]*vcross[0]*sgn/dt/area*nktv2p/ftm2v;
+ values_local[ibin][m+1] += rmass[i]*vcross[1]*sgn/dt/area*nktv2p/ftm2v;
+ values_local[ibin][m+2] += rmass[i]*vcross[2]*sgn/dt/area*nktv2p/ftm2v;
+ } else {
+ double vcross[3];
+ vcross[0] = vi[0]-fi[0]/mass[itype]/2*ftm2v*dt;
+ vcross[1] = vi[1]-fi[1]/mass[itype]/2*ftm2v*dt;
+ vcross[2] = vi[2]-fi[2]/mass[itype]/2*ftm2v*dt;
+
+ values_local[ibin][m] += mass[itype]*vcross[0]*sgn/dt/area*nktv2p/ftm2v;
+ values_local[ibin][m+1] += mass[itype]*vcross[1]*sgn/dt/area*nktv2p/ftm2v;
+ values_local[ibin][m+2] += mass[itype]*vcross[2]*sgn/dt/area*nktv2p/ftm2v;
+ }
}
}
}
--
GitLab