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

Merge pull request #2504 from jrgissing/bond/react-add-logic-to-constraints

Browse Source 
Bond/react: add logic to constraints
This commit is contained in:
Axel Kohlmeyer
2020-12-23 15:39:30 -05:00
committed by GitHub
parent 06d7e5ab02 a7fa4739f5
commit 41d6648ac9
3 changed files with 262 additions and 189 deletions
Download Patch File Download Diff File Expand all files Collapse all files

21
doc/src/fix_bond_react.rst
View File

@ -457,6 +457,23 @@ example, the molecule fragment could consist of only the backbone
atoms of a polymer chain. This constraint can be used to enforce a
specific relative position and orientation between reacting molecules.
By default, all constraints must be satisfied for the reaction to
occur. In other words, constraints are evaluated as a series of
logical values using the logical AND operator "&&". More complex logic
can be achieved by explicitly adding the logical AND operator "&&" or
the logical OR operator "||" after a given constraint command. If a
logical operator is specified after a constraint, it must be placed
after all constraint parameters, on the same line as the constraint
(one per line). Similarly, parentheses can be used to group
constraints. The expression that results from concatenating all
constraints should be a valid logical expression that can be read by
the :doc:`variable <variable>` command after converting each
constraint to a logical value. Because exactly one constraint is
allowed per line, having a valid logical expression implies that left
parentheses "(" should only appear before a constraint, and right
parentheses ")" should only appear after a constraint and before any
logical operator.
Once a reaction site has been successfully identified, data structures
within LAMMPS that store bond topology are updated to reflect the
post-reacted molecule template. All force fields with fixed bonds,
@ -599,8 +616,8 @@ reset_mol_ids = yes, custom_charges = no, molecule = off
.. _Gissinger:
**(Gissinger)** Gissinger, Jensen and Wise, Polymer, 128, 211 (2017).
**(Gissinger)** Gissinger, Jensen and Wise, Polymer, 128, 211-217 (2017).
.. _Gissinger2020:
**(Gissinger)** Gissinger, Jensen and Wise, Macromolecules (2020, in press).
**(Gissinger)** Gissinger, Jensen and Wise, Macromolecules, 53, 22, 9953–9961 (2020).

410
src/USER-REACTION/fix_bond_react.cpp
View File

