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

src/fix.cpp

@@ -17,6 +17,7 @@
 #include "atom.h"
 #include "group.h"
 #include "force.h"
+#include "comm.h"
 #include "atom_masks.h"
 #include "memory.h"
 #include "error.h"
@@ -58,6 +59,7 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   box_change_size = box_change_shape = box_change_domain = 0;
   thermo_energy = 0;
   rigid_flag = 0;
+  peatom_flag = 0;
   virial_flag = 0;
   no_change_box = 0;
   time_integrate = 0;
@@ -89,7 +91,8 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   // set vflag_atom = 0 b/c some fixes grow vatom in grow_arrays()
   //   which may occur outside of timestepping
 
-  maxvatom = 0;
+  maxeatom = maxvatom = 0;
+  eatom = NULL;
   vatom = NULL;
   vflag_atom = 0;
 
@@ -113,6 +116,7 @@ Fix::~Fix()
 
   delete [] id;
   delete [] style;
+  memory->destroy(eatom);
   memory->destroy(vatom);
 }
 
@@ -148,9 +152,65 @@ void Fix::modify_params(int narg, char **arg)
   }
 }
 
+/* ----------------------------------------------------------------------
+   setup for energy, virial computation
+   see integrate::ev_set() for values of eflag (0-3) and vflag (0-6)
+   fixes call this if use ev_tally()
+------------------------------------------------------------------------- */
+
+void Fix::ev_setup(int eflag, int vflag)
+{
+  int i,n;
+
+  evflag = 1;
+
+  eflag_either = eflag;
+  eflag_global = eflag % 2;
+  eflag_atom = eflag / 2;
+
+  vflag_either = vflag;
+  vflag_global = vflag % 4;
+  vflag_atom = vflag / 4;
+
+  // reallocate per-atom arrays if necessary
+
+  if (eflag_atom && atom->nlocal > maxeatom) {
+    maxeatom = atom->nmax;
+    memory->destroy(eatom);
+    memory->create(eatom,maxeatom,"fix:eatom");
+  }
+  if (vflag_atom && atom->nlocal > maxvatom) {
+    maxvatom = atom->nmax;
+    memory->destroy(vatom);
+    memory->create(vatom,maxvatom,6,"fix:vatom");
+  }
+
+  // zero accumulators
+  // no global energy variable to zero (unlike pair,bond,angle,etc)
+  // fixes tally it individually via fix_modify energy yes and compute_scalar()
+
+  if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0;
+  if (eflag_atom) {
+    n = atom->nlocal;
+    for (i = 0; i < n; i++) eatom[i] = 0.0;
+  }
+  if (vflag_atom) {
+    n = atom->nlocal;
+    for (i = 0; i < n; i++) {
+      vatom[i][0] = 0.0;
+      vatom[i][1] = 0.0;
+      vatom[i][2] = 0.0;
+      vatom[i][3] = 0.0;
+      vatom[i][4] = 0.0;
+      vatom[i][5] = 0.0;
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    setup for virial computation
    see integrate::ev_set() for values of vflag (0-6)
+   fixes call this if use v_tally()
 ------------------------------------------------------------------------- */
 
 void Fix::v_setup(int vflag)
@@ -187,12 +247,42 @@ void Fix::v_setup(int vflag)
 }
 
 /* ----------------------------------------------------------------------
-   tally virial into global and per-atom accumulators
+   tally per-atom energy and global/per-atom virial into accumulators
+   n = # of local owned atoms involved, with local indices in list
+   eng = total energy for the interaction involving total atoms
+   v = total virial for the interaction involving total atoms
+   increment per-atom energy of each atom in list by 1/total fraction
+   v_tally tallies virial
+   this method can be used when fix computes energy/forces in post_force()
+     e.g. fix cmap: compute energy and virial only on owned atoms
+       whether newton_bond is on or off
+     other procs will tally left-over fractions for atoms they own
+------------------------------------------------------------------------- */
+
+void Fix::ev_tally(int n, int *list, double total, double eng, double *v)
+{
+  int m;
+
+  if (eflag_atom) {
+    double fraction = eng/total;
+    for (int i = 0; i < n; i++)
+      eatom[list[i]] += fraction;
+  }
+
+  v_tally(n,list,total,v);
+}
+
+
+/* ----------------------------------------------------------------------
+   tally virial into global and per-atom accumulators
+   n = # of local owned atoms involved, with local indices in list
    v = total virial for the interaction involving total atoms
-   n = # of local atoms involved, with local indices in list
    increment global virial by n/total fraction
    increment per-atom virial of each atom in list by 1/total fraction
-   assumes other procs will tally left-over fractions
+   this method can be used when fix computes forces in post_force()
+     e.g. fix shake, fix rigid: compute virial only on owned atoms
+       whether newton_bond is on or off
+     other procs will tally left-over fractions for atoms they own
 ------------------------------------------------------------------------- */
 
 void Fix::v_tally(int n, int *list, double total, double *v)

src/fix.h

@@ -38,6 +38,7 @@ class Fix : protected Pointers {
   int thermo_energy;             // 1 if fix_modify enabled ThEng, 0 if not
   int nevery;                    // how often to call an end_of_step fix
   int rigid_flag;                // 1 if Fix integrates rigid bodies, 0 if not
+  int peatom_flag;               // 1 if Fix contributes per-atom eng, 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
   int time_integrate;            // 1 if fix performs time integration, 0 if no
@@ -89,7 +90,7 @@ class Fix : protected Pointers {
   int comm_border;               // size of border communication (0 if none)
 
   double virial[6];              // accumlated virial
-  double **vatom;                // accumulated per-atom virial
+  double *eatom,**vatom;         // accumulated per-atom energy/virial
 
   int restart_reset;             // 1 if restart just re-initialized fix
 
@@ -215,12 +216,15 @@ class Fix : protected Pointers {
   int instance_me;        // which Fix class instantiation I am
 
   int evflag;
-  int vflag_global,vflag_atom;
+  int eflag_either,eflag_global,eflag_atom;
-  int maxvatom;
+  int vflag_either,vflag_global,vflag_atom;
+  int maxeatom,maxvatom;
 
   int copymode;   // if set, do not deallocate during destruction
                   // required when classes are used as functors by Kokkos
 
+  void ev_setup(int, int);
+  void ev_tally(int, int *, double, double, double *);
   void v_setup(int);
   void v_tally(int, int *, double, double *);