ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #3792 from evoyiatzis/master

Browse Source 
Inclusion of bond & angle contributions to "compute stress/mop"
This commit is contained in:
Axel Kohlmeyer
2023-06-07 06:09:58 -04:00
committed by GitHub
parent 74a9cfc110 53b1af7720
commit f406aa1492
6 changed files with 702 additions and 180 deletions
Download Patch File Download Diff File Expand all files Collapse all files

118
doc/src/compute_stress_mop.rst
View File

@ -18,7 +18,7 @@ Syntax
* style = *stress/mop* or *stress/mop/profile*
* dir = *x* or *y* or *z* is the direction normal to the plane
* args = argument specific to the compute style
* keywords = *kin* or *conf* or *total* (one of more can be specified)
* keywords = *kin* or *conf* or *total* or *pair* or *bond* or *angle* (one or more can be specified)
.. parsed-literal::
@ -45,85 +45,107 @@ Examples
Description
"""""""""""
Compute *stress/mop* and compute *stress/mop/profile* define computations that
calculate components of the local stress tensor using the method of
planes :ref:`(Todd) <mop-todd>`. Specifically in compute *stress/mop* calculates 3
components are computed in directions *dir*,\ *x*\ ; *dir*,\ *y*\ ; and
*dir*,\ *z*\ ; where *dir* is the direction normal to the plane, while
in compute *stress/mop/profile* the profile of the stress is computed.
Compute *stress/mop* and compute *stress/mop/profile* define
computations that calculate components of the local stress tensor using
the method of planes :ref:`(Todd) <mop-todd>`. Specifically in compute
*stress/mop* calculates 3 components are computed in directions *dir*,\
*x*\ ; *dir*,\ *y*\ ; and *dir*,\ *z*\ ; where *dir* is the direction
normal to the plane, while in compute *stress/mop/profile* the profile
of the stress is computed.
Contrary to methods based on histograms of atomic stress (i.e., using
:doc:`compute stress/atom <compute_stress_atom>`), the method of planes is
compatible with mechanical balance in heterogeneous systems and at
:doc:`compute stress/atom <compute_stress_atom>`), the method of planes
is compatible with mechanical balance in heterogeneous systems and at
interfaces :ref:`(Todd) <mop-todd>`.
The stress tensor is the sum of a kinetic term and a configurational
term, which are given respectively by Eq. (21) and Eq. (16) in
:ref:`(Todd) <mop-todd>`. For the kinetic part, the algorithm considers that
atoms have crossed the plane if their positions at times :math:`t-\Delta t`
and :math:`t` are one on either side of the plane, and uses the velocity at
time :math:`t-\Delta t/2` given by the velocity Verlet algorithm.
:ref:`(Todd) <mop-todd>`. For the kinetic part, the algorithm considers
that atoms have crossed the plane if their positions at times
:math:`t-\Delta t` and :math:`t` are one on either side of the plane,
and uses the velocity at time :math:`t-\Delta t/2` given by the velocity
Verlet algorithm.
Between one and three keywords can be used to indicate which
contributions to the stress must be computed: kinetic stress (kin),
configurational stress (conf), and/or total stress (total).
.. versionadded:: TBD
NOTE 1: The configurational stress is computed considering all pairs of atoms where at least one atom belongs to group group-ID.
contributions from bond and angle potentials
Between one and six keywords can be used to indicate which contributions
to the stress must be computed: total stress (total), kinetic stress
(kin), configurational stress (conf), stress due to bond stretching
(bond), stress due to angle bending (angle) and/or due to pairwise
non-bonded interactions (pair). The angle keyword is currently
available only for the *stress/mop* command and **not** the
*stress/mop/profile* command.
NOTE 1: The configurational stress is computed considering all pairs of
atoms where at least one atom belongs to group group-ID.
NOTE 2: The local stress does not include any Lennard-Jones tail
corrections to the stress added by the :doc:`pair_modify tail yes <pair_modify>`
command, since those are contributions to the global system pressure.
corrections to the stress added by the :doc:`pair_modify tail yes
<pair_modify>` command, since those are contributions to the global
system pressure.
NOTE 3: The local stress profile generated by compute *stress/mop/profile*
is similar to that obtained by compute
:doc:`stress/cartesian <compute_stress_profile>`.
A key difference is that compute *stress/mop/profile* considers particles
crossing a set of planes, while compute *stress/cartesian* computes averages
for a set of small volumes. More information
on the similarities and differences can be found in
:ref:`(Ikeshoji)<Ikeshoji2>`.
NOTE 3: The local stress profile generated by compute
*stress/mop/profile* is similar to that obtained by compute
:doc:`stress/cartesian <compute_stress_profile>`. A key difference is
that compute *stress/mop/profile* considers particles crossing a set of
planes, while compute *stress/cartesian* computes averages for a set of
small volumes. More information on the similarities and differences can
be found in :ref:`(Ikeshoji)<Ikeshoji2>`.
Output info
"""""""""""
Compute *stress/mop* calculates a global vector (indices starting at 1), with 3
values for each declared keyword (in the order the keywords have been
declared). For each keyword, the stress tensor components are ordered as
follows: stress_dir,x, stress_dir,y, and stress_dir,z.
Compute *stress/mop* calculates a global vector (indices starting at 1),
with 3 values for each declared keyword (in the order the keywords have
been declared). For each keyword, the stress tensor components are
ordered as follows: stress_dir,x, stress_dir,y, and stress_dir,z.
Compute *stress/mop/profile* instead calculates a global array, with 1 column
giving the position of the planes where the stress tensor was computed,
and with 3 columns of values for each declared keyword (in the order the
keywords have been declared). For each keyword, the profiles of stress
tensor components are ordered as follows: stress_dir,x; stress_dir,y;
and stress_dir,z.
Compute *stress/mop/profile* instead calculates a global array, with 1
column giving the position of the planes where the stress tensor was
computed, and with 3 columns of values for each declared keyword (in the
order the keywords have been declared). For each keyword, the profiles
of stress tensor components are ordered as follows: stress_dir,x;
stress_dir,y; and stress_dir,z.
The values are in pressure :doc:`units <units>`.
The values produced by this compute can be accessed by various :doc:`output commands <Howto_output>`.
For instance, the results can be written to a file using the
:doc:`fix ave/time <fix_ave_time>` command. Please see the example
in the examples/PACKAGES/mop folder.
The values produced by this compute can be accessed by various
:doc:`output commands <Howto_output>`. For instance, the results can be
written to a file using the :doc:`fix ave/time <fix_ave_time>`
command. Please see the example in the examples/PACKAGES/mop folder.
Restrictions
""""""""""""
These styles are part of the EXTRA-COMPUTE package. They are only enabled if
LAMMPS is built with that package. See the :doc:`Build package <Build_package>`
doc page on for more info.
These styles are part of the EXTRA-COMPUTE package. They are only
enabled if LAMMPS is built with that package. See the :doc:`Build
package <Build_package>` doc page on for more info.
The method is only implemented for 3d orthogonal simulation boxes whose
size does not change in time, and axis-aligned planes.
The method only works with two-body pair interactions, because it
requires the class method pair->single() to be implemented. In
particular, it does not work with more than two-body pair interactions,
intra-molecular interactions, and long range (kspace) interactions.
requires the class method ``Pair::single()`` to be implemented, which is
not possible for manybody potentials. In particular, compute
*stress/mop/profile* does not work with more than two-body pair
interactions, long range (kspace) interactions and
angle/dihedral/improper intramolecular interactions. Similarly, compute
*stress/mop* does not work with more than two-body pair interactions,
long range (kspace) interactions and dihedral/improper intramolecular
interactions but works with all bond interactions with the class method
single() implemented and all angle interactions with the class method
born_matrix() implemented.
Related commands
""""""""""""""""
:doc:`compute stress/atom <compute_stress_atom>`, :doc:`compute pressure <compute_pressure>`, :doc:`compute stress/cartesian <compute_stress_profile>`, :doc:`compute stress/cylinder <compute_stress_profile>`, :doc:`compute stress/spherical <compute_stress_profile>`
:doc:`compute stress/atom <compute_stress_atom>`,
:doc:`compute pressure <compute_pressure>`,
:doc:`compute stress/cartesian <compute_stress_profile>`,
:doc:`compute stress/cylinder <compute_stress_profile>`,
:doc:`compute stress/spherical <compute_stress_profile>`
Default
"""""""

492
src/EXTRA-COMPUTE/compute_stress_mop.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -14,15 +13,21 @@
/*------------------------------------------------------------------------
Contributing Authors : Romain Vermorel (LFCR), Laurent Joly (ULyon)
Support for bonds and angles added by : Evangelos Voyiatzis (NovaMechanics)
--------------------------------------------------------------------------*/
#include "compute_stress_mop.h"
#include "angle.h"
#include "atom.h"
#include "atom_vec.h"
#include "bond.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "molecule.h"
#include "neigh_list.h"
#include "neighbor.h"
#include "pair.h"
@ -33,40 +38,52 @@
using namespace LAMMPS_NS;
enum{X,Y,Z};
enum{TOTAL,CONF,KIN};
#define BIG 1000000000
enum { X, Y, Z };
enum { TOTAL, CONF, KIN, PAIR, BOND, ANGLE };
// clang-format off
/* ---------------------------------------------------------------------- */
ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg < 6) error->all(FLERR,"Illegal compute stress/mop command");
if (narg < 6) utils::missing_cmd_args(FLERR, "compute stress/mop", error);
MPI_Comm_rank(world,&me);
bondflag = 0;
angleflag = 0;
// set compute mode and direction of plane(s) for pressure calculation
if (strcmp(arg[3],"x")==0) {
if (strcmp(arg[3],"x") == 0) {
dir = X;
} else if (strcmp(arg[3],"y")==0) {
} else if (strcmp(arg[3],"y") == 0) {
dir = Y;
} else if (strcmp(arg[3],"z")==0) {
} else if (strcmp(arg[3],"z") == 0) {
dir = Z;
} else error->all(FLERR,"Illegal compute stress/mop command");
// Position of the plane
if (strcmp(arg[4],"lower")==0) {
if (strcmp(arg[4],"lower") == 0) {
pos = domain->boxlo[dir];
} else if (strcmp(arg[4],"upper")==0) {
} else if (strcmp(arg[4],"upper") == 0) {
pos = domain->boxhi[dir];
} else if (strcmp(arg[4],"center")==0) {
} else if (strcmp(arg[4],"center") == 0) {
pos = 0.5*(domain->boxlo[dir]+domain->boxhi[dir]);
} else pos = utils::numeric(FLERR,arg[4],false,lmp);
// plane inside the box
if ((pos > domain->boxhi[dir]) || (pos < domain->boxlo[dir])) {
error->warning(FLERR, "The specified initial plane lies outside of the simulation box");
double dx[3] = {0.0, 0.0, 0.0};
dx[dir] = pos - 0.5*(domain->boxhi[dir] + domain->boxlo[dir]);
domain->minimum_image(dx[0], dx[1], dx[2]);
pos = 0.5*(domain->boxhi[dir] + domain->boxlo[dir]) + dx[dir];
if ((pos > domain->boxhi[dir]) || (pos < domain->boxlo[dir]))
error->all(FLERR, "Plane for compute stress/mop is out of bounds");
}
if (pos < (domain->boxlo[dir]+domain->prd_half[dir])) {
pos1 = pos + domain->prd[dir];
} else {
@ -96,34 +113,53 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) :
which[nvalues] = TOTAL;
nvalues++;
}
} else if (strcmp(arg[iarg],"pair") == 0) {
for (i=0; i<3; i++) {
which[nvalues] = PAIR;
nvalues++;
}
} else if (strcmp(arg[iarg],"bond") == 0) {
for (i=0; i<3; i++) {
which[nvalues] = BOND;
nvalues++;
}
} else if (strcmp(arg[iarg],"angle") == 0) {
for (i=0; i<3; i++) {
which[nvalues] = ANGLE;
nvalues++;
}
} else error->all(FLERR, "Illegal compute stress/mop command"); //break;
iarg++;
}
// Error check
// Error checks
// 3D only
if (domain->dimension < 3)
error->all(FLERR, "Compute stress/mop incompatible with simulation dimension");
if (domain->dimension != 3)
error->all(FLERR, "Compute stress/mop requires a 3d system");
// orthogonal simulation box
if (domain->triclinic != 0)
error->all(FLERR, "Compute stress/mop incompatible with triclinic simulation box");
// plane inside the box
if (pos >domain->boxhi[dir] || pos <domain->boxlo[dir])
error->all(FLERR, "Plane for compute stress/mop is out of bounds");
error->all(FLERR, "Compute stress/mop is incompatible with triclinic simulation box");
// Initialize some variables
values_local = values_global = vector = nullptr;
bond_local = nullptr;
bond_global = nullptr;
angle_local = nullptr;
angle_global = nullptr;
// this fix produces a global vector
memory->create(vector,nvalues,"stress/mop:vector");
memory->create(values_local,nvalues,"stress/mop/spatial:values_local");
memory->create(values_global,nvalues,"stress/mop/spatial:values_global");
memory->create(values_local,nvalues,"stress/mop:values_local");
memory->create(values_global,nvalues,"stress/mop:values_global");
memory->create(bond_local,nvalues,"stress/mop:bond_local");
memory->create(bond_global,nvalues,"stress/mop:bond_global");
memory->create(angle_local,nvalues,"stress/mop:angle_local");
memory->create(angle_global,nvalues,"stress/mop:angle_global");
size_vector = nvalues;
vector_flag = 1;
@ -135,11 +171,14 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) :
ComputeStressMop::~ComputeStressMop()
{
delete [] which;
delete[] which;
memory->destroy(values_local);
memory->destroy(values_global);
memory->destroy(bond_local);
memory->destroy(bond_global);
memory->destroy(angle_local);
memory->destroy(angle_global);
memory->destroy(vector);
}
@ -169,31 +208,39 @@ void ComputeStressMop::init()
// 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");
error->all(FLERR, "Compute stress/mop requires a fixed size simulation box");
// This compute requires a pair style with pair_single method implemented
if (force->pair == nullptr)
if (!force->pair)
error->all(FLERR,"No pair style is defined for compute stress/mop");
if (force->pair->single_enable == 0)
error->all(FLERR,"Pair style does not support compute stress/mop");
// Warnings
// Errors
if (me==0) {
if (comm->me == 0) {
//Compute stress/mop only accounts for pair interactions.
// issue a warning if any intramolecular potential or Kspace is defined.
// issue an error for unimplemented intramolecular potentials or Kspace.
if (force->bond!=nullptr)
error->warning(FLERR,"compute stress/mop does not account for bond potentials");
if (force->angle!=nullptr)
error->warning(FLERR,"compute stress/mop does not account for angle potentials");
if (force->dihedral!=nullptr)
error->warning(FLERR,"compute stress/mop does not account for dihedral potentials");
if (force->improper!=nullptr)
error->warning(FLERR,"compute stress/mop does not account for improper potentials");
if (force->kspace!=nullptr)
if (force->bond) bondflag = 1;
if (force->angle) {
if (force->angle->born_matrix_enable == 0) {
if ((strcmp(force->angle_style, "zero") != 0) && (strcmp(force->angle_style, "none") != 0))
error->all(FLERR,"compute stress/mop does not account for angle potentials");
} else {
angleflag = 1;
}
}
if (force->dihedral) {
if ((strcmp(force->dihedral_style, "zero") != 0) && (strcmp(force->dihedral_style, "none") != 0))
error->all(FLERR,"compute stress/mop does not account for dihedral potentials");
}
if (force->improper) {
if ((strcmp(force->improper_style, "zero") != 0) && (strcmp(force->improper_style, "none") != 0))
error->all(FLERR,"compute stress/mop does not account for improper potentials");
}
if (force->kspace)
error->warning(FLERR,"compute stress/mop does not account for kspace contributions");
}
@ -221,14 +268,31 @@ void ComputeStressMop::compute_vector()
compute_pairs();
// sum pressure contributions over all procs
MPI_Allreduce(values_local,values_global,nvalues,
MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(values_local,values_global,nvalues,MPI_DOUBLE,MPI_SUM,world);
int m;
for (m=0; m<nvalues; m++) {
vector[m] = values_global[m];
if (bondflag) {
//Compute bond contribution on separate procs
compute_bonds();
} else {
for (int i=0; i<nvalues; i++) bond_local[i] = 0.0;
}
// sum bond contribution over all procs
MPI_Allreduce(bond_local,bond_global,nvalues,MPI_DOUBLE,MPI_SUM,world);
if (angleflag) {
//Compute angle contribution on separate procs
compute_angles();
} else {
for (int i=0; i<nvalues; i++) angle_local[i] = 0.0;
}
// sum angle contribution over all procs
MPI_Allreduce(angle_local,angle_global,nvalues,MPI_DOUBLE,MPI_SUM,world);
for (int m=0; m<nvalues; m++) {
vector[m] = values_global[m] + bond_global[m] + angle_global[m];
}
}
@ -280,8 +344,8 @@ void ComputeStressMop::compute_pairs()
double xj[3];
m = 0;
while (m<nvalues) {
if (which[m] == CONF || which[m] == TOTAL) {
while (m < nvalues) {
if ((which[m] == CONF) || (which[m] == TOTAL) || (which[m] == PAIR)) {
//Compute configurational contribution to pressure
for (ii = 0; ii < inum; ii++) {
@ -316,47 +380,34 @@ void ComputeStressMop::compute_pairs()
if (newton_pair || j < nlocal) {
//check if ij pair is across plane, add contribution to pressure
if (((xi[dir]>pos) && (xj[dir]<pos)) || ((xi[dir]>pos1) && (xj[dir]<pos1))) {
if (((xi[dir] > pos) && (xj[dir] < pos)) || ((xi[dir] > pos1) && (xj[dir] < pos1))) {
pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);
values_local[m] += fpair*(xi[0]-xj[0])/area*nktv2p;
values_local[m+1] += fpair*(xi[1]-xj[1])/area*nktv2p;
values_local[m+2] += fpair*(xi[2]-xj[2])/area*nktv2p;
}
else if (((xi[dir]<pos) && (xj[dir]>pos)) || ((xi[dir]<pos1) && (xj[dir]>pos1))) {
} else if (((xi[dir] < pos) && (xj[dir] > pos)) || ((xi[dir] < pos1) && (xj[dir] > pos1))) {
pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);
values_local[m] -= fpair*(xi[0]-xj[0])/area*nktv2p;
values_local[m+1] -= fpair*(xi[1]-xj[1])/area*nktv2p;
values_local[m+2] -= fpair*(xi[2]-xj[2])/area*nktv2p;
}
} else {
if (((xi[dir]>pos) && (xj[dir]<pos)) || ((xi[dir]>pos1) && (xj[dir]<pos1))) {
if (((xi[dir] > pos) && (xj[dir] < pos)) || ((xi[dir] > pos1) && (xj[dir] < pos1))) {
pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);
values_local[m] += fpair*(xi[0]-xj[0])/area*nktv2p;
values_local[m+1] += fpair*(xi[1]-xj[1])/area*nktv2p;
values_local[m+2] += fpair*(xi[2]-xj[2])/area*nktv2p;
}
}
}
}
}
// Compute kinetic contribution to pressure
// counts local particles transfers across the plane
if (which[m] == KIN || which[m] == TOTAL) {
if ((which[m] == KIN) || (which[m] == TOTAL)) {
double sgn;
for (int i = 0; i < nlocal; i++) {
@ -383,13 +434,13 @@ void ComputeStressMop::compute_pairs()
//coordinates at t-dt (based on Velocity-Verlet alg.)
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;
xj[0] = xi[0]-vi[0]*dt+fi[0]/2.0/rmass[i]*dt*dt*ftm2v;
xj[1] = xi[1]-vi[1]*dt+fi[1]/2.0/rmass[i]*dt*dt*ftm2v;
xj[2] = xi[2]-vi[2]*dt+fi[2]/2.0/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;
xj[0] = xi[0]-vi[0]*dt+fi[0]/2.0/mass[itype]*dt*dt*ftm2v;
xj[1] = xi[1]-vi[1]*dt+fi[1]/2.0/mass[itype]*dt*dt*ftm2v;
xj[2] = xi[2]-vi[2]*dt+fi[2]/2.0/mass[itype]*dt*dt*ftm2v;
}
// because LAMMPS does not put atoms back in the box
@ -399,21 +450,21 @@ void ComputeStressMop::compute_pairs()
double pos_temp = pos+copysign(1.0,domain->prd_half[dir]-pos)*domain->prd[dir];
if (fabs(xi[dir]-pos)<fabs(xi[dir]-pos_temp)) pos_temp = pos;
if (((xi[dir]-pos_temp)*(xj[dir]-pos_temp)<0)) {
if (((xi[dir]-pos_temp)*(xj[dir]-pos_temp)) < 0) {
//sgn = copysign(1.0,vi[dir]-vcm[dir]);
// sgn = copysign(1.0,vi[dir]-vcm[dir]);
sgn = copysign(1.0,vi[dir]);
//approximate crossing velocity by v(t-dt/2) (based on Velocity-Verlet alg.)
// approximate crossing velocity by v(t-dt/2) (based on Velocity-Verlet alg.)
double vcross[3];
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;
vcross[0] = vi[0]-fi[0]/rmass[i]/2.0*ftm2v*dt;
vcross[1] = vi[1]-fi[1]/rmass[i]/2.0*ftm2v*dt;
vcross[2] = vi[2]-fi[2]/rmass[i]/2.0*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;
vcross[0] = vi[0]-fi[0]/mass[itype]/2.0*ftm2v*dt;
vcross[1] = vi[1]-fi[1]/mass[itype]/2.0*ftm2v*dt;
vcross[2] = vi[2]-fi[2]/mass[itype]/2.0*ftm2v*dt;
}
values_local[m] += mass[itype]*vcross[0]*sgn/dt/area*nktv2p/ftm2v;
@ -423,7 +474,288 @@ void ComputeStressMop::compute_pairs()
}
}
}
m+=3;
m += 3;
}
}
/*------------------------------------------------------------------------
* compute bond contribution to pressure of local proc
*-------------------------------------------------------------------------*/
void ComputeStressMop::compute_bonds()
{
int i, nb, atom1, atom2, imol, iatom, btype;
tagint tagprev;
double rsq, fpair;
double **x = atom->x;
tagint *tag = atom->tag;
int *num_bond = atom->num_bond;
tagint **bond_atom = atom->bond_atom;
int **bond_type = atom->bond_type;
int *mask = atom->mask;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
int molecular = atom->molecular;
Bond *bond = force->bond;
double dx[3] = {0.0, 0.0, 0.0};
double x_bond_1[3] = {0.0, 0.0, 0.0};
double x_bond_2[3] = {0.0, 0.0, 0.0};
double local_contribution[3] = {0.0, 0.0, 0.0};
// initialization
for (int i = 0; i < nvalues; i++) bond_local[i] = 0.0;
// loop over all bonded atoms in the current proc
for (atom1 = 0; atom1 < nlocal; atom1++) {
if (!(mask[atom1] & groupbit)) continue;
if (molecular == 1)
nb = num_bond[atom1];
else {
if (molindex[atom1] < 0) continue;
imol = molindex[atom1];
iatom = molatom[atom1];
nb = onemols[imol]->num_bond[iatom];
}
for (i = 0; i < nb; i++) {
if (molecular == 1) {
btype = bond_type[atom1][i];
atom2 = atom->map(bond_atom[atom1][i]);
} else {
tagprev = tag[atom1] - iatom - 1;
btype = onemols[imol]->bond_type[iatom][i];
atom2 = atom->map(onemols[imol]->bond_atom[iatom][i] + tagprev);
}
if (atom2 < 0 || !(mask[atom2] & groupbit)) continue;
if (newton_bond == 0 && tag[atom1] > tag[atom2]) continue;
if (btype <= 0) continue;
// minimum image of atom1 with respect to the plane of interest
dx[0] = x[atom1][0];
dx[1] = x[atom1][1];
dx[2] = x[atom1][2];
dx[dir] -= pos;
domain->minimum_image(dx[0], dx[1], dx[2]);
x_bond_1[0] = dx[0];
x_bond_1[1] = dx[1];
x_bond_1[2] = dx[2];
x_bond_1[dir] += pos;
// minimum image of atom2 with respect to atom1
dx[0] = x[atom2][0] - x_bond_1[0];
dx[1] = x[atom2][1] - x_bond_1[1];
dx[2] = x[atom2][2] - x_bond_1[2];
domain->minimum_image(dx[0], dx[1], dx[2]);
x_bond_2[0] = x_bond_1[0] + dx[0];
x_bond_2[1] = x_bond_1[1] + dx[1];
x_bond_2[2] = x_bond_1[2] + dx[2];
// check if the bond vector crosses the plane of interest
double tau = (x_bond_1[dir] - pos) / (x_bond_1[dir] - x_bond_2[dir]);
if ((tau <= 1) && (tau >= 0)) {
dx[0] = x_bond_1[0] - x_bond_2[0];
dx[1] = x_bond_1[1] - x_bond_2[1];
dx[2] = x_bond_1[2] - x_bond_2[2];
rsq = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
bond->single(btype, rsq, atom1, atom2, fpair);
double sgn = copysign(1.0, x_bond_1[dir] - pos);
local_contribution[0] += sgn*fpair*dx[0]/area*nktv2p;
local_contribution[1] += sgn*fpair*dx[1]/area*nktv2p;
local_contribution[2] += sgn*fpair*dx[2]/area*nktv2p;
}
}
}
// loop over the keywords and if necessary add the bond contribution
int m = 0;
while (m<nvalues) {
if (which[m] == CONF || which[m] == TOTAL || which[m] == BOND) {
bond_local[m] = local_contribution[0];
bond_local[m+1] = local_contribution[1];
bond_local[m+2] = local_contribution[2];
}
m += 3;
}
}
/*------------------------------------------------------------------------
compute angle contribution to pressure of local proc
-------------------------------------------------------------------------*/
void ComputeStressMop::compute_angles()
{
int na, atom1, atom2, atom3, imol, iatom, atype;
tagint tagprev;
double r1, r2, cos_theta;
double **x = atom->x;
tagint *tag = atom->tag;
int *num_angle = atom->num_angle;
tagint **angle_atom1 = atom->angle_atom1;
tagint **angle_atom2 = atom->angle_atom2;
tagint **angle_atom3 = atom->angle_atom3;
int **angle_type = atom->angle_type;
int *mask = atom->mask;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int nlocal = atom->nlocal;
int molecular = atom->molecular;
// loop over all atoms and their angles
Angle *angle = force->angle;
double duang, du2ang;
double dx[3] = {0.0, 0.0, 0.0};
double dx_left[3] = {0.0, 0.0, 0.0};
double dx_right[3] = {0.0, 0.0, 0.0};
double x_angle_left[3] = {0.0, 0.0, 0.0};
double x_angle_middle[3] = {0.0, 0.0, 0.0};
double x_angle_right[3] = {0.0, 0.0, 0.0};
double dcos_theta[3] = {0.0, 0.0, 0.0};
double local_contribution[3] = {0.0, 0.0, 0.0};
// initialization
for (int i = 0; i < nvalues; i++) angle_local[i] = 0.0;
for (atom2 = 0; atom2 < nlocal; atom2++) {
if (!(mask[atom2] & groupbit)) continue;
if (molecular == 1)
na = num_angle[atom2];
else {
if (molindex[atom2] < 0) continue;
imol = molindex[atom2];
iatom = molatom[atom2];
na = onemols[imol]->num_angle[iatom];
}
for (int i = 0; i < na; i++) {
if (molecular == 1) {
if (tag[atom2] != angle_atom2[atom2][i]) continue;
atype = angle_type[atom2][i];
atom1 = atom->map(angle_atom1[atom2][i]);
atom3 = atom->map(angle_atom3[atom2][i]);
} else {
if (tag[atom2] != onemols[imol]->angle_atom2[atom2][i]) continue;
atype = onemols[imol]->angle_type[atom2][i];
tagprev = tag[atom2] - iatom - 1;
atom1 = atom->map(onemols[imol]->angle_atom1[atom2][i] + tagprev);
atom3 = atom->map(onemols[imol]->angle_atom3[atom2][i] + tagprev);
}
if (atom1 < 0 || !(mask[atom1] & groupbit)) continue;
if (atom3 < 0 || !(mask[atom3] & groupbit)) continue;
if (atype <= 0) continue;
// minimum image of atom1 with respect to the plane of interest
dx[0] = x[atom1][0];
dx[1] = x[atom1][1];
dx[2] = x[atom1][2];
dx[dir] -= pos;
domain->minimum_image(dx[0], dx[1], dx[2]);
x_angle_left[0] = dx[0];
x_angle_left[1] = dx[1];
x_angle_left[2] = dx[2];
x_angle_left[dir] += pos;
// minimum image of atom2 with respect to atom1
dx_left[0] = x[atom2][0] - x_angle_left[0];
dx_left[1] = x[atom2][1] - x_angle_left[1];
dx_left[2] = x[atom2][2] - x_angle_left[2];
domain->minimum_image(dx_left[0], dx_left[1], dx_left[2]);
x_angle_middle[0] = x_angle_left[0] + dx_left[0];
x_angle_middle[1] = x_angle_left[1] + dx_left[1];
x_angle_middle[2] = x_angle_left[2] + dx_left[2];
// minimum image of atom3 with respect to atom2
dx_right[0] = x[atom3][0] - x_angle_middle[0];
dx_right[1] = x[atom3][1] - x_angle_middle[1];
dx_right[2] = x[atom3][2] - x_angle_middle[2];
domain->minimum_image(dx_right[0], dx_right[1], dx_right[2]);
x_angle_right[0] = x_angle_middle[0] + dx_right[0];
x_angle_right[1] = x_angle_middle[1] + dx_right[1];
x_angle_right[2] = x_angle_middle[2] + dx_right[2];
// check if any bond vector crosses the plane of interest
double tau_right = (x_angle_right[dir] - pos) / (x_angle_right[dir] - x_angle_middle[dir]);
double tau_left = (x_angle_middle[dir] - pos) / (x_angle_middle[dir] - x_angle_left[dir]);
bool right_cross = ((tau_right >=0) && (tau_right <= 1));
bool left_cross = ((tau_left >=0) && (tau_left <= 1));
// no bonds crossing the plane
if (!right_cross && !left_cross) continue;
// compute the cos(theta) of the angle
r1 = sqrt(dx_left[0]*dx_left[0] + dx_left[1]*dx_left[1] + dx_left[2]*dx_left[2]);
r2 = sqrt(dx_right[0]*dx_right[0] + dx_right[1]*dx_right[1] + dx_right[2]*dx_right[2]);
cos_theta = -(dx_right[0]*dx_left[0] + dx_right[1]*dx_left[1] + dx_right[2]*dx_left[2])/(r1*r2);
if (cos_theta > 1.0) cos_theta = 1.0;
if (cos_theta < -1.0) cos_theta = -1.0;
// The method returns derivative with regards to cos(theta)
angle->born_matrix(atype, atom1, atom2, atom3, duang, du2ang);
// only right bond crossing the plane
if (right_cross && !left_cross)
{
double sgn = copysign(1.0, x_angle_right[dir] - pos);
dcos_theta[0] = sgn*(dx_right[0]*cos_theta/r2 + dx_left[0]/r1)/r2;
dcos_theta[1] = sgn*(dx_right[1]*cos_theta/r2 + dx_left[1]/r1)/r2;
dcos_theta[2] = sgn*(dx_right[2]*cos_theta/r2 + dx_left[2]/r1)/r2;
}
// only left bond crossing the plane
if (!right_cross && left_cross)
{
double sgn = copysign(1.0, x_angle_left[dir] - pos);
dcos_theta[0] = -sgn*(dx_left[0]*cos_theta/r1 + dx_right[0]/r2)/r1;
dcos_theta[1] = -sgn*(dx_left[1]*cos_theta/r1 + dx_right[1]/r2)/r1;
dcos_theta[2] = -sgn*(dx_left[2]*cos_theta/r1 + dx_right[2]/r2)/r1;
}
// both bonds crossing the plane
if (right_cross && left_cross)
{
// due to right bond
double sgn = copysign(1.0, x_angle_middle[dir] - pos);
dcos_theta[0] = -sgn*(dx_right[0]*cos_theta/r2 + dx_left[0]/r1)/r2;
dcos_theta[1] = -sgn*(dx_right[1]*cos_theta/r2 + dx_left[1]/r1)/r2;
dcos_theta[2] = -sgn*(dx_right[2]*cos_theta/r2 + dx_left[2]/r1)/r2;
// due to left bond
dcos_theta[0] += sgn*(dx_left[0]*cos_theta/r1 + dx_right[0]/r2)/r1;
dcos_theta[1] += sgn*(dx_left[1]*cos_theta/r1 + dx_right[1]/r2)/r1;
dcos_theta[2] += sgn*(dx_left[2]*cos_theta/r1 + dx_right[2]/r2)/r1;
}
// final contribution of the given angle term
local_contribution[0] += duang*dcos_theta[0]/area*nktv2p;
local_contribution[1] += duang*dcos_theta[1]/area*nktv2p;
local_contribution[2] += duang*dcos_theta[2]/area*nktv2p;
}
}
// loop over the keywords and if necessary add the angle contribution
int m = 0;
while (m < nvalues) {
if (which[m] == CONF || which[m] == TOTAL || which[m] == ANGLE) {
angle_local[m] = local_contribution[0];
angle_local[m+1] = local_contribution[1];
angle_local[m+2] = local_contribution[2];
}
m += 3;
}
}

8
src/EXTRA-COMPUTE/compute_stress_mop.h
View File

@ -38,11 +38,17 @@ class ComputeStressMop : public Compute {
private:
void compute_pairs();
void compute_bonds();
void compute_angles();
int me, nvalues, dir;
int nvalues, dir;
int *which;
int bondflag, angleflag;
double *values_local, *values_global;
double *bond_local, *bond_global;
double *angle_local, *angle_global;
double pos, pos1, dt, nktv2p, ftm2v;
double area;
class NeighList *list;

249
src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -14,47 +13,51 @@
/*------------------------------------------------------------------------
Contributing Authors : Romain Vermorel (LFCR), Laurent Joly (ULyon)
Support for bonds added by : Evangelos Voyiatzis (NovaMechanics)
--------------------------------------------------------------------------*/
#include "compute_stress_mop_profile.h"
#include "atom.h"
#include "update.h"
#include "atom_vec.h"
#include "bond.h"
#include "comm.h"
#include "domain.h"
#include "neighbor.h"
#include "force.h"
#include "pair.h"
#include "neigh_list.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "molecule.h"
#include "neigh_list.h"
#include "neighbor.h"
#include "pair.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
enum{X,Y,Z};
enum{LOWER,CENTER,UPPER,COORD};
enum{TOTAL,CONF,KIN};
#define BIG 1000000000
enum { X, Y, Z };
enum { LOWER, CENTER, UPPER, COORD };
enum { TOTAL, CONF, KIN, PAIR, BOND };
// clang-format off
/* ---------------------------------------------------------------------- */
ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg < 7) error->all(FLERR,"Illegal compute stress/mop/profile command");
if (narg < 7) utils::missing_cmd_args(FLERR, "compute stress/mop/profile", error);
MPI_Comm_rank(world,&me);
bondflag = 0;
// set compute mode and direction of plane(s) for pressure calculation
if (strcmp(arg[3],"x")==0) {
if (strcmp(arg[3],"x") == 0) {
dir = X;
} else if (strcmp(arg[3],"y")==0) {
} else if (strcmp(arg[3],"y") == 0) {
dir = Y;
} else if (strcmp(arg[3],"z")==0) {
} else if (strcmp(arg[3],"z") == 0) {
dir = Z;
} else error->all(FLERR,"Illegal compute stress/mop/profile command");
@ -64,8 +67,7 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a
else if (strcmp(arg[4],"center") == 0) originflag = CENTER;
else if (strcmp(arg[4],"upper") == 0) originflag = UPPER;
else originflag = COORD;
if (originflag == COORD)
origin = utils::numeric(FLERR,arg[4],false,lmp);
if (originflag == COORD) origin = utils::numeric(FLERR,arg[4],false,lmp);
delta = utils::numeric(FLERR,arg[5],false,lmp);
invdelta = 1.0/delta;
@ -92,6 +94,16 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a
which[nvalues] = TOTAL;
nvalues++;
}
} else if (strcmp(arg[iarg],"pair") == 0) {
for (i=0; i<3; i++) {
which[nvalues] = PAIR;
nvalues++;
}
} else if (strcmp(arg[iarg],"bond") == 0) {
for (i=0; i<3; i++) {
which[nvalues] = BOND;
nvalues++;
}
} else error->all(FLERR, "Illegal compute stress/mop/profile command"); //break;
iarg++;
@ -114,6 +126,9 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a
nbins = 0;
coord = coordp = nullptr;
values_local = values_global = array = nullptr;
bond_local = nullptr;
bond_global = nullptr;
local_contribution = nullptr;
// bin setup
@ -133,13 +148,15 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a
ComputeStressMopProfile::~ComputeStressMopProfile()
{
delete [] which;
delete[] which;
memory->destroy(coord);
memory->destroy(coordp);
memory->destroy(values_local);
memory->destroy(values_global);
memory->destroy(bond_local);
memory->destroy(bond_global);
memory->destroy(local_contribution);
memory->destroy(array);
}
@ -147,7 +164,6 @@ ComputeStressMopProfile::~ComputeStressMopProfile()
void ComputeStressMopProfile::init()
{
// conversion constants
nktv2p = force->nktv2p;
@ -173,27 +189,30 @@ void ComputeStressMopProfile::init()
//This compute requires a pair style with pair_single method implemented
if (force->pair == nullptr)
if (!force->pair)
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
// Errors
if (me==0) {
if (comm->me == 0) {
//Compute stress/mop/profile only accounts for pair interactions.
// issue a warning if any intramolecular potential or Kspace is defined.
// issue an error if any intramolecular potential or Kspace is defined.
if (force->bond!=nullptr)
error->warning(FLERR,"compute stress/mop/profile does not account for bond potentials");
if (force->angle!=nullptr)
error->warning(FLERR,"compute stress/mop/profile does not account for angle potentials");
if (force->dihedral!=nullptr)
error->warning(FLERR,"compute stress/mop/profile does not account for dihedral potentials");
if (force->improper!=nullptr)
error->warning(FLERR,"compute stress/mop/profile does not account for improper potentials");
if (force->kspace!=nullptr)
if (force->bond) bondflag = 1;
if (force->angle)
if ((strcmp(force->angle_style, "zero") != 0) && (strcmp(force->angle_style, "none") != 0))
error->all(FLERR,"compute stress/mop/profile does not account for angle potentials");
if (force->dihedral)
if ((strcmp(force->dihedral_style, "zero") != 0) && (strcmp(force->dihedral_style, "none") != 0))
error->all(FLERR,"compute stress/mop/profile does not account for dihedral potentials");
if (force->improper)
if ((strcmp(force->improper_style, "zero") != 0) && (strcmp(force->improper_style, "none") != 0))
error->all(FLERR,"compute stress/mop/profile does not account for improper potentials");
if (force->kspace)
error->warning(FLERR,"compute stress/mop/profile does not account for kspace contributions");
}
@ -222,17 +241,29 @@ void ComputeStressMopProfile::compute_array()
compute_pairs();
// sum pressure contributions over all procs
MPI_Allreduce(&values_local[0][0],&values_global[0][0],nbins*nvalues,
MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&values_local[0][0],&values_global[0][0],nbins*nvalues,MPI_DOUBLE,MPI_SUM,world);
int ibin,m,mo;
for (ibin=0; ibin<nbins; ibin++) {
if (bondflag) {
//Compute bond contribution on separate procs
compute_bonds();
} else {
for (int m = 0; m < nbins; m++) {
for (int i = 0; i < nvalues; i++) {
bond_local[m][i] = 0.0;
}
}
}
// sum bond contribution over all procs
MPI_Allreduce(&bond_local[0][0],&bond_global[0][0],nbins*nvalues,MPI_DOUBLE,MPI_SUM,world);
for (int ibin=0; ibin<nbins; ibin++) {
array[ibin][0] = coord[ibin][0];
mo=1;
m = 0;
while (m<nvalues) {
array[ibin][m+mo] = values_global[ibin][m];
int mo = 1;
int m = 0;
while (m < nvalues) {
array[ibin][m+mo] = values_global[ibin][m] + bond_global[ibin][m];
m++;
}
}
@ -289,8 +320,8 @@ void ComputeStressMopProfile::compute_pairs()
double xj[3];
m = 0;
while (m<nvalues) {
if (which[m] == CONF || which[m] == TOTAL) {
while (m < nvalues) {
if ((which[m] == CONF) || (which[m] == TOTAL) || (which[m] == PAIR)) {
// Compute configurational contribution to pressure
@ -377,7 +408,7 @@ void ComputeStressMopProfile::compute_pairs()
// compute kinetic contribution to pressure
// counts local particles transfers across the plane
if (which[m] == KIN || which[m] == TOTAL) {
if ((which[m] == KIN) || (which[m] == TOTAL)) {
double sgn;
@ -415,7 +446,7 @@ void ComputeStressMopProfile::compute_pairs()
xj[2] = xi[2]-vi[2]*dt+fi[2]/2/mass[itype]*dt*dt*ftm2v;
}
for (ibin=0;ibin<nbins;ibin++) {
for (ibin = 0; ibin < nbins; ibin++) {
pos = coord[ibin][0];
pos1 = coordp[ibin][0];
@ -452,6 +483,129 @@ void ComputeStressMopProfile::compute_pairs()
}
}
/*------------------------------------------------------------------------
compute bond contribution to pressure of local proc
-------------------------------------------------------------------------*/
void ComputeStressMopProfile::compute_bonds()
{
int i, nb, atom1, atom2, imol, iatom, btype;
tagint tagprev;
double rsq, fpair;
double **x = atom->x;
tagint *tag = atom->tag;
int *num_bond = atom->num_bond;
tagint **bond_atom = atom->bond_atom;
int **bond_type = atom->bond_type;
int *mask = atom->mask;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
int molecular = atom->molecular;
Bond *bond = force->bond;
double dx[3] = {0.0, 0.0, 0.0};
double x_bond_1[3] = {0.0, 0.0, 0.0};
double x_bond_2[3] = {0.0, 0.0, 0.0};
// initialization
for (int m = 0; m < nbins; m++) {
for (int i = 0; i < nvalues; i++) {
bond_local[m][i] = 0.0;
}
local_contribution[m][0] = 0.0;
local_contribution[m][1] = 0.0;
local_contribution[m][2] = 0.0;
}
// loop over all bonded atoms in the current proc
for (atom1 = 0; atom1 < nlocal; atom1++) {
if (!(mask[atom1] & groupbit)) continue;
if (molecular == 1)
nb = num_bond[atom1];
else {
if (molindex[atom1] < 0) continue;
imol = molindex[atom1];
iatom = molatom[atom1];
nb = onemols[imol]->num_bond[iatom];
}
for (i = 0; i < nb; i++) {
if (molecular == 1) {
btype = bond_type[atom1][i];
atom2 = atom->map(bond_atom[atom1][i]);
} else {
tagprev = tag[atom1] - iatom - 1;
btype = onemols[imol]->bond_type[iatom][i];
atom2 = atom->map(onemols[imol]->bond_atom[iatom][i] + tagprev);
}
if (atom2 < 0 || !(mask[atom2] & groupbit)) continue;
if (newton_bond == 0 && tag[atom1] > tag[atom2]) continue;
if (btype <= 0) continue;
for (int ibin = 0; ibin<nbins; ibin++) {
double pos = coord[ibin][0];
// minimum image of atom1 with respect to the plane of interest
dx[0] = x[atom1][0];
dx[1] = x[atom1][1];
dx[2] = x[atom1][2];
dx[dir] -= pos;
domain->minimum_image(dx[0], dx[1], dx[2]);
x_bond_1[0] = dx[0];
x_bond_1[1] = dx[1];
x_bond_1[2] = dx[2];
x_bond_1[dir] += pos;
// minimum image of atom2 with respect to atom1
dx[0] = x[atom2][0] - x_bond_1[0];
dx[1] = x[atom2][1] - x_bond_1[1];
dx[2] = x[atom2][2] - x_bond_1[2];
domain->minimum_image(dx[0], dx[1], dx[2]);
x_bond_2[0] = x_bond_1[0] + dx[0];
x_bond_2[1] = x_bond_1[1] + dx[1];
x_bond_2[2] = x_bond_1[2] + dx[2];
// check if the bond vector crosses the plane of interest
double tau = (x_bond_1[dir] - pos) / (x_bond_1[dir] - x_bond_2[dir]);
if ((tau <= 1) && (tau >= 0)) {
dx[0] = x_bond_1[0] - x_bond_2[0];
dx[1] = x_bond_1[1] - x_bond_2[1];
dx[2] = x_bond_1[2] - x_bond_2[2];
rsq = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
bond->single(btype, rsq, atom1, atom2, fpair);
double sgn = copysign(1.0, x_bond_1[dir] - pos);
local_contribution[ibin][0] += sgn*fpair*dx[0]/area*nktv2p;
local_contribution[ibin][1] += sgn*fpair*dx[1]/area*nktv2p;
local_contribution[ibin][2] += sgn*fpair*dx[2]/area*nktv2p;
}
}
}
}
// loop over the keywords and if necessary add the bond contribution
int m = 0;
while (m < nvalues) {
if ((which[m] == CONF) || (which[m] == TOTAL) || (which[m] == BOND)) {
for (int ibin = 0; ibin < nbins; ibin++) {
bond_local[ibin][m] = local_contribution[ibin][0];
bond_local[ibin][m+1] = local_contribution[ibin][1];
bond_local[ibin][m+2] = local_contribution[ibin][2];
}
}
m += 3;
}
}
/* ----------------------------------------------------------------------
setup 1d bins and their extent and coordinates
called at init()
@ -492,6 +646,9 @@ void ComputeStressMopProfile::setup_bins()
memory->create(coordp,nbins,1,"stress/mop/profile:coordp");
memory->create(values_local,nbins,nvalues,"stress/mop/profile:values_local");
memory->create(values_global,nbins,nvalues,"stress/mop/profile:values_global");
memory->create(bond_local,nbins,nvalues,"stress/mop/profile:bond_local");
memory->create(bond_global,nbins,nvalues,"stress/mop/profile:bond_global");
memory->create(local_contribution,nbins,3,"stress/mop/profile:local_contribution");
// set bin coordinates
for (i = 0; i < nbins; i++) {

7
src/EXTRA-COMPUTE/compute_stress_mop_profile.h
View File

@ -38,16 +38,21 @@ class ComputeStressMopProfile : public Compute {
private:
void compute_pairs();
void compute_bonds();
void setup_bins();
int me, nvalues, dir;
int nvalues, dir;
int *which;
int bondflag;
int originflag;
double origin, delta, offset, invdelta;
int nbins;
double **coord, **coordp;
double **values_local, **values_global;
double **bond_local, **bond_global;
double **local_contribution;
double dt, nktv2p, ftm2v;
double area;

8
unittest/commands/test_compute_global.cpp
View File

@ -169,7 +169,7 @@ TEST_F(ComputeGlobalTest, Geometry)
command("compute mom1 all momentum");
command("compute mom2 allwater momentum");
command("compute mop1 all stress/mop x 0.0 total");
command("compute mop2 all stress/mop/profile z lower 0.5 kin conf");
command("compute mop2 all stress/mop/profile z lower 0.5 kin pair");
thermo_style += " c_mu1 c_mu2 c_mop1[*] c_mop2[1][1]";
}
@ -225,9 +225,9 @@ TEST_F(ComputeGlobalTest, Geometry)
EXPECT_DOUBLE_EQ(mom2[0], -0.022332069630161717);
EXPECT_DOUBLE_EQ(mom2[1], -0.056896553865696115);
EXPECT_DOUBLE_EQ(mom2[2], 0.069179891052881484);
EXPECT_DOUBLE_EQ(mop1[0], 3522311.3572200728);
EXPECT_DOUBLE_EQ(mop1[1], 2871104.9055934539);
EXPECT_DOUBLE_EQ(mop1[2], -4136077.5224247416);
EXPECT_DOUBLE_EQ(mop1[0], 3536584.0880458541);
EXPECT_DOUBLE_EQ(mop1[1], 2887485.033995091);
EXPECT_DOUBLE_EQ(mop1[2], -4154145.8952306858);
EXPECT_DOUBLE_EQ(mop2[0][0], -8.0869239999999998);
EXPECT_DOUBLE_EQ(mop2[0][1], 0.0);
EXPECT_DOUBLE_EQ(mop2[0][2], 0.0);