From 0ca04bc6b35037d74505fcc84403dd7d4c8bc3fe Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 8 May 2013 15:24:36 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9834 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-REAXC/compute_spec_atom.cpp | 696 --------------------------- src/USER-REAXC/compute_spec_atom.h | 107 ---- src/USER-REAXC/fix_reaxc_bonds.cpp | 25 +- src/USER-REAXC/fix_species.cpp | 682 -------------------------- src/USER-REAXC/fix_species.h | 76 --- src/USER-REAXC/pair_reax_c.cpp | 80 +-- 6 files changed, 14 insertions(+), 1652 deletions(-) delete mode 100644 src/USER-REAXC/compute_spec_atom.cpp delete mode 100644 src/USER-REAXC/compute_spec_atom.h delete mode 100644 src/USER-REAXC/fix_species.cpp delete mode 100644 src/USER-REAXC/fix_species.h diff --git a/src/USER-REAXC/compute_spec_atom.cpp b/src/USER-REAXC/compute_spec_atom.cpp deleted file mode 100644 index e7a988640a..0000000000 --- a/src/USER-REAXC/compute_spec_atom.cpp +++ /dev/null @@ -1,696 +0,0 @@ -/* ---------------------------------------------------------------------- - 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],"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 deleted file mode 100644 index 77f3bfc302..0000000000 --- a/src/USER-REAXC/compute_spec_atom.h +++ /dev/null @@ -1,107 +0,0 @@ -/* ---------------------------------------------------------------------- - 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_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; - struct _reax_system *system; - struct _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_reaxc_bonds.cpp b/src/USER-REAXC/fix_reaxc_bonds.cpp index 9c2fc62d99..3b9cce52d2 100644 --- a/src/USER-REAXC/fix_reaxc_bonds.cpp +++ b/src/USER-REAXC/fix_reaxc_bonds.cpp @@ -145,14 +145,15 @@ void FixReaxCBonds::Output_ReaxC_Bonds(bigint ntimestep, FILE *fp) } for (i = 0; i < nmax; i++) { - numneigh[i] = 0.0; + numneigh[i] = 0; for (j = 0; j < MAXREAXBOND; j++) { neighid[i][j] = 0; abo[i][j] = 0.0; } } + numbonds = 0; - FindBond( system, lists, numbonds); + FindBond(lists, numbonds); // allocate a temporary buffer for the snapshot info MPI_Allreduce(&numbonds,&numbonds_max,1,MPI_INT,MPI_MAX,world); @@ -162,10 +163,10 @@ void FixReaxCBonds::Output_ReaxC_Bonds(bigint ntimestep, FILE *fp) for (i = 0; i < nbuf; i ++) buf[i] = 0.0; // Pass information to buffer - PassBuffer( system, buf, nbuf_local); + PassBuffer(buf, nbuf_local); // Receive information from buffer for output - RecvBuffer( system, buf, nbuf, nbuf_local, nlocal_tot, numbonds_max); + RecvBuffer(buf, nbuf, nbuf_local, nlocal_tot, numbonds_max); memory->destroy(buf); @@ -173,8 +174,7 @@ void FixReaxCBonds::Output_ReaxC_Bonds(bigint ntimestep, FILE *fp) /* ---------------------------------------------------------------------- */ -void FixReaxCBonds::FindBond( struct _reax_system *system, - struct _reax_list *lists, int &numbonds) +void FixReaxCBonds::FindBond(struct _reax_list *lists, int &numbonds) { int *ilist, i, ii, inum; int j, pj, nj, jtag, jtype; @@ -192,7 +192,7 @@ void FixReaxCBonds::FindBond( struct _reax_system *system, for( pj = Start_Index(i, reaxc->lists); pj < End_Index(i, reaxc->lists); ++pj ) { bo_ij = &( reaxc->lists->select.bond_list[pj] ); j = bo_ij->nbr; - jtag = reaxc->system->my_atoms[j].orig_id; + jtag = atom->tag[j]; bo_tmp = bo_ij->bo_data.BO; if (bo_tmp > bo_cut) { @@ -208,8 +208,7 @@ void FixReaxCBonds::FindBond( struct _reax_system *system, } /* ---------------------------------------------------------------------- */ -void FixReaxCBonds::PassBuffer( struct _reax_system *system, double *buf, - int &nbuf_local) +void FixReaxCBonds::PassBuffer(double *buf, int &nbuf_local) { int i, j, k, jtag, numbonds; int nlocal = atom->nlocal; @@ -244,13 +243,13 @@ void FixReaxCBonds::PassBuffer( struct _reax_system *system, double *buf, /* ---------------------------------------------------------------------- */ -void FixReaxCBonds::RecvBuffer( struct _reax_system *system, double *buf, - int nbuf, int nbuf_local, int natoms, int maxnum) +void FixReaxCBonds::RecvBuffer(double *buf, int nbuf, int nbuf_local, + int natoms, int maxnum) { int i, j, k, l, itype, jtype, itag, jtag; int inode, nlocal_tmp, numbonds, molid; int nlocal = atom->nlocal; - int ntimestep = update->ntimestep; + bigint ntimestep = update->ntimestep; double sbotmp, nlptmp, avqtmp, abotmp; double cutof3 = reaxc->control->bg_cut; @@ -348,7 +347,7 @@ double FixReaxCBonds::memory_usage() { double bytes; - bytes += 3.0*nmax*sizeof(double); + bytes = 3.0*nmax*sizeof(double); bytes += nmax*sizeof(int); bytes += 1.0*nmax*MAXREAXBOND*sizeof(double); bytes += 1.0*nmax*MAXREAXBOND*sizeof(int); diff --git a/src/USER-REAXC/fix_species.cpp b/src/USER-REAXC/fix_species.cpp deleted file mode 100644 index 815c73b393..0000000000 --- a/src/USER-REAXC/fix_species.cpp +++ /dev/null @@ -1,682 +0,0 @@ -/* ---------------------------------------------------------------------- - 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" -#include "reaxc_defs.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 deleted file mode 100644 index c3267ce8ea..0000000000 --- a/src/USER-REAXC/fix_species.h +++ /dev/null @@ -1,76 +0,0 @@ -/* ---------------------------------------------------------------------- - 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" /* universal defines inside namespace */ - -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; - -}; -} - -#endif -#endif diff --git a/src/USER-REAXC/pair_reax_c.cpp b/src/USER-REAXC/pair_reax_c.cpp index 746d68c11b..d4945afee6 100644 --- a/src/USER-REAXC/pair_reax_c.cpp +++ b/src/USER-REAXC/pair_reax_c.cpp @@ -20,7 +20,7 @@ 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 + Fix reax/c/bonds and fix reax/c/species for pair_style reax/c added by Ray Shan (Sandia) ------------------------------------------------------------------------- */ @@ -104,10 +104,6 @@ PairReaxC::PairReaxC(LAMMPS *lmp) : Pair(lmp) setup_flag = 0; - fixbond_flag = fixspecies_flag = 0; - tmpr = NULL; - tmpid = NULL; - tmpbo = NULL; nmax = 0; } @@ -152,9 +148,6 @@ PairReaxC::~PairReaxC() delete [] eta; delete [] gamma; } - memory->destroy(tmpr); - memory->destroy(tmpid); - memory->destroy(tmpbo); delete [] pvector; @@ -197,6 +190,7 @@ void PairReaxC::settings(int narg, char **arg) control->hbond_cut = 7.50; control->thb_cut = 0.001; control->thb_cutsq = 0.00001; + control->bg_cut = 0.3; out_control->write_steps = 0; out_control->traj_method = 0; @@ -280,17 +274,6 @@ void PairReaxC::coeff( int nargs, char **args ) map[i-2] = -1; continue; } - - /* - itmp = atoi(args[i]) - 1; - map[i-2] = itmp; - - // error check - - if (itmp < 0 || itmp >= nreax_types) - error->all(FLERR,"Non-existent ReaxFF type"); - */ - } int n = atom->ntypes; @@ -308,7 +291,6 @@ void PairReaxC::coeff( int nargs, char **args ) if (itmp != n) error->all(FLERR,"Non-existent ReaxFF type"); - int count = 0; for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) { @@ -542,33 +524,6 @@ void PairReaxC::compute(int eflag, int vflag) Output_Results( system, control, data, &lists, out_control, mpi_data ); - if(fixbond_flag) - fixbond( 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,MAXSPECBOND,"pair:tmpr"); - memory->create(tmpid,nmax,MAXSPECBOND,"pair:tmpid"); - memory->create(tmpbo,nmax,MAXSPECBOND,"pair:tmpbo"); - } - - for (i = 0; i < system->N; i ++) - for (j = 0; j < MAXSPECBOND; j ++) { - tmpr[i][j] = tmpbo[i][j] = 0.0; - tmpid[i][j] = -1; - } - - FindBond(); - - } - } /* ---------------------------------------------------------------------- */ @@ -771,34 +726,3 @@ 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 > MAXSPECBOND) error->all(FLERR,"Increase MAXSPECBOND in fix_reaxc_species.h"); - } - } - } -} - -/* ---------------------------------------------------------------------- */