git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@1119 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/compute.cpp

@@ -18,10 +18,13 @@
 #include "group.h"
 #include "domain.h"
 #include "lattice.h"
+#include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
+#define DELTA 4
+
 /* ---------------------------------------------------------------------- */
 
 Compute::Compute(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
@@ -48,7 +51,9 @@ Compute::Compute(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   vector_atom = NULL;
   
   scalar_flag = vector_flag = peratom_flag = 0;
-  tempflag = pressflag = 0;
+  tempflag = pressflag = peflag = 0;
+  timeflag = 0;
+  invoked = 0;
   npre = 0;
   id_pre = NULL;
   comm_forward = comm_reverse = 0;
@@ -57,6 +62,11 @@ Compute::Compute(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
 
   extra_dof = 3;
   dynamic = 0;
+
+  // setup list of timesteps
+
+  ntime = maxtime = 0;
+  tlist = NULL;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -68,6 +78,8 @@ Compute::~Compute()
 
   for (int i = 0; i < npre; i++) delete [] id_pre[i];
   delete [] id_pre;
+
+  memory->sfree(tlist);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -76,7 +88,8 @@ void Compute::modify_params(int narg, char **arg)
 {
   if (narg == 0) error->all("Illegal compute_modify command");
 
-  int iarg = 0;  while (iarg < narg) {
+  int iarg = 0;
+  while (iarg < narg) {
     if (strcmp(arg[iarg],"extra") == 0) {
       if (iarg+2 > narg) error->all("Illegal compute_modify command");
       extra_dof = atoi(arg[iarg+1]);
@@ -87,6 +100,60 @@ void Compute::modify_params(int narg, char **arg)
       else if (strcmp(arg[iarg+1],"yes") == 0) dynamic = 1;
       else error->all("Illegal compute_modify command");
       iarg += 2;
+    } else if (strcmp(arg[iarg],"thermo") == 0) {
+      if (iarg+2 > narg) error->all("Illegal compute_modify command");
+      if (strcmp(arg[iarg+1],"no") == 0) thermoflag = 0;
+      else if (strcmp(arg[iarg+1],"yes") == 0) thermoflag = 1;
+      else error->all("Illegal compute_modify command");
+      iarg += 2;
     } else error->all("Illegal compute_modify command");
   }
 }
+
+/* ----------------------------------------------------------------------
+   add ntimestep to list of timesteps the compute will be called on
+   do not add if already in list
+   search from top downward, since list of times is in decreasing order
+------------------------------------------------------------------------- */
+
+void Compute::add_step(int ntimestep)
+{
+  // i = location in list to insert ntimestep
+
+  int i;
+  for (i = ntime-1; i >= 0; i--) {
+    if (ntimestep == tlist[i]) return;
+    if (ntimestep < tlist[i]) break;
+  }
+  i++;
+
+  // extend list as needed
+
+  if (ntime == maxtime) {
+    maxtime += DELTA;
+    tlist = (int *)
+      memory->srealloc(tlist,maxtime*sizeof(int),"compute:tlist");
+  }
+
+  // move remainder of list upward and insert ntimestep
+
+  for (int j = ntime-1; j >= i; j--) tlist[j+1] = tlist[j];
+  tlist[i] = ntimestep;
+  ntime++;
+}
+
+/* ----------------------------------------------------------------------
+   return 1/0 if ntimestep is or is not in list of calling timesteps
+   if value(s) on top of list are less than ntimestep, delete them
+   search from top downward, since list of times is in decreasing order
+------------------------------------------------------------------------- */
+
+int Compute::match_step(int ntimestep)
+{
+  for (int i = ntime-1; i >= 0; i--) {
+    if (ntimestep < tlist[i]) return 0;
+    if (ntimestep == tlist[i]) return 1;
+    if (ntimestep > tlist[i]) ntime--;
+  }
+  return 0;
+}

src/compute.h

@@ -39,6 +39,14 @@ class Compute : protected Pointers {
                       // must have both compute_scalar, compute_vector
   int pressflag;      // 1 if Compute can be used as pressure (uses virial)
                       // must have both compute_scalar, compute_vector
+  int peflag;         // 1 if Compute calculates PE (uses Force energies)
+
+  int timeflag;       // 1 if Compute stores list of timesteps it's called on
+  int ntime;          // # of entries in time list
+  int maxtime;        // max # of entries time list can hold
+  int *tlist;         // time list of steps the Compute is called on
+  int invoked;        // 1 if Compute was invoked (e.g. by a variable)
+
   double dof;         // degrees-of-freedom for temperature
 
   int npre;           // # of computes to compute before this one
@@ -61,10 +69,15 @@ class Compute : protected Pointers {
   virtual int pack_reverse_comm(int, int, double *) {return 0;}
   virtual void unpack_reverse_comm(int, int *, double *) {}
 
+  void add_step(int);
+  int match_step(int);
+
   virtual double memory_usage() {return 0.0;}
 
  protected:
-  int extra_dof,dynamic;
+  int extra_dof;               // extra DOF for temperature computes
+  int dynamic;                 // recount atoms for temperature computes
+  int thermoflag;              // 1 if include fix PE for PE computes
 };
 
 }

src/compute_pe.cpp

@@ -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.
+------------------------------------------------------------------------- */
+
+#include "mpi.h"
+#include "string.h"
+#include "compute_pe.h"
+#include "atom.h"
+#include "force.h"
+#include "pair.h"
+#include "bond.h"
+#include "angle.h"
+#include "dihedral.h"
+#include "improper.h"
+#include "kspace.h"
+#include "modify.h"
+#include "domain.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputePE::ComputePE(LAMMPS *lmp, int narg, char **arg) : 
+  Compute(lmp, narg, arg)
+{
+  if (narg < 3) error->all("Illegal compute pe command");
+  if (igroup) error->all("Compute pe must use group all");
+
+  if (narg == 3) {
+    pairflag = 1;
+    bondflag = angleflag = dihedralflag = improperflag = 1;
+    kspaceflag = 1;
+    thermoflag = 1;
+  } else {
+    pairflag = 0;
+    bondflag = angleflag = dihedralflag = improperflag = 0;
+    kspaceflag = 0;
+    thermoflag = 0;
+    int iarg = 3;
+    while (iarg < narg) {
+      if (strcmp(arg[iarg],"pair") == 0) pairflag = 1;
+      else if (strcmp(arg[iarg],"bond") == 0) bondflag = 1;
+      else if (strcmp(arg[iarg],"angle") == 0) angleflag = 1;
+      else if (strcmp(arg[iarg],"dihedral") == 0) dihedralflag = 1;
+      else if (strcmp(arg[iarg],"improper") == 0) improperflag = 1;
+      else if (strcmp(arg[iarg],"kspace") == 0) kspaceflag = 1;
+      else error->all("Illegal compute pe command");
+      iarg++;
+    }
+  }
+
+  // settings
+
+  scalar_flag = 1;
+  extensive = 1;
+  peflag = 1;
+  timeflag = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double ComputePE::compute_scalar()
+{
+  invoked = 1;
+  double one = 0.0;
+
+  if (pairflag) {
+    if (force->pair) one += force->pair->eng_vdwl + force->pair->eng_coul;
+    if (force->bond) one += force->bond->eng_vdwl;
+    if (force->dihedral) 
+      one += force->dihedral->eng_vdwl + force->dihedral->eng_coul;
+  }
+
+  if (atom->molecular) {
+    if (bondflag && force->bond) one += force->bond->energy;
+    if (angleflag && force->angle) one += force->angle->energy;
+    if (dihedralflag && force->dihedral) one += force->dihedral->energy;
+    if (improperflag && force->improper) one += force->improper->energy;
+  }
+
+  MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
+
+  if (kspaceflag && force->kspace) scalar += force->kspace->energy;
+  if (pairflag && force->pair && force->pair->tail_flag) {
+    double volume = domain->xprd * domain->yprd * domain->zprd;
+    scalar += force->pair->etail / volume;
+  }
+
+  if (thermoflag && modify->n_thermo_energy) scalar += modify->thermo_energy();
+
+  return scalar;
+}

src/compute_pe.h

@@ -0,0 +1,34 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#ifndef COMPUTE_PE_H
+#define COMPUTE_PE_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputePE : public Compute {
+ public:
+  ComputePE(class LAMMPS *, int, char **);
+  ~ComputePE() {}
+  void init() {}
+  double compute_scalar();
+
+ private:
+  int pairflag,bondflag,angleflag,dihedralflag,improperflag,kspaceflag;
+};
+
+}
+
+#endif

src/compute_pressure.cpp

@@ -36,7 +36,6 @@ ComputePressure::ComputePressure(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg)
 {
   if (narg != 4) error->all("Illegal compute pressure command");
-
   if (igroup) error->all("Compute pressure must use group all");
 
   // store temperature ID used by pressure computation
@@ -54,10 +53,13 @@ ComputePressure::ComputePressure(LAMMPS *lmp, int narg, char **arg) :
   if (modify->compute[icompute]->tempflag == 0)
     error->all("Compute pressure temp ID does not compute temperature");
 
+  // settings
+
   scalar_flag = vector_flag = 1;
   size_vector = 6;
   extensive = 0;
   pressflag = 1;
+  timeflag = 1;
 
   vector = new double[6];
   nvirial = 0;
@@ -175,6 +177,7 @@ void ComputePressure::virial_compute(int n)
   int i,j;
   double v[6],*vcomponent;
 
+  invoked = 1;
   for (i = 0; i < n; i++) v[i] = 0.0;
 
   // sum contributions to virial from forces and fixes

src/fix.cpp

@@ -38,7 +38,6 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   force_reneighbor = 0;
   box_change = 0;
   thermo_energy = 0;
-  pressure_every = 0;
   rigid_flag = 0;
   virial_flag = 0;
   no_change_box = 0;

src/fix.h

@@ -30,7 +30,6 @@ class Fix : protected Pointers {
   int next_reneighbor;           // next timestep to force a reneighboring
   int thermo_energy;             // 1 if ThEng enabled via fix_modify, 0 if not
   int nevery;                    // how often to call an end_of_step fix
-  int pressure_every;            // how often fix needs virial computed
   int rigid_flag;                // 1 if Fix integrates rigid bodies, 0 if not
   int virial_flag;               // 1 if Fix contributes to virial, 0 if not
   int no_change_box;             // 1 if cannot swap ortho <-> triclinic

src/fix_ave_atom.cpp

@@ -51,8 +51,6 @@ FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) :
   if (modify->compute[icompute]->peratom_flag == 0)
     error->all("Fix ave/atom compute does not calculate per-atom info");
 
-  if (modify->compute[icompute]->pressflag) pressure_every = nevery;
-
   peratom_flag = 1;
 
   // setup list of computes to call, including pre-computes

src/fix_ave_time.cpp

@@ -128,9 +128,6 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
       error->all("Fix ave/time fix does not calculate a vector");
   }
 
-  if (which == COMPUTE &&
-      modify->compute[icompute]->pressflag) pressure_every = nevery;
-
   // setup list of computes to call, including pre-computes
 
   compute = NULL;
@@ -293,6 +290,8 @@ void FixAveTime::end_of_step()
   // accumulate results of compute or fix to local copy
   
   if (which == COMPUTE) {
+    modify->clearstep_compute();
+
     if (sflag) {
       double value;
       for (i = 0; i < ncompute; i++) value = compute[i]->compute_scalar();
@@ -312,10 +311,18 @@ void FixAveTime::end_of_step()
   }
 
   // done if irepeat < nrepeat
+  // else reset irepeat and nvalid
 
   irepeat++;
-  nvalid += nevery;
+  if (irepeat < nrepeat) {
-  if (irepeat < nrepeat) return;
+    nvalid += nevery;
+    if (which == COMPUTE) modify->addstep_compute(nvalid);
+    return;
+  }
+
+  irepeat = 0;
+  nvalid = update->ntimestep+nfreq - (nrepeat-1)*nevery;
+  if (which == COMPUTE) modify->addstep_compute(nvalid);
 
   // average the final result for the Nfreq timestep
 
@@ -323,11 +330,6 @@ void FixAveTime::end_of_step()
   if (sflag) scalar /= repeat;
   if (vflag) for (i = 0; i < size_vector; i++) vector[i] /= repeat;
 
-  // reset irepeat and nvalid
-
-  irepeat = 0;
-  nvalid = update->ntimestep+nfreq - (nrepeat-1)*nevery;
-
   // if ave = ONE, only single Nfreq timestep value is needed
   // if ave = RUNNING, combine with all previous Nfreq timestep values
   // if ave = WINDOW, comine with nwindow most recent Nfreq timestep values

src/fix_nph.cpp

@@ -45,7 +45,6 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) :
   if (narg < 4) error->all("Illegal fix nph command");
 
   restart_global = 1;
-  pressure_every = 1;
   box_change = 1;
   scalar_flag = 1;
   scalar_vector_freq = 1;
@@ -332,6 +331,10 @@ void FixNPH::setup()
     pressure->compute_vector();
   }
   couple();
+
+  // trigger virial computation on next timestep
+
+  pressure->add_step(update->ntimestep+1);
 }
 
 /* ----------------------------------------------------------------------
@@ -444,6 +447,10 @@ void FixNPH::final_integrate()
   }
   couple();
 
+  // trigger virial computation on next timestep
+
+  pressure->add_step(update->ntimestep+1);
+
   // update omega_dot
   // for non-varying dims, p_freq is 0.0, so omega_dot doesn't change
 

src/fix_npt.cpp

@@ -44,7 +44,6 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) :
   if (narg < 7) error->all("Illegal fix npt command");
 
   restart_global = 1;
-  pressure_every = 1;
   box_change = 1;
   scalar_flag = 1;
   scalar_vector_freq = 1;
@@ -334,6 +333,10 @@ void FixNPT::setup()
     pressure->compute_vector();
   }
   couple();
+
+  // trigger virial computation on next timestep
+
+  pressure->add_step(update->ntimestep+1);
 }
 
 /* ----------------------------------------------------------------------
@@ -454,6 +457,10 @@ void FixNPT::final_integrate()
   }
   couple();
 
+  // trigger virial computation on next timestep
+
+  pressure->add_step(update->ntimestep+1);
+
   // update eta_dot
 
   f_eta = t_freq*t_freq * (t_current/t_target - 1.0);

src/fix_print.cpp

@@ -16,6 +16,7 @@
 #include "fix_print.h"
 #include "update.h"
 #include "input.h"
+#include "modify.h"
 #include "variable.h"
 #include "error.h"
 
@@ -67,11 +68,16 @@ void FixPrint::end_of_step()
   // make a copy of line to work on
   // substitute for $ variables (no printing)
   // append a newline and print final copy
+  // variable evaluation may invoke a compute that affects Verlet::eflag,vflag
+
+  modify->clearstep_compute();
 
   strcpy(copy,line);
   input->substitute(copy,0);
   strcat(copy,"\n");
 
+  modify->addstep_compute(update->ntimestep + nevery);
+
   if (me == 0) {
     if (screen) fprintf(screen,copy);
     if (logfile) fprintf(logfile,copy);

src/integrate.cpp

@@ -0,0 +1,53 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#include "stdlib.h"
+#include "integrate.h"
+#include "compute.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+Integrate::Integrate(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
+{
+  elist = vlist = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+Integrate::~Integrate()
+{
+  delete [] elist;
+  delete [] vlist;
+}
+
+/* ----------------------------------------------------------------------
+   set eflag if a pe compute is called this timestep
+   set vflag if a pressure compute is called this timestep
+------------------------------------------------------------------------- */
+
+void Integrate::ev_set(int ntimestep)
+{
+  int i;
+
+  eflag = 0;
+  for (i = 0; i < nelist; i++)
+    if (elist[i]->match_step(ntimestep)) break;
+  if (i < nelist) eflag = 1;
+
+  vflag = 0;
+  for (i = 0; i < nvlist; i++)
+    if (vlist[i]->match_step(ntimestep)) break;
+  if (i < nvlist) vflag = virial_style;
+}

src/integrate.h

@@ -20,14 +20,24 @@ namespace LAMMPS_NS {
 
 class Integrate : protected Pointers {
  public:
-  Integrate(class LAMMPS *lmp, int, char **) : Pointers(lmp) {}
+  Integrate(class LAMMPS *, int, char **);
-  virtual ~Integrate() {}
+  virtual ~Integrate();
   virtual void init() = 0;
   virtual void setup() = 0;
   virtual void iterate(int) = 0;
   virtual void cleanup() {}
-  virtual double memory_usage() {return 0.0;}
   virtual void reset_dt() {}
+  virtual double memory_usage() {return 0.0;}
+
+ protected:
+  int eflag,vflag;                  // flags for energy/virial computation
+  int virial_style;                 // compute virial explicitly or implicitly
+
+  int nelist,nvlist;                // # of PE,virial coputes for eflag,vflag
+  class Compute **elist;            // list of Computes to check
+  class Compute **vlist;
+
+  void ev_set(int);
 };
 
 }

src/min_cg.cpp

@@ -36,7 +36,9 @@
 #include "thermo.h"
 #include "update.h"
 #include "modify.h"
+#include "compute.h"
 #include "fix_minimize.h"
+#include "thermo.h"
 #include "timer.h"
 #include "memory.h"
 #include "error.h"
@@ -55,7 +57,17 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-MinCG::MinCG(LAMMPS *lmp) : Min(lmp) {}
+MinCG::MinCG(LAMMPS *lmp) : Min(lmp)
+{
+  vlist = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+MinCG::~MinCG()
+{
+  delete [] vlist;
+}
 
 /* ---------------------------------------------------------------------- */
 
@@ -77,11 +89,27 @@ void MinCG::init()
   setup_vectors();
   for (int i = 0; i < ndof; i++) h[i] = g[i] = 0.0;
 
-  // virial_thermo is how virial should be computed on thermo timesteps
+  // virial_style:
-  // 1 = computed explicity by pair, 2 = computed implicitly by pair
+  // 1 if computed explicitly by pair->compute via sum over pair interactions
+  // 2 if computed implicitly by pair->virial_compute via sum over ghost atoms
 
-  if (force->newton_pair) virial_thermo = 2;
+  if (force->newton_pair) virial_style = 2;
-  else virial_thermo = 1;
+  else virial_style = 1;
+
+  // vlist = list of computes for pressure
+
+  delete [] vlist;
+  vlist = NULL;
+
+  nvlist = 0;
+  for (int i = 0; i < modify->ncompute; i++)
+    if (modify->compute[i]->pressflag) nvlist++;
+
+  if (nvlist) vlist = new Compute*[nvlist];
+
+  nvlist = 0;
+  for (int i = 0; i < modify->ncompute; i++)
+    if (modify->compute[i]->pressflag) vlist[nvlist++] = modify->compute[i];
 
   // set flags for what arrays to clear in force_clear()
   // need to clear torques if array exists
@@ -123,11 +151,17 @@ void MinCG::run()
   double tmp,*f;
 
   // set initial force & energy
+  // normalize energy if thermo PE does
 
   setup();
   setup_vectors();
-  output->thermo->compute_pe();
+
-  ecurrent = output->thermo->potential_energy;
+  int id = modify->find_compute("thermo_pe");
+  if (id < 0) error->all("Minimization could not find thermo_pe compute");
+  pe_compute = modify->compute[id];
+
+  ecurrent = pe_compute->compute_scalar();
+  if (output->thermo->normflag) ecurrent /= atom->natoms;
 
   // stats for Finish to print
 	
@@ -221,8 +255,7 @@ void MinCG::setup()
 
   // compute all forces
 
-  int eflag = 1;
+  ev_set(update->ntimestep);
-  int vflag = virial_thermo;
   force_clear(vflag);
   
   if (force->pair) force->pair->compute(eflag,vflag);
@@ -252,7 +285,7 @@ void MinCG::setup()
 
 void MinCG::iterate(int n)
 {
-  int i,gradsearch,fail;
+  int i,gradsearch,fail,ntimestep;
   double alpha,beta,gg,dot[2],dotall[2];
   double *f;
 
@@ -268,7 +301,7 @@ void MinCG::iterate(int n)
 
   for (niter = 0; niter < n; niter++) {
 
-    update->ntimestep++;
+    ntimestep = ++update->ntimestep;
 
     // line minimization along direction h from current atom->x
 
@@ -317,9 +350,9 @@ void MinCG::iterate(int n)
 
     // output for thermo, dump, restart files
 
-    if (output->next == update->ntimestep) {
+    if (output->next == ntimestep) {
       timer->stamp();
-      output->write(update->ntimestep);
+      output->write(ntimestep);
       timer->stamp(TIME_OUTPUT);
     }
   }
@@ -375,11 +408,7 @@ void MinCG::eng_force(int *pndof, double **px, double **ph, double *peng)
     setup_vectors();
   }
 
-  // eflag is always set, since minimizer needs potential energy
+  ev_set(update->ntimestep);
-
-  int eflag = 1;
-  int vflag = 0;
-  if (output->next_thermo == update->ntimestep) vflag = virial_thermo;
   force_clear(vflag);
 
   timer->stamp();
@@ -411,10 +440,11 @@ void MinCG::eng_force(int *pndof, double **px, double **ph, double *peng)
 
   if (modify->n_min_post_force) modify->min_post_force(vflag);
 
-  // compute potential energy of system via Thermo
+  // compute potential energy of system
+  // normalize if thermo PE does
 
-  output->thermo->compute_pe();
+  ecurrent = pe_compute->compute_scalar();
-  ecurrent = output->thermo->potential_energy;
+  if (output->thermo->normflag) ecurrent /= atom->natoms;
 
   // return updated ptrs to caller since atoms may have migrated
 
@@ -424,6 +454,22 @@ void MinCG::eng_force(int *pndof, double **px, double **ph, double *peng)
   *peng = ecurrent;
 }
 
+/* ----------------------------------------------------------------------
+   eflag is always set, since minimizer needs potential energy
+   set vflag if a pressure compute is called this timestep
+------------------------------------------------------------------------- */
+
+void MinCG::ev_set(int ntimestep)
+{
+  int i;
+
+  eflag = 1;
+  vflag = 0;
+  for (i = 0; i < nvlist; i++)
+    if (vlist[i]->match_step(ntimestep)) break;
+  if (i < nvlist) vflag = virial_style;
+}
+
 /* ----------------------------------------------------------------------
    clear force on own & ghost atoms
    setup and clear other arrays as needed

src/min_cg.h

@@ -21,21 +21,26 @@ namespace LAMMPS_NS {
 class MinCG : public Min {
  public:
   MinCG(class LAMMPS *);
+  virtual ~MinCG();
   void init();
   void run();
-
   virtual void iterate(int);
 
  protected:
-  int virial_thermo;          // what vflag should be on thermo steps (1,2)
+  int eflag,vflag;            // flags for energy/virial computation
+  int virial_style;           // compute virial explicitly or implicitly
   int pairflag,torqueflag;
   int neigh_every,neigh_delay,neigh_dist_check;   // copies of reneigh criteria
   int triclinic;              // 0 if domain is orthog, 1 if triclinic
 
+  int nvlist;                 // # of PE,virial coputes for eflag,vflag
+  class Compute **vlist;      // list of Computes to check
+
   int maxpair;                // copies of Update quantities
   double **f_pair;
 
   class FixMinimize *fix_minimize;  // fix that stores gradient vecs
+  class Compute *pe_compute;  // compute for potential energy
   double ecurrent;            // current potential energy
   double mindist,maxdist;     // min/max dist for coord delta in line search
 
@@ -54,6 +59,7 @@ class MinCG : public Min {
   void setup();
   void setup_vectors();
   void eng_force(int *, double **, double **, double *);
+  void ev_set(int);
   void force_clear(int);
 };
 

src/modify.cpp

@@ -154,8 +154,12 @@ void Modify::init()
   list_init(MIN_ENERGY,n_min_energy,list_min_energy);
 
   // init each compute
+  // notify relevant computes they may be called on this timestep
 
-  for (i = 0; i < ncompute; i++) compute[i]->init();
+  for (i = 0; i < ncompute; i++) {
+    compute[i]->init();
+    if (compute[i]->timeflag) compute[i]->add_step(update->ntimestep);
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -568,6 +572,28 @@ int Modify::find_compute(char *id)
   return icompute;
 }
 
+/* ----------------------------------------------------------------------
+   clear the invoked flag of computes that stores their next invocation
+   called by classes that are invoking computes
+------------------------------------------------------------------------- */
+
+void Modify::clearstep_compute()
+{
+  for (int icompute = 0; icompute < ncompute; icompute++)
+    if (compute[icompute]->timeflag) compute[icompute]->invoked = 0;
+}
+
+/* ----------------------------------------------------------------------
+   schedule the next invocation of computes that were invoked
+   called by classes that invoked computes to schedule the next invocation
+------------------------------------------------------------------------- */
+
+void Modify::addstep_compute(int ntimestep)
+{
+  for (int icompute = 0; icompute < ncompute; icompute++)
+    if (compute[icompute]->invoked) compute[icompute]->add_step(ntimestep);
+}
+
 /* ----------------------------------------------------------------------
    write to restart file for all Fixes with restart info
    (1) fixes that have global state

src/modify.h

@@ -65,6 +65,8 @@ class Modify : protected Pointers {
   void modify_compute(int, char **);
   void delete_compute(char *);
   int find_compute(char *);
+  void clearstep_compute();
+  void addstep_compute(int);
 
   void write_restart(FILE *);
   int read_restart(FILE *);

src/output.cpp

@@ -46,18 +46,25 @@ using namespace LAMMPS_NS;
 
 Output::Output(LAMMPS *lmp) : Pointers(lmp)
 {
-  // create 2 default computes for temp and pressure
+  // create default computes for temp,pressure,pe
 
   char **newarg = new char*[4];
   newarg[0] = (char *) "thermo_temp";
   newarg[1] = (char *) "all";
   newarg[2] = (char *) "temp";
   modify->add_compute(3,newarg);
+
   newarg[0] = (char *) "thermo_pressure";
   newarg[1] = (char *) "all";
   newarg[2] = (char *) "pressure";
   newarg[3] = (char *) "thermo_temp";
   modify->add_compute(4,newarg);
+
+  newarg[0] = (char *) "thermo_pe";
+  newarg[1] = (char *) "all";
+  newarg[2] = (char *) "pe";
+  modify->add_compute(3,newarg);
+
   delete [] newarg;
 
   // create default Thermo class
@@ -157,6 +164,8 @@ void Output::setup(int flag)
   // set next_thermo to multiple of every or last step of run (if smaller)
   // if every = 0, set next_thermo to last step of run
 
+  modify->clearstep_compute();
+
   thermo->header();
   thermo->compute(0);
   last_thermo = ntimestep;
@@ -166,6 +175,8 @@ void Output::setup(int flag)
     next_thermo = MYMIN(next_thermo,update->laststep);
   } else next_thermo = update->laststep;
 
+  modify->addstep_compute(next_thermo);
+
   // next = next timestep any output will be done
 
   next = MYMIN(next_dump_any,next_restart);
@@ -220,10 +231,12 @@ void Output::write(int ntimestep)
   // insure next_thermo forces output on last step of run
 
   if (next_thermo == ntimestep && last_thermo != ntimestep) {
+    modify->clearstep_compute();
     thermo->compute(1);
     last_thermo = ntimestep;
     next_thermo += thermo_every;
     next_thermo = MYMIN(next_thermo,update->laststep);
+    modify->addstep_compute(next_thermo);
   }
 
   // next = next timestep any output will be done

src/respa.cpp

@@ -33,6 +33,7 @@
 #include "output.h"
 #include "update.h"
 #include "modify.h"
+#include "compute.h"
 #include "fix_respa.h"
 #include "memory.h"
 #include "timer.h"
@@ -40,9 +41,6 @@
 
 using namespace LAMMPS_NS;
 
-#define MIN(A,B) ((A) < (B)) ? (A) : (B)
-#define MAX(A,B) ((A) > (B)) ? (A) : (B)
-
 /* ---------------------------------------------------------------------- */
 
 Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg)
@@ -225,7 +223,7 @@ Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg)
     cutoff[3] = cutoff[1];
   }
 
