Merge branch 'master' into kim-v2-update

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})
 
 ######################################

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
 

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

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) {

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;
 

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

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];

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");
 

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");
 

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");
 

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);