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

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;
+}

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

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;
+}

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);
 };
 
 }