tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
sjplimp
2015-02-13 16:53:48 +00:00
parent 4682db87f2
commit 8f818715f4
31 changed files with 5531 additions and 1742 deletions
Download Patch File Download Diff File Expand all files Collapse all files

103
src/compute.cpp
View File

@ -95,9 +95,7 @@ Compute::Compute(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
ntime = maxtime = 0;
tlist = NULL;
// setup map for molecule IDs
molmap = NULL;
// data masks
datamask = ALL_MASK;
datamask_ext = ALL_MASK;
@ -114,9 +112,7 @@ Compute::~Compute()
{
delete [] id;
delete [] style;
memory->destroy(tlist);
memory->destroy(molmap);
}
/* ---------------------------------------------------------------------- */
@ -230,100 +226,3 @@ void Compute::clearstep()
{
ntime = 0;
}
/* ----------------------------------------------------------------------
identify molecule IDs with atoms in group
warn if any atom in group has molecule ID = 0
warn if any molecule has only some atoms in group
return Ncount = # of molecules with atoms in group
set molmap to NULL if molecule IDs include all in range from 1 to Ncount
else: molecule IDs range from idlo to idhi
set molmap to vector of length idhi-idlo+1
molmap[id-idlo] = index from 0 to Ncount-1
return idlo and idhi
------------------------------------------------------------------------- */
int Compute::molecules_in_group(tagint &idlo, tagint &idhi)
{
int i;
memory->destroy(molmap);
molmap = NULL;
// find lo/hi molecule ID for any atom in group
// warn if atom in group has ID = 0
tagint *molecule = atom->molecule;
int *mask = atom->mask;
int nlocal = atom->nlocal;
tagint lo = BIG;
tagint hi = -BIG;
int flag = 0;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (molecule[i] == 0) flag = 1;
lo = MIN(lo,molecule[i]);
hi = MAX(hi,molecule[i]);
}
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
if (flagall && comm->me == 0)
error->warning(FLERR,"Atom with molecule ID = 0 included in "
"compute molecule group");
MPI_Allreduce(&lo,&idlo,1,MPI_LMP_TAGINT,MPI_MIN,world);
MPI_Allreduce(&hi,&idhi,1,MPI_LMP_TAGINT,MPI_MAX,world);
if (idlo == BIG) return 0;
// molmap = vector of length nlen
// set to 1 for IDs that appear in group across all procs, else 0
tagint nlen_tag = idhi-idlo+1;
if (nlen_tag > MAXSMALLINT)
error->all(FLERR,"Too many molecules for compute");
int nlen = (int) nlen_tag;
memory->create(molmap,nlen,"compute:molmap");
for (i = 0; i < nlen; i++) molmap[i] = 0;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
molmap[molecule[i]-idlo] = 1;
int *molmapall;
memory->create(molmapall,nlen,"compute:molmapall");
MPI_Allreduce(molmap,molmapall,nlen,MPI_INT,MPI_MAX,world);
// nmolecules = # of non-zero IDs in molmap
// molmap[i] = index of molecule, skipping molecules not in group with -1
int nmolecules = 0;
for (i = 0; i < nlen; i++)
if (molmapall[i]) molmap[i] = nmolecules++;
else molmap[i] = -1;
memory->destroy(molmapall);
// warn if any molecule has some atoms in group and some not in group
flag = 0;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) continue;
if (molecule[i] < idlo || molecule[i] > idhi) continue;
if (molmap[molecule[i]-idlo] >= 0) flag = 1;
}
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
if (flagall && comm->me == 0)
error->warning(FLERR,
"One or more compute molecules has atoms not in group");
// if molmap simply stores 1 to Nmolecules, then free it
if (idlo == 1 && idhi == nmolecules && nlen == nmolecules) {
memory->destroy(molmap);
molmap = NULL;
}
return nmolecules;
}

4
src/compute.h
View File