@ -82,6 +82,9 @@ enum{ACCEPT,REJECT,PROCEED,CONTINUE,GUESSFAIL,RESTORE};
// types of available reaction constraints
enum{DISTANCE,ANGLE,DIHEDRAL,ARRHENIUS,RMSD};
// ID type used by constraint
enum{ATOM,FRAG};
// keyword values that accept variables as input
enum{NEVERY,RMIN,RMAX,PROB};
@ -115,7 +118,6 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
global_freq = 1;
extvector = 0;
rxnID = 0;
nconstraints = 0;
narrhenius = 0;
status = PROCEED;
@ -206,7 +208,9 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
memory->create(stabilize_steps_flag,nreacts,"bond/react:stabilize_steps_flag");
memory->create(custom_charges_fragid,nreacts,"bond/react:custom_charges_fragid");
memory->create(molecule_keyword,nreacts,"bond/react:molecule_keyword");
memory->create(constraints,1,MAXCONARGS,"bond/react:constraints");
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");
@ -227,6 +231,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
stabilize_steps_flag[i] = 0;
custom_charges_fragid[i] = -1;
molecule_keyword[i] = OFF;
nconstraints[i] = 0;
// set default limit duration to 60 timesteps
limit_duration[i] = 60;
reaction_count[i] = 0;
@ -440,9 +445,11 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
rrhandom = new class RanMars*[narrhenius];
int tmp = 0;
for (int i = 0; i < nconstraints; i++) {
if (constraints[i][1] == ARRHENIUS) {
rrhandom[tmp++] = new RanMars(lmp,(int) constraints[i][6] + me);
for (int i = 0; i < nreacts; i++) {
for (int j = 0; j < nconstraints[i]; j++) {
if (constraints[j][i].type == ARRHENIUS) {
rrhandom[tmp++] = new RanMars(lmp,(int) constraints[j][i].par[4] + me);
}
}
}
@ -562,6 +569,9 @@ FixBondReact::~FixBondReact()
memory->destroy(stabilize_steps_flag);
memory->destroy(custom_charges_fragid);
memory->destroy(molecule_keyword);
memory->destroy(constraints);
memory->destroy(nconstraints);
// need to delete rxn_name and constraintstr
memory->destroy(iatomtype);
memory->destroy(jatomtype);
@ -1795,163 +1805,178 @@ int FixBondReact::check_constraints()
tagint atom1,atom2;
double **x = atom->x;
for (int i = 0; i < nconstraints; i++) {
if (constraints[i][0] == rxnID) {
if (constraints[i][1] == DISTANCE) {
get_IDcoords((int) constraints[i][2], (int) constraints[i][3], x1);
get_IDcoords((int) constraints[i][4], (int) constraints[i][5], x2);
delx = x1[0] - x2[0];
dely = x1[1] - x2[1];
delz = x1[2] - x2[2];
domain->minimum_image(delx,dely,delz); // ghost location fix
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < constraints[i][6] || rsq > constraints[i][7]) return 0;
} else if (constraints[i][1] == ANGLE) {
get_IDcoords((int) constraints[i][2], (int) constraints[i][3], x1);
get_IDcoords((int) constraints[i][4], (int) constraints[i][5], x2);
get_IDcoords((int) constraints[i][6], (int) constraints[i][7], x3);
int *satisfied;
memory->create(satisfied,nconstraints[rxnID],"bond/react:satisfied");
for (int i = 0; i < nconstraints[rxnID]; i++)
satisfied[i] = 1;
// 1st bond
delx1 = x1[0] - x2[0];
dely1 = x1[1] - x2[1];
delz1 = x1[2] - x2[2];
domain->minimum_image(delx1,dely1,delz1); // ghost location fix
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
r1 = sqrt(rsq1);
for (int i = 0; i < nconstraints[rxnID]; i++) {
if (constraints[i][rxnID].type == DISTANCE) {
get_IDcoords(constraints[i][rxnID].idtype[0], constraints[i][rxnID].id[0], x1);
get_IDcoords(constraints[i][rxnID].idtype[1], constraints[i][rxnID].id[1], x2);
delx = x1[0] - x2[0];
dely = x1[1] - x2[1];
delz = x1[2] - x2[2];
domain->minimum_image(delx,dely,delz); // ghost location fix
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < constraints[i][rxnID].par[0] || rsq > constraints[i][rxnID].par[1]) satisfied[i] = 0;
} else if (constraints[i][rxnID].type == ANGLE) {
get_IDcoords(constraints[i][rxnID].idtype[0], constraints[i][rxnID].id[0], x1);
get_IDcoords(constraints[i][rxnID].idtype[1], constraints[i][rxnID].id[1], x2);
get_IDcoords(constraints[i][rxnID].idtype[2], constraints[i][rxnID].id[2], x3);
// 2nd bond
delx2 = x3[0] - x2[0];
dely2 = x3[1] - x2[1];
delz2 = x3[2] - x2[2];
domain->minimum_image(delx2,dely2,delz2); // ghost location fix
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
r2 = sqrt(rsq2);
// 1st bond
delx1 = x1[0] - x2[0];
dely1 = x1[1] - x2[1];
delz1 = x1[2] - x2[2];
domain->minimum_image(delx1,dely1,delz1); // ghost location fix
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
r1 = sqrt(rsq1);
// angle (cos and sin)
c = delx1*delx2 + dely1*dely2 + delz1*delz2;
c /= r1*r2;
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
if (acos(c) < constraints[i][8] || acos(c) > constraints[i][9]) return 0;
} else if (constraints[i][1] == DIHEDRAL) {
// phi calculation from dihedral style harmonic
get_IDcoords((int) constraints[i][2], (int) constraints[i][3], x1);
get_IDcoords((int) constraints[i][4], (int) constraints[i][5], x2);
get_IDcoords((int) constraints[i][6], (int) constraints[i][7], x3);
get_IDcoords((int) constraints[i][8], (int) constraints[i][9], x4);
// 2nd bond
delx2 = x3[0] - x2[0];
dely2 = x3[1] - x2[1];
delz2 = x3[2] - x2[2];
domain->minimum_image(delx2,dely2,delz2); // ghost location fix
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
r2 = sqrt(rsq2);
vb1x = x1[0] - x2[0];
vb1y = x1[1] - x2[1];
vb1z = x1[2] - x2[2];
domain->minimum_image(vb1x,vb1y,vb1z);
// angle (cos and sin)
c = delx1*delx2 + dely1*dely2 + delz1*delz2;
c /= r1*r2;
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
if (acos(c) < constraints[i][rxnID].par[0] || acos(c) > constraints[i][rxnID].par[1]) satisfied[i] = 0;
} else if (constraints[i][rxnID].type == DIHEDRAL) {
// phi calculation from dihedral style harmonic
get_IDcoords(constraints[i][rxnID].idtype[0], constraints[i][rxnID].id[0], x1);
get_IDcoords(constraints[i][rxnID].idtype[1], constraints[i][rxnID].id[1], x2);
get_IDcoords(constraints[i][rxnID].idtype[2], constraints[i][rxnID].id[2], x3);
get_IDcoords(constraints[i][rxnID].idtype[3], constraints[i][rxnID].id[3], x4);
vb2x = x3[0] - x2[0];
vb2y = x3[1] - x2[1];
vb2z = x3[2] - x2[2];
domain->minimum_image(vb2x,vb2y,vb2z);
vb1x = x1[0] - x2[0];
vb1y = x1[1] - x2[1];
vb1z = x1[2] - x2[2];
domain->minimum_image(vb1x,vb1y,vb1z);
vb2xm = -vb2x;
vb2ym = -vb2y;
vb2zm = -vb2z;
domain->minimum_image(vb2xm,vb2ym,vb2zm);
vb2x = x3[0] - x2[0];
vb2y = x3[1] - x2[1];
vb2z = x3[2] - x2[2];
domain->minimum_image(vb2x,vb2y,vb2z);
vb3x = x4[0] - x3[0];
vb3y = x4[1] - x3[1];
vb3z = x4[2] - x3[2];
domain->minimum_image(vb3x,vb3y,vb3z);
vb2xm = -vb2x;
vb2ym = -vb2y;
vb2zm = -vb2z;
domain->minimum_image(vb2xm,vb2ym,vb2zm);
ax = vb1y*vb2zm - vb1z*vb2ym;
ay = vb1z*vb2xm - vb1x*vb2zm;
az = vb1x*vb2ym - vb1y*vb2xm;
bx = vb3y*vb2zm - vb3z*vb2ym;
by = vb3z*vb2xm - vb3x*vb2zm;
bz = vb3x*vb2ym - vb3y*vb2xm;
vb3x = x4[0] - x3[0];
vb3y = x4[1] - x3[1];
vb3z = x4[2] - x3[2];
domain->minimum_image(vb3x,vb3y,vb3z);
rasq = ax*ax + ay*ay + az*az;
rbsq = bx*bx + by*by + bz*bz;
rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm;
rg = sqrt(rgsq);
ax = vb1y*vb2zm - vb1z*vb2ym;
ay = vb1z*vb2xm - vb1x*vb2zm;
az = vb1x*vb2ym - vb1y*vb2xm;
bx = vb3y*vb2zm - vb3z*vb2ym;
by = vb3z*vb2xm - vb3x*vb2zm;
bz = vb3x*vb2ym - vb3y*vb2xm;
ra2inv = rb2inv = 0.0;
if (rasq > 0) ra2inv = 1.0/rasq;
if (rbsq > 0) rb2inv = 1.0/rbsq;
rabinv = sqrt(ra2inv*rb2inv);
rasq = ax*ax + ay*ay + az*az;
rbsq = bx*bx + by*by + bz*bz;
rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm;
rg = sqrt(rgsq);
c = (ax*bx + ay*by + az*bz)*rabinv;
s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z);
ra2inv = rb2inv = 0.0;
if (rasq > 0) ra2inv = 1.0/rasq;
if (rbsq > 0) rb2inv = 1.0/rbsq;
rabinv = sqrt(ra2inv*rb2inv);
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
phi = atan2(s,c);
c = (ax*bx + ay*by + az*bz)*rabinv;
s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z);
ANDgate = 0;
if (constraints[i][10] < constraints[i][11]) {
if (phi > constraints[i][10] && phi < constraints[i][11]) ANDgate = 1;
} else {
if (phi > constraints[i][10] || phi < constraints[i][11]) ANDgate = 1;
}
if (constraints[i][12] < constraints[i][13]) {
if (phi > constraints[i][12] && phi < constraints[i][13]) ANDgate = 1;
} else {
if (phi > constraints[i][12] || phi < constraints[i][13]) ANDgate = 1;
}
if (ANDgate != 1) return 0;
} else if (constraints[i][1] == ARRHENIUS) {
t = get_temperature();
prrhob = constraints[i][3]*pow(t,constraints[i][4])*
exp(-constraints[i][5]/(force->boltz*t));
if (prrhob < rrhandom[(int) constraints[i][2]]->uniform()) return 0;
} else if (constraints[i][1] == RMSD) {
// call superpose
int iatom;
int iref = -1; // choose first atom as reference
int n2superpose = 0;
double **xfrozen; // coordinates for the "frozen" target molecule
double **xmobile; // coordinates for the "mobile" molecule
int ifragment = constraints[i][3];
if (ifragment >= 0) {
for (int j = 0; j < onemol->natoms; j++)
if (onemol->fragmentmask[ifragment][j]) n2superpose++;
memory->create(xfrozen,n2superpose,3,"bond/react:xfrozen");
memory->create(xmobile,n2superpose,3,"bond/react:xmobile");
int myincr = 0;
for (int j = 0; j < onemol->natoms; j++) {
if (onemol->fragmentmask[ifragment][j]) {
iatom = atom->map(glove[j][1]);
if (iref == -1) iref = iatom;
iatom = domain->closest_image(iref,iatom);
for (int k = 0; k < 3; k++) {
xfrozen[myincr][k] = x[iatom][k];
xmobile[myincr][k] = onemol->x[j][k];
}
myincr++;
}
}
} else {
int iatom;
int iref = -1; // choose first atom as reference
n2superpose = onemol->natoms;
memory->create(xfrozen,n2superpose,3,"bond/react:xfrozen");
memory->create(xmobile,n2superpose,3,"bond/react:xmobile");
for (int j = 0; j < n2superpose; j++) {
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
phi = atan2(s,c);
ANDgate = 0;
if (constraints[i][rxnID].par[0] < constraints[i][rxnID].par[1]) {
if (phi > constraints[i][rxnID].par[0] && phi < constraints[i][rxnID].par[1]) ANDgate = 1;
} else {
if (phi > constraints[i][rxnID].par[0] || phi < constraints[i][rxnID].par[1]) ANDgate = 1;
}
if (constraints[i][rxnID].par[2] < constraints[i][rxnID].par[3]) {
if (phi > constraints[i][rxnID].par[2] && phi < constraints[i][rxnID].par[3]) ANDgate = 1;
} else {
if (phi > constraints[i][rxnID].par[2] || phi < constraints[i][rxnID].par[3]) ANDgate = 1;
}
if (ANDgate != 1) satisfied[i] = 0;
} else if (constraints[i][rxnID].type == ARRHENIUS) {
t = get_temperature();
prrhob = constraints[i][rxnID].par[1]*pow(t,constraints[i][rxnID].par[2])*
exp(-constraints[i][rxnID].par[3]/(force->boltz*t));
if (prrhob < rrhandom[(int) constraints[i][rxnID].par[0]]->uniform()) satisfied[i] = 0;
} else if (constraints[i][rxnID].type == RMSD) {
// call superpose
int iatom;
int iref = -1; // choose first atom as reference
int n2superpose = 0;
double **xfrozen; // coordinates for the "frozen" target molecule
double **xmobile; // coordinates for the "mobile" molecule
int ifragment = constraints[i][rxnID].id[0];
if (ifragment >= 0) {
for (int j = 0; j < onemol->natoms; j++)
if (onemol->fragmentmask[ifragment][j]) n2superpose++;
memory->create(xfrozen,n2superpose,3,"bond/react:xfrozen");
memory->create(xmobile,n2superpose,3,"bond/react:xmobile");
int myincr = 0;
for (int j = 0; j < onemol->natoms; j++) {
if (onemol->fragmentmask[ifragment][j]) {
iatom = atom->map(glove[j][1]);
if (iref == -1) iref = iatom;
iatom = domain->closest_image(iref,iatom);
for (int k = 0; k < 3; k++) {
xfrozen[j][k] = x[iatom][k];
xmobile[j][k] = onemol->x[j][k];
xfrozen[myincr][k] = x[iatom][k];
xmobile[myincr][k] = onemol->x[j][k];
}
myincr++;
}
}
} else {
int iatom;
int iref = -1; // choose first atom as reference
n2superpose = onemol->natoms;
memory->create(xfrozen,n2superpose,3,"bond/react:xfrozen");
memory->create(xmobile,n2superpose,3,"bond/react:xmobile");
for (int j = 0; j < n2superpose; j++) {
iatom = atom->map(glove[j][1]);
if (iref == -1) iref = iatom;
iatom = domain->closest_image(iref,iatom);
for (int k = 0; k < 3; k++) {
xfrozen[j][k] = x[iatom][k];
xmobile[j][k] = onemol->x[j][k];
}
}
Superpose3D<double, double **> superposer(n2superpose);
double rmsd = superposer.Superpose(xfrozen, xmobile);
if (rmsd > constraints[i][2]) return 0;
memory->destroy(xfrozen);
memory->destroy(xmobile);
}
Superpose3D<double, double **> superposer(n2superpose);
double rmsd = superposer.Superpose(xfrozen, xmobile);
memory->destroy(xfrozen);
memory->destroy(xmobile);
if (rmsd > constraints[i][rxnID].par[0]) satisfied[i] = 0;
}
}
if (nconstraints[rxnID] > 0) {
char evalstr[MAXLINE],*ptr,valstr;
strcpy(evalstr,constraintstr[rxnID]);
for (int i = 0; i < nconstraints[rxnID]; i++) {
sprintf(&valstr,"%d", satisfied[i]);
ptr = strchr(evalstr,'C');
*ptr = valstr;
}
double verdict = input->variable->evaluate_boolean(evalstr);
if (verdict == 0.0) return 0;
}
// let's also check chirality within 'check_constraint'
for (int i = 0; i < onemol->natoms; i++) {
if (chiral_atoms[i][0][rxnID] == 1) {
@ -1975,6 +2000,7 @@ int FixBondReact::check_constraints()
}
}
memory->destroy(satisfied);
return 1;
}
@ -1987,7 +2013,7 @@ fragment: given pre-reacted molID (onemol) and fragID,
void FixBondReact::get_IDcoords(int mode, int myID, double *center)
{
double **x = atom->x;
if (mode == 1) {
if (mode == ATOM) {
int iatom = atom->map(glove[myID-1][1]);
for (int i = 0; i < 3; i++)
center[i] = x[iatom][i];
@ -3216,8 +3242,8 @@ void FixBondReact::read(int myrxn)
else if (strstr(line,"deleteIDs")) sscanf(line,"%d",&ndelete);
else if (strstr(line,"chiralIDs")) sscanf(line,"%d",&nchiral);
else if (strstr(line,"constraints")) {
sscanf(line,"%d",&nconstr);
memory->grow(constraints,nconstraints+nconstr,MAXCONARGS,"bond/react:constraints");
sscanf(line,"%d",&nconstraints[myrxn]);
memory->grow(constraints,nconstraints[myrxn],nreacts,"bond/react:constraints");
} else break;
}
@ -3252,7 +3278,7 @@ void FixBondReact::read(int myrxn)
} else if (strcmp(keyword,"ChiralIDs") == 0) {
ChiralCenters(line, myrxn);
} else if (strcmp(keyword,"Constraints") == 0) {
Constraints(line, myrxn);
ReadConstraints(line, myrxn);
} else error->one(FLERR,"Bond/react: Unknown section in map file");
parse_keyword(1,line,keyword);
@ -3350,69 +3376,89 @@ void FixBondReact::ChiralCenters(char *line, int myrxn)
}
}
void FixBondReact::Constraints(char *line, int myrxn)
void FixBondReact::ReadConstraints(char *line, int myrxn)
{
double tmp[MAXCONARGS];
char **strargs;
char **strargs,*ptr;
memory->create(strargs,MAXCONARGS,MAXLINE,"bond/react:strargs");
char *constraint_type = new char[MAXLINE];
for (int i = 0; i < nconstr; i++) {
strcpy(constraintstr[myrxn],"("); // string for boolean constraint logic
for (int i = 0; i < nconstraints[myrxn]; i++) {
readline(line);
if (ptr = strrchr(line,'(')) { // reverse char search
strncat(constraintstr[myrxn],line,ptr-line+1);
line = ptr + 1;
}
// 'C' indicates where to sub in next constraint
strcat(constraintstr[myrxn],"C");
if (ptr = strchr(line,')')) {
strncat(constraintstr[myrxn],ptr,strrchr(line,')')-ptr+1);
}
if (ptr = strstr(line,"&&")) {
strcat(constraintstr[myrxn],"&&");
*ptr = '\0';
} else if (ptr = strstr(line,"||")) {
strcat(constraintstr[myrxn],"||");
*ptr = '\0';
} else if (i+1 < nconstraints[myrxn]){
strcat(constraintstr[myrxn],"&&");
}
if (ptr = strchr(line,')'))
*ptr = '\0';
sscanf(line,"%s",constraint_type);
constraints[nconstraints][0] = myrxn;
if (strcmp(constraint_type,"distance") == 0) {
constraints[nconstraints][1] = DISTANCE;
constraints[i][myrxn].type = DISTANCE;
sscanf(line,"%*s %s %s %lg %lg",strargs[0],strargs[1],&tmp[0],&tmp[1]);
readID(strargs[0], nconstraints, 2, 3);
readID(strargs[1], nconstraints, 4, 5);
readID(strargs[0], i, myrxn, 0);
readID(strargs[1], i, myrxn, 1);
// cutoffs
constraints[nconstraints][6] = tmp[0]*tmp[0]; // using square of distance
constraints[nconstraints][7] = tmp[1]*tmp[1];
constraints[i][myrxn].par[0] = tmp[0]*tmp[0]; // using square of distance
constraints[i][myrxn].par[1] = tmp[1]*tmp[1];
} else if (strcmp(constraint_type,"angle") == 0) {
constraints[nconstraints][1] = ANGLE;
constraints[i][myrxn].type = ANGLE;
sscanf(line,"%*s %s %s %s %lg %lg",strargs[0],strargs[1],strargs[2],&tmp[0],&tmp[1]);
readID(strargs[0], nconstraints, 2, 3);
readID(strargs[1], nconstraints, 4, 5);
readID(strargs[2], nconstraints, 6, 7);
constraints[nconstraints][8] = tmp[0]/180.0 * MY_PI;
constraints[nconstraints][9] = tmp[1]/180.0 * MY_PI;
readID(strargs[0], i, myrxn, 0);
readID(strargs[1], i, myrxn, 1);
readID(strargs[2], i, myrxn, 2);
constraints[i][myrxn].par[0] = tmp[0]/180.0 * MY_PI;
constraints[i][myrxn].par[1] = tmp[1]/180.0 * MY_PI;
} else if (strcmp(constraint_type,"dihedral") == 0) {
constraints[nconstraints][1] = DIHEDRAL;
constraints[i][myrxn].type = DIHEDRAL;
tmp[2] = 181.0; // impossible range
tmp[3] = 182.0;
sscanf(line,"%*s %s %s %s %s %lg %lg %lg %lg",strargs[0],strargs[1],
strargs[2],strargs[3],&tmp[0],&tmp[1],&tmp[2],&tmp[3]);
readID(strargs[0], nconstraints, 2, 3);
readID(strargs[1], nconstraints, 4, 5);
readID(strargs[2], nconstraints, 6, 7);
readID(strargs[3], nconstraints, 8, 9);
constraints[nconstraints][10] = tmp[0]/180.0 * MY_PI;
constraints[nconstraints][11] = tmp[1]/180.0 * MY_PI;
constraints[nconstraints][12] = tmp[2]/180.0 * MY_PI;
constraints[nconstraints][13] = tmp[3]/180.0 * MY_PI;
readID(strargs[0], i, myrxn, 0);
readID(strargs[1], i, myrxn, 1);
readID(strargs[2], i, myrxn, 2);
readID(strargs[3], i, myrxn, 3);
constraints[i][myrxn].par[0] = tmp[0]/180.0 * MY_PI;
constraints[i][myrxn].par[1] = tmp[1]/180.0 * MY_PI;
constraints[i][myrxn].par[2] = tmp[2]/180.0 * MY_PI;
constraints[i][myrxn].par[3] = tmp[3]/180.0 * MY_PI;
} else if (strcmp(constraint_type,"arrhenius") == 0) {
constraints[nconstraints][1] = ARRHENIUS;
constraints[nconstraints][2] = narrhenius++;
constraints[i][myrxn].type = ARRHENIUS;
constraints[i][myrxn].par[0] = narrhenius++;
sscanf(line,"%*s %lg %lg %lg %lg",&tmp[0],&tmp[1],&tmp[2],&tmp[3]);
constraints[nconstraints][3] = tmp[0];
constraints[nconstraints][4] = tmp[1];
constraints[nconstraints][5] = tmp[2];
constraints[nconstraints][6] = tmp[3];
constraints[i][myrxn].par[1] = tmp[0];
constraints[i][myrxn].par[2] = tmp[1];
constraints[i][myrxn].par[3] = tmp[2];
constraints[i][myrxn].par[4] = tmp[3];
} else if (strcmp(constraint_type,"rmsd") == 0) {
constraints[nconstraints][1] = RMSD;
constraints[i][myrxn].type = RMSD;
strcpy(strargs[0],"0");
sscanf(line,"%*s %lg %s",&tmp[0],strargs[0]);
constraints[nconstraints][2] = tmp[0]; // RMSDmax
constraints[nconstraints][3] = -1; // optional molecule fragment
constraints[i][myrxn].par[0] = tmp[0]; // RMSDmax
constraints[i][myrxn].id[0] = -1; // optional molecule fragment
if (isalpha(strargs[0][0])) {
int ifragment = onemol->findfragment(strargs[0]);
if (ifragment < 0) error->one(FLERR,"Bond/react: Molecule fragment does not exist");
else constraints[nconstraints][3] = ifragment;
else constraints[i][myrxn].id[0] = ifragment;
}
} else
error->one(FLERR,"Bond/react: Illegal constraint type in 'Constraints' section of map file");
nconstraints++;
}
strcat(constraintstr[myrxn],")"); // close boolean constraint logic string
delete [] constraint_type;
memory->destroy(strargs);
}
@ -3422,18 +3468,18 @@ if ID starts with character, assume it is a pre-reaction molecule fragment ID
otherwise, it is a pre-reaction atom ID
---------------------------------------------------------------------- */
void FixBondReact::readID(char *strarg, int iconstr, int mode, int myID)
void FixBondReact::readID(char *strarg, int iconstr, int myrxn, int i)
{
if (isalpha(strarg[0])) {
constraints[iconstr][mode] = 0; // fragment vs. atom ID flag
constraints[iconstr][myrxn].idtype[i] = FRAG; // fragment vs. atom ID flag
int ifragment = onemol->findfragment(strarg);
if (ifragment < 0) error->one(FLERR,"Bond/react: Molecule fragment does not exist");
constraints[iconstr][myID] = ifragment;
constraints[iconstr][myrxn].id[i] = ifragment;
} else {
constraints[iconstr][mode] = 1; // fragment vs. atom ID flag
constraints[iconstr][myrxn].idtype[i] = ATOM; // fragment vs. atom ID flag
int iatom = atoi(strarg);
if (iatom > onemol->natoms) error->one(FLERR,"Bond/react: Invalid template atom ID in map file");
constraints[iconstr][myID] = iatom;
constraints[iconstr][myrxn].id[i] = iatom;
}
}

20
src/USER-REACTION/fix_bond_react.h
View File

@ -31,7 +31,9 @@ namespace LAMMPS_NS {
class FixBondReact : public Fix {
public:
enum {MAXLINE=256};
enum {MAXLINE=256}; // max length of line read from files
enum {MAXCONIDS=4}; // max # of IDs used by any constraint
enum {MAXCONPAR=4}; // max # of constraint parameters
FixBondReact(class LAMMPS *, int, char **);
~FixBondReact();
@ -66,9 +68,9 @@ class FixBondReact : public Fix {
int *stabilize_steps_flag;
int *custom_charges_fragid;
int *molecule_keyword;
int nconstraints;
int *nconstraints;
char **constraintstr;
int narrhenius;
double **constraints;
int **var_flag,**var_id; // for keyword values with variable inputs
int status;
int *groupbits;
@ -114,7 +116,7 @@ class FixBondReact : public Fix {
int *ibonding,*jbonding;
int *closeneigh; // indicates if bonding atoms of a rxn are 1-2, 1-3, or 1-4 neighbors
int nedge,nequivalent,ndelete,nchiral,nconstr; // # edge, equivalent atoms in mapping file
int nedge,nequivalent,ndelete,nchiral; // # edge, equivalent atoms in mapping file
int attempted_rxn; // there was an attempt!
int *local_rxn_count;
int *ghostly_rxn_count;
@ -155,7 +157,7 @@ class FixBondReact : public Fix {
void DeleteAtoms(char *, int);
void CustomCharges(int, int);
void ChiralCenters(char *, int);
void Constraints(char *, int);
void ReadConstraints(char *, int);
void readID(char *, int, int, int);
void make_a_guess ();
@ -195,6 +197,14 @@ class FixBondReact : public Fix {
};
Set *set;
struct Constraint {
int type;
int id[MAXCONIDS];
int idtype[MAXCONIDS];
double par[MAXCONPAR];
};
Constraint **constraints;
// DEBUG
void print_bb();