change energy tally during minimize

src/RIGID/fix_shake.cpp

@@ -67,6 +67,11 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
   centroidstressflag = CENTROID_AVAIL;
   next_output = -1;
 
+  // to avoid uninitialized access
+  vflag_post_force = 0;
+  eflag_pre_reverse = 0;
+  ebond = 0.0;
+
   // error check
 
   molecular = atom->molecular;
@@ -324,6 +329,7 @@ int FixShake::setmask()
   mask |= PRE_NEIGHBOR;
   mask |= POST_FORCE;
   mask |= POST_FORCE_RESPA;
+  mask |= MIN_PRE_REVERSE;
   mask |= MIN_POST_FORCE;
   mask |= POST_RUN;
   return mask;
@@ -505,6 +511,13 @@ void FixShake::min_setup(int vflag)
   min_post_force(vflag);
 }
 
+/* --------------------------------------------------------------------- */
+
+void FixShake::setup_pre_reverse(int eflag, int vflag)
+{
+  min_pre_reverse(eflag,vflag);
+}
+
 /* ----------------------------------------------------------------------
    build list of SHAKE clusters to constrain
    if one or more atoms in cluster are on this proc,
@@ -513,6 +526,7 @@ void FixShake::min_setup(int vflag)
 
 void FixShake::pre_neighbor()
 {
+  ebond = 0.0;
   int atom1,atom2,atom3,atom4;
 
   // local copies of atom quantities
@@ -653,6 +667,15 @@ void FixShake::post_force_respa(int vflag, int ilevel, int iloop)
   vflag_post_force = vflag;
 }
 
+/* ----------------------------------------------------------------------
+   store eflag so it can be used in min_post_force
+------------------------------------------------------------------------- */
+
+void FixShake::min_pre_reverse(int eflag, int /*vflag*/)
+{
+  eflag_pre_reverse = eflag;
+}
+
 /* ----------------------------------------------------------------------
    substitute shake constraints with very strong bonds
 ------------------------------------------------------------------------- */
@@ -668,6 +691,9 @@ void FixShake::min_post_force(int vflag)
       next_output = (ntimestep/output_every)*output_every + output_every;
   } else next_output = -1;
 
+  int eflag = eflag_pre_reverse;
+  ev_init(eflag, vflag);
+
   x = atom->x;
   f = atom->f;
   nlocal = atom->nlocal;
@@ -2522,18 +2548,33 @@ void FixShake::bond_force(tagint id1, tagint id2, double length)
   const double dr = r - length;
   const double rk = kbond * dr;
   const double fbond = (r > 0.0) ? -2.0 * rk / r : 0.0;
+  const double eb = rk*dr;
+  int list[2];
+  int nlist = 0;
 
   if (i1 < nlocal) {
     f[i1][0] += delx * fbond;
     f[i1][1] += dely * fbond;
     f[i1][2] += delz * fbond;
-    ebond += 0.5*rk*dr;
+    list[nlist++] = i1;
+    ebond += 0.5*eb;
   }
   if (i2 < nlocal) {
     f[i2][0] -= delx * fbond;
     f[i2][1] -= dely * fbond;
     f[i2][2] -= delz * fbond;
-    ebond += 0.5*rk*dr;
+    list[nlist++] = i2;
+    ebond += 0.5*eb;
+  }
+  if (evflag) {
+    double v[6];
+    v[0] = 0.5 * delx * delx * fbond;
+    v[1] = 0.5 * dely * dely * fbond;
+    v[2] = 0.5 * delz * delz * fbond;
+    v[3] = 0.5 * delx * dely * fbond;
+    v[4] = 0.5 * delx * delz * fbond;
+    v[5] = 0.5 * dely * delz * fbond;
+    ev_tally(nlist, list, 2.0, eb, v);
   }
 }
 

src/RIGID/fix_shake.h