@ -148,10 +148,6 @@ class Compute : protected Pointers {
double **vbiasall; // stored velocity bias for all atoms
int maxbias; // size of vbiasall array
int *molmap; // convert molecule ID to local index
int molecules_in_group(tagint &, tagint &);
inline int sbmask(int j) {
return j >> SBBITS & 3;
}

347
src/compute_atom_molecule.cpp
View File

@ -1,347 +0,0 @@
/* ----------------------------------------------------------------------
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 "string.h"
#include "stdlib.h"
#include "compute_atom_molecule.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "compute.h"
#include "fix.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
enum{COMPUTE,FIX,VARIABLE};
#define INVOKED_PERATOM 8
/* ---------------------------------------------------------------------- */
ComputeAtomMolecule::
ComputeAtomMolecule(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg < 4) error->all(FLERR,"Illegal compute atom/molecule command");
if (atom->molecular == 0)
error->all(FLERR,"Compute atom/molecule requires molecular atom style");
// parse args
which = new int[narg-3];
argindex = new int[narg-3];
ids = new char*[narg-3];
value2index = new int[narg-3];
nvalues = 0;
int iarg = 3;
while (iarg < narg) {
if (strncmp(arg[iarg],"c_",2) == 0 ||
strncmp(arg[iarg],"f_",2) == 0 ||
strncmp(arg[iarg],"v_",2) == 0) {
if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE;
else if (arg[iarg][0] == 'f') which[nvalues] = FIX;
else if (arg[iarg][0] == 'v') which[nvalues] = VARIABLE;
int n = strlen(arg[iarg]);
char *suffix = new char[n];
strcpy(suffix,&arg[iarg][2]);
char *ptr = strchr(suffix,'[');
if (ptr) {
if (suffix[strlen(suffix)-1] != ']')
error->all(FLERR,"Illegal compute reduce command");
argindex[nvalues] = atoi(ptr+1);
*ptr = '\0';
} else argindex[nvalues] = 0;
n = strlen(suffix) + 1;
ids[nvalues] = new char[n];
strcpy(ids[nvalues],suffix);
nvalues++;
delete [] suffix;
} else error->all(FLERR,"Illegal compute atom/molecule command");
iarg++;
}
// setup and error check
for (int i = 0; i < nvalues; i++) {
if (which[i] == COMPUTE) {
int icompute = modify->find_compute(ids[i]);
if (icompute < 0)
error->all(FLERR,"Compute ID for compute atom/molecule does not exist");
if (modify->compute[icompute]->peratom_flag == 0)
error->all(FLERR,"Compute atom/molecule compute does not "
"calculate per-atom values");
if (argindex[i] == 0 &&
modify->compute[icompute]->size_peratom_cols != 0)
error->all(FLERR,"Compute atom/molecule compute does not "
"calculate a per-atom vector");
if (argindex[i] && modify->compute[icompute]->size_peratom_cols == 0)
error->all(FLERR,"Compute atom/molecule compute does not "
"calculate a per-atom array");
if (argindex[i] &&
argindex[i] > modify->compute[icompute]->size_peratom_cols)
error->all(FLERR,"Compute atom/molecule 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 compute atom/molecule does not exist");
if (modify->fix[ifix]->peratom_flag == 0)
error->all(FLERR,"Compute atom/molecule fix does not "
"calculate per-atom values");
if (argindex[i] == 0 &&
modify->fix[ifix]->size_peratom_cols != 0)
error->all(FLERR,"Compute atom/molecule fix does not "
"calculate a per-atom vector");
if (argindex[i] && modify->fix[ifix]->size_peratom_cols == 0)
error->all(FLERR,"Compute atom/molecule fix does not "
"calculate a per-atom array");
if (argindex[i] &&
argindex[i] > modify->fix[ifix]->size_peratom_cols)
error->all(FLERR,
"Compute atom/molecule fix array is accessed out-of-range");
} else if (which[i] == VARIABLE) {
int ivariable = input->variable->find(ids[i]);
if (ivariable < 0)
error->all(FLERR,
"Variable name for compute atom/molecule does not exist");
if (input->variable->atomstyle(ivariable) == 0)
error->all(FLERR,"Compute atom/molecule variable is not "
"atom-style variable");
}
}
// setup molecule-based data
nmolecules = molecules_in_group(idlo,idhi);
vone = vector = NULL;
aone = array = NULL;
if (nvalues == 1) {
memory->create(vone,nmolecules,"atom/molecule:vone");
memory->create(vector,nmolecules,"atom/molecule:vector");
vector_flag = 1;
size_vector = nmolecules;
extvector = 0;
} else {
memory->create(aone,nmolecules,nvalues,"atom/molecule:aone");
memory->create(array,nmolecules,nvalues,"atom/molecule:array");
array_flag = 1;
size_array_rows = nmolecules;
size_array_cols = nvalues;
extarray = 0;
}
maxatom = 0;
scratch = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeAtomMolecule::~ComputeAtomMolecule()
{
delete [] which;
delete [] argindex;
for (int m = 0; m < nvalues; m++) delete [] ids[m];
delete [] ids;
delete [] value2index;
memory->destroy(vone);
memory->destroy(vector);
memory->destroy(aone);
memory->destroy(array);
memory->destroy(scratch);
}
/* ---------------------------------------------------------------------- */
void ComputeAtomMolecule::init()
{
int ntmp = molecules_in_group(idlo,idhi);
if (ntmp != nmolecules)
error->all(FLERR,"Molecule count changed in compute atom/molecule");
// set indices and check validity of all computes,fixes,variables
for (int m = 0; m < nvalues; m++) {
if (which[m] == COMPUTE) {
int icompute = modify->find_compute(ids[m]);
if (icompute < 0)
error->all(FLERR,"Compute ID for compute atom/molecule does not exist");
value2index[m] = icompute;
} else if (which[m] == FIX) {
int ifix = modify->find_fix(ids[m]);
if (ifix < 0)
error->all(FLERR,"Fix ID for compute atom/molecule does not exist");
value2index[m] = ifix;
} else if (which[m] == VARIABLE) {
int ivariable = input->variable->find(ids[m]);
if (ivariable < 0)
error->all(FLERR,
"Variable name for compute atom/molecule does not exist");
value2index[m] = ivariable;
} else value2index[m] = -1;
}
}
/* ---------------------------------------------------------------------- */
void ComputeAtomMolecule::compute_vector()
{
int i,j,n;
tagint imol;
invoked_vector = update->ntimestep;
for (n = 0; n < nmolecules; n++) vone[n] = 0.0;
compute_one(0);
tagint *molecule = atom->molecule;
int *mask = atom->mask;
int nlocal = atom->nlocal;
j = 0;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
vone[imol] += peratom[j];
}
j += nstride;
}
int me;
MPI_Comm_rank(world,&me);
MPI_Allreduce(vone,vector,nmolecules,MPI_DOUBLE,MPI_SUM,world);
}
/* ---------------------------------------------------------------------- */
void ComputeAtomMolecule::compute_array()
{
int i,j,m,n;
tagint imol;
invoked_array = update->ntimestep;
tagint *molecule = atom->molecule;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (m = 0; m < nvalues; m++) {
for (n = 0; n < nmolecules; n++) aone[n][m] = 0.0;
compute_one(m);
j = 0;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
aone[imol][m] += peratom[j];
}
j += nstride;
}
}
if (array)
MPI_Allreduce(&aone[0][0],&array[0][0],nvalues*nmolecules,
MPI_DOUBLE,MPI_SUM,world);
}
/* ----------------------------------------------------------------------
calculate per-atom values for one input M
invoke the appropriate compute,fix,variable
reallocate scratch if necessary for per-atom variable scratch space
------------------------------------------------------------------------- */
void ComputeAtomMolecule::compute_one(int m)
{
int vidx = value2index[m];
int aidx = argindex[m];
// invoke compute if not previously invoked
if (which[m] == COMPUTE) {
Compute *compute = modify->compute[vidx];
if (!(compute->invoked_flag & INVOKED_PERATOM)) {
compute->compute_peratom();
compute->invoked_flag |= INVOKED_PERATOM;
}
if (aidx == 0) {
peratom = compute->vector_atom;
nstride = 1;
} else {
peratom = &compute->array_atom[0][aidx-1];
nstride = compute->size_peratom_cols;
}
// access fix fields, check if fix frequency is a match
} else if (which[m] == FIX) {
if (update->ntimestep % modify->fix[vidx]->peratom_freq)
error->all(FLERR,"Fix used in compute atom/molecule not computed "
"at compatible time");
Fix *fix = modify->fix[vidx];
if (aidx == 0) {
peratom = fix->vector_atom;
nstride = 1;
} else {
peratom = &fix->array_atom[0][aidx-1];
nstride = fix->size_peratom_cols;
}
// evaluate atom-style variable
} else if (which[m] == VARIABLE) {
if (atom->nlocal > maxatom) {
maxatom = atom->nmax;
memory->destroy(scratch);
memory->create(scratch,maxatom,"atom/molecule:scratch");
}
peratom = scratch;
input->variable->compute_atom(vidx,igroup,peratom,1,0);
nstride = 1;
}
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputeAtomMolecule::memory_usage()
{
double bytes = (bigint) nmolecules * 2*nvalues * sizeof(double);
if (molmap) bytes += (idhi-idlo+1) * sizeof(int);
bytes += maxatom * sizeof(double);
return bytes;
}

126
src/compute_atom_molecule.h
View File

@ -1,126 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
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 COMPUTE_CLASS
ComputeStyle(atom/molecule,ComputeAtomMolecule)
#else
#ifndef LMP_COMPUTE_ATOM_MOLECULE_H
#define LMP_COMPUTE_ATOM_MOLECULE_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeAtomMolecule : public Compute {
public:
ComputeAtomMolecule(class LAMMPS *, int, char **);
~ComputeAtomMolecule();
void init();
void compute_vector();
void compute_array();
double memory_usage();
private:
int nvalues,nmolecules;
tagint idlo,idhi;
int *which,*argindex,*value2index;
char **ids;
int nstride,maxatom;
double *vone;
double **aone;
double *scratch;
double *peratom;
void compute_one(int);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Compute atom/molecule requires molecular atom style
Self-explanatory.
E: Compute ID for compute atom/molecule does not exist
Self-explanatory.
E: Compute atom/molecule compute does not calculate per-atom values
Self-explanatory.
E: Compute atom/molecule compute does not calculate a per-atom vector
Self-explanatory.
E: Compute atom/molecule compute does not calculate a per-atom array
Self-explanatory.
E: Compute atom/molecule compute array is accessed out-of-range
Self-explanatory.
E: Fix ID for compute atom/molecule does not exist
Self-explanatory.
E: Compute atom/molecule fix does not calculate per-atom values
Self-explanatory.
E: Compute atom/molecule fix does not calculate a per-atom vector
Self-explanatory.
E: Compute atom/molecule fix does not calculate a per-atom array
Self-explanatory.
E: Compute atom/molecule fix array is accessed out-of-range
Self-explanatory.
E: Variable name for compute atom/molecule does not exist
Self-explanatory.
E: Compute atom/molecule variable is not atom-style variable
Self-explanatory.
E: Molecule count changed in compute atom/molecule
Number of molecules must remain constant over time.
E: Fix used in compute atom/molecule not computed at compatible time
The fix must produce per-atom quantities on timesteps that the compute
needs them.
*/

1557
src/compute_chunk_atom.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

125
src/compute_chunk_atom.h Normal file
View File

@ -0,0 +1,125 @@
/* -*- c++ -*- ----------------------------------------------------------
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 COMPUTE_CLASS
ComputeStyle(chunk/atom,ComputeChunkAtom)
#else
#ifndef LMP_COMPUTE_CHUNK_ATOM_H
#define LMP_COMPUTE_CHUNK_ATOM_H
#include "compute.h"
#include <map>
namespace LAMMPS_NS {
class ComputeChunkAtom : public Compute {
public:
int nchunk,ncoord,compress,idsflag,lockcount;
int computeflag; // 1 if this compute invokes other computes
double chunk_volume_scalar;
double *chunk_volume_vec;
double **coord;
int *ichunk,*chunkID;
ComputeChunkAtom(class LAMMPS *, int, char **);
~ComputeChunkAtom();
void init();
void setup();
void compute_peratom();
void set_arrays(int);
double memory_usage();
void lock(class Fix *, bigint, bigint);
void unlock(class Fix *);
int setup_chunks();
void compute_ichunk();
private:
int which,binflag;
int regionflag,nchunksetflag,nchunkflag,discard;
int limit,limitstyle,limitfirst;
int scaleflag;
double xscale,yscale,zscale;
int argindex;
char *cfvid;
int ndim;
int dim[3],originflag[3],nlayers[3];
int minflag[3],maxflag[3];
double origin[3],delta[3];
double offset[3],invdelta[3];
double minvalue[3],maxvalue[3];
char *idregion;
class Region *region;
class Compute *cchunk;
class Fix *fchunk;
int vchunk;
int maxvar;
double *varatom;
char *id_fix;
class FixStore *fixstore;
class Fix *lockfix; // ptr to FixAveChunk that is locking out setups
// NULL if no lock currently in place
bigint lockstart,lockstop; // timesteps for start and stop of locking
bigint invoked_setup; // last timestep setup_chunks and nchunk calculated
bigint invoked_ichunk; // last timestep ichunk values calculated
int nmax,nmaxint;
double *chunk;
int molcheck; // one-time check if all molecule atoms in chunk
int *exclude; // 1 if atom is not assigned to any chunk
std::map<int,int> *hash; // store original chunks IDs before compression
// static variable for ring communication callback to access class data
// callback functions for ring communication
static ComputeChunkAtom *cptr;
static void idring(int, char *);
void assign_chunk_ids();
void compress_chunk_ids();
void check_molecules();
int setup_bins();
void bin_volumes();
void atom2bin1d();
void atom2bin2d();
void atom2bin3d();
void readdim(int, char **, int, int);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
W: More than one compute ke/atom
It is not efficient to use compute ke/atom more than once.
*/

242
src/compute_com_chunk.cpp Normal file
View File

@ -0,0 +1,242 @@
/* ----------------------------------------------------------------------
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 "string.h"
#include "compute_com_chunk.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "compute_chunk_atom.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
enum{ONCE,NFREQ,EVERY};
/* ---------------------------------------------------------------------- */
ComputeCOMChunk::ComputeCOMChunk(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 4) error->all(FLERR,"Illegal compute com/chunk command");
array_flag = 1;
size_array_cols = 3;
size_array_rows = 0;
size_array_rows_variable = 1;
extarray = 0;
// ID of compute chunk/atom
int n = strlen(arg[3]) + 1;
idchunk = new char[n];
strcpy(idchunk,arg[3]);
init();
// chunk-based data
nchunk = 1;
maxchunk = 0;
massproc = masstotal = NULL;
com = comall = NULL;
allocate();
firstflag = massneed = 1;
}
/* ---------------------------------------------------------------------- */
ComputeCOMChunk::~ComputeCOMChunk()
{
delete [] idchunk;
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(com);
memory->destroy(comall);
}
/* ---------------------------------------------------------------------- */
void ComputeCOMChunk::init()
{
int icompute = modify->find_compute(idchunk);
if (icompute < 0)
error->all(FLERR,"Chunk/atom compute does not exist for compute com/chunk");
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
if (strcmp(cchunk->style,"chunk/atom") != 0)
error->all(FLERR,"Compute com/chunk does not use chunk/atom compute");
}
/* ---------------------------------------------------------------------- */
void ComputeCOMChunk::setup()
{
// one-time calculation of per-chunk mass
// done in setup, so that ComputeChunkAtom::setup() is already called
if (firstflag && cchunk->idsflag == ONCE) {
firstflag = massneed = 0;
compute_array();
}
}
/* ---------------------------------------------------------------------- */
void ComputeCOMChunk::compute_array()
{
int index;
double massone;
double unwrap[3];
invoked_array = update->ntimestep;
// compute chunk/atom assigns atoms to chunk IDs
// extract ichunk index vector from compute
// ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms
nchunk = cchunk->setup_chunks();
cchunk->compute_ichunk();
int *ichunk = cchunk->ichunk;
if (nchunk > maxchunk) allocate();
// zero local per-chunk values
for (int i = 0; i < nchunk; i++)
com[i][0] = com[i][1] = com[i][2] = 0.0;
if (massneed)
for (int i = 0; i < nchunk; i++) massproc[i] = 0.0;
// compute COM for each chunk
double **x = atom->x;
int *mask = atom->mask;
int *type = atom->type;
imageint *image = atom->image;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
index = ichunk[i]-1;
if (index < 0) continue;
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
domain->unmap(x[i],image[i],unwrap);
com[index][0] += unwrap[0] * massone;
com[index][1] += unwrap[1] * massone;
com[index][2] += unwrap[2] * massone;
if (massneed) massproc[index] += massone;
}
MPI_Allreduce(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world);
if (massneed)
MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world);
for (int i = 0; i < nchunk; i++) {
if (masstotal[i] > 0.0) {
comall[i][0] /= masstotal[i];
comall[i][1] /= masstotal[i];
comall[i][2] /= masstotal[i];
} else comall[i][0] = comall[i][1] = comall[i][2] = 0.0;
}
}
/* ----------------------------------------------------------------------
lock methods: called by fix ave/time
these methods insure vector/array size is locked for Nfreq epoch
by passing lock info along to compute chunk/atom
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
increment lock counter
------------------------------------------------------------------------- */
void ComputeCOMChunk::lock_enable()
{
cchunk->lockcount++;
}
/* ----------------------------------------------------------------------
decrement lock counter in compute chunk/atom, it if still exists
------------------------------------------------------------------------- */
void ComputeCOMChunk::lock_disable()
{
int icompute = modify->find_compute(idchunk);
if (icompute >= 0) {
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
cchunk->lockcount--;
}
}
/* ----------------------------------------------------------------------
calculate and return # of chunks = length of vector/array
------------------------------------------------------------------------- */
int ComputeCOMChunk::lock_length()
{
nchunk = cchunk->setup_chunks();
return nchunk;
}
/* ----------------------------------------------------------------------
set the lock from startstep to stopstep
------------------------------------------------------------------------- */
void ComputeCOMChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep)
{
cchunk->lock(fixptr,startstep,stopstep);
}
/* ----------------------------------------------------------------------
unset the lock
------------------------------------------------------------------------- */
void ComputeCOMChunk::unlock(Fix *fixptr)
{
cchunk->unlock(fixptr);
}
/* ----------------------------------------------------------------------
free and reallocate per-chunk arrays
------------------------------------------------------------------------- */
void ComputeCOMChunk::allocate()
{
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(com);
memory->destroy(comall);
size_array_rows = maxchunk = nchunk;
memory->create(massproc,maxchunk,"com/chunk:massproc");
memory->create(masstotal,maxchunk,"com/chunk:masstotal");
memory->create(com,maxchunk,3,"com/chunk:com");
memory->create(comall,maxchunk,3,"com/chunk:comall");
array = comall;
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputeCOMChunk::memory_usage()
{
double bytes = (bigint) maxchunk * 2 * sizeof(double);
bytes += (bigint) maxchunk * 2*3 * sizeof(double);
return bytes;
}

28
src/compute_com_molecule.h → src/compute_com_chunk.h
View File

@ -13,31 +13,43 @@
#ifdef COMPUTE_CLASS
ComputeStyle(com/molecule,ComputeCOMMolecule)
ComputeStyle(com/chunk,ComputeCOMChunk)
#else
#ifndef LMP_COMPUTE_COM_MOLECULE_H
#define LMP_COMPUTE_COM_MOLECULE_H
#ifndef LMP_COMPUTE_COM_CHUNK_H
#define LMP_COMPUTE_COM_CHUNK_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeCOMMolecule : public Compute {
class ComputeCOMChunk : public Compute {
public:
ComputeCOMMolecule(class LAMMPS *, int, char **);
~ComputeCOMMolecule();
ComputeCOMChunk(class LAMMPS *, int, char **);
~ComputeCOMChunk();
void init();
void setup();
void compute_array();
void lock_enable();
void lock_disable();
int lock_length();
void lock(class Fix *, bigint, bigint);
void unlock(class Fix *);
double memory_usage();
private:
int nmolecules;
tagint idlo,idhi;
int nchunk,maxchunk;
int firstflag,massneed;
char *idchunk;
class ComputeChunkAtom *cchunk;
double *massproc,*masstotal;
double **com,**comall;
void allocate();
};
}

148
src/compute_com_molecule.cpp
View File

@ -1,148 +0,0 @@
/* ----------------------------------------------------------------------
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 "compute_com_molecule.h"
#include "atom.h"
#include "update.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeCOMMolecule::ComputeCOMMolecule(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 3) error->all(FLERR,"Illegal compute com/molecule command");
if (atom->molecular == 0)
error->all(FLERR,"Compute com/molecule requires molecular atom style");
array_flag = 1;
size_array_cols = 3;
extarray = 0;
// setup molecule-based data
nmolecules = molecules_in_group(idlo,idhi);
size_array_rows = nmolecules;
memory->create(massproc,nmolecules,"com/molecule:massproc");
memory->create(masstotal,nmolecules,"com/molecule:masstotal");
memory->create(com,nmolecules,3,"com/molecule:com");
memory->create(comall,nmolecules,3,"com/molecule:comall");
array = comall;
// compute masstotal for each molecule
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *type = atom->type;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
tagint imol;
double massone;
for (int i = 0; i < nmolecules; i++) massproc[i] = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
massproc[imol] += massone;
}
MPI_Allreduce(massproc,masstotal,nmolecules,MPI_DOUBLE,MPI_SUM,world);
}
/* ---------------------------------------------------------------------- */
ComputeCOMMolecule::~ComputeCOMMolecule()
{
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(com);
memory->destroy(comall);
}
/* ---------------------------------------------------------------------- */
void ComputeCOMMolecule::init()
{
int ntmp = molecules_in_group(idlo,idhi);
if (ntmp != nmolecules)
error->all(FLERR,"Molecule count changed in compute com/molecule");
}
/* ---------------------------------------------------------------------- */
void ComputeCOMMolecule::compute_array()
{
tagint imol;
double massone;
double unwrap[3];
invoked_array = update->ntimestep;
for (int i = 0; i < nmolecules; i++)
com[i][0] = com[i][1] = com[i][2] = 0.0;
double **x = atom->x;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *type = atom->type;
imageint *image = atom->image;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
domain->unmap(x[i],image[i],unwrap);
com[imol][0] += unwrap[0] * massone;
com[imol][1] += unwrap[1] * massone;
com[imol][2] += unwrap[2] * massone;
}
MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules,
MPI_DOUBLE,MPI_SUM,world);
for (int i = 0; i < nmolecules; i++) {
comall[i][0] /= masstotal[i];
comall[i][1] /= masstotal[i];
comall[i][2] /= masstotal[i];
}
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputeCOMMolecule::memory_usage()
{
double bytes = (bigint) nmolecules * 2 * sizeof(double);
if (molmap) bytes += (idhi-idlo+1) * sizeof(int);
bytes += (bigint) nmolecules * 2*3 * sizeof(double);
return bytes;
}

357
src/compute_gyration_chunk.cpp Normal file
View File

