diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst
index 4ad2261bb0..5662fc76d4 100644
--- a/doc/src/compute_stress_mop.rst
+++ b/doc/src/compute_stress_mop.rst
@@ -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
 """""""
diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp
index 60f2d76e06..98e3bf7043 100644
--- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp
+++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp
@@ -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;
+  }
+}
diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.h b/src/EXTRA-COMPUTE/compute_stress_mop.h
index 5357d36371..86140dc278 100644
--- a/src/EXTRA-COMPUTE/compute_stress_mop.h
+++ b/src/EXTRA-COMPUTE/compute_stress_mop.h
@@ -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;
diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp
index 709f109cb1..1f1dc30733 100644
--- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp
+++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp
@@ -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++) {
diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h
index b58f762c12..2b0ffef0f8 100644
--- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h
+++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h
@@ -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;
diff --git a/unittest/commands/test_compute_global.cpp b/unittest/commands/test_compute_global.cpp
index 6c365c8c2b..04913f51b5 100644
--- a/unittest/commands/test_compute_global.cpp
+++ b/unittest/commands/test_compute_global.cpp
@@ -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);