@@ -34,10 +34,12 @@ class FixShake : public Fix {
   int setmask() override;
   void init() override;
   void setup(int) override;
+  void setup_pre_reverse(int, int) override;
   void min_setup(int) override;
   void pre_neighbor() override;
   void post_force(int) override;
   void post_force_respa(int, int, int) override;
+  void min_pre_reverse(int, int) override;
   void min_post_force(int) override;
   void post_run() override;
 
@@ -64,6 +66,7 @@ class FixShake : public Fix {
 
  protected:
   int vflag_post_force;    // store the vflag of last post_force call
+  int eflag_pre_reverse;   // store the eflag of last pre_reverse call
   int respa;               // 0 = vel. Verlet, 1 = respa
   int rattle;              // 0 = SHAKE, 1 = RATTLE
   double tolerance;        // SHAKE tolerance

src/compute_pe_atom.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -13,30 +12,31 @@
 ------------------------------------------------------------------------- */
 
 #include "compute_pe_atom.h"
-#include <cstring>
+
-#include "atom.h"
-#include "update.h"
-#include "comm.h"
-#include "force.h"
-#include "pair.h"
-#include "bond.h"
 #include "angle.h"
+#include "atom.h"
+#include "bond.h"
+#include "comm.h"
 #include "dihedral.h"
+#include "error.h"
+#include "force.h"
 #include "improper.h"
 #include "kspace.h"
-#include "modify.h"
 #include "memory.h"
-#include "error.h"
+#include "modify.h"
+#include "pair.h"
+#include "update.h"
+
+#include <cstring>
 
 using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
 ComputePEAtom::ComputePEAtom(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg),
+    Compute(lmp, narg, arg), energy(nullptr)
-  energy(nullptr)
 {
-  if (narg < 3) error->all(FLERR,"Illegal compute pe/atom command");
+  if (narg < 3) error->all(FLERR, "Illegal compute pe/atom command");
 
   peratom_flag = 1;
   size_peratom_cols = 0;
@@ -56,14 +56,22 @@ ComputePEAtom::ComputePEAtom(LAMMPS *lmp, int narg, char **arg) :
     fixflag = 0;
     int iarg = 3;
     while (iarg < narg) {
-      if (strcmp(arg[iarg],"pair") == 0) pairflag = 1;
+      if (strcmp(arg[iarg], "pair") == 0)
-      else if (strcmp(arg[iarg],"bond") == 0) bondflag = 1;
+        pairflag = 1;
-      else if (strcmp(arg[iarg],"angle") == 0) angleflag = 1;
+      else if (strcmp(arg[iarg], "bond") == 0)
-      else if (strcmp(arg[iarg],"dihedral") == 0) dihedralflag = 1;
+        bondflag = 1;
-      else if (strcmp(arg[iarg],"improper") == 0) improperflag = 1;
+      else if (strcmp(arg[iarg], "angle") == 0)
-      else if (strcmp(arg[iarg],"kspace") == 0) kspaceflag = 1;
+        angleflag = 1;
-      else if (strcmp(arg[iarg],"fix") == 0) fixflag = 1;
+      else if (strcmp(arg[iarg], "dihedral") == 0)
-      else error->all(FLERR,"Illegal compute pe/atom command");
+        dihedralflag = 1;
+      else if (strcmp(arg[iarg], "improper") == 0)
+        improperflag = 1;
+      else if (strcmp(arg[iarg], "kspace") == 0)
+        kspaceflag = 1;
+      else if (strcmp(arg[iarg], "fix") == 0)
+        fixflag = 1;
+      else
+        error->all(FLERR, "Illegal compute pe/atom command");
       iarg++;
     }
   }
@@ -86,7 +94,7 @@ void ComputePEAtom::compute_peratom()
 
   invoked_peratom = update->ntimestep;
   if (update->eflag_atom != invoked_peratom)
-    error->all(FLERR,"Per-atom energy was not tallied on needed timestep");
+    error->all(FLERR, "Per-atom energy was not tallied on needed timestep");
 
   // grow local energy array if necessary
   // needs to be atom->nmax in length
@@ -94,7 +102,7 @@ void ComputePEAtom::compute_peratom()
   if (atom->nmax > nmax) {
     memory->destroy(energy);
     nmax = atom->nmax;
-    memory->create(energy,nmax,"pe/atom:energy");
+    memory->create(energy, nmax, "pe/atom:energy");
     vector_atom = energy;
   }
 
@@ -153,13 +161,11 @@ void ComputePEAtom::compute_peratom()
   // add in per-atom contributions from relevant fixes
   // always only for owned atoms, not ghost
 
-  if (fixflag && modify->n_energy_atom)
+  if (fixflag && modify->n_energy_atom) modify->energy_atom(nlocal, energy);
-    modify->energy_atom(nlocal,energy);
 
   // communicate ghost energy between neighbor procs
 
-  if (force->newton || (force->kspace && force->kspace->tip4pflag))
+  if (force->newton || (force->kspace && force->kspace->tip4pflag)) comm->reverse_comm(this);
-    comm->reverse_comm(this);
 
   // zero energy of atoms not in group
   // only do this after comm since ghost contributions must be included
@@ -174,7 +180,7 @@ void ComputePEAtom::compute_peratom()
 
 int ComputePEAtom::pack_reverse_comm(int n, int first, double *buf)
 {
-  int i,m,last;
+  int i, m, last;
 
   m = 0;
   last = first + n;
@@ -186,7 +192,7 @@ int ComputePEAtom::pack_reverse_comm(int n, int first, double *buf)
 
 void ComputePEAtom::unpack_reverse_comm(int n, int *list, double *buf)
 {
-  int i,j,m;
+  int i, j, m;
 
   m = 0;
   for (i = 0; i < n; i++) {
@@ -201,6 +207,6 @@ void ComputePEAtom::unpack_reverse_comm(int n, int *list, double *buf)
 
 double ComputePEAtom::memory_usage()
 {
-  double bytes = (double)nmax * sizeof(double);
+  double bytes = (double) nmax * sizeof(double);
   return bytes;
 }