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

src/REPLICA/tad.cpp

@@ -71,7 +71,6 @@ TAD::~TAD()
 
 void TAD::command(int narg, char **arg)
 {
-
   fix_event_list = NULL;
   n_event_list = 0;
   nmax_event_list = 0;

src/compute_displace_atom.cpp

@@ -20,6 +20,7 @@
 #include "domain.h"
 #include "modify.h"
 #include "fix.h"
+#include "fix_store.h"
 #include "memory.h"
 #include "error.h"
 
@@ -35,25 +36,43 @@ ComputeDisplaceAtom::ComputeDisplaceAtom(LAMMPS *lmp, int narg, char **arg) :
   peratom_flag = 1;
   size_peratom_cols = 4;
 
-  // create a new fix store/state style
+  // create a new fix STORE style
-  // id = compute-ID + store_state, fix group = compute group
+  // id = compute-ID + COMPUTE_STORE, fix group = compute group
 
-  int n = strlen(id) + strlen("_store_state") + 1;
+  int n = strlen(id) + strlen("_COMPUTE_STORE") + 1;
   id_fix = new char[n];
   strcpy(id_fix,id);
-  strcat(id_fix,"_store_state");
+  strcat(id_fix,"_COMPUTE_STORE");
 
-  char **newarg = new char*[7];
+  char **newarg = new char*[5];
   newarg[0] = id_fix;
   newarg[1] = group->names[igroup];
-  newarg[2] = (char *) "store/state";
+  newarg[2] = (char *) "STORE";
-  newarg[3] = (char *) "0";
+  newarg[3] = (char *) "1";
-  newarg[4] = (char *) "xu";
+  newarg[4] = (char *) "3";
-  newarg[5] = (char *) "yu";
+  modify->add_fix(5,newarg);
-  newarg[6] = (char *) "zu";
+  fix = (FixStore *) modify->fix[modify->nfix-1];
-  modify->add_fix(7,newarg);
   delete [] newarg;
 
+  // calculate xu,yu,zu for fix store array
+  // skip if reset from restart file
+
+  if (fix->restart_reset) fix->restart_reset = 0;
+  else {
+    double **xoriginal = fix->astore;
+
+    double **x = atom->x;
+    int *mask = atom->mask;
+    tagint *image = atom->image;
+    int nlocal = atom->nlocal;
+
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) domain->unmap(x[i],image[i],xoriginal[i]);
+      else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0;
+  }
+
+  // per-atom displacement array
+
   nmax = 0;
   displace = NULL;
 }
@@ -78,7 +97,7 @@ void ComputeDisplaceAtom::init()
 
   int ifix = modify->find_fix(id_fix);
   if (ifix < 0) error->all(FLERR,"Could not find compute displace/atom fix ID");
-  fix = modify->fix[ifix];
+  fix = (FixStore *) modify->fix[ifix];
 }
 
 /* ---------------------------------------------------------------------- */
@@ -100,7 +119,7 @@ void ComputeDisplaceAtom::compute_peratom()
   // original unwrapped position is stored by fix
   // for triclinic, need to unwrap current atom coord via h matrix
 
-  double **xoriginal = fix->array_atom;
+  double **xoriginal = fix->astore;
 
   double **x = atom->x;
   int *mask = atom->mask;

src/compute_displace_atom.h

@@ -36,7 +36,7 @@ class ComputeDisplaceAtom : public Compute {
   int nmax;
   double **displace;
   char *id_fix;
-  class Fix *fix;
+  class FixStore *fix;
 };
 
 }

src/compute_msd.cpp

@@ -19,6 +19,7 @@
 #include "domain.h"
 #include "modify.h"
 #include "fix.h"
+#include "fix_store.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
@@ -49,28 +50,57 @@ ComputeMSD::ComputeMSD(LAMMPS *lmp, int narg, char **arg) :
     } else error->all(FLERR,"Illegal compute msd command");
   }
 
-  // create a new fix store/state style with or without com keyword
+  // create a new fix STORE style
-  // id = compute-ID + store_state, fix group = compute group
+  // id = compute-ID + COMPUTE_STORE, fix group = compute group
 
-  int n = strlen(id) + strlen("_store_state") + 1;
+  int n = strlen(id) + strlen("_COMPUTE_STORE") + 1;
   id_fix = new char[n];
   strcpy(id_fix,id);
-  strcat(id_fix,"_store_state");
+  strcat(id_fix,"_COMPUTE_STORE");
 