@ -0,0 +1,357 @@
/* ----------------------------------------------------------------------
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 "math.h"
#include "string.h"
#include "compute_gyration_chunk.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "compute_chunk_atom.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeGyrationChunk::ComputeGyrationChunk(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg < 4) error->all(FLERR,"Illegal compute gyration/chunk command");
// ID of compute chunk/atom
int n = strlen(arg[6]) + 1;
idchunk = new char[n];
strcpy(idchunk,arg[6]);
init();
// optional args
tensor = 0;
int iarg = 3;
while (iarg < narg) {
if (strcmp(arg[iarg],"tensor") == 0) {
tensor = 1;
iarg++;
} else error->all(FLERR,"Illegal compute gyration/chunk command");
}
if (tensor) {
array_flag = 1;
size_array_cols = 6;
size_array_rows = 0;
size_array_rows_variable = 1;
extarray = 0;
} else {
vector_flag = 1;
size_vector = 0;
size_vector_variable = 1;
extvector = 0;
}
// chunk-based data
nchunk = 1;
maxchunk = 0;
massproc = masstotal = NULL;
com = comall = NULL;
rg = rgall = NULL;
rgt = rgtall = NULL;
allocate();
}
/* ---------------------------------------------------------------------- */
ComputeGyrationChunk::~ComputeGyrationChunk()
{
delete [] idchunk;
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(com);
memory->destroy(comall);
memory->destroy(rg);
memory->destroy(rgall);
memory->destroy(rgt);
memory->destroy(rgtall);
}
/* ---------------------------------------------------------------------- */
void ComputeGyrationChunk::init()
{
int icompute = modify->find_compute(idchunk);
if (icompute < 0)
error->all(FLERR,"Chunk/atom compute does not exist for "
"compute gyration/chunk");
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
if (strcmp(cchunk->style,"chunk/atom") != 0)
error->all(FLERR,"Compute gyration/chunk does not use chunk/atom compute");
}
/* ---------------------------------------------------------------------- */
void ComputeGyrationChunk::compute_vector()
{
int i,index;
double dx,dy,dz,massone;
double unwrap[3];
invoked_array = update->ntimestep;
com_chunk();
int *ichunk = cchunk->ichunk;
for (i = 0; i < nchunk; i++) rg[i] = 0.0;
// compute Rg for each chunk
double **x = atom->x;
int *mask = atom->mask;
int *type = atom->type;
imageint *image = atom->image;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
index = ichunk[i]-1;
if (index < 0) continue;
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - comall[index][0];
dy = unwrap[1] - comall[index][1];
dz = unwrap[2] - comall[index][2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
rg[index] += (dx*dx + dy*dy + dz*dz) * massone;
}
MPI_Allreduce(rg,rgall,nchunk,MPI_DOUBLE,MPI_SUM,world);
for (int i = 0; i < nchunk; i++)
rgall[i] = sqrt(rgall[i]/masstotal[i]);
}
/* ---------------------------------------------------------------------- */
void ComputeGyrationChunk::compute_array()
{
int i,j,index;
double dx,dy,dz,massone;
double unwrap[3];
invoked_array = update->ntimestep;
com_chunk();
int *ichunk = cchunk->ichunk;
for (i = 0; i < nchunk; i++)
for (j = 0; j < 6; j++) rgt[i][j] = 0.0;
double **x = atom->x;
int *mask = atom->mask;
int *type = atom->type;
imageint *image = atom->image;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
index = ichunk[i]-1;
if (index < 0) continue;
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - comall[index][0];
dy = unwrap[1] - comall[index][1];
dz = unwrap[2] - comall[index][2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
rgt[index][0] += dx*dx * massone;
rgt[index][1] += dy*dy * massone;
rgt[index][2] += dz*dz * massone;
rgt[index][3] += dx*dy * massone;
rgt[index][4] += dx*dz * massone;
rgt[index][5] += dy*dz * massone;
}
if (nchunk)
MPI_Allreduce(&rgt[0][0],&rgtall[0][0],nchunk*6,MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < nchunk; i++)
for (j = 0; j < 6; j++)
rgtall[i][j] = rgtall[i][j]/masstotal[i];
}
/* ----------------------------------------------------------------------
calculate per-chunk COM, used by both scalar and tensor
------------------------------------------------------------------------- */
void ComputeGyrationChunk::com_chunk()
{
int index;
double massone;
double unwrap[3];
// compute chunk/atom assigns atoms to chunk IDs
// extract ichunk index vector from compute
// ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms
nchunk = cchunk->setup_chunks();
cchunk->compute_ichunk();
int *ichunk = cchunk->ichunk;
if (nchunk > maxchunk) allocate();
// zero local per-chunk values
for (int i = 0; i < nchunk; i++) {
massproc[i] = 0.0;
com[i][0] = com[i][1] = com[i][2] = 0.0;
}
// compute COM for each chunk
double **x = atom->x;
int *mask = atom->mask;
int *type = atom->type;
imageint *image = atom->image;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
index = ichunk[i]-1;
if (index < 0) continue;
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
domain->unmap(x[i],image[i],unwrap);
massproc[index] += massone;
com[index][0] += unwrap[0] * massone;
com[index][1] += unwrap[1] * massone;
com[index][2] += unwrap[2] * massone;
}
MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world);
for (int i = 0; i < nchunk; i++) {
comall[i][0] /= masstotal[i];
comall[i][1] /= masstotal[i];
comall[i][2] /= masstotal[i];
}
}
/* ----------------------------------------------------------------------
lock methods: called by fix ave/time
these methods insure vector/array size is locked for Nfreq epoch
by passing lock info along to compute chunk/atom
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
increment lock counter
------------------------------------------------------------------------- */
void ComputeGyrationChunk::lock_enable()
{
cchunk->lockcount++;
}
/* ----------------------------------------------------------------------
decrement lock counter in compute chunk/atom, it if still exists
------------------------------------------------------------------------- */
void ComputeGyrationChunk::lock_disable()
{
int icompute = modify->find_compute(idchunk);
if (icompute >= 0) {
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
cchunk->lockcount--;
}
}
/* ----------------------------------------------------------------------
calculate and return # of chunks = length of vector/array
------------------------------------------------------------------------- */
int ComputeGyrationChunk::lock_length()
{
nchunk = cchunk->setup_chunks();
return nchunk;
}
/* ----------------------------------------------------------------------
set the lock from startstep to stopstep
------------------------------------------------------------------------- */
void ComputeGyrationChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep)
{
cchunk->lock(fixptr,startstep,stopstep);
}
/* ----------------------------------------------------------------------
unset the lock
------------------------------------------------------------------------- */
void ComputeGyrationChunk::unlock(Fix *fixptr)
{
cchunk->unlock(fixptr);
}
/* ----------------------------------------------------------------------
free and reallocate per-chunk arrays
------------------------------------------------------------------------- */
void ComputeGyrationChunk::allocate()
{
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(com);
memory->destroy(comall);
memory->destroy(rg);
memory->destroy(rgall);
memory->destroy(rgt);
memory->destroy(rgtall);
size_array_rows = maxchunk = nchunk;
memory->create(massproc,maxchunk,"gyration/chunk:massproc");
memory->create(masstotal,maxchunk,"gyration/chunk:masstotal");
memory->create(com,maxchunk,3,"gyration/chunk:com");
memory->create(comall,maxchunk,3,"gyration/chunk:comall");
if (tensor) {
memory->create(rgt,maxchunk,6,"gyration/chunk:rgt");
memory->create(rgtall,maxchunk,6,"gyration/chunk:rgtall");
array = rgtall;
} else {
memory->create(rg,maxchunk,"gyration/chunk:rg");
memory->create(rgall,maxchunk,"gyration/chunk:rgall");
vector = rgall;
}
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputeGyrationChunk::memory_usage()
{
double bytes = (bigint) maxchunk * 2 * sizeof(double);
bytes += (bigint) maxchunk * 2*3 * sizeof(double);
if (tensor) bytes += (bigint) maxchunk * 2*6 * sizeof(double);
else bytes += (bigint) maxchunk * 2 * sizeof(double);
return bytes;
}

32
src/compute_gyration_molecule.h → src/compute_gyration_chunk.h
View File

@ -13,37 +13,47 @@
#ifdef COMPUTE_CLASS
ComputeStyle(gyration/molecule,ComputeGyrationMolecule)
ComputeStyle(gyration/chunk,ComputeGyrationChunk)
#else
#ifndef LMP_COMPUTE_GYRATION_MOLECULE_H
#define LMP_COMPUTE_GYRATION_MOLECULE_H
#ifndef LMP_COMPUTE_GYRATION_CHUNK_H
#define LMP_COMPUTE_GYRATION_CHUNK_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeGyrationMolecule : public Compute {
class ComputeGyrationChunk : public Compute {
public:
ComputeGyrationMolecule(class LAMMPS *, int, char **);
~ComputeGyrationMolecule();
ComputeGyrationChunk(class LAMMPS *, int, char **);
~ComputeGyrationChunk();
void init();
void compute_vector();
void compute_array();
void lock_enable();
void lock_disable();
int lock_length();
void lock(class Fix *, bigint, bigint);
void unlock(class Fix *);
double memory_usage();
private:
int nchunk,maxchunk;
char *idchunk;
class ComputeChunkAtom *cchunk;
int tensor;
int nmolecules;
tagint idlo,idhi;
double *massproc,*masstotal;
double **com,**comall;
double *rg;
double **rgt;
double *rg,*rgall;
double **rgt,**rgtall;
void molcom();
void com_chunk();
void allocate();
};
}

275
src/compute_gyration_molecule.cpp
View File

@ -1,275 +0,0 @@
/* ----------------------------------------------------------------------
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 "math.h"
#include "string.h"
#include "compute_gyration_molecule.h"
#include "atom.h"
#include "update.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeGyrationMolecule::ComputeGyrationMolecule(LAMMPS *lmp,
int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg < 3) error->all(FLERR,"Illegal compute gyration/molecule command");
if (atom->molecular == 0)
error->all(FLERR,"Compute gyration/molecule requires molecular atom style");
tensor = 0;
int iarg = 3;
while (iarg < narg) {
if (strcmp(arg[iarg],"tensor") == 0) {
tensor = 1;
iarg++;
} else error->all(FLERR,"Illegal compute gyration/molecule command");
}
// setup molecule-based data
nmolecules = molecules_in_group(idlo,idhi);
memory->create(massproc,nmolecules,"gyration/molecule:massproc");
memory->create(masstotal,nmolecules,"gyration/molecule:masstotal");
memory->create(com,nmolecules,3,"gyration/molecule:com");
memory->create(comall,nmolecules,3,"gyration/molecule:comall");
rg = vector = NULL;
rgt = array = NULL;
if (tensor) {
memory->create(rgt,nmolecules,6,"gyration/molecule:rgt");
memory->create(array,nmolecules,6,"gyration/molecule:array");
array_flag = 1;
size_array_rows = nmolecules;
size_array_cols = 6;
extarray = 0;
} else {
memory->create(rg,nmolecules,"gyration/molecule:rg");
memory->create(vector,nmolecules,"gyration/molecule:vector");
vector_flag = 1;
size_vector = nmolecules;
extvector = 0;
}
// compute masstotal for each molecule
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *type = atom->type;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
tagint imol;
double massone;
for (int i = 0; i < nmolecules; i++) massproc[i] = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
massproc[imol] += massone;
}
MPI_Allreduce(massproc,masstotal,nmolecules,MPI_DOUBLE,MPI_SUM,world);
}
/* ---------------------------------------------------------------------- */
ComputeGyrationMolecule::~ComputeGyrationMolecule()
{
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(com);
memory->destroy(comall);
memory->destroy(rg);
memory->destroy(rgt);
}
/* ---------------------------------------------------------------------- */
void ComputeGyrationMolecule::init()
{
int ntmp = molecules_in_group(idlo,idhi);
if (ntmp != nmolecules)
error->all(FLERR,"Molecule count changed in compute gyration/molecule");
}
/* ---------------------------------------------------------------------- */
void ComputeGyrationMolecule::compute_vector()
{
tagint imol;
double dx,dy,dz,massone;
double unwrap[3];
invoked_array = update->ntimestep;
molcom();
for (int i = 0; i < nmolecules; i++) rg[i] = 0.0;
double **x = atom->x;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *type = atom->type;
imageint *image = atom->image;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - comall[imol][0];
dy = unwrap[1] - comall[imol][1];
dz = unwrap[2] - comall[imol][2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
rg[imol] += (dx*dx + dy*dy + dz*dz) * massone;
}
MPI_Allreduce(rg,vector,nmolecules,MPI_DOUBLE,MPI_SUM,world);
for (int i = 0; i < nmolecules; i++)
vector[i] = sqrt(vector[i]/masstotal[i]);
}
/* ---------------------------------------------------------------------- */
void ComputeGyrationMolecule::compute_array()
{
int i,j;
tagint imol;
double dx,dy,dz,massone;
double unwrap[3];
invoked_array = update->ntimestep;
molcom();
for (i = 0; i < nmolecules; i++)
for (j = 0; j < 6; j++)
rgt[i][j] = 0.0;
double **x = atom->x;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *type = atom->type;
imageint *image = atom->image;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - comall[imol][0];
dy = unwrap[1] - comall[imol][1];
dz = unwrap[2] - comall[imol][2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
rgt[imol][0] += dx*dx * massone;
rgt[imol][1] += dy*dy * massone;
rgt[imol][2] += dz*dz * massone;
rgt[imol][3] += dx*dy * massone;
rgt[imol][4] += dx*dz * massone;
rgt[imol][5] += dy*dz * massone;
}
if (nmolecules)
MPI_Allreduce(&rgt[0][0],&array[0][0],nmolecules*6,
MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < nmolecules; i++)
for (j = 0; j < 6; j++)
array[i][j] = array[i][j]/masstotal[i];
}
/* ----------------------------------------------------------------------
calculate per-molecule COM
------------------------------------------------------------------------- */
void ComputeGyrationMolecule::molcom()
{
tagint imol;
double massone;
double unwrap[3];
for (int i = 0; i < nmolecules; i++)
com[i][0] = com[i][1] = com[i][2] = 0.0;
double **x = atom->x;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *type = atom->type;
imageint *image = atom->image;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
domain->unmap(x[i],image[i],unwrap);
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
com[imol][0] += unwrap[0] * massone;
com[imol][1] += unwrap[1] * massone;
com[imol][2] += unwrap[2] * massone;
}
MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules,
MPI_DOUBLE,MPI_SUM,world);
for (int i = 0; i < nmolecules; i++) {
comall[i][0] /= masstotal[i];
comall[i][1] /= masstotal[i];
comall[i][2] /= masstotal[i];
}
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputeGyrationMolecule::memory_usage()
{
double bytes = (bigint) nmolecules * 2 * sizeof(double);
if (molmap) bytes += (idhi-idlo+1) * sizeof(int);
bytes += (bigint) nmolecules * 2*3 * sizeof(double);
if (tensor) bytes += (bigint) nmolecules * 2*6 * sizeof(double);
else bytes += (bigint) nmolecules * 2 * sizeof(double);
return bytes;
}

