bond/react RMSD constraint: manual rebase (bc49e854) · Commits · 郑智淋 / lammps

src/USER-REACTION/fix_bond_react.cpp

+58 −73

Original line number	Diff line number	Diff line
		@@ -19,6 +19,7 @@ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu)
		#include <mpi.h>
		#include <cmath>
		#include <cstring>
		#include <string>
		#include "update.h"
		#include "modify.h"
		#include "respa.h"
		@@ -32,6 +33,7 @@ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu)
		#include "neigh_list.h"
		#include "neigh_request.h"
		#include "random_mars.h"
		#include "reset_mol_ids.h"
		#include "molecule.h"
		#include "group.h"
		#include "citeme.h"
		@@ -41,6 +43,7 @@ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu)
		#include "error.h"
		#include "input.h"
		#include "variable.h"
		#include "fmt/format.h"
		#include "superpose3d.h"

		#include <algorithm>
		@@ -91,6 +94,7 @@ FixBondReact::FixBondReact(LAMMPS lmp, int narg, char *arg) :
		fix1 = NULL;
		fix2 = NULL;
		fix3 = NULL;
		reset_mol_ids = NULL;

		if (narg < 8) error->all(FLERR,"Illegal fix bond/react command: "
		"too few arguments");
		@@ -143,7 +147,8 @@ FixBondReact::FixBondReact(LAMMPS lmp, int narg, char *arg) :

		int iarg = 3;
		stabilization_flag = 0;
		int num_common_keywords = 1;
		reset_mol_ids_flag = 1;
		int num_common_keywords = 2;
		for (int m = 0; m < num_common_keywords; m++) {
		if (strcmp(arg[iarg],"stabilization") == 0) {
		if (strcmp(arg[iarg+1],"no") == 0) {
		@@ -161,11 +166,23 @@ FixBondReact::FixBondReact(LAMMPS lmp, int narg, char *arg) :
		nve_limit_xmax = arg[iarg+3];
		iarg += 4;
		}
		} else if (strcmp(arg[iarg],"reset_mol_ids") == 0) {
		if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: "
		"'reset_mol_ids' keyword has too few arguments");
		if (strcmp(arg[iarg+1],"yes") == 0) ; // default
		if (strcmp(arg[iarg+1],"no") == 0) reset_mol_ids_flag = 0;
		iarg += 2;
		} else if (strcmp(arg[iarg],"react") == 0) {
		break;
		} else error->all(FLERR,"Illegal fix bond/react command: unknown keyword");
		}

		if (reset_mol_ids_flag) {
		delete reset_mol_ids;
		reset_mol_ids = new ResetMolIDs(lmp);
		reset_mol_ids->create_computes(id,group->names[igroup]);
		}

		// set up common variables as vectors of length 'nreacts'
		// nevery, cutoff, onemol, twomol, superimpose file

		@@ -228,11 +245,11 @@ FixBondReact::FixBondReact(LAMMPS lmp, int narg, char *arg) :

		int n = strlen(arg[iarg]) + 1;
		if (n > MAXLINE) error->all(FLERR,"Reaction name (react-ID) is too long (limit: 256 characters)");
		strncpy(rxn_name[rxn],arg[iarg++],n);
		strcpy(rxn_name[rxn],arg[iarg++]);

		int igroup = group->find(arg[iarg++]);
		if (igroup == -1) error->all(FLERR,"Could not find fix group ID");
		groupbits[rxn] = group->bitmask[igroup];
		int groupid = group->find(arg[iarg++]);
		if (groupid == -1) error->all(FLERR,"Could not find fix group ID");
		groupbits[rxn] = group->bitmask[groupid];

		if (strncmp(arg[iarg],"v_",2) == 0) {
		n = strlen(&arg[iarg][2]) + 1;
		@@ -498,6 +515,8 @@ FixBondReact::~FixBondReact()
		}
		delete [] random;

		delete reset_mol_ids;

		memory->destroy(partner);
		memory->destroy(finalpartner);
		memory->destroy(ncreate);
		@@ -576,17 +595,11 @@ FixBondReact::~FixBondReact()
		delete [] set;

		if (group) {
		char **newarg;
		newarg = new char*[2];
		newarg[0] = master_group;
		newarg[1] = (char *) "delete";
		group->assign(2,newarg);
		group->assign(std::string(master_group) + " delete");
		if (stabilization_flag == 1) {
		newarg[0] = exclude_group;
		group->assign(2,newarg);
		group->assign(std::string(exclude_group) + " delete");
		delete [] exclude_group;
		}
		delete [] newarg;
		}
		}

		@@ -609,59 +622,38 @@ it will have the name 'i_limit_tags' and will be intitialized to 0 (not in group

		void FixBondReact::post_constructor()
		{
		int len;
		// let's add the limit_tags per-atom property fix
		int len = strlen("bond_react_props_internal") + 1;
		id_fix2 = new char[len];
		strcpy(id_fix2,"bond_react_props_internal");
		std::string cmd = std::string("bond_react_props_internal");
		id_fix2 = new char[cmd.size()+1];
		strcpy(id_fix2,cmd.c_str());

		int ifix = modify->find_fix(id_fix2);
		if (ifix == -1) {
		char *newarg = new char[7];
		newarg[0] = (char *) "bond_react_props_internal";
		newarg[1] = (char *) "all"; // group ID is ignored
		newarg[2] = (char *) "property/atom";
		newarg[3] = (char *) "i_limit_tags";
		newarg[4] = (char *) "i_react_tags";
		newarg[5] = (char *) "ghost";
		newarg[6] = (char *) "yes";
		modify->add_fix(7,newarg);
		delete [] newarg;
		cmd += std::string(" all property/atom i_limit_tags i_react_tags ghost yes");
		modify->add_fix(cmd);
		}

		// create master_group if not already existing
		// NOTE: limit_tags and react_tags automaticaly intitialized to zero (unless read from restart)
		group->find_or_create(master_group);
		char **newarg;
		newarg = new char*[5];
		newarg[0] = master_group;
		newarg[1] = (char *) "dynamic";
		newarg[2] = (char *) "all";
		newarg[3] = (char *) "property";
		newarg[4] = (char *) "limit_tags";
		group->assign(5,newarg);
		delete [] newarg;
		cmd = fmt::format("{} dynamic all property limit_tags",master_group);
		group->assign(cmd);

		if (stabilization_flag == 1) {
		int igroup = group->find(exclude_group);
		int groupid = group->find(exclude_group);
		// create exclude_group if not already existing, or use as parent group if static
		if (igroup == -1 \|\| group->dynamic[igroup] == 0) {
		if (groupid == -1 \|\| group->dynamic[groupid] == 0) {
		// create stabilization per-atom property
		len = strlen("bond_react_stabilization_internal") + 1;
		id_fix3 = new char[len];
		strcpy(id_fix3,"bond_react_stabilization_internal");
		cmd = std::string("bond_react_stabilization_internal");
		id_fix3 = new char[cmd.size()+1];
		strcpy(id_fix3,cmd.c_str());

		ifix = modify->find_fix(id_fix3);
		if (ifix == -1) {
		char *newarg = new char[6];
		newarg[0] = (char *) id_fix3;
		newarg[1] = (char *) "all"; // group ID is ignored
		newarg[2] = (char *) "property/atom";
		newarg[3] = (char *) "i_statted_tags";
		newarg[4] = (char *) "ghost";
		newarg[5] = (char *) "yes";
		modify->add_fix(6,newarg);
		cmd += std::string(" all property/atom i_statted_tags ghost yes");
		modify->add_fix(cmd);
		fix3 = modify->fix[modify->nfix-1];
		delete [] newarg;
		}

		len = strlen("statted_tags") + 1;
		@@ -681,16 +673,11 @@ void FixBondReact::post_constructor()
		strcat(exclude_group,"_REACT");

		group->find_or_create(exclude_group);
		char **newarg;
		newarg = new char*[5];
		newarg[0] = exclude_group;
		newarg[1] = (char *) "dynamic";
		if (igroup == -1) newarg[2] = (char *) "all";
		else newarg[2] = (char *) exclude_PARENT_group;
		newarg[3] = (char *) "property";
		newarg[4] = (char *) "statted_tags";
		group->assign(5,newarg);
		delete [] newarg;
		if (groupid == -1)
		cmd = fmt::format("{} dynamic all property statted_tags",exclude_group);
		else
		cmd = fmt::format("{} dynamic {} property statted_tags",exclude_group,exclude_PARENT_group);
		group->assign(cmd);
		delete [] exclude_PARENT_group;

		// on to statted_tags (system-wide thermostat)
		@@ -738,21 +725,16 @@ void FixBondReact::post_constructor()


		// let's create a new nve/limit fix to limit newly reacted atoms
		len = strlen("bond_react_MASTER_nve_limit") + 1;
		id_fix1 = new char[len];
		strcpy(id_fix1,"bond_react_MASTER_nve_limit");
		cmd = std::string("bond_react_MASTER_nve_limit");
		id_fix1 = new char[cmd.size()+1];
		strcpy(id_fix1,cmd.c_str());

		ifix = modify->find_fix(id_fix1);

		if (ifix == -1) {
		char *newarg = new char[4];
		newarg[0] = id_fix1;
		newarg[1] = master_group;
		newarg[2] = (char *) "nve/limit";
		newarg[3] = nve_limit_xmax;
		modify->add_fix(4,newarg);
		cmd += fmt::format(" {} nve/limit {}",master_group,nve_limit_xmax);
		modify->add_fix(cmd);
		fix1 = modify->fix[modify->nfix-1];
		delete [] newarg;
		}
		}
		}
		@@ -2136,7 +2118,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
		for (int i = 0; i < twomol->natoms; i++) {
		if (twomol->type[i] != onemol->type[equivalences[i][1][myrxn]-1] && landlocked_atoms[i][myrxn] == 0) {
		char str[128];
		snprintf(str,128,"Bond/react: Atom affected by reaction %s too close to template edge",rxn_name[myrxn]);
		snprintf(str,128,"Bond/react: Atom type affected by reaction %s too close to template edge",rxn_name[myrxn]);
		error->all(FLERR,str);
		}
		}
		@@ -2155,7 +2137,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
		if (onemol_batom == equivalences[twomol_atomj-1][1][myrxn]) {
		if (twomol->bond_type[i][j] != onemol->bond_type[onemol_atomi-1][m]) {
		char str[128];
		snprintf(str,128,"Bond/react: Atom affected by reaction %s too close to template edge",rxn_name[myrxn]);
		snprintf(str,128,"Bond/react: Bond type affected by reaction %s too close to template edge",rxn_name[myrxn]);
		error->all(FLERR,str);
		}
		}
		@@ -2167,7 +2149,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
		if (onemol_batom == equivalences[i][1][myrxn]) {
		if (twomol->bond_type[i][j] != onemol->bond_type[onemol_atomj-1][m]) {
		char str[128];
		snprintf(str,128,"Bond/react: Atom affected by reaction %s too close to template edge",rxn_name[myrxn]);
		snprintf(str,128,"Bond/react: Bond type affected by reaction %s too close to template edge",rxn_name[myrxn]);
		error->all(FLERR,str);
		}
		}
		@@ -2578,7 +2560,7 @@ void FixBondReact::ghost_glovecast()
		}

		/* ----------------------------------------------------------------------
		update charges, types, special lists and all topology
		update molecule IDs, charges, types, special lists and all topology
		------------------------------------------------------------------------- */

		void FixBondReact::update_everything()
		@@ -3117,6 +3099,9 @@ void FixBondReact::update_everything()
		atom->natoms -= ndel;
		// done deleting atoms

		// reset mol ids
		if (reset_mol_ids_flag) reset_mol_ids->reset();

		// something to think about: this could done much more concisely if
		// all atom-level info (bond,angles, etc...) were kinda inherited from a common data struct --JG

src/USER-REACTION/fix_bond_react.h

+5 −1

Original line number	Diff line number	Diff line
		@@ -61,6 +61,7 @@ class FixBondReact : public Fix {
		int max_rxn,nlocalskips,*nghostlyskips;
		tagint lastcheck;
		int stabilization_flag;
		int reset_mol_ids_flag;
		int custom_exclude_flag;
		int *stabilize_steps_flag;
		int *update_edges_flag;
		@@ -93,6 +94,7 @@ class FixBondReact : public Fix {
		class RanMars **random; // random number for 'prob' keyword
		class RanMars **rrhandom; // random number for Arrhenius constraint
		class NeighList *list;
		class ResetMolIDs *reset_mol_ids; // class for resetting mol IDs

		int reacted_mol,unreacted_mol;
		int *limit_duration; // indicates how long to relax
		@@ -240,8 +242,10 @@ E: Bond/react: Invalid template atom ID in map file
		Atom IDs in molecule templates range from 1 to the number of atoms in the template.

		E or W: Bond/react: Atom affected by reaction %s too close to template edge
		Bond/react: Atom type affected by reaction %s too close to template edge
		Bond/react: Bond type affected by reaction %s too close to template edge

		This means an atom which changes type or connectivity during the
		This means an atom (or bond) that changes type or connectivity during the
		reaction is too close to an 'edge' atom defined in the map file. This
		could cause incorrect assignment of bonds, angle, etc. Generally, this
		means you must include more atoms in your templates, such that there

Admin message