diff --git a/src/compute_pe_atom.cpp b/src/compute_pe_atom.cpp new file mode 100755 index 0000000000..dcea56ddb0 --- /dev/null +++ b/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; +} diff --git a/src/compute_pe_atom.h b/src/compute_pe_atom.h new file mode 100755 index 0000000000..13065aaca6 --- /dev/null +++ b/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