256
src/compute_inertia_chunk.cpp Normal file
View File

@ -0,0 +1,256 @@
/* ----------------------------------------------------------------------
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 "string.h"
#include "compute_inertia_chunk.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "compute_chunk_atom.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeInertiaChunk::ComputeInertiaChunk(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 4) error->all(FLERR,"Illegal compute inertia/chunk command");
array_flag = 1;
size_array_cols = 6;
size_array_rows = 0;
size_array_rows_variable = 1;
extarray = 0;
// ID of compute chunk/atom
int n = strlen(arg[6]) + 1;
idchunk = new char[n];
strcpy(idchunk,arg[6]);
init();
// chunk-based data
nchunk = 1;
maxchunk = 0;
massproc = masstotal = NULL;
com = comall = NULL;
inertia = inertiaall = NULL;
allocate();
}
/* ---------------------------------------------------------------------- */
ComputeInertiaChunk::~ComputeInertiaChunk()
{
delete [] idchunk;
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(com);
memory->destroy(comall);
memory->destroy(inertia);
memory->destroy(inertiaall);
}
/* ---------------------------------------------------------------------- */
void ComputeInertiaChunk::init()
{
int icompute = modify->find_compute(idchunk);
if (icompute < 0)
error->all(FLERR,"Chunk/atom compute does not exist for "
"compute inertia/chunk");
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
if (strcmp(cchunk->style,"chunk/atom") != 0)
error->all(FLERR,"Compute inertia/chunk does not use chunk/atom compute");
}
/* ---------------------------------------------------------------------- */
void ComputeInertiaChunk::compute_array()
{
int i,j,index;
double dx,dy,dz,massone;
double unwrap[3];
invoked_array = update->ntimestep;
// compute chunk/atom assigns atoms to chunk IDs
// extract ichunk index vector from compute
// ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms
nchunk = cchunk->setup_chunks();
cchunk->compute_ichunk();
int *ichunk = cchunk->ichunk;
if (nchunk > maxchunk) allocate();
// zero local per-chunk values
for (int i = 0; i < nchunk; i++) {
massproc[i] = 0.0;
com[i][0] = com[i][1] = com[i][2] = 0.0;
for (j = 0; j < 6; j++) inertia[i][j] = 0.0;
}
// compute COM for each chunk
double **x = atom->x;
int *mask = atom->mask;
int *type = atom->type;
imageint *image = atom->image;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
index = ichunk[i]-1;
if (index < 0) continue;
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
domain->unmap(x[i],image[i],unwrap);
massproc[index] += massone;
com[index][0] += unwrap[0] * massone;
com[index][1] += unwrap[1] * massone;
com[index][2] += unwrap[2] * massone;
}
MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world);
for (int i = 0; i < nchunk; i++) {
comall[i][0] /= masstotal[i];
comall[i][1] /= masstotal[i];
comall[i][2] /= masstotal[i];
}
// compute inertia tensor for each chunk
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
index = ichunk[i]-1;
if (index < 0) continue;
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - comall[index][0];
dy = unwrap[1] - comall[index][1];
dz = unwrap[2] - comall[index][2];
inertia[index][0] += massone * (dy*dy + dz*dz);
inertia[index][1] += massone * (dx*dx + dz*dz);
inertia[index][2] += massone * (dx*dx + dy*dy);
inertia[index][3] -= massone * dx*dy;
inertia[index][4] -= massone * dy*dz;
inertia[index][5] -= massone * dx*dz;
}
MPI_Allreduce(&inertia[0][0],&inertiaall[0][0],6*nchunk,
MPI_DOUBLE,MPI_SUM,world);
}
/* ----------------------------------------------------------------------
lock methods: called by fix ave/time
these methods insure vector/array size is locked for Nfreq epoch
by passing lock info along to compute chunk/atom
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
increment lock counter
------------------------------------------------------------------------- */
void ComputeInertiaChunk::lock_enable()
{
cchunk->lockcount++;
}
/* ----------------------------------------------------------------------
decrement lock counter in compute chunk/atom, it if still exists
------------------------------------------------------------------------- */
void ComputeInertiaChunk::lock_disable()
{
int icompute = modify->find_compute(idchunk);
if (icompute >= 0) {
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
cchunk->lockcount--;
}
}
/* ----------------------------------------------------------------------
calculate and return # of chunks = length of vector/array
------------------------------------------------------------------------- */
int ComputeInertiaChunk::lock_length()
{
nchunk = cchunk->setup_chunks();
return nchunk;
}
/* ----------------------------------------------------------------------
set the lock from startstep to stopstep
------------------------------------------------------------------------- */
void ComputeInertiaChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep)
{
cchunk->lock(fixptr,startstep,stopstep);
}
/* ----------------------------------------------------------------------
unset the lock
------------------------------------------------------------------------- */
void ComputeInertiaChunk::unlock(Fix *fixptr)
{
cchunk->unlock(fixptr);
}
/* ----------------------------------------------------------------------
free and reallocate per-chunk arrays
------------------------------------------------------------------------- */
void ComputeInertiaChunk::allocate()
{
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(com);
memory->destroy(comall);
memory->destroy(inertia);
memory->destroy(inertiaall);
size_array_rows = maxchunk = nchunk;
memory->create(massproc,maxchunk,"inertia/chunk:massproc");
memory->create(masstotal,maxchunk,"inertia/chunk:masstotal");
memory->create(com,maxchunk,3,"inertia/chunk:com");
memory->create(comall,maxchunk,3,"inertia/chunk:comall");
memory->create(inertia,maxchunk,6,"inertia/chunk:inertia");
memory->create(inertiaall,maxchunk,6,"inertia/chunk:inertiaall");
array = inertiaall;
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputeInertiaChunk::memory_usage()
{
double bytes = (bigint) maxchunk * 2 * sizeof(double);
bytes += (bigint) maxchunk * 2*3 * sizeof(double);
bytes += (bigint) maxchunk * 2*6 * sizeof(double);
return bytes;
}

26
src/compute_inertia_molecule.h → src/compute_inertia_chunk.h
View File

@ -13,32 +13,42 @@
#ifdef COMPUTE_CLASS
ComputeStyle(inertia/molecule,ComputeInertiaMolecule)
ComputeStyle(inertia/chunk,ComputeInertiaChunk)
#else
#ifndef LMP_COMPUTE_INERTIA_MOLECULE_H
#define LMP_COMPUTE_INERTIA_MOLECULE_H
#ifndef LMP_COMPUTE_INERTIA_CHUNK_H
#define LMP_COMPUTE_INERTIA_CHUNK_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeInertiaMolecule : public Compute {
class ComputeInertiaChunk : public Compute {
public:
ComputeInertiaMolecule(class LAMMPS *, int, char **);
~ComputeInertiaMolecule();
ComputeInertiaChunk(class LAMMPS *, int, char **);
~ComputeInertiaChunk();
void init();
void compute_array();
void lock_enable();
void lock_disable();
int lock_length();
void lock(class Fix *, bigint, bigint);
void unlock(class Fix *);
double memory_usage();
private:
int nmolecules;
tagint idlo,idhi;
int nchunk,maxchunk;
char *idchunk;
class ComputeChunkAtom *cchunk;
double *massproc,*masstotal;
double **com,**comall;
double **inertia,**inertiaall;
void allocate();
};
}

185
src/compute_inertia_molecule.cpp
View File

@ -1,185 +0,0 @@
/* ----------------------------------------------------------------------
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 "compute_inertia_molecule.h"
#include "atom.h"
#include "update.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeInertiaMolecule::
ComputeInertiaMolecule(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 3) error->all(FLERR,"Illegal compute inertia/molecule command");
if (atom->molecular == 0)
error->all(FLERR,"Compute inertia/molecule requires molecular atom style");
array_flag = 1;
size_array_cols = 6;
extarray = 0;
// setup molecule-based data
nmolecules = molecules_in_group(idlo,idhi);
size_array_rows = nmolecules;
memory->create(massproc,nmolecules,"inertia/molecule:massproc");
memory->create(masstotal,nmolecules,"inertia/molecule:masstotal");
memory->create(com,nmolecules,3,"inertia/molecule:com");
memory->create(comall,nmolecules,3,"inertia/molecule:comall");
memory->create(inertia,nmolecules,6,"inertia/molecule:inertia");
memory->create(inertiaall,nmolecules,6,"inertia/molecule:inertiaall");
array = inertiaall;
// compute masstotal for each molecule
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *type = atom->type;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
tagint imol;
double massone;
for (int i = 0; i < nmolecules; i++) massproc[i] = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
massproc[imol] += massone;
}
MPI_Allreduce(massproc,masstotal,nmolecules,MPI_DOUBLE,MPI_SUM,world);
}
/* ---------------------------------------------------------------------- */
ComputeInertiaMolecule::~ComputeInertiaMolecule()
{
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(com);
memory->destroy(comall);
memory->destroy(inertia);
memory->destroy(inertiaall);
}
/* ---------------------------------------------------------------------- */
void ComputeInertiaMolecule::init()
{
int ntmp = molecules_in_group(idlo,idhi);
if (ntmp != nmolecules)
error->all(FLERR,"Molecule count changed in compute inertia/molecule");
}
/* ---------------------------------------------------------------------- */
void ComputeInertiaMolecule::compute_array()
{
int i,j;
tagint imol;
double dx,dy,dz,massone;
double unwrap[3];
invoked_array = update->ntimestep;
double **x = atom->x;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *type = atom->type;
imageint *image = atom->image;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
// center-of-mass for each molecule
for (i = 0; i < nmolecules; i++)
com[i][0] = com[i][1] = com[i][2] = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
domain->unmap(x[i],image[i],unwrap);
com[imol][0] += unwrap[0] * massone;
com[imol][1] += unwrap[1] * massone;
com[imol][2] += unwrap[2] * massone;
}
MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules,
MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < nmolecules; i++) {
comall[i][0] /= masstotal[i];
comall[i][1] /= masstotal[i];
comall[i][2] /= masstotal[i];
}
// inertia tensor for each molecule
for (i = 0; i < nmolecules; i++)
for (j = 0; j < 6; j++)
inertia[i][j] = 0.0;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - comall[imol][0];
dy = unwrap[1] - comall[imol][1];
dz = unwrap[2] - comall[imol][2];
inertia[imol][0] += massone * (dy*dy + dz*dz);
inertia[imol][1] += massone * (dx*dx + dz*dz);
inertia[imol][2] += massone * (dx*dx + dy*dy);
inertia[imol][3] -= massone * dx*dy;
inertia[imol][4] -= massone * dy*dz;
inertia[imol][5] -= massone * dx*dz;
}
MPI_Allreduce(&inertia[0][0],&inertiaall[0][0],6*nmolecules,
MPI_DOUBLE,MPI_SUM,world);
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputeInertiaMolecule::memory_usage()
{
double bytes = (bigint) nmolecules * 2 * sizeof(double);
if (molmap) bytes += (idhi-idlo+1) * sizeof(int);
bytes += (bigint) nmolecules * 2*3 * sizeof(double);
bytes += (bigint) nmolecules * 2*6 * sizeof(double);
return bytes;
}

261
src/compute_msd_chunk.cpp Normal file
View File

