Merge pull request #1477 from jrgissing/bond/react-reaction_topology_overflow (47c6cd86) · Commits · 郑智淋 / lammps

doc/src/Errors_messages.txt

+56 −0

Original line number	Diff line number	Diff line
		@@ -610,6 +610,62 @@ This means there is something invalid about the topology definitions. :dd

		The data file header lists bonds but no bond types. :dd

		{Bond/react: Cannot use fix bond/react with non-molecular systems} :dt

		Only systems with bonds that can be changed can be used. Atom_style
		template does not qualify. :dd

		{Bond/react: Rmax cutoff is longer than pairwise cutoff} :dt

		This is not allowed because bond creation is done using the pairwise
		neighbor list. :dd

		{Bond/react: Molecule template ID for fix bond/react does not exist} :dt

		A valid molecule template must have been created with the molecule
		command. :dd

		{Bond/react: Reaction templates must contain the same number of atoms} :dt

		There should be a one-to-one correspondence between atoms in the
		pre-reacted and post-reacted templates, as specified by the map file. :dd

		{Bond/react: Unknown section in map file} :dt

		Please ensure reaction map files are properly formatted. :dd

		{Bond/react: Atom affected by reaction too close to template edge} :dt

		This means an atom which changes type during the reaction is too close
		to an 'edge' atom defined in the superimpose file. This could cause
		incorrect assignment of bonds, angle, etc. Generally, this means you
		must include more atoms in your templates, such that there are at
		least two atoms between each atom involved in the reaction and an edge
		atom. :dd

		{Bond/react: Fix bond/react needs ghost atoms from farther away} :dt

		This is because a processor needs to superimpose the entire unreacted
		molecule template onto simulation atoms it knows about. The
		comm_modify cutoff command can be used to extend the communication
		range. :dd

		{Bond/react: A deleted atom cannot remain bonded to an atom that is not deleted} :dt

		Self-explanatory. :dd

		{Bond/react special bond generation overflow} :dt

		The number of special bonds per-atom created by a reaction exceeds the
		system setting. See the read_data or create_box command for how to
		specify this value. :dd

		{Bond/react topology/atom exceed system topology/atom} :dt

		The number of bonds, angles etc per-atom created by a reaction exceeds
		the system setting. See the read_data or create_box command for how to
		specify this value. :dd

		{Both restart files must use % or neither} :dt

		Self-explanatory. :dd

doc/src/Errors_warnings.txt

+5 −0

Original line number	Diff line number	Diff line
		@@ -82,6 +82,11 @@ bond/angle/dihedral. LAMMPS computes this by taking the maximum bond
		length, multiplying by the number of bonds in the interaction (e.g. 3
		for a dihedral) and adding a small amount of stretch. :dd

		{Bond/react: An atom in 'react #%d' changes bond connectivity but not atom type} :dt

		You may want to double-check that all atom types are properly assigned
		in the post-reaction template. :dd

		{Both groups in compute group/group have a net charge; the Kspace boundary correction to energy will be non-zero} :dt

		Self-explanatory. :dd

doc/src/fix_bond_react.txt

+41 −30

