From 35eb67c68ef8fa77544f93bebc7e9af854c736ca Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 1 Oct 2013 14:12:06 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10842 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-REAXC/compute_spec_atom.cpp | 646 +++++++++++++++++++++++++++ src/USER-REAXC/compute_spec_atom.h | 104 +++++ 2 files changed, 750 insertions(+) create mode 100644 src/USER-REAXC/compute_spec_atom.cpp create mode 100644 src/USER-REAXC/compute_spec_atom.h diff --git a/src/USER-REAXC/compute_spec_atom.cpp b/src/USER-REAXC/compute_spec_atom.cpp new file mode 100644 index 0000000000..0e7cf1ff8d --- /dev/null +++ b/src/USER-REAXC/compute_spec_atom.cpp @@ -0,0 +1,646 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Labo0ratories + 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_spec_atom.h" +#include "math_extra.h" +#include "atom.h" +#include "update.h" +#include "force.h" +#include "domain.h" +#include "memory.h" +#include "error.h" + +#include "reaxc_defs.h" +#include "reaxc_types.h" +#include "pair_reax_c.h" + +using namespace LAMMPS_NS; + +enum{KEYWORD,COMPUTE,FIX,VARIABLE}; + +/* ---------------------------------------------------------------------- */ + +ComputeSpecAtom::ComputeSpecAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg < 4) error->all(FLERR,"Illegal compute reax/c/atom command"); + + peratom_flag = 1; + nvalues = narg - 3; + if (nvalues == 1) size_peratom_cols = 0; + else size_peratom_cols = nvalues; + + // Initiate reaxc + reaxc = (PairReaxC *) force->pair_match("reax/c",1); + + pack_choice = new FnPtrPack[nvalues]; + + int i; + for (int iarg = 3; iarg < narg; iarg++) { + i = iarg-3; + + // standard lammps attributes + if (strcmp(arg[iarg],"q") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_q; + + } else if (strcmp(arg[iarg],"x") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_x; + } else if (strcmp(arg[iarg],"y") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_y; + } else if (strcmp(arg[iarg],"z") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_z; + + } else if (strcmp(arg[iarg],"vx") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_vx; + } else if (strcmp(arg[iarg],"vy") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_vy; + } else if (strcmp(arg[iarg],"vz") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_vz; + + // from pair_reax_c + } else if (strcmp(arg[iarg],"abo01") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo01; + } else if (strcmp(arg[iarg],"abo02") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo02; + } else if (strcmp(arg[iarg],"abo03") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo03; + } else if (strcmp(arg[iarg],"abo04") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo04; + } else if (strcmp(arg[iarg],"abo05") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo05; + } else if (strcmp(arg[iarg],"abo06") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo06; + } else if (strcmp(arg[iarg],"abo07") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo07; + } else if (strcmp(arg[iarg],"abo08") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo08; + } else if (strcmp(arg[iarg],"abo09") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo09; + } else if (strcmp(arg[iarg],"abo10") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo10; + } else if (strcmp(arg[iarg],"abo11") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo11; + } else if (strcmp(arg[iarg],"abo12") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo12; + } else if (strcmp(arg[iarg],"abo13") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo13; + } else if (strcmp(arg[iarg],"abo14") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo14; + } else if (strcmp(arg[iarg],"abo15") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo15; + } else if (strcmp(arg[iarg],"abo16") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo16; + } else if (strcmp(arg[iarg],"abo17") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo17; + } else if (strcmp(arg[iarg],"abo18") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo18; + } else if (strcmp(arg[iarg],"abo19") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo19; + } else if (strcmp(arg[iarg],"abo20") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo20; + } else if (strcmp(arg[iarg],"abo21") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo21; + } else if (strcmp(arg[iarg],"abo22") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo22; + } else if (strcmp(arg[iarg],"abo23") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo23; + } else if (strcmp(arg[iarg],"abo24") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_abo24; + + } else error->all(FLERR,"Invalid keyword in compute reax/c/atom command"); + } + + nmax = 0; + vector = NULL; + array = NULL; + +} + +/* ---------------------------------------------------------------------- */ + +ComputeSpecAtom::~ComputeSpecAtom() +{ + delete [] pack_choice; + memory->destroy(vector); + memory->destroy(array); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::compute_peratom() +{ + invoked_peratom = update->ntimestep; + + // grow vector or array if necessary + + if (atom->nlocal > nmax) { + nmax = atom->nmax; + if (nvalues == 1) { + memory->destroy(vector); + memory->create(vector,nmax,"property/atom:vector"); + vector_atom = vector; + } else { + memory->destroy(array); + memory->create(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 { + if (nmax > 0) { + buf = &array[0][0]; + for (int n = 0; n < nvalues; n++) + (this->*pack_choice[n])(n); + } + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double ComputeSpecAtom::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 ComputeSpecAtom::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 ComputeSpecAtom::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 ComputeSpecAtom::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 ComputeSpecAtom::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 ComputeSpecAtom::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 ComputeSpecAtom::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 ComputeSpecAtom::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 ComputeSpecAtom::pack_abo01(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][0]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo02(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][1]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo03(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][2]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo04(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][3]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo05(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][4]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo06(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][5]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo07(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][6]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo08(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][7]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo09(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][8]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo10(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][9]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo11(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][10]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo12(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][11]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo13(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][12]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo14(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][13]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo15(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][14]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo16(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][15]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo17(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][16]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo18(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][17]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo19(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][18]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo20(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][19]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo21(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][20]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo22(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][21]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo23(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][22]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_abo24(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpbo[i][23]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ diff --git a/src/USER-REAXC/compute_spec_atom.h b/src/USER-REAXC/compute_spec_atom.h new file mode 100644 index 0000000000..50e4643167 --- /dev/null +++ b/src/USER-REAXC/compute_spec_atom.h @@ -0,0 +1,104 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Labo0ratories + 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(SPEC/ATOM,ComputeSpecAtom) + +#else + +#ifndef LMP_COMPUTE_SPEC_ATOM_H +#define LMP_COMPUTE_SPEC_ATOM_H + +#include "compute.h" +#include "pointers.h" + +namespace LAMMPS_NS { + +class ComputeSpecAtom : public Compute { + public: + ComputeSpecAtom(class LAMMPS *, int, char **); + ~ComputeSpecAtom(); + void init() {} + void compute_peratom(); + double memory_usage(); + + private: + int nvalues; + int nmax; + double *vector; + double **array; + double *buf; + double *vbuf; + + typedef void (ComputeSpecAtom::*FnPtrPack)(int); + FnPtrPack *pack_choice; + + void pack_q(int); + void pack_x(int); + void pack_y(int); + void pack_z(int); + void pack_vx(int); + void pack_vy(int); + void pack_vz(int); + + void pack_abo01(int); + void pack_abo02(int); + void pack_abo03(int); + void pack_abo04(int); + void pack_abo05(int); + void pack_abo06(int); + void pack_abo07(int); + void pack_abo08(int); + void pack_abo09(int); + void pack_abo10(int); + void pack_abo11(int); + void pack_abo12(int); + void pack_abo13(int); + void pack_abo14(int); + void pack_abo15(int); + void pack_abo16(int); + void pack_abo17(int); + void pack_abo18(int); + void pack_abo19(int); + void pack_abo20(int); + void pack_abo21(int); + void pack_abo22(int); + void pack_abo23(int); + void pack_abo24(int); + + class PairReaxC *reaxc; +}; + +} + +#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. + +E: Compute reaxc/atom for atom reaxc that isn't allocated + +Self-explanatory. + +E: Invalid keyword in compute reaxc/atom command + +Self-explanatory. + +*/