@ -0,0 +1,261 @@
/* ----------------------------------------------------------------------
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 "string.h"
#include "compute_msd_chunk.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "compute_chunk_atom.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeMSDChunk::ComputeMSDChunk(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 4) error->all(FLERR,"Illegal compute msd/chunk command");
array_flag = 1;
size_array_cols = 4;
size_array_rows = 0;
size_array_rows_variable = 1;
extarray = 0;
// ID of compute chunk/atom
int n = strlen(arg[6]) + 1;
idchunk = new char[n];
strcpy(idchunk,arg[6]);
init();
massproc = masstotal = NULL;
com = comall = cominit = NULL;
msd = NULL;
firstflag = 1;
}
/* ---------------------------------------------------------------------- */
ComputeMSDChunk::~ComputeMSDChunk()
{
delete [] idchunk;
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(com);
memory->destroy(comall);
memory->destroy(cominit);
memory->destroy(msd);
}
/* ---------------------------------------------------------------------- */
void ComputeMSDChunk::init()
{
int icompute = modify->find_compute(idchunk);
if (icompute < 0)
error->all(FLERR,"Chunk/atom compute does not exist for compute msd/chunk");
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
if (strcmp(cchunk->style,"chunk/atom") != 0)
error->all(FLERR,"Compute msd/chunk does not use chunk/atom compute");
}
/* ----------------------------------------------------------------------
compute initial COM for each chunk
only once on timestep compute is defined, when firstflag = 1
------------------------------------------------------------------------- */
void ComputeMSDChunk::setup()
{
if (!firstflag) return;
firstflag = 0;
compute_array();
for (int i = 0; i < nchunk; i++) {
cominit[i][0] = comall[i][0];
cominit[i][1] = comall[i][1];
cominit[i][2] = comall[i][2];
}
}
/* ---------------------------------------------------------------------- */
void ComputeMSDChunk::compute_array()
{
int index;
double massone;
double unwrap[3];
invoked_array = update->ntimestep;
// compute chunk/atom assigns atoms to chunk IDs
// extract ichunk index vector from compute
// ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms
int n = cchunk->setup_chunks();
cchunk->compute_ichunk();
int *ichunk = cchunk->ichunk;
// first time call, allocate per-chunk arrays
// thereafter, require nchunk remain the same
if (firstflag) allocate(n);
else if (n != nchunk)
error->all(FLERR,"Compute msd/chunk nchunk is not static");
// zero local per-chunk values
for (int i = 0; i < nchunk; i++) {
massproc[i] = 0.0;
com[i][0] = com[i][1] = com[i][2] = 0.0;
}
// compute current COM for each chunk
double **x = atom->x;
int *mask = atom->mask;
int *type = atom->type;
imageint *image = atom->image;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
index = ichunk[i]-1;
if (index < 0) continue;
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
domain->unmap(x[i],image[i],unwrap);
massproc[index] += massone;
com[index][0] += unwrap[0] * massone;
com[index][1] += unwrap[1] * massone;
com[index][2] += unwrap[2] * massone;
}
MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world);
for (int i = 0; i < nchunk; i++) {
comall[i][0] /= masstotal[i];
comall[i][1] /= masstotal[i];
comall[i][2] /= masstotal[i];
}
// MSD is difference between current and initial COM
// cominit does not yet exist when called from constructor
if (firstflag) return;
double dx,dy,dz;
for (int i = 0; i < nchunk; i++) {
dx = comall[i][0] - cominit[i][0];
dy = comall[i][1] - cominit[i][1];
dz = comall[i][2] - cominit[i][2];
msd[i][0] = dx*dx;
msd[i][1] = dy*dy;
msd[i][2] = dz*dz;
msd[i][3] = dx*dx + dy*dy + dz*dz;
}
}
/* ----------------------------------------------------------------------
lock methods: called by fix ave/time
these methods insure vector/array size is locked for Nfreq epoch
by passing lock info along to compute chunk/atom
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
increment lock counter
------------------------------------------------------------------------- */
void ComputeMSDChunk::lock_enable()
{
cchunk->lockcount++;
}
/* ----------------------------------------------------------------------
decrement lock counter in compute chunk/atom, it if still exists
------------------------------------------------------------------------- */
void ComputeMSDChunk::lock_disable()
{
int icompute = modify->find_compute(idchunk);
if (icompute >= 0) {
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
cchunk->lockcount--;
}
}
/* ----------------------------------------------------------------------
calculate and return # of chunks = length of vector/array
------------------------------------------------------------------------- */
int ComputeMSDChunk::lock_length()
{
nchunk = cchunk->setup_chunks();
return nchunk;
}
/* ----------------------------------------------------------------------
set the lock from startstep to stopstep
------------------------------------------------------------------------- */
void ComputeMSDChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep)
{
cchunk->lock(fixptr,startstep,stopstep);
}
/* ----------------------------------------------------------------------
unset the lock
------------------------------------------------------------------------- */
void ComputeMSDChunk::unlock(Fix *fixptr)
{
cchunk->unlock(fixptr);
}
/* ----------------------------------------------------------------------
one-time allocate of per-chunk arrays
------------------------------------------------------------------------- */
void ComputeMSDChunk::allocate(int n)
{
size_array_rows = nchunk = n;
memory->create(massproc,nchunk,"msd/chunk:massproc");
memory->create(masstotal,nchunk,"msd/chunk:masstotal");
memory->create(com,nchunk,3,"msd/chunk:com");
memory->create(comall,nchunk,3,"msd/chunk:comall");
memory->create(cominit,nchunk,3,"msd/chunk:cominit");
memory->create(msd,nchunk,4,"msd/chunk:msd");
array = msd;
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputeMSDChunk::memory_usage()
{
double bytes = (bigint) nchunk * 2 * sizeof(double);
bytes += (bigint) nchunk * 3*3 * sizeof(double);
bytes += (bigint) nchunk * 4 * sizeof(double);
return bytes;
}

34
src/compute_msd_molecule.h → src/compute_msd_chunk.h
View File

@ -13,33 +13,45 @@
#ifdef COMPUTE_CLASS
ComputeStyle(msd/molecule,ComputeMSDMolecule)
ComputeStyle(msd/chunk,ComputeMSDChunk)
#else
#ifndef LMP_COMPUTE_MSD_MOLECULE_H
#define LMP_COMPUTE_MSD_MOLECULE_H
#ifndef LMP_COMPUTE_MSD_CHUNK_H
#define LMP_COMPUTE_MSD_CHUNK_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeMSDMolecule : public Compute {
class ComputeMSDChunk : public Compute {
public:
ComputeMSDMolecule(class LAMMPS *, int, char **);
~ComputeMSDMolecule();
ComputeMSDChunk(class LAMMPS *, int, char **);
~ComputeMSDChunk();
void init();
void setup();
void compute_array();
void lock_enable();
void lock_disable();
int lock_length();
void lock(class Fix *, bigint, bigint);
void unlock(class Fix *);
double memory_usage();
private:
int nmolecules;
tagint idlo,idhi;
int firstflag;
int nchunk;
char *idchunk;
class ComputeChunkAtom *cchunk;
double *massproc,*masstotal;
double **com,**comall,**cominit;
double **msd;
int firstflag;
void allocate(int);
};
}
@ -55,11 +67,11 @@ Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Compute msd/molecule requires molecular atom style
E: Compute com/molecule requires molecular atom style
Self-explanatory.
E: Molecule count changed in compute msd/molecule
E: Molecule count changed in compute com/molecule
Number of molecules must remain constant over time.

181
src/compute_msd_molecule.cpp
View File

@ -1,181 +0,0 @@
/* ----------------------------------------------------------------------
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 "compute_msd_molecule.h"
#include "atom.h"
#include "update.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeMSDMolecule::ComputeMSDMolecule(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 3) error->all(FLERR,"Illegal compute msd/molecule command");
if (atom->molecular == 0)
error->all(FLERR,"Compute msd/molecule requires molecular atom style");
array_flag = 1;
size_array_cols = 4;
extarray = 0;
// setup molecule-based data and initial COM positions
nmolecules = molecules_in_group(idlo,idhi);
size_array_rows = nmolecules;
memory->create(massproc,nmolecules,"msd/molecule:massproc");
memory->create(masstotal,nmolecules,"msd/molecule:masstotal");
memory->create(com,nmolecules,3,"msd/molecule:com");
memory->create(comall,nmolecules,3,"msd/molecule:comall");
memory->create(cominit,nmolecules,3,"msd/molecule:cominit");
memory->create(msd,nmolecules,4,"msd/molecule:msd");
array = msd;
// compute masstotal for each molecule
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *type = atom->type;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
tagint imol;
double massone;
for (int i = 0; i < nmolecules; i++) massproc[i] = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
massproc[imol] += massone;
}
MPI_Allreduce(massproc,masstotal,nmolecules,MPI_DOUBLE,MPI_SUM,world);
// compute initial COM for each molecule
firstflag = 1;
compute_array();
for (int i = 0; i < nmolecules; i++) {
cominit[i][0] = comall[i][0];
cominit[i][1] = comall[i][1];
cominit[i][2] = comall[i][2];
}
firstflag = 0;
}
/* ---------------------------------------------------------------------- */
ComputeMSDMolecule::~ComputeMSDMolecule()
{
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(com);
memory->destroy(comall);
memory->destroy(cominit);
memory->destroy(msd);
}
/* ---------------------------------------------------------------------- */
void ComputeMSDMolecule::init()
{
int ntmp = molecules_in_group(idlo,idhi);
if (ntmp != nmolecules)
error->all(FLERR,"Molecule count changed in compute msd/molecule");
}
/* ---------------------------------------------------------------------- */
void ComputeMSDMolecule::compute_array()
{
tagint imol;
double dx,dy,dz,massone;
double unwrap[3];
invoked_array = update->ntimestep;
// compute current COM positions
for (int i = 0; i < nmolecules; i++)
com[i][0] = com[i][1] = com[i][2] = 0.0;
double **x = atom->x;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *type = atom->type;
imageint *image = atom->image;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
domain->unmap(x[i],image[i],unwrap);
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
com[imol][0] += unwrap[0] * massone;
com[imol][1] += unwrap[1] * massone;
com[imol][2] += unwrap[2] * massone;
}
MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules,
MPI_DOUBLE,MPI_SUM,world);
for (int i = 0; i < nmolecules; i++) {
comall[i][0] /= masstotal[i];
comall[i][1] /= masstotal[i];
comall[i][2] /= masstotal[i];
}
// MSD is difference between current and initial COM
// cominit does not yet exist when called from constructor
if (firstflag) return;
for (int i = 0; i < nmolecules; i++) {
dx = comall[i][0] - cominit[i][0];
dy = comall[i][1] - cominit[i][1];
dz = comall[i][2] - cominit[i][2];
msd[i][0] = dx*dx;
msd[i][1] = dy*dy;
msd[i][2] = dz*dz;
msd[i][3] = dx*dx + dy*dy + dz*dz;
}
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputeMSDMolecule::memory_usage()
{
double bytes = (bigint) nmolecules * 2 * sizeof(double);
if (molmap) bytes += (idhi-idlo+1) * sizeof(int);
bytes += (bigint) nmolecules * 2*3 * sizeof(double);
bytes += (bigint) nmolecules * 4 * sizeof(double);
return bytes;
}

335
src/compute_property_chunk.cpp Normal file
View File

