Merge pull request #1734 from jrgissing/arrhenius_constraint (3d7881b3) · Commits · 郑智淋 / lammps

doc/src/Eqs/fix_bond_react.jpg

0 → 100644

+2.37 KiB

Loading image diff...

doc/src/Eqs/fix_bond_react.tex

0 → 100644

+9 −0

Original line number	Diff line number	Diff line
		\documentstyle[12pt]{article}
		\pagestyle{empty}
		\begin{document}

		\begin{eqnarray*}
		k = AT^{n}e^{\frac{-E_{a}}{k_{B}T}}
		\end{eqnarray*}

		\end{document}

doc/src/fix_bond_react.txt

+23 −2

Original line number	Diff line number	Diff line
		@@ -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 to 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,

doc/utils/sphinx-config/false_positives.txt

+1 −0

Original line number	Diff line number	Diff line
		@@ -108,6 +108,7 @@ Archlinux
		arcsin
		arg
		args
		arrhenius
		Arun
		arXiv
		asin

src/USER-MISC/fix_bond_react.cpp

+64 −3

Original line number	Diff line number	Diff line
		@@ -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++;

Admin message