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

Merge pull request #2887 from jrgissing/bond-react-fixes

Browse Source 
Bond react fixes
This commit is contained in:
Axel Kohlmeyer
2021-08-27 16:58:33 -04:00
committed by GitHub
parent cc69431dd3 4089d7757d
commit 3e058b734b
3 changed files with 170 additions and 43 deletions
Download Patch File Download Diff File Expand all files Collapse all files

187
src/REACTION/fix_bond_react.cpp
View File

@ -140,6 +140,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
rxnfunclist[0] = "rxnsum";
rxnfunclist[1] = "rxnave";
nvvec = 0;
ncustomvars = 0;
vvec = nullptr;
nxspecial = nullptr;
@ -232,7 +233,6 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
memory->create(molecule_keyword,nreacts,"bond/react:molecule_keyword");
memory->create(nconstraints,nreacts,"bond/react:nconstraints");
memory->create(constraintstr,nreacts,MAXLINE,"bond/react:constraintstr");
memory->create(constraints,0,nreacts,"bond/react:constraints");
memory->create(var_flag,NUMVARVALS,nreacts,"bond/react:var_flag");
memory->create(var_id,NUMVARVALS,nreacts,"bond/react:var_id");
memory->create(iatomtype,nreacts,"bond/react:iatomtype");
@ -488,6 +488,9 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
if (custom_charges_fragid[i] >= 0) CustomCharges(custom_charges_fragid[i],i);
}
// get the names of per-atom variables needed by 'rxn' functions of custom constraint
customvarnames();
// initialize Marsaglia RNG with processor-unique seed (Arrhenius prob)
rrhandom = new RanMars*[narrhenius];
@ -617,7 +620,6 @@ FixBondReact::~FixBondReact()
memory->destroy(stabilize_steps_flag);
memory->destroy(custom_charges_fragid);
memory->destroy(molecule_keyword);
memory->destroy(constraints);
memory->destroy(nconstraints);
memory->destroy(constraintstr);
memory->destroy(create_atoms_flag);
@ -1020,6 +1022,11 @@ void FixBondReact::post_integrate()
return;
}
// evaluate custom constraint variable values here and forward_comm
get_customvars();
commflag = 1;
comm->forward_comm_fix(this,ncustomvars);
// run through the superimpose algorithm
// this checks if simulation topology matches unreacted mol template
superimpose_algorithm();
@ -2086,6 +2093,100 @@ double FixBondReact::get_temperature(tagint **myglove, int row_offset, int col)
return t;
}
/* ----------------------------------------------------------------------
get per-atom variable names used by custom constraint
------------------------------------------------------------------------- */
void FixBondReact::customvarnames()
{
int pos,pos1,pos2,pos3,prev3;
std::string varstr,argstr,varid;
// search all constraints' varstr for special 'rxn' functions
// add variable names to customvarstrs
// add values to customvars
for (rxnID = 0; rxnID < nreacts; rxnID++) {
for (int i = 0; i < nconstraints[rxnID]; i++) {
if (constraints[i][rxnID].type == CUSTOM) {
varstr = constraints[i][rxnID].str;
prev3 = -1;
while (true) {
// find next reaction special function occurrence
pos1 = INT_MAX;
for (int i = 0; i < nrxnfunction; i++) {
pos = varstr.find(rxnfunclist[i],prev3+1);
if (pos == std::string::npos) continue;
if (pos < pos1) pos1 = pos;
}
if (pos1 == INT_MAX) break;
pos2 = varstr.find("(",pos1);
pos3 = varstr.find(")",pos2);
if (pos2 == std::string::npos || pos3 == std::string::npos)
error->all(FLERR,"Bond/react: Illegal rxn function syntax\n");
prev3 = pos3;
argstr = varstr.substr(pos2+1,pos3-pos2-1);
argstr.erase(remove_if(argstr.begin(), argstr.end(), isspace), argstr.end()); // remove whitespace
pos2 = argstr.find(",");
if (pos2 != std::string::npos) varid = argstr.substr(0,pos2);
else varid = argstr;
// check if we already know about this variable
int varidflag = 0;
for (int j = 0; j < ncustomvars; j++) {
if (customvarstrs[j] == varid) {
varidflag = 1;
break;
}
}
if (!varidflag) {
customvarstrs.resize(ncustomvars+1);
customvarstrs[ncustomvars++] = varid;
}
}
}
}
}
}
/* ----------------------------------------------------------------------
evaluate per-atom variables needed for custom constraint
------------------------------------------------------------------------- */
void FixBondReact::get_customvars()
{
double *tempvvec;
std::string varid;
int nall = atom->nlocal + atom->nghost;
memory->create(tempvvec,nall,"bond/react:tempvvec");
if (vvec == nullptr) {
memory->create(vvec,nall,ncustomvars,"bond/react:vvec");
nvvec = nall;
}
if (nvvec < nall) {
memory->grow(vvec,nall,ncustomvars,"bond/react:vvec");
nvvec = nall;
}
for (int i = 0; i < ncustomvars; i++) {
varid = customvarstrs[i];
if (varid.substr(0,2) != "v_") error->all(FLERR,"Bond/react: Reaction special function variable "
"name should begin with 'v_'");
varid = varid.substr(2);
int ivar = input->variable->find(varid.c_str());
if (ivar < 0)
error->all(FLERR,"Bond/react: Reaction special function variable "
"name does not exist");
if (!input->variable->atomstyle(ivar))
error->all(FLERR,"Bond/react: Reaction special function must "
"reference an atom-style variable");
input->variable->compute_atom(ivar,igroup,tempvvec,1,0);
for (int j = 0; j < nall; j++) vvec[j][i] = tempvvec[j];
}
memory->destroy(tempvvec);
}
/* ----------------------------------------------------------------------
evaulate expression for variable constraint
------------------------------------------------------------------------- */
@ -2120,6 +2221,7 @@ double FixBondReact::custom_constraint(std::string varstr)
evlstr.push_back(varstr.substr(prev3+1,pos1-(prev3+1)));
prev3 = pos3;
argstr = varstr.substr(pos2+1,pos3-pos2-1);
argstr.erase(remove_if(argstr.begin(), argstr.end(), isspace), argstr.end()); // remove whitespace
pos2 = argstr.find(",");
if (pos2 != std::string::npos) {
varid = argstr.substr(0,pos2);
@ -2145,25 +2247,19 @@ currently two 'rxn' functions: rxnsum and rxnave
double FixBondReact::rxnfunction(std::string rxnfunc, std::string varid,
std::string fragid)
{
if (varid.substr(0,2) != "v_") error->one(FLERR,"Bond/react: Reaction special function variable "
"name should begin with 'v_'");
varid = varid.substr(2);
int ivar = input->variable->find(varid.c_str());
int ivar = -1;
for (int i = 0; i < ncustomvars; i++) {
if (varid == customvarstrs[i]) {
ivar = i;
break;
}
}
// variable name should always be found, at this point
// however, let's double check for completeness
if (ivar < 0)
error->one(FLERR,"Bond/react: Reaction special function variable "
"name does not exist");
if (!input->variable->atomstyle(ivar))
error->one(FLERR,"Bond/react: Reaction special function must "
"reference an atom-style variable");
if (vvec == nullptr) {
memory->create(vvec,atom->nlocal,"bond/react:vvec");
nvvec = atom->nlocal;
}
if (nvvec < atom->nlocal) {
memory->grow(vvec,atom->nlocal,"bond/react:vvec");
nvvec = atom->nlocal;
}
input->variable->compute_atom(ivar,igroup,vvec,1,0);
int ifrag = -1;
if (fragid != "all") {
ifrag = onemol->findfragment(fragid.c_str());
@ -2178,14 +2274,14 @@ double FixBondReact::rxnfunction(std::string rxnfunc, std::string varid,
if (fragid == "all") {
for (int i = 0; i < onemol->natoms; i++) {
iatom = atom->map(glove[i][1]);
sumvvec += vvec[iatom];
sumvvec += vvec[iatom][ivar];
}
nsum = onemol->natoms;
} else {
for (int i = 0; i < onemol->natoms; i++) {
if (onemol->fragmentmask[ifrag][i]) {
iatom = atom->map(glove[i][1]);
sumvvec += vvec[iatom];
sumvvec += vvec[iatom][ivar];
nsum++;
}
}
@ -2830,6 +2926,19 @@ void FixBondReact::update_everything()
rxnID = local_mega_glove[0][i];
// reactions already shuffled from dedup procedure, so can skip first N
if (iskip[rxnID]++ < nlocalskips[rxnID]) continue;
// atoms inserted here for serial MPI_STUBS build only
if (create_atoms_flag[rxnID] == 1) {
onemol = atom->molecules[unreacted_mol[rxnID]];
twomol = atom->molecules[reacted_mol[rxnID]];
if (insert_atoms(local_mega_glove,i)) {
inserted_atoms_flag = 1;
} else { // create aborted
reaction_count_total[rxnID]--;
continue;
}
}
for (int j = 0; j < max_natoms+1; j++)
update_mega_glove[j][update_num_mega] = local_mega_glove[j][i];
update_num_mega++;
@ -2842,7 +2951,7 @@ void FixBondReact::update_everything()
// we can insert atoms here, now that reactions are finalized
// can't do it any earlier, due to skipped reactions (max_rxn)
// reactions that create atoms are always treated as 'global'
// for MPI build, reactions that create atoms are always treated as 'global'
if (create_atoms_flag[rxnID] == 1) {
onemol = atom->molecules[unreacted_mol[rxnID]];
twomol = atom->molecules[reacted_mol[rxnID]];
@ -2858,17 +2967,18 @@ void FixBondReact::update_everything()
update_mega_glove[j][update_num_mega] = global_mega_glove[j][i];
update_num_mega++;
}
// if inserted atoms and global map exists, reset map now instead
// of waiting for comm since other pre-exchange fixes may use it
// invoke map_init() b/c atom count has grown
// do this once after all atom insertions
if (inserted_atoms_flag == 1 && atom->map_style != Atom::MAP_NONE) {
atom->map_init();
atom->map_set();
}
}
delete [] iskip;
// if inserted atoms and global map exists, reset map now instead
// of waiting for comm since other pre-exchange fixes may use it
// invoke map_init() b/c atom count has grown
// do this once after all atom insertions
if (inserted_atoms_flag == 1 && atom->map_style != Atom::MAP_NONE) {
atom->map_init();
atom->map_set();
}
// mark to-delete atoms
nlocal = atom->nlocal;
if (nlocal > nmark) {
@ -3668,7 +3778,7 @@ void FixBondReact::read(int myrxn)
else if (strstr(line,"constraints")) {
sscanf(line,"%d",&nconstraints[myrxn]);
if (maxnconstraints < nconstraints[myrxn]) maxnconstraints = nconstraints[myrxn];
memory->grow(constraints,maxnconstraints,nreacts,"bond/react:constraints");
constraints.resize(maxnconstraints, std::vector<Constraint>(nreacts));
} else break;
}
@ -4027,6 +4137,15 @@ int FixBondReact::pack_forward_comm(int n, int *list, double *buf,
m = 0;
if (commflag == 1) {
for (i = 0; i < n; i++) {
j = list[i];
for (k = 0; k < ncustomvars; k++)
buf[m++] = vvec[j][k];
}
return m;
}
if (commflag == 2) {
for (i = 0; i < n; i++) {
j = list[i];
@ -4051,12 +4170,16 @@ int FixBondReact::pack_forward_comm(int n, int *list, double *buf,
void FixBondReact::unpack_forward_comm(int n, int first, double *buf)
{
int i,j,m,ns,last;
int i,j,k,m,ns,last;
m = 0;
last = first + n;
if (commflag == 2) {
if (commflag == 1) {
for (i = first; i < last; i++)
for (k = 0; k < ncustomvars; k++)
vvec[i][k] = buf[m++];
} else if (commflag == 2) {
for (i = first; i < last; i++)
partner[i] = (tagint) ubuf(buf[m++]).i;
} else {

20
src/REACTION/fix_bond_react.h
View File

@ -139,8 +139,6 @@ class FixBondReact : public Fix {
int **delete_atoms; // atoms in pre-reacted templates to delete
int **create_atoms; // atoms in post-reacted templates to create
int ***chiral_atoms; // pre-react chiral atoms. 1) flag 2) orientation 3-4) ordered atom types
int nvvec;
double *vvec; // per-atom vector to store variable constraint atom-style variable values
int **nxspecial, **onemol_nxspecial, **twomol_nxspecial; // full number of 1-4 neighbors
tagint **xspecial, **onemol_xspecial, **twomol_xspecial; // full 1-4 neighbor list
@ -151,14 +149,12 @@ class FixBondReact : public Fix {
// for all mega_gloves and global_mega_glove: first row is the ID of bond/react
tagint **local_mega_glove; // consolidation local of reaction instances
tagint **ghostly_mega_glove; // consolidation nonlocal of reaction instances
tagint *
*global_mega_glove; // consolidation (inter-processor) of gloves containing nonlocal atoms
tagint **global_mega_glove; // consolidation (inter-processor) of gloves containing nonlocal atoms
int *localsendlist; // indicates ghosts of other procs
int local_num_mega; // num of local reaction instances
int ghostly_num_mega; // num of ghostly reaction instances
int global_megasize; // num of reaction instances in global_mega_glove
int *
pioneers; // during Superimpose Algorithm, atoms which have been assigned, but whose first neighbors haven't
int *pioneers; // during Superimpose Algorithm, atoms which have been assigned, but whose first neighbors haven't
int glove_counter; // used to determine when to terminate Superimpose Algorithm
void read(int);
@ -180,8 +176,10 @@ class FixBondReact : public Fix {
int check_constraints();
void get_IDcoords(int, int, double *);
double get_temperature(tagint **, int, int);
double custom_constraint(std::string);
double rxnfunction(std::string, std::string, std::string);
void customvarnames(); // get per-atom variables names used by custom constraint
void get_customvars(); // evaluate local values for variables names used by custom constraint
double custom_constraint(std::string); // evaulate expression for custom constraint
double rxnfunction(std::string, std::string, std::string); // eval rxn_sum and rxn_ave
int get_chirality(double[12]); // get handedness given an ordered set of coordinates
void open(char *);
@ -216,7 +214,11 @@ class FixBondReact : public Fix {
double par[MAXCONPAR];
std::string str;
};
Constraint **constraints;
int ncustomvars;
std::vector<std::string> customvarstrs;
int nvvec;
double **vvec; // per-atom vector to store variable constraint atom-style variable values
std::vector<std::vector<Constraint>> constraints;
// DEBUG

6
src/molecule.cpp
View File

@ -148,13 +148,15 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) :
if (me == 0)
utils::logmesg(lmp,"Read molecule template {}:\n {} molecules\n"
" {} fragments\n"
" {} atoms with max type {}\n"
" {} bonds with max type {}\n"
" {} angles with max type {}\n"
" {} dihedrals with max type {}\n"
" {} impropers with max type {}\n", id,nmolecules,
natoms,ntypes,nbonds,nbondtypes,nangles,nangletypes,
ndihedrals,ndihedraltypes,nimpropers,nimpropertypes);
nfragments,natoms,ntypes,nbonds,nbondtypes,nangles,
nangletypes,ndihedrals,ndihedraltypes,nimpropers,
nimpropertypes);
}
/* ---------------------------------------------------------------------- */