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

src/compute_group_group.cpp

@@ -65,7 +65,9 @@ ComputeGroupGroup::~ComputeGroupGroup()
 
 void ComputeGroupGroup::init()
 {
-  if (force->pair == NULL || force->pair->single_enable == 0)
+  if (force->pair == NULL)
+    error->all("Pair style does not support compute group/group");
+  if (force->pair->single_enable == 0)
     error->all("Pair style does not support compute group/group");
 
   pair = force->pair;

src/compute_pair_local.cpp

@@ -0,0 +1,236 @@
+/* ----------------------------------------------------------------------
+   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_pair_local.h"
+#include "atom.h"
+#include "update.h"
+#include "force.h"
+#include "pair.h"
+#include "neighbor.h"
+#include "neigh_request.h"
+#include "neigh_list.h"
+#include "group.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define DELTA 10000
+
+/* ---------------------------------------------------------------------- */
+
+ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg)
+{
+  if (narg < 4) error->all("Illegal compute pair/local command");
+
+  local_flag = 1;
+  nvalues = narg - 3;
+  if (nvalues == 1) size_local_cols = 0;
+  else size_local_cols = nvalues;
+
+  dflag = eflag = fflag = -1;
+  nvalues = 0;
+
+  int i;
+  for (int iarg = 3; iarg < narg; iarg++) {
+    i = iarg-3;
+    if (strcmp(arg[iarg],"dist") == 0) dflag = nvalues++;
+    else if (strcmp(arg[iarg],"eng") == 0) eflag = nvalues++;
+    else if (strcmp(arg[iarg],"force") == 0) fflag = nvalues++;
+    else error->all("Invalid keyword in compute pair/local command");
+  }
+
+  nmax = 0;
+  vector = NULL;
+  array = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputePairLocal::~ComputePairLocal()
+{
+  memory->sfree(vector);
+  memory->destroy_2d_double_array(array);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePairLocal::init()
+{
+  if (force->pair == NULL) 
+    error->all("No pair style is defined for compute pair/local");
+  if (force->pair->single_enable == 0)
+    error->all("Pair style does not support compute pair/local");
+
+  // need an occasional half neighbor list
+
+  int irequest = neighbor->request((void *) this);
+  neighbor->requests[irequest]->pair = 0;
+  neighbor->requests[irequest]->compute = 1;
+  neighbor->requests[irequest]->occasional = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePairLocal::init_list(int id, NeighList *ptr)
+{
+  list = ptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePairLocal::compute_local()
+{
+  invoked_local = update->ntimestep;
+
+  // count local entries and compute pair info
+
+  ncount = compute_pairs(0);
+  if (ncount > nmax) reallocate(ncount);
+  size_local_rows = ncount;
+  ncount = compute_pairs(1);
+}
+
+/* ----------------------------------------------------------------------
+   count pairs and compute pair info on this proc
+   only count pair once if newton_pair is off
+   both atom I,J must be in group
+   if flag is set, compute requested info about pair
+------------------------------------------------------------------------- */
+
+int ComputePairLocal::compute_pairs(int flag)
+{
+  int i,j,m,n,ii,jj,inum,jnum,itype,jtype;
+  double xtmp,ytmp,ztmp,delx,dely,delz;
+  double rsq,eng,fpair,factor_coul,factor_lj;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  double *dbuf,*ebuf,*fbuf;
+
+  double **x = atom->x;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  int nall = nlocal + atom->nghost;
+  double *special_coul = force->special_coul;
+  double *special_lj = force->special_lj;
+  int newton_pair = force->newton_pair;
+
+  // invoke half neighbor list (will copy or build if necessary)
+
+  if (flag == 0) neighbor->build_one(list->index);
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  // loop over neighbors of my atoms
+  // skip if I or J are not in group
+
+  if (flag) {
+    if (nvalues == 1) {
+      if (dflag >= 0) dbuf = vector;
+      if (eflag >= 0) ebuf = vector;
+      if (fflag >= 0) fbuf = vector;
+    } else {
+      if (dflag >= 0) dbuf = &array[0][dflag];
+      if (eflag >= 0) ebuf = &array[0][eflag];
+      if (fflag >= 0) fbuf = &array[0][fflag];
+    }
+  }
+
+  Pair *pair = force->pair;
+  double **cutsq = force->pair->cutsq;
+
+  m = n = 0;
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    if (!(mask[i] & groupbit)) continue;
+
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+
+      if (j < nall) factor_coul = factor_lj = 1.0;
+      else {
+	factor_coul = special_coul[j/nall];
+	factor_lj = special_lj[j/nall];
+	j %= nall;
+      }
+
+      if (!(mask[j] & groupbit)) continue;
+      if (newton_pair == 0 && j >= nlocal) continue;
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx*delx + dely*dely + delz*delz;
+      jtype = type[j];
+      if (rsq >= cutsq[itype][jtype]) continue;
+	
+      if (flag) {
+	if (dflag >= 0) dbuf[n] = sqrt(rsq);
+	if (eflag >= 0 || fflag >= 0) {
+	  eng = pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);
+	  if (eflag >= 0) ebuf[n] = eng;
+	  if (fflag >= 0) fbuf[n] = fpair;
+	}
+	n += nvalues;
+      }
+
+      m++;
+    }
+  }
+
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePairLocal::reallocate(int n)
+{
+  // grow vector or array and indices array
+
+  while (nmax < n) nmax += DELTA;
+
+  if (nvalues == 1) {
+    memory->sfree(vector);
+    vector = (double *) memory->smalloc(nmax*sizeof(double),
+					"pair/local:vector");
+    vector_local = vector;
+  } else {
+    memory->destroy_2d_double_array(array);
+    array = memory->create_2d_double_array(nmax,nvalues,
+					   "pair/local:array");
+    array_local = array;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local data
+------------------------------------------------------------------------- */
+
+double ComputePairLocal::memory_usage()
+{
+  double bytes = nmax*nvalues * sizeof(double);
+  return bytes;
+}

src/compute_pair_local.h

@@ -0,0 +1,46 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#ifndef COMPUTE_PAIR_LOCAL_H
+#define COMPUTE_PAIR_LOCAL_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputePairLocal : public Compute {
+ public:
+  ComputePairLocal(class LAMMPS *, int, char **);
+  ~ComputePairLocal();
+  void init();
+  void init_list(int, class NeighList *);
+  void compute_local();
+  double memory_usage();
+
+ private:
+  int nvalues,dflag,eflag,fflag;
+  int ncount;
+
+  int nmax;
+  double *vector;
+  double **array;
+
+  class NeighList *list;
+
+  int compute_pairs(int);
+  void reallocate(int);
+};
+
+}
+
+#endif

src/style.h

@@ -94,6 +94,7 @@ CommandStyle(write_restart,WriteRestart)
 #include "compute_ke_atom.h"
 #include "compute_msd.h"
 #include "compute_msd_molecule.h"
+#include "compute_pair_local.h"
 #include "compute_pe.h"
 #include "compute_pe_atom.h"
 #include "compute_pressure.h"
@@ -134,6 +135,7 @@ ComputeStyle(ke,ComputeKE)
 ComputeStyle(ke/atom,ComputeKEAtom)
 ComputeStyle(msd,ComputeMSD)
 ComputeStyle(msd/molecule,ComputeMSDMolecule)
+ComputeStyle(pair/local,ComputePairLocal)
 ComputeStyle(pe,ComputePE)
 ComputeStyle(pe/atom,ComputePEAtom)
 ComputeStyle(pressure,ComputePressure)