From afac9a0f14a9c3a447bd7fb471dc759c36f3709d Mon Sep 17 00:00:00 2001 From: sjplimp Date: Sat, 6 Apr 2013 23:08:09 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9754 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/REPLICA/tad.cpp | 1 - src/compute_displace_atom.cpp | 45 +++++--- src/compute_displace_atom.h | 2 +- src/compute_msd.cpp | 62 ++++++++--- src/compute_msd.h | 2 +- src/fix.cpp | 1 + src/fix.h | 1 + src/fix_store.cpp | 189 ++++++++++++++++++++++++++++++++++ src/fix_store.h | 54 ++++++++++ src/fix_store_state.cpp | 66 ++++++++---- src/modify.cpp | 1 + src/my_pool.h | 2 +- 12 files changed, 371 insertions(+), 55 deletions(-) create mode 100644 src/fix_store.cpp create mode 100644 src/fix_store.h diff --git a/src/REPLICA/tad.cpp b/src/REPLICA/tad.cpp index 23633af5dc..ab93a5906e 100644 --- a/src/REPLICA/tad.cpp +++ b/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; diff --git a/src/compute_displace_atom.cpp b/src/compute_displace_atom.cpp index b89acd7f56..998c3e82b3 100644 --- a/src/compute_displace_atom.cpp +++ b/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 - // id = compute-ID + store_state, fix group = compute group + // create a new fix STORE style + // 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[3] = (char *) "0"; - newarg[4] = (char *) "xu"; - newarg[5] = (char *) "yu"; - newarg[6] = (char *) "zu"; - modify->add_fix(7,newarg); + newarg[2] = (char *) "STORE"; + newarg[3] = (char *) "1"; + newarg[4] = (char *) "3"; + modify->add_fix(5,newarg); + fix = (FixStore *) modify->fix[modify->nfix-1]; 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; diff --git a/src/compute_displace_atom.h b/src/compute_displace_atom.h index 74e41cce45..be3fbe5074 100644 --- a/src/compute_displace_atom.h +++ b/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; }; } diff --git a/src/compute_msd.cpp b/src/compute_msd.cpp index c77da886ba..05fa21e091 100644 --- a/src/compute_msd.cpp +++ b/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 - // id = compute-ID + store_state, fix group = compute group + // create a new fix STORE style + // 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[3] = (char *) "0"; - newarg[4] = (char *) "xu"; - newarg[5] = (char *) "yu"; - newarg[6] = (char *) "zu"; - newarg[7] = (char *) "com"; - newarg[8] = (char *) "yes"; - if (comflag) modify->add_fix(9,newarg); - else modify->add_fix(7,newarg); + newarg[2] = (char *) "STORE"; + newarg[3] = (char *) "1"; + newarg[4] = (char *) "3"; + modify->add_fix(5,newarg); + fix = (FixStore *) modify->fix[modify->nfix-1]; 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; diff --git a/src/compute_msd.h b/src/compute_msd.h index d2a35a4d62..54c044a7c2 100644 --- a/src/compute_msd.h +++ b/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; }; } diff --git a/src/fix.cpp b/src/fix.cpp index 6ba167e221..56045e8573 100644 --- a/src/fix.cpp +++ b/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; diff --git a/src/fix.h b/src/fix.h index eb198cae6c..570b7f005f 100644 --- a/src/fix.h +++ b/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; diff --git a/src/fix_store.cpp b/src/fix_store.cpp new file mode 100644 index 0000000000..be80e0ebf2 --- /dev/null +++ b/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 (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; +} diff --git a/src/fix_store.h b/src/fix_store.h new file mode 100644 index 0000000000..f3b140152c --- /dev/null +++ b/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 diff --git a/src/fix_store_state.cpp b/src/fix_store_state.cpp index 0fab57d7f9..94353500a7 100644 --- a/src/fix_store_state.cpp +++ b/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]); diff --git a/src/modify.cpp b/src/modify.cpp index 75bcf1325b..8c7e296cd0 100644 --- a/src/modify.cpp +++ b/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"); diff --git a/src/my_pool.h b/src/my_pool.h index 4b384e5525..83864a6426 100644 --- a/src/my_pool.h +++ b/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