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

add option to enforce atom chirality

This commit is contained in:
jrgissing 2019-12-02 12:27:57 -05:00
parent 36e102516f
commit 819fe9ec56
2 changed files with 109 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -36,6 +36,7 @@ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu)
#include "group.h"
#include "citeme.h"
#include "math_const.h"
#include "math_extra.h"
#include "memory.h"
#include "error.h"
@ -297,6 +298,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
memory->create(landlocked_atoms,max_natoms,nreacts,"bond/react:landlocked_atoms");
memory->create(custom_edges,max_natoms,nreacts,"bond/react:custom_edges");
memory->create(delete_atoms,max_natoms,nreacts,"bond/react:delete_atoms");
memory->create(chiral_atoms,max_natoms,6,nreacts,"bond/react:chiral_atoms");
for (int j = 0; j < nreacts; j++)
for (int i = 0; i < max_natoms; i++) {
@ -304,6 +306,9 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
if (update_edges_flag[j] == 1) custom_edges[i][j] = 1;
else custom_edges[i][j] = 0;
delete_atoms[i][j] = 0;
for (int k = 0; k < 6; k++) {
chiral_atoms[i][k][j] = 0;
}
}
// read all map files afterward
@ -430,6 +435,7 @@ FixBondReact::~FixBondReact()
memory->destroy(landlocked_atoms);
memory->destroy(custom_edges);
memory->destroy(delete_atoms);
memory->destroy(chiral_atoms);
memory->destroy(nevery);
memory->destroy(cutsq);
@ -1649,10 +1655,11 @@ evaluate constraints: return 0 if any aren't satisfied
int FixBondReact::check_constraints()
{
tagint atom1,atom2,atom3;
double delx,dely,delz,rsq;
double unwrap[3],delx,dely,delz,rsq;
double delx1,dely1,delz1,delx2,dely2,delz2;
double rsq1,rsq2,r1,r2,c,t,prrhob;
imageint *image = atom->image;
double **x = atom->x;
for (int i = 0; i < nconstraints; i++) {
@ -1701,6 +1708,30 @@ int FixBondReact::check_constraints()
}
}
}
// let's also check chirality within 'check_constraint'
for (int i = 0; i < onemol->natoms; i++) {
if (chiral_atoms[i][0][rxnID] == 1) {
double my4coords[12];
// already ensured, by transitive property, that chiral simulation atom has four neighs
for (int j = 0; j < 4; j++) {
atom1 = atom->map(glove[i][1]);
// loop over known types involved in chiral center
for (int jj = 0; jj < 4; jj++) {
if (atom->type[atom->map(xspecial[atom1][j])] == chiral_atoms[i][jj+2][rxnID]) {
atom2 = atom->map(xspecial[atom1][j]);
domain->unmap(x[atom2],image[atom2],unwrap);
for (int k = 0; k < 3; k++) {
my4coords[3*jj+k] = unwrap[k];
}
break;
}
}
}
if (get_chirality(my4coords) != chiral_atoms[i][1][rxnID]) return 0;
}
}
return 1;
}
@ -1741,6 +1772,33 @@ double FixBondReact::get_temperature()
return t;
}
/* ----------------------------------------------------------------------
return handedness (1 or -1) of a chiral center, given ordered set of coordinates
------------------------------------------------------------------------- */
int FixBondReact::get_chirality(double four_coords[12])
{
// define oriented plane with first three coordinates
double vec1[3],vec2[3],vec3[3],vec4[3],mean3[3],dot;
for (int i = 0; i < 3; i++) {
vec1[i] = four_coords[i]-four_coords[i+3];
vec2[i] = four_coords[i+3]-four_coords[i+6];
}
MathExtra::cross3(vec1,vec2,vec3);
for (int i = 0; i < 3; i++) {
mean3[i] = (four_coords[i] + four_coords[i+3] +
four_coords[i+6])/3;
vec4[i] = four_coords[i+9] - mean3[i];
}
dot = MathExtra::dot3(vec3,vec4);
dot = dot/fabs(dot);
return (int) dot;
}
/* ----------------------------------------------------------------------
Get xspecials for current molecule templates
------------------------------------------------------------------------- */
@ -2869,6 +2927,7 @@ void FixBondReact::read(int myrxn)
}
else if (strstr(line,"customIDs")) sscanf(line,"%d",&ncustom);
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");
@ -2900,6 +2959,8 @@ void FixBondReact::read(int myrxn)
CustomEdges(line, myrxn);
} else if (strcmp(keyword,"DeleteIDs") == 0) {
DeleteAtoms(line, myrxn);
} else if (strcmp(keyword,"ChiralIDs") == 0) {
ChiralCenters(line, myrxn);
} else if (strcmp(keyword,"Constraints") == 0) {
Constraints(line, myrxn);
} else error->one(FLERR,"Bond/react: Unknown section in map file");
@ -2977,6 +3038,37 @@ void FixBondReact::DeleteAtoms(char *line, int myrxn)
}
}
void FixBondReact::ChiralCenters(char *line, int myrxn)
{
int tmp;
for (int i = 0; i < nchiral; i++) {
readline(line);
sscanf(line,"%d",&tmp);
chiral_atoms[tmp-1][0][myrxn] = 1;
if (onemol->xflag == 0)
error->one(FLERR,"Bond/react: Molecule template 'Coords' section required for chiralIDs keyword");
if ((int) onemol_nxspecial[tmp-1][0] != 4)
error->one(FLERR,"Bond/react: Chiral atoms must have exactly four first neighbors");
for (int j = 0; j < 4; j++) {
for (int k = j+1; k < 4; k++) {
if (onemol->type[onemol_xspecial[tmp-1][j]-1] ==
onemol->type[onemol_xspecial[tmp-1][k]-1])
error->one(FLERR,"Bond/react: First neighbors of chiral atoms must be of mutually different types");
}
}
// record order of atom types, and coords
double my4coords[12];
for (int j = 0; j < 4; j++) {
chiral_atoms[tmp-1][j+2][myrxn] = onemol->type[onemol_xspecial[tmp-1][j]-1];
for (int k = 0; k < 3; k++) {
my4coords[3*j+k] = onemol->x[onemol_xspecial[tmp-1][j]-1][k];
}
}
// get orientation
chiral_atoms[tmp-1][1][myrxn] = get_chirality(my4coords);
}
}
void FixBondReact::Constraints(char *line, int myrxn)
{
double tmp[MAXCONARGS];

17
src/USER-MISC/fix_bond_react.h
View File

@ -110,7 +110,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,ncustom,ndelete,nconstr; // # edge, equivalent, custom atoms in mapping file
int nedge,nequivalent,ncustom,ndelete,nchiral,nconstr; // # edge, equivalent, custom atoms in mapping file
int attempted_rxn; // there was an attempt!
int *local_rxn_count;
int *ghostly_rxn_count;
@ -126,6 +126,7 @@ class FixBondReact : public Fix {
int **landlocked_atoms; // all atoms at least three bonds away from edge atoms
int **custom_edges; // atoms in molecule templates with incorrect valences
int **delete_atoms; // atoms in pre-reacted templates to delete
int ***chiral_atoms; // pre-react chiral atoms. 1) flag 2) orientation 3-4) ordered atom types
int **nxspecial,**onemol_nxspecial,**twomol_nxspecial; // full number of 1-4 neighbors
tagint **xspecial,**onemol_xspecial,**twomol_xspecial; // full 1-4 neighbor list
@ -149,6 +150,7 @@ class FixBondReact : public Fix {
void Equivalences(char *,int);
void CustomEdges(char *,int);
void DeleteAtoms(char *,int);
void ChiralCenters(char *,int);
void Constraints(char *,int);
void make_a_guess ();
@ -159,6 +161,7 @@ class FixBondReact : public Fix {
void ring_check();
int check_constraints();
double get_temperature();
int get_chirality(double[12]); // get handedness given an ordered set of coordinates
void open(char *);
void readline(char *);
@ -249,6 +252,18 @@ E: Bond/react: A deleted atom cannot remain bonded to an atom that is not delete
Self-explanatory.
E: Bond/react: First neighbors of chiral atoms must be of mutually different types
Self-explanatory.
E: Bond/react: Chiral atoms must have exactly four first neighbors
Self-explanatory.
E: Bond/react: Molecule template 'Coords' section required for chiralIDs keyword
The coordinates of atoms in the pre-reacted template are used to determine chirality.
E: Bond/react special bond generation overflow
The number of special bonds per-atom created by a reaction exceeds the