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

Merge branch 'master' into collected-small-changes

This commit is contained in:
Axel Kohlmeyer 2020-08-26 11:49:29 -04:00
parent 14b66d1f84 e2fd95d8ed
commit c8af729701
No known key found for this signature in database
GPG Key ID: D9B44E93BF0C375A
6 changed files with 160 additions and 70 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
doc/src/Errors_messages.rst
View File

@ -502,7 +502,8 @@ Doc page with :doc:`WARNING messages <Errors_warnings>`
*Bond/react: Unknown section in map file*
Please ensure reaction map files are properly formatted.
*Bond/react: Atom affected by reaction too close to template edge*
*Bond/react: Atom type affected by reaction too close to template edge*
*Bond/react: Bond type affected by reaction too close to template edge*
This means an atom which 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.

31
doc/src/fix_bond_react.rst
View File

@ -14,19 +14,22 @@ Syntax
react react-ID react-group-ID Nevery Rmin Rmax template-ID(pre-reacted) template-ID(post-reacted) map_file individual_keyword values ...
...
* ID, group-ID are documented in :doc:`fix <fix>` command. Group-ID is ignored.
* ID, group-ID are documented in :doc:`fix <fix>` command.
* bond/react = style name of this fix command
* the common keyword/values may be appended directly after 'bond/react'
* this applies to all reaction specifications (below)
* common_keyword = *stabilization*
* common_keyword = *stabilization* or *reset_mol_ids*
.. parsed-literal::
*stabilization* values = *no* or *yes* *group-ID* *xmax*
*no* = no reaction site stabilization
*no* = no reaction site stabilization (default)
*yes* = perform reaction site stabilization
*group-ID* = user-assigned prefix for the dynamic group of atoms not currently involved in a reaction
*xmax* = xmax value that is used by an internally-created :doc:`nve/limit <fix_nve_limit>` integrator
*reset_mol_ids* values = *yes* or *no*
*yes* = update molecule IDs based on new global topology (default)
*no* = do not update molecule IDs
* react = mandatory argument indicating new reaction specification
* react-ID = user-assigned name for the reaction
@ -50,9 +53,9 @@ Syntax
*stabilize_steps* value = timesteps
timesteps = number of timesteps to apply the internally-created :doc:`nve/limit <fix_nve_limit>` fix to reacting atoms
*update_edges* value = *none* or *charges* or *custom*
none = do not update topology near the edges of reaction templates
charges = update atomic charges of all atoms in reaction templates
custom = force the update of user-specified atomic charges
*none* = do not update topology near the edges of reaction templates
*charges* = update atomic charges of all atoms in reaction templates
*custom* = force the update of user-specified atomic charges
Examples
""""""""
@ -154,6 +157,13 @@ due to the internal dynamic grouping performed by fix bond/react.
If the group-ID is an existing static group, react-group-IDs
should also be specified as this static group, or a subset.
The *reset_mol_ids* keyword invokes the :doc:`reset_mol_ids <reset_mol_ids>`
command after a reaction occurs, to ensure that molecule IDs are
consistent with the new bond topology. The group-ID used for
:doc:`reset_mol_ids <reset_mol_ids>` is the group-ID for this fix.
Resetting molecule IDs is necessarily a global operation, and so can
be slow for very large systems.
The following comments pertain to each *react* argument (in other
words, can be customized for each reaction, or reaction step):
@ -203,9 +213,10 @@ surrounding topology. As described below, the bonding atom pairs of
the pre-reacted template are specified by atom ID in the map file. The
pre-reacted molecule template should contain as few atoms as possible
while still completely describing the topology of all atoms affected
by the reaction. For example, if the force field contains dihedrals,
the pre-reacted template should contain any atom within three bonds of
reacting atoms.
by the reaction (which includes all atoms that change atom type or
connectivity, and all bonds that change bond type). For example, if
the force field contains dihedrals, the pre-reacted template should
contain any atom within three bonds of reacting atoms.
Some atoms in the pre-reacted template that are not reacting may have
missing topology with respect to the simulation. For example, the
@ -554,7 +565,7 @@ Default
"""""""
The option defaults are stabilization = no, prob = 1.0, stabilize_steps = 60,
update_edges = none
reset_mol_ids = yes, update_edges = none
----------

44
src/USER-REACTION/fix_bond_react.cpp
View File

@ -33,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"
@ -92,6 +93,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");
@ -144,7 +146,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) {
@ -162,11 +165,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
@ -229,11 +244,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;
@ -499,6 +514,8 @@ FixBondReact::~FixBondReact()
}
delete [] random;
delete reset_mol_ids;
memory->destroy(partner);
memory->destroy(finalpartner);
memory->destroy(ncreate);
@ -623,9 +640,9 @@ void FixBondReact::post_constructor()
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
cmd = std::string("bond_react_stabilization_internal");
id_fix3 = new char[cmd.size()+1];
@ -655,7 +672,7 @@ void FixBondReact::post_constructor()
strcat(exclude_group,"_REACT");
group->find_or_create(exclude_group);
if (igroup == -1)
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);
@ -2061,7 +2078,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);
}
}
@ -2080,7 +2097,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);
}
}
@ -2092,7 +2109,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);
}
}
@ -2503,7 +2520,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()
@ -3042,6 +3059,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

6
src/USER-REACTION/fix_bond_react.h
View File

@ -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

131
src/reset_mol_ids.cpp
View File

@ -17,7 +17,6 @@
#include "reset_mol_ids.h"
#include <mpi.h>
#include <string>
#include "atom.h"
#include "domain.h"
#include "comm.h"
@ -33,10 +32,29 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ResetMolIDs::ResetMolIDs(LAMMPS *lmp) : Pointers(lmp) {}
ResetMolIDs::ResetMolIDs(LAMMPS *lmp) : Pointers(lmp) {
cfa = NULL;
cca = NULL;
// default settings
compressflag = 1;
singleflag = 0;
offset = -1;
}
/* ---------------------------------------------------------------------- */
ResetMolIDs::~ResetMolIDs()
{
modify->delete_compute(idfrag);
if (compressflag) modify->delete_compute(idchunk);
}
/* ----------------------------------------------------------------------
called as reset_mol_ids command in input script
------------------------------------------------------------------------- */
void ResetMolIDs::command(int narg, char **arg)
{
if (domain->box_exist == 0)
@ -49,13 +67,7 @@ void ResetMolIDs::command(int narg, char **arg)
// process args
if (narg < 1) error->all(FLERR,"Illegal reset_mol_ids command");
int igroup = group->find(arg[0]);
if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID");
int groupbit = group->bitmask[igroup];
int compressflag = 1;
int singleflag = 0;
tagint offset = -1;
char *groupid = arg[0];
int iarg = 1;
while (iarg < narg) {
@ -86,20 +98,6 @@ void ResetMolIDs::command(int narg, char **arg)
MPI_Barrier(world);
double time1 = MPI_Wtime();
// create instances of compute fragment/atom, compute reduce (if needed),
// and compute chunk/atom. all use the group-ID for this command
const std::string idfrag = "reset_mol_ids_FRAGMENT_ATOM";
if (singleflag)
modify->add_compute(fmt::format("{} {} fragment/atom single yes",idfrag,arg[0]));
else
modify->add_compute(fmt::format("{} {} fragment/atom single no",idfrag,arg[0]));
const std::string idchunk = "reset_mol_ids_CHUNK_ATOM";
if (compressflag)
modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes",
idchunk,arg[0]));
// initialize system since comm->borders() will be invoked
lmp->init();
@ -116,12 +114,73 @@ void ResetMolIDs::command(int narg, char **arg)
comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
// create computes
create_computes((char *) "COMMAND",groupid);
// reset molecule IDs
reset();
// total time
MPI_Barrier(world);
if (comm->me == 0) {
if (nchunk < 0)
utils::logmesg(lmp,fmt::format(" number of new molecule IDs = unknown\n"));
else
utils::logmesg(lmp,fmt::format(" number of new molecule IDs = {}\n",nchunk));
utils::logmesg(lmp,fmt::format(" reset_mol_ids CPU = {:.3f} seconds\n",
MPI_Wtime()-time1));
}
}
/* ----------------------------------------------------------------------
create computes used by reset_mol_ids
------------------------------------------------------------------------- */
void ResetMolIDs::create_computes(char *fixid, char *groupid)
{
int igroup = group->find(groupid);
if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID");
groupbit = group->bitmask[igroup];
// create instances of compute fragment/atom, compute reduce (if needed),
// and compute chunk/atom. all use the group-ID for this command.
// 'fixid' allows for creating independent instances of the computes
idfrag = fmt::format("{}_reset_mol_ids_FRAGMENT_ATOM",fixid);
if (singleflag)
modify->add_compute(fmt::format("{} {} fragment/atom single yes",idfrag,groupid));
else
modify->add_compute(fmt::format("{} {} fragment/atom single no",idfrag,groupid));
idchunk = fmt::format("{}_reset_mol_ids_CHUNK_ATOM",fixid);
if (compressflag)
modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes",
idchunk,groupid));
int icompute = modify->find_compute(idfrag);
cfa = (ComputeFragmentAtom *) modify->compute[icompute];
if (compressflag) {
icompute = modify->find_compute(idchunk);
cca = (ComputeChunkAtom *) modify->compute[icompute];
}
}
/* ----------------------------------------------------------------------
called from command() and directly from fixes that update molecules
------------------------------------------------------------------------- */
void ResetMolIDs::reset()
{
// invoke peratom method of compute fragment/atom
// walks bond connectivity and assigns each atom a fragment ID
// if singleflag = 0, atoms w/out bonds will be assigned fragID = 0
int icompute = modify->find_compute(idfrag);
ComputeFragmentAtom *cfa = (ComputeFragmentAtom *) modify->compute[icompute];
cfa->compute_peratom();
double *fragIDs = cfa->vector_atom;
@ -138,15 +197,13 @@ void ResetMolIDs::command(int narg, char **arg)
// if compressflag = 0, done
// set nchunk = -1 since cannot easily determine # of new molecule IDs
int nchunk = -1;
nchunk = -1;
// if compressflag = 1, invoke peratom method of compute chunk/atom
// will compress new molecule IDs to be contiguous 1 to Nmol
// NOTE: use of compute chunk/atom limits Nmol to a 32-bit int
if (compressflag) {
icompute = modify->find_compute(idchunk);
ComputeChunkAtom *cca = (ComputeChunkAtom *) modify->compute[icompute];
cca->compute_peratom();
double *chunkIDs = cca->vector_atom;
nchunk = cca->nchunk;
@ -193,22 +250,4 @@ void ResetMolIDs::command(int narg, char **arg)
}
}
}
// clean up
modify->delete_compute(idfrag);
if (compressflag) modify->delete_compute(idchunk);
// total time
MPI_Barrier(world);
if (comm->me == 0) {
if (nchunk < 0)
utils::logmesg(lmp,fmt::format(" number of new molecule IDs = unknown\n"));
else
utils::logmesg(lmp,fmt::format(" number of new molecule IDs = {}\n",nchunk));
utils::logmesg(lmp,fmt::format(" reset_mol_ids CPU = {:.3f} seconds\n",
MPI_Wtime()-time1));
}
}

15
src/reset_mol_ids.h
View File

@ -21,13 +21,28 @@ CommandStyle(reset_mol_ids,ResetMolIDs)
#define LMP_RESET_MOL_IDS_H
#include "pointers.h"
#include <string>
namespace LAMMPS_NS {
class ResetMolIDs : protected Pointers {
public:
ResetMolIDs(class LAMMPS *);
~ResetMolIDs();
void command(int, char **);
void create_computes(char *, char *);
void reset();
private:
std::string idfrag, idchunk;
int nchunk;
int groupbit;
int compressflag; // 1 = contiguous values for new IDs
int singleflag; // 0 = mol IDs of single atoms set to 0
tagint offset; // offset for contiguous mol ID values
class ComputeFragmentAtom *cfa;
class ComputeChunkAtom *cca;
};
}