From 2472cdc89d61feb4fbfd673bb58d5ca3215dba0c Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 16 Dec 2009 16:56:56 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@3544 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/compute_property_atom.cpp | 1027 +++++++++++++++++++++++++++++++++ src/compute_property_atom.h | 91 +++ src/dump_custom.cpp | 81 +-- src/style.h | 2 + src/variable.cpp | 6 +- 5 files changed, 1164 insertions(+), 43 deletions(-) create mode 100644 src/compute_property_atom.cpp create mode 100644 src/compute_property_atom.h diff --git a/src/compute_property_atom.cpp b/src/compute_property_atom.cpp new file mode 100644 index 0000000000..94fe9b91d3 --- /dev/null +++ b/src/compute_property_atom.cpp @@ -0,0 +1,1027 @@ +/* ---------------------------------------------------------------------- + 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_property_atom.h" +#include "atom.h" +#include "update.h" +#include "domain.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +// customize by adding keyword +// same as in dump_custom.cpp + +enum{ID,MOL,TYPE,MASS, + X,Y,Z,XS,YS,ZS,XSTRI,YSTRI,ZSTRI,XU,YU,ZU,XUTRI,YUTRI,ZUTRI,IX,IY,IZ, + VX,VY,VZ,FX,FY,FZ, + Q,MUX,MUY,MUZ,RADIUS,OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ, + QUATW,QUATI,QUATJ,QUATK,TQX,TQY,TQZ}; + +/* ---------------------------------------------------------------------- */ + +ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg < 4) error->all("Illegal compute property/atom command"); + + peratom_flag = 1; + nvalues = narg - 3; + if (nvalues == 1) size_peratom_cols = 0; + else size_peratom_cols = nvalues; + + pack_choice = new FnPtrPack[nvalues]; + + int i; + for (int iarg = 3; iarg < narg; iarg++) { + i = iarg-3; + + if (strcmp(arg[iarg],"id") == 0) { + pack_choice[i] = &ComputePropertyAtom::pack_id; + } else if (strcmp(arg[iarg],"mol") == 0) { + if (!atom->molecule_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_molecule; + } else if (strcmp(arg[iarg],"type") == 0) { + pack_choice[i] = &ComputePropertyAtom::pack_type; + } else if (strcmp(arg[iarg],"mass") == 0) { + pack_choice[i] = &ComputePropertyAtom::pack_mass; + + } else if (strcmp(arg[iarg],"x") == 0) { + pack_choice[i] = &ComputePropertyAtom::pack_x; + } else if (strcmp(arg[iarg],"y") == 0) { + pack_choice[i] = &ComputePropertyAtom::pack_y; + } else if (strcmp(arg[iarg],"z") == 0) { + pack_choice[i] = &ComputePropertyAtom::pack_z; + } else if (strcmp(arg[iarg],"xs") == 0) { + if (domain->triclinic) + pack_choice[i] = &ComputePropertyAtom::pack_xs_triclinic; + else pack_choice[i] = &ComputePropertyAtom::pack_xs; + } else if (strcmp(arg[iarg],"ys") == 0) { + if (domain->triclinic) + pack_choice[i] = &ComputePropertyAtom::pack_ys_triclinic; + else pack_choice[i] = &ComputePropertyAtom::pack_ys; + } else if (strcmp(arg[iarg],"zs") == 0) { + if (domain->triclinic) + pack_choice[i] = &ComputePropertyAtom::pack_zs_triclinic; + else pack_choice[i] = &ComputePropertyAtom::pack_zs; + } else if (strcmp(arg[iarg],"xu") == 0) { + if (domain->triclinic) + pack_choice[i] = &ComputePropertyAtom::pack_xu_triclinic; + else pack_choice[i] = &ComputePropertyAtom::pack_xu; + } else if (strcmp(arg[iarg],"yu") == 0) { + if (domain->triclinic) + pack_choice[i] = &ComputePropertyAtom::pack_yu_triclinic; + else pack_choice[i] = &ComputePropertyAtom::pack_yu; + } else if (strcmp(arg[iarg],"zu") == 0) { + if (domain->triclinic) + pack_choice[i] = &ComputePropertyAtom::pack_zu_triclinic; + else pack_choice[i] = &ComputePropertyAtom::pack_zu; + } else if (strcmp(arg[iarg],"ix") == 0) { + pack_choice[i] = &ComputePropertyAtom::pack_ix; + } else if (strcmp(arg[iarg],"iy") == 0) { + pack_choice[i] = &ComputePropertyAtom::pack_iy; + } else if (strcmp(arg[iarg],"iz") == 0) { + pack_choice[i] = &ComputePropertyAtom::pack_iz; + + } else if (strcmp(arg[iarg],"vx") == 0) { + pack_choice[i] = &ComputePropertyAtom::pack_vx; + } else if (strcmp(arg[iarg],"vy") == 0) { + pack_choice[i] = &ComputePropertyAtom::pack_vy; + } else if (strcmp(arg[iarg],"vz") == 0) { + pack_choice[i] = &ComputePropertyAtom::pack_vz; + } else if (strcmp(arg[iarg],"fx") == 0) { + pack_choice[i] = &ComputePropertyAtom::pack_fx; + } else if (strcmp(arg[iarg],"fy") == 0) { + pack_choice[i] = &ComputePropertyAtom::pack_fy; + } else if (strcmp(arg[iarg],"fz") == 0) { + pack_choice[i] = &ComputePropertyAtom::pack_fz; + + } else if (strcmp(arg[iarg],"q") == 0) { + if (!atom->q_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_q; + } else if (strcmp(arg[iarg],"mux") == 0) { + if (!atom->mu_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_mux; + } else if (strcmp(arg[iarg],"muy") == 0) { + if (!atom->mu_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_muy; + } else if (strcmp(arg[iarg],"muz") == 0) { + if (!atom->mu_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_muz; + + } else if (strcmp(arg[iarg],"radius") == 0) { + if (!atom->radius_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_radius; + } else if (strcmp(arg[iarg],"omegax") == 0) { + if (!atom->omega_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_omegax; + } else if (strcmp(arg[iarg],"omegay") == 0) { + if (!atom->omega_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_omegay; + } else if (strcmp(arg[iarg],"omegaz") == 0) { + if (!atom->omega_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_omegaz; + } else if (strcmp(arg[iarg],"angmomx") == 0) { + if (!atom->angmom_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_angmomx; + } else if (strcmp(arg[iarg],"angmomy") == 0) { + if (!atom->angmom_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_angmomy; + } else if (strcmp(arg[iarg],"angmomz") == 0) { + if (!atom->angmom_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_angmomz; + + } else if (strcmp(arg[iarg],"quatw") == 0) { + if (!atom->quat_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_quatw; + } else if (strcmp(arg[iarg],"quati") == 0) { + if (!atom->quat_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_quati; + } else if (strcmp(arg[iarg],"quatj") == 0) { + if (!atom->quat_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_quatj; + } else if (strcmp(arg[iarg],"quatk") == 0) { + if (!atom->quat_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_quatk; + } else if (strcmp(arg[iarg],"tqx") == 0) { + if (!atom->torque_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_tqx; + } else if (strcmp(arg[iarg],"tqy") == 0) { + if (!atom->torque_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_tqy; + } else if (strcmp(arg[iarg],"tqz") == 0) { + if (!atom->torque_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_tqz; + + } else error->all("Invalid keyword in compute property/atom command"); + } + + nmax = 0; + vector = NULL; + array = NULL; +} + +/* ---------------------------------------------------------------------- */ + +ComputePropertyAtom::~ComputePropertyAtom() +{ + delete [] pack_choice; + memory->sfree(vector); + memory->destroy_2d_double_array(array); +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::compute_peratom() +{ + invoked_peratom = update->ntimestep; + + // grow vector or array if necessary + + if (atom->nlocal > nmax) { + nmax = atom->nmax; + if (nvalues == 1) { + memory->sfree(vector); + vector = (double *) memory->smalloc(nmax*sizeof(double), + "property/atom:vector"); + vector_atom = vector; + } else { + memory->destroy_2d_double_array(array); + array = memory->create_2d_double_array(nmax,nvalues, + "property/atom:array"); + array_atom = array; + } + } + + // fill vector or array with per-atom values + + if (nvalues == 1) { + buf = vector; + (this->*pack_choice[0])(0); + } else { + buf = array[0]; + for (int n = 0; n < nvalues; n++) + (this->*pack_choice[n])(n); + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double ComputePropertyAtom::memory_usage() +{ + double bytes = nmax*nvalues * sizeof(double); + return bytes; +} + +/* ---------------------------------------------------------------------- + one method for every keyword compute property/atom can output + the atom property is packed into buf starting at n with stride nvalues + customize a new keyword by adding a method +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_id(int n) +{ + int *tag = atom->tag; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = tag[i]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_molecule(int n) +{ + int *molecule = atom->molecule; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = molecule[i]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_type(int n) +{ + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = type[i]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_mass(int n) +{ + int *type = atom->type; + double *mass = atom->mass; + double *rmass = atom->rmass; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = rmass[i]; + else buf[n] = 0.0; + n += nvalues; + } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = mass[type[i]]; + else buf[n] = 0.0; + n += nvalues; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_x(int n) +{ + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = x[i][0]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_y(int n) +{ + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = x[i][1]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_z(int n) +{ + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = x[i][2]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_xs(int n) +{ + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double boxxlo = domain->boxlo[0]; + double invxprd = 1.0/domain->xprd; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = (x[i][0] - boxxlo) * invxprd; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_ys(int n) +{ + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double boxylo = domain->boxlo[1]; + double invyprd = 1.0/domain->yprd; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = (x[i][1] - boxylo) * invyprd; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_zs(int n) +{ + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double boxzlo = domain->boxlo[2]; + double invzprd = 1.0/domain->zprd; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = (x[i][2] - boxzlo) * invzprd; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_xs_triclinic(int n) +{ + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double *boxlo = domain->boxlo; + double *h_inv = domain->h_inv; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = h_inv[0]*(x[i][0]-boxlo[0]) + + h_inv[5]*(x[i][1]-boxlo[1]) + h_inv[4]*(x[i][2]-boxlo[2]); + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_ys_triclinic(int n) +{ + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double *boxlo = domain->boxlo; + double *h_inv = domain->h_inv; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = h_inv[1]*(x[i][1]-boxlo[1]) + h_inv[3]*(x[i][2]-boxlo[2]); + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_zs_triclinic(int n) +{ + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double *boxlo = domain->boxlo; + double *h_inv = domain->h_inv; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = h_inv[2]*(x[i][2]-boxlo[2]); + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_xu(int n) +{ + double **x = atom->x; + int *image = atom->image; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double xprd = domain->xprd; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = x[i][0] + ((image[i] & 1023) - 512) * xprd; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_yu(int n) +{ + double **x = atom->x; + int *image = atom->image; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double yprd = domain->yprd; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = x[i][1] + ((image[i] >> 10 & 1023) - 512) * yprd; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_zu(int n) +{ + double **x = atom->x; + int *image = atom->image; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double zprd = domain->zprd; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = x[i][2] + ((image[i] >> 20) - 512) * zprd; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_xu_triclinic(int n) +{ + double **x = atom->x; + int *image = atom->image; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double *h = domain->h; + int xbox,ybox,zbox; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + xbox = (image[i] & 1023) - 512; + ybox = (image[i] >> 10 & 1023) - 512; + zbox = (image[i] >> 20) - 512; + buf[n] = x[i][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox; + } else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_yu_triclinic(int n) +{ + double **x = atom->x; + int *image = atom->image; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double *h = domain->h; + int ybox,zbox; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + ybox = (image[i] >> 10 & 1023) - 512; + zbox = (image[i] >> 20) - 512; + buf[n] = x[i][1] + h[1]*ybox + h[3]*zbox; + } else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_zu_triclinic(int n) +{ + double **x = atom->x; + int *image = atom->image; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double *h = domain->h; + int zbox; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + zbox = (image[i] >> 20) - 512; + buf[n] = x[i][2] + h[2]*zbox; + } else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_ix(int n) +{ + int *image = atom->image; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = (image[i] & 1023) - 512; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_iy(int n) +{ + int *image = atom->image; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = (image[i] >> 10 & 1023) - 512; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_iz(int n) +{ + int *image = atom->image; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = (image[i] >> 20) - 512; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_vx(int n) +{ + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = v[i][0]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_vy(int n) +{ + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = v[i][1]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_vz(int n) +{ + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = v[i][2]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_fx(int n) +{ + double **f = atom->f; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = f[i][0]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_fy(int n) +{ + double **f = atom->f; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = f[i][1]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_fz(int n) +{ + double **f = atom->f; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = f[i][2]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_q(int n) +{ + double *q = atom->q; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = q[i]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_mux(int n) +{ + double **mu = atom->mu; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = mu[i][0]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_muy(int n) +{ + double **mu = atom->mu; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = mu[i][1]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_muz(int n) +{ + double **mu = atom->mu; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = mu[i][2]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_radius(int n) +{ + double *radius = atom->radius; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = radius[i]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_omegax(int n) +{ + double **omega = atom->omega; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = omega[i][0]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_omegay(int n) +{ + double **omega = atom->omega; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = omega[i][1]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_omegaz(int n) +{ + double **omega = atom->omega; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = omega[i][2]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_angmomx(int n) +{ + double **angmom = atom->angmom; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = angmom[i][0]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_angmomy(int n) +{ + double **angmom = atom->angmom; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = angmom[i][1]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_angmomz(int n) +{ + double **angmom = atom->angmom; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = angmom[i][2]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_quatw(int n) +{ + double **quat = atom->quat; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = quat[i][0]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_quati(int n) +{ + double **quat = atom->quat; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = quat[i][1]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_quatj(int n) +{ + double **quat = atom->quat; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = quat[i][2]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_quatk(int n) +{ + double **quat = atom->quat; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = quat[i][3]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_tqx(int n) +{ + double **torque = atom->torque; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = torque[i][0]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_tqy(int n) +{ + double **torque = atom->torque; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = torque[i][1]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_tqz(int n) +{ + double **torque = atom->torque; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = torque[i][2]; + else buf[n] = 0.0; + n += nvalues; + } +} diff --git a/src/compute_property_atom.h b/src/compute_property_atom.h new file mode 100644 index 0000000000..43a38bd891 --- /dev/null +++ b/src/compute_property_atom.h @@ -0,0 +1,91 @@ +/* ---------------------------------------------------------------------- + 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_PROPERTY_ATOM_H +#define COMPUTE_PROPERTY_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputePropertyAtom : public Compute { + public: + ComputePropertyAtom(class LAMMPS *, int, char **); + ~ComputePropertyAtom(); + void init() {} + void compute_peratom(); + double memory_usage(); + + private: + int nvalues; + int nmax; + double *vector; + double **array; + double *buf; + + typedef void (ComputePropertyAtom::*FnPtrPack)(int); + FnPtrPack *pack_choice; // ptrs to pack functions + + void pack_id(int); + void pack_molecule(int); + void pack_type(int); + void pack_mass(int); + + void pack_x(int); + void pack_y(int); + void pack_z(int); + void pack_xs(int); + void pack_ys(int); + void pack_zs(int); + void pack_xs_triclinic(int); + void pack_ys_triclinic(int); + void pack_zs_triclinic(int); + void pack_xu(int); + void pack_yu(int); + void pack_zu(int); + void pack_xu_triclinic(int); + void pack_yu_triclinic(int); + void pack_zu_triclinic(int); + void pack_ix(int); + void pack_iy(int); + void pack_iz(int); + + void pack_vx(int); + void pack_vy(int); + void pack_vz(int); + void pack_fx(int); + void pack_fy(int); + void pack_fz(int); + void pack_q(int); + void pack_mux(int); + void pack_muy(int); + void pack_muz(int); + void pack_radius(int); + void pack_omegax(int); + void pack_omegay(int); + void pack_omegaz(int); + void pack_angmomx(int); + void pack_angmomy(int); + void pack_angmomz(int); + void pack_quatw(int); + void pack_quati(int); + void pack_quatj(int); + void pack_quatk(int); + void pack_tqx(int); + void pack_tqy(int); + void pack_tqz(int); +}; + +} + +#endif diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 07b97fb055..a079d8d76f 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -30,7 +30,8 @@ using namespace LAMMPS_NS; -// customize by adding keyword to 1st enum +// customize by adding keyword +// same as in compute_property.cpp, also customize that command enum{ID,MOL,TYPE,MASS, X,Y,Z,XS,YS,ZS,XSTRI,YSTRI,ZSTRI,XU,YU,ZU,XUTRI,YUTRI,ZUTRI,IX,IY,IZ, @@ -392,7 +393,7 @@ int DumpCustom::count() nstride = 1; } else if (thresh_array[ithresh] == MOL) { if (!atom->molecule_flag) - error->all("Threshhold for an atom quantity that isn't allocated"); + error->all("Threshhold for an atom property that isn't allocated"); int *molecule = atom->molecule; for (i = 0; i < nlocal; i++) dchoose[i] = molecule[i]; ptr = dchoose; @@ -577,92 +578,92 @@ int DumpCustom::count() nstride = 3; } else if (thresh_array[ithresh] == Q) { if (!atom->q_flag) - error->all("Threshhold for an atom quantity that isn't allocated"); + error->all("Threshhold for an atom property that isn't allocated"); ptr = atom->q; nstride = 1; } else if (thresh_array[ithresh] == MUX) { if (!atom->mu_flag) - error->all("Threshhold for an atom quantity that isn't allocated"); + error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->mu[0][0]; nstride = 3; } else if (thresh_array[ithresh] == MUY) { if (!atom->mu_flag) - error->all("Threshhold for an atom quantity that isn't allocated"); + error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->mu[0][1]; nstride = 3; } else if (thresh_array[ithresh] == MUZ) { if (!atom->mu_flag) - error->all("Threshhold for an atom quantity that isn't allocated"); + error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->mu[0][2]; nstride = 3; } else if (thresh_array[ithresh] == RADIUS) { if (!atom->radius_flag) - error->all("Threshhold for an atom quantity that isn't allocated"); + error->all("Threshhold for an atom property that isn't allocated"); ptr = atom->radius; nstride = 1; } else if (thresh_array[ithresh] == OMEGAX) { if (!atom->omega_flag) - error->all("Threshhold for an atom quantity that isn't allocated"); + error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->omega[0][0]; nstride = 3; } else if (thresh_array[ithresh] == OMEGAY) { if (!atom->omega_flag) - error->all("Threshhold for an atom quantity that isn't allocated"); + error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->omega[0][1]; nstride = 3; } else if (thresh_array[ithresh] == OMEGAZ) { if (!atom->omega_flag) - error->all("Threshhold for an atom quantity that isn't allocated"); + error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->omega[0][2]; nstride = 3; } else if (thresh_array[ithresh] == ANGMOMX) { if (!atom->angmom_flag) - error->all("Threshhold for an atom quantity that isn't allocated"); + error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->angmom[0][0]; nstride = 3; } else if (thresh_array[ithresh] == ANGMOMY) { if (!atom->angmom_flag) - error->all("Threshhold for an atom quantity that isn't allocated"); + error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->angmom[0][1]; nstride = 3; } else if (thresh_array[ithresh] == ANGMOMZ) { if (!atom->angmom_flag) - error->all("Threshhold for an atom quantity that isn't allocated"); + error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->angmom[0][2]; nstride = 3; } else if (thresh_array[ithresh] == QUATW) { if (!atom->quat_flag) - error->all("Threshhold for an atom quantity that isn't allocated"); + error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->quat[0][0]; nstride = 4; } else if (thresh_array[ithresh] == QUATI) { if (!atom->quat_flag) - error->all("Threshhold for an atom quantity that isn't allocated"); + error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->quat[0][1]; nstride = 4; } else if (thresh_array[ithresh] == QUATJ) { if (!atom->quat_flag) - error->all("Threshhold for an atom quantity that isn't allocated"); + error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->quat[0][2]; nstride = 4; } else if (thresh_array[ithresh] == QUATK) { if (!atom->quat_flag) - error->all("Threshhold for an atom quantity that isn't allocated"); + error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->quat[0][3]; nstride = 4; } else if (thresh_array[ithresh] == TQX) { if (!atom->torque_flag) - error->all("Threshhold for an atom quantity that isn't allocated"); + error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->torque[0][0]; nstride = 3; } else if (thresh_array[ithresh] == TQY) { if (!atom->torque_flag) - error->all("Threshhold for an atom quantity that isn't allocated"); + error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->torque[0][1]; nstride = 3; } else if (thresh_array[ithresh] == TQZ) { if (!atom->torque_flag) - error->all("Threshhold for an atom quantity that isn't allocated"); + error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->torque[0][2]; nstride = 3; @@ -795,7 +796,7 @@ void DumpCustom::parse_fields(int narg, char **arg) vtype[i] = INT; } else if (strcmp(arg[iarg],"mol") == 0) { if (!atom->molecule_flag) - error->all("Dumping an atom quantity that isn't allocated"); + error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_molecule; vtype[i] = INT; } else if (strcmp(arg[iarg],"type") == 0) { @@ -869,94 +870,94 @@ void DumpCustom::parse_fields(int narg, char **arg) } else if (strcmp(arg[iarg],"q") == 0) { if (!atom->q_flag) - error->all("Dumping an atom quantity that isn't allocated"); + error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_q; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"mux") == 0) { if (!atom->mu_flag) - error->all("Dumping an atom quantity that isn't allocated"); + error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_mux; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"muy") == 0) { if (!atom->mu_flag) - error->all("Dumping an atom quantity that isn't allocated"); + error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_muy; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"muz") == 0) { if (!atom->mu_flag) - error->all("Dumping an atom quantity that isn't allocated"); + error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_muz; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"radius") == 0) { if (!atom->radius_flag) - error->all("Dumping an atom quantity that isn't allocated"); + error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_radius; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"omegax") == 0) { if (!atom->omega_flag) - error->all("Dumping an atom quantity that isn't allocated"); + error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_omegax; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"omegay") == 0) { if (!atom->omega_flag) - error->all("Dumping an atom quantity that isn't allocated"); + error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_omegay; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"omegaz") == 0) { if (!atom->omega_flag) - error->all("Dumping an atom quantity that isn't allocated"); + error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_omegaz; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"angmomx") == 0) { if (!atom->angmom_flag) - error->all("Dumping an atom quantity that isn't allocated"); + error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_angmomx; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"angmomy") == 0) { if (!atom->angmom_flag) - error->all("Dumping an atom quantity that isn't allocated"); + error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_angmomy; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"angmomz") == 0) { if (!atom->angmom_flag) - error->all("Dumping an atom quantity that isn't allocated"); + error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_angmomz; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"quatw") == 0) { if (!atom->quat_flag) - error->all("Dumping an atom quantity that isn't allocated"); + error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_quatw; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"quati") == 0) { if (!atom->quat_flag) - error->all("Dumping an atom quantity that isn't allocated"); + error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_quati; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"quatj") == 0) { if (!atom->quat_flag) - error->all("Dumping an atom quantity that isn't allocated"); + error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_quatj; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"quatk") == 0) { if (!atom->quat_flag) - error->all("Dumping an atom quantity that isn't allocated"); + error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_quatk; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"tqx") == 0) { if (!atom->torque_flag) - error->all("Dumping an atom quantity that isn't allocated"); + error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_tqx; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"tqy") == 0) { if (!atom->torque_flag) - error->all("Dumping an atom quantity that isn't allocated"); + error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_tqy; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"tqz") == 0) { if (!atom->torque_flag) - error->all("Dumping an atom quantity that isn't allocated"); + error->all("Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_tqz; vtype[i] = DOUBLE; @@ -1458,7 +1459,7 @@ void DumpCustom::pack_variable(int n) /* ---------------------------------------------------------------------- one method for every keyword dump custom can output - the atom quantity is packed into buf starting at n with stride size_one + the atom property is packed into buf starting at n with stride size_one customize a new keyword by adding a method ------------------------------------------------------------------------- */ diff --git a/src/style.h b/src/style.h index 0087a83da8..1484271809 100644 --- a/src/style.h +++ b/src/style.h @@ -93,6 +93,7 @@ CommandStyle(write_restart,WriteRestart) #include "compute_reduce.h" #include "compute_reduce_region.h" #include "compute_erotate_sphere.h" +#include "compute_property_atom.h" #include "compute_stress_atom.h" #include "compute_temp.h" #include "compute_temp_com.h" @@ -122,6 +123,7 @@ ComputeStyle(pressure,ComputePressure) ComputeStyle(reduce,ComputeReduce) ComputeStyle(reduce/region,ComputeReduceRegion) ComputeStyle(erotate/sphere,ComputeERotateSphere) +ComputeStyle(property/atom,ComputePropertyAtom) ComputeStyle(stress/atom,ComputeStressAtom) ComputeStyle(temp,ComputeTemp) ComputeStyle(temp/com,ComputeTempCOM) diff --git a/src/variable.cpp b/src/variable.cpp index 28b1c63b35..af35c248bc 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -1024,7 +1024,7 @@ double Variable::evaluate(char *str, Tree **tree) i = ptr-str+1; } - // v_name = non atom-style variable = global value + // v_name = scalar from non atom-style global scalar if (nbracket == 0 && style[ivar] != ATOM) { @@ -1039,7 +1039,7 @@ double Variable::evaluate(char *str, Tree **tree) treestack[ntreestack++] = newtree; } else argstack[nargstack++] = atof(var); - // v_name = atom-style variable + // v_name = vector from atom-style per-atom vector } else if (nbracket == 0 && style[ivar] == ATOM) { @@ -1049,7 +1049,7 @@ double Variable::evaluate(char *str, Tree **tree) double tmp = evaluate(data[ivar][0],&newtree); treestack[ntreestack++] = newtree; - // v_name[N] = global value from atom-style variable + // v_name[N] = scalar from atom-style per-atom vector // compute the per-atom variable in result // use peratom2global to extract single value from result