From d429143589a3fd14fb42470ceb04db27ef3efb30 Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sat, 19 Oct 2019 22:15:52 -0600 Subject: [PATCH] bond/react: Arrhenius constraint --- doc/src/Eqs/fix_bond_react.jpg | Bin 0 -> 2427 bytes doc/src/Eqs/fix_bond_react.tex | 9 +++++ doc/src/fix_bond_react.txt | 23 ++++++++++- src/USER-MISC/fix_bond_react.cpp | 67 +++++++++++++++++++++++++++++-- src/USER-MISC/fix_bond_react.h | 5 ++- 5 files changed, 99 insertions(+), 5 deletions(-) create mode 100644 doc/src/Eqs/fix_bond_react.jpg create mode 100644 doc/src/Eqs/fix_bond_react.tex diff --git a/doc/src/Eqs/fix_bond_react.jpg b/doc/src/Eqs/fix_bond_react.jpg new file mode 100644 index 0000000000000000000000000000000000000000..d63b983230c0c26a7bf9c50b38593bdd65196a35 GIT binary patch literal 2427 zcmbuBS5(u962|`_^m0KG+A2jwA-KRI3km{WDUlE$h%_lm7ZPbJXy{582&j~>w9q0= zq(lf!su+j}QbGiz2umkYgAz&>-|qYU&C56EJ7?yc?_ticzp!Tk5feB94uC))ApE2N z_9S5NtNXuTcK`ysz;8eh2y_Zyzdum`2l!;e|G>+`$Ir!m;=x2tJOB)UI5{{@#0gy7 zU``$o2QR<}fr^QWpOWO4u*;&eCjj1GBd!B8l(&_- zm=EMhkUClXbJy-6;|rs&CmR+|h0P;O$1+?V-MFNd(KK+SYPZE*t3zD-jK05i%vbG> zOtm6G>nwG2U&r>2XN=mNm+_cDcanZcL#aKC^>-J26rd#TQ}OQB!VB4Y*-f>nBfVd* zs(>eKO*?NWCtV__6&qXQA#y>dXH$kp1-lhrPqTsG4+j3JXgBF!+3mzpyVNZB;OCAO zS9Q$H_g|=-qUKV|LvNNH z=P4i4V99&&CPKuJE7zQMhPUjC%m@#W)46qr5a}BKa`pN5UZ{*~3sk{=caMsC$UMl? zKB`LrsUM4yDhl0LRnKBfo9>J#JPE^H%vmRJQdbdMHag1G6}V+Nn8tT1%^hmd6fxxC)9w-d zgbBmru*o|%Qf>_EY&aVz*R5s)e{0Pgzp`KhjCRp>Z`k*;ermPT|uUmkdCBWTh4ah-hsq*l|PNm z&Fh3>*TP5#^P}13(FmG>%-&3GGZrN@wljYOv-eh*ZH+V|>irjy4Ca-V)s*zJGh*&! z@B}tq@E;*m2`c5EOwuVm#wh7EKxC<7M0hncjHnRSI0T!ev$@paRg zdBPgMN85_h&~a@c{z&}RacfKDegup}m=wQmsp)Wq4b;RQMvmKrwki3U*<;#rl%_3G zlE0Hu@10hk2?rzGK6OuxB9clz)U+Es|67wdl~%_wE;q0c^7V|48C61*gSjzcNpbRV~l+lJ{wNT44x>AhK9dEn3o*`Y3%(UsVJu_l|0{|2V#}=)?OV z<<480q0GdR(eqIyr2I!|5v=79={Do`3pZCb<{l6LVd?xM)b$H8Z(!R?(?0kNG zZzpc4aWf}RYo1uqJs9-Te21pxKy(OEDR>nzdzr_YZZ>1*ihGCN@{dM1-hGfDbs-@e zG~Q(I9a9T$aBkl9Tz>opS#Y@|hDquX7^|d~vT$_GO}|XtXKR6iE=uBp2lBKZ4MWeU zg<=V4pT1cK8JCPN9jA6TGO}9Ea*)l{(l^3lQPVYN7!#q7e0E%wI=57<5*e|LELr3# z#0cd3Vv3l~Q_SQEdKa9#{o=X%!$R!YKVNU~^jozwkgQsp#@#U(K%4Igc;xr4QRmDT zBjf_CApGeHvMpPQ4VWA_HeWOkwBL#N>%5`G-GMs6U51#Lg@8)jM`hKON{5ev)2Z+o zTM!o(iRNxDw9I1)sM2I-amuDigJnlKWugd-V@Z92+p%B2y)pLNf#2^BY7$v0`NumM zspr5INj~Se>H|4LZF+^hZpovRUZ%Q2SA?z)s!PZqLq_%TWg@nT#hL!)UD8}Art}bnnPq{%~_hps=cn8ntAJu&;{&bl<gLkyem zj9d)+)8AeGDTksZRNgY`^N{NYT2#Wf;<7lRO$I--f+B_vEFVi$VxsAZS1)y*2{4j1 zW(BQLs(*$>>6$p#{m43!3jyV zXT$K*Zy!7Zi+darnzS6~yvWJ=mfN4EbJ}o|akPhqtTC~s>Q&8W+oXnNofo76O;jX) z2(;<%>a^bk^WD?(6PHNg{Ex>X>g%us1-#~zSSove>Eu+V9MkHqhTjZ#yFs=eu~%#k zv2~K4H+!~KGG)_ssd4Ho=EpcD?GhS7z9A$c?;mQUZJ+CLD!Rh*CI4>aP4TtkV7!Ig z#BB!~scwEZ2e&&qT(Q{ARgca*Zb$3dg=i^J{nJo>ewcnwqWsOf)wS literal 0 HcmV?d00001 diff --git a/doc/src/Eqs/fix_bond_react.tex b/doc/src/Eqs/fix_bond_react.tex new file mode 100644 index 0000000000..9400656038 --- /dev/null +++ b/doc/src/Eqs/fix_bond_react.tex @@ -0,0 +1,9 @@ +\documentstyle[12pt]{article} +\pagestyle{empty} +\begin{document} + +\begin{eqnarray*} + k = AT^{n}e^{\frac{-E_{a}}{k_{B}T}} +\end{eqnarray*} + +\end{document} diff --git a/doc/src/fix_bond_react.txt b/doc/src/fix_bond_react.txt index fb8ed95afb..d101a20da5 100644 --- a/doc/src/fix_bond_react.txt +++ b/doc/src/fix_bond_react.txt @@ -266,7 +266,7 @@ either 'none' or 'charges.' Further details are provided in the discussion of the 'update_edges' keyword. The fourth optional section begins with the keyword 'Constraints' and lists additional criteria that must be satisfied in order for the reaction to occur. Currently, -there are two types of constraints available, as discussed below. +there are three types of constraints available, as discussed below. A sample map file is given below: @@ -320,6 +320,27 @@ the central atom). Angles must be specified in degrees. This constraint can be used to enforce a certain orientation between reacting molecules. +The constraint of type 'arrhenius' imposes an additional reaction +probability according to the temperature-dependent Arrhenius equation: + +:c,image(Eqs/fix_bond_react.jpg) + +The Arrhenius constraint has the following syntax: + +arrhenius {A} {n} {E_a} {seed} :pre + +where 'arrhenius' is the required keyword, {A} is the pre-exponential +factor, {n} is the exponent of the temperature dependence, {E_a} is +the activation energy ("units"_units.html of energy), and {seed} is a +random number seed. The temperature is defined as the instantaneous +temperature averaged over all atoms in the reaction site, and is +calculated in the same manner as for example +"compute_temp_chunk"_compute_temp_chunk.html. Currently, there are no +options for additional temperature averaging or velocity-biased +temperature calculations. A uniform random number between 0 and 1 is +generated using {seed}; if this number is less than the result of the +Arrhenius equation above, the reaction is permitted occur. + 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, diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp index c34a7494d9..f52227ce26 100644 --- a/src/USER-MISC/fix_bond_react.cpp +++ b/src/USER-MISC/fix_bond_react.cpp @@ -71,7 +71,7 @@ static const char cite_fix_bond_react[] = enum{ACCEPT,REJECT,PROCEED,CONTINUE,GUESSFAIL,RESTORE}; // types of available reaction constraints -enum{DISTANCE,ANGLE}; +enum{DISTANCE,ANGLE,ARRHENIUS}; /* ---------------------------------------------------------------------- */ @@ -100,6 +100,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : extvector = 0; rxnID = 0; nconstraints = 0; + narrhenius = 0; status = PROCEED; nxspecial = NULL; @@ -322,6 +323,16 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : find_landlocked_atoms(i); } + // initialize Marsaglia RNG with processor-unique seed (Arrhenius prob) + + 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++) { delete [] files[i]; } @@ -348,7 +359,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : } } - // initialize Marsaglia RNG with processor-unique seed + // initialize Marsaglia RNG with processor-unique seed ('prob' keyword) random = new class RanMars*[nreacts]; for (int i = 0; i < nreacts; i++) { @@ -1640,7 +1651,7 @@ int FixBondReact::check_constraints() tagint atom1,atom2,atom3; double delx,dely,delz,rsq; double delx1,dely1,delz1,delx2,dely2,delz2; - double rsq1,rsq2,r1,r2,c; + double rsq1,rsq2,r1,r2,c,t,prrhob; double **x = atom->x; @@ -1680,12 +1691,54 @@ int FixBondReact::check_constraints() if (c > 1.0) c = 1.0; if (c < -1.0) c = -1.0; if (acos(c) < constraints[i][5] || acos(c) > constraints[i][6]) 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; } } } return 1; } +/* ---------------------------------------------------------------------- +compute local temperature: average over all atoms in reaction template +------------------------------------------------------------------------- */ + +double FixBondReact::get_temperature() +{ + int i,ilocal; + double adof = domain->dimension; + + double **v = atom->v; + double *mass = atom->mass; + double *rmass = atom->rmass; + int *type = atom->type; + + double t = 0.0; + + if (rmass) { + for (i = 0; i < onemol->natoms; i++) { + ilocal = atom->map(glove[i][1]); + t += (v[ilocal][0]*v[ilocal][0] + v[ilocal][1]*v[ilocal][1] + + v[ilocal][2]*v[ilocal][2]) * rmass[ilocal]; + } + } else { + for (i = 0; i < onemol->natoms; i++) { + ilocal = atom->map(glove[i][1]); + t += (v[ilocal][0]*v[ilocal][0] + v[ilocal][1]*v[ilocal][1] + + v[ilocal][2]*v[ilocal][2]) * mass[type[ilocal]]; + } + } + + // final temperature + double dof = adof*onemol->natoms; + double tfactor = force->mvv2e / (dof * force->boltz); + t *= tfactor; + return t; +} + /* ---------------------------------------------------------------------- Get xspecials for current molecule templates ------------------------------------------------------------------------- */ @@ -2947,6 +3000,14 @@ void FixBondReact::Constraints(char *line, int myrxn) constraints[nconstraints][4] = tmp[2]; constraints[nconstraints][5] = tmp[3]/180.0 * MY_PI; constraints[nconstraints][6] = tmp[4]/180.0 * MY_PI; + } else if (strcmp(constraint_type,"arrhenius") == 0) { + constraints[nconstraints][1] = ARRHENIUS; + constraints[nconstraints][2] = 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]; } else error->one(FLERR,"Bond/react: Illegal constraint type in 'Constraints' section of map file"); nconstraints++; diff --git a/src/USER-MISC/fix_bond_react.h b/src/USER-MISC/fix_bond_react.h index eda26f129d..f59e051ed1 100644 --- a/src/USER-MISC/fix_bond_react.h +++ b/src/USER-MISC/fix_bond_react.h @@ -65,6 +65,7 @@ class FixBondReact : public Fix { int *stabilize_steps_flag; int *update_edges_flag; int nconstraints; + int narrhenius; double **constraints; int status; int *groupbits; @@ -88,7 +89,8 @@ class FixBondReact : public Fix { Fix *fix2; // properties/atom used to indicate 1) relaxing atoms // 2) to which 'react' atom belongs Fix *fix3; // property/atom used for system-wide thermostat - class RanMars **random; + class RanMars **random; // random number for 'prob' keyword + class RanMars **rrhandom; // random number for Arrhenius constraint class NeighList *list; int *reacted_mol,*unreacted_mol; @@ -156,6 +158,7 @@ class FixBondReact : public Fix { void inner_crosscheck_loop(); void ring_check(); int check_constraints(); + double get_temperature(); void open(char *); void readline(char *);