jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #1067 from jrgissing/bond/react-better-template-checks 

Bond/react: better template validity checks
This commit is contained in:
Steve Plimpton 2018-08-22 08:36:44 -06:00 committed by GitHub
parent 25bc74eb21 2baa0b3452
commit 316f265e3d
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
2 changed files with 41 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
doc/src/fix_bond_react.txt
View File

@ -125,13 +125,13 @@ template. If both these conditions are met, the reaction site is
modified to match the post-reaction template.
A bonding atom pair will be identified if several conditions are met.
First, a pair of atoms within the specified react-group-ID of type
typei and typej must separated by a distance between {Rmin} and
First, a pair of atoms I,J within the specified react-group-ID of type
itype and jtype must separated by a distance between {Rmin} and
{Rmax}. It is possible that multiple bonding atom pairs are
identified: if the bonding atoms in the pre-reacted template are not
1-2, 1-3, or 1-4 neighbors, the closest bonding atom partner is set as
its bonding partner; otherwise, the farthest potential partner is
chosen. Then, if both an atomi and atomj have each other as their
chosen. Then, if both an atom I and atom J have each other as their
nearest bonding partners, these two atoms are identified as the
bonding atom pair of the reaction site. Once this unique bonding atom
pair is identified for each reaction, there could two or more
@ -323,7 +323,7 @@ bond/break"_fix_bond_break.html, "fix bond/swap"_fix_bond_swap.html,
[Default:]
The option defaults are stabilization = no, stabilize_steps = 60
The option defaults are stabilization = no, prob = 1.0, stabilize_steps = 60
:line

40
src/USER-MISC/fix_bond_react.cpp
View File

@ -1063,8 +1063,8 @@ void FixBondReact::superimpose_algorithm()
hang_catch++;
// let's go ahead and catch the simplest of hangs
//if (hang_catch > onemol->natoms*4)
if (hang_catch > atom->nlocal*3) {
error->all(FLERR,"Excessive iteration of superimpose algorithm");
if (hang_catch > atom->nlocal*30) {
error->one(FLERR,"Excessive iteration of superimpose algorithm");
}
}
}
@ -1554,6 +1554,40 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
error->one(FLERR,"Atom affected by reaction too close to template edge");
}
// additionally, if a bond changes type, but neither involved atom is landlocked, bad
// would someone want to change an angle type but not bond or atom types? (etc.) ...hopefully not yet
for (int i = 0; i < twomol->natoms; i++) {
if (landlocked_atoms[i][myrxn] == 0) {
int twomol_atomi = i+1;
for (int j = 0; j < twomol->num_bond[i]; j++) {
int twomol_atomj = twomol->bond_atom[i][j];
if (landlocked_atoms[twomol_atomj-1][myrxn] == 0) {
int onemol_atomi = equivalences[i][1][myrxn];
int onemol_batom;
for (int m = 0; m < onemol->num_bond[onemol_atomi-1]; m++) {
onemol_batom = onemol->bond_atom[onemol_atomi-1][m];
if (onemol_batom == equivalences[twomol_atomj-1][1][myrxn]) {
if (twomol->bond_type[i][j] != onemol->bond_type[onemol_atomi-1][m]) {
error->one(FLERR,"Bond type affected by reaction too close to template edge");
}
}
}
if (newton_bond) {
int onemol_atomj = equivalences[twomol_atomj-1][1][myrxn];
for (int m = 0; m < onemol->num_bond[onemol_atomj-1]; m++) {
onemol_batom = onemol->bond_atom[onemol_atomj-1][m];
if (onemol_batom == equivalences[i][1][myrxn]) {
if (twomol->bond_type[i][j] != onemol->bond_type[onemol_atomj-1][m]) {
error->one(FLERR,"Bond type affected by reaction too close to template edge");
}
}
}
}
}
}
}
}
// also, if atoms change number of bonds, but aren't landlocked, that could be bad
if (me == 0)
for (int i = 0; i < twomol->natoms; i++) {
@ -2067,7 +2101,7 @@ void FixBondReact::update_everything()
nspecial[atom->map(update_mega_glove[jj+1][i])][1]++;
nspecial[atom->map(update_mega_glove[jj+1][i])][2]++;
}
for (int n = nspecial[atom->map(update_mega_glove[jj+1][i])][2]; n > insert_num; n--) {
for (int n = nspecial[atom->map(update_mega_glove[jj+1][i])][2]-1; n > insert_num; n--) {
special[atom->map(update_mega_glove[jj+1][i])][n] = special[atom->map(update_mega_glove[jj+1][i])][n-1];
}
special[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->special[j][k]-1][1][rxnID]][i];