From 93145d50a4e7735780c49a3b0b9e64f13023ccfd Mon Sep 17 00:00:00 2001 From: sjplimp Date: Mon, 18 Oct 2010 18:33:01 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5064 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/compute_atom_molecule.cpp | 336 ++++++++++++++++++++++++++++++ src/compute_atom_molecule.h | 54 +++++ src/compute_property_molecule.cpp | 45 +++- src/compute_property_molecule.h | 1 + 4 files changed, 430 insertions(+), 6 deletions(-) create mode 100644 src/compute_atom_molecule.cpp create mode 100644 src/compute_atom_molecule.h diff --git a/src/compute_atom_molecule.cpp b/src/compute_atom_molecule.cpp new file mode 100644 index 0000000000..08dc7b5944 --- /dev/null +++ b/src/compute_atom_molecule.cpp @@ -0,0 +1,336 @@ +/* ---------------------------------------------------------------------- + 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("Illegal compute atom/molecule command"); + + if (atom->molecular == 0) + error->all("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("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("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("Compute ID for compute atom/molecule does not exist"); + if (modify->compute[icompute]->peratom_flag == 0) + error->all("Compute atom/molecule compute does not " + "calculate per-atom values"); + if (argindex[i] == 0 && + modify->compute[icompute]->size_peratom_cols != 0) + error->all("Compute atom/molecule compute does not " + "calculate a per-atom vector"); + if (argindex[i] && modify->compute[icompute]->size_peratom_cols == 0) + error->all("Compute atom/molecule compute does not " + "calculate a per-atom array"); + if (argindex[i] && + argindex[i] > modify->compute[icompute]->size_peratom_cols) + error->all("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("Fix ID for compute atom/molecule does not exist"); + if (modify->fix[ifix]->peratom_flag) + error->all("Compute atom/molecule fix does not " + "calculate per-atom values"); + if (argindex[i] == 0 && + modify->fix[ifix]->size_peratom_cols != 0) + error->all("Compute atom/molecule fix does not " + "calculate a per-atom vector"); + if (argindex[i] && modify->fix[ifix]->size_peratom_cols == 0) + error->all("Compute atom/molecule fix does not " + "calculate a per-atom array"); + if (argindex[i] && + argindex[i] > modify->fix[ifix]->size_peratom_cols) + error->all("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("Variable name for compute atom/molecule does not exist"); + if (input->variable->atomstyle(ivariable) == 0) + error->all("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) { + vone = (double *) memory->smalloc(nmolecules*sizeof(double), + "atom/molecule:vone"); + vector = (double *) memory->smalloc(nmolecules*sizeof(double), + "atom/molecule:vector"); + vector_flag = 1; + size_vector = nmolecules; + extvector = 0; + } else { + aone = memory->create_2d_double_array(nmolecules,nvalues, + "atom/molecule:aone"); + array = memory->create_2d_double_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() +{ + memory->sfree(vone); + memory->sfree(vector); + memory->sfree(aone); + memory->destroy_2d_double_array(array); + memory->sfree(scratch); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeAtomMolecule::init() +{ + int ntmp = molecules_in_group(idlo,idhi); + if (ntmp != nmolecules) + error->all("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("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("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("Variable name for compute atom/molecule does not exist"); + value2index[m] = ivariable; + + } else value2index[m] = -1; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeAtomMolecule::compute_vector() +{ + int i,j,m,imol; + + invoked_vector = update->ntimestep; + + for (m = 0; m < nmolecules; m++) vone[m] = 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] = scratch[j]; + } + j += nstride; + } + + 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 (n = 0; n < nvalues; n++) { + for (m = 0; m < nmolecules; m++) aone[n][m] = 0.0; + compute_one(n); + + j = 0; + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + imol = molecule[i]; + if (molmap) imol = molmap[imol-idlo]; + else imol--; + aone[n][imol] = scratch[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 vector if necessary +------------------------------------------------------------------------- */ + +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) { + scratch = compute->vector_atom; + nstride = 1; + } else { + scratch = &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("Fix used in compute atom/molecule not computed " + "at compatible time"); + Fix *fix = modify->fix[vidx]; + + if (aidx == 0) { + scratch = fix->vector_atom; + nstride = 1; + } else { + scratch = &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->sfree(scratch); + scratch = (double *) + memory->smalloc(maxatom*sizeof(double),"atom/molecule:scratch"); + } + + input->variable->compute_atom(vidx,igroup,scratch,1,0); + } +} + +/* ---------------------------------------------------------------------- + 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; +} diff --git a/src/compute_atom_molecule.h b/src/compute_atom_molecule.h new file mode 100644 index 0000000000..322c1a92ec --- /dev/null +++ b/src/compute_atom_molecule.h @@ -0,0 +1,54 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef 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; + int idlo,idhi; + + int *which,*argindex,*flavor,*value2index; + char **ids; + + int nstride,maxatom; + double *vone; + double **aone; + double *scratch; + + void compute_one(int); +}; + +} + +#endif +#endif diff --git a/src/compute_property_molecule.cpp b/src/compute_property_molecule.cpp index ca50dcaf66..16c30ab855 100644 --- a/src/compute_property_molecule.cpp +++ b/src/compute_property_molecule.cpp @@ -32,7 +32,6 @@ ComputePropertyMolecule(LAMMPS *lmp, int narg, char **arg) : error->all("Compute property/molecule requires molecular atom style"); nvalues = narg - 3; - if (nvalues == 1) size_array_cols = 0; pack_choice = new FnPtrPack[nvalues]; @@ -42,6 +41,8 @@ ComputePropertyMolecule(LAMMPS *lmp, int narg, char **arg) : 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("Invalid keyword in compute property/molecule command"); } @@ -128,13 +129,45 @@ double ComputePropertyMolecule::memory_usage() customize a new keyword by adding a method ------------------------------------------------------------------------- */ -/* ---------------------------------------------------------------------- */ - void ComputePropertyMolecule::pack_mol(int n) { - for (int i = idlo; i <= idhi; i++) - if (molmap == NULL || molmap[i-idlo] >= 0) { - buf[n] = i; + for (int 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,imol; + + int *count_one = new int[nmolecules]; + for (m = 0; m < nmolecules; m++) count_one[m] = 0; + + int *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; +} diff --git a/src/compute_property_molecule.h b/src/compute_property_molecule.h index ce0f0fcae3..4317c67f6b 100644 --- a/src/compute_property_molecule.h +++ b/src/compute_property_molecule.h @@ -43,6 +43,7 @@ class ComputePropertyMolecule : public Compute { FnPtrPack *pack_choice; // ptrs to pack functions void pack_mol(int); + void pack_count(int); }; }