@ -0,0 +1,335 @@
/* ----------------------------------------------------------------------
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 "string.h"
#include "compute_property_chunk.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "compute_chunk_atom.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputePropertyChunk::ComputePropertyChunk(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg < 5) error->all(FLERR,"Illegal compute property/chunk command");
// ID of compute chunk/atom
int n = strlen(arg[3]) + 1;
idchunk = new char[n];
strcpy(idchunk,arg[3]);
init();
// parse values
nvalues = narg - 4;
pack_choice = new FnPtrPack[nvalues];
countflag = 0;
int i;
for (int iarg = 4; iarg < narg; iarg++) {
i = iarg-4;
if (strcmp(arg[iarg],"count") == 0) {
pack_choice[i] = &ComputePropertyChunk::pack_count;
countflag = 1;
} else if (strcmp(arg[iarg],"id") == 0) {
if (!cchunk->compress)
error->all(FLERR,"Compute chunk/atom stores no IDs for "
"compute property/chunk");
pack_choice[i] = &ComputePropertyChunk::pack_id;
} else if (strcmp(arg[iarg],"coord1") == 0) {
if (cchunk->ncoord < 1)
error->all(FLERR,"Compute chunk/atom stores no coord1 for "
"compute property/chunk");
pack_choice[i] = &ComputePropertyChunk::pack_coord1;
} else if (strcmp(arg[iarg],"coord2") == 0) {
if (cchunk->ncoord < 2)
error->all(FLERR,"Compute chunk/atom stores no coord2 for "
"compute property/chunk");
pack_choice[i] = &ComputePropertyChunk::pack_coord2;
} else if (strcmp(arg[iarg],"coord3") == 0) {
if (cchunk->ncoord < 3)
error->all(FLERR,"Compute chunk/atom stores no coord3 for "
"compute property/chunk");
pack_choice[i] = &ComputePropertyChunk::pack_coord3;
} else error->all(FLERR,
"Invalid keyword in compute property/chunk command");
}
// initialization
nchunk = 1;
maxchunk = 0;
vector = NULL;
array = NULL;
count_one = count_all = NULL;
allocate();
if (nvalues == 1) {
vector_flag = 1;
size_vector = 0;
size_vector_variable = 1;
extvector = 0;
} else {
array_flag = 1;
size_array_cols = nvalues;
size_array_rows = 0;
size_array_rows_variable = 1;
extarray = 0;
}
}
/* ---------------------------------------------------------------------- */
ComputePropertyChunk::~ComputePropertyChunk()
{
delete [] idchunk;
delete [] pack_choice;
memory->destroy(vector);
memory->destroy(array);
memory->destroy(count_one);
memory->destroy(count_all);
}
/* ---------------------------------------------------------------------- */
void ComputePropertyChunk::init()
{
int icompute = modify->find_compute(idchunk);
if (icompute < 0)
error->all(FLERR,"Chunk/atom compute does not exist for "
"compute property/chunk");
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
if (strcmp(cchunk->style,"chunk/atom") != 0)
error->all(FLERR,"Compute property/chunk does not use chunk/atom compute");
}
/* ---------------------------------------------------------------------- */
void ComputePropertyChunk::compute_vector()
{
invoked_vector = update->ntimestep;
// compute chunk/atom assigns atoms to chunk IDs
// if need count, extract ichunk index vector from compute
// ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms
nchunk = cchunk->setup_chunks();
if (nchunk > maxchunk) allocate();
if (countflag) {
cchunk->compute_ichunk();
ichunk = cchunk->ichunk;
}
// fill vector
buf = vector;
(this->*pack_choice[0])(0);
}
/* ---------------------------------------------------------------------- */
void ComputePropertyChunk::compute_array()
{
invoked_array = update->ntimestep;
// compute chunk/atom assigns atoms to chunk IDs
// if need count, extract ichunk index vector from compute
// ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms
nchunk = cchunk->setup_chunks();
if (nchunk > maxchunk) allocate();
if (countflag) {
cchunk->compute_ichunk();
ichunk = cchunk->ichunk;
}
// fill array
if (array) buf = &array[0][0];
for (int n = 0; n < nvalues; n++)
(this->*pack_choice[n])(n);
}
/* ----------------------------------------------------------------------
lock methods: called by fix ave/time
these methods insure vector/array size is locked for Nfreq epoch
by passing lock info along to compute chunk/atom
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
increment lock counter
------------------------------------------------------------------------- */
void ComputePropertyChunk::lock_enable()
{
cchunk->lockcount++;
}
/* ----------------------------------------------------------------------
decrement lock counter in compute chunk/atom, it if still exists
------------------------------------------------------------------------- */
void ComputePropertyChunk::lock_disable()
{
int icompute = modify->find_compute(idchunk);
if (icompute >= 0) {
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
cchunk->lockcount--;
}
}
/* ----------------------------------------------------------------------
calculate and return # of chunks = length of vector/array
------------------------------------------------------------------------- */
int ComputePropertyChunk::lock_length()
{
nchunk = cchunk->setup_chunks();
return nchunk;
}
/* ----------------------------------------------------------------------
set the lock from startstep to stopstep
------------------------------------------------------------------------- */
void ComputePropertyChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep)
{
cchunk->lock(fixptr,startstep,stopstep);
}
/* ----------------------------------------------------------------------
unset the lock
------------------------------------------------------------------------- */
void ComputePropertyChunk::unlock(Fix *fixptr)
{
cchunk->unlock(fixptr);
}
/* ----------------------------------------------------------------------
free and reallocate per-chunk arrays
------------------------------------------------------------------------- */
void ComputePropertyChunk::allocate()
{
memory->destroy(vector);
memory->destroy(array);
memory->destroy(count_one);
memory->destroy(count_all);
size_array_rows = maxchunk = nchunk;
if (nvalues == 1) memory->create(vector,maxchunk,"property/chunk:vector");
else memory->create(array,maxchunk,nvalues,"property/chunk:array");
if (countflag) {
memory->create(count_one,maxchunk,"property/chunk:count_one");
memory->create(count_one,maxchunk,"property/chunk:count_all");
}
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputePropertyChunk::memory_usage()
{
double bytes = (bigint) nchunk * nvalues * sizeof(double);
if (countflag) bytes += (bigint) nchunk * 2 * sizeof(int);
return bytes;
}
/* ----------------------------------------------------------------------
one method for every keyword compute property/chunk can output
the property is packed into buf starting at n with stride nvalues
customize a new keyword by adding a method
------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
void ComputePropertyChunk::pack_count(int n)
{
int index;
for (int m = 0; m < nchunk; m++) count_one[m] = 0;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
index = ichunk[i]-1;
if (index < 0) continue;
count_one[index]++;
}
}
MPI_Allreduce(count_one,count_all,nchunk,MPI_INT,MPI_SUM,world);
for (int m = 0; m < nchunk; m++) {
buf[n] = count_all[m];
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyChunk::pack_id(int n)
{
int *origID = cchunk->chunkID;
for (int m = 0; m < nchunk; m++) {
buf[n] = origID[m];
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyChunk::pack_coord1(int n)
{
double **coord = cchunk->coord;
for (int m = 0; m < nchunk; m++) {
buf[n] = coord[m][0];
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyChunk::pack_coord2(int n)
{
double **coord = cchunk->coord;
for (int m = 0; m < nchunk; m++) {
buf[n] = coord[m][1];
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyChunk::pack_coord3(int n)
{
double **coord = cchunk->coord;
for (int m = 0; m < nchunk; m++) {
buf[n] = coord[m][2];
n += nvalues;
}
}

36
src/compute_property_molecule.h → src/compute_property_chunk.h
View File

@ -13,37 +13,53 @@
#ifdef COMPUTE_CLASS
ComputeStyle(property/molecule,ComputePropertyMolecule)
ComputeStyle(property/chunk,ComputePropertyChunk)
#else
#ifndef LMP_COMPUTE_PROPERTY_MOLECULE_H
#define LMP_COMPUTE_PROPERTY_MOLECULE_H
#ifndef LMP_COMPUTE_CHUNK_MOLECULE_H
#define LMP_COMPUTE_CHUNK_MOLECULE_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputePropertyMolecule : public Compute {
class ComputePropertyChunk : public Compute {
public:
ComputePropertyMolecule(class LAMMPS *, int, char **);
~ComputePropertyMolecule();
ComputePropertyChunk(class LAMMPS *, int, char **);
~ComputePropertyChunk();
void init();
void compute_vector();
void compute_array();
void lock_enable();
void lock_disable();
int lock_length();
void lock(class Fix *, bigint, bigint);
void unlock(class Fix *);
double memory_usage();
private:
int nvalues,nmolecules;
tagint idlo,idhi;
int nchunk,maxchunk;
char *idchunk;
class ComputeChunkAtom *cchunk;
int *ichunk;
int nvalues,countflag;
double *buf;
int *count_one,*count_all;
typedef void (ComputePropertyMolecule::*FnPtrPack)(int);
void allocate();
typedef void (ComputePropertyChunk::*FnPtrPack)(int);
FnPtrPack *pack_choice; // ptrs to pack functions
void pack_mol(int);
void pack_count(int);
void pack_id(int);
void pack_coord1(int);
void pack_coord2(int);
void pack_coord3(int);
};
}

173
src/compute_property_molecule.cpp
View File

@ -1,173 +0,0 @@
/* ----------------------------------------------------------------------
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 "string.h"
#include "compute_property_molecule.h"
#include "atom.h"
#include "update.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputePropertyMolecule::
ComputePropertyMolecule(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg < 4) error->all(FLERR,"Illegal compute property/molecule command");
if (atom->molecular == 0)
error->all(FLERR,"Compute property/molecule requires molecular atom style");
nvalues = narg - 3;
pack_choice = new FnPtrPack[nvalues];
int i;
for (int iarg = 3; iarg < narg; iarg++) {
i = iarg-3;
if (strcmp(arg[iarg],"mol") == 0)
pack_choice[i] = &ComputePropertyMolecule::pack_mol;
else if (strcmp(arg[iarg],"count") == 0)
pack_choice[i] = &ComputePropertyMolecule::pack_count;
else error->all(FLERR,
"Invalid keyword in compute property/molecule command");
}
// setup molecule-based data
nmolecules = molecules_in_group(idlo,idhi);
vector = NULL;
array = NULL;
if (nvalues == 1) {
memory->create(vector,nmolecules,"property/molecule:vector");
vector_flag = 1;
size_vector = nmolecules;
extvector = 0;
} else {
memory->create(array,nmolecules,nvalues,"property/molecule:array");
array_flag = 1;
size_array_rows = nmolecules;
size_array_cols = nvalues;
extarray = 0;
}
// fill vector or array with molecule values
if (nvalues == 1) {
buf = vector;
(this->*pack_choice[0])(0);
} else {
if (array) buf = &array[0][0];
for (int n = 0; n < nvalues; n++)
(this->*pack_choice[n])(n);
}
}
/* ---------------------------------------------------------------------- */
ComputePropertyMolecule::~ComputePropertyMolecule()
{
delete [] pack_choice;
memory->destroy(vector);
memory->destroy(array);
}
/* ---------------------------------------------------------------------- */
void ComputePropertyMolecule::init()
{
int ntmp = molecules_in_group(idlo,idhi);
if (ntmp != nmolecules)
error->all(FLERR,"Molecule count changed in compute property/molecule");
}
/* ---------------------------------------------------------------------- */
void ComputePropertyMolecule::compute_vector()
{
invoked_vector = update->ntimestep;
}
/* ---------------------------------------------------------------------- */
void ComputePropertyMolecule::compute_array()
{
invoked_array = update->ntimestep;
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputePropertyMolecule::memory_usage()
{
double bytes = (bigint) nmolecules * nvalues * sizeof(double);
if (molmap) bytes += (idhi-idlo+1) * sizeof(int);
return bytes;
}
/* ----------------------------------------------------------------------
one method for every keyword compute property/molecule can output
the atom property is packed into buf starting at n with stride nvalues
customize a new keyword by adding a method
------------------------------------------------------------------------- */
void ComputePropertyMolecule::pack_mol(int n)
{
for (tagint m = idlo; m <= idhi; m++)
if (molmap == NULL || molmap[m-idlo] >= 0) {
buf[n] = m;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyMolecule::pack_count(int n)
{
int i,m;
tagint imol;
int *count_one = new int[nmolecules];
for (m = 0; m < nmolecules; m++) count_one[m] = 0;
tagint *molecule = atom->molecule;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
count_one[imol]++;
}
int *count_all = new int[nmolecules];
MPI_Allreduce(count_one,count_all,nmolecules,MPI_INT,MPI_SUM,world);
for (m = 0; m < nmolecules; m++)
if (molmap == NULL || molmap[m] >= 0) {
buf[n] = count_all[m];
n += nvalues;
}
delete [] count_one;
delete [] count_all;
}

407
src/compute_temp_chunk.cpp Normal file
View File

@ -0,0 +1,407 @@
/* ----------------------------------------------------------------------
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 "string.h"
#include "compute_temp_chunk.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "modify.h"
#include "compute_chunk_atom.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeTempChunk::ComputeTempChunk(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg < 4) error->all(FLERR,"Illegal compute temp/chunk command");
vector_flag = 1;
size_vector = 1;
size_vector_variable = 1;
extvector = 0;
// ID of compute chunk/atom
int n = strlen(arg[3]) + 1;
idchunk = new char[n];
strcpy(idchunk,arg[3]);
biasflag = 0;
init();
// optional args
comflag = 0;
biasflag = 0;
id_bias = NULL;
adof = domain->dimension;
cdof = 0.0;
int iarg = 4;
while (iarg < narg) {
if (strcmp(arg[iarg],"com") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute temp/chunk command");
if (strcmp(arg[iarg+1],"yes") == 0) comflag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) comflag = 0;
else error->all(FLERR,"Illegal compute temp/chunk command");
iarg += 2;
} else if (strcmp(arg[iarg],"bias") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute temp/chunk command");
biasflag = 1;
int n = strlen(arg[iarg+1]) + 1;
id_bias = new char[n];
strcpy(id_bias,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"adof") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute temp/chunk command");
adof = force->numeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"cdof") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute temp/chunk command");
cdof = force->numeric(FLERR,arg[iarg+1]);
iarg += 2;
} else error->all(FLERR,"Illegal compute temp/chunk command");
}
if (comflag && biasflag)
error->all(FLERR,"Cannot use both com and bias with compute temp/chunk");
// error check on bias compute
if (biasflag) {
int i = modify->find_compute(id_bias);
if (i < 0)
error->all(FLERR,"Could not find compute ID for temperature bias");
tbias = modify->compute[i];
if (tbias->tempflag == 0)
error->all(FLERR,"Bias compute does not calculate temperature");
if (tbias->tempbias == 0)
error->all(FLERR,"Bias compute does not calculate a velocity bias");
}
// chunk-based data
nchunk = 1;
maxchunk = 0;
ke = keall = NULL;
count = countall = NULL;
massproc = masstotal = NULL;
vcm = vcmall = NULL;
allocate();
}
/* ---------------------------------------------------------------------- */
ComputeTempChunk::~ComputeTempChunk()
{
delete [] idchunk;
delete [] id_bias;
memory->destroy(ke);
memory->destroy(keall);
memory->destroy(count);
memory->destroy(countall);
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(vcm);
memory->destroy(vcmall);
}
/* ---------------------------------------------------------------------- */
void ComputeTempChunk::init()
{
int icompute = modify->find_compute(idchunk);
if (icompute < 0)
error->all(FLERR,"Chunk/atom compute does not exist for "
"compute temp/chunk");
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
if (strcmp(cchunk->style,"chunk/atom") != 0)
error->all(FLERR,"Compute temp/chunk does not use chunk/atom compute");
if (biasflag) {
int i = modify->find_compute(id_bias);
if (i < 0)
error->all(FLERR,"Could not find compute ID for temperature bias");
tbias = modify->compute[i];
}
}
/* ---------------------------------------------------------------------- */
void ComputeTempChunk::compute_vector()
{
int index;
invoked_vector = update->ntimestep;
// compute chunk/atom assigns atoms to chunk IDs
// extract ichunk index vector from compute
// ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms
nchunk = cchunk->setup_chunks();
cchunk->compute_ichunk();
int *ichunk = cchunk->ichunk;
if (nchunk > maxchunk) allocate();
// calculate COM velocity for each chunk
if (comflag) vcm_compute();
// remove velocity bias
if (biasflag) {
if (tbias->invoked_scalar != update->ntimestep) tbias->compute_scalar();
tbias->remove_bias_all();
}
// zero local per-chunk values
for (int i = 0; i < nchunk; i++) {
count[i] = 0;
ke[i] = 0.0;
}
// compute temperature for each chunk
// option for removing COM velocity
double **v = atom->v;
double *mass = atom->mass;
double *rmass = atom->rmass;
int *mask = atom->mask;
int *type = atom->type;
int nlocal = atom->nlocal;
if (!comflag) {
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
index = ichunk[i]-1;
if (index < 0) continue;
ke[index] += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
rmass[i];
count[index]++;
}
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
index = ichunk[i]-1;
if (index < 0) continue;
ke[index] += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
mass[type[i]];
count[index]++;
}
}
} else {
double vx,vy,vz;
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
index = ichunk[i]-1;
if (index < 0) continue;
vx = v[i][0] - vcmall[index][0];
vy = v[i][1] - vcmall[index][1];
vz = v[i][2] - vcmall[index][2];
ke[index] += (vx*vx + vy*vy + vz*vz) * rmass[i];
count[index]++;
}
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
index = ichunk[i]-1;
if (index < 0) continue;
vx = v[i][0] - vcmall[index][0];
vy = v[i][1] - vcmall[index][1];
vz = v[i][2] - vcmall[index][2];
ke[index] += (vx*vx + vy*vy + vz*vz) * mass[type[i]];
count[index]++;
}
}
}
MPI_Allreduce(ke,keall,nchunk,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(count,countall,nchunk,MPI_INT,MPI_SUM,world);
// restore velocity bias
if (biasflag) tbias->restore_bias_all();
// normalize temperatures by per-chunk DOF
double dof,tfactor;
double mvv2e = force->mvv2e;
double boltz = force->boltz;
for (int i = 0; i < nchunk; i++) {
dof = cdof + adof*countall[i];
if (dof > 0.0) tfactor = mvv2e / (dof * boltz);
else tfactor = 0.0;
keall[i] *= tfactor;
}
}
/* ----------------------------------------------------------------------
calculate velocity of COM for each chunk
------------------------------------------------------------------------- */
void ComputeTempChunk::vcm_compute()
{
int index;
double massone;
int *ichunk = cchunk->ichunk;
for (int i = 0; i < nchunk; i++) {
vcm[i][0] = vcm[i][1] = vcm[i][2] = 0.0;
massproc[i] = 0.0;
}
double **v = atom->v;
int *mask = atom->mask;
int *type = atom->type;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
index = ichunk[i]-1;
if (index < 0) continue;
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
vcm[index][0] += v[i][0] * massone;
vcm[index][1] += v[i][1] * massone;
vcm[index][2] += v[i][2] * massone;
massproc[index] += massone;
}
MPI_Allreduce(&vcm[0][0],&vcmall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world);
for (int i = 0; i < nchunk; i++) {
vcmall[i][0] /= masstotal[i];
vcmall[i][1] /= masstotal[i];
vcmall[i][2] /= masstotal[i];
}
}
/* ----------------------------------------------------------------------
lock methods: called by fix ave/time
these methods insure vector/array size is locked for Nfreq epoch
by passing lock info along to compute chunk/atom
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
increment lock counter
------------------------------------------------------------------------- */
void ComputeTempChunk::lock_enable()
{
cchunk->lockcount++;
}
/* ----------------------------------------------------------------------
decrement lock counter in compute chunk/atom, it if still exists
------------------------------------------------------------------------- */
void ComputeTempChunk::lock_disable()
{
int icompute = modify->find_compute(idchunk);
if (icompute >= 0) {
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
cchunk->lockcount--;
}
}
/* ----------------------------------------------------------------------
calculate and return # of chunks = length of vector/array
------------------------------------------------------------------------- */
int ComputeTempChunk::lock_length()
{
nchunk = cchunk->setup_chunks();
return nchunk;
}
/* ----------------------------------------------------------------------
set the lock from startstep to stopstep
------------------------------------------------------------------------- */
void ComputeTempChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep)
{
cchunk->lock(fixptr,startstep,stopstep);
}
/* ----------------------------------------------------------------------
unset the lock
------------------------------------------------------------------------- */
void ComputeTempChunk::unlock(Fix *fixptr)
{
cchunk->unlock(fixptr);
}
/* ----------------------------------------------------------------------
free and reallocate per-chunk arrays
------------------------------------------------------------------------- */
void ComputeTempChunk::allocate()
{
memory->destroy(ke);
memory->destroy(keall);
memory->destroy(count);
memory->destroy(countall);
size_vector = maxchunk = nchunk;
memory->create(ke,maxchunk,"temp/chunk:ke");
memory->create(keall,maxchunk,"temp/chunk:keall");
memory->create(count,maxchunk,"temp/chunk:count");
memory->create(countall,maxchunk,"temp/chunk:countall");
vector = keall;
if (comflag) {
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(vcm);
memory->destroy(vcmall);
memory->create(massproc,maxchunk,"vcm/chunk:massproc");
memory->create(masstotal,maxchunk,"vcm/chunk:masstotal");
memory->create(vcm,maxchunk,3,"vcm/chunk:vcm");
memory->create(vcmall,maxchunk,3,"vcm/chunk:vcmall");
}
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputeTempChunk::memory_usage()
{
double bytes = (bigint) maxchunk * 2 * sizeof(double);
bytes = (bigint) maxchunk * 2 * sizeof(int);
if (comflag) {
bytes += (bigint) maxchunk * 2 * sizeof(double);
bytes += (bigint) maxchunk * 2*3 * sizeof(double);
}
return bytes;
}

80
src/compute_temp_chunk.h Normal file
View File

@ -0,0 +1,80 @@
/* -*- c++ -*- ----------------------------------------------------------
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 COMPUTE_CLASS
ComputeStyle(temp/chunk,ComputeTempChunk)
#else
#ifndef LMP_COMPUTE_TEMP_CHUNK_H
#define LMP_COMPUTE_TEMP_CHUNK_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeTempChunk : public Compute {
public:
ComputeTempChunk(class LAMMPS *, int, char **);
~ComputeTempChunk();
void init();
void compute_vector();
void lock_enable();
void lock_disable();
int lock_length();
void lock(class Fix *, bigint, bigint);
void unlock(class Fix *);
double memory_usage();
private:
int nchunk,maxchunk,comflag,biasflag;
char *idchunk;
class ComputeChunkAtom *cchunk;
double adof,cdof;
char *id_bias;
class Compute *tbias; // ptr to additional bias compute
double *ke,*keall;
int *count,*countall;
double *massproc,*masstotal;
double **vcm,**vcmall;
void vcm_compute();
void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Compute com/molecule requires molecular atom style
Self-explanatory.
E: Molecule count changed in compute com/molecule
Number of molecules must remain constant over time.
*/

252
src/compute_torque_chunk.cpp Normal file
View File

@ -0,0 +1,252 @@
/* ----------------------------------------------------------------------
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 "string.h"
#include "compute_torque_chunk.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "compute_chunk_atom.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeTorqueChunk::ComputeTorqueChunk(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 4) error->all(FLERR,"Illegal compute inertia/chunk command");
array_flag = 1;
size_array_cols = 3;
size_array_rows = 0;
size_array_rows_variable = 1;
extarray = 0;
// ID of compute chunk/atom
int n = strlen(arg[6]) + 1;
idchunk = new char[n];
strcpy(idchunk,arg[6]);
init();
// chunk-based data
nchunk = 1;
maxchunk = 0;
massproc = masstotal = NULL;
com = comall = NULL;
torque = torqueall = NULL;
allocate();
}
/* ---------------------------------------------------------------------- */
ComputeTorqueChunk::~ComputeTorqueChunk()
{
delete [] idchunk;
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(com);
memory->destroy(comall);
memory->destroy(torque);
memory->destroy(torqueall);
}
/* ---------------------------------------------------------------------- */
void ComputeTorqueChunk::init()
{
int icompute = modify->find_compute(idchunk);
if (icompute < 0)
error->all(FLERR,"Chunk/atom compute does not exist for "
"compute torque/chunk");
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
if (strcmp(cchunk->style,"chunk/atom") != 0)
error->all(FLERR,"Compute torque/chunk does not use chunk/atom compute");
}
/* ---------------------------------------------------------------------- */
void ComputeTorqueChunk::compute_array()
{
int i,j,index;
double dx,dy,dz,massone;
double unwrap[3];
invoked_array = update->ntimestep;
// compute chunk/atom assigns atoms to chunk IDs
// extract ichunk index vector from compute
// ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms
nchunk = cchunk->setup_chunks();
cchunk->compute_ichunk();
int *ichunk = cchunk->ichunk;
if (nchunk > maxchunk) allocate();
// zero local per-chunk values
for (int i = 0; i < nchunk; i++) {
massproc[i] = 0.0;
com[i][0] = com[i][1] = com[i][2] = 0.0;
torque[i][0] = torque[i][1] = torque[i][2] = 0.0;
}
// compute COM for each chunk
double **x = atom->x;
int *mask = atom->mask;
int *type = atom->type;
imageint *image = atom->image;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
index = ichunk[i]-1;
if (index < 0) continue;
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
domain->unmap(x[i],image[i],unwrap);
massproc[index] += massone;
com[index][0] += unwrap[0] * massone;
com[index][1] += unwrap[1] * massone;
com[index][2] += unwrap[2] * massone;
}
MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world);
for (int i = 0; i < nchunk; i++) {
comall[i][0] /= masstotal[i];
comall[i][1] /= masstotal[i];
comall[i][2] /= masstotal[i];
}
// compute torque on each chunk
double **f = atom->f;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
index = ichunk[i]-1;
if (index < 0) continue;
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - comall[index][0];
dy = unwrap[1] - comall[index][1];
dz = unwrap[2] - comall[index][2];
torque[i][0] += dy*f[i][2] - dz*f[i][1];
torque[i][1] += dz*f[i][0] - dx*f[i][2];
torque[i][2] += dx*f[i][1] - dy*f[i][0];
}
MPI_Allreduce(&torque[0][0],&torqueall[0][0],3*nchunk,
MPI_DOUBLE,MPI_SUM,world);
}
/* ----------------------------------------------------------------------
lock methods: called by fix ave/time
these methods insure vector/array size is locked for Nfreq epoch
by passing lock info along to compute chunk/atom
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
increment lock counter
------------------------------------------------------------------------- */
void ComputeTorqueChunk::lock_enable()
{
cchunk->lockcount++;
}
/* ----------------------------------------------------------------------
decrement lock counter in compute chunk/atom, it if still exists
------------------------------------------------------------------------- */
void ComputeTorqueChunk::lock_disable()
{
int icompute = modify->find_compute(idchunk);
if (icompute >= 0) {
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
cchunk->lockcount--;
}
}
/* ----------------------------------------------------------------------
calculate and return # of chunks = length of vector/array
------------------------------------------------------------------------- */
int ComputeTorqueChunk::lock_length()
{
nchunk = cchunk->setup_chunks();
return nchunk;
}
/* ----------------------------------------------------------------------
set the lock from startstep to stopstep
------------------------------------------------------------------------- */
void ComputeTorqueChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep)
{
cchunk->lock(fixptr,startstep,stopstep);
}
/* ----------------------------------------------------------------------
unset the lock
------------------------------------------------------------------------- */
void ComputeTorqueChunk::unlock(Fix *fixptr)
{
cchunk->unlock(fixptr);
}
/* ----------------------------------------------------------------------
free and reallocate per-chunk arrays
------------------------------------------------------------------------- */
void ComputeTorqueChunk::allocate()
{
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(com);
memory->destroy(comall);
memory->destroy(torque);
memory->destroy(torqueall);
size_array_rows = maxchunk = nchunk;
memory->create(massproc,maxchunk,"torque/chunk:massproc");
memory->create(masstotal,maxchunk,"torque/chunk:masstotal");
memory->create(com,maxchunk,3,"torque/chunk:com");
memory->create(comall,maxchunk,3,"torque/chunk:comall");
memory->create(torque,maxchunk,6,"torque/chunk:torque");
memory->create(torqueall,maxchunk,6,"torque/chunk:torqueall");
array = torqueall;
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputeTorqueChunk::memory_usage()
{
double bytes = (bigint) maxchunk * 2 * sizeof(double);
bytes += (bigint) maxchunk * 2*3 * sizeof(double);
bytes += (bigint) maxchunk * 2*3 * sizeof(double);
return bytes;
}

75
src/compute_torque_chunk.h Normal file
View File

@ -0,0 +1,75 @@
/* -*- c++ -*- ----------------------------------------------------------
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 COMPUTE_CLASS
ComputeStyle(torque/chunk,ComputeTorqueChunk)
#else
#ifndef LMP_COMPUTE_TORQUE_CHUNK_H
#define LMP_COMPUTE_TORQUE_CHUNK_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeTorqueChunk : public Compute {
public:
ComputeTorqueChunk(class LAMMPS *, int, char **);
~ComputeTorqueChunk();
void init();
void compute_array();
void lock_enable();
void lock_disable();
int lock_length();
void lock(class Fix *, bigint, bigint);
void unlock(class Fix *);
double memory_usage();
private:
int nchunk,maxchunk;
char *idchunk;
class ComputeChunkAtom *cchunk;
double *massproc,*masstotal;
double **com,**comall;
double **torque,**torqueall;
void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Compute inertia/molecule requires molecular atom style
Self-explanatory.
E: Molecule count changed in compute inertia/molecule
Number of molecules must remain constant over time.
*/

237
src/compute_vcm_chunk.cpp Normal file
View File

@ -0,0 +1,237 @@
/* ----------------------------------------------------------------------
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 "string.h"
#include "compute_vcm_chunk.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "compute_chunk_atom.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
enum{ONCE,NFREQ,EVERY};
/* ---------------------------------------------------------------------- */
ComputeVCMChunk::ComputeVCMChunk(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 4) error->all(FLERR,"Illegal compute vcm/chunk command");
array_flag = 1;
size_array_cols = 3;
size_array_rows = 0;
size_array_rows_variable = 1;
extarray = 0;
// ID of compute chunk/atom
int n = strlen(arg[6]) + 1;
idchunk = new char[n];
strcpy(idchunk,arg[6]);
init();
// chunk-based data
nchunk = 1;
maxchunk = 0;
massproc = masstotal = NULL;
vcm = vcmall = NULL;
allocate();
firstflag = massneed = 1;
}
/* ---------------------------------------------------------------------- */
ComputeVCMChunk::~ComputeVCMChunk()
{
delete [] idchunk;
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(vcm);
memory->destroy(vcmall);
}
/* ---------------------------------------------------------------------- */
void ComputeVCMChunk::init()
{
int icompute = modify->find_compute(idchunk);
if (icompute < 0)
error->all(FLERR,"Chunk/atom compute does not exist for compute vcm/chunk");
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
if (strcmp(cchunk->style,"chunk/atom") != 0)
error->all(FLERR,"Compute vcm/chunk does not use chunk/atom compute");
}
/* ---------------------------------------------------------------------- */
void ComputeVCMChunk::setup()
{
// one-time calculation of per-chunk mass
// done in setup, so that ComputeChunkAtom::setup() is already called
if (firstflag && cchunk->idsflag == ONCE) {
firstflag = massneed = 0;
compute_array();
}
}
/* ---------------------------------------------------------------------- */
void ComputeVCMChunk::compute_array()
{
int index;
double massone;
invoked_array = update->ntimestep;
// compute chunk/atom assigns atoms to chunk IDs
// extract ichunk index vector from compute
// ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms
nchunk = cchunk->setup_chunks();
cchunk->compute_ichunk();
int *ichunk = cchunk->ichunk;
if (nchunk > maxchunk) allocate();
// zero local per-chunk values
for (int i = 0; i < nchunk; i++)
vcm[i][0] = vcm[i][1] = vcm[i][2] = 0.0;
if (massneed)
for (int i = 0; i < nchunk; i++) massproc[i] = 0.0;
// compute VCM for each chunk
double **v = atom->v;
int *mask = atom->mask;
int *type = atom->type;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
index = ichunk[i]-1;
if (index < 0) continue;
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
vcm[index][0] += v[i][0] * massone;
vcm[index][1] += v[i][1] * massone;
vcm[index][2] += v[i][2] * massone;
if (massneed) massproc[index] += massone;
}
MPI_Allreduce(&vcm[0][0],&vcmall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world);
if (massneed)
MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world);
for (int i = 0; i < nchunk; i++) {
vcmall[i][0] /= masstotal[i];
vcmall[i][1] /= masstotal[i];
vcmall[i][2] /= masstotal[i];
}
}
/* ----------------------------------------------------------------------
lock methods: called by fix ave/time
these methods insure vector/array size is locked for Nfreq epoch
by passing lock info along to compute chunk/atom
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
increment lock counter
------------------------------------------------------------------------- */
void ComputeVCMChunk::lock_enable()
{
cchunk->lockcount++;
}
/* ----------------------------------------------------------------------
decrement lock counter in compute chunk/atom, it if still exists
------------------------------------------------------------------------- */
void ComputeVCMChunk::lock_disable()
{
int icompute = modify->find_compute(idchunk);
if (icompute >= 0) {
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
cchunk->lockcount--;
}
}
/* ----------------------------------------------------------------------
calculate and return # of chunks = length of vector/array
------------------------------------------------------------------------- */
int ComputeVCMChunk::lock_length()
{
nchunk = cchunk->setup_chunks();
return nchunk;
}
/* ----------------------------------------------------------------------
set the lock from startstep to stopstep
------------------------------------------------------------------------- */
void ComputeVCMChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep)
{
cchunk->lock(fixptr,startstep,stopstep);
}
/* ----------------------------------------------------------------------
unset the lock
------------------------------------------------------------------------- */
void ComputeVCMChunk::unlock(Fix *fixptr)
{
cchunk->unlock(fixptr);
}
/* ----------------------------------------------------------------------
free and reallocate per-chunk arrays
------------------------------------------------------------------------- */
void ComputeVCMChunk::allocate()
{
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(vcm);
memory->destroy(vcmall);
size_array_rows = maxchunk = nchunk;
memory->create(massproc,maxchunk,"vcm/chunk:massproc");
memory->create(masstotal,maxchunk,"vcm/chunk:masstotal");
memory->create(vcm,maxchunk,3,"vcm/chunk:vcm");
memory->create(vcmall,maxchunk,3,"vcm/chunk:vcmall");
array = vcmall;
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputeVCMChunk::memory_usage()
{
double bytes = (bigint) maxchunk * 2 * sizeof(double);
bytes += (bigint) maxchunk * 2*3 * sizeof(double);
return bytes;
}

28
src/compute_vcm_molecule.h → src/compute_vcm_chunk.h
View File

@ -13,31 +13,43 @@
#ifdef COMPUTE_CLASS
ComputeStyle(vcm/molecule,ComputeVCMMolecule)
ComputeStyle(vcm/chunk,ComputeVCMChunk)
#else
#ifndef LMP_COMPUTE_VCM_MOLECULE_H
#define LMP_COMPUTE_VCM_MOLECULE_H
#ifndef LMP_COMPUTE_VCM_CHUNK_H
#define LMP_COMPUTE_VCM_CHUNK_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeVCMMolecule : public Compute {
class ComputeVCMChunk : public Compute {
public:
ComputeVCMMolecule(class LAMMPS *, int, char **);
~ComputeVCMMolecule();
ComputeVCMChunk(class LAMMPS *, int, char **);
~ComputeVCMChunk();
void init();
void setup();
void compute_array();
void lock_enable();
void lock_disable();
int lock_length();
void lock(class Fix *, bigint, bigint);
void unlock(class Fix *);
double memory_usage();
private:
int nmolecules;
tagint idlo,idhi;
int nchunk,maxchunk;
int firstflag,massneed;
char *idchunk;
class ComputeChunkAtom *cchunk;
double *massproc,*masstotal;
double **vcm,**vcmall;
void allocate();
};
}

