From a797a0d193cacfad1ca0329927ddb52f19799f58 Mon Sep 17 00:00:00 2001
From: RomainVermorel <romain.vermorel@univ-pau.fr>
Date: Tue, 4 Sep 2018 14:02:19 +0200
Subject: [PATCH] changed computes names to stress/mop and stress/mop/profile
---
doc/src/compute_stress_mop.txt | 111 +++++
examples/USER/mop/in.compute_stress_mop | 62 +++
src/USER-MOP/README | 10 +-
src/USER-MOP/compute_stress_mop.cpp | 424 +++++++++++++++++
src/USER-MOP/compute_stress_mop.h | 102 ++++
src/USER-MOP/compute_stress_mop_profile.cpp | 496 ++++++++++++++++++++
src/USER-MOP/compute_stress_mop_profile.h | 113 +++++
7 files changed, 1313 insertions(+), 5 deletions(-)
create mode 100644 doc/src/compute_stress_mop.txt
create mode 100755 examples/USER/mop/in.compute_stress_mop
create mode 100644 src/USER-MOP/compute_stress_mop.cpp
create mode 100644 src/USER-MOP/compute_stress_mop.h
create mode 100644 src/USER-MOP/compute_stress_mop_profile.cpp
create mode 100644 src/USER-MOP/compute_stress_mop_profile.h
diff --git a/doc/src/compute_stress_mop.txt b/doc/src/compute_stress_mop.txt
new file mode 100644
index 0000000000..8d95ef42e2
--- /dev/null
+++ b/doc/src/compute_stress_mop.txt
@@ -0,0 +1,111 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+compute stress/mop command :h3
+compute stress/mop/profile command :h3
+
+
+[Syntax:]
+
+compute ID group-ID style dir args keywords ... :pre
+
+ID, group-ID are documented in "compute"_compute.html command
+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) :ul
+ {stress/mop} args = pos
+ pos = {lower} or {center} or {upper} or coordinate value (distance units) is the position of the plane
+ {stress/mop/profile} args = origin delta
+ origin = {lower} or {center} or {upper} or coordinate value (distance units) is the position of the first plane
+ delta = value (distance units) is the distance between planes :pre
+
+compute 1 all stress/mop x lower total
+compute 1 liquid stress/mop z 0.0 kin conf
+fix 1 all ave/time 10 1000 10000 c_1\[*\] file mop.time
+fix 1 all ave/time 10 1000 10000 c_1\[2\] file mop.time :pre
+
+compute 1 all stress/mop/profile x lower 0.1 total
+compute 1 liquid stress/mop/profile z 0.0 0.25 kin conf
+fix 1 all ave/time 500 20 10000 c_1\[*\] ave running overwrite file mopp.time mode vector :pre
+
+
+[Description:]
+
+Compute {stress/mop} and compute {stress/mop/profile} define computations that
+calculate components of the local stress tensor using the method of
+planes "(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
+"compute stress/atom"_compute_stress_atom.html), the method of planes is
+compatible with mechanical balance in heterogeneous systems and at
+interfaces "(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
+"(Todd)"_#mop-todd. For the kinetic part, the algorithm considers that
+atoms have crossed the plane if their positions at times t-dt and t are
+one on either side of the plane, and uses the velocity at time t-dt/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).
+
+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 pressure added by the "pair_modify tail
+yes"_pair_modify.html command, since those are contributions to the global system pressure.
+
+[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/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 "units"_units.html.
+
+The values produced by this compute can be accessed by various "output commands"_Howto_output.html. For instance, the results can be written to a file using the "fix ave/time"_fix_ave_time.html command. Please see the example in the examples/USER/mop folder.
+
+[Restrictions:]
+
+This style is part of the USER-MOP package. It is only enabled if LAMMPS
+is built with that package. See the "Build package"_Build_package.html
+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.
+
+[Related commands:]
+
+"compute stress/atom"_compute_stress_atom.html
+
+[Default:] none
+
+:line
+
+:link(mop-todd)
+[(Todd)] B. D. Todd, Denis J. Evans, and Peter J. Daivis: "Pressure tensor for inhomogeneous fluids",
+Phys. Rev. E 52, 1627 (1995).
diff --git a/examples/USER/mop/in.compute_stress_mop b/examples/USER/mop/in.compute_stress_mop
new file mode 100755
index 0000000000..cfa97f5dfe
--- /dev/null
+++ b/examples/USER/mop/in.compute_stress_mop
@@ -0,0 +1,62 @@
+variable T equal 0.8
+variable p_solid equal 0.05
+
+lattice fcc 1.0
+region box block 0.0 6.0 0.0 6.0 -2.0 12.0
+create_box 2 box
+
+mass * 1.0
+pair_style lj/cut 2.5
+pair_coeff * * 1.0 1.0
+pair_coeff 1 2 0.5 1.0
+pair_coeff 2 2 0.0 0.0
+neigh_modify delay 0
+
+region solid_bottom block INF INF INF INF -1.1 0.1
+region liquid block INF INF INF INF 1.1 8.9
+region solid_up block INF INF INF INF 9.9 11.1
+
+create_atoms 1 region liquid
+delete_atoms porosity liquid 0.26 88765
+group liquid region liquid
+
+create_atoms 2 region solid_bottom
+group solid_bottom region solid_bottom
+create_atoms 2 region solid_up
+group solid_up region solid_up
+group solid union solid_bottom solid_up
+
+variable faSolid equal ${p_solid}*lx*ly/count(solid_up)
+fix piston_up solid_up aveforce NULL NULL -${faSolid}
+fix freeze_up solid_up setforce 0.0 0.0 NULL
+fix freeze_bottom solid_bottom setforce 0.0 0.0 0.0
+fix nvesol solid nve
+compute Tliq liquid temp
+fix nvtliq liquid nvt temp $T $T 0.5
+fix_modify nvtliq temp Tliq
+
+thermo 10000
+thermo_modify flush yes temp Tliq
+
+dump 1 all atom 10000 dump.lammpstrj
+
+fix fxbal all balance 1000 1.05 shift z 10 1.05
+velocity liquid create $T 47298 dist gaussian rot yes
+run 50000
+undump 1
+reset_timestep 0
+
+compute bin_z liquid chunk/atom bin/1d z 0.0 0.1 units box
+
+compute liquidStress_ke liquid stress/atom NULL ke
+compute liquidStress_vir liquid stress/atom NULL virial
+fix profile_z liquid ave/chunk 10 1000 10000 bin_z density/number temp c_liquidStress_ke[1] c_liquidStress_ke[2] c_liquidStress_ke[3] c_liquidStress_ke[4] c_liquidStress_ke[5] c_liquidStress_ke[6] c_liquidStress_vir[1] c_liquidStress_vir[2] c_liquidStress_vir[3] c_liquidStress_vir[4] c_liquidStress_vir[5] c_liquidStress_vir[6] ave running overwrite file profile.z
+
+compute mopz0 all stress/mop z center kin conf
+fix mopz0t all ave/time 10 1000 10000 c_mopz0[*] file mopz0.time
+
+compute moppz liquid stress/mop/profile z 0.0 0.1 kin conf
+fix moppzt all ave/time 100 100 10000 c_moppz[*] ave running overwrite file moppz.time mode vector
+
+run 40000
+
diff --git a/src/USER-MOP/README b/src/USER-MOP/README
index d8e7d2f732..15fcb3701a 100644
--- a/src/USER-MOP/README
+++ b/src/USER-MOP/README
@@ -1,21 +1,21 @@
USER-MOP
============
-This package provides the mop and mop/profile compute styles.
-See the doc page for compute mop or compute mop/profile command for how to use
+This package provides the stress/mop and stress/mop/profile compute styles.
+See the doc page for compute stress/mop or compute stress/mop/profile command for how to use
them.
PACKAGE DESCRIPTION
-------------------
-This is a LAMMPS (http://lammps.sandia.gov/) compute style that calculates
+These are LAMMPS (http://lammps.sandia.gov/) compute styles that calculate
components of the local stress tensor using the method of planes as described in
the paper by Todd et al. (B. D. Todd, Denis J. Evans, and Peter J. Daivis: "Pressure
tensor for inhomogeneous fluids", Phys. Rev. E 52, 1627 (1995)).
This package contains the source files of two compute styles:
-* compute mop calculates components of the local stress tensor using the method of planes, applied on one plane.
-* compute mop/profile calculates profiles of components of the local stress tensor using the method of planes, applied on an array of regularly spaced planes.
+* compute stress/mop calculates components of the local stress tensor using the method of planes, applied on one plane.
+* compute stress/mop/profile calculates profiles of components of the local stress tensor using the method of planes, applied on an array of regularly spaced planes.
The persons who created these files are Laurent Joly at University of Lyon 1 and Romain Vermorel at University of Pau and Pays de l'Adour.
Contact them directly if you have questions.
diff --git a/src/USER-MOP/compute_stress_mop.cpp b/src/USER-MOP/compute_stress_mop.cpp
new file mode 100644
index 0000000000..f764e0ee7e
--- /dev/null
+++ b/src/USER-MOP/compute_stress_mop.cpp
@@ -0,0 +1,424 @@
+/* ----------------------------------------------------------------------
+ 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 Authors : Romain Vermorel (LFCR), Laurent Joly (ULyon)
+ --------------------------------------------------------------------------*/
+
+#include <mpi.h>
+#include <cmath>
+#include <cstring>
+#include <cstdlib>
+
+#include "compute_stress_mop.h"
+#include "atom.h"
+#include "update.h"
+#include "domain.h"
+#include "group.h"
+#include "modify.h"
+#include "fix.h"
+#include "neighbor.h"
+#include "force.h"
+#include "pair.h"
+#include "neigh_request.h"
+#include "neigh_list.h"
+#include "error.h"
+#include "memory.h"
+
+using namespace LAMMPS_NS;
+
+enum{X,Y,Z};
+enum{TOTAL,CONF,KIN};
+
+#define BIG 1000000000
+
+/* ---------------------------------------------------------------------- */
+
+ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) :
+ Compute(lmp, narg, arg)
+{
+ if (narg < 6) error->all(FLERR,"Illegal compute stress/mop command");
+
+ MPI_Comm_rank(world,&me);
+
+ // set compute mode and direction of plane(s) for pressure calculation
+
+ if (strcmp(arg[3],"x")==0) {
+ dir = X;
+ } else if (strcmp(arg[3],"y")==0) {
+ dir = Y;
+ } 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) {
+ pos = domain->boxlo[dir];
+ } else if (strcmp(arg[4],"upper")==0) {
+ pos = domain->boxhi[dir];
+ } 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 {
+ pos1 = pos - domain->prd[dir];
+ }
+
+ // parse values until one isn't recognized
+
+ which = new int[3*(narg-5)];
+ nvalues = 0;
+ int i;
+
+ int iarg=5;
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"conf") == 0) {
+ for (i=0; i<3; i++) {
+ which[nvalues] = CONF;
+ nvalues++;
+ }
+ } else if (strcmp(arg[iarg],"kin") == 0) {
+ for (i=0; i<3; i++) {
+ which[nvalues] = KIN;
+ nvalues++;
+ }
+ } else if (strcmp(arg[iarg],"total") == 0) {
+ for (i=0; i<3; i++) {
+ which[nvalues] = TOTAL;
+ nvalues++;
+ }
+ } else error->all(FLERR, "Illegal compute stress/mop command"); //break;
+
+ iarg++;
+ }
+
+ // Error check
+ // 3D only
+
+ if (domain->dimension < 3)
+ error->all(FLERR, "Compute stress/mop incompatible with simulation dimension");
+
+ // 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");
+
+
+ // Initialize some variables
+
+ values_local = values_global = vector = NULL;
+
+ // 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");
+ size_vector = nvalues;
+
+ vector_flag = 1;
+ extvector = 0;
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeStressMop::~ComputeStressMop()
+{
+
+ delete [] which;
+
+ memory->destroy(values_local);
+ memory->destroy(values_global);
+ memory->destroy(vector);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeStressMop::init()
+{
+
+ // Conversion constants
+
+ nktv2p = force->nktv2p;
+ ftm2v = force->ftm2v;
+
+ // Plane area
+
+ area = 1;
+ int i;
+ for (i=0; i<3; i++){
+ if (i!=dir) area = area*domain->prd[i];
+ }
+
+ // 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)
+ 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
+
+ if (me==0){
+
+ //Compute stress/mop only accounts for pair interactions.
+ // issue a warning if any intramolecular potential or Kspace is defined.
+
+ if (force->bond!=NULL)
+ error->warning(FLERR,"compute stress/mop does not account for bond potentials");
+ if (force->angle!=NULL)
+ error->warning(FLERR,"compute stress/mop does not account for angle potentials");
+ if (force->dihedral!=NULL)
+ error->warning(FLERR,"compute stress/mop does not account for dihedral potentials");
+ if (force->improper!=NULL)
+ error->warning(FLERR,"compute stress/mop does not account for improper potentials");
+ if (force->kspace!=NULL)
+ error->warning(FLERR,"compute stress/mop does not account for kspace contributions");
+ }
+
+ // need an occasional half neighbor list
+ int irequest = neighbor->request((void *) this);
+ neighbor->requests[irequest]->pair = 0;
+ neighbor->requests[irequest]->compute = 1;
+ neighbor->requests[irequest]->occasional = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeStressMop::init_list(int id, NeighList *ptr)
+{
+ list = ptr;
+}
+
+
+/* ----------------------------------------------------------------------
+ compute output vector
+ ------------------------------------------------------------------------- */
+
+void ComputeStressMop::compute_vector()
+{
+ invoked_array = update->ntimestep;
+
+ //Compute pressures on separate procs
+ compute_pairs();
+
+ // sum pressure contributions over all procs
+ 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];
+ }
+
+}
+
+
+/*------------------------------------------------------------------------
+ compute pressure contribution of local proc
+ -------------------------------------------------------------------------*/
+
+void ComputeStressMop::compute_pairs()
+
+{
+ int i,j,m,n,ii,jj,inum,jnum,itype,jtype;
+ double delx,dely,delz;
+ double rsq,eng,fpair,factor_coul,factor_lj;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ double *mass = atom->mass;
+ int *type = atom->type;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+ double *special_coul = force->special_coul;
+ double *special_lj = force->special_lj;
+ int newton_pair = force->newton_pair;
+
+
+ // zero out arrays for one sample
+
+ for (i = 0; i < nvalues; i++) values_local[i] = 0.0;
+
+ // invoke half neighbor list (will copy or build if necessary)
+
+ neighbor->build_one(list);
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // loop over neighbors of my atoms
+
+ Pair *pair = force->pair;
+ double **cutsq = force->pair->cutsq;
+
+ // Parse values
+
+ double xi[3];
+ double vi[3];
+ double fi[3];
+ double xj[3];
+
+ m = 0;
+ while (m<nvalues) {
+ if (which[m] == CONF || which[m] == TOTAL) {
+
+ //Compute configurational contribution to pressure
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+
+ xi[0] = atom->x[i][0];
+ xi[1] = atom->x[i][1];
+ xi[2] = atom->x[i][2];
+ itype = type[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ factor_lj = special_lj[sbmask(j)];
+ factor_coul = special_coul[sbmask(j)];
+ j &= NEIGHMASK;
+
+ // skip if neither I nor J are in group
+ if (!(mask[i] & groupbit || mask[j] & groupbit)) continue;
+
+ xj[0] = atom->x[j][0];
+ xj[1] = atom->x[j][1];
+ xj[2] = atom->x[j][2];
+ delx = xi[0] - xj[0];
+ dely = xi[1] - xj[1];
+ delz = xi[2] - xj[2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ jtype = type[j];
+ if (rsq >= cutsq[itype][jtype]) continue;
+
+ if (newton_pair || j < nlocal) {
+
+ //check if ij pair is accross plane, add contribution to pressure
+ 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)) ){
+
+ 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)) ) {
+
+ 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){
+ double vcm[3];
+ double masstotal,sgn;
+
+ for (int i = 0; i < nlocal; i++){
+
+ // skip if I is not in group
+ if (mask[i] & groupbit){
+
+ itype = type[i];
+
+ //coordinates at t
+ xi[0] = atom->x[i][0];
+ xi[1] = atom->x[i][1];
+ xi[2] = atom->x[i][2];
+
+ //velocities at t
+ vi[0] = atom->v[i][0];
+ vi[1] = atom->v[i][1];
+ vi[2] = atom->v[i][2];
+
+ //forces at t
+ fi[0] = atom->f[i][0];
+ fi[1] = atom->f[i][1];
+ 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;
+
+ // because LAMMPS does not put atoms back in the box
+ // at each timestep, must check atoms going through the
+ // image of the plane that is closest to the box
+
+ double pos_temp = pos+copysign(1,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)){
+
+ //sgn = copysign(1,vi[dir]-vcm[dir]);
+ sgn = copysign(1,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[m] += mass[itype]*vcross[0]*sgn/dt/area*nktv2p/ftm2v;
+ values_local[m+1] += mass[itype]*vcross[1]*sgn/dt/area*nktv2p/ftm2v;
+ values_local[m+2] += mass[itype]*vcross[2]*sgn/dt/area*nktv2p/ftm2v;
+ }
+ }
+ }
+ }
+ m+=3;
+ }
+}
+
diff --git a/src/USER-MOP/compute_stress_mop.h b/src/USER-MOP/compute_stress_mop.h
new file mode 100644
index 0000000000..2047d1d54f
--- /dev/null
+++ b/src/USER-MOP/compute_stress_mop.h
@@ -0,0 +1,102 @@
+/* ----------------------------------------------------------------------
+ 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 Authors : Romain Vermorel (LFCR), Laurent Joly (ULyon)
+ --------------------------------------------------------------------------*/
+
+#ifdef COMPUTE_CLASS
+
+ComputeStyle(stress/mop,ComputeStressMop)
+
+#else
+
+#ifndef LMP_COMPUTE_STRESS_MOP_H
+#define LMP_COMPUTE_STRESS_MOP_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+ class ComputeStressMop : public Compute {
+ public:
+ ComputeStressMop(class LAMMPS *, int, char **);
+ virtual ~ComputeStressMop();
+ void init();
+ void init_list(int, class NeighList *);
+ void compute_vector();
+
+ private:
+
+ void compute_pairs();
+
+ int me,nvalues,dir;
+ int *which;
+
+ double *values_local,*values_global;
+ double pos,pos1,dt,nktv2p,ftm2v;
+ double area;
+ class NeighList *list;
+
+ };
+
+}
+
+#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: Compute stress/mop incompatible with simulation dimension
+
+ Compute stress/mop only works with 3D simulations.
+
+ E: Compute stress/mop incompatible with triclinic simulation box
+
+ Self-explanatory.
+
+ E: Compute stress/mop requires a fixed simulation box
+
+ Compute stress/mop is not compatible with any change of volume or shape
+ or boundary conditions of the simulation box.
+
+ E: No pair style is defined for compute stress/mop
+
+ Self-explanatory. Compute stress/mop requires the definition of a pair style.
+
+ E: Pair style does not support compute stress/mop
+
+ The pair style does not have a single() function, so it can
+ not be invoked by compute stress/mop.
+
+ W: compute stress/mop does not account for bond potentials
+
+ W: compute stress/mop does not account for angle potentials
+
+ W: compute stress/mop does not account for dihedral potentials
+
+ W: compute stress/mop does not account for improper potentials
+
+ W: compute stress/mop does not account for kspace contributions
+
+ Compute stress/mop only accounts for pairwise additive interactions for
+ the computation of local stress tensor components.
+
+ */
+
diff --git a/src/USER-MOP/compute_stress_mop_profile.cpp b/src/USER-MOP/compute_stress_mop_profile.cpp
new file mode 100644
index 0000000000..7ba0622aa4
--- /dev/null
+++ b/src/USER-MOP/compute_stress_mop_profile.cpp
@@ -0,0 +1,496 @@
+/* ----------------------------------------------------------------------
+ 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 Authors : Romain Vermorel (LFCR), Laurent Joly (ULyon)
+ --------------------------------------------------------------------------*/
+
+#include <mpi.h>
+#include <cmath>
+#include <cstring>
+#include <cstdlib>
+
+#include "compute_stress_mop_profile.h"
+#include "atom.h"
+#include "update.h"
+#include "domain.h"
+#include "group.h"
+#include "modify.h"
+#include "fix.h"
+#include "neighbor.h"
+#include "force.h"
+#include "pair.h"
+#include "neigh_request.h"
+#include "neigh_list.h"
+#include "error.h"
+#include "memory.h"
+
+using namespace LAMMPS_NS;
+
+enum{X,Y,Z};
+enum{LOWER,CENTER,UPPER,COORD};
+enum{TOTAL,CONF,KIN};
+
+#define BIG 1000000000
+
+/* ---------------------------------------------------------------------- */
+
+ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **arg) :
+ Compute(lmp, narg, arg)
+{
+ if (narg < 7) error->all(FLERR,"Illegal compute stress/mop/profile command");
+
+ MPI_Comm_rank(world,&me);
+
+ // set compute mode and direction of plane(s) for pressure calculation
+
+ if (strcmp(arg[3],"x")==0) {
+ dir = X;
+ } else if (strcmp(arg[3],"y")==0) {
+ dir = Y;
+ } else if (strcmp(arg[3],"z")==0) {
+ dir = Z;
+ } else error->all(FLERR,"Illegal compute stress/mop/profile command");
+
+ // bin parameters
+
+ if (strcmp(arg[4],"lower") == 0) originflag = LOWER;
+ 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 = force->numeric(FLERR,arg[4]);
+ delta = force->numeric(FLERR,arg[5]);
+ invdelta = 1.0/delta;
+
+ // parse values until one isn't recognized
+
+ which = new int[3*(narg-6)];
+ nvalues = 0;
+ int i;
+
+ int iarg=6;
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"conf") == 0) {
+ for (i=0; i<3; i++) {
+ which[nvalues] = CONF;
+ nvalues++;
+ }
+ } else if (strcmp(arg[iarg],"kin") == 0) {
+ for (i=0; i<3; i++) {
+ which[nvalues] = KIN;
+ nvalues++;
+ }
+ } else if (strcmp(arg[iarg],"total") == 0) {
+ for (i=0; i<3; i++) {
+ which[nvalues] = TOTAL;
+ nvalues++;
+ }
+ } else error->all(FLERR, "Illegal compute stress/mop/profile command"); //break;
+
+ iarg++;
+ }
+
+ // check domain related errors
+
+ // 3D only
+
+ if (domain->dimension < 3)
+ error->all(FLERR, "Compute stress/mop/profile incompatible with simulation dimension");
+
+ // orthogonal simulation box
+
+ if (domain->triclinic != 0)
+ error->all(FLERR, "Compute stress/mop/profile incompatible with triclinic simulation box");
+
+ // initialize some variables
+
+ nbins = 0;
+ coord = coordp = NULL;
+ values_local = values_global = array = NULL;
+
+ // bin setup
+
+ setup_bins();
+
+ // this fix produces a global array
+
+ memory->create(array,nbins,1+nvalues,"stress/mop/profile:array");
+ size_array_rows = nbins;
+ size_array_cols = 1 + nvalues;
+
+ array_flag = 1;
+ extarray = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeStressMopProfile::~ComputeStressMopProfile()
+{
+
+ delete [] which;
+
+ memory->destroy(coord);
+ memory->destroy(coordp);
+ memory->destroy(values_local);
+ memory->destroy(values_global);
+ memory->destroy(array);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeStressMopProfile::init()
+{
+
+ // conversion constants
+
+ nktv2p = force->nktv2p;
+ ftm2v = force->ftm2v;
+
+ // plane area
+
+ area = 1;
+ int i;
+ for (i=0; i<3; i++) {
+ if (i!=dir) area = area*domain->prd[i];
+ }
+
+ // 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){
+
+ //Compute stress/mop/profile only accounts for pair interactions.
+ // issue a warning if any intramolecular potential or Kspace is defined.
+
+ if (force->bond!=NULL)
+ error->warning(FLERR,"compute stress/mop/profile does not account for bond potentials");
+ if (force->angle!=NULL)
+ error->warning(FLERR,"compute stress/mop/profile does not account for angle potentials");
+ if (force->dihedral!=NULL)
+ error->warning(FLERR,"compute stress/mop/profile does not account for dihedral potentials");
+ if (force->improper!=NULL)
+ error->warning(FLERR,"compute stress/mop/profile does not account for improper potentials");
+ if (force->kspace!=NULL)
+ error->warning(FLERR,"compute stress/mop/profile does not account for kspace contributions");
+ }
+
+ // need an occasional half neighbor list
+
+ int irequest = neighbor->request((void *) this);
+ neighbor->requests[irequest]->pair = 0;
+ neighbor->requests[irequest]->compute = 1;
+ neighbor->requests[irequest]->occasional = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeStressMopProfile::init_list(int id, NeighList *ptr)
+{
+ list = ptr;
+}
+
+
+/* ----------------------------------------------------------------------
+ compute output array
+ ------------------------------------------------------------------------- */
+
+void ComputeStressMopProfile::compute_array()
+{
+ invoked_array = update->ntimestep;
+
+ //Compute pressures on separate procs
+ 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);
+
+ int ibin,m,mo;
+ for (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];
+ m++;
+ }
+ }
+}
+
+
+/*------------------------------------------------------------------------
+ compute pressure contribution of local proc
+ -------------------------------------------------------------------------*/
+
+void ComputeStressMopProfile::compute_pairs()
+
+{
+ int i,j,m,n,ii,jj,inum,jnum,itype,jtype,ibin;
+ double delx,dely,delz;
+ double rsq,eng,fpair,factor_coul,factor_lj;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+ double pos,pos1,pos_temp;
+
+ double *mass = atom->mass;
+ int *type = atom->type;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+ double *special_coul = force->special_coul;
+ double *special_lj = force->special_lj;
+ int newton_pair = force->newton_pair;
+
+
+ // zero out arrays for one sample
+ for (m = 0; m < nbins; m++) {
+ for (i = 0; i < nvalues; i++) values_local[m][i] = 0.0;
+ }
+
+ // invoke half neighbor list (will copy or build if necessary)
+
+ neighbor->build_one(list);
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // loop over neighbors of my atoms
+
+ Pair *pair = force->pair;
+ double **cutsq = force->pair->cutsq;
+
+ // parse values
+
+ double xi[3];
+ double vi[3];
+ double fi[3];
+ double xj[3];
+
+ m = 0;
+ while (m<nvalues) {
+ if (which[m] == CONF || which[m] == TOTAL) {
+
+ // Compute configurational contribution to pressure
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+
+ xi[0] = atom->x[i][0];
+ xi[1] = atom->x[i][1];
+ xi[2] = atom->x[i][2];
+ itype = type[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ factor_lj = special_lj[sbmask(j)];
+ factor_coul = special_coul[sbmask(j)];
+ j &= NEIGHMASK;
+
+ // skip if neither I nor J are in group
+
+ if (!(mask[i] & groupbit || mask[j] & groupbit)) continue;
+
+ xj[0] = atom->x[j][0];
+ xj[1] = atom->x[j][1];
+ xj[2] = atom->x[j][2];
+ delx = xi[0] - xj[0];
+ dely = xi[1] - xj[1];
+ delz = xi[2] - xj[2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ jtype = type[j];
+ if (rsq >= cutsq[itype][jtype]) continue;
+
+ if (newton_pair || j < nlocal) {
+
+ for (ibin=0;ibin<nbins;ibin++) {
+ pos = coord[ibin][0];
+ pos1 = coordp[ibin][0];
+
+ //check if ij pair is accross plane, add contribution to pressure
+
+ 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[ibin][m] += fpair*(xi[0]-xj[0])/area*nktv2p;
+ values_local[ibin][m+1] += fpair*(xi[1]-xj[1])/area*nktv2p;
+ values_local[ibin][m+2] += fpair*(xi[2]-xj[2])/area*nktv2p;
+
+ } 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[ibin][m] -= fpair*(xi[0]-xj[0])/area*nktv2p;
+ values_local[ibin][m+1] -= fpair*(xi[1]-xj[1])/area*nktv2p;
+ values_local[ibin][m+2] -= fpair*(xi[2]-xj[2])/area*nktv2p;
+ }
+ }
+ } else {
+
+ for (ibin=0;ibin<nbins;ibin++) {
+ pos = coord[ibin][0];
+ pos1 = coordp[ibin][0];
+
+ //check if ij pair is accross plane, add contribution to pressure
+
+ 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[ibin][m] += fpair*(xi[0]-xj[0])/area*nktv2p;
+ values_local[ibin][m+1] += fpair*(xi[1]-xj[1])/area*nktv2p;
+ values_local[ibin][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){
+
+ double vcm[3];
+ double masstotal,sgn;
+
+ for (int i = 0; i < nlocal; i++){
+
+ // skip if I is not in group
+
+ if (mask[i] & groupbit){
+
+ itype = type[i];
+
+ //coordinates at t
+ xi[0] = atom->x[i][0];
+ xi[1] = atom->x[i][1];
+ xi[2] = atom->x[i][2];
+
+ //velocities at t
+ vi[0] = atom->v[i][0];
+ vi[1] = atom->v[i][1];
+ vi[2] = atom->v[i][2];
+
+ //forces at t
+ fi[0] = atom->f[i][0];
+ fi[1] = atom->f[i][1];
+ 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;
+
+ for (ibin=0;ibin<nbins;ibin++) {
+ pos = coord[ibin][0];
+ pos1 = coordp[ibin][0];
+
+ if (((xi[dir]-pos)*(xj[dir]-pos)*(xi[dir]-pos1)*(xj[dir]-pos1)<0)){
+
+ sgn = copysign(1,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;
+ }
+ }
+ }
+ }
+ }
+ m+=3;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ setup 1d bins and their extent and coordinates
+ called at init()
+ ------------------------------------------------------------------------- */
+
+void ComputeStressMopProfile::setup_bins()
+{
+ int i,j,k,m,n;
+ double lo,hi,coord1,coord2;
+
+ double *boxlo,*boxhi,*prd;
+ boxlo = domain->boxlo;
+ boxhi = domain->boxhi;
+ prd = domain->prd;
+
+ if (originflag == LOWER) origin = boxlo[dir];
+ else if (originflag == UPPER) origin = boxhi[dir];
+ else if (originflag == CENTER)
+ origin = 0.5 * (boxlo[dir] + boxhi[dir]);
+
+ if (origin < boxlo[dir]) {
+ error->all(FLERR,"Origin of bins for compute stress/mop/profile is out of bounds" );
+ } else {
+ n = static_cast<int> ((origin - boxlo[dir]) * invdelta);
+ lo = origin - n*delta;
+ }
+ if (origin < boxhi[dir]) {
+ n = static_cast<int> ((boxhi[dir] - origin) * invdelta);
+ hi = origin + n*delta;
+ } else {
+ error->all(FLERR,"Origin of bins for compute stress/mop/profile is out of bounds" );
+ }
+
+ offset = lo;
+ nbins = static_cast<int> ((hi-lo) * invdelta + 1.5);
+
+ //allocate bin arrays
+ memory->create(coord,nbins,1,"stress/mop/profile:coord");
+ 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");
+
+ // set bin coordinates
+ for (i = 0; i < nbins; i++) {
+ coord[i][0] = offset + i*delta;
+ if ( coord[i][0] < (domain->boxlo[dir]+domain->prd_half[dir]) ) {
+ coordp[i][0] = coord[i][0] + domain->prd[dir];
+ } else {
+ coordp[i][0] = coord[i][0] - domain->prd[dir];
+ }
+ }
+}
diff --git a/src/USER-MOP/compute_stress_mop_profile.h b/src/USER-MOP/compute_stress_mop_profile.h
new file mode 100644
index 0000000000..648d86cc7f
--- /dev/null
+++ b/src/USER-MOP/compute_stress_mop_profile.h
@@ -0,0 +1,113 @@
+/* ----------------------------------------------------------------------
+ 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 Authors : Romain Vermorel (LFCR), Laurent Joly (ULyon)
+ --------------------------------------------------------------------------*/
+
+#ifdef COMPUTE_CLASS
+
+ComputeStyle(stress/mop/profile,ComputeStressMopProfile)
+
+#else
+
+#ifndef LMP_COMPUTE_STRESS_MOP_PROFILE_H
+#define LMP_COMPUTE_STRESS_MOP_PROFILE_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+ class ComputeStressMopProfile : public Compute {
+ public:
+ ComputeStressMopProfile(class LAMMPS *, int, char **);
+ virtual ~ComputeStressMopProfile();
+ void init();
+ void init_list(int, class NeighList *);
+ void compute_array();
+
+ private:
+
+ void compute_pairs();
+ void setup_bins();
+
+ int me,nvalues,dir;
+ int *which;
+
+ int originflag;
+ double origin,delta,offset,invdelta;
+ int nbins;
+ double **coord,**coordp;
+ double **values_local,**values_global;
+
+ int ndim;
+ double dt,nktv2p,ftm2v;
+ double area;
+ class NeighList *list;
+
+ };
+
+}
+
+#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: Compute stress/mop/profile incompatible with simulation dimension
+
+ Compute stress/mop/profile only works with 3D simulations.
+
+ E: Compute stress/mop/profile incompatible with triclinic simulation box
+
+ Self-explanatory.
+
+ E: Compute stress/mop/profile requires a fixed simulation box
+
+ Compute stress/mop/profile is not compatible with any change of volume or shape
+ or boundary conditions of the simulation box.
+
+ E: No pair style is defined for compute stress/mop/profile
+
+ Self-explanatory. Compute stress/mop/profile requires the definition of a pair style.
+
+ E: Pair style does not support compute stress/mop/profile
+
+ The pair style does not have a single() function, so it can
+ not be invoked by compute stress/mop/profile.
+
+ E: Origin of bins for compute stress/mop/profile is out of bounds
+
+ Self-explanatory.
+
+ W: compute stress/mop/profile does not account for bond potentials
+
+ W: compute stress/mop/profile does not account for angle potentials
+
+ W: compute stress/mop/profile does not account for dihedral potentials
+
+ W: compute stress/mop/profile does not account for improper potentials
+
+ W: compute stress/mop/profile does not account for kspace contributions
+
+ Compute stress/mop/profile only accounts for pairwise additive interactions for
+ the computation of local stress tensor components.
+
+*/
+
--
GitLab