Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F70062458
compute_atom_molecule.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Jul 5, 02:31
Size
9 KB
Mime Type
text/x-c
Expires
Sun, Jul 7, 02:31 (2 d)
Engine
blob
Format
Raw Data
Handle
18495116
Attached To
rLAMMPS lammps
compute_atom_molecule.cpp
View Options
/* ----------------------------------------------------------------------
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)
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,imol;
invoked_vector = update->ntimestep;
for (n = 0; n < nmolecules; n++) vone[n] = 0.0;
compute_one(0);
int *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,imol;
invoked_array = update->ntimestep;
int *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;
}
}
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_array_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_array_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 = 2*nmolecules*nvalues * sizeof(double);
if (molmap) bytes += (idhi-idlo+1) * sizeof(int);
bytes += maxatom * sizeof(double);
return bytes;
}
Event Timeline
Log In to Comment