145
src/compute_vcm_molecule.cpp
View File

@ -1,145 +0,0 @@
/* ----------------------------------------------------------------------
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 "compute_vcm_molecule.h"
#include "atom.h"
#include "update.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeVCMMolecule::ComputeVCMMolecule(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 3) error->all(FLERR,"Illegal compute vcm/molecule command");
if (atom->molecular == 0)
error->all(FLERR,"Compute vcm/molecule requires molecular atom style");
array_flag = 1;
size_array_cols = 3;
extarray = 0;
// setup molecule-based data
nmolecules = molecules_in_group(idlo,idhi);
size_array_rows = nmolecules;
memory->create(massproc,nmolecules,"vcm/molecule:massproc");
memory->create(masstotal,nmolecules,"vcm/molecule:masstotal");
memory->create(vcm,nmolecules,3,"vcm/molecule:vcm");
memory->create(vcmall,nmolecules,3,"vcm/molecule:vcmall");
array = vcmall;
// compute masstotal for each molecule
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *type = atom->type;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
tagint imol;
double massone;
for (int i = 0; i < nmolecules; i++) massproc[i] = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
massproc[imol] += massone;
}
MPI_Allreduce(massproc,masstotal,nmolecules,MPI_DOUBLE,MPI_SUM,world);
}
/* ---------------------------------------------------------------------- */
ComputeVCMMolecule::~ComputeVCMMolecule()
{
memory->destroy(massproc);
memory->destroy(masstotal);
memory->destroy(vcm);
memory->destroy(vcmall);
}
/* ---------------------------------------------------------------------- */
void ComputeVCMMolecule::init()
{
int ntmp = molecules_in_group(idlo,idhi);
if (ntmp != nmolecules)
error->all(FLERR,"Molecule count changed in compute vcm/molecule");
}
/* ---------------------------------------------------------------------- */
void ComputeVCMMolecule::compute_array()
{
tagint imol;
double massone;
invoked_array = update->ntimestep;
for (int i = 0; i < nmolecules; i++)
vcm[i][0] = vcm[i][1] = vcm[i][2] = 0.0;
double **v = atom->v;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *type = atom->type;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
vcm[imol][0] += v[i][0] * massone;
vcm[imol][1] += v[i][1] * massone;
vcm[imol][2] += v[i][2] * massone;
}
MPI_Allreduce(&vcm[0][0],&vcmall[0][0],3*nmolecules,
MPI_DOUBLE,MPI_SUM,world);
for (int i = 0; i < nmolecules; i++) {
vcmall[i][0] /= masstotal[i];
vcmall[i][1] /= masstotal[i];
vcmall[i][2] /= masstotal[i];
}
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputeVCMMolecule::memory_usage()
{
double bytes = (bigint) nmolecules * 2 * sizeof(double);
if (molmap) bytes += (idhi-idlo+1) * sizeof(int);
bytes += (bigint) nmolecules * 2*3 * sizeof(double);
return bytes;
}

