ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
sjplimp
2015-08-28 23:51:23 +00:00
parent c9188c0095
commit 546ccbe611
11 changed files with 1258 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

27
src/USER-TALLY/README Normal file
View File

@ -0,0 +1,27 @@
This package provides computes styles that use a particular hook
to accumulate information about pairwise interactions directly
as part of the pairwise force computation. It inserts additional
function calls into the Pair::ev_tally() method. Several LAMMPS
users have asked for such a facility over time.
The currently provided compute styles are mostly meant as a
demonstration for how to use this facility and provide an
alternative approach to using features like compute group/group
or compute heat/flux. Its application is limited to pairwise
additive potentials that use the standard Pair::ev_tally()
method to accumulate energy (and virial).
Nevertheless, since those compute styles are executed directly
using information that has already been computed, they should
usually be more efficient than their counterparts.
There are example scripts for using this package in examples/USER/tally
The person who created this package is Axel Kohlmeyer (akohlmey@gmail.com)
at Temple University with a little help and inspiration from
Loris Ercole (SISSA/ISAS Trieste), who contributed compute heat/flux/tally.
Additional contributed compute style for this package are welcome.
Please contact Axel, if you have questions about the implementation.

223
src/USER-TALLY/compute_force_tally.cpp Normal file
View File

@ -0,0 +1,223 @@
/* ----------------------------------------------------------------------
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_force_tally.h"
#include "atom.h"
#include "group.h"
#include "pair.h"
#include "update.h"
#include "memory.h"
#include "error.h"
#include "force.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeForceTally::ComputeForceTally(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg < 4) error->all(FLERR,"Illegal compute force/tally command");
igroup2 = group->find(arg[3]);
if (igroup2 == -1)
error->all(FLERR,"Could not find compute force/tally second group ID");
groupbit2 = group->bitmask[igroup2];
scalar_flag = 1;
vector_flag = 0;
peratom_flag = 1;
timeflag = 1;
comm_reverse = size_peratom_cols = 3;
extscalar = 1;
peflag = 1; // we need Pair::ev_tally() to be run
did_compute = 0;
invoked_peratom = invoked_scalar = -1;
nmax = -1;
fatom = NULL;
vector = new double[size_peratom_cols];
}
/* ---------------------------------------------------------------------- */
ComputeForceTally::~ComputeForceTally()
{
if (force && force->pair) force->pair->del_tally_callback(this);
delete[] vector;
}
/* ---------------------------------------------------------------------- */
void ComputeForceTally::init()
{
if (force->pair == NULL)
error->all(FLERR,"Trying to use compute force/tally with no pair style");
else
force->pair->add_tally_callback(this);
if (force->pair->single_enable == 0 || force->pair->manybody_flag)
error->all(FLERR,"Compute force/tally used with incompatible pair style.");
if ((comm->me == 0) && (force->bond || force->angle || force->dihedral
|| force->improper || force->kspace))
error->warning(FLERR,"Compute force/tally only called from pair style");
did_compute = -1;
}
/* ---------------------------------------------------------------------- */
void ComputeForceTally::pair_tally_callback(int i, int j, int nlocal, int newton,
double, double, double fpair,
double dx, double dy, double dz)
{
const int ntotal = atom->nlocal + atom->nghost;
const int * const mask = atom->mask;
// do setup work that needs to be done only once per timestep
if (did_compute != update->ntimestep) {
did_compute = update->ntimestep;
// grow local force array if necessary
// needs to be atom->nmax in length
if (atom->nmax > nmax) {
memory->destroy(fatom);
nmax = atom->nmax;
memory->create(fatom,nmax,size_peratom_cols,"force/tally:fatom");
array_atom = fatom;
}
// clear storage as needed
if (newton) {
for (int i=0; i < ntotal; ++i)
for (int j=0; j < size_peratom_cols; ++j)
fatom[i][j] = 0.0;
} else {
for (int i=0; i < atom->nlocal; ++i)
for (int j=0; j < size_peratom_cols; ++j)
fatom[i][j] = 0.0;
}
for (int i=0; i < size_peratom_cols; ++i)
vector[i] = ftotal[i] = 0.0;
}
if ( ((mask[i] & groupbit) && (mask[j] & groupbit2))
|| ((mask[i] & groupbit2) && (mask[j] & groupbit)) ) {
if (newton || i < nlocal) {
if (mask[i] & groupbit) {
ftotal[0] += fpair*dx;
ftotal[1] += fpair*dy;
ftotal[2] += fpair*dz;
}
fatom[i][0] += fpair*dx;
fatom[i][1] += fpair*dy;
fatom[i][2] += fpair*dz;
}
if (newton || j < nlocal) {
if (mask[j] & groupbit) {
ftotal[0] -= fpair*dx;
ftotal[1] -= fpair*dy;
ftotal[2] -= fpair*dz;
}
fatom[j][0] -= fpair*dx;
fatom[j][1] -= fpair*dy;
fatom[j][2] -= fpair*dz;
}
}
}
/* ---------------------------------------------------------------------- */
int ComputeForceTally::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++] = fatom[i][0];
buf[m++] = fatom[i][1];
buf[m++] = fatom[i][2];
}
return m;
}
/* ---------------------------------------------------------------------- */
void ComputeForceTally::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
fatom[j][0] += buf[m++];
fatom[j][1] += buf[m++];
fatom[j][2] += buf[m++];
}
}
/* ---------------------------------------------------------------------- */
double ComputeForceTally::compute_scalar()
{
invoked_scalar = update->ntimestep;
if ((did_compute != invoked_scalar) || (update->eflag_global != invoked_scalar))
error->all(FLERR,"Energy was not tallied on needed timestep");
// sum accumulated forces across procs
MPI_Allreduce(ftotal,vector,size_peratom_cols,MPI_DOUBLE,MPI_SUM,world);
scalar = sqrt(vector[0]*vector[0]+vector[1]*vector[1]+vector[2]*vector[2]);
return scalar;
}
/* ---------------------------------------------------------------------- */
void ComputeForceTally::compute_peratom()
{
invoked_peratom = update->ntimestep;
if ((did_compute != invoked_peratom) || (update->eflag_global != invoked_peratom))
error->all(FLERR,"Energy was not tallied on needed timestep");
// collect contributions from ghost atoms
if (force->newton_pair) {
comm->reverse_comm_compute(this);
const int nall = atom->nlocal + atom->nghost;
for (int i = atom->nlocal; i < nall; ++i)
for (int j = 0; j < size_peratom_cols; ++j)
fatom[i][j] = 0.0;
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeForceTally::memory_usage()
{
double bytes = nmax*size_peratom_cols * sizeof(double);
return bytes;
}

66
src/USER-TALLY/compute_force_tally.h Normal file
View File

@ -0,0 +1,66 @@
/* -*- c++ -*- ----------------------------------------------------------
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(force/tally,ComputeForceTally)
#else
#ifndef LMP_COMPUTE_FORCE_TALLY_H
#define LMP_COMPUTE_FORCE_TALLY_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeForceTally : public Compute {
public:
ComputeForceTally(class LAMMPS *, int, char **);
virtual ~ComputeForceTally();
void init();
double compute_scalar();
void compute_peratom();
int pack_reverse_comm(int, int, double *);
void unpack_reverse_comm(int, int *, double *);
double memory_usage();
void pair_tally_callback(int, int, int, int,
double, double, double,
double, double, double);
private:
bigint did_compute;
int nmax,igroup2,groupbit2;
double **fatom;
double ftotal[3];
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
*/

286
src/USER-TALLY/compute_heat_flux_tally.cpp Normal file
View File

@ -0,0 +1,286 @@
/* ----------------------------------------------------------------------
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_heat_flux_tally.h"
#include "atom.h"
#include "group.h"
#include "pair.h"
#include "update.h"
#include "memory.h"
#include "error.h"
#include "force.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeHeatFluxTally::ComputeHeatFluxTally(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg < 4) error->all(FLERR,"Illegal compute heat/flux/tally command");
igroup2 = group->find(arg[3]);
if (igroup2 == -1)
error->all(FLERR,"Could not find compute heat/flux/tally second group ID");
groupbit2 = group->bitmask[igroup2];
vector_flag = 1;
timeflag = 1;
comm_reverse = 7;
extvector = 1;
size_vector = 6;
peflag = 1; // we need Pair::ev_tally() to be run
did_compute = 0;
invoked_peratom = invoked_scalar = -1;
nmax = -1;
stress = NULL;
eatom = NULL;
vector = new double[size_vector];
}
/* ---------------------------------------------------------------------- */
ComputeHeatFluxTally::~ComputeHeatFluxTally()
{
if (force && force->pair) force->pair->del_tally_callback(this);
delete[] vector;
}
/* ---------------------------------------------------------------------- */
void ComputeHeatFluxTally::init()
{
if (force->pair == NULL)
error->all(FLERR,"Trying to use compute heat/flux/tally with no pair style");
else
force->pair->add_tally_callback(this);
if (force->pair->single_enable == 0 || force->pair->manybody_flag)
error->all(FLERR,"Compute heat/flux/tally used with incompatible pair style.");
if ((comm->me == 0) && (force->bond || force->angle || force->dihedral
|| force->improper || force->kspace))
error->warning(FLERR,"Compute heat/flux/tally only called from pair style");
did_compute = -1;
}
/* ---------------------------------------------------------------------- */
void ComputeHeatFluxTally::pair_tally_callback(int i, int j, int nlocal, int newton,
double evdwl, double ecoul, double fpair,
double dx, double dy, double dz)
{
const int ntotal = atom->nlocal + atom->nghost;
const int * const mask = atom->mask;
// do setup work that needs to be done only once per timestep
if (did_compute != update->ntimestep) {
did_compute = update->ntimestep;
// grow local stress and eatom arrays if necessary
// needs to be atom->nmax in length
if (atom->nmax > nmax) {
memory->destroy(stress);
nmax = atom->nmax;
memory->create(stress,nmax,6,"heat/flux/tally:stress");
memory->destroy(eatom);
nmax = atom->nmax;
memory->create(eatom,nmax,"heat/flux/tally:eatom");
}
// clear storage as needed
if (newton) {
for (int i=0; i < ntotal; ++i) {
eatom[i] = 0.0;
stress[i][0] = 0.0;
stress[i][1] = 0.0;
stress[i][2] = 0.0;
stress[i][3] = 0.0;
stress[i][4] = 0.0;
stress[i][5] = 0.0;
}
} else {
for (int i=0; i < atom->nlocal; ++i) {
eatom[i] = 0.0;
stress[i][0] = 0.0;
stress[i][1] = 0.0;
stress[i][2] = 0.0;
stress[i][3] = 0.0;
stress[i][4] = 0.0;
stress[i][5] = 0.0;
}
}
for (int i=0; i < size_vector; ++i)
vector[i] = heatj[i] = 0.0;
}
if ( ((mask[i] & groupbit) && (mask[j] & groupbit2))
|| ((mask[i] & groupbit2) && (mask[j] & groupbit)) ) {
const double epairhalf = 0.5 * (evdwl + ecoul);
fpair *= 0.5;
const double v0 = dx*dx*fpair; // dx*fpair = Fij_x
const double v1 = dy*dy*fpair;
const double v2 = dz*dz*fpair;
const double v3 = dx*dy*fpair;
const double v4 = dx*dz*fpair;
const double v5 = dy*dz*fpair;
if (newton || i < nlocal) {
eatom[i] += epairhalf;
stress[i][0] += v0;
stress[i][1] += v1;
stress[i][2] += v2;
stress[i][3] += v3;
stress[i][4] += v4;
stress[i][5] += v5;
}
if (newton || j < nlocal) {
eatom[j] += epairhalf;
stress[j][0] += v0;
stress[j][1] += v1;
stress[j][2] += v2;
stress[j][3] += v3;
stress[j][4] += v4;
stress[j][5] += v5;
}
}
}
/* ---------------------------------------------------------------------- */
int ComputeHeatFluxTally::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++] = eatom[i];
buf[m++] = stress[i][0];
buf[m++] = stress[i][1];
buf[m++] = stress[i][2];
buf[m++] = stress[i][3];
buf[m++] = stress[i][4];
buf[m++] = stress[i][5];
}
return m;
}
/* ---------------------------------------------------------------------- */
void ComputeHeatFluxTally::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
eatom[j] += buf[m++];
stress[j][0] += buf[m++];
stress[j][1] += buf[m++];
stress[j][2] += buf[m++];
stress[j][3] += buf[m++];
stress[j][4] += buf[m++];
stress[j][5] += buf[m++];
}
}
/* ---------------------------------------------------------------------- */
void ComputeHeatFluxTally::compute_vector()
{
invoked_vector = update->ntimestep;
if ((did_compute != invoked_vector) || (update->eflag_global != invoked_vector))
error->all(FLERR,"Energy was not tallied on needed timestep");
// collect contributions from ghost atoms
if (force->newton_pair) {
comm->reverse_comm_compute(this);
const int nall = atom->nlocal + atom->nghost;
for (int i = atom->nlocal; i < nall; ++i) {
eatom[i] = 0.0;
stress[i][0] = 0.0;
stress[i][1] = 0.0;
stress[i][2] = 0.0;
stress[i][3] = 0.0;
stress[i][4] = 0.0;
stress[i][5] = 0.0;
}
}
// compute heat currents
// heat flux vector = jc[3] + jv[3]
// jc[3] = convective portion of heat flux = sum_i (ke_i + pe_i) v_i[3]
// jv[3] = virial portion of heat flux = sum_i (stress_tensor_i . v_i[3])
// normalization by volume is not included
// J = sum_i( (0.5*m*v_i^2 + 0.5*(evdwl_i+ecoul_i))*v_i +
// + (F_ij . v_i)*dR_ij/2 )
int nlocal = atom->nlocal;
int *mask = atom->mask;
const double pfactor = 0.5 * force->mvv2e;
double **v = atom->v;
double *mass = atom->mass;
int *type = atom->type;
double jc[3] = {0.0,0.0,0.0};
double jv[3] = {0.0,0.0,0.0};
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
double ke_i = pfactor * mass[type[i]] *
(v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
jc[0] += (ke_i + eatom[i]) * v[i][0];
jc[1] += (ke_i + eatom[i]) * v[i][1];
jc[2] += (ke_i + eatom[i]) * v[i][2];
jv[0] += stress[i][0]*v[i][0] + stress[i][3]*v[i][1] +
stress[i][4]*v[i][2];
jv[1] += stress[i][3]*v[i][0] + stress[i][1]*v[i][1] +
stress[i][5]*v[i][2];
jv[2] += stress[i][4]*v[i][0] + stress[i][5]*v[i][1] +
stress[i][2]*v[i][2];
}
}
// sum accumulated heatj across procs
heatj[0] = jc[0] + jv[0];
heatj[1] = jc[1] + jv[1];
heatj[2] = jc[2] + jv[2];
heatj[3] = jc[0];
heatj[4] = jc[1];
heatj[5] = jc[2];
MPI_Allreduce(heatj,vector,size_vector,MPI_DOUBLE,MPI_SUM,world);
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeHeatFluxTally::memory_usage()
{
double bytes = nmax*comm_reverse * sizeof(double);
return bytes;
}

