diff --git a/src/USER-REAXC/Install.sh b/src/USER-REAXC/Install.sh index 279a3f1349..c4a5051cb8 100755 --- a/src/USER-REAXC/Install.sh +++ b/src/USER-REAXC/Install.sh @@ -6,11 +6,15 @@ if (test $1 = 1) then cp -p fix_qeq_reax.cpp .. cp -p fix_reax_c.cpp .. cp -p fix_reaxc_bonds.cpp .. + cp -p fix_species.cpp .. + cp -p compute_spec_atom.cpp .. cp -p pair_reax_c.h .. cp -p fix_qeq_reax.h .. cp -p fix_reax_c.h .. cp -p fix_reaxc_bonds.h .. + cp -p fix_species.h .. + cp -p compute_spec_atom.h .. cp -p reaxc_allocate.cpp .. cp -p reaxc_basic_comm.cpp .. @@ -64,11 +68,15 @@ elif (test $1 = 0) then rm -f ../fix_qeq_reax.cpp rm -f ../fix_reax_c.cpp rm -f ../fix_reaxc_bonds.cpp + rm -f ../fix_species.cpp + rm -f ../compute_spec_atom.cpp rm -f ../pair_reax_c.h rm -f ../fix_qeq_reax.h rm -f ../fix_reax_c.h rm -f ../fix_reaxc_bonds.h + rm -f ../fix_species.h + rm -f ../compute_spec_atom.h rm -f ../reaxc_allocate.cpp rm -f ../reaxc_basic_comm.cpp diff --git a/src/USER-REAXC/compute_spec_atom.cpp b/src/USER-REAXC/compute_spec_atom.cpp new file mode 100644 index 0000000000..fb867f41e2 --- /dev/null +++ b/src/USER-REAXC/compute_spec_atom.cpp @@ -0,0 +1,692 @@ +/* ---------------------------------------------------------------------- + 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" + +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],"jid01") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_jid01; + } else if (strcmp(arg[iarg],"jid02") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_jid02; + } else if (strcmp(arg[iarg],"jid03") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_jid03; + } else if (strcmp(arg[iarg],"jid04") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_jid04; + } else if (strcmp(arg[iarg],"jid05") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_jid05; + } else if (strcmp(arg[iarg],"jid06") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_jid06; + } else if (strcmp(arg[iarg],"jid07") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_jid07; + } else if (strcmp(arg[iarg],"jid08") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_jid08; + } else if (strcmp(arg[iarg],"jid09") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_jid09; + } else if (strcmp(arg[iarg],"jid10") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_jid10; + } else if (strcmp(arg[iarg],"jid11") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_jid11; + } else if (strcmp(arg[iarg],"jid12") == 0) { + pack_choice[i] = &ComputeSpecAtom::pack_jid12; + } 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 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_jid01(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpid[i][0]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_jid02(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpid[i][1]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_jid03(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpid[i][2]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_jid04(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpid[i][3]; + else buf[n] = 0.0; + n += nvalues; + } +} + + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_jid05(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpid[i][4]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_jid06(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpid[i][5]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_jid07(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpid[i][6]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_jid08(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpid[i][7]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_jid09(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpid[i][8]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_jid10(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpid[i][9]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSpecAtom::pack_jid11(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpid[i][10]; + 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_jid12(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = reaxc->tmpid[i][11]; + 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; + } +} + +/* ---------------------------------------------------------------------- */ + diff --git a/src/USER-REAXC/compute_spec_atom.h b/src/USER-REAXC/compute_spec_atom.h new file mode 100644 index 0000000000..c7d554cc5f --- /dev/null +++ b/src/USER-REAXC/compute_spec_atom.h @@ -0,0 +1,110 @@ +/* ---------------------------------------------------------------------- + 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 "pair_reax_c.h" +#include "reaxc_types.h" +#include "reaxc_defs.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_jid01(int); + void pack_jid02(int); + void pack_jid03(int); + void pack_jid04(int); + void pack_jid05(int); + void pack_jid06(int); + void pack_jid07(int); + void pack_jid08(int); + void pack_jid09(int); + void pack_jid10(int); + void pack_jid11(int); + void pack_jid12(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); + + class PairReaxC *reaxc; + reax_system *system; + reax_atom *my_atoms; +}; + +} + +#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. + +*/ diff --git a/src/USER-REAXC/fix_species.cpp b/src/USER-REAXC/fix_species.cpp new file mode 100644 index 0000000000..d170606017 --- /dev/null +++ b/src/USER-REAXC/fix_species.cpp @@ -0,0 +1,681 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Ray Shan (Sandia, tnshan@sandia.gov) +------------------------------------------------------------------------- */ + +#include "lmptype.h" +#include "stdlib.h" +#include "math.h" +#include "atom.h" +#include "string.h" +#include "fix_ave_atom.h" +#include "fix_species.h" +#include "domain.h" +#include "update.h" +#include "pair_reax_c.h" +#include "modify.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "comm.h" +#include "force.h" +#include "compute.h" +#include "input.h" +#include "variable.h" +#include "memory.h" +#include "error.h" +#include "reaxc_list.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +FixSpecies::FixSpecies(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg < 7) error->all(FLERR,"Illegal fix species command"); + + MPI_Comm_rank(world,&me); + MPI_Comm_size(world,&nprocs); + ntypes = atom->ntypes; + + nevery = atoi(arg[3]); + nrepeat = atoi(arg[4]); + global_freq = nfreq = atoi(arg[5]); + + comm_forward = 1; + + if (nevery <= 0 || nrepeat <= 0 || nfreq <= 0) + error->all(FLERR,"Illegal fix species command"); + if (nfreq % nevery || (nrepeat-1)*nevery >= nfreq) + error->all(FLERR,"Illegal fix species command"); + + tmparg = NULL; + memory->create(tmparg,4,4,"species:tmparg"); + strcpy(tmparg[0],arg[3]); + strcpy(tmparg[1],arg[4]); + strcpy(tmparg[2],arg[5]); + + if (me == 0) { + fp = fopen(arg[6],"w"); + if (fp == NULL) { + char str[128]; + sprintf(str,"Cannot open fix species file %s",arg[6]); + error->one(FLERR,str); + } + } + + BOCut = NULL; + clusterID = NULL; + + Name = NULL; + MolName = NULL; + MolType = NULL; + NMol = NULL; + nd = NULL; + molmap = NULL; + + nmax = 0; + setupflag = 0; + + // set default bond order cutoff + int n, i, j, itype, jtype; + double bo_cut; + bg_cut = 0.30; + n = ntypes+1; + memory->create(BOCut,n,n,"/species:BOCut"); + for (i = 1; i < n; i ++) + for (j = 1; j < n; j ++) + BOCut[i][j] = bg_cut; + + // optional args + eletype = NULL; + ele = posspec = filepos = NULL; + eleflag = posflag = padflag = 0; + + int iarg = 7; + while (iarg < narg) { + + // set BO cutoff + if (strcmp(arg[iarg],"cutoff") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix species command"); + itype = atoi(arg[iarg+1]); + jtype = atoi(arg[iarg+2]); + bo_cut = atof(arg[iarg+3]); + if (itype > ntypes || jtype > ntypes) + error->all(FLERR,"Illegal fix species command"); + if (itype <= 0 || jtype <= 0) + error->all(FLERR,"Illegal fix species command"); + if (bo_cut > 1.0 || bo_cut < 0.0) + error->all(FLERR,"Illegal fix species command"); + + BOCut[itype][jtype] = bo_cut; + BOCut[jtype][itype] = bo_cut; + iarg += 4; + + // modify element type names + } else if (strcmp(arg[iarg],"element") == 0) { + if (iarg+ntypes+1 > narg) error->all(FLERR,"Illegal fix species command"); + + int nchar = 2; + eletype = (char**) malloc(ntypes*sizeof(char*)); + for (int i = 0; i < ntypes; i ++) { + if (strlen(arg[iarg+1+i]) > nchar) + error->all(FLERR,"Illegal fix species command"); + eletype[i] = (char*) malloc(nchar*sizeof(char)); + strcpy(eletype[i],arg[iarg+1+i]); + } + eleflag = 1; + iarg += ntypes + 1; + + } else error->all(FLERR,"Illegal fix species command"); + } + + if (!eleflag) { + memory->create(ele,ntypes+1,"species:ele"); + ele[0]='C'; + ele[1]='H'; + ele[2]='O'; + ele[3]='N'; + } + +} + +/* ---------------------------------------------------------------------- */ + +FixSpecies::~FixSpecies() +{ + memory->destroy(ele); + memory->destroy(BOCut); + memory->destroy(clusterID); + + memory->destroy(nd); + memory->destroy(Name); + memory->destroy(NMol); + memory->destroy(MolType); + memory->destroy(MolName); + + if (me == 0) fclose(fp); + + modify->delete_compute("SPECATOM"); + modify->delete_fix("SPECBOND"); +} + +/* ---------------------------------------------------------------------- */ + +int FixSpecies::setmask() +{ + int mask = 0; + mask |= END_OF_STEP; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixSpecies::setup(int vflag) +{ + ntotal = static_cast (atom->natoms); + memory->create(Name,ntypes,"species:Name"); + end_of_step(); +} + +/* ---------------------------------------------------------------------- */ + +void FixSpecies::init() +{ + if (atom->tag_enable == 0) + error->all(FLERR,"Cannot use fix specis unless atoms have IDs"); + + reaxc = (PairReaxC *) force->pair_match("reax/c",1); + + if (reaxc == NULL) error->all(FLERR,"Cannot use fix species without " + "pair_style reax/c"); + + reaxc->fixspecies_flag = 1; + + nvalid = update->ntimestep+nfreq; + + // request neighbor list + int irequest = neighbor->request((void *) this); + neighbor->requests[irequest]->pair = 0; + neighbor->requests[irequest]->fix = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; + neighbor->requests[irequest]->occasional = 1; + + // check if this fix has been called twice + int count = 0; + for (int i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"species") == 0) count++; + if (count > 1 && comm->me == 0) + error->warning(FLERR,"More than one fix species"); + + if (!setupflag) { + // create a compute to store properties + create_compute(); + + // create a fix to point to fix_ave_atom for averaging stored properties + create_fix(); + + setupflag = 1; + } + +} + +/* ---------------------------------------------------------------------- */ + +void FixSpecies::create_compute() +{ + int narg; + char **args; + + narg = 22; + args = new char*[narg]; + args[0] = (char *) "SPECATOM"; + args[1] = (char *) "all"; + args[2] = (char *) "spec/atom"; + args[3] = (char *) "q"; + args[4] = (char *) "x"; + args[5] = (char *) "y"; + args[6] = (char *) "z"; + args[7] = (char *) "vx"; + args[8] = (char *) "vy"; + args[9] = (char *) "vz"; + args[10] = (char *) "abo01"; + args[11] = (char *) "abo02"; + args[12] = (char *) "abo03"; + args[13] = (char *) "abo04"; + args[14] = (char *) "abo05"; + args[15] = (char *) "abo06"; + args[16] = (char *) "abo07"; + args[17] = (char *) "abo08"; + args[18] = (char *) "abo09"; + args[19] = (char *) "abo10"; + args[20] = (char *) "abo11"; + args[21] = (char *) "abo12"; + modify->add_compute(narg,args); + delete [] args; +} + +/* ---------------------------------------------------------------------- */ + +void FixSpecies::create_fix() +{ + int narg; + char **args; + + narg = 25; + args = new char*[narg]; + args[0] = (char *) "SPECBOND"; + args[1] = (char *) "all"; + args[2] = (char *) "ave/atom"; + args[3] = tmparg[0]; + args[4] = tmparg[1]; + args[5] = tmparg[2]; + args[6] = (char *) "c_SPECATOM[1]"; // q, array_atoms[i][0] + args[7] = (char *) "c_SPECATOM[2]"; // x, 1 + args[8] = (char *) "c_SPECATOM[3]"; // y, 2 + args[9] = (char *) "c_SPECATOM[4]"; // z, 3 + args[10] = (char *) "c_SPECATOM[5]"; // vx, 4 + args[11] = (char *) "c_SPECATOM[6]"; // vy, 5 + args[12] = (char *) "c_SPECATOM[7]"; // vz, 6 + args[13] = (char *) "c_SPECATOM[8]"; // abo01, 7 + args[14] = (char *) "c_SPECATOM[9]"; + args[15] = (char *) "c_SPECATOM[10]"; + args[16] = (char *) "c_SPECATOM[11]"; + args[17] = (char *) "c_SPECATOM[12]"; + args[18] = (char *) "c_SPECATOM[13]"; + args[19] = (char *) "c_SPECATOM[14]"; + args[20] = (char *) "c_SPECATOM[15]"; + args[21] = (char *) "c_SPECATOM[16]"; + args[22] = (char *) "c_SPECATOM[17]"; + args[23] = (char *) "c_SPECATOM[18]"; + args[24] = (char *) "c_SPECATOM[19]"; // abo12, 18 + modify->add_fix(narg,args); + f_SPECBOND = (FixAveAtom *) modify->fix[modify->nfix-1]; + delete [] args; + +} + +/* ---------------------------------------------------------------------- */ + +void FixSpecies::init_list(int id, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void FixSpecies::end_of_step() +{ + Output_ReaxC_Bonds(update->ntimestep,fp); + if (me == 0) fflush(fp); +} + +/* ---------------------------------------------------------------------- */ + +void FixSpecies::Output_ReaxC_Bonds(bigint ntimestep, FILE *fp) + +{ + int i, j, k, itype, jtype, itag, jtag; + int b, nbuf, nbuf_local, inode; + int nlocal_max, numbonds, numbonds_max, count; + int Nmole, Nspec; + double *bbuf; + int *ilist,*jlist,*numneigh,**firstneigh; + int ii, jj, inum, jnum; + + // point to fix_ave_atom + f_SPECBOND->end_of_step(); + + if (ntimestep != nvalid) return; + + nlocal = atom->nlocal; + nghost = atom->nghost; + nall = nlocal + nghost; + if (atom->nmax > nmax) { + nmax = atom->nmax; + memory->destroy(clusterID); + memory->create(clusterID,nmax,"species:clusterID"); + } + + Nmole = Nspec = count = 0; + + FindMolecule(); + + SortMolecule (Nmole); + + FindSpecies(Nmole, Nspec); + + if (me == 0) + WriteFormulas (Nmole, Nspec); + + nvalid += nfreq; +} + +/* ---------------------------------------------------------------------- */ + +void FixSpecies::FindMolecule () +{ + int i,j,ii,jj,inum,jnum,itype,jtype,jtag,k,ktag,itag,loop,ntot; + int change,done,anychange; + int *mask = atom->mask; + int *ilist, *jlist, *numneigh, **firstneigh; + double bo_tmp,bo_cut; + double **spec_atom = f_SPECBOND->array_atom; + + neighbor->build_one(list->index); + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (mask[i] & groupbit) clusterID[i] = atom->tag[i]; + else clusterID[i] = 0; + } + + loop = 0; + while (1) { + comm->forward_comm_fix(this); + loop ++; + + change = 0; + while (1) { + done = 1; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (!(mask[i] & groupbit)) continue; + + itype = atom->type[i]; + + for (jj = 0; jj < MAXSPECBOND; jj++) { + j = reaxc->tmpid[i][jj]; + + if (j < i) continue; + if (!(mask[j] & groupbit)) continue; + if (clusterID[i] == clusterID[j]) continue; + + jtype = atom->type[j]; + bo_cut = BOCut[itype][jtype]; + bo_tmp = spec_atom[i][jj+7]; + + if (bo_tmp > bo_cut) { + clusterID[i] = clusterID[j] = MIN(clusterID[i],clusterID[j]); + done = 0; + } + } + } + if (!done) change = 1; + if (done) break; + } + + MPI_Allreduce(&change,&anychange,1,MPI_INT,MPI_MAX,world); + if (!anychange) break; + + MPI_Allreduce(&loop,&ntot,1,MPI_INT,MPI_SUM,world); + if (ntot >= 20*nprocs) break; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixSpecies::SortMolecule(int &Nmole) +{ + memory->destroy(molmap); + molmap = NULL; + + int m, n, idlo, idhi; + int *mask =atom->mask; + int lo = ntotal; + int hi = -ntotal; + int flag = 0; + for (n = 0; n < nlocal; n++) { + if (!(mask[n] & groupbit)) continue; + if (clusterID[n] == 0) flag = 1; + lo = MIN(lo,nint(clusterID[n])); + hi = MAX(hi,nint(clusterID[n])); + } + int flagall; + MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); + if (flagall && me == 0) + error->warning(FLERR,"Atom with cluster ID = 0 included in " + "fix species group"); + MPI_Allreduce(&lo,&idlo,1,MPI_INT,MPI_MIN,world); + MPI_Allreduce(&hi,&idhi,1,MPI_INT,MPI_MAX,world); + if (idlo == ntotal) + if (me == 0) + error->warning(FLERR,"Atom with cluster ID = maxmol " + "included in fix species group"); + + int nlen = idhi - idlo + 1; + memory->create(molmap,nlen,"species:molmap"); + for (n = 0; n < nlen; n++) molmap[n] = 0; + + for (n = 0; n < nlocal; n++) { + if (!(mask[n] & groupbit)) continue; + molmap[nint(clusterID[n])-idlo] = 1; + } + + int *molmapall; + memory->create(molmapall,nlen,"species:molmapall"); + MPI_Allreduce(molmap,molmapall,nlen,MPI_INT,MPI_MAX,world); + + Nmole = 0; + for (n = 0; n < nlen; n++) { + if (molmapall[n]) molmap[n] = Nmole++; + else molmap[n] = -1; + } + memory->destroy(molmapall); + + flag = 0; + for (n = 0; n < nlocal; n++) { + if (mask[n] & groupbit) continue; + if (nint(clusterID[n]) < idlo || nint(clusterID[n]) > idhi) continue; + if (molmap[nint(clusterID[n])-idlo] >= 0) flag = 1; + } + + MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); + if (flagall && comm->me == 0) + error->warning(FLERR,"One or more cluster has atoms not in group"); + + for (n = 0; n < nlocal; n++) { + if (!(mask[n] & groupbit)) continue; + clusterID[n] = molmap[nint(clusterID[n])-idlo]+1; + } + memory->destroy(molmap); + molmap = NULL; + +} + +/* ---------------------------------------------------------------------- */ + +void FixSpecies::FindSpecies(int Nmole, int &Nspec) +{ + int inum, *ilist; + int i, j, k, l, m, n, itype, cid; + int flag_identity, flag_mol, flag_spec; + int flag_tmp; + int *mask =atom->mask; + int *Nameall, *NMolall; + + memory->destroy(MolName); + MolName = NULL; + memory->create(MolName,Nmole*(ntypes+1),"species:MolName"); + + memory->destroy(NMol); + NMol = NULL; + memory->create(NMol,Nmole,"species:NMol"); + for (m = 0; m < Nmole; m ++) + NMol[m] = 1; + + memory->create(Nameall,ntypes,"species:Nameall"); + memory->create(NMolall,Nmole,"species:NMolall"); + + for (m = 1, Nspec = 0; m <= Nmole; m ++) { + for (n = 0; n < ntypes; n ++) Name[n] = 0; + for (n = 0, flag_mol = 0; n < nlocal; n ++) { + if (!(mask[n] & groupbit)) continue; + cid = nint(clusterID[n]); + if (cid == m) { + itype = atom->type[n]-1; + Name[itype] ++; + flag_mol = 1; + } + } + MPI_Allreduce(&flag_mol,&flag_tmp,1,MPI_INT,MPI_MAX,world); + flag_mol = flag_tmp; + + MPI_Allreduce(Name,Nameall,ntypes,MPI_INT,MPI_SUM,world); + for (n = 0; n < ntypes; n++) Name[n] = Nameall[n]; + + if (flag_mol == 1) { + flag_identity = 1; + for (k = 0; k < Nspec; k ++) { + flag_spec=0; + for (l = 0; l < ntypes; l ++) + if (MolName[ntypes*k+l] != Name[l]) flag_spec = 1; + if (flag_spec == 0) NMol[k] ++; + flag_identity *= flag_spec; + } + if (Nspec == 0 || flag_identity == 1) { + for (l = 0; l < ntypes; l ++) + MolName[ntypes*Nspec+l] = Name[l]; + Nspec ++; + } + } + } + memory->destroy(NMolall); + memory->destroy(Nameall); + + memory->destroy(nd); + nd = NULL; + memory->create(nd,Nspec,"species:nd"); + + memory->destroy(MolType); + MolType = NULL; + memory->create(MolType,Nspec*(ntypes+2),"species:MolType"); +} + +/* ---------------------------------------------------------------------- */ + +int FixSpecies::CheckExistence(int id, int ntypes) +{ + int i, j, molid, flag, num; + + for (i = 0; i < Nmoltype; i ++) { + flag = 0; + for (j = 0; j < ntypes; j ++) { + molid = MolType[ntypes * i + j]; + if (molid != MolName[ntypes * id + j]) flag = 1; + } + if (flag == 0) return i; + } + for (i = 0; i < ntypes; i ++) + MolType[ntypes * Nmoltype + i] = MolName[ntypes *id + i]; + + Nmoltype ++; + return Nmoltype - 1; +} + +/* ---------------------------------------------------------------------- */ + +void FixSpecies::WriteFormulas(int Nmole, int Nspec) +{ + int i, j, k, l, jj, itemp; + int inode; + bigint ntimestep = update->ntimestep; + + fprintf(fp,"# Timestep No_Moles No_Specs "); + + Nmoltype = 0; + + for (i = 0; i < Nspec; i ++) + nd[i] = CheckExistence(i, ntypes); + + for (i = 0; i < Nmoltype; i ++) { + for (j = 0;j < ntypes; j ++) { + itemp = MolType[ntypes * i + j]; + if (itemp != 0) { + if (eletype) fprintf(fp,"%s",eletype[j]); + else fprintf(fp,"%c",ele[j]); + if (itemp != 1) fprintf(fp,"%d",itemp); + } + } + fprintf(fp,"\t"); + } + fprintf(fp,"\n"); + + fprintf(fp,"%11d",ntimestep); + fprintf(fp,"%11d%11d\t",Nmole,Nspec); + + for (i = 0; i < Nmoltype; i ++) + fprintf(fp," %d\t",NMol[i]); + fprintf(fp,"\n"); + +} + +/* ---------------------------------------------------------------------- */ + +int FixSpecies::nint(const double &r) +{ + int i = 0; + if (r>0.0) i = static_cast(r+0.5); + else if (r<0.0) i = static_cast(r-0.5); + return i; +} + +/* ---------------------------------------------------------------------- */ + +int FixSpecies::pack_comm(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = clusterID[j]; + } + return 1; +} + +/* ---------------------------------------------------------------------- */ + +void FixSpecies::unpack_comm(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) clusterID[i] = nint(buf[m++]); +} + +/* ---------------------------------------------------------------------- */ + +double FixSpecies::memory_usage() +{ + double bytes; + + bytes += nmax*sizeof(double); + bytes += nmax*nall*sizeof(double); + + return bytes; +} + +/* ---------------------------------------------------------------------- */ diff --git a/src/USER-REAXC/fix_species.h b/src/USER-REAXC/fix_species.h new file mode 100644 index 0000000000..d34980bf00 --- /dev/null +++ b/src/USER-REAXC/fix_species.h @@ -0,0 +1,85 @@ +/* ---------------------------------------------------------------------- + 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 FIX_CLASS + +FixStyle(species,FixSpecies) + +#else + +#ifndef LMP_FIX_SPECIES_H +#define LMP_FIX_SPECIES_H + +#include "fix.h" +#include "pointers.h" + +#include "pair_reax_c.h" +#include "reaxc_types.h" +#include "reaxc_defs.h" +// #include "pair_foo.h" + +#define MAXSPECBOND 12 +#define BUFLEN 1000 + +namespace LAMMPS_NS { + +class FixSpecies : public Fix { + public: + FixSpecies(class LAMMPS *, int, char **); + ~FixSpecies(); + int setmask(); + void init(); + void init_list(int, class NeighList *); + void setup(int); + void end_of_step(); + + private: + int eleflag, posflag, multipos, padflag, setupflag; + int me, nprocs, nmax, nlocal, nghost, nall, ntypes, ntotal; + int nrepeat, nfreq, Nmoltype; + int *Name, *MolName, *NMol, *nd, *MolType, *molmap, *clusterID; + + double bg_cut, **BOCut; + + char *ele, **eletype; + char **tmparg; + char *posspec, *filepos; + + FILE *fp, *pos; + + void Output_ReaxC_Bonds(bigint, FILE *); + void create_compute(); + void create_fix(); + void FindMolecule(); + void SortMolecule(int &); + void FindSpecies(int, int &); + void WriteFormulas(int, int); + int CheckExistence(int, int); + + int nint(const double &); + int pack_comm(int, int *, double *, int, int *); + void unpack_comm(int, int, double *); + double memory_usage(); + + bigint nvalid; + + class NeighList *list; + class FixAveAtom *f_SPECBOND; + class PairReaxC *reaxc; + // class PairFoo *foo; + +}; +} + +#endif +#endif diff --git a/src/USER-REAXC/pair_reax_c.cpp b/src/USER-REAXC/pair_reax_c.cpp index 4f42028f41..67058da7d7 100644 --- a/src/USER-REAXC/pair_reax_c.cpp +++ b/src/USER-REAXC/pair_reax_c.cpp @@ -14,11 +14,14 @@ /* ---------------------------------------------------------------------- Contributing author: Hasan Metin Aktulga, Purdue University (now at Lawrence Berkeley National Laboratory, hmaktulga@lbl.gov) - Please cite the related publication: H. M. Aktulga, J. C. Fogarty, S. A. Pandit, A. Y. Grama, "Parallel Reactive Molecular Dynamics: Numerical Methods and Algorithmic Techniques", Parallel Computing, in press. + --- + Per-atom energy/virial added by Ray Shan (Sandia) + Fix reax/c/bonds and fix species for pair_style reax/c added by + Ray Shan (Sandia) ------------------------------------------------------------------------- */ #include "pair_reax_c.h" @@ -102,6 +105,10 @@ PairReaxC::PairReaxC(LAMMPS *lmp) : Pair(lmp) setup_flag = 0; fixbond_flag = fixspecies_flag = 0; + tmpr = NULL; + tmpid = NULL; + tmpbo = NULL; + nmax = 0; } /* ---------------------------------------------------------------------- */ @@ -121,13 +128,10 @@ PairReaxC::~PairReaxC() Delete_List( lists+BONDS, world ); Delete_List( lists+THREE_BODIES, world ); Delete_List( lists+FAR_NBRS, world ); - // fprintf( stderr, "3\n" ); - DeAllocate_Workspace( control, workspace ); DeAllocate_System( system ); } - //fprintf( stderr, "4\n" ); memory->destroy( system ); memory->destroy( control ); @@ -136,7 +140,6 @@ PairReaxC::~PairReaxC() memory->destroy( lists ); memory->destroy( out_control ); memory->destroy( mpi_data ); - //fprintf( stderr, "5\n" ); // deallocate interface storage if( allocated ) { @@ -149,10 +152,12 @@ PairReaxC::~PairReaxC() delete [] eta; delete [] gamma; } + memory->destroy(tmpr); + memory->destroy(tmpid); + memory->destroy(tmpbo); delete [] pvector; - //fprintf( stderr, "6\n" ); } /* ---------------------------------------------------------------------- */ @@ -438,13 +443,6 @@ void PairReaxC::compute(int eflag, int vflag) if (eflag || vflag) ev_setup(eflag,vflag); else ev_unset(); -/* if ((eflag_atom || vflag_atom) && firstwarn) { - firstwarn = 0; - if (comm->me == 0) - error->warning(FLERR,"Pair reax/c cannot yet compute " - "per-atom energy or stress"); - } */ - if (vflag_global) control->virial = 1; else control->virial = 0; @@ -522,13 +520,6 @@ void PairReaxC::compute(int eflag, int vflag) if (vflag_fdotr) virial_fdotr_compute(); -// #if defined(LOG_PERFORMANCE) -// if( comm->me == 0 && fix_qeq != NULL ) { -// data->timing.s_matvecs += fix_qeq->matvecs; -// data->timing.qEq += fix_qeq->qeq_time; -// } -// #endif - // Set internal timestep counter to that of LAMMPS data->step = update->ntimestep; @@ -538,8 +529,29 @@ void PairReaxC::compute(int eflag, int vflag) if(fixbond_flag) fixbond( system, control, data, &lists, out_control, mpi_data ); - if(fixspecies_flag) - fixspecies( system, control, data, &lists, out_control, mpi_data ); + // populate tmpr and tmpbo arrays for fix reax/c/species + int i, j; + + if(fixspecies_flag) { + if (system->N > nmax) { + memory->destroy(tmpr); + memory->destroy(tmpid); + memory->destroy(tmpbo); + nmax = system->N; + memory->create(tmpr,nmax,MAXBOND,"pair:tmpr"); + memory->create(tmpid,nmax,MAXBOND,"pair:tmpid"); + memory->create(tmpbo,nmax,MAXBOND,"pair:tmpbo"); + } + + for (i = 0; i < system->N; i ++) + for (j = 0; j < MAXBOND; j ++) { + tmpr[i][j] = tmpbo[i][j] = 0.0; + tmpid[i][j] = -1; + } + + FindBond(); + + } } @@ -743,3 +755,34 @@ void *PairReaxC::extract(const char *str, int &dim) } /* ---------------------------------------------------------------------- */ + +void PairReaxC::FindBond() +{ + int i, ii, j, pj, jtag, nj, jtmp, jj; + double bo_tmp, bo_cut, rij, rsq, r_tmp; + + bond_data *bo_ij; + bo_cut = 0.20; + + for (i = 0; i < system->n; i++) { + nj = 0; + for( pj = Start_Index(i, lists); pj < End_Index(i, lists); ++pj ) { + bo_ij = &( lists->select.bond_list[pj] ); + j = bo_ij->nbr; + if (j < i) continue; + + bo_tmp = bo_ij->bo_data.BO; + r_tmp = bo_ij->d; + + if (bo_tmp >= bo_cut ) { + tmpr[i][nj] = r_tmp; + tmpid[i][nj] = j; + tmpbo[i][nj] = bo_tmp; + nj ++; + if (nj > MAXBOND) error->all(FLERR,"Increase MAXSPECBOND in fix_reaxc_species.h"); + } + } + } +} + +/* ---------------------------------------------------------------------- */ diff --git a/src/USER-REAXC/pair_reax_c.h b/src/USER-REAXC/pair_reax_c.h index 30b1d4e14a..060ef2ce8e 100644 --- a/src/USER-REAXC/pair_reax_c.h +++ b/src/USER-REAXC/pair_reax_c.h @@ -46,6 +46,8 @@ class PairReaxC : public Pair { double init_one(int, int); void *extract(const char *, int &); int fixbond_flag, fixspecies_flag; + int **tmpid; + double ** tmpbo, **tmpr; control_params *control; reax_system *system; @@ -73,6 +75,9 @@ class PairReaxC : public Pair { int write_reax_lists(); void read_reax_forces(); void setup(); + + int nmax; + void FindBond(); }; }