1038
src/fix_ave_chunk.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

173
src/fix_ave_chunk.h Normal file
View File

@ -0,0 +1,173 @@
/* -*- c++ -*- ----------------------------------------------------------
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(ave/chunk,FixAveChunk)
#else
#ifndef LMP_FIX_AVE_CHUNK_H
#define LMP_FIX_AVE_CHUNK_H
#include "stdio.h"
#include "fix.h"
namespace LAMMPS_NS {
class FixAveChunk : public Fix {
public:
FixAveChunk(class LAMMPS *, int, char **);
~FixAveChunk();
int setmask();
void init();
void setup(int);
void end_of_step();
double compute_array(int,int);
double memory_usage();
void reset_timestep(bigint);
private:
int me,nvalues;
int nrepeat,nfreq,irepeat;
int normflag,scaleflag,overwrite,biasflag,colextra;
bigint nvalid;
double adof,cdof;
char *tstring,*sstring,*id_bias;
int *which,*argindex,*value2index;
char **ids;
class Compute *tbias; // ptr to additional bias compute
FILE *fp;
int ave,nwindow;
int normcount,iwindow,window_limit;
int nchunk,maxchunk;
char *idchunk;
class ComputeChunkAtom *cchunk;
long filepos;
int maxvar;
double *varatom;
// one,many,sum vecs/arrays are used with a single Nfreq epoch
// total,list vecs/arrays are used across epochs
double *count_one,*count_many,*count_sum;
double **values_one,**values_many,**values_sum;
double *count_total,**count_list;
double **values_total,***values_list;
void allocate();
bigint nextvalid();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Cannot use fix ave/spatial z for 2 dimensional model
Self-explanatory.
E: Same dimension twice in fix ave/spatial
Self-explanatory.
E: Region ID for fix ave/spatial does not exist
Self-explanatory.
E: Cannot open fix ave/spatial file %s
The specified file cannot be opened. Check that the path and name are
correct.
E: Compute ID for fix ave/spatial does not exist
Self-explanatory.
E: Fix ave/spatial compute does not calculate per-atom values
A compute used by fix ave/spatial must generate per-atom values.
E: Fix ave/spatial compute does not calculate a per-atom vector
A compute used by fix ave/spatial must generate per-atom values.
E: Fix ave/spatial compute does not calculate a per-atom array
Self-explanatory.
E: Fix ave/spatial compute vector is accessed out-of-range
The index for the vector is out of bounds.
E: Fix ID for fix ave/spatial does not exist
Self-explanatory.
E: Fix ave/spatial fix does not calculate per-atom values
A fix used by fix ave/spatial must generate per-atom values.
E: Fix ave/spatial fix does not calculate a per-atom vector
A fix used by fix ave/spatial must generate per-atom values.
E: Fix ave/spatial fix does not calculate a per-atom array
Self-explanatory.
E: Fix ave/spatial fix vector is accessed out-of-range
The index for the vector is out of bounds.
E: Variable name for fix ave/spatial does not exist
Self-explanatory.
E: Fix ave/spatial variable is not atom-style variable
A variable used by fix ave/spatial must generate per-atom values.
E: Fix ave/spatial for triclinic boxes requires units reduced
Self-explanatory.
E: Fix ave/spatial settings invalid with changing box size
If the box size changes, only the units reduced option can be
used.
E: Fix for fix ave/spatial not computed at compatible time
Fixes generate their values on specific timesteps. Fix ave/spatial is
requesting a value on a non-allowed timestep.
E: Fix ave/spatial missed timestep
You cannot reset the timestep to a value beyond where the fix
expects to next perform averaging.
*/

7
src/fix_ave_spatial.cpp
View File

@ -29,6 +29,7 @@
#include "compute.h"
#include "input.h"
#include "variable.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
@ -49,6 +50,12 @@ enum{NODISCARD,MIXED,YESDISCARD};
FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (comm->me == 0)
error->warning(FLERR,"The fix ave/spatial command has been replaced "
"by the more flexible fix ave/chunk and compute chunk/atom "
"commands -- fix ave/spatial will be removed in "
"the summer of 2015");
if (narg < 6) error->all(FLERR,"Illegal fix ave/spatial command");
MPI_Comm_rank(world,&me);