-  char **newarg = new char*[9];
+  char **newarg = new char*[5];
   newarg[0] = id_fix;
   newarg[1] = group->names[igroup];
-  newarg[2] = (char *) "store/state";
+  newarg[2] = (char *) "STORE";
-  newarg[3] = (char *) "0";
+  newarg[3] = (char *) "1";
-  newarg[4] = (char *) "xu";
+  newarg[4] = (char *) "3";
-  newarg[5] = (char *) "yu";
+  modify->add_fix(5,newarg);
-  newarg[6] = (char *) "zu";
+  fix = (FixStore *) modify->fix[modify->nfix-1];
-  newarg[7] = (char *) "com";
-  newarg[8] = (char *) "yes";
-  if (comflag) modify->add_fix(9,newarg);
-  else modify->add_fix(7,newarg);
   delete [] newarg;
 
+  // calculate xu,yu,zu for fix store array
+  // skip if reset from restart file
+
+  if (fix->restart_reset) fix->restart_reset = 0;
+  else {
+    double **xoriginal = fix->astore;
+
+    double **x = atom->x;
+    int *mask = atom->mask;
+    tagint *image = atom->image;
+    int nlocal = atom->nlocal;
+
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) domain->unmap(x[i],image[i],xoriginal[i]);
+      else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0;
+
+    // adjust for COM if requested
+    
+    if (comflag) {
+      double cm[3];
+      masstotal = group->mass(igroup);
+      group->xcm(igroup,masstotal,cm);
+      for (int i = 0; i < nlocal; i++)
+        if (mask[i] & groupbit) {
+          xoriginal[i][0] -= cm[0];
+          xoriginal[i][1] -= cm[1];
+          xoriginal[i][2] -= cm[2];
+        }
+    }
+  }
+
+  // displacement vector
+
   vector = new double[4];
 }
 
@@ -94,7 +124,7 @@ void ComputeMSD::init()
 
   int ifix = modify->find_fix(id_fix);
   if (ifix < 0) error->all(FLERR,"Could not find compute msd fix ID");
-  fix = modify->fix[ifix];
+  fix = (FixStore *) modify->fix[ifix];
 
   // nmsd = # of atoms in group
 
@@ -129,7 +159,7 @@ void ComputeMSD::compute_vector()
   // relative to center of mass if comflag is set
   // for triclinic, need to unwrap current atom coord via h matrix
 
-  double **xoriginal = fix->array_atom;
+  double **xoriginal = fix->astore;
 
   double **x = atom->x;
   int *mask = atom->mask;

src/compute_msd.h

@@ -35,7 +35,7 @@ class ComputeMSD : public Compute {
   int comflag,nmsd;
   double masstotal;
   char *id_fix;
-  class Fix *fix;
+  class FixStore *fix;
 };
 
 }

src/fix.cpp

@@ -66,6 +66,7 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   peratom_flag = local_flag = 0;
 
   comm_forward = comm_reverse = 0;
+  restart_reset = 0;
 
   maxvatom = 0;
   vatom = NULL;

src/fix.h