Original line number	Diff line number	Diff line
		@@ -18,8 +18,8 @@ fix ID group-ID bond/react common_keyword values ...

		ID, group-ID are documented in "fix"_fix.html command. Group-ID is ignored. :ulb,l
		bond/react = style name of this fix command :l
		zero or more common keyword/value pairs may be appended directly after 'bond/react' :l
		these apply to all reaction specifications (below) :l
		the common keyword/values may be appended directly after 'bond/react' :l
		this applies to all reaction specifications (below) :l
		common_keyword = {stabilization} :l
		{stabilization} values = {no} or {yes} {group-ID} {xmax}
		{no} = no reaction site stabilization
		@@ -136,10 +136,12 @@ words, can be customized for each reaction, or reaction step):
		A check for possible new reaction sites is performed every {Nevery}
		timesteps.

		Two conditions must be met for a reaction to occur. First a bonding
		atom pair must be identified. Second, the topology surrounding the
		bonding atom pair must match the topology of the pre-reaction
		template. If both these conditions are met, the reaction site is
		Three physical conditions must be met for a reaction to occur. First,
		a bonding atom pair must be identified within the reaction distance
		cutoffs. Second, the topology surrounding the bonding atom pair must
		match the topology of the pre-reaction template. Finally, any reaction
		constraints listed in the map file (see below) must be satisfied. If
		all of these conditions are met, the reaction site is eligible to be
		modified to match the post-reaction template.

		A bonding atom pair will be identified if several conditions are met.
		@@ -203,14 +205,24 @@ new types must also be defined during the setup of a given simulation.
		A discussion of correctly handling this is also provided on the
		"molecule"_molecule.html command page.

		NOTE: When a reaction occurs, it is possible that the resulting
		topology/atom (e.g. special bonds, dihedrals, etc.) exceeds that of
		the existing system and reaction templates. As when inserting
		molecules, enough space for this increased topology/atom must be
		reserved by using the relevant "extra" keywords to the
		"read_data"_read_data.html or "create_box"_create_box.html commands.

		The map file is a text document with the following format:

		A map file has a header and a body. The header of map file the
		contains one mandatory keyword and three optional keywords. The
		mandatory keyword is 'equivalences' and the optional keywords are
		'edgeIDs' and 'deleteIDs' and 'customIDs':
		contains one mandatory keyword and four optional keywords. The
		mandatory keyword is 'equivalences':

		N {equivalences} = # of atoms N in the reaction molecule templates :pre

		The optional keywords are 'edgeIDs', 'deleteIDs', 'customIDs' and
		'constraints':

		N {equivalences} = # of atoms N in the reaction molecule templates
		N {edgeIDs} = # of edge atoms N in the pre-reacted molecule template
		N {deleteIDs} = # of atoms N that are specified for deletion
		N {customIDs} = # of atoms N that are specified for a custom update
		@@ -244,8 +256,8 @@ A sample map file is given below:

		# this is a map file :pre

		2 edgeIDs
		7 equivalences :pre
		7 equivalences
		2 edgeIDs :pre

		BondingIDs :pre

		@@ -297,26 +309,25 @@ can allow for the possibility of one or more reverse reactions.

		The optional keywords deal with the probability of a given reaction
		occurring as well as the stable equilibration of each reaction site as
		it occurs.
		it occurs:

		The {prob} keyword can affect whether an eligible reaction actually
		occurs. The fraction setting must be a value between 0.0 and 1.0. A
		uniform random number between 0.0 and 1.0 is generated and the
		The {prob} keyword can affect whether or not an eligible reaction
		actually occurs. The fraction setting must be a value between 0.0 and
		1.0. A uniform random number between 0.0 and 1.0 is generated and the
		eligible reaction only occurs if the random number is less than the
		fraction. Up to N reactions are permitted to occur, as optionally
		specified by the {max_rxn} keyword.

		The {stabilize_steps} keyword allows for the specification of how many
		timesteps a reaction site is stabilized before being returned to the
		overall system thermostat.

		In order to produce the most physical behavior, this 'reaction site
		equilibration time' should be tuned to be as small as possible while
		retaining stability for a given system or reaction step. After a
		limited number of case studies, this number has been set to a default
		of 60 timesteps. Ideally, it should be individually tuned for each fix
		reaction step. Note that in some situations, decreasing rather than
		increasing this parameter will result in an increase in stability.
		overall system thermostat. In order to produce the most physical
		behavior, this 'reaction site equilibration time' should be tuned to
		be as small as possible while retaining stability for a given system
		or reaction step. After a limited number of case studies, this number
		has been set to a default of 60 timesteps. Ideally, it should be
		individually tuned for each fix reaction step. Note that in some
		situations, decreasing rather than increasing this parameter will
		result in an increase in stability.

		The {update_edges} keyword can increase the number of atoms whose
		atomic charges are updated, when the pre-reaction template contains
		@@ -324,11 +335,11 @@ edge atoms. When the value is set to 'charges,' all atoms' atomic
		charges are updated to those specified by the post-reaction template,
		including atoms near the edge of reaction templates. When the value is
		set to 'custom,' an additional section must be included in the map
		file that specifies whether to update charges, on a per-atom basis.
		The format of this section is detailed above. Listing a pre-reaction
		atom ID with a value of 'charges' will force the update of the atom's
		charge, even if it is near a template edge. Atoms not near a template
		edge are unaffected by this setting.
		file that specifies whether or not to update charges, on a per-atom
		basis. The format of this section is detailed above. Listing a
		pre-reaction atom ID with a value of 'charges' will force the update
		of the atom's charge, even if it is near a template edge. Atoms not
		near a template edge are unaffected by this setting.

		A few other considerations:

doc/utils/sphinx-config/false_positives.txt

+1 −0

Original line number	Diff line number	Diff line
		@@ -2844,6 +2844,7 @@ unoptimized
		unpadded
		unphysical
		unphysically
		unreacted
		unscaled
		unsets
		unsmoothed

src/USER-MISC/fix_bond_react.cpp

+45 −34

Original line number	Diff line number	Diff line
		@@ -308,8 +308,7 @@ FixBondReact::FixBondReact(LAMMPS lmp, int narg, char *arg) :
		onemol->check_attributes(0);
		twomol->check_attributes(0);
		if (onemol->natoms != twomol->natoms)
		error->all(FLERR,"Post-reacted template must contain the same "
		"number of atoms as the pre-reacted template");
		error->all(FLERR,"Bond/react: Reaction templates must contain the same number of atoms");
		get_molxspecials();
		read(i);
		fclose(fp);
		@@ -324,7 +323,7 @@ FixBondReact::FixBondReact(LAMMPS lmp, int narg, char *arg) :
		delete [] files;

		if (atom->molecular != 1)
		error->all(FLERR,"Cannot use fix bond/react with non-molecular systems");
		error->all(FLERR,"Bond/react: Cannot use fix bond/react with non-molecular systems");

		// check if bonding atoms are 1-2, 1-3, or 1-4 bonded neighbors
		// if so, we don't need non-bonded neighbor list
		@@ -472,6 +471,7 @@ FixBondReact::~FixBondReact()
		delete [] guess_branch;
		delete [] pioneer_count;

		if (group) {
		char **newarg;
		newarg = new char*[2];
		newarg[0] = master_group;
		@@ -484,6 +484,7 @@ FixBondReact::~FixBondReact()
		}
		delete [] newarg;
		}
		}

		/* ---------------------------------------------------------------------- */

		@@ -663,7 +664,7 @@ void FixBondReact::init()
		// check cutoff for iatomtype,jatomtype
		for (int i = 0; i < nreacts; i++) {
		if (force->pair == NULL \|\| cutsq[i][1] > force->pair->cutsq[iatomtype[i]][jatomtype[i]])
		error->all(FLERR,"Fix bond/react cutoff is longer than pairwise cutoff");
		error->all(FLERR,"Bond/react: Fix bond/react cutoff is longer than pairwise cutoff");
		}

		// need a half neighbor list, built every Nevery steps
		@@ -1172,7 +1173,7 @@ void FixBondReact::superimpose_algorithm()
		// let's go ahead and catch the simplest of hangs
		//if (hang_catch > onemol->natoms*4)
		if (hang_catch > atom->nlocal*30) {
		error->one(FLERR,"Excessive iteration of superimpose algorithm");
		error->one(FLERR,"Bond/react: Excessive iteration of superimpose algorithm");
		}
		}
		}
		@@ -1285,7 +1286,7 @@ void FixBondReact::make_a_guess()

		for (int i = 0; i < nxspecial[atom->map(glove[pion][1])][0]; i++) {
		if (atom->map(xspecial[atom->map(glove[pion][1])][i]) < 0) {
		error->all(FLERR,"Fix bond/react needs ghost atoms from further away1"); // parallel issues.
		error->all(FLERR,"Bond/react: Fix bond/react needs ghost atoms from further away1"); // parallel issues.
		}
		if (i_limit_tags[(int)atom->map(xspecial[atom->map(glove[pion][1])][i])] != 0) {
		status = GUESSFAIL;
		@@ -1396,7 +1397,7 @@ void FixBondReact::check_a_neighbor()

		//another check for ghost atoms. perhaps remove the one in make_a_guess
		if (atom->map(glove[(int)onemol_xspecial[pion][neigh]-1][1]) < 0) {
		error->all(FLERR,"Fix bond/react needs ghost atoms from further away2");
		error->all(FLERR,"Bond/react: Fix bond/react needs ghost atoms from further away2");
		}

		for (int j = 0; j < onemol_nxspecial[onemol_xspecial[pion][neigh]-1][0]; j++) {
		@@ -1448,7 +1449,7 @@ void FixBondReact::check_a_neighbor()

		//another check for ghost atoms. perhaps remove the one in make_a_guess
		if (atom->map(glove[(int)onemol_xspecial[pion][neigh]-1][1]) < 0) {
		error->all(FLERR,"Fix bond/react needs ghost atoms from further away3");
		error->all(FLERR,"Bond/react: Fix bond/react needs ghost atoms from further away3");
		}

		for (int ii = 0; ii < onemol_nxspecial[onemol_xspecial[pion][neigh]-1][0]; ii++) {
		@@ -1490,7 +1491,7 @@ void FixBondReact::crosscheck_the_neighbor()
		glove[onemol_xspecial[pion][trace]-1][0] == 0) {

		if (avail_guesses == MAXGUESS) {
		error->warning(FLERR,"Fix bond/react failed because MAXGUESS set too small. ask developer for info");
		error->warning(FLERR,"Bond/react: Fix bond/react failed because MAXGUESS set too small. ask developer for info");
		status = GUESSFAIL;
		return;
		}
		@@ -1559,7 +1560,7 @@ void FixBondReact::inner_crosscheck_loop()

		//another check for ghost atoms. perhaps remove the one in make_a_guess
		if (atom->map(glove[(int)onemol_xspecial[pion][neigh]-1][1]) < 0) {
		error->all(FLERR,"Fix bond/react needs ghost atoms from further away4");
		error->all(FLERR,"Bond/react: Fix bond/react needs ghost atoms from further away4");
		}

		if (guess_branch[avail_guesses-1] == 0) avail_guesses--;
		@@ -1720,7 +1721,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
		// if atoms change types, but aren't landlocked, that's bad
		for (int i = 0; i < twomol->natoms; i++) {
		if (twomol->type[i] != onemol->type[equivalences[i][1][myrxn]-1] && landlocked_atoms[i][myrxn] == 0)
		error->one(FLERR,"Atom affected by reaction too close to template edge");
		error->one(FLERR,"Bond/react: Atom affected by reaction too close to template edge");
		}

		// additionally, if a bond changes type, but neither involved atom is landlocked, bad
		@@ -1736,7 +1737,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
		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");
		error->one(FLERR,"Bond/react: Bond type affected by reaction too close to template edge");
		}
		}
		}
		@@ -1746,7 +1747,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
		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");
		error->one(FLERR,"Bond/react: Bond type affected by reaction too close to template edge");
		}
		}
		}
		@@ -1762,7 +1763,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
		int ii = reverse_equiv[i][1][myrxn] - 1;
		for (int j = 0; j < twomol_nxspecial[ii][0]; j++) {
		if (delete_atoms[equivalences[twomol_xspecial[ii][j]-1][1][myrxn]-1][myrxn] == 0) {
		error->one(FLERR,"A deleted atom cannot remain bonded to an atom that is not deleted");
		error->one(FLERR,"Bond/react: A deleted atom cannot remain bonded to an atom that is not deleted");
		}
		}
		}
		@@ -1773,7 +1774,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
		for (int i = 0; i < twomol->natoms; i++) {
		if (twomol_nxspecial[i][0] != onemol_nxspecial[equivalences[i][1][myrxn]-1][0] && landlocked_atoms[i][myrxn] == 0) {
		char str[128];
		sprintf(str,"An atom in 'react #%d' changes bond connectivity but not atom type",myrxn+1);
		sprintf(str,"Bond/react: An atom in 'react #%d' changes bond connectivity but not atom type",myrxn+1);
		error->warning(FLERR,str);
		break;
		}
		@@ -2261,7 +2262,7 @@ void FixBondReact::update_everything()
		if (landlocked_atoms[j][rxnID] == 1) {
		for (int k = 0; k < nspecial[atom->map(update_mega_glove[jj+1][i])][2]; k++) {
		if (atom->map(special[atom->map(update_mega_glove[jj+1][i])][k]) < 0) {
		error->all(FLERR,"Fix bond/react needs ghost atoms from further away - most likely too many processors");
		error->all(FLERR,"Bond/react: Fix bond/react needs ghost atoms from further away - most likely too many processors");
		}
		}
		}
		@@ -2322,6 +2323,8 @@ void FixBondReact::update_everything()
		nspecial[atom->map(update_mega_glove[jj+1][i])][1]++;
		nspecial[atom->map(update_mega_glove[jj+1][i])][2]++;
		}
		if (nspecial[atom->map(update_mega_glove[jj+1][i])][2] > atom->maxspecial)
		error->one(FLERR,"Bond/react special bond generation overflow");
		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];
		}
		@@ -2383,6 +2386,8 @@ void FixBondReact::update_everything()
		bond_type[atom->map(update_mega_glove[jj+1][i])][insert_num] = twomol->bond_type[j][p];
		bond_atom[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->bond_atom[j][p]-1][1][rxnID]][i];
		num_bond[atom->map(update_mega_glove[jj+1][i])]++;
		if (num_bond[atom->map(update_mega_glove[jj+1][i])] > atom->bond_per_atom)
		error->one(FLERR,"Bond/react topology/atom exceed system topology/atom");
		delta_bonds++;
		}
		}
		@@ -2457,6 +2462,8 @@ void FixBondReact::update_everything()
		angle_atom2[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->angle_atom2[j][p]-1][1][rxnID]][i];
		angle_atom3[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->angle_atom3[j][p]-1][1][rxnID]][i];
		num_angle[atom->map(update_mega_glove[jj+1][i])]++;
		if (num_angle[atom->map(update_mega_glove[jj+1][i])] > atom->angle_per_atom)
		error->one(FLERR,"Bond/react topology/atom exceed system topology/atom");
		delta_angle++;
		}
		}
		@@ -2538,6 +2545,8 @@ void FixBondReact::update_everything()
		dihedral_atom3[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->dihedral_atom3[j][p]-1][1][rxnID]][i];
		dihedral_atom4[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->dihedral_atom4[j][p]-1][1][rxnID]][i];
		num_dihedral[atom->map(update_mega_glove[jj+1][i])]++;
		if (num_dihedral[atom->map(update_mega_glove[jj+1][i])] > atom->dihedral_per_atom)
		error->one(FLERR,"Bond/react topology/atom exceed system topology/atom");
		delta_dihed++;
		}
		}
		@@ -2619,6 +2628,8 @@ void FixBondReact::update_everything()
		improper_atom3[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->improper_atom3[j][p]-1][1][rxnID]][i];
		improper_atom4[atom->map(update_mega_glove[jj+1][i])][insert_num] = update_mega_glove[equivalences[twomol->improper_atom4[j][p]-1][1][rxnID]][i];
		num_improper[atom->map(update_mega_glove[jj+1][i])]++;
		if (num_improper[atom->map(update_mega_glove[jj+1][i])] > atom->improper_per_atom)
		error->one(FLERR,"Bond/react topology/atom exceed system topology/atom");
		delta_imprp++;
		}
		}
		@@ -2716,7 +2727,7 @@ void FixBondReact::read(int myrxn)

		// skip 1st line of file
		eof = fgets(line,MAXLINE,fp);
		if (eof == NULL) error->one(FLERR,"Unexpected end of superimpose file");
		if (eof == NULL) error->one(FLERR,"Bond/react: Unexpected end of superimpose file");

		// read header lines
		// skip blank lines or lines that start with "#"
		@@ -2770,7 +2781,7 @@ void FixBondReact::read(int myrxn)
		DeleteAtoms(line, myrxn);
		} else if (strcmp(keyword,"Constraints") == 0) {
		Constraints(line, myrxn);
		} else error->one(FLERR,"Unknown section in superimpose file");
		} else error->one(FLERR,"Bond/react: Unknown section in map file");

		parse_keyword(1,line,keyword);

		@@ -2778,13 +2789,13 @@ void FixBondReact::read(int myrxn)

		// error check
		if (bondflag == 0 \|\| equivflag == 0)
		error->all(FLERR,"Superimpose file missing BondingIDs or Equivalences section\n");
		error->all(FLERR,"Bond/react: Map file missing BondingIDs or Equivalences section\n");

		if (update_edges_flag[myrxn] == 2 && customedgesflag == 0)
		error->all(FLERR,"Map file must have a Custom Edges section when using 'update_edges custom'\n");
		error->all(FLERR,"Bond/react: Map file must have a Custom Edges section when using 'update_edges custom'\n");

		if (update_edges_flag[myrxn] != 2 && customedgesflag == 1)
		error->all(FLERR,"Specify 'update_edges custom' to include Custom Edges section in map file\n");
		error->all(FLERR,"Bond/react: Specify 'update_edges custom' to include Custom Edges section in map file\n");
		}

		void FixBondReact::EdgeIDs(char *line, int myrxn)
		@@ -2830,7 +2841,7 @@ void FixBondReact::CustomEdges(char *line, int myrxn)
		else if (strcmp(edgemode,"charges") == 0)
		custom_edges[tmp-1][myrxn] = 1;
		else
		error->one(FLERR,"Illegal value in 'Custom Edges' section of map file");
		error->one(FLERR,"Bond/react: Illegal value in 'Custom Edges' section of map file");
		}
		delete [] edgemode;
		}
		@@ -2861,7 +2872,7 @@ void FixBondReact::Constraints(char *line, int myrxn)
		constraints[myrxn][3] = tmp[2]*tmp[2]; // using square of distance
		constraints[myrxn][4] = tmp[3]*tmp[3];
		} else
		error->one(FLERR,"Illegal constraint type in 'Constraints' section of map file");
		error->one(FLERR,"Bond/react: Illegal constraint type in 'Constraints' section of map file");
		}
		delete [] constraint_type;
		}
		@@ -2871,7 +2882,7 @@ void FixBondReact::open(char *file)
		fp = fopen(file,"r");
		if (fp == NULL) {
		char str[128];
		snprintf(str,128,"Cannot open superimpose file %s",file);
		snprintf(str,128,"Bond/react: Cannot open map file %s",file);
		error->one(FLERR,str);
		}
		}
		@@ -2884,7 +2895,7 @@ void FixBondReact::readline(char *line)
		else n = strlen(line) + 1;
		}
		MPI_Bcast(&n,1,MPI_INT,0,world);
		if (n == 0) error->all(FLERR,"Unexpected end of superimpose file");
		if (n == 0) error->all(FLERR,"Bond/react: Unexpected end of map file");
		MPI_Bcast(line,n,MPI_CHAR,0,world);
		}

Admin message