-  // allocate other needed level arrays
+  // allocate other needed arrays
 
   newton = new int[nlevels];
   step = new double[nlevels];
@@ -271,28 +269,29 @@ void Respa::init()
     if (force->pair && force->pair->respa_enable == 0)
       error->all("Pair style does not support rRESPA inner/middle/outer");
 
-  // setup virial computations for timestepping
   // virial_style = 1 (explicit) since never computed implicitly like Verlet
-  // virial_every = 1 if computed every timestep (NPT,NPH)
-  // fix arrays store info on fixes that need virial computed occasionally
 
   virial_style = 1;
 
-  virial_every = 0;
+  // elist,vlist = list of computes for PE and pressure
-  nfix_virial = 0;
-  for (int i = 0; i < modify->nfix; i++)
-    if (modify->fix[i]->pressure_every == 1) virial_every = 1;
-    else if (modify->fix[i]->pressure_every > 1) nfix_virial++;
 
-  if (nfix_virial) {
+  delete [] elist;
-    delete [] fix_virial_every;
+  delete [] vlist;
-    delete [] next_fix_virial;
+  elist = vlist = NULL;
-    fix_virial_every = new int[nfix_virial];
+
-    next_fix_virial = new int[nfix_virial];
+  nelist = nvlist = 0;
-    nfix_virial = 0;
+  for (int i = 0; i < modify->ncompute; i++) {
-    for (int i = 0; i < modify->nfix; i++)
+    if (modify->compute[i]->peflag) nelist++;
-      if (modify->fix[i]->pressure_every > 1)
+    if (modify->compute[i]->pressflag) nvlist++;
-	fix_virial_every[nfix_virial++] = modify->fix[i]->pressure_every;
+  }
+
+  if (nelist) elist = new Compute*[nelist];
+  if (nvlist) vlist = new Compute*[nvlist];
+
+  nelist = nvlist = 0;
+  for (int i = 0; i < modify->ncompute; i++) {
+    if (modify->compute[i]->peflag) elist[nelist++] = modify->compute[i];
+    if (modify->compute[i]->pressflag) vlist[nvlist++] = modify->compute[i];
   }
 
   // step[] = timestep for each level
@@ -347,8 +346,7 @@ void Respa::setup()
 
   // compute all forces
 
-  eflag = 1;
+  ev_set(update->ntimestep);
-  vflag = virial_style;
 
   for (int ilevel = 0; ilevel < nlevels; ilevel++) {
     force_clear(newton[ilevel]);
@@ -379,21 +377,6 @@ void Respa::setup()
   modify->setup();
   sum_flevel_f();
   output->setup(1);
-
-  // setup virial computations for timestepping
-
-  int ntimestep = update->ntimestep;
-  next_virial = 0;
-  if (virial_every) next_virial = ntimestep + 1;
-  else {
-    for (int ivirial = 0; ivirial < nfix_virial; ivirial++) {
-      next_fix_virial[ivirial] = 
-	(ntimestep/fix_virial_every[ivirial])*fix_virial_every[ivirial] + 
-	fix_virial_every[ivirial];
-      if (ivirial) next_virial = MIN(next_virial,next_fix_virial[ivirial]);
-      else next_virial = next_fix_virial[0];
-    }
-  }
 }
 
 /* ----------------------------------------------------------------------
@@ -408,16 +391,7 @@ void Respa::iterate(int n)
 
     ntimestep = ++update->ntimestep;
 
-    // eflag/vflag = 0/1 for energy/virial computation
+    ev_set(ntimestep);
-
-    if (ntimestep == output->next_thermo) eflag = 1;
-    else eflag = 0;
-
-    if (ntimestep == output->next_thermo || ntimestep == next_virial) {
-      vflag = virial_style;
-      if (virial_every) next_virial++;
-      else next_virial = fix_virial(ntimestep);
-    } else vflag = 0;
 
     recurse(nlevels-1);
 
@@ -541,7 +515,8 @@ void Respa::recurse(int ilevel)
       timer->stamp(TIME_COMM);
     }
   
-    if (modify->n_post_force_respa) modify->post_force_respa(vflag,ilevel,iloop);
+    if (modify->n_post_force_respa)
+      modify->post_force_respa(vflag,ilevel,iloop);
     modify->final_integrate_respa(ilevel);
   }
 
@@ -623,20 +598,3 @@ void Respa::sum_flevel_f()
     }
   }
 }
-
-/* ----------------------------------------------------------------------
-   return next timestep virial should be computed
-   based on one or more fixes that need virial computed periodically
-------------------------------------------------------------------------- */
-
-int Respa::fix_virial(int ntimestep)
-{
-  int next;
-  for (int ivirial = 0; ivirial < nfix_virial; ivirial++) {
-    if (ntimestep == next_fix_virial[ivirial])
-      next_fix_virial[ivirial] += fix_virial_every[ivirial];
-    if (ivirial) next = MIN(next,next_fix_virial[ivirial]);
-    else next = next_fix_virial[0];
-  }
-  return next;
-}

src/respa.h

@@ -45,13 +45,6 @@ class Respa : public Integrate {
   void copy_flevel_f(int);
 
  private:
-  int eflag,vflag;                  // flags for energy/virial computation
-  int virial_style;                 // compute virial explicitly (not implicit)
-  int virial_every;                 // 1 if virial computed every step
-  int next_virial;                  // next timestep to compute virial
-  int nfix_virial;                  // # of fixes that need virial occasionally
-  int *fix_virial_every;            // frequency they require it
-  int *next_fix_virial;             // next timestep they need it
   int triclinic;                    // 0 if domain is orthog, 1 if triclinic
 
   int *newton;                      // newton flag at each level
@@ -60,7 +53,6 @@ class Respa : public Integrate {
   void recurse(int);
   void force_clear(int);
   void sum_flevel_f();
-  int fix_virial(int);
 };
 
 }

src/run.cpp

@@ -17,6 +17,7 @@
 #include "domain.h"
 #include "update.h"
 #include "integrate.h"
+#include "modify.h"
 #include "output.h"
 #include "finish.h"
 #include "input.h"
@@ -182,7 +183,13 @@ void Run::command(int narg, char **arg)
       if (postflag || nleft == 0) finish.end(1);
       else finish.end(0);
 
-      if (commandstr) char *command = input->one(commandstr);
+      // command string may invoke a compute that affects Verlet::eflag,vflag
+
+      if (commandstr) {
+	modify->clearstep_compute();
+	char *command = input->one(commandstr);
+	modify->addstep_compute(update->ntimestep + nevery);
+      }
 
       nleft -= nsteps;
       iter++;

src/style.h

@@ -82,6 +82,7 @@ CommandStyle(write_restart,WriteRestart)
 #include "compute_ebond_atom.h"
 #include "compute_epair_atom.h"
 #include "compute_ke_atom.h"
+#include "compute_pe.h"
 #include "compute_pressure.h"
 #include "compute_rotate_dipole.h"
 #include "compute_rotate_gran.h"
@@ -104,6 +105,7 @@ ComputeStyle(coord/atom,ComputeCoordAtom)
 ComputeStyle(ebond/atom,ComputeEbondAtom)
 ComputeStyle(epair/atom,ComputeEpairAtom)
 ComputeStyle(ke/atom,ComputeKEAtom)
+ComputeStyle(pe,ComputePE)
 ComputeStyle(pressure,ComputePressure)
 ComputeStyle(rotate/dipole,ComputeRotateDipole)
 ComputeStyle(rotate/gran,ComputeRotateGran)

src/temper.cpp

@@ -21,9 +21,11 @@
 #include "temper.h"
 #include "universe.h"
 #include "domain.h"
+#include "atom.h"
 #include "update.h"
 #include "integrate.h"
 #include "modify.h"
+#include "compute.h"
 #include "force.h"
 #include "output.h"
 #include "thermo.h"
@@ -93,15 +95,6 @@ void Temper::command(int narg, char **arg)
   if (nswaps*nevery != nsteps) 
     error->universe_all("Non integer # of swaps in temper command");
 
-  // thermodynamics must be computed on swap steps
-  // potential energy must be computed by thermo
-
-  if (nevery % output->thermo_every)
-    error->universe_all("Thermodynamics not computed on tempering swap steps");
-
-  if (output->thermo->peflag == 0)
-    error->universe_all("Thermodynamics must compute PE for temper command");
-
   // fix style must be appropriate for temperature control
 
   if (strcmp(modify->fix[whichfix]->style,"nvt") == 0) fixstyle = NVT;
@@ -126,6 +119,14 @@ void Temper::command(int narg, char **arg)
   iworld = universe->iworld;
   boltz = force->boltz;
 
+  // pe_compute = ptr to thermo_pe compute
+  // notify compute it will be called at first swap
+
+  int id = modify->find_compute("thermo_pe");
+  if (id < 0) error->all("Tempering could not find thermo_pe compute");
+  Compute *pe_compute = modify->compute[id];
+  pe_compute->add_step(update->ntimestep + nevery);
+
   // create MPI communicator for root proc from each world
 
   int color;
@@ -203,6 +204,13 @@ void Temper::command(int narg, char **arg)
 
     update->integrate->iterate(nevery);
 
+    // compute PE, normalize if thermo PE does
+    // notify compute it will be called at next swap
+
+    pe = pe_compute->compute_scalar();
+    if (output->thermo->normflag) pe /= atom->natoms;
+    pe_compute->add_step(update->ntimestep + nevery);
+
     // which = which of 2 kinds of swaps to do (0,1)
 
     if (!ranswap) which = iswap % 2;
@@ -235,7 +243,6 @@ void Temper::command(int narg, char **arg)
 
     swap = 0;
     if (partner != -1) {
-      pe = output->thermo->potential_energy;
       if (me_universe > partner) 
 	MPI_Send(&pe,1,MPI_DOUBLE,partner,0,universe->uworld);
       else

src/thermo.cpp

@@ -111,14 +111,16 @@ Thermo::Thermo(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
 
   temperature = NULL;
   pressure = NULL;
+  pe = NULL;
   rotate_dipole = NULL;
   rotate_gran = NULL;
 
-  index_temp = index_press = index_drot = index_grot = -1;
+  index_temp = index_press = index_pe = index_drot = index_grot = -1;
   internal_drot = internal_grot = 0;
 
   id_temp = (char *) "thermo_temp";
   id_press = (char *) "thermo_pressure";
+  id_pe = (char *) "thermo_pe";
   id_drot = (char *) "thermo_rotate_dipole";
   id_grot = (char *) "thermo_rotate_gran";
 
@@ -131,7 +133,7 @@ Thermo::Thermo(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   parse_fields(line);
 
   // create the requested compute styles
-  // temperature and pressure always exist b/c Output class created them
+  // temperature,pressure,pe always exist b/c Output class created them
 
   if (index_drot >= 0) {
     create_compute(id_drot,(char *) "rotate/dipole",NULL);
@@ -247,6 +249,7 @@ void Thermo::init()
 
   if (index_temp >= 0) temperature = computes[index_temp];
   if (index_press >= 0) pressure = computes[index_press];
+  if (index_pe >= 0) pe = computes[index_pe];
   if (index_drot >= 0) rotate_dipole = computes[index_drot];
   if (index_grot >= 0) rotate_gran = computes[index_grot];
 }
@@ -593,13 +596,12 @@ void Thermo::deallocate()
 
 /* ----------------------------------------------------------------------
    parse list of thermo keywords from str
-   set compute flags (temp, press, etc)
+   set compute flags (temp, press, pe, etc)
 ------------------------------------------------------------------------- */
 
 void Thermo::parse_fields(char *str)
 {
   nfield = 0;
-  peflag = 0;
 
   // customize a new keyword by adding to if statement
 
@@ -622,39 +624,50 @@ void Thermo::parse_fields(char *str)
       index_press = add_compute(id_press,0);
     } else if (strcmp(word,"pe") == 0) {
       addfield("PotEng",&Thermo::compute_pe,FLOAT);
-      peflag = 1;
+      index_pe = add_compute(id_pe,0);
     } else if (strcmp(word,"ke") == 0) {
       addfield("KinEng",&Thermo::compute_ke,FLOAT);
       index_temp = add_compute(id_temp,0);
     } else if (strcmp(word,"etotal") == 0) {
       addfield("TotEng",&Thermo::compute_etotal,FLOAT);
-      peflag = 1;
       index_temp = add_compute(id_temp,0);
+      index_pe = add_compute(id_pe,0);
     } else if (strcmp(word,"enthalpy") == 0) {
       addfield("Enthalpy",&Thermo::compute_enthalpy,FLOAT);
       index_temp = add_compute(id_temp,0);
       index_press = add_compute(id_press,0);
+      index_pe = add_compute(id_pe,0);
 
     } else if (strcmp(word,"evdwl") == 0) {
       addfield("E_vdwl",&Thermo::compute_evdwl,FLOAT);
+      index_pe = add_compute(id_pe,0);
     } else if (strcmp(word,"ecoul") == 0) {
       addfield("E_coul",&Thermo::compute_ecoul,FLOAT);
+      index_pe = add_compute(id_pe,0);
     } else if (strcmp(word,"epair") == 0) {
       addfield("E_pair",&Thermo::compute_epair,FLOAT);
+      index_pe = add_compute(id_pe,0);
     } else if (strcmp(word,"ebond") == 0) {
       addfield("E_bond",&Thermo::compute_ebond,FLOAT);
+      index_pe = add_compute(id_pe,0);
     } else if (strcmp(word,"eangle") == 0) {
       addfield("E_angle",&Thermo::compute_eangle,FLOAT);
+      index_pe = add_compute(id_pe,0);
     } else if (strcmp(word,"edihed") == 0) {
       addfield("E_dihed",&Thermo::compute_edihed,FLOAT);
+      index_pe = add_compute(id_pe,0);
     } else if (strcmp(word,"eimp") == 0) {
       addfield("E_impro",&Thermo::compute_eimp,FLOAT);
+      index_pe = add_compute(id_pe,0);
     } else if (strcmp(word,"emol") == 0) {
       addfield("E_mol",&Thermo::compute_emol,FLOAT);
+      index_pe = add_compute(id_pe,0);
     } else if (strcmp(word,"elong") == 0) {
       addfield("E_long",&Thermo::compute_elong,FLOAT);
+      index_pe = add_compute(id_pe,0);
     } else if (strcmp(word,"etail") == 0) {
       addfield("E_tail",&Thermo::compute_etail,FLOAT);
+      index_pe = add_compute(id_pe,0);
 
     } else if (strcmp(word,"vol") == 0) {
       addfield("Volume",&Thermo::compute_vol,FLOAT);
@@ -876,34 +889,78 @@ void Thermo::create_compute(char *id, char *cstyle, char *extra)
 }
 
 /* ----------------------------------------------------------------------
-   compute a single thermodyanmic value matching word to custom list
+   compute a single thermodyanmic value
-   called when a variable is evaluated in input script
+   word is any supported keyword in custom list
+   called when a variable is evaluated by Variable class
    return value as double in answer
    return 0 if OK, 1 if str is invalid
-   customize a new keyword by adding to if statement
+   customize a new keyword by adding to if statement and error tests
 ------------------------------------------------------------------------- */
 
 int Thermo::evaluate_keyword(char *word, double *answer)
 {
-  // could do tests to insure user does not invoke variable at wrong time
+  // don't allow use of thermo keyword before first run
+  // since system is not setup (e.g. no forces have been called for energy)
 
-  // all keywords:
+  if (update->first_update == 0)
-  //   if (domain->box_exist == 0)
+    error->all("Variable equal thermo keyword used before initial run");
-  //   error->all("Variable equal keyword used before simulation box defined");
-  // keywords except step,atoms,vol,ly,ly,lz
-  //   if (update->first_update == 0)
-  //     error->all("Variable equal keyword used before initial run");
-  //   if (update->ntimestep != output->next_thermo && me == 0)
-  //     error->warning("Variable equal keyword used with non-current thermo");
-  // keywords that require Compute objects like T,P
-  //   ptrs like temperature will not be set if init() hasn't been called
-  //   Compute objects will not exist if thermo style isn't using them
 
-  // toggle thermoflag off/on so inidividual compute routines don't assume
+  // check if Compute pointers exist for keywords that need them
-  //   they are being called from compute() which pre-calls Compute objects
+  // error if thermo style does not use the Compute
+  // all energy-related keywords require PE, even if they don't call it
+  // this is b/c Verlet::eflag is triggered by compute_pe invocation schedule
+
+  if (strcmp(word,"temp") == 0 || strcmp(word,"press") == 0 ||
+      strcmp(word,"ke") == 0 || strcmp(word,"etotal") == 0 ||
+      strcmp(word,"pxx") == 0 || strcmp(word,"pyy") == 0 ||
+      strcmp(word,"pzz") == 0 || strcmp(word,"pxy") == 0 ||
+      strcmp(word,"pxz") == 0 || strcmp(word,"pyz") == 0)
+    if (!temperature)
+      error->all("Variable uses compute via thermo keyword unused by thermo");
+
+  if (strcmp(word,"press") == 0 ||
+      strcmp(word,"pxx") == 0 || strcmp(word,"pyy") == 0 ||
+      strcmp(word,"pzz") == 0 || strcmp(word,"pxy") == 0 ||
+      strcmp(word,"pxz") == 0 || strcmp(word,"pyz") == 0)
+    if (!pressure)
+      error->all("Variable uses compute via thermo keyword unused by thermo");
+
+  if (strcmp(word,"pe") == 0 || strcmp(word,"etotal") == 0 ||
+      strcmp(word,"enthalpy") == 0 || strcmp(word,"evdwl") == 0 ||
+      strcmp(word,"ecoul") == 0 || strcmp(word,"epair") == 0 ||
+      strcmp(word,"ebond") == 0 || strcmp(word,"eangle") == 0 ||
+      strcmp(word,"edihed") == 0 || strcmp(word,"eimp") == 0 ||
+      strcmp(word,"emol") == 0 || strcmp(word,"elong") == 0 ||
+      strcmp(word,"etail") == 0)
+    if (!pe)
+      error->all("Variable uses compute via thermo keyword unused by thermo");
+
+  if (strcmp(word,"drot") == 0)
+    if (!rotate_dipole)
+      error->all("Variable uses compute via thermo keyword unused by thermo");
+
+  if (strcmp(word,"grot") == 0)
+    if (!rotate_gran)
+      error->all("Variable uses compute via thermo keyword unused by thermo");
+
+  // set compute_pe invocation flag for keywords that use energy
+  // but don't call compute_pe explicitly
+
+  if (strcmp(word,"evdwl") == 0 || strcmp(word,"ecoul") == 0 ||
+      strcmp(word,"epair") == 0 || strcmp(word,"ebond") == 0 ||
+      strcmp(word,"eangle") == 0 || strcmp(word,"edihed") == 0 ||
+      strcmp(word,"eimp") == 0 || strcmp(word,"emol") == 0 ||
+      strcmp(word,"elong") == 0 || strcmp(word,"etail") == 0)
+    pe->invoked = 1;
+
+  // toggle thermoflag off/on
+  // so individual compute routines know they are not being called from
+  // Thermo::compute() which pre-calls Computes
 
   thermoflag = 0;
 
+  // invoke the lo-level thermo routine to compute the variable value
+
   if (strcmp(word,"step") == 0) {
     compute_step();
     dvalue = ivalue;
@@ -914,12 +971,13 @@ int Thermo::evaluate_keyword(char *word, double *answer)
   else if (strcmp(word,"cpu") == 0) compute_cpu();
 
   else if (strcmp(word,"temp") == 0) compute_temp();
+
   else if (strcmp(word,"press") == 0) compute_press();
   else if (strcmp(word,"pe") == 0) compute_pe();
   else if (strcmp(word,"ke") == 0) compute_ke();
   else if (strcmp(word,"etotal") == 0) compute_etotal();
   else if (strcmp(word,"enthalpy") == 0) compute_enthalpy();
-
+  
   else if (strcmp(word,"evdwl") == 0) compute_evdwl();
   else if (strcmp(word,"ecoul") == 0) compute_ecoul();
   else if (strcmp(word,"epair") == 0) compute_epair();
@@ -930,7 +988,7 @@ int Thermo::evaluate_keyword(char *word, double *answer)
   else if (strcmp(word,"emol") == 0) compute_emol();
   else if (strcmp(word,"elong") == 0) compute_elong();
   else if (strcmp(word,"etail") == 0) compute_etail();
-
+  
   else if (strcmp(word,"vol") == 0) compute_vol();
   else if (strcmp(word,"lx") == 0) compute_lx();
   else if (strcmp(word,"ly") == 0) compute_ly();
@@ -949,7 +1007,7 @@ int Thermo::evaluate_keyword(char *word, double *answer)
   else if (strcmp(word,"pxy") == 0) compute_pxy();
   else if (strcmp(word,"pxz") == 0) compute_pxz();
   else if (strcmp(word,"pyz") == 0) compute_pyz();
-
+  
   else if (strcmp(word,"drot") == 0) compute_drot();
   else if (strcmp(word,"grot") == 0) compute_grot();
 
@@ -1049,33 +1107,9 @@ void Thermo::compute_press()
 
 void Thermo::compute_pe()
 {
-  double tmp = 0.0;
+  if (thermoflag) dvalue = pe->scalar;
-  if (force->pair) tmp += force->pair->eng_vdwl + force->pair->eng_coul;
+  else dvalue = pe->compute_scalar();
-  if (force->bond) tmp += force->bond->eng_vdwl;
-  if (force->dihedral) 
-    tmp += force->dihedral->eng_vdwl + force->dihedral->eng_coul;
-
-  if (atom->molecular) {
-    if (force->bond) tmp += force->bond->energy;
-    if (force->angle) tmp += force->angle->energy;
-    if (force->dihedral) tmp += force->dihedral->energy;
-    if (force->improper) tmp += force->improper->energy;
-  }
-
-  MPI_Allreduce(&tmp,&dvalue,1,MPI_DOUBLE,MPI_SUM,world);
-
-  if (force->kspace) dvalue += force->kspace->energy;
-  if (force->pair && force->pair->tail_flag) {
-    double volume = domain->xprd * domain->yprd * domain->zprd;
-    dvalue += force->pair->etail / volume;
-  }
-  if (modify->n_thermo_energy) dvalue += modify->thermo_energy();
-
   if (normflag) dvalue /= natoms;
-
-  // also store PE in potential_energy so other classes can grab it
-
-  potential_energy = dvalue;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1093,6 +1127,7 @@ void Thermo::compute_ke()
 void Thermo::compute_etotal()
 {
   compute_pe();
+
   double ke;
   if (thermoflag) ke = temperature->scalar;
   else ke = temperature->compute_scalar();

src/thermo.h

@@ -24,9 +24,7 @@ class Thermo : protected Pointers {
 
  public:
   char *style;
-  int peflag;
   int normflag;          // 0 if do not normalize by atoms, 1 if normalize
-  double potential_energy;
 
   Thermo(class LAMMPS *, int, char **);
   ~Thermo();
@@ -71,10 +69,10 @@ class Thermo : protected Pointers {
                          // internal = 1/0 if Thermo created them or not
                          // id = ID of Compute objects
                          // Compute * = ptrs to the Compute objects
-  int index_temp,index_press,index_drot,index_grot;
+  int index_temp,index_press,index_pe,index_drot,index_grot;
   int internal_drot,internal_grot;
-  char *id_temp,*id_press,*id_drot,*id_grot;
+  char *id_temp,*id_press,*id_pe,*id_drot,*id_grot;
-  class Compute *temperature,*pressure,*rotate_dipole,*rotate_gran;
+  class Compute *temperature,*pressure,*pe,*rotate_dipole,*rotate_gran;
 
   int ncompute;          // # of Compute objects called by thermo
   char **id_compute;     // their IDs

src/verlet.cpp

@@ -27,6 +27,7 @@
 #include "output.h"
 #include "update.h"
 #include "modify.h"
+#include "compute.h"
 #include "fix.h"
 #include "timer.h"
 #include "memory.h"
@@ -34,25 +35,10 @@
 
 using namespace LAMMPS_NS;
 
-#define MIN(A,B) ((A) < (B)) ? (A) : (B)
-#define MAX(A,B) ((A) > (B)) ? (A) : (B)
-
 /* ---------------------------------------------------------------------- */
 
 Verlet::Verlet(LAMMPS *lmp, int narg, char **arg) :
-  Integrate(lmp, narg, arg)
+  Integrate(lmp, narg, arg) {}
-{
-  fix_virial_every = NULL;
-  next_fix_virial = NULL;
-}
-
-/* ---------------------------------------------------------------------- */
-
-Verlet::~Verlet()
-{
-  delete [] fix_virial_every;
-  delete [] next_fix_virial;
-}
 
 /* ----------------------------------------------------------------------
    initialization before run
@@ -65,30 +51,32 @@ void Verlet::init()
   if (modify->nfix == 0)
     error->warning("No fixes defined, atoms won't move");
 
-  // setup virial computations for timestepping
+  // virial_style:
-  // virial_style = 1 if computed explicitly by pair
+  // 1 if computed explicitly by pair->compute via sum over pair interactions
-  //                2 if computed implicitly by pair (sum over ghost atoms)
+  // 2 if computed implicitly by pair->virial_compute via sum over ghost atoms
-  // virial_every = 1 if computed every timestep (NPT,NPH)
-  // fix arrays store info on fixes that need virial computed occasionally
 
   if (force->newton_pair) virial_style = 2;
   else virial_style = 1;
 
-  virial_every = 0;
+  // elist,vlist = list of computes for PE and pressure
-  nfix_virial = 0;
-  for (int i = 0; i < modify->nfix; i++)
-    if (modify->fix[i]->pressure_every == 1) virial_every = 1;
-    else if (modify->fix[i]->pressure_every > 1) nfix_virial++;
 
-  if (nfix_virial) {
+  delete [] elist;
-    delete [] fix_virial_every;
+  delete [] vlist;
-    delete [] next_fix_virial;
+  elist = vlist = NULL;
-    fix_virial_every = new int[nfix_virial];
+
-    next_fix_virial = new int[nfix_virial];
+  nelist = nvlist = 0;
-    nfix_virial = 0;
+  for (int i = 0; i < modify->ncompute; i++) {
-    for (int i = 0; i < modify->nfix; i++)
+    if (modify->compute[i]->peflag) nelist++;
-      if (modify->fix[i]->pressure_every > 1)
+    if (modify->compute[i]->pressflag) nvlist++;
-	fix_virial_every[nfix_virial++] = modify->fix[i]->pressure_every;
+  }
+
+  if (nelist) elist = new Compute*[nelist];
+  if (nvlist) vlist = new Compute*[nvlist];
+
+  nelist = nvlist = 0;
+  for (int i = 0; i < modify->ncompute; i++) {
+    if (modify->compute[i]->peflag) elist[nelist++] = modify->compute[i];
+    if (modify->compute[i]->pressflag) vlist[nvlist++] = modify->compute[i];
   }
 
   // set flags for what arrays to clear in force_clear()
@@ -127,8 +115,7 @@ void Verlet::setup()
 
   // compute all forces
 
-  int eflag = 1;
+  ev_set(update->ntimestep);
-  int vflag = virial_style;
   force_clear(vflag);
   
   if (force->pair) force->pair->compute(eflag,vflag);
@@ -149,21 +136,6 @@ void Verlet::setup()
 
   modify->setup();
   output->setup(1);
-
-  // setup virial computations for timestepping
-
-  int ntimestep = update->ntimestep;
-  next_virial = 0;
-  if (virial_every) next_virial = ntimestep + 1;
-  else {
-    for (int ivirial = 0; ivirial < nfix_virial; ivirial++) {
-      next_fix_virial[ivirial] = 
-	(ntimestep/fix_virial_every[ivirial])*fix_virial_every[ivirial] + 
-	fix_virial_every[ivirial];
-      if (ivirial) next_virial = MIN(next_virial,next_fix_virial[ivirial]);
-      else next_virial = next_fix_virial[0];
-    }
-  }
 }
 
 /* ----------------------------------------------------------------------
@@ -172,7 +144,7 @@ void Verlet::setup()
 
 void Verlet::iterate(int n)
 {
-  int eflag,vflag,nflag,ntimestep;
+  int nflag,ntimestep;
 
   for (int i = 0; i < n; i++) {
 
@@ -209,19 +181,10 @@ void Verlet::iterate(int n)
       timer->stamp(TIME_NEIGHBOR);
     }
 
-    // eflag/vflag = 0/1/2 for energy/virial computation
-
-    if (ntimestep == output->next_thermo) eflag = 1;
-    else eflag = 0;
-
-    if (ntimestep == output->next_thermo || ntimestep == next_virial) {
-      vflag = virial_style;
-      if (virial_every) next_virial++;
-      else next_virial = fix_virial(ntimestep);
-    } else vflag = 0;
-
     // force computations
 
+    ev_set(ntimestep);
+    ///printf("AAA %d %d %d\n",ntimestep,eflag,vflag);
     force_clear(vflag);
 
     timer->stamp();
@@ -299,20 +262,3 @@ void Verlet::force_clear(int vflag)
     }
   }
 }
-
-/* ----------------------------------------------------------------------
-   return next timestep virial should be computed
-   based on one or more fixes that need virial computed periodically
-------------------------------------------------------------------------- */
-
-int Verlet::fix_virial(int ntimestep)
-{
-  int next;
-  for (int ivirial = 0; ivirial < nfix_virial; ivirial++) {
-    if (ntimestep == next_fix_virial[ivirial])
-      next_fix_virial[ivirial] += fix_virial_every[ivirial];
-    if (ivirial) next = MIN(next,next_fix_virial[ivirial]);
-    else next = next_fix_virial[0];
-  }
-  return next;
-}

src/verlet.h

@@ -21,24 +21,16 @@ namespace LAMMPS_NS {
 class Verlet : public Integrate {
  public:
   Verlet(class LAMMPS *, int, char **);
-  ~Verlet();
+  ~Verlet() {}
   void init();
   void setup();
   void iterate(int);
 
  private:
-  int virial_style;                 // compute virial explicitly or implicitly
-  int virial_every;                 // 1 if virial computed every step
-  int next_virial;                  // next timestep to compute virial
-  int nfix_virial;                  // # of fixes that need virial occasionally
-  int *fix_virial_every;            // frequency they require it
-  int *next_fix_virial;             // next timestep they need it
   int triclinic;                    // 0 if domain is orthog, 1 if triclinic
-
   int torqueflag;                   // arrays to zero out every step
 
   void force_clear(int);
-  int fix_virial(int);
 };
 
 }