@@ -76,6 +76,7 @@ class Fix : protected Pointers {
   double virial[6];              // accumlated virial
   double **vatom;                // accumulated per-atom virial
 
+  int restart_reset;             // 1 if restart just re-initialized fix
   unsigned int datamask;
   unsigned int datamask_ext;
 

src/fix_store.cpp

@@ -0,0 +1,189 @@
+/* ----------------------------------------------------------------------
+   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 "string.h"
+#include "fix_store.h"
+#include "atom.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+/* ---------------------------------------------------------------------- */
+
+FixStore::FixStore(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
+{
+  if (narg != 5) error->all(FLERR,"Illegal fix store command");
+
+  // syntax: id group style 0/1 nvalue
+
+  restart_peratom = atoi(arg[3]);
+  nvalues = atoi(arg[4]);
+
+  vecflag = 0;
+  if (nvalues == 1) vecflag = 1;
+
+  // perform initial allocation of atom-based array
+  // register with Atom class
+
+  vstore = NULL;
+  astore = NULL;
+  grow_arrays(atom->nmax);
+  atom->add_callback(0);
+  if (restart_peratom) atom->add_callback(1);
+
+  // zero the storage
+  // since may be exchanged before filled by caller, e.g. fix store/force
+
+  int nlocal = atom->nlocal;
+  if (vecflag)
+    for (int i = 0; i < nlocal; i++)
+      vstore[i] = 0.0;
+  else 
+    for (int i = 0; i < nlocal; i++)
+      for (int j = 0; j < nvalues; j++)
+	astore[i][j] = 0.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixStore::~FixStore()
+{
+  // unregister callbacks to this fix from Atom class
+
+  atom->delete_callback(id,0);
+  if (restart_peratom) atom->delete_callback(id,1);
+
+  memory->destroy(vstore);
+  memory->destroy(astore);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixStore::setmask()
+{
+  int mask = 0;
+  return mask;
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double FixStore::memory_usage()
+{
+  double bytes = atom->nmax*nvalues * sizeof(double);
+  return bytes;
+}
+
+/* ----------------------------------------------------------------------
+   allocate atom-based array
+------------------------------------------------------------------------- */
+
+void FixStore::grow_arrays(int nmax)
+{
+  if (vecflag) memory->grow(vstore,nmax,"store:vstore");
+  else memory->grow(astore,nmax,nvalues,"store:astore");
+}
+
+/* ----------------------------------------------------------------------
+   copy values within local atom-based array
+------------------------------------------------------------------------- */
+
+void FixStore::copy_arrays(int i, int j, int delflag)
+{
+  if (vecflag) vstore[j] = vstore[i];
+  else
+    for (int m = 0; m < nvalues; m++)
+      astore[j][m] = astore[i][m];
+}
+
+/* ----------------------------------------------------------------------
+   pack values in local atom-based array for exchange with another proc
+------------------------------------------------------------------------- */
+
+int FixStore::pack_exchange(int i, double *buf)
+{
+  if (vecflag) buf[0] = vstore[i];
+  else
+    for (int m = 0; m < nvalues; m++)
+      buf[m] = astore[i][m];
+  return nvalues;
+}
+
+/* ----------------------------------------------------------------------
+   unpack values in local atom-based array from exchange with another proc
+------------------------------------------------------------------------- */
+
+int FixStore::unpack_exchange(int nlocal, double *buf)
+{
+  if (vecflag) vstore[nlocal] = buf[0];
+  else 
+    for (int m = 0; m < nvalues; m++)
+      astore[nlocal][m] = buf[m];
+  return nvalues;
+}
+
+/* ----------------------------------------------------------------------
+   pack values in local atom-based arrays for restart file
+------------------------------------------------------------------------- */
+
+int FixStore::pack_restart(int i, double *buf)
+{
+  buf[0] = nvalues+1;
+  if (vecflag) buf[1] = vstore[i];
+  else
+    for (int m = 0; m < nvalues; m++)
+      buf[m+1] = astore[i][m];
+  return nvalues+1;
+}
+
+/* ----------------------------------------------------------------------
+   unpack values from atom->extra array to restart the fix
+------------------------------------------------------------------------- */
+
+void FixStore::unpack_restart(int nlocal, int nth)
+{
+  double **extra = atom->extra;
+
+  // skip to Nth set of extra values
+
+  int m = 0;
+  for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
+  m++;
+
+  if (vecflag) vstore[nlocal] = extra[nlocal][m];
+  else
+    for (int i = 0; i < nvalues; i++) 
+      astore[nlocal][i] = extra[nlocal][m++];
+}
+
+/* ----------------------------------------------------------------------
+   maxsize of any atom's restart data
+------------------------------------------------------------------------- */
+
+int FixStore::maxsize_restart()
+{
+  return nvalues+1;
+}
+
+/* ----------------------------------------------------------------------
+   size of atom nlocal's restart data
+------------------------------------------------------------------------- */
+
+int FixStore::size_restart(int nlocal)
+{
+  return nvalues+1;
+}

src/fix_store.h

@@ -0,0 +1,54 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(STORE,FixStore)
+
+#else
+
+#ifndef LMP_FIX_STORE_H
+#define LMP_FIX_STORE_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixStore : public Fix {
+ public:
+  double *vstore;        // vector storage if nvalues = 1
+  double **astore;       // array storage if nvalues > 1
+
+  FixStore(class LAMMPS *, int, char **);
+  ~FixStore();
+  int setmask();
+
+  double memory_usage();
+  void grow_arrays(int);
+  void copy_arrays(int, int, int);
+  int pack_exchange(int, double *);
+  int unpack_exchange(int, double *);
+  int pack_restart(int, double *);
+  void unpack_restart(int, int);
+  int size_restart(int);
+  int maxsize_restart();
+
+ private:
+  int nvalues;                  // total # of values per atom
+  int vecflag;                  // 1 if nvalues = 1
+};
+
+}
+
+#endif
+#endif

src/fix_store_state.cpp

@@ -126,60 +126,74 @@ FixStoreState::FixStoreState(LAMMPS *lmp, int narg, char **arg) :
 
     } else if (strcmp(arg[iarg],"q") == 0) {
       if (!atom->q_flag)
-        error->all(FLERR,"Fix store/state for atom property that isn't allocated");
+        error->all(FLERR,
+                   "Fix store/state for atom property that isn't allocated");
       pack_choice[nvalues++] = &FixStoreState::pack_q;
     } else if (strcmp(arg[iarg],"mux") == 0) {
       if (!atom->mu_flag)
-        error->all(FLERR,"Fix store/state for atom property that isn't allocated");
+        error->all(FLERR,
+                   "Fix store/state for atom property that isn't allocated");
       pack_choice[nvalues++] = &FixStoreState::pack_mux;
     } else if (strcmp(arg[iarg],"muy") == 0) {
       if (!atom->mu_flag)
-        error->all(FLERR,"Fix store/state for atom property that isn't allocated");
+        error->all(FLERR,
+                   "Fix store/state for atom property that isn't allocated");
       pack_choice[nvalues++] = &FixStoreState::pack_muy;
     } else if (strcmp(arg[iarg],"muz") == 0) {
       if (!atom->mu_flag)
-        error->all(FLERR,"Fix store/state for atom property that isn't allocated");
+        error->all(FLERR,
+                   "Fix store/state for atom property that isn't allocated");
       pack_choice[nvalues++] = &FixStoreState::pack_muz;
 
     } else if (strcmp(arg[iarg],"radius") == 0) {
       if (!atom->radius_flag)
-        error->all(FLERR,"Fix store/state for atom property that isn't allocated");
+        error->all(FLERR,
+                   "Fix store/state for atom property that isn't allocated");
       pack_choice[nvalues++] = &FixStoreState::pack_radius;
     } else if (strcmp(arg[iarg],"omegax") == 0) {
       if (!atom->omega_flag)
-        error->all(FLERR,"Fix store/state for atom property that isn't allocated");
+        error->all(FLERR,
+                   "Fix store/state for atom property that isn't allocated");
       pack_choice[nvalues++] = &FixStoreState::pack_omegax;
     } else if (strcmp(arg[iarg],"omegay") == 0) {
       if (!atom->omega_flag)
-        error->all(FLERR,"Fix store/state for atom property that isn't allocated");
+        error->all(FLERR,
+                   "Fix store/state for atom property that isn't allocated");
       pack_choice[nvalues++] = &FixStoreState::pack_omegay;
     } else if (strcmp(arg[iarg],"omegaz") == 0) {
       if (!atom->omega_flag)
-        error->all(FLERR,"Fix store/state for atom property that isn't allocated");
+        error->all(FLERR,
+                   "Fix store/state for atom property that isn't allocated");
       pack_choice[nvalues++] = &FixStoreState::pack_omegaz;
     } else if (strcmp(arg[iarg],"angmomx") == 0) {
       if (!atom->angmom_flag)
-        error->all(FLERR,"Fix store/state for atom property that isn't allocated");
+        error->all(FLERR,
+                   "Fix store/state for atom property that isn't allocated");
       pack_choice[nvalues++] = &FixStoreState::pack_angmomx;
     } else if (strcmp(arg[iarg],"angmomy") == 0) {
       if (!atom->angmom_flag)
-        error->all(FLERR,"Fix store/state for atom property that isn't allocated");
+        error->all(FLERR,
+                   "Fix store/state for atom property that isn't allocated");
       pack_choice[nvalues++] = &FixStoreState::pack_angmomy;
     } else if (strcmp(arg[iarg],"angmomz") == 0) {
       if (!atom->angmom_flag)
-        error->all(FLERR,"Fix store/state for atom property that isn't allocated");
+        error->all(FLERR,
+                   "Fix store/state for atom property that isn't allocated");
       pack_choice[nvalues++] = &FixStoreState::pack_angmomz;
     } else if (strcmp(arg[iarg],"tqx") == 0) {
       if (!atom->torque_flag)
-        error->all(FLERR,"Fix store/state for atom property that isn't allocated");
+        error->all(FLERR,
+                   "Fix store/state for atom property that isn't allocated");
       pack_choice[nvalues++] = &FixStoreState::pack_tqx;
     } else if (strcmp(arg[iarg],"tqy") == 0) {
       if (!atom->torque_flag)
-        error->all(FLERR,"Fix store/state for atom property that isn't allocated");
+        error->all(FLERR,
+                   "Fix store/state for atom property that isn't allocated");
       pack_choice[nvalues++] = &FixStoreState::pack_tqy;
     } else if (strcmp(arg[iarg],"tqz") == 0) {
       if (!atom->torque_flag)
-        error->all(FLERR,"Fix store/state for atom property that isn't allocated");
+        error->all(FLERR,
+                   "Fix store/state for atom property that isn't allocated");
       pack_choice[nvalues++] = &FixStoreState::pack_tqz;
 
     } else if (strncmp(arg[iarg],"c_",2) == 0 ||
@@ -242,26 +256,34 @@ FixStoreState::FixStoreState(LAMMPS *lmp, int narg, char **arg) :
         error->all(FLERR,"Fix store/state compute does not "
                    "calculate a per-atom vector");
       if (argindex[i] && modify->compute[icompute]->size_peratom_cols == 0)
-        error->all(FLERR,"Fix store/state compute does not "
+        error->all(FLERR,
+                   "Fix store/state compute does not "
                    "calculate a per-atom array");
       if (argindex[i] &&
           argindex[i] > modify->compute[icompute]->size_peratom_cols)
-        error->all(FLERR,"Fix store/state compute array is accessed out-of-range");
+        error->all(FLERR,
+                   "Fix store/state compute array is accessed out-of-range");
 
     } else if (which[i] == FIX) {
       int ifix = modify->find_fix(ids[i]);
       if (ifix < 0)
-        error->all(FLERR,"Fix ID for fix store/state does not exist");
+        error->all(FLERR,
+                   "Fix ID for fix store/state does not exist");
       if (modify->fix[ifix]->peratom_flag == 0)
-        error->all(FLERR,"Fix store/state fix does not calculate per-atom values");
+        error->all(FLERR,
+                   "Fix store/state fix does not calculate per-atom values");
       if (argindex[i] == 0 && modify->fix[ifix]->size_peratom_cols != 0)
-        error->all(FLERR,"Fix store/state fix does not calculate a per-atom vector");
+        error->all(FLERR,
+                   "Fix store/state fix does not calculate a per-atom vector");
       if (argindex[i] && modify->fix[ifix]->size_peratom_cols == 0)
-        error->all(FLERR,"Fix store/state fix does not calculate a per-atom array");
+        error->all(FLERR,
+                   "Fix store/state fix does not calculate a per-atom array");
       if (argindex[i] && argindex[i] > modify->fix[ifix]->size_peratom_cols)
-        error->all(FLERR,"Fix store/state fix array is accessed out-of-range");
+        error->all(FLERR,
+                   "Fix store/state fix array is accessed out-of-range");
       if (nevery % modify->fix[ifix]->peratom_freq)
-        error->all(FLERR,"Fix for fix store/state not computed at compatible time");
+        error->all(FLERR,
+                   "Fix for fix store/state not computed at compatible time");
 
     } else if (which[i] == VARIABLE) {
       int ivariable = input->variable->find(ids[i]);

src/modify.cpp

@@ -728,6 +728,7 @@ void Modify::add_fix(int narg, char **arg, char *suffix)
         strcmp(style_restart_peratom[i],fix[ifix]->style) == 0) {
       for (int j = 0; j < atom->nlocal; j++)
         fix[ifix]->unpack_restart(j,index_restart_peratom[i]);
+      fix[ifix]->restart_reset = 1;
       if (comm->me == 0) {
         char *str = (char *) ("Resetting per-atom state of Fix %s Style %s "
                      "from restart file info\n");

src/my_pool.h

@@ -32,7 +32,7 @@ methods:
    T *get(index) = return ptr/index to unused chunk of size maxchunk
    T *get(N,index) = return ptr/index to unused chunk of size N
                      minchunk < N < maxchunk required
-   put(index) = return indexed chunk to pool
+   put(index) = return indexed chunk to pool (same index returned by get)
    int size() = return total size of allocated pages in bytes
 public varaibles:
    ndatum = total # of stored datums