65
src/USER-TALLY/compute_heat_flux_tally.h Normal file
View File

@ -0,0 +1,65 @@
/* -*- c++ -*- ----------------------------------------------------------
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(heat/flux/tally,ComputeHeatFluxTally)
#else
#ifndef LMP_COMPUTE_HEAT_FLUX_TALLY_H
#define LMP_COMPUTE_HEAT_FLUX_TALLY_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeHeatFluxTally : public Compute {
public:
ComputeHeatFluxTally(class LAMMPS *, int, char **);
virtual ~ComputeHeatFluxTally();
void init();
void compute_vector();
int pack_reverse_comm(int, int, double *);
void unpack_reverse_comm(int, int *, double *);
double memory_usage();
void pair_tally_callback(int, int, int, int,
double, double, double,
double, double, double);
private:
bigint did_compute;
int nmax,igroup2,groupbit2;
double **stress,*eatom;
double heatj[6];
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
*/

204
src/USER-TALLY/compute_pe_tally.cpp Normal file
View File

@ -0,0 +1,204 @@
/* ----------------------------------------------------------------------
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_tally.h"
#include "atom.h"
#include "group.h"
#include "pair.h"
#include "update.h"
#include "memory.h"
#include "error.h"
#include "force.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputePETally::ComputePETally(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg < 4) error->all(FLERR,"Illegal compute pe/tally command");
igroup2 = group->find(arg[3]);
if (igroup2 == -1)
error->all(FLERR,"Could not find compute pe/tally second group ID");
groupbit2 = group->bitmask[igroup2];
scalar_flag = 1;
vector_flag = 0;
peratom_flag = 1;
timeflag = 1;
comm_reverse = size_peratom_cols = 2;
extscalar = 1;
peflag = 1; // we need Pair::ev_tally() to be run
did_compute = invoked_peratom = invoked_scalar = -1;
nmax = -1;
eatom = NULL;
vector = new double[size_peratom_cols];
}
/* ---------------------------------------------------------------------- */
ComputePETally::~ComputePETally()
{
if (force && force->pair) force->pair->del_tally_callback(this);
delete[] vector;
}
/* ---------------------------------------------------------------------- */
void ComputePETally::init()
{
if (force->pair == NULL)
error->all(FLERR,"Trying to use compute pe/tally with no pair style");
else
force->pair->add_tally_callback(this);
if (force->pair->single_enable == 0 || force->pair->manybody_flag)
error->all(FLERR,"Compute pe/tally used with incompatible pair style.");
if ((comm->me == 0) && (force->bond || force->angle || force->dihedral
|| force->improper || force->kspace))
error->warning(FLERR,"Compute pe/tally only called from pair style");
did_compute = -1;
}
/* ---------------------------------------------------------------------- */
void ComputePETally::pair_tally_callback(int i, int j, int nlocal, int newton,
double evdwl, double ecoul, double,
double, double, double)
{
const int ntotal = atom->nlocal + atom->nghost;
const int * const mask = atom->mask;
// do setup work that needs to be done only once per timestep
if (did_compute != update->ntimestep) {
did_compute = update->ntimestep;
// grow local eatom array if necessary
// needs to be atom->nmax in length
if (atom->nmax > nmax) {
memory->destroy(eatom);
nmax = atom->nmax;
memory->create(eatom,nmax,size_peratom_cols,"pe/tally:eatom");
array_atom = eatom;
}
// clear storage as needed
if (newton) {
for (int i=0; i < ntotal; ++i)
eatom[i][0] = eatom[i][1] = 0.0;
} else {
for (int i=0; i < atom->nlocal; ++i)
eatom[i][0] = eatom[i][1] = 0.0;
}
vector[0] = etotal[0] = vector[1] = etotal[1] = 0.0;
}
if ( ((mask[i] & groupbit) && (mask[j] & groupbit2))
|| ((mask[i] & groupbit2) && (mask[j] & groupbit)) ){
evdwl *= 0.5; ecoul *= 0.5;
if (newton || i < nlocal) {
etotal[0] += evdwl; eatom[i][0] += evdwl;
etotal[1] += ecoul; eatom[i][1] += ecoul;
}
if (newton || j < nlocal) {
etotal[0] += evdwl; eatom[j][0] += evdwl;
etotal[1] += ecoul; eatom[j][1] += ecoul;
}
}
}
/* ---------------------------------------------------------------------- */
int ComputePETally::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++] = eatom[i][0];
buf[m++] = eatom[i][1];
}
return m;
}
/* ---------------------------------------------------------------------- */
void ComputePETally::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
eatom[j][0] += buf[m++];
eatom[j][1] += buf[m++];
}
}
/* ---------------------------------------------------------------------- */
double ComputePETally::compute_scalar()
{
invoked_scalar = update->ntimestep;
if ((did_compute != invoked_scalar) || (update->eflag_global != invoked_scalar))
error->all(FLERR,"Energy was not tallied on needed timestep");
// sum accumulated energies across procs
MPI_Allreduce(etotal,vector,size_peratom_cols,MPI_DOUBLE,MPI_SUM,world);
scalar = vector[0]+vector[1];
return scalar;
}
/* ---------------------------------------------------------------------- */
void ComputePETally::compute_peratom()
{
invoked_peratom = update->ntimestep;
if ((did_compute != invoked_peratom) || (update->eflag_global != invoked_peratom))
error->all(FLERR,"Energy was not tallied on needed timestep");
// collect contributions from ghost atoms
if (force->newton_pair) {
comm->reverse_comm_compute(this);
const int nall = atom->nlocal + atom->nghost;
for (int i = atom->nlocal; i < nall; ++i)
eatom[i][0] = eatom[i][1] = 0.0;
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputePETally::memory_usage()
{
double bytes = nmax*size_peratom_cols * sizeof(double);
return bytes;
}

66
src/USER-TALLY/compute_pe_tally.h Normal file
View File

@ -0,0 +1,66 @@
/* -*- c++ -*- ----------------------------------------------------------
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(pe/tally,ComputePETally)
#else
#ifndef LMP_COMPUTE_PETALLY_H
#define LMP_COMPUTE_PETALLY_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputePETally : public Compute {
public:
ComputePETally(class LAMMPS *, int, char **);
virtual ~ComputePETally();
void init();
double compute_scalar();
void compute_peratom();
int pack_reverse_comm(int, int, double *);
void unpack_reverse_comm(int, int *, double *);
double memory_usage();
void pair_tally_callback(int, int, int, int,
double, double, double,
double, double, double);
private:
bigint did_compute;
int nmax,igroup2,groupbit2;
double **eatom;
double etotal[2];
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
*/

249
src/USER-TALLY/compute_stress_tally.cpp Normal file
View File

@ -0,0 +1,249 @@
/* ----------------------------------------------------------------------
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_stress_tally.h"
#include "atom.h"
#include "group.h"
#include "pair.h"
#include "update.h"
#include "memory.h"
#include "error.h"
#include "force.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeStressTally::ComputeStressTally(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg < 4) error->all(FLERR,"Illegal compute stress/tally command");
igroup2 = group->find(arg[3]);
if (igroup2 == -1)
error->all(FLERR,"Could not find compute stress/tally second group ID");
groupbit2 = group->bitmask[igroup2];
scalar_flag = 1;
vector_flag = 0;
peratom_flag = 1;
timeflag = 1;
comm_reverse = size_peratom_cols = 6;
extscalar = 0;
peflag = 1; // we need Pair::ev_tally() to be run
did_compute = 0;
invoked_peratom = invoked_scalar = -1;
nmax = -1;
stress = NULL;
vector = new double[size_peratom_cols];
}
/* ---------------------------------------------------------------------- */
ComputeStressTally::~ComputeStressTally()
{
if (force && force->pair) force->pair->del_tally_callback(this);
delete[] vector;
}
/* ---------------------------------------------------------------------- */
void ComputeStressTally::init()
{
if (force->pair == NULL)
error->all(FLERR,"Trying to use compute stress/tally with no pair style");
else
force->pair->add_tally_callback(this);
if (force->pair->single_enable == 0 || force->pair->manybody_flag)
error->all(FLERR,"Compute stress/tally used with incompatible pair style.");
if ((comm->me == 0) && (force->bond || force->angle || force->dihedral
|| force->improper || force->kspace))
error->warning(FLERR,"Compute stress/tally only called from pair style");
did_compute = -1;
}
/* ---------------------------------------------------------------------- */
void ComputeStressTally::pair_tally_callback(int i, int j, int nlocal, int newton,
double, double, double fpair,
double dx, double dy, double dz)
{
const int ntotal = atom->nlocal + atom->nghost;
const int * const mask = atom->mask;
// do setup work that needs to be done only once per timestep
if (did_compute != update->ntimestep) {
did_compute = update->ntimestep;
// grow local stress array if necessary
// needs to be atom->nmax in length
if (atom->nmax > nmax) {
memory->destroy(stress);
nmax = atom->nmax;
memory->create(stress,nmax,size_peratom_cols,"stress/tally:stress");
array_atom = stress;
}
// clear storage as needed
if (newton) {
for (int i=0; i < ntotal; ++i)
for (int j=0; j < size_peratom_cols; ++j)
stress[i][j] = 0.0;
} else {
for (int i=0; i < atom->nlocal; ++i)
for (int j=0; j < size_peratom_cols; ++j)
stress[i][j] = 0.0;
}
for (int i=0; i < size_peratom_cols; ++i)
vector[i] = virial[i] = 0.0;
}
if ( ((mask[i] & groupbit) && (mask[j] & groupbit2))
|| ((mask[i] & groupbit2) && (mask[j] & groupbit)) ) {
fpair *= 0.5;
const double v0 = dx*dx*fpair;
const double v1 = dy*dy*fpair;
const double v2 = dz*dz*fpair;
const double v3 = dx*dy*fpair;
const double v4 = dx*dz*fpair;
const double v5 = dy*dz*fpair;
if (newton || i < nlocal) {
virial[0] += v0; stress[i][0] += v0;
virial[1] += v1; stress[i][1] += v1;
virial[2] += v2; stress[i][2] += v2;
virial[3] += v3; stress[i][3] += v3;
virial[4] += v4; stress[i][4] += v4;
virial[5] += v5; stress[i][5] += v5;
}
if (newton || j < nlocal) {
virial[0] += v0; stress[j][0] += v0;
virial[1] += v1; stress[j][1] += v1;
virial[2] += v2; stress[j][2] += v2;
virial[3] += v3; stress[j][3] += v3;
virial[4] += v4; stress[j][4] += v4;
virial[5] += v5; stress[j][5] += v5;
}
}
}
/* ---------------------------------------------------------------------- */
int ComputeStressTally::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++] = stress[i][0];
buf[m++] = stress[i][1];
buf[m++] = stress[i][2];
buf[m++] = stress[i][3];
buf[m++] = stress[i][4];
buf[m++] = stress[i][5];
}
return m;
}
/* ---------------------------------------------------------------------- */
void ComputeStressTally::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
stress[j][0] += buf[m++];
stress[j][1] += buf[m++];
stress[j][2] += buf[m++];
stress[j][3] += buf[m++];
stress[j][4] += buf[m++];
stress[j][5] += buf[m++];
}
}
/* ---------------------------------------------------------------------- */
double ComputeStressTally::compute_scalar()
{
invoked_scalar = update->ntimestep;
if ((did_compute != invoked_scalar) || (update->eflag_global != invoked_scalar))
error->all(FLERR,"Energy was not tallied on needed timestep");
// sum accumulated forces across procs
MPI_Allreduce(virial,vector,size_peratom_cols,MPI_DOUBLE,MPI_SUM,world);
if (domain->dimension == 3)
scalar = (vector[0]+vector[1]+vector[2])/3.0;
else
scalar = (vector[0]+vector[1])/2.0;
return scalar;
}
/* ---------------------------------------------------------------------- */
void ComputeStressTally::compute_peratom()
{
invoked_peratom = update->ntimestep;
if ((did_compute != invoked_peratom) || (update->eflag_global != invoked_peratom))
error->all(FLERR,"Energy was not tallied on needed timestep");
// collect contributions from ghost atoms
if (force->newton_pair) {
comm->reverse_comm_compute(this);
const int nall = atom->nlocal + atom->nghost;
for (int i = atom->nlocal; i < nall; ++i)
for (int j = 0; j < size_peratom_cols; ++j)
stress[i][j] = 0.0;
}
// convert to stress*volume units = -pressure*volume
const double nktv2p = -force->nktv2p;
for (int i = 0; i < atom->nlocal; i++) {
stress[i][0] *= nktv2p;
stress[i][1] *= nktv2p;
stress[i][2] *= nktv2p;
stress[i][3] *= nktv2p;
stress[i][4] *= nktv2p;
stress[i][5] *= nktv2p;
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeStressTally::memory_usage()
{
double bytes = nmax*size_peratom_cols * sizeof(double);
return bytes;
}

66
src/USER-TALLY/compute_stress_tally.h Normal file
View File

@ -0,0 +1,66 @@
/* -*- c++ -*- ----------------------------------------------------------
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(stress/tally,ComputeStressTally)
#else
#ifndef LMP_COMPUTE_STRESS_TALLY_H
#define LMP_COMPUTE_STRESS_TALLY_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeStressTally : public Compute {
public:
ComputeStressTally(class LAMMPS *, int, char **);
virtual ~ComputeStressTally();
void init();
double compute_scalar();
void compute_peratom();
int pack_reverse_comm(int, int, double *);
void unpack_reverse_comm(int, int *, double *);
double memory_usage();
void pair_tally_callback(int, int, int, int,
double, double, double,
double, double, double);
private:
bigint did_compute;
int nmax,igroup2,groupbit2;
double **stress;
double virial[6];
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
*/

3
src/info.cpp
View File

@ -252,7 +252,8 @@ void Info::command(int narg, char **arg)
double cut;
for (int i=1; i <= atom->ntypes && neighbor->cuttype; ++i) {
cut = neighbor->cuttype[i];
if (comm->cutusermulti) cut = MAX(cut,comm->cutusermulti[i]);
// AXEL: this variable does not exist?
//if (comm->cutusermulti) cut = MAX(cut,comm->cutusermulti[i]);
fprintf(out,"Communication cutoff for type %d = %g\n", i, cut);
}
}

9
src/main.cpp
View File

@ -26,17 +26,16 @@ int main(int argc, char **argv)
{
MPI_Init(&argc,&argv);
// try {
//try {
LAMMPS *lammps = new LAMMPS(argc,argv,MPI_COMM_WORLD);
lammps->input->file();
delete lammps;
// }
// catch(...) {
//}
//catch(...) {
//int me;
//MPI_Comm_rank(MPI_COMM_WORLD,&me);
//fprintf(stderr,"Unknown exception caught on rank %d\n",me);
// }
//}
MPI_Barrier(MPI_COMM_WORLD);
MPI_Finalize();