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

Added fix_species

Browse Source 
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9515 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
athomps
2013-02-15 21:28:35 +00:00
parent 7ebdff8bbd
commit d01f873370
7 changed files with 1646 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
src/USER-REAXC/Install.sh
View File

@ -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

692
src/USER-REAXC/compute_spec_atom.cpp Normal file
View File

@ -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;
}
}
/* ---------------------------------------------------------------------- */

110
src/USER-REAXC/compute_spec_atom.h Normal file
View File

@ -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.
*/

681
src/USER-REAXC/fix_species.cpp Normal file
View File

@ -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<int> (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<int>(r+0.5);
else if (r<0.0) i = static_cast<int>(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;
}
/* ---------------------------------------------------------------------- */

85
src/USER-REAXC/fix_species.h Normal file
View File

@ -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

87
src/USER-REAXC/pair_reax_c.cpp
View File

@ -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");
}
}
}
}
/* ---------------------------------------------------------------------- */

5
src/USER-REAXC/pair_reax_c.h
View File

@ -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();
};
}