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

src/compute_pe_atom.cpp

@@ -0,0 +1,179 @@
+/* ----------------------------------------------------------------------
+   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_pe_atom.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "pair.h"
+#include "bond.h"
+#include "angle.h"
+#include "dihedral.h"
+#include "improper.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputePEAtom::ComputePEAtom(LAMMPS *lmp, int narg, char **arg) : 
+  Compute(lmp, narg, arg)
+{
+  if (narg < 3) error->all("Illegal compute pe/atom command");
+
+  peratom_flag = 1;
+  size_peratom = 0;
+  peatomflag = 1;
+  timeflag = 1;
+  comm_reverse = 1;
+
+  if (narg == 3) {
+    pairflag = 1;
+    bondflag = angleflag = dihedralflag = improperflag = 1;
+  } else {
+    pairflag = 0;
+    bondflag = angleflag = dihedralflag = improperflag = 0;
+    int iarg = 3;
+    while (iarg < narg) {
+      if (strcmp(arg[iarg],"pair") == 0) pairflag = 1;
+      else if (strcmp(arg[iarg],"bond") == 0) bondflag = 1;
+      else if (strcmp(arg[iarg],"angle") == 0) angleflag = 1;
+      else if (strcmp(arg[iarg],"dihedral") == 0) dihedralflag = 1;
+      else if (strcmp(arg[iarg],"improper") == 0) improperflag = 1;
+      else error->all("Illegal compute pe/atom command");
+      iarg++;
+    }
+  }
+
+  nmax = 0;
+  energy = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputePEAtom::~ComputePEAtom()
+{
+  memory->sfree(energy);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePEAtom::compute_peratom()
+{
+  int i;
+
+  invoked = 1;
+
+  // grow local energy array if necessary
+
+  if (atom->nmax > nmax) {
+    memory->sfree(energy);
+    nmax = atom->nmax;
+    energy = (double *) 
+      memory->smalloc(nmax*sizeof(double),"compute/pe/atom:energy");
+    scalar_atom = energy;
+  }
+
+  // npair includes ghosts if either newton flag is set
+  //   b/c some bonds/dihedrals call pair::ev_tally with pairwise info
+  // nbond includes ghosts if newton_bond is set
+  // ntotal includes ghosts if either newton flag is set
+
+  int nlocal = atom->nlocal;
+  int npair = nlocal;
+  int nbond = nlocal;
+  int ntotal = nlocal;
+  if (force->newton) npair += atom->nghost;
+  if (force->newton_bond) nbond += atom->nghost;
+  if (force->newton) ntotal += atom->nghost;
+
+  // clear local energy array
+
+  for (i = 0; i < ntotal; i++) energy[i] = 0.0;
+
+  // add in per-atom contributions from each force
+  
+  if (pairflag && force->pair) {
+    double *eatom = force->pair->eatom;
+    for (i = 0; i < npair; i++) energy[i] += eatom[i];
+  }
+
+  if (bondflag && force->bond) {
+    double *eatom = force->bond->eatom;
+    for (i = 0; i < nbond; i++) energy[i] += eatom[i];
+  }
+
+  if (angleflag && force->angle) {
+    double *eatom = force->angle->eatom;
+    for (i = 0; i < nbond; i++) energy[i] += eatom[i];
+  }
+
+  if (dihedralflag && force->dihedral) {
+    double *eatom = force->dihedral->eatom;
+    for (i = 0; i < nbond; i++) energy[i] += eatom[i];
+  }
+
+  if (improperflag && force->improper) {
+    double *eatom = force->improper->eatom;
+    for (i = 0; i < nbond; i++) energy[i] += eatom[i];
+  }
+
+  // communicate ghost energy between neighbor procs
+
+  if (force->newton) comm->reverse_comm_compute(this);
+
+  // zero energy of atoms not in group
+  // only do this after comm since ghost contributions must be included
+
+  int *mask = atom->mask;
+
+  for (i = 0; i < nlocal; i++)
+    if (!(mask[i] & groupbit)) energy[i] = 0.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int ComputePEAtom::pack_reverse_comm(int n, int first, double *buf)
+{
+  int i,m,last;
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) buf[m++] = energy[i];
+  return 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePEAtom::unpack_reverse_comm(int n, int *list, double *buf)
+{
+  int i,j,m;
+
+  m = 0;
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    energy[j] += buf[m++];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double ComputePEAtom::memory_usage()
+{
+  double bytes = nmax * sizeof(double);
+  return bytes;
+}

src/compute_pe_atom.h

@@ -0,0 +1,39 @@
+/* ----------------------------------------------------------------------
+   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_PE_ATOM_H
+#define COMPUTE_PE_ATOM_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputePEAtom : public Compute {
+ public:
+  ComputePEAtom(class LAMMPS *, int, char **);
+  ~ComputePEAtom();
+  void init() {}
+  void compute_peratom();
+  int pack_reverse_comm(int, int, double *);
+  void unpack_reverse_comm(int, int *, double *);
+  double memory_usage();
+
+ private:
+  int pairflag,bondflag,angleflag,dihedralflag,improperflag;
+  int nmax;
+  double *energy;
+};
+
+}
+
+#endif