diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index b3c0a70e4a..418bdd0dba 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -193,6 +193,8 @@ if(PKG_MEAM OR PKG_USER-H5MD OR PKG_USER-QMMM OR PKG_USER-SCAFACOS) enable_language(C) endif() +include_directories(${LAMMPS_SOURCE_DIR}) + # do MPI detection after language activation, if MPI for these language is required find_package(MPI QUIET) option(BUILD_MPI "Build MPI version" ${MPI_FOUND}) @@ -1166,7 +1168,6 @@ set(LAMMPS_STYLE_HEADERS_DIR ${CMAKE_CURRENT_BINARY_DIR}/styles) GenerateStyleHeaders(${LAMMPS_STYLE_HEADERS_DIR}) -include_directories(${LAMMPS_SOURCE_DIR}) include_directories(${LAMMPS_STYLE_HEADERS_DIR}) ###################################### diff --git a/doc/src/fix_bond_react.txt b/doc/src/fix_bond_react.txt index 8e7cb1bdae..f806da3324 100644 --- a/doc/src/fix_bond_react.txt +++ b/doc/src/fix_bond_react.txt @@ -41,7 +41,10 @@ react = mandatory argument indicating new reaction specification :l fraction = initiate reaction with this probability if otherwise eligible seed = random number seed (positive integer) {stabilize_steps} value = timesteps - timesteps = number of timesteps to apply internally created nve/limit.html :pre + timesteps = number of timesteps to apply internally created nve/limit.html + {update_edges} value = {none} or {charges} :l + none = do not update topology near the edges of reaction templates + charges = update atomic charges of all atoms in reaction templates :pre :ule [Examples:] @@ -155,7 +158,17 @@ Some atoms in the pre-reacted template that are not reacting may have missing topology with respect to the simulation. For example, the pre-reacted template may contain an atom that would connect to the rest of a long polymer chain. These are referred to as edge atoms, and -are also specified in the map file. +are also specified in the map file. When the pre-reaction template +contains edge atoms, not all atoms, bonds, charges, etc. specified in +the reaction templates will be updated. Specifically, topology that +involves only atoms that are 'too near' to template edges will not be +updated. The definition of 'too near the edge' depends on which +interactions are defined in the simulation. If the simulation has +defined dihedrals, atoms within two bonds of edge atoms are considered +'too near the edge.' If the simulation defines angles, but not +dihedrals, atoms within one bond of edge atoms are considered 'too +near the edge.' If just bonds are defined, only edge atoms are +considered 'too near the edge.' Note that some care must be taken when a building a molecule template for a given simulation. All atom types in the pre-reacted template @@ -255,6 +268,12 @@ 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. +The {update_edges} keyword can increase the number of atoms whose +atomic charges are updated, when the pre-reaction template contains +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. + 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 @@ -323,7 +342,7 @@ bond/break"_fix_bond_break.html, "fix bond/swap"_fix_bond_swap.html, [Default:] -The option defaults are stabilization = no, prob = 1.0, stabilize_steps = 60 +The option defaults are stabilization = no, prob = 1.0, stabilize_steps = 60, update_edges = none :line diff --git a/lib/latte/Install.py b/lib/latte/Install.py index 3b211858dd..1e1f3040c2 100644 --- a/lib/latte/Install.py +++ b/lib/latte/Install.py @@ -78,7 +78,6 @@ def geturl(url,fname): if which('curl') != None: cmd = 'curl -L -o "%s" %s' % (fname,url) - print(cmd) try: subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True) success = True @@ -87,7 +86,6 @@ def geturl(url,fname): if not success and which('wget') != None: cmd = 'wget -O "%s" %s' % (fname,url) - print(cmd) try: subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True) success = True diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp index b30f1b36c6..705bb8ff70 100644 --- a/src/USER-MISC/fix_bond_react.cpp +++ b/src/USER-MISC/fix_bond_react.cpp @@ -163,6 +163,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : memory->create(seed,nreacts,"bond/react:seed"); memory->create(limit_duration,nreacts,"bond/react:limit_duration"); memory->create(stabilize_steps_flag,nreacts,"bond/react:stabilize_steps_flag"); + memory->create(update_edges_flag,nreacts,"bond/react:update_edges_flag"); memory->create(iatomtype,nreacts,"bond/react:iatomtype"); memory->create(jatomtype,nreacts,"bond/react:jatomtype"); memory->create(ibonding,nreacts,"bond/react:ibonding"); @@ -178,6 +179,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : fraction[i] = 1; seed[i] = 12345; stabilize_steps_flag[i] = 0; + update_edges_flag[i] = 0; // set default limit duration to 60 timesteps limit_duration[i] = 60; reaction_count[i] = 0; @@ -249,6 +251,11 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : limit_duration[rxn] = force->numeric(FLERR,arg[iarg+1]); stabilize_steps_flag[rxn] = 1; iarg += 2; + } else if (strcmp(arg[iarg],"update_edges") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: " + "'update_edges' has too few arguments"); + if (strcmp(arg[iarg+1],"charges") == 0) update_edges_flag[rxn] = 1; + iarg += 2; } else error->all(FLERR,"Illegal fix bond/react command: unknown keyword"); } } @@ -379,6 +386,7 @@ FixBondReact::~FixBondReact() memory->destroy(seed); memory->destroy(limit_duration); memory->destroy(stabilize_steps_flag); + memory->destroy(update_edges_flag); memory->destroy(iatomtype); memory->destroy(jatomtype); @@ -2022,13 +2030,13 @@ void FixBondReact::update_everything() } // update charges and types of landlocked atoms - // here, add check for charge instead of requiring it for (int i = 0; i < update_num_mega; i++) { rxnID = update_mega_glove[0][i]; twomol = atom->molecules[reacted_mol[rxnID]]; for (int j = 0; j < twomol->natoms; j++) { int jj = equivalences[j][1][rxnID]-1; - if (landlocked_atoms[j][rxnID] == 1 && atom->map(update_mega_glove[jj+1][i]) >= 0 && + if ((landlocked_atoms[j][rxnID] == 1 || update_edges_flag[rxnID] == 1) && + atom->map(update_mega_glove[jj+1][i]) >= 0 && atom->map(update_mega_glove[jj+1][i]) < nlocal) { type[atom->map(update_mega_glove[jj+1][i])] = twomol->type[j]; if (twomol->qflag && atom->q_flag) { diff --git a/src/USER-MISC/fix_bond_react.h b/src/USER-MISC/fix_bond_react.h index 08756f8131..b960ec8e73 100644 --- a/src/USER-MISC/fix_bond_react.h +++ b/src/USER-MISC/fix_bond_react.h @@ -58,6 +58,7 @@ class FixBondReact : public Fix { tagint lastcheck; int stabilization_flag; int *stabilize_steps_flag; + int *update_edges_flag; int status; int *groupbits; diff --git a/src/comm_brick.cpp b/src/comm_brick.cpp index e96b0f38e3..bf7f3e2f8f 100644 --- a/src/comm_brick.cpp +++ b/src/comm_brick.cpp @@ -1477,8 +1477,10 @@ void CommBrick::free_multi() void *CommBrick::extract(const char *str, int &dim) { + dim = 0; if (strcmp(str,"localsendlist") == 0) { int i, iswap, isend; + dim = 1; if (!localsendlist) memory->create(localsendlist,atom->nlocal,"comm:localsendlist"); else diff --git a/src/compute_adf.cpp b/src/compute_adf.cpp index 6e6239f076..885a051b5b 100644 --- a/src/compute_adf.cpp +++ b/src/compute_adf.cpp @@ -46,14 +46,17 @@ enum{DEGREE, RADIAN, COSINE}; ComputeADF::ComputeADF(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), ilo(NULL), ihi(NULL), jlo(NULL), jhi(NULL), klo(NULL), khi(NULL), - iatomcount(NULL), iatomcountall(NULL), hist(NULL), histall(NULL), - iatomflag(NULL), - jatomflag(NULL), rcutinnerj(NULL), rcutouterj(NULL), - katomflag(NULL), rcutinnerk(NULL), rcutouterk(NULL), - maxjatom(NULL), numjatom(NULL), neighjatom(NULL), - maxkatom(NULL), numkatom(NULL), neighkatom(NULL), - maxjkatom(NULL), numjkatom(NULL), neighjkatom(NULL), bothjkatom(NULL) + rcutinnerj(NULL), rcutinnerk(NULL), + rcutouterj(NULL), rcutouterk(NULL), + list(NULL), + iatomcount(NULL), iatomcountall(NULL), iatomflag(NULL), + maxjatom(NULL), maxkatom(NULL), + numjatom(NULL), numkatom(NULL), + neighjatom(NULL),neighkatom(NULL), + jatomflag(NULL), katomflag(NULL), + maxjkatom(NULL), numjkatom(NULL), + neighjkatom(NULL), bothjkatom(NULL), delrjkatom(NULL) { int nargsperadf = 7; @@ -358,9 +361,9 @@ void ComputeADF::init_list(int /*id*/, NeighList *ptr) void ComputeADF::compute_array() { int i,j,k,m,ii,jj,jatom,katom,jk,jjj,kkk; - int inum,jnum,itype,jtype,ibin,ihisto; + int inum,jnum,itype,jtype,ibin; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; - int *ilist,*jlist,*klist,*numneigh,**firstneigh; + int *ilist,*jlist,*numneigh,**firstneigh; double factor_lj,factor_coul; double delr1[3],delr2[3],rinv1,rinv2,rinv12,cs,theta; @@ -394,11 +397,9 @@ void ComputeADF::compute_array() double **x = atom->x; int *type = atom->type; int *mask = atom->mask; - int nlocal = atom->nlocal; double *special_coul = force->special_coul; double *special_lj = force->special_lj; - int newton_pair = force->newton_pair; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; diff --git a/src/compute_angle_local.cpp b/src/compute_angle_local.cpp index 95faad54ad..cf5c9e808d 100644 --- a/src/compute_angle_local.cpp +++ b/src/compute_angle_local.cpp @@ -38,7 +38,7 @@ enum{THETA,ENG,VARIABLE}; ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), - bstyle(NULL), vstr(NULL), vvar(NULL), tstr(NULL), vlocal(NULL), alocal(NULL) + bstyle(NULL), vvar(NULL), tstr(NULL), vstr(NULL), vlocal(NULL), alocal(NULL) { if (narg < 4) error->all(FLERR,"Illegal compute angle/local command"); diff --git a/src/compute_bond_local.cpp b/src/compute_bond_local.cpp index 6a179cf1b4..72ce2b5f8d 100644 --- a/src/compute_bond_local.cpp +++ b/src/compute_bond_local.cpp @@ -39,7 +39,7 @@ enum{DIST,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE,VARIABLE}; ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), - bstyle(NULL), vstr(NULL), vvar(NULL), dstr(NULL), vlocal(NULL), alocal(NULL) + bstyle(NULL), vvar(NULL), dstr(NULL), vstr(NULL), vlocal(NULL), alocal(NULL) { if (narg < 4) error->all(FLERR,"Illegal compute bond/local command"); diff --git a/src/compute_dihedral_local.cpp b/src/compute_dihedral_local.cpp index 7444630090..919081236c 100644 --- a/src/compute_dihedral_local.cpp +++ b/src/compute_dihedral_local.cpp @@ -40,7 +40,7 @@ enum{PHI,VARIABLE}; ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), - bstyle(NULL), vstr(NULL), vvar(NULL), pstr(NULL), vlocal(NULL), alocal(NULL) + bstyle(NULL), vvar(NULL), pstr(NULL), vstr(NULL), vlocal(NULL), alocal(NULL) { if (narg < 4) error->all(FLERR,"Illegal compute dihedral/local command"); diff --git a/src/deprecated.cpp b/src/deprecated.cpp index 2f2282f07c..b937482669 100644 --- a/src/deprecated.cpp +++ b/src/deprecated.cpp @@ -36,7 +36,7 @@ static void writemsg(LAMMPS *lmp, const char *msg, int abend=1) /* ---------------------------------------------------------------------- */ -void Deprecated::command(int narg, char **arg) +void Deprecated::command(int /* narg */, char ** /* arg */) { if (strcmp(input->command,"DEPRECATED") == 0) { writemsg(lmp,"\nCommand 'DEPRECATED' is a dummy command\n\n",0);