adjustments to balancing weights and factors, also XOR op for formulas, if, dump_modify thresh

Merge pull request #193 from akohlmey/eam-bugfix
bugfix for eam/alloy/omp and eam/fs/omp
2016-10-05 15:46:20 -06:00 · 2016-10-05 10:54:36 -06:00 · 2016-10-05 10:54:08 -06:00 · 2016-10-05 10:33:39 -06:00 · 2016-10-04 14:33:26 -04:00 · 2016-10-04 13:57:21 -04:00
83 changed files with 7943 additions and 8349 deletions
--- a/doc/src/Manual.txt
+++ b/doc/src/Manual.txt
@ -1,7 +1,7 @@
 <!-- HTML_ONLY -->
 <HEAD>
 <TITLE>LAMMPS Users Manual</TITLE>
-<META NAME="docnumber" CONTENT="29 Sep 2016 version">
+<META NAME="docnumber" CONTENT="5 Oct 2016 version">
 <META NAME="author" CONTENT="http://lammps.sandia.gov - Sandia National Laboratories">
 <META NAME="copyright" CONTENT="Copyright (2003) Sandia Corporation.  This software and manual is distributed under the GNU General Public License.">
 </HEAD>
@ -21,7 +21,7 @@
 <H1></H1>

 LAMMPS Documentation :c,h3
-29 Sep 2016 version :c,h4
+5 Oct 2016 version :c,h4

 Version info: :h4

--- a/doc/src/Section_commands.txt
+++ b/doc/src/Section_commands.txt
@ -896,7 +896,7 @@ KOKKOS, o = USER-OMP, t = OPT.
 "lubricate/poly (o)"_pair_lubricate.html,
 "lubricateU"_pair_lubricateU.html,
 "lubricateU/poly"_pair_lubricateU.html,
-"meam (o)"_pair_meam.html,
+"meam"_pair_meam.html,
 "mie/cut (o)"_pair_mie.html,
 "morse (got)"_pair_morse.html,
 "nb3b/harmonic (o)"_pair_nb3b_harmonic.html,
@ -956,7 +956,7 @@ package"_Section_start.html#start_3.
 "lj/sdk/coul/long (go)"_pair_sdk.html,
 "lj/sdk/coul/msm (o)"_pair_sdk.html,
 "lj/sf (o)"_pair_lj_sf.html,
-"meam/spline"_pair_meam_spline.html,
+"meam/spline (o)"_pair_meam_spline.html,
 "meam/sw/spline"_pair_meam_sw_spline.html,
 "mgpt"_pair_mgpt.html,
 "morse/smooth/linear"_pair_morse.html,
--- a/doc/src/balance.txt
+++ b/doc/src/balance.txt
@ -319,14 +319,16 @@ accurately would be impractical and slow down the computation.
 Instead the {weight} keyword implements several ways to influence the
 per-particle weights empirically by properties readily available or
 using the user's knowledge of the system.  Note that the absolute
-value of the weights are not important; their ratio is what is used to
-assign particles to processors.  A particle with a weight of 2.5 is
-assumed to require 5x more computational than a particle with a weight
-of 0.5.
+value of the weights are not important; only their relative ratios
+affect which particle is assigned to which processor.  A particle with
+a weight of 2.5 is assumed to require 5x more computational than a
+particle with a weight of 0.5.  For all the options below the weight
+assigned to a particle must be a positive value; an error will be be
+generated if a weight is <= 0.0.

 Below is a list of possible weight options with a short description of
 their usage and some example scenarios where they might be applicable.
-It is possible to apply multiple weight flags and the weightins they
+It is possible to apply multiple weight flags and the weightings they
 induce will be combined through multiplication.  Most of the time,
 however, it is sufficient to use just one method.

@ -346,13 +348,24 @@ the computational cost for each group remains constant over time.
 This is a purely empirical weighting, so a series test runs to tune
 the assigned weight factors for optimal performance is recommended.

-The {neigh} weight style assigns a weight to each particle equal to
-its number of neighbors divided by the avergage number of neighbors
-for all particles.  The {factor} setting is then appied as an overall
-scale factor to all the {neigh} weights which allows tuning of the
-impact of this style.  A {factor} smaller than 1.0 (e.g. 0.8) often
-results in the best performance, since the number of neighbors is
-likely to overestimate the ideal weight.
+The {neigh} weight style assigns the same weight to each particle
+owned by a processor based on the total count of neighbors in the
+neighbor list owned by that processor.  The motivation is that more
+neighbors means a higher computational cost.  The style does not use
+neighbors per atom to assign a unique weight to each atom, because
+that value can vary depending on how the neighbor list is built.
+
+The {factor} setting is applied as an overall scale factor to the
+{neigh} weights which allows adjustment of their impact on the
+balancing operation.  The specified {factor} value must be positive.
+A value > 1.0 will increase the weights so that the ratio of max
+weight to min weight increases by {factor}.  A value < 1.0 will
+decrease the weights so that the ratio of max weight to min weight
+decreases by {factor}.  In both cases the intermediate weight values
+increase/decrease proportionally as well.  A value = 1.0 has no effect
+on the {neigh} weights.  As a rule of thumb, we have found a {factor}
+of about 0.8 often results in the best performance, since the number
+of neighbors is likely to overestimate the ideal weight.

 This weight style is useful for systems where there are different
 cutoffs used for different pairs of interations, or the density
@ -368,35 +381,48 @@ weights are computed.  Inserting a "run 0 post no"_run.html command
 before issuing the {balance} command, may be a workaround for this
 case, as it will induce the neighbor list to be built.

-The {time} weight style uses "timer data"_timer.html to estimate a
-weight for each particle.  It uses the same information as is used for
-the "MPI task timing breakdown"_Section_start.html#start_8, namely,
-the timings for sections {Pair}, {Bond}, {Kspace}, and {Neigh}.  The
-time spent in these sections of the timestep are measured for each MPI
-rank, summed up, then converted into a cost for each MPI rank relative
-to the average cost over all MPI ranks for the same sections.  That
-cost then evenly distributed over all the particles owned by that
-rank.  Finally, the {factor} setting is then appied as an overall
-scale factor to all the {time} weights as a way to fine tune the
-impact of this weight style.  Good {factor} values to use are
-typically between 0.5 and 1.2.
+The {time} weight style uses "timer data"_timer.html to estimate
+weights.  It assigns the same weight to each particle owned by a
+processor based on the total computational time spent by that
+processor.  See details below on what time window is used.  It uses
+the same timing information as is used for the "MPI task timing
+breakdown"_Section_start.html#start_8, namely, for sections {Pair},
+{Bond}, {Kspace}, and {Neigh}.  The time spent in those portions of
+the timestep are measured for each MPI rank, summed, then divided by
+the number of particles owned by that processor.  I.e. the weight is
+an effective CPU time/particle averaged over the particles on that
+processor.

-For the {balance} command the timing data is taken from the preceding
-run command, i.e. the timings are for the entire previous run.  For
-the {fix balance} command the timing data is for only the timesteps
-since the last balancing operation was performed.  If timing
-information for the required sections is not available, e.g. at the
-beginning of a run, or when the "timer"_timer.html command is set to
-either {loop} or {off}, a warning is issued.  In this case no weights
-are computed.
+The {factor} setting is applied as an overall scale factor to the
+{time} weights which allows adjustment of their impact on the
+balancing operation.  The specified {factor} value must be positive.
+A value > 1.0 will increase the weights so that the ratio of max
+weight to min weight increases by {factor}.  A value < 1.0 will
+decrease the weights so that the ratio of max weight to min weight
+decreases by {factor}.  In both cases the intermediate weight values
+increase/decrease proportionally as well.  A value = 1.0 has no effect
+on the {time} weights.  As a rule of thumb, effective values to use
+are typicall between 0.5 and 1.2.  Note that the timer quantities
+mentioned above can be affected by communication which occurs in the
+middle of the operations, e.g. pair styles with intermediate exchange
+of data witin the force computation, and likewise for KSpace solves.

-This weight style is the most generic one, and should be tried first,
-if neither the {group} or {neigh} styles are easily applicable.
-However, since the computed cost function is averaged over all local
-particles this weight style may not be highly accurate.  This style
-can also be effective as a secondary weight in combination with either
-{group} or {neigh} to offset some of inaccuracies in either of those
-heuristics.
+When using the {time} weight style with the {balance} command, the
+timing data is taken from the preceding run command, i.e. the timings
+are for the entire previous run.  For the {fix balance} command the
+timing data is for only the timesteps since the last balancing
+operation was performed.  If timing information for the required
+sections is not available, e.g. at the beginning of a run, or when the
+"timer"_timer.html command is set to either {loop} or {off}, a warning
+is issued.  In this case no weights are computed.
+
+NOTE: The {time} weight style is the most generic option, and should
+be tried first, unless the {group} style is easily applicable.
+However, since the computed cost function is averaged over all
+particles on a processor, the weights may not be highly accurate.
+This style can also be effective as a secondary weight in combination
+with either {group} or {neigh} to offset some of inaccuracies in
+either of those heuristics.

 The {var} weight style assigns per-particle weights by evaluating an
 "atom-style variable"_variable.html specified by {name}.  This is
--- a/doc/src/dump_modify.txt
+++ b/doc/src/dump_modify.txt
@ -49,8 +49,8 @@ keyword = {append} or {buffer} or {element} or {every} or {fileper} or {first} o
     -N = sort per-atom lines in descending order by the Nth column
  {thresh} args = attribute operation value
    attribute = same attributes (x,fy,etotal,sxx,etc) used by dump custom style
-    operation = "<" or "<=" or ">" or ">=" or "==" or "!="
-    value = numeric value to compare to
+    operation = "<" or "<=" or ">" or ">=" or "==" or "!=" or "|^"
+    value = numeric value to compare to, or LAST
    these 3 args can be replaced by the word "none" to turn off thresholding
  {unwrap} arg = {yes} or {no} :pre
 these keywords apply only to the {image} and {movie} "styles"_dump_image.html :l
@ -458,16 +458,59 @@ as well as memory, versus unsorted output.

 The {thresh} keyword only applies to the dump {custom}, {cfg},
 {image}, and {movie} styles.  Multiple thresholds can be specified.
-Specifying "none" turns off all threshold criteria.  If thresholds are
+Specifying {none} turns off all threshold criteria.  If thresholds are
 specified, only atoms whose attributes meet all the threshold criteria
 are written to the dump file or included in the image.  The possible
 attributes that can be tested for are the same as those that can be
 specified in the "dump custom"_dump.html command, with the exception
 of the {element} attribute, since it is not a numeric value.  Note
-that different attributes can be output by the dump custom command
-than are used as threshold criteria by the dump_modify command.
-E.g. you can output the coordinates and stress of atoms whose energy
-is above some threshold.
+that a different attributes can be used than those output by the "dump
+custom"_dump.html command.  E.g. you can output the coordinates and
+stress of atoms whose energy is above some threshold.
+
+If an atom-style variable is used as the attribute, then it can
+produce continuous numeric values or effective Boolean 0/1 values
+which may be useful for the comparision operation.  Boolean values can
+be generated by variable formulas that use comparison or Boolean math
+operators or special functions like gmask() and rmask() and grmask().
+See the "variable"_variable.html command doc page for details.
+
+NOTE: The LAST option, discussed below, is not yet implemented.  It
+will be soon.
+
+The specified value must be a simple numeric value or the word LAST.
+If LAST is used, it refers to the value of the attribute the last time
+the dump command was invoked to produce a snapshot.  This is a way to
+only dump atoms whose attribute has changed (or not changed).
+Three examples follow.
+
+dump_modify ... thresh ix != LAST :pre
+
+This will dump atoms which have crossed the periodic x boundary of the
+simulation box since the last dump.  (Note that atoms that crossed
+once and then crossed back between the two dump timesteps would not be
+included.)
+
+region foo sphere 10 20 10 15 
+variable inregion atom rmask(foo)
+dump_modify ... thresh v_inregion |^ LAST
+
+This will dump atoms which crossed the boundary of the spherical
+region since the last dump.
+
+variable charge atom "(q > 0.5) || (q < -0.5)"
+dump_modify ... thresh v_charge |^ LAST
+
+This will dump atoms whose charge has changed from an absolute value
+less than 1/2 to greater than 1/2 (or vice versa) since the last dump.
+E.g. due to reactions and subsequent charge equilibration in a
+reactive force field.
+
+The choice of operations are the usual comparison operators.  The XOR
+operation (exclusive or) is also included as "|^".  In this context,
+XOR means that if either the attribute or value is 0.0 and the other
+is non-zero, then the result is "true" and the threshold criterion is
+met.  Otherwise it is not met.

 :line

--- a/doc/src/if.txt
+++ b/doc/src/if.txt
@ -139,7 +139,7 @@ InP, myString, a123, ab_23_cd, etc :pre

 and Boolean operators:

-A == B, A != B, A < B, A <= B, A > B, A >= B, A && B, A || B, !A :pre
+A == B, A != B, A < B, A <= B, A > B, A >= B, A && B, A || B, A |^ B, !A :pre

 Each A and B is a number or string or a variable reference like $a or
 $\{abc\}, or A or B can be another Boolean expression.
@ -155,9 +155,10 @@ precedence: the unary logical NOT operator "!" has the highest
 precedence, the 4 relational operators "<", "<=", ">", and ">=" are
 next; the two remaining relational operators "==" and "!=" are next;
 then the logical AND operator "&&"; and finally the logical OR
-operator "||" has the lowest precedence.  Parenthesis can be used to
-group one or more portions of an expression and/or enforce a different
-order of evaluation than what would occur with the default precedence.
+operator "||" and logical XOR (exclusive or) operator "|^" have the
+lowest precedence.  Parenthesis can be used to group one or more
+portions of an expression and/or enforce a different order of
+evaluation than what would occur with the default precedence.

 When the 6 relational operators (first 6 in list above) compare 2
 numbers, they return either a 1.0 or 0.0 depending on whether the
@ -171,9 +172,11 @@ relationship between A and B is TRUE or FALSE (or just A).  The
 logical AND operator will return 1.0 if both its arguments are
 non-zero, else it returns 0.0.  The logical OR operator will return
 1.0 if either of its arguments is non-zero, else it returns 0.0.  The
-logical NOT operator returns 1.0 if its argument is 0.0, else it
-returns 0.0.  The 3 logical operators can only be used to operate on
-numbers, not on strings.
+logical XOR operator will return 1.0 if one of its arguments is zero
+and the other non-zero, else it returns 0.0.  The logical NOT operator
+returns 1.0 if its argument is 0.0, else it returns 0.0.  The 3
+logical operators can only be used to operate on numbers, not on
+strings.

 The overall Boolean expression produces a TRUE result if the result is
 non-zero.  If the result is zero, the expression result is FALSE.
--- a/doc/src/variable.txt
+++ b/doc/src/variable.txt
@ -47,7 +47,7 @@ style = {delete} or {index} or {loop} or {world} or {universe} or {uloop} or {st
    constants = PI, version, on, off, true, false, yes, no
    thermo keywords = vol, ke, press, etc from "thermo_style"_thermo_style.html
    math operators = (), -x, x+y, x-y, x*y, x/y, x^y, x%y,
-                     x == y, x != y, x < y, x <= y, x > y, x >= y, x && y, x || y, !x
+                     x == y, x != y, x < y, x <= y, x > y, x >= y, x && y, x || y, x |^ y, !x
    math functions = sqrt(x), exp(x), ln(x), log(x), abs(x),
                     sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x),
                     random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x)
@ -450,7 +450,7 @@ Number: 0.2, 100, 1.0e20, -15.4, etc
 Constant: PI, version, on, off, true, false, yes, no
 Thermo keywords: vol, pe, ebond, etc
 Math operators: (), -x, x+y, x-y, x*y, x/y, x^y, x%y, \
-     x == y, x != y, x < y, x <= y, x > y, x >= y, x && y, x || y, !x
+     x == y, x != y, x < y, x <= y, x > y, x >= y, x && y, x || y, x |^ y, !x
 Math functions: sqrt(x), exp(x), ln(x), log(x), abs(x), \
     sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x), \
     random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x), \
@ -551,9 +551,10 @@ division and the modulo operator "%" are next; addition and
 subtraction are next; the 4 relational operators "<", "<=", ">", and
 ">=" are next; the two remaining relational operators "==" and "!="
 are next; then the logical AND operator "&&"; and finally the logical
-OR operator "||" has the lowest precedence.  Parenthesis can be used
-to group one or more portions of a formula and/or enforce a different
-order of evaluation than what would occur with the default precedence.
+OR operator "||" and logical XOR (exclusive or) operator "|^" have the
+lowest precedence.  Parenthesis can be used to group one or more
+portions of a formula and/or enforce a different order of evaluation
+than what would occur with the default precedence.

 NOTE: Because a unary minus is higher precedence than exponentiation,
 the formula "-2^2" will evaluate to 4, not -4.  This convention is
@ -568,8 +569,10 @@ return 1.0 for all atoms whose x-coordinate is less than 10.0, and 0.0
 for the others.  The logical AND operator will return 1.0 if both its
 arguments are non-zero, else it returns 0.0.  The logical OR operator
 will return 1.0 if either of its arguments is non-zero, else it
-returns 0.0.  The logical NOT operator returns 1.0 if its argument is
-0.0, else it returns 0.0.
+returns 0.0.  The logical XOR operator will return 1.0 if one of its
+arguments is zero and the other non-zero, else it returns 0.0.  The
+logical NOT operator returns 1.0 if its argument is 0.0, else it
+returns 0.0.

 These relational and logical operators can be used as a masking or
 selection operation in a formula.  For example, the number of atoms
--- a/examples/balance/in.balance.bond.fast
+++ b/examples/balance/in.balance.bond.fast
@ -43,7 +43,7 @@ fix             2 all wall/lj93 xlo 0.0 1 1 2.5 xhi $x 1 1 2.5
 fix             3 all wall/lj93 ylo 0.0 1 1 2.5 yhi $y 1 1 2.5

 comm_style      tiled
-comm_modify     cutoff 7.5
+comm_modify     cutoff 10.0   # because bonds stretch a long ways
 fix             10 all balance 50 0.9 rcb

 #compute         1 all property/atom proc
--- a/examples/balance/in.balance.clock.dynamic
+++ b/examples/balance/in.balance.clock.dynamic
@ -52,4 +52,3 @@ fix		0 all balance 50 1.0 shift x 5 1.0 &
                weight neigh 0.5 weight time 0.66 weight store WEIGHT
 run             500
 run             500
-
--- a/examples/body/in.body
+++ b/examples/body/in.body
@ -11,13 +11,19 @@ velocity	all create 1.44 87287 loop geom
 pair_style	body 5.0
 pair_coeff	* * 1.0 1.0

-neighbor	0.3 bin
+neighbor	0.5 bin
+neigh_modify    every 1 delay 0 check yes

 fix		1 all nve/body
+#fix		1 all nvt/body temp 1.44 1.44 1.0
 fix		2 all enforce2d

 #compute         1 all body/local type 1 2 3
 #dump            1 all local 100 dump.body index c_1[1] c_1[2] c_1[3] c_1[4]

-thermo          500
+#dump		2 all image 1000 image.*.jpg type type &
+#		zoom 1.6 adiam 1.5 body type 1.0 0
+#dump_modify    2 pad 5
+
+thermo          100
 run		10000
--- a/lib/linalg/Makefile.gfortran
+++ b/lib/linalg/Makefile.gfortran
--- a/lib/linalg/Makefile.mingw32-cross
+++ b/lib/linalg/Makefile.mingw32-cross
--- a/lib/linalg/Makefile.mingw64-cross
+++ b/lib/linalg/Makefile.mingw64-cross
--- a/lib/meam/Makefile.tbird
+++ b/lib/meam/Makefile.tbird
--- a/lib/meam/meam_data.F
+++ b/lib/meam/meam_data.F
--- a/lib/meam/meam_dens_final.F
+++ b/lib/meam/meam_dens_final.F
--- a/lib/meam/meam_dens_init.F
+++ b/lib/meam/meam_dens_init.F
--- a/lib/meam/meam_force.F
+++ b/lib/meam/meam_force.F
--- a/lib/meam/meam_setup_done.F
+++ b/lib/meam/meam_setup_done.F
@ -183,6 +183,16 @@ c
      real*8, external :: zbl
      real*8, external :: compute_phi

+c check for previously allocated arrays and free them
+      if(allocated(phir)) deallocate(phir)
+      if(allocated(phirar)) deallocate(phirar)
+      if(allocated(phirar1)) deallocate(phirar1)
+      if(allocated(phirar2)) deallocate(phirar2)
+      if(allocated(phirar3)) deallocate(phirar3)
+      if(allocated(phirar4)) deallocate(phirar4)
+      if(allocated(phirar5)) deallocate(phirar5)
+      if(allocated(phirar6)) deallocate(phirar6)
+
 c allocate memory for array that defines the potential
      allocate(phir(nr,(neltypes*(neltypes+1))/2))

--- a/lib/meam/meam_setup_global.F
+++ b/lib/meam/meam_setup_global.F
--- a/lib/meam/meam_setup_param.F
+++ b/lib/meam/meam_setup_param.F
--- a/lib/reax/Makefile.g95
+++ b/lib/reax/Makefile.g95
--- a/lib/reax/Makefile.gfortran
+++ b/lib/reax/Makefile.gfortran
--- a/python/lammps.py
+++ b/python/lammps.py
@ -53,6 +53,7 @@ class lammps(object):

  def __init__(self,name="",cmdargs=None,ptr=None,comm=None):
    self.comm = comm
+    self.opened = 0

    # determine module location

@ -133,21 +134,20 @@ class lammps(object):
          # self.lmp = self.lib.lammps_open_no_mpi(0,None)

    else:
-      if isinstance(ptr,lammps):
-        # magic to convert ptr to ctypes ptr
-        pythonapi.PyCObject_AsVoidPtr.restype = c_void_p
-        pythonapi.PyCObject_AsVoidPtr.argtypes = [py_object]
-        self.lmp = c_void_p(pythonapi.PyCObject_AsVoidPtr(ptr))
-      else:
-        self.lmp = None
-        raise TypeError('Unsupported type passed as "ptr"')
+      # magic to convert ptr to ctypes ptr
+      pythonapi.PyCObject_AsVoidPtr.restype = c_void_p
+      pythonapi.PyCObject_AsVoidPtr.argtypes = [py_object]
+      self.lmp = c_void_p(pythonapi.PyCObject_AsVoidPtr(ptr))

  def __del__(self):
-    if self.lmp and self.opened: self.lib.lammps_close(self.lmp)
+    if self.lmp and self.opened:
+      self.lib.lammps_close(self.lmp)
+      self.opened = 0

  def close(self):
    if self.opened: self.lib.lammps_close(self.lmp)
    self.lmp = None
+    self.opened = 0

  def version(self):
    return self.lib.lammps_version(self.lmp)
@ -507,8 +507,7 @@ class PyLammps(object):
      elif isinstance(ptr,lammps):
        self.lmp = ptr
      else:
-        self.lmp = None
-        raise TypeError('Unsupported type passed as "ptr"')
+        self.lmp = lammps(name=name,cmdargs=cmdargs,ptr=ptr,comm=comm)
    else:
      self.lmp = lammps(name=name,cmdargs=cmdargs,ptr=None,comm=comm)
    print("LAMMPS output is captured by PyLammps wrapper")
--- a/src/MANYBODY/pair_vashishta_table.cpp
+++ b/src/MANYBODY/pair_vashishta_table.cpp
--- a/src/MANYBODY/pair_vashishta_table.h
+++ b/src/MANYBODY/pair_vashishta_table.h
--- a/src/MOLECULE/fix_cmap.cpp
+++ b/src/MOLECULE/fix_cmap.cpp
@ -89,10 +89,9 @@ FixCMAP::FixCMAP(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
  memory->create(d2cmapgrid,6,CMAPDIM,CMAPDIM,"cmap:d2grid");
  memory->create(d12cmapgrid,6,CMAPDIM,CMAPDIM,"cmap:d12grid");

-  // read, initialize, broadcast cmapgrid
+  // read and setup CMAP data

-  if (me == 0) read_grid_map(arg[3]);
-  MPI_Bcast(&cmapgrid[0][0][0],6*CMAPDIM*CMAPDIM,MPI_DOUBLE,0,world);
+  read_grid_map(arg[3]);

  // perform initial allocation of atom-based arrays
  // register with Atom class
--- a/src/USER-OMP/pair_eam_alloy_omp.cpp
+++ b/src/USER-OMP/pair_eam_alloy_omp.cpp
@ -101,6 +101,7 @@ void PairEAMAlloyOMP::coeff(int narg, char **arg)
        if (i == j) atom->set_mass(i,setfl->mass[map[i]]);
        count++;
      }
+      scale[i][j] = 1.0;
    }
  }

--- a/src/USER-OMP/pair_eam_fs_omp.cpp
+++ b/src/USER-OMP/pair_eam_fs_omp.cpp
@ -101,6 +101,7 @@ void PairEAMFSOMP::coeff(int narg, char **arg)
        if (i == j) atom->set_mass(i,fs->mass[map[i]]);
        count++;
      }
+      scale[i][j] = 1.0;
    }
  }

--- a/src/angle.cpp
+++ b/src/angle.cpp
@ -40,9 +40,6 @@ Angle::Angle(LAMMPS *lmp) : Pointers(lmp)
  vatom = NULL;
  setflag = NULL;

-  datamask = ALL_MASK;
-  datamask_ext = ALL_MASK;
-
  execution_space = Host;
  datamask_read = ALL_MASK;
  datamask_modify = ALL_MASK;
--- a/src/angle.h
+++ b/src/angle.h
@ -29,10 +29,9 @@ class Angle : protected Pointers {
  double energy;                  // accumulated energies
  double virial[6];               // accumlated virial
  double *eatom,**vatom;          // accumulated per-atom energy/virial
-  unsigned int datamask;
-  unsigned int datamask_ext;

  // KOKKOS host/device flag and data masks
+
  ExecutionSpace execution_space;
  unsigned int datamask_read,datamask_modify;
  int copymode;
@ -51,9 +50,6 @@ class Angle : protected Pointers {
  virtual double single(int, int, int, int) = 0;
  virtual double memory_usage();

-  virtual unsigned int data_mask() {return datamask;}
-  virtual unsigned int data_mask_ext() {return datamask_ext;}
-
 protected:
  int suffix_flag;             // suffix compatibility flag

--- a/src/atom.cpp
+++ b/src/atom.cpp
@ -208,9 +208,6 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
  atom_style = NULL;
  avec = NULL;

-  datamask = ALL_MASK;
-  datamask_ext = ALL_MASK;
-
  avec_map = new AtomVecCreatorMap();

 #define ATOM_CLASS
--- a/src/atom.h
+++ b/src/atom.h
@ -124,11 +124,6 @@ class Atom : protected Pointers {
  char **iname,**dname;
  int nivector,ndvector;

-  // used by USER-CUDA to flag used per-atom arrays
-
-  unsigned int datamask;
-  unsigned int datamask_ext;
-
  // atom style and per-atom array existence flags
  // customize by adding new flag

--- a/src/atom_vec.h
+++ b/src/atom_vec.h
@ -156,10 +156,6 @@ E: Invalid atom_style command

 Self-explanatory.

-E: USER-CUDA package requires a cuda enabled atom_style
-
-Self-explanatory.
-
 E: KOKKOS package requires a kokkos enabled atom_style

 Self-explanatory.
--- a/src/balance.cpp
+++ b/src/balance.cpp
@ -272,7 +272,7 @@ void Balance::command(int narg, char **arg)
  // imbinit = initial imbalance

  double maxinit;
-  init_imbalance();
+  init_imbalance(0);
  set_weights();
  double imbinit = imbalance_factor(maxinit);

@ -543,12 +543,13 @@ void Balance::weight_storage(char *prefix)

 /* ----------------------------------------------------------------------
   invoke init() for each Imbalance class
+   flag = 0 for call from Balance, 1 for call from FixBalance
 ------------------------------------------------------------------------- */

-void Balance::init_imbalance()
+void Balance::init_imbalance(int flag)
 {
  if (!wtflag) return;
-  for (int n = 0; n < nimbalance; n++) imbalances[n]->init();
+  for (int n = 0; n < nimbalance; n++) imbalances[n]->init(flag);
 }

 /* ----------------------------------------------------------------------
--- a/src/balance.h
+++ b/src/balance.h
@ -38,7 +38,7 @@ class Balance : protected Pointers {
  void command(int, char **);
  void options(int, int, char **);
  void weight_storage(char *);
-  void init_imbalance();
+  void init_imbalance(int);
  void set_weights();
  double imbalance_factor(double &);
  void shift_setup(char *, int, double);
--- a/src/bond.cpp
+++ b/src/bond.cpp
@ -44,9 +44,6 @@ Bond::Bond(LAMMPS *lmp) : Pointers(lmp)
  vatom = NULL;
  setflag = NULL;

-  datamask = ALL_MASK;
-  datamask_ext = ALL_MASK;
-
  execution_space = Host;
  datamask_read = ALL_MASK;
  datamask_modify = ALL_MASK;
--- a/src/bond.h
+++ b/src/bond.h
@ -29,10 +29,9 @@ class Bond : protected Pointers {
  double energy;                  // accumulated energies
  double virial[6];               // accumlated virial
  double *eatom,**vatom;          // accumulated per-atom energy/virial
-  unsigned int datamask;
-  unsigned int datamask_ext;

  // KOKKOS host/device flag and data masks
+
  ExecutionSpace execution_space;
  unsigned int datamask_read,datamask_modify;
  int copymode;
@ -51,9 +50,6 @@ class Bond : protected Pointers {
  virtual double single(int, double, int, int, double &) = 0;
  virtual double memory_usage();

-  virtual unsigned int data_mask() {return datamask;}
-  virtual unsigned int data_mask_ext() {return datamask_ext;}
-
  void write_file(int, char**);

 protected:
--- a/src/comm_tiled.h
+++ b/src/comm_tiled.h
@ -155,10 +155,6 @@ class CommTiled : public Comm {

 /* ERROR/WARNING messages:

-E: USER-CUDA package does not yet support comm_style tiled
-
-Self-explanatory.
-
 E: KOKKOS package does not yet support comm_style tiled

 Self-explanatory.
--- a/src/compute.cpp
+++ b/src/compute.cpp
@ -99,9 +99,6 @@ Compute::Compute(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp),
  
  // data masks

-  datamask = ALL_MASK;
-  datamask_ext = ALL_MASK;
-
  execution_space = Host;
  datamask_read = ALL_MASK;
  datamask_modify = ALL_MASK;
--- a/src/compute.h
+++ b/src/compute.h
@ -84,9 +84,6 @@ class Compute : protected Pointers {
  int comm_reverse;         // size of reverse communication (0 if none)
  int dynamic_group_allow;  // 1 if can be used with dynamic group, else 0

-  unsigned int datamask;
-  unsigned int datamask_ext;
-
  // KOKKOS host/device flag and data masks

  ExecutionSpace execution_space;
@ -140,9 +137,6 @@ class Compute : protected Pointers {
                                   double, double, double,
                                   double, double, double) {}

-  virtual int unsigned data_mask() {return datamask;}
-  virtual int unsigned data_mask_ext() {return datamask_ext;}
-
 protected:
  int instance_me;             // which Compute class instantiation I am

--- a/src/dihedral.cpp
+++ b/src/dihedral.cpp
@ -41,9 +41,6 @@ Dihedral::Dihedral(LAMMPS *lmp) : Pointers(lmp)
  vatom = NULL;
  setflag = NULL;

-  datamask = ALL_MASK;
-  datamask_ext = ALL_MASK;
-
  execution_space = Host;
  datamask_read = ALL_MASK;
  datamask_modify = ALL_MASK;
--- a/src/dihedral.h
+++ b/src/dihedral.h
@ -29,10 +29,9 @@ class Dihedral : protected Pointers {
  double energy;                     // accumulated energy
  double virial[6];                  // accumlated virial
  double *eatom,**vatom;             // accumulated per-atom energy/virial
-  unsigned int datamask;
-  unsigned int datamask_ext;

  // KOKKOS host/device flag and data masks
+
  ExecutionSpace execution_space;
  unsigned int datamask_read,datamask_modify;
  int copymode;
@ -49,9 +48,6 @@ class Dihedral : protected Pointers {
  virtual void write_data(FILE *) {}
  virtual double memory_usage();

-  virtual unsigned int data_mask() {return datamask;}
-  virtual unsigned int data_mask_ext() {return datamask_ext;}
-
 protected:
  int suffix_flag;             // suffix compatibility flag

--- a/src/dump_custom.cpp
+++ b/src/dump_custom.cpp
@ -43,7 +43,7 @@ enum{ID,MOL,PROC,PROCP1,TYPE,ELEMENT,MASS,
     OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ,
     TQX,TQY,TQZ,
     COMPUTE,FIX,VARIABLE,INAME,DNAME};
-enum{LT,LE,GT,GE,EQ,NEQ};
+enum{LT,LE,GT,GE,EQ,NEQ,XOR};
 enum{INT,DOUBLE,STRING,BIGINT};    // same as in DumpCFG

 #define INVOKED_PERATOM 8
@ -947,6 +947,11 @@ int DumpCustom::count()
      } else if (thresh_op[ithresh] == NEQ) {
        for (i = 0; i < nlocal; i++, ptr += nstride)
          if (choose[i] && *ptr == value) choose[i] = 0;
+      } else if (thresh_op[ithresh] == XOR) {
+        for (i = 0; i < nlocal; i++, ptr += nstride)
+          if (choose[i] && (*ptr == 0.0 && value == 0.0) || 
+              (*ptr != 0.0 && value != 0.0))
+            choose[i] = 0;
      }
    }
  }
@ -1835,6 +1840,7 @@ int DumpCustom::modify_param(int narg, char **arg)
    else if (strcmp(arg[2],">=") == 0) thresh_op[nthresh] = GE;
    else if (strcmp(arg[2],"==") == 0) thresh_op[nthresh] = EQ;
    else if (strcmp(arg[2],"!=") == 0) thresh_op[nthresh] = NEQ;
+    else if (strcmp(arg[2],"|^") == 0) thresh_op[nthresh] = XOR;
    else error->all(FLERR,"Invalid dump_modify threshold operator");

    // set threshold value
--- a/src/fix.cpp
+++ b/src/fix.cpp
@ -95,10 +95,7 @@ id(NULL), style(NULL), eatom(NULL), vatom(NULL)
  maxeatom = maxvatom = 0;
  vflag_atom = 0;

-  // CUDA and KOKKOS per-fix data masks
-
-  datamask = ALL_MASK;
-  datamask_ext = ALL_MASK;
+  // KOKKOS per-fix data masks

  execution_space = Host;
  datamask_read = ALL_MASK;
--- a/src/fix.h
+++ b/src/fix.h
@ -99,11 +99,6 @@ class Fix : protected Pointers {
  ExecutionSpace execution_space;
  unsigned int datamask_read,datamask_modify;

-  // USER-CUDA per-fix data masks
-
-  unsigned int datamask;
-  unsigned int datamask_ext;
-
  Fix(class LAMMPS *, int, char **);
  virtual ~Fix();
  void modify_params(int, char **);
@ -211,9 +206,6 @@ class Fix : protected Pointers {

  virtual double memory_usage() {return 0.0;}

-  virtual unsigned int data_mask() {return datamask;}
-  virtual unsigned int data_mask_ext() {return datamask_ext;}
-
 protected:
  int instance_me;        // which Fix class instantiation I am

--- a/src/fix_balance.cpp
+++ b/src/fix_balance.cpp
@ -26,6 +26,7 @@
 #include "modify.h"
 #include "fix_store.h"
 #include "rcb.h"
+#include "timer.h"
 #include "error.h"

 using namespace LAMMPS_NS;
@ -113,7 +114,8 @@ FixBalance::FixBalance(LAMMPS *lmp, int narg, char **arg) :
  // only force reneighboring if nevery > 0

  if (nevery) force_reneighbor = 1;
-
+  lastbalance = -1;
+  
  // compute initial outputs

  itercount = 0;
@ -153,7 +155,7 @@ void FixBalance::init()
  if (force->kspace) kspace_flag = 1;
  else kspace_flag = 0;

-  balance->init_imbalance();
+  balance->init_imbalance(1);
 }

 /* ---------------------------------------------------------------------- */
@ -170,6 +172,12 @@ void FixBalance::setup(int vflag)

 void FixBalance::setup_pre_exchange()
 {
+  // do not allow rebalancing twice on same timestep
+  // even if wanted to, can mess up elapsed time in ImbalanceTime
+
+  if (update->ntimestep == lastbalance) return;
+  lastbalance = update->ntimestep;
+
  // insure atoms are in current box & update box via shrink-wrap
  // has to be be done before rebalance() invokes Irregular::migrate_atoms()
  //   since it requires atoms be inside simulation box
@ -202,6 +210,12 @@ void FixBalance::pre_exchange()

  if (nevery && update->ntimestep < next_reneighbor) return;

+  // do not allow rebalancing twice on same timestep
+  // even if wanted to, can mess up elapsed time in ImbalanceTime
+
+  if (update->ntimestep == lastbalance) return;
+  lastbalance = update->ntimestep;
+
  // insure atoms are in current box & update box via shrink-wrap
  // no exchange() since doesn't matter if atoms are assigned to correct procs

@ -284,7 +298,7 @@ void FixBalance::rebalance()
  if (kspace_flag) force->kspace->setup_grid();

  // pending triggers pre_neighbor() to compute final imbalance factor
-  // can only be done after atoms migrate in caller's comm->exchange()
+  // can only be done after atoms migrate in comm->exchange()

  pending = 1;
 }
--- a/src/fix_balance.h
+++ b/src/fix_balance.h
@ -53,6 +53,7 @@ class FixBalance : public Fix {
  int itercount;                // iteration count of last call to Balance
  int kspace_flag;              // 1 if KSpace solver defined
  int pending;
+  bigint lastbalance;           // last timestep balancing was attempted

  class Balance *balance;
  class Irregular *irregular;
--- a/src/fix_tmd.cpp
+++ b/src/fix_tmd.cpp
@ -460,9 +460,9 @@ void FixTMD::readfile(char *file)
      else
        n = sscanf(bufptr,TAGINT_FORMAT " %lg %lg %lg",&itag,&x,&y,&z);

-      if (n < 0 && me == 0) {
-        error->warning(FLERR,"Ignoring empty or incorrectly formatted"
-                             " line in target file");
+      if (n < 0) {
+        if (me == 0) error->warning(FLERR,"Ignoring empty or incorrectly"
+                                    " formatted line in target file");
        bufptr = next + 1;
        continue;
      }
--- a/src/imbalance.h
+++ b/src/imbalance.h
@ -25,13 +25,13 @@ class Imbalance : protected Pointers {
  virtual ~Imbalance() {};

  // parse options. return number of arguments consumed (required)
-  virtual int options(int narg, char **arg) = 0;
+  virtual int options(int, char **) = 0;
  // reinitialize internal data (needed for fix balance) (optional)
-  virtual void init() {};
+  virtual void init(int) {};
  // compute and apply weight factors to local atom array (required)
-  virtual void compute(double *weights) = 0;
+  virtual void compute(double *) = 0;
  // print information about the state of this imbalance compute (required)
-  virtual void info(FILE *fp) = 0;
+  virtual void info(FILE *) = 0;

  // disallow default and copy constructor, assignment operator
  // private:
--- a/src/imbalance_group.cpp
+++ b/src/imbalance_group.cpp
@ -49,6 +49,7 @@ int ImbalanceGroup::options(int narg, char **arg)
    if (id[i] < 0)
      error->all(FLERR,"Unknown group in balance weight command");
    factor[i] = force->numeric(FLERR,arg[2*i+2]);
+    if (factor[i] <= 0.0) error->all(FLERR,"Illegal balance weight command");
  }
  return 2*num+1;
 }
@ -65,12 +66,10 @@ void ImbalanceGroup::compute(double *weight)

  for (int i = 0; i < nlocal; ++i) {
    const int imask = mask[i];
-    double iweight = weight[i];
    for (int j = 0; j < num; ++j) {
      if (imask & bitmask[id[j]])
-        iweight *= factor[j];
+        weight[i] *= factor[j];
    }
-    weight[i] = iweight;
  }
 }

--- a/src/imbalance_neigh.cpp
+++ b/src/imbalance_neigh.cpp
@ -23,6 +23,8 @@

 using namespace LAMMPS_NS;

+#define BIG 1.0e20
+
 /* -------------------------------------------------------------------- */

 ImbalanceNeigh::ImbalanceNeigh(LAMMPS *lmp) : Imbalance(lmp)
@ -36,7 +38,7 @@ int ImbalanceNeigh::options(int narg, char **arg)
 {
  if (narg < 1) error->all(FLERR,"Illegal balance weight command");
  factor = force->numeric(FLERR,arg[0]);
-  if (factor < 0.0) error->all(FLERR,"Illegal balance weight command");
+  if (factor <= 0.0) error->all(FLERR,"Illegal balance weight command");
  return 1;
 }

@ -49,7 +51,7 @@ void ImbalanceNeigh::compute(double *weight)
  if (factor == 0.0) return;

  // find suitable neighbor list
-  // we can only make use of certain (conventional) neighbor lists
+  // can only use certain conventional neighbor lists

  for (req = 0; req < neighbor->old_nrequest; ++req) {
    if ((neighbor->old_requests[req]->half ||
@ -62,36 +64,46 @@ void ImbalanceNeigh::compute(double *weight)

  if (req >= neighbor->old_nrequest || neighbor->ago < 0) {
    if (comm->me == 0 && !did_warn)
-      error->warning(FLERR,"No suitable neighbor list found. "
-                     "Neighbor weighted balancing skipped");
+      error->warning(FLERR,"Balance weight neigh skipped b/c no list found");
    did_warn = 1;
    return;
  }

+  // neighsum = total neigh count for atoms on this proc
+  // localwt = weight assigned to each owned atom
+
  NeighList *list = neighbor->lists[req];
-  bigint neighsum = 0;
-  
  const int inum = list->inum;
  const int * const ilist = list->ilist;
  const int * const numneigh = list->numneigh;
+  int nlocal = atom->nlocal;

-  // first pass: get local number of neighbors
-
+  bigint neighsum = 0;
  for (int i = 0; i < inum; ++i) neighsum += numneigh[ilist[i]];
+  double localwt = 0.0;
+  if (nlocal) localwt = 1.0*neighsum/nlocal;

-  double allatoms = static_cast <double>(atom->natoms);
-  if (allatoms == 0.0) allatoms = 1.0;
-  double allavg;
-  double myavg = static_cast<double>(neighsum)/allatoms;
-  MPI_Allreduce(&myavg,&allavg,1,MPI_DOUBLE,MPI_SUM,world);
-  
-  // second pass: compute and apply weights
+  if (nlocal && localwt <= 0.0) error->one(FLERR,"Balance weight <= 0.0");

-  double scale = 1.0/allavg;
-  for (int ii = 0; ii < inum; ++ii) {
-    const int i = ilist[ii];
-    weight[i] *= (1.0-factor) + factor*scale*numneigh[i];
+  // apply factor if specified != 1.0
+  // wtlo,wthi = lo/hi values excluding 0.0 due to no atoms on this proc
+  // lo value does not change
+  // newhi = new hi value to give hi/lo ratio factor times larger/smaller
+  // expand/contract all localwt values from lo->hi to lo->newhi
+
+  if (factor != 1.0) {
+    double wtlo,wthi;
+    if (localwt == 0.0) localwt = BIG;
+    MPI_Allreduce(&localwt,&wtlo,1,MPI_DOUBLE,MPI_MIN,world);
+    if (localwt == BIG) localwt = 0.0;
+    MPI_Allreduce(&localwt,&wthi,1,MPI_DOUBLE,MPI_MAX,world);
+    if (wtlo == wthi) return;
+
+    double newhi = wthi*factor;
+    localwt = wtlo + ((localwt-wtlo)/(wthi-wtlo)) * (newhi-wtlo);
  }
+
+  for (int i = 0; i < nlocal; i++) weight[i] *= localwt;
 }

 /* -------------------------------------------------------------------- */
--- a/src/imbalance_time.cpp
+++ b/src/imbalance_time.cpp
@ -21,6 +21,8 @@

 using namespace LAMMPS_NS;

+#define BIG 1.0e20
+
 /* -------------------------------------------------------------------- */

 ImbalanceTime::ImbalanceTime(LAMMPS *lmp) : Imbalance(lmp) {}
@ -31,51 +33,74 @@ int ImbalanceTime::options(int narg, char **arg)
 {
  if (narg < 1) error->all(FLERR,"Illegal balance weight command");
  factor = force->numeric(FLERR,arg[0]);
-  if (factor < 0.0) error->all(FLERR,"Illegal balance weight command");
+  if (factor <= 0.0) error->all(FLERR,"Illegal balance weight command");
  return 1;
 }

 /* ----------------------------------------------------------------------
-   reset last, needed for fix balance caller
+   reset last and timers if necessary
 ------------------------------------------------------------------------- */

-void ImbalanceTime::init()
+void ImbalanceTime::init(int flag)
 {
  last = 0.0;
+
+  // flag = 1 if called from FixBalance at start of run
+  //   init Timer, so accumulated time not carried over from previous run
+  // should NOT init Timer if called from Balance, it uses time from last run
+
+  if (flag) timer->init();
 }

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

 void ImbalanceTime::compute(double *weight)
 {
-  const int nlocal = atom->nlocal;
-  const bigint natoms = atom->natoms;
+  if (!timer->has_normal()) return;

-  if (factor == 0.0) return;
+  // cost = CPU time for relevant timers since last invocation
+  // localwt = weight assigned to each owned atom
+  // just return if no time yet tallied

-  // compute the cost function of based on relevant timers
-  
-  if (timer->has_normal()) {
-      double cost = -last;
-      cost += timer->get_wall(Timer::PAIR);
-      cost += timer->get_wall(Timer::NEIGH);
-      cost += timer->get_wall(Timer::BOND);
-      cost += timer->get_wall(Timer::KSPACE);
+  double cost = -last;
+  cost += timer->get_wall(Timer::PAIR);
+  cost += timer->get_wall(Timer::NEIGH);
+  cost += timer->get_wall(Timer::BOND);
+  cost += timer->get_wall(Timer::KSPACE);

-      double allcost;
-      MPI_Allreduce(&cost,&allcost,1,MPI_DOUBLE,MPI_SUM,world);
+  double maxcost;
+  MPI_Allreduce(&cost,&maxcost,1,MPI_DOUBLE,MPI_MAX,world);
+  if (maxcost <= 0.0) return;

-      if ((allcost > 0.0) && (nlocal > 0)) {
-        const double avgcost = allcost/natoms;
-        const double localcost = cost/nlocal;
-        const double scale = (1.0-factor) + factor*localcost/avgcost;
-        for (int i = 0; i < nlocal; ++i) weight[i] *= scale;
-      }
+  int nlocal = atom->nlocal;
+  double localwt = 0.0;
+  if (nlocal) localwt = cost/nlocal;

-      // record time up to this point
+  if (nlocal && localwt <= 0.0) error->one(FLERR,"Balance weight <= 0.0");

-      last += cost;
+  // apply factor if specified != 1.0
+  // wtlo,wthi = lo/hi values excluding 0.0 due to no atoms on this proc
+  // lo value does not change
+  // newhi = new hi value to give hi/lo ratio factor times larger/smaller
+  // expand/contract all localwt values from lo->hi to lo->newhi
+
+  if (factor != 1.0) {
+    double wtlo,wthi;
+    if (localwt == 0.0) localwt = BIG;
+    MPI_Allreduce(&localwt,&wtlo,1,MPI_DOUBLE,MPI_MIN,world);
+    if (localwt == BIG) localwt = 0.0;
+    MPI_Allreduce(&localwt,&wthi,1,MPI_DOUBLE,MPI_MAX,world);
+    if (wtlo == wthi) return;
+
+    double newhi = wthi*factor;
+    localwt = wtlo + ((localwt-wtlo)/(wthi-wtlo)) * (newhi-wtlo);
  }
+
+  for (int i = 0; i < nlocal; i++) weight[i] *= localwt;
+  
+  // record time up to this point
+  
+  last += cost;
 }

 /* -------------------------------------------------------------------- */
--- a/src/imbalance_time.h
+++ b/src/imbalance_time.h
@ -27,7 +27,7 @@ class ImbalanceTime : public Imbalance {
  // parse options, return number of arguments consumed
  virtual int options(int, char **);
  // reinitialize internal data
-  virtual void init();
+  virtual void init(int);
  // compute and apply weight factors to local atom array
  virtual void compute(double *);
  // print information about the state of this imbalance compute
--- a/src/imbalance_var.cpp
+++ b/src/imbalance_var.cpp
@ -45,14 +45,14 @@ int ImbalanceVar::options(int narg, char **arg)
  int len = strlen(arg[0]) + 1;
  name = new char[len];
  memcpy(name,arg[0],len);
-  init();
+  init(0);

  return 1;
 }

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

-void ImbalanceVar::init()
+void ImbalanceVar::init(int flag)
 {
  id = input->variable->find(name);
  if (id < 0) {
@ -75,7 +75,15 @@ void ImbalanceVar::compute(double *weight)
  memory->create(values,nlocal,"imbalance:values");

  input->variable->compute_atom(id,all,values,1,0);
-  for (int i = 0; i < nlocal; ++i) weight[i] *= values[i];
+
+  int flag = 0;
+  for (int i = 0; i < nlocal; i++)
+    if (values[i] <= 0.0) flag = 1;
+  int flagall;
+  MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
+  if (flagall) error->one(FLERR,"Balance weight <= 0.0");
+
+  for (int i = 0; i < nlocal; i++) weight[i] *= values[i];

  memory->destroy(values);
 }
--- a/src/imbalance_var.h
+++ b/src/imbalance_var.h
@ -27,7 +27,7 @@ class ImbalanceVar : public Imbalance {
  // parse options. return number of arguments consumed.
  virtual int options(int, char **);
  // re-initialize internal data, e.g. variable ID
-  virtual void init();
+  virtual void init(int);
  // compute per-atom imbalance and apply to weight array
  virtual void compute(double *);
  // print information about the state of this imbalance compute (required)
--- a/src/improper.cpp
+++ b/src/improper.cpp
@ -38,9 +38,6 @@ Improper::Improper(LAMMPS *lmp) : Pointers(lmp)
  vatom = NULL;
  setflag = NULL;

-  datamask = ALL_MASK;
-  datamask_ext = ALL_MASK;
-
  execution_space = Host;
  datamask_read = ALL_MASK;
  datamask_modify = ALL_MASK;
--- a/src/improper.h
+++ b/src/improper.h
@ -29,10 +29,9 @@ class Improper : protected Pointers {
  double energy;                  // accumulated energies
  double virial[6];               // accumlated virial
  double *eatom,**vatom;          // accumulated per-atom energy/virial
-  unsigned int datamask;
-  unsigned int datamask_ext;

  // KOKKOS host/device flag and data masks
+
  ExecutionSpace execution_space;
  unsigned int datamask_read,datamask_modify;
  int copymode;
@ -49,9 +48,6 @@ class Improper : protected Pointers {
  virtual void write_data(FILE *) {}
  virtual double memory_usage();

-  virtual unsigned int data_mask() {return datamask;}
-  virtual unsigned int data_mask_ext() {return datamask_ext;}
-
 protected:
  int suffix_flag;             // suffix compatibility flag

--- a/src/input.h
+++ b/src/input.h
@ -328,12 +328,6 @@ E: Package command after simulation box is defined
 The package command cannot be used afer a read_data, read_restart, or
 create_box command.

-E: Package cuda command without USER-CUDA package enabled
-
-The USER-CUDA package must be installed via "make yes-user-cuda"
-before LAMMPS is built, and the "-c on" must be used to enable the
-package.
-
 E: Package gpu command without GPU package installed

 The GPU package must be installed via "make yes-gpu" before LAMMPS is
--- a/src/kspace.cpp
+++ b/src/kspace.cpp
@ -88,9 +88,6 @@ KSpace::KSpace(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
  eatom = NULL;
  vatom = NULL;

-  datamask = ALL_MASK;
-  datamask_ext = ALL_MASK;
-
  execution_space = Host;
  datamask_read = ALL_MASK;
  datamask_modify = ALL_MASK;
--- a/src/kspace.h
+++ b/src/kspace.h
@ -80,10 +80,8 @@ class KSpace : protected Pointers {

  int group_group_enable;         // 1 if style supports group/group calculation

-  unsigned int datamask;
-  unsigned int datamask_ext;
-
  // KOKKOS host/device flag and data masks
+
  ExecutionSpace execution_space;
  unsigned int datamask_read,datamask_modify;
  int copymode;
--- a/src/lammps.h
+++ b/src/lammps.h
@ -168,19 +168,10 @@ E: Cannot use -cuda on and -kokkos on together

 This is not allowed since both packages can use GPUs.

-E: Cannot use -cuda on without USER-CUDA installed
-
-The USER-CUDA package must be installed via "make yes-user-cuda"
-before LAMMPS is built.
-
 E: Cannot use -kokkos on without KOKKOS installed

 Self-explanatory.

-E: Using suffix cuda without USER-CUDA package enabled
-
-Self-explanatory.
-
 E: Using suffix gpu without GPU package installed

 Self-explanatory.
--- a/src/pair.cpp
+++ b/src/pair.cpp
@ -100,9 +100,6 @@ Pair::Pair(LAMMPS *lmp) : Pointers(lmp)

  // KOKKOS per-fix data masks

-  datamask = ALL_MASK;
-  datamask_ext = ALL_MASK;
-
  execution_space = Host;
  datamask_read = ALL_MASK;
  datamask_modify = ALL_MASK;
--- a/src/pair.h
+++ b/src/pair.h
@ -97,9 +97,6 @@ class Pair : protected Pointers {
  class NeighList *listmiddle;
  class NeighList *listouter;

-  unsigned int datamask;
-  unsigned int datamask_ext;
-
  int allocated;                 // 0/1 = whether arrays are allocated
                                 //       public so external driver can check
  int compute_flag;              // 0 if skip compute()
@ -191,9 +188,6 @@ class Pair : protected Pointers {
  virtual void min_xf_get(int) {}
  virtual void min_x_set(int) {}

-  virtual unsigned int data_mask() {return datamask;}
-  virtual unsigned int data_mask_ext() {return datamask_ext;}
-
  // management of callbacks to be run from ev_tally()

 protected:
--- a/src/suffix.h
+++ b/src/suffix.h
@ -20,9 +20,8 @@ namespace Suffix {
  static const int NONE = 0;
  static const int OPT  = 1<<0;
  static const int GPU  = 1<<1;
-  static const int CUDA = 1<<2;
-  static const int OMP  = 1<<3;
-  static const int INTEL  = 1<<4;
+  static const int OMP  = 1<<2;
+  static const int INTEL  = 1<<3;
 }

 }
--- a/src/update.h
+++ b/src/update.h
@ -81,14 +81,6 @@ class Update : protected Pointers {

 /* ERROR/WARNING messages:

-E: USER-CUDA mode requires CUDA variant of run style
-
-CUDA mode is enabled, so the run style must include a cuda suffix.
-
-E: USER-CUDA mode requires CUDA variant of min style
-
-CUDA mode is enabled, so the min style must include a cuda suffix.
-
 E: Illegal ... command

 Self-explanatory.  Check the input script syntax and compare to the
--- a/src/variable.cpp
+++ b/src/variable.cpp
@ -56,11 +56,11 @@ enum{INDEX,LOOP,WORLD,UNIVERSE,ULOOP,STRING,GETENV,
 enum{ARG,OP};

 // customize by adding a function
-// if add before OR,
+// if add before XOR:
 // also set precedence level in constructor and precedence length in *.h

 enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,MODULO,UNARY,
-     NOT,EQ,NE,LT,LE,GT,GE,AND,OR,
+     NOT,EQ,NE,LT,LE,GT,GE,AND,OR,XOR,
     SQRT,EXP,LN,LOG,ABS,SIN,COS,TAN,ASIN,ACOS,ATAN,ATAN2,
     RANDOM,NORMAL,CEIL,FLOOR,ROUND,RAMP,STAGGER,LOGFREQ,LOGFREQ2,
     STRIDE,STRIDE2,VDISPLACE,SWIGGLE,CWIGGLE,GMASK,RMASK,GRMASK,
@ -103,7 +103,7 @@ Variable::Variable(LAMMPS *lmp) : Pointers(lmp)
  // customize by assigning a precedence level

  precedence[DONE] = 0;
-  precedence[OR] = 1;
+  precedence[OR] = precedence[XOR] = 1;
  precedence[AND] = 2;
  precedence[EQ] = precedence[NE] = 3;
  precedence[LT] = precedence[LE] = precedence[GT] = precedence[GE] = 4;
@ -2088,9 +2088,9 @@ double Variable::evaluate(char *str, Tree **tree)
        op = AND;
        i++;
      } else if (onechar == '|') {
-        if (str[i+1] != '|')
-          error->all(FLERR,"Invalid syntax in variable formula");
-        op = OR;
+        if (str[i+1] == '|') op = OR;
+        else if (str[i+1] == '^') op = XOR;
+        else error->all(FLERR,"Invalid syntax in variable formula");
        i++;
      } else op = DONE;

@ -2180,6 +2180,10 @@ double Variable::evaluate(char *str, Tree **tree)
          } else if (opprevious == OR) {
            if (value1 != 0.0 || value2 != 0.0) argstack[nargstack++] = 1.0;
            else argstack[nargstack++] = 0.0;
+          } else if (opprevious == XOR) {
+            if ((value1 == 0.0 && value2 != 0.0) || 
+                (value1 != 0.0 && value2 == 0.0)) argstack[nargstack++] = 1.0;
+            else argstack[nargstack++] = 0.0;
          }
        }
      }
@ -2390,6 +2394,17 @@ double Variable::collapse_tree(Tree *tree)
    return tree->value;
  }

+  if (tree->type == XOR) {
+    arg1 = collapse_tree(tree->first);
+    arg2 = collapse_tree(tree->second);
+    if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
+    tree->type = VALUE;
+    if ((arg1 == 0.0 && arg2 != 0.0) || (arg1 != 0.0 && arg2 == 0.0))
+      tree->value = 1.0;
+    else tree->value = 0.0;
+    return tree->value;
+  }
+
  if (tree->type == SQRT) {
    arg1 = collapse_tree(tree->first);
    if (tree->first->type != VALUE) return 0.0;
@ -2809,6 +2824,13 @@ double Variable::eval_tree(Tree *tree, int i)
      return 1.0;
    else return 0.0;
  }
+  if (tree->type == XOR) {
+    if ((eval_tree(tree->first,i) == 0.0 && eval_tree(tree->second,i) != 0.0) 
+        ||
+        (eval_tree(tree->first,i) != 0.0 && eval_tree(tree->second,i) == 0.0))
+      return 1.0;
+    else return 0.0;
+  }

  if (tree->type == SQRT) {
    arg1 = eval_tree(tree->first,i);
@ -4678,9 +4700,9 @@ double Variable::evaluate_boolean(char *str)
        op = AND;
        i++;
      } else if (onechar == '|') {
-        if (str[i+1] != '|')
-          error->all(FLERR,"Invalid Boolean syntax in if command");
-        op = OR;
+        if (str[i+1] == '|') op = OR;
+        else if (str[i+1] == '^') op = XOR;
+        else error->all(FLERR,"Invalid Boolean syntax in if command");
        i++;
      } else op = DONE;

@ -4764,6 +4786,12 @@ double Variable::evaluate_boolean(char *str)
          if (flag2) error->all(FLERR,"Invalid Boolean syntax in if command");
          if (value1 != 0.0 || value2 != 0.0) argstack[nargstack].value = 1.0;
          else argstack[nargstack].value = 0.0;
+        } else if (opprevious == XOR) {
+          if (flag2) error->all(FLERR,"Invalid Boolean syntax in if command");
+          if ((value1 == 0.0 && value2 != 0.0) ||
+              (value1 != 0.0 && value2 == 0.0))
+            argstack[nargstack].value = 1.0;
+          else argstack[nargstack].value = 0.0;
        }

        argstack[nargstack++].flag = 0;
@ -4785,72 +4813,6 @@ double Variable::evaluate_boolean(char *str)
  return argstack[0].value;
 }

-/* ---------------------------------------------------------------------- */
-
-unsigned int Variable::data_mask(int ivar)
-{
-  if (eval_in_progress[ivar]) return EMPTY_MASK;
-  eval_in_progress[ivar] = 1;
-  unsigned int datamask = data_mask(data[ivar][0]);
-  eval_in_progress[ivar] = 0;
-  return datamask;
-}
-
-/* ---------------------------------------------------------------------- */
-
-unsigned int Variable::data_mask(char *str)
-{
-  unsigned int datamask = EMPTY_MASK;
-
-  for (unsigned int i = 0; i < strlen(str)-2; i++) {
-    int istart = i;
-    while (isalnum(str[i]) || str[i] == '_') i++;
-    int istop = i-1;
-
-    int n = istop - istart + 1;
-    char *word = new char[n+1];
-    strncpy(word,&str[istart],n);
-    word[n] = '\0';
-
-    // ----------------
-    // compute
-    // ----------------
-
-    if ((strncmp(word,"c_",2) == 0) && (i>0) && (!isalnum(str[i-1]))) {
-      if (domain->box_exist == 0)
-        error->all(FLERR,
-                   "Variable evaluation before simulation box is defined");
-
-      int icompute = modify->find_compute(word+2);
-      if (icompute < 0)
-        error->all(FLERR,"Invalid compute ID in variable formula");
-
-      datamask &= modify->compute[icompute]->data_mask();
-    }
-
-    if ((strncmp(word,"f_",2) == 0) && (i>0) && (!isalnum(str[i-1]))) {
-      if (domain->box_exist == 0)
-        error->all(FLERR,
-                   "Variable evaluation before simulation box is defined");
-
-      int ifix = modify->find_fix(word+2);
-      if (ifix < 0) error->all(FLERR,"Invalid fix ID in variable formula");
-
-      datamask &= modify->fix[ifix]->data_mask();
-    }
-
-    if ((strncmp(word,"v_",2) == 0) && (i>0) && (!isalnum(str[i-1]))) {
-      int ivar = find(word+2);
-      if (ivar < 0) error->all(FLERR,"Invalid variable name in variable formula");
-      datamask &= data_mask(ivar);
-    }
-
-    delete [] word;
-  }
-
-  return datamask;
-}
-
 /* ----------------------------------------------------------------------
   class to read variable values from a file
   for flag = SCALARFILE, reads one value per line
--- a/src/variable.h
+++ b/src/variable.h
@ -49,9 +49,6 @@ class Variable : protected Pointers {
  tagint int_between_brackets(char *&, int);
  double evaluate_boolean(char *);

-  unsigned int data_mask(int ivar);
-  unsigned int data_mask(char *str);
-
 private:
  int me;
  int nvar;                // # of defined variables
@ -78,8 +75,8 @@ class Variable : protected Pointers {
  class RanMars *randomequal;   // random number generator for equal-style vars
  class RanMars *randomatom;    // random number generator for atom-style vars

-  int precedence[17];      // precedence level of math operators
-                           // set length to include up to OR in enum
+  int precedence[18];      // precedence level of math operators
+                           // set length to include up to XOR in enum

  class Python *python;    // ptr to embedded Python interpreter

--- a/src/version.h
+++ b/src/version.h
@ -1 +1 @@
-#define LAMMPS_VERSION "29 Sep 2016"
+#define LAMMPS_VERSION "5 Oct 2016"
--- a/tools/ch2lmp/README
+++ b/tools/ch2lmp/README
@ -41,9 +41,9 @@ should be -option=value (e.g. -border=5).
 -border      add border to all sides of simulation box [default: 0 A]
 -ax          rotation around x-axis
 -ay          rotation around y-axis
-az          rotation around z-axis           
+-az          rotation around z-axis
 -cd          correction for dihedral for carbohydrate systems
-cmap        add CMAP section to data file and fix cmap command lines in 
+-cmap        add CMAP section to data file and fix cmap command lines in
               input script" (NOTE: requires use of *.pdb file)

 In the "example" folder, you will find example files that were created
@ -69,26 +69,26 @@ file has to the corresponding names in the charmm FF files. You'll
 need to add a "pdbalias residue x xnew" line for each change that
 needs to be made.  The *.pgn should contain something like this:

-package require psfgen   
+package require psfgen
 topology top_all27_na.rtf
-pdbalias residue A ADE 
-pdbalias residue T THY 
-pdbalias residue G GUA 
-pdbalias residue C CYT 
+pdbalias residue A ADE
+pdbalias residue T THY
+pdbalias residue G GUA
+pdbalias residue C CYT
 .
 .
 .
-segment A {pdb 1ac7_pared.pdb}   
-coordpdb 1ac7_pared.pdb A   
-guesscoord   
-writepdb 1ac7.pdb   
-writepsf charmm 1ac7.psf  
-exit 
+segment A {pdb 1ac7_pared.pdb}
+coordpdb 1ac7_pared.pdb A
+guesscoord
+writepdb 1ac7.pdb
+writepsf charmm 1ac7.psf
+exit

 5) Type "vmd -e 1ac7.pgn" to build the 1ac7.psf file, and the new
   1ac7.pdb file.

-6) Run charmm2lammps.pl by typing: 
+6) Run charmm2lammps.pl by typing:
 "perl charmm2lammps.pl all27_na 1ac7 -charmm -border=1 -pdb_ctrl -water -ions"

 7) Run lammps by typing: "lmp < 1ac7.in"
@ -105,7 +105,7 @@ molecule. The -pdb_ctrl option produces the 1ac7_ctrl.pdb file that
 can be visualized in a standard visualization package such as VMD. The
 -charmm option put comments into the LAMMPS data file (everything
 after the # sign is a comment) for user convenience in tracking atom
-types etc. according to CHARMM nomenclature.   
+types etc. according to CHARMM nomenclature.

 The example molecule provided above (i.e., 1ac7) is a DNA fragment.
 If instead, a peptide longer than 2 amino acid residues or a protein
--- a/tools/ch2lmp/charmm2lammps.pl
+++ b/tools/ch2lmp/charmm2lammps.pl
--- a/tools/ch2lmp/example/1ac7.data
+++ b/tools/ch2lmp/example/1ac7.data
--- a/tools/ch2lmp/example/1ac7.in
+++ b/tools/ch2lmp/example/1ac7.in
@ -1,4 +1,4 @@
-# Created by charmm2lammps v1.6.2 on Fri Jun  3 14:07:02 MDT 2005
+# Created by charmm2lammps v1.9.0 on Mon Oct  3 06:55:39 EDT 2016

 units           real
 neigh_modify    delay 2 every 1
@ -9,21 +9,23 @@ angle_style     charmm
 dihedral_style  charmm
 improper_style  harmonic

-pair_style      lj/charmm/coul/long 8 10
+pair_style      lj/charmm/coul/long 8 12
 pair_modify     mix arithmetic
-kspace_style    pppm 1e-4
+kspace_style    pppm 1e-6

 read_data       1ac7.data

 special_bonds   charmm
-fix             1 all nve
-fix             2 all shake 1e-6 500 0 m 1.0 a 93
-velocity        all create 0.0 12345678 dist uniform
-
 thermo          1
 thermo_style    multi
 timestep        0.5

+minimize 0.0 0.0 1000 10000
+
+fix             1 all nve
+fix             2 all shake 1e-6 500 0 m 1.0 a 93
+velocity        all create 0.0 12345678 dist uniform
+
 restart         10 1ac7.restart1 1ac7.restart2
 dump            1 all atom 10 1ac7.dump
 dump_modify     1 image yes scale yes
--- a/tools/ch2lmp/example/1ac7.pdb
+++ b/tools/ch2lmp/example/1ac7.pdb
--- a/tools/ch2lmp/example/1ac7.pdb1
+++ b/tools/ch2lmp/example/1ac7.pdb1
--- a/tools/ch2lmp/example/1ac7.pgn
+++ b/tools/ch2lmp/example/1ac7.pgn
--- a/tools/ch2lmp/example/1ac7.psf
+++ b/tools/ch2lmp/example/1ac7.psf
--- a/tools/ch2lmp/example/1ac7_ctrl.pdb
+++ b/tools/ch2lmp/example/1ac7_ctrl.pdb
--- a/tools/ch2lmp/example/1ac7_ctrl.psf
+++ b/tools/ch2lmp/example/1ac7_ctrl.psf
--- a/tools/ch2lmp/example/1ac7_pared.pdb
+++ b/tools/ch2lmp/example/1ac7_pared.pdb
--- a/tools/ch2lmp/example/par_all27_na.prm
+++ b/tools/ch2lmp/example/par_all27_na.prm
--- a/tools/ch2lmp/example/top_all27_na.rtf
+++ b/tools/ch2lmp/example/top_all27_na.rtf
--- a/tools/ch2lmp/lammps2pdb.pl
+++ b/tools/ch2lmp/lammps2pdb.pl
@ -8,9 +8,9 @@
 #  		conjunction with vmd
 #
 #  Notes:	Copyright by author for Sandia National Laboratories
-#    20040903	Conception date of v1.0: rudimentary script for collagen 
+#    20040903	Conception date of v1.0: rudimentary script for collagen
 #    		project.
-#    20050423	Conception date of v2.0: 
+#    20050423	Conception date of v2.0:
 #    		- changed data access through indexing data directly on disk;
 #    		- added all command line options
 #    20050425	Corrected focussing to use a target molecule's moment of
@ -22,21 +22,21 @@
 #    		the data stream.

 # subroutines
-  
+
  sub test
  {
    my $name		= shift(@_);
-    
+
    printf("Error: file %s not found\n", $name) if (!scalar(stat($name)));
    return !scalar(stat($name));
  }

-  
+
  sub initialize
  {
    my $k		= 0;
-    my @options		= ("-help", "-nstart", "-dn", "-cut", "-repair", 
-			   "-units", "-quiet", "-nodetect", "-data", "-pbc", 
+    my @options		= ("-help", "-nstart", "-dn", "-cut", "-repair",
+			   "-units", "-quiet", "-nodetect", "-data", "-pbc",
 			   "-focus", "-center", "-exclude");
    my @remarks		= ("display this message",
 			   "starting frame [-1]",
@ -61,7 +61,7 @@
      "* Expected files in current directory: project.data, project.dump",
      "* Generated output files: project_trj.psf, project_trj.pdb",
      "* Uses project_ctrl.psf if available");
-    
+

    $program		= "lammps2pdb";
    $version		= "2.2.5";
@ -96,7 +96,7 @@
 	$info		= $switch ? 0 : 1 if (!$k--);
 	$detect		= $switch ? 0 : 1 if (!$k--);
 	$data_name	= $arg[1] if (!$k--);
-	if (!$k--) { 
+	if (!$k--) {
 	  if ($switch) { $pbc{ALL} = 1; }
 	  else { foreach (split(",",$arg[1])) { $pbc{uc($_)} = 1; }}}
 	if (!$k--) { foreach (split(",",$arg[1])) { $focus{uc($_)} = uc($_);}}
@ -120,7 +120,7 @@
      printf("\n");
      exit(-1);
    }
-    printf("%s v%s (c)%s\n\n", $program, $version, $year) if ($info);   
+    printf("%s v%s (c)%s\n\n", $program, $version, $year) if ($info);

    $data_name		= $project.".data" if ($data_name eq "");
    $traject_name	= $project.".dump";
@ -133,7 +133,7 @@
    my $flag		= test($data_name);
    printf("Conversion aborted\n\n") if ($flag);
    exit(-1) if ($flag);
-    
+
    # data input

    create_atom_ids();
@ -145,9 +145,9 @@

    open(PSF, ">".$psf_name) if (!$psf_ctrl);
    open(PDB, ">".$pdb_name);
-    
+
    # align center with focus
-    
+
    %center		= %focus if (scalar(%focus));
  }

@ -166,24 +166,24 @@
    return @_;
  }

-  
+
  sub V_Add					# V_Add(@a, @b) = @a + @b;
  {
    return (@_[0]+@_[3], @_[1]+@_[4], @_[2]+@_[5]);
  }
-  
-    
+
+
  sub V_Subtr					# V_Subtr(@a, @b) = @a - @b;
  {
    return (@_[0]-@_[3], @_[1]-@_[4], @_[2]-@_[5]);
  }
-  
-    
+
+
  sub V_Dot					# V_Dot(@a, @b) = @a . @b;
  {
    return (@_[0]*@_[3]+@_[1]*@_[4]+@_[2]*@_[5]);
  }
-  
+

  sub V_Cross					# V_Cross(@a, @b) = @a x @b;
  {
@ -191,7 +191,7 @@
      @_[0]*@_[4]-@_[1]*@_[3]);
  }

-  
+
  sub V_Mult					# V_Mult($f, @a) = $f * @a;
  {
    return (@_[0]*@_[1], @_[0]*@_[2], @_[0]*@_[3]);
@ -203,11 +203,11 @@
    return V_Mult(1/sqrt(V_Dot(@_[0,1,2],@_[0,1,2])), @_[0,1,2]);
  }

-  
+
  sub pbc					# periodic -0.5*L <= x < 0.5*L
  {
    my $x		= @_[0]/@_[1]+0.5;
-    
+
    return @_[0]-@_[1]*($x<0 ? int($x)-1 : int($x));
  }

@ -216,8 +216,8 @@
  {
    return (pbc(@_[0], @_[3]), pbc(@_[1], @_[4]), pbc(@_[2], @_[5]));
  }
-  
-  
+
+
  sub V_Cmp					# V_Cmp(abs(@a), abs(@b))
  {
    return -1 if (abs($_[0])<abs($_[3]));
@ -228,12 +228,12 @@
    return  1 if (abs($_[2])>abs($_[5]));
    return 0;
  }
-  
-    
+
+
  sub V_Sort					# sort on descending absolute
  {						# value
    my @v		= @_;
-    
+
    for (my $i=0; $i<scalar(@v)-3; $i+=3)
    {
      for (my $j=$i+3; $j<scalar(@v); $j+=3)
@ -247,13 +247,13 @@
    return @v;
  }

-  
+
 # Matrix routines

  sub M_String					# M_String(@A)
  {
    my @M;
-    
+
    for (my $i=0; $i<3; ++$i) { push(@M, V_String(splice(@_, 0, 3))); }
    return "{".join(", ", @M)."}";
  }
@ -270,7 +270,7 @@
    return (@_[0], @_[3], @_[6], @_[1], @_[4], @_[7], @_[2], @_[5], @_[8]);
  }

- 
+
  sub M_Dot					# M_Dot(@A, @B) = @A . @B;
  {
    return (
@ -288,7 +288,7 @@
    return V_Dot(@_[0,1,2], V_Cross(@_[3,4,5], @_[6,7,8]));
  }

-  
+
  sub M_Mult					# M_Mult($a, @A) = $a * @A
  {
    return (
@ -296,12 +296,12 @@
      @_[0]*@_[4], @_[0]*@_[5], @_[0]*@_[6],
      @_[0]*@_[7], @_[0]*@_[8], @_[0]*@_[9]);
  }
-  
+

  sub M_Unit { return (1,0,0, 0,1,0, 0,0,1); }

  sub PI { return 4*atan2(1,1); }
-  
+
  sub M_Rotate					# M_Rotate($n, $alpha) = @R_$n;
  {						# vmd convention
    my $n		= shift(@_);
@ -316,7 +316,7 @@
    return ($cos,-$sin,0, $sin,$cos,0, 0,0,1) if ($n==2); # around z-axis
    return M_Unit();
  }
-    
+

  sub M_RotateNormal				# returns R.(1,0,0) = @a/|@a|;
  {
@ -324,7 +324,7 @@
    my @n		= V_Mult(1.0/sqrt(V_Dot(@_[0,1,2], @_)), @_);
    my $sina		= $n[1]<0 ? -sqrt($n[1]*$n[1]+$n[2]*$n[2]) :
 				     sqrt($n[1]*$n[1]+$n[2]*$n[2]);
-    
+
    if ($sina)
    {
      my $cosa		= $n[0];
@ -337,8 +337,8 @@
    }
    return @R;
  }
-  
-  
+
+
  sub MV_Dot					# MV_Dot(@A, @b) = @A . @b;
  {
    return (V_Dot(@_[0,1,2], @_[9,10,11]), V_Dot(@_[3,4,5], @_[9,10,11]),
@ -374,14 +374,14 @@
  {
    return (@_[0]+@_[2], @_[1]+@_[3]);
  }
-  
-    
+
+
  sub C_Subtr					# z = z1 - z2
  {
    return (@_[0]-@_[2], @_[1]-@_[3]);
  }
-  
-    
+
+
  sub C_Mult					# z = z1 * z2
  {
    return (@_[0]*@_[2]-@_[1]*@_[3], @_[0]*@_[3]+@_[2]*@_[1]);
@ -404,7 +404,7 @@
    return ($r*cos($a), $r*sin($a));
  }

-  
+
  sub C_Correct
  {
    return (abs(@_[0])<1e-14 ? 0 : @_[0], abs(@_[1])<1e-14 ? 0 : @_[1]);
@ -422,7 +422,7 @@
    return (@_[0], -@_[1]);
  }

-  
+
  sub C_String
  {
    return @_[0]." + ".@_[1]."i";
@ -434,12 +434,12 @@
  sub R_Sort
  {
    my $n		= scalar(@_);
-    
+
    for (my $i=0; $i<$n-2; $i+=2)
    {
      for (my $j=$i+2; $j<$n; $j+=2)
      {
-	if (@_[$j]<@_[$i]) { 
+	if (@_[$j]<@_[$i]) {
 	  my @t = @_[$i,$i+1]; @_[$i,$i+1] = @_[$j,$j+1]; @_[$j,$j+1] = @t; }
 	else { if ((@_[$j]==@_[$i])&&(@_[$j+1]<@_[$i+1])) {
 	  my @t = @_[$i,$i+1]; @_[$i,$i+1] = @_[$j,$j+1]; @_[$j,$j+1] = @t; } }
@ -448,7 +448,7 @@
    return @_;
  }

-  
+
  sub R_First
  {
    return (0, 0) if (abs(@_[1])<1e-14);
@ -476,16 +476,16 @@
    my @B		= (0, 0);
    my @z1		= (0.5, 0.5*sqrt(3));
    my @z2		= C_Conj(@z1);
-   
+
    if (abs($f1)<1e-3) { 				# limit f1 -> 0
      @A		= ($f2<0 ? abs(2*$f2)**(1/3) : 0, 0); }
-    else { 
+    else {
      if (abs($f2)<1e-14) {				# limit f2 -> 0
 	my $f		= sqrt(abs($f1))/$c3;
 	@A		= $f1<0 ? (-$f*$z1[1], 0.5*$f) : ($f, 0);
 	@B		= $f1<0 ? (-$A[0], $A[1]) : ($f, 0); }
      else {
-	@B		= C_Pow(C_Add(($f2, 0), 
+	@B		= C_Pow(C_Add(($f2, 0),
 			    C_Pow(($f1*$f1*$f1+$f2*$f2, 0), 1/2)), 1/3);
 	@A		= C_Div(($f1/$c3, 0), @B);
 	@B		= ($B[0]/$c3, $B[1]/$c3); } }
@ -503,21 +503,21 @@
    my $input		= shift;
    my $dlines		= shift;
    my $read;
-    
+
    while ($dlines--) { $read = <$input>; }
    return $read;
  }

-  
+
  sub rewind_read
  {
    my $input		= shift;
-    
+
    seek($input, 0, SEEK_SET);
  }


-# create crossreference tables 
+# create crossreference tables

  sub create_psf_index					# make an index of id
  {							# locations
@ -538,13 +538,13 @@
    }
  }

- 
+
  sub psf_goto						# goto $ident in <PSF>
  {
    create_psf_index() if (!scalar(%PSFIndex));
    my $id		= shift(@_);
    my @n 		= split(" ", $PSFIndex{$id});
-    
+
    @PSFBuffer		= ();
    if (!scalar(@n))
    {
@ -571,7 +571,7 @@
    my $id;
    my %hash;
    my %size;
-    
+
    foreach ((masses,atoms,bonds,angles,dihedrals,impropers)) { $hash{$_}=$_; }
    open(DATA, "<".$data_name);
    for (my $i=0; $i<2; ++$i) { my $tmp = <DATA>; }	# skip first two lines
@ -604,7 +604,7 @@
    }
  }

-  
+
  sub goto_data
  {
    create_data_index() if (!scalar(%DATAIndex));
@ -623,7 +623,7 @@


 # create atom and residue identifiers
-  
+
  sub create_names
  {
    return if (scalar(@names));
@ -636,7 +636,7 @@
    my @letter		= ("X", "Y", "Z", "P", "Q", "R", "S", "A", "B", "C");
    my $k		= 0;
    my %atom;
-    
+
    $names[0]		= "";
    foreach (@mass) { $atom{$_} = shift(@name); }
    for (my $i=1; $i<=$n; ++$i)
@ -658,7 +658,7 @@
    my @data		= @_;
    my $p		= $data[1]." ".$data[2];
    my $k;
-    
+
    for ($k=0; ($k<scalar(@data))&&(substr($data[$k],0,1) ne "#"); ++$k) { }
    @data		= splice(@data, $k-($k<8 ? 3 : 6), $k<8 ? 3 : 6);
    foreach (@L)
@ -669,8 +669,8 @@
    }
    return $p;
  }
-  
-  
+
+
  sub create_atom_ids
  {
    my $res		= 0;
@ -682,7 +682,7 @@
    my $id;
    my %link;
    my %special;
-    
+
    printf("Info: creating atom ids\n") if ($info);
    create_names();
    $n			= goto_data(atoms);
@ -696,8 +696,8 @@
      {
 	if ((($tmp = $link{$id = join(" ", sort(split(" ", $id)))}) eq "")&&
 	    (($tmp = $special{$id}) eq ""))
-	{  
-	  $link{$id}	= 
+	{
+	  $link{$id}	=
 	      $tmp 	= "R".($res<10 ? "0" : "").$res;
 	  ++$res;
 	}
@ -712,12 +712,12 @@
    }
  }

-  
+
  sub crossover
  {
    my @d		= V_Subtr((split(" ", $position[@_[0]]))[0,1,2],
 				  (split(" ", $position[@_[1]]))[0,1,2]);
-  
+
    $d[0]		/= $l[3];
    $d[1]		/= $l[4];
    $d[2]		/= $l[5];
@ -744,7 +744,7 @@
  {
    my $n		= scalar(@bonds);
    my $i		= 0;
-    
+
    printf("Info: deleting excluded bonds\n") if ($info);
    while ($i<$n)
    {
@ -754,12 +754,12 @@
      else { ++$i; }
    }
  }
-  
-  
+
+
  sub create_bonds
  {
    my $n		= goto_data(bonds);
-    
+
    printf("Info: creating bonds\n") if ($info);
    for (my $i=0; $i<$n; ++$i)
    {
@ -779,20 +779,20 @@
    while (!eof(TRAJECT)&&(substr(lc(join(" ", split(" ", <TRAJECT>))),
 				  0,length($subject)) ne $subject)) {}
  }
-  
-  
+
+
  sub read_traject
  {
    my @box;
    my @l;
-    
+
    advance_traject("timestep");
    my $timestep	= (split(" ", <TRAJECT>))[0];
    advance_traject("number of atoms");
    my $n		= (split(" ", <TRAJECT>))[0];
    advance_traject("box bounds");
    for (my $i=0; $i<3; ++$i)
-    { 
+    {
      my @data		= split(" ", <TRAJECT>);
      $box[$i]		= $data[0];			# box edge
      $l[$i]		= $data[1]-$data[0];		# box length
@ -806,9 +806,9 @@
    return ($timestep, $n, @box, @l);
  }

-  
+
 # pdb format
- 
+
  sub eigen_vector				# eigen_vector(@A, $l)
  {
    my @A		= splice(@_,0,9);
@ -826,8 +826,8 @@
    return (0,0,1) if ($A[8]==1);
    return (0,0,0);
  }
-  
-    
+
+
  sub pdb_inertia
  {
    my @s		= (
@ -852,13 +852,13 @@
 			    M_Transpose(M_RotateNormal(MV_Dot(@A,@b))));
    return M_Dot(@B, @A);
  }
-    
-      
+
+

  sub pdb_focus		# using moment of inertia
  {
    my @R		= pdb_inertia(@_);
-    
+
    printf("Info: focussing\n") if ($info);
    foreach (@position)
    {
@ -867,7 +867,7 @@
    }
  }

-  
+
  sub pdb_center	
  {
    my @c		= splice(@_, 0, 3);
@ -881,7 +881,7 @@
    }
  }

-  
+
  sub pdb_pbc
  {
    printf("Info: applying periodicity\n") if ($info);
@ -909,7 +909,7 @@
    }
  }

-  
+
  sub pdb_positions
  {
    my @m		= (0,0,0,0,0,0,0,0,0);
@ -918,7 +918,7 @@
    my $mass;
    my @p;
    my $d;
-  
+
    foreach (@traject)
    {
      my @arg		= split(" ");
@ -974,7 +974,7 @@
    printf(PDB "%-11.11s", "P 1");
    printf(PDB "%3.3s\n", "1");
  }
-  
+

  sub pdb_atoms
  {
@ -998,15 +998,15 @@
      foreach (@p) { $_ = 0 if (abs($_)<1e-4); }
      printf(PDB "ATOM  ");				# pdb command
      printf(PDB "%6.6s ", ++$n);			# atom number
-      printf(PDB "%-3.3s ", 
+      printf(PDB "%-3.3s ",
 	$psf_ctrl ? $psf[4] : $names[$p[4]]);		# atom name
-      printf(PDB "%-3.3s ", 
+      printf(PDB "%-3.3s ",
 	$psf_ctrl ? $psf[3] : $residue[$nres]);		# residue name
      printf(PDB "%5.5s ", $nres);			# residue number
      printf(PDB "%3.3s ", "");				# empty placeholder
      printf(PDB "%7.7s %7.7s %7.7s ",
       	$p[0], $p[1], $p[2]);				# positions
-      printf(PDB "%5.5s %5.5s %4.4s ", 
+      printf(PDB "%5.5s %5.5s %4.4s ",
 	"1.00", "0.00", "");				# trailing variables
      printf(PDB "%-4.4s\n",
       	$psf_ctrl ? $psf[1] : $cluster[$nres]);		# cluster name
@ -1015,7 +1015,7 @@
    };
    printf(PDB "END\n");
  }
- 
+

  sub pdb_timestep
  {
@ -1045,7 +1045,7 @@
    my $l		= 0;
    my $k		= 0;
    my @extra;
-    
+
    for (my $i=0; $i<$n; ++$i)
    {
      my @arg		= split(" ", <DATA>);
@ -1073,13 +1073,13 @@
    }
    printf(PSF "\n");
  }
-  
+

  sub psf_bonds
  {
    my $npairs		= 0;
-  
-    delete_exclude() if (scalar(%exclude)>0); 
+
+    delete_exclude() if (scalar(%exclude)>0);
    delete_crossovers() if ($cut);
    printf(PSF "%8.8s !NBOND\n", scalar(@bonds));
    foreach (@bonds)
@ -1098,7 +1098,7 @@
  initialize();

  # create .pdb file
- 
+
  $ncurrent		= -1;
  while ($traject_flag&&!eof(TRAJECT))
  {
@ -1134,4 +1134,4 @@

  close(TRAJECT);
  close(DATA);
-  
+
--- a/tools/ch2lmp/other/in_mkpdb
+++ b/tools/ch2lmp/other/in_mkpdb
Author	SHA1	Message	Date
Steve Plimpton	02bfa898ee	adjustments to balancing weights and factors, also XOR op for formulas, if, dump_modify thresh	2016-10-05 15:46:20 -06:00
sjplimp	030df745bc	Merge pull request #193 from akohlmey/eam-bugfix bugfix for eam/alloy/omp and eam/fs/omp	2016-10-05 10:54:36 -06:00
sjplimp	6a97211932	Merge pull request #192 from rbberger/python-interface-bugfix Revert type checking commit from July	2016-10-05 10:54:08 -06:00
Steve Plimpton	c46be7db62	changes to imbalance weight factors	2016-10-05 10:33:39 -06:00
Axel Kohlmeyer	4381db846b	correct the bug discovered by stan due to uninitialized scale factors for eam/alloy/omp and eam/fs/omp	2016-10-04 14:33:26 -04:00
Richard Berger	e2caf5c105	Fix code path which allows passing a C++ ptr to PyLammps	2016-10-04 13:57:21 -04:00
Steve Plimpton	11c2892e54	Merge branch 'restrict-weights-and-weight-factors' of https://github.com/akohlmey/lammps	2016-10-04 09:49:09 -06:00
Richard Berger	91be47a0d0	Revert type checking commit from July `0aebb2eabe`	2016-10-04 11:43:12 -04:00
sjplimp	ab92529b19	Merge pull request #191 from akohlmey/updated-charmm2lammps Updated charmm2lammps	2016-10-03 17:59:21 -06:00
sjplimp	e079362776	Merge pull request #190 from akohlmey/small-bufixes-and-enhancements Small bufixes and enhancements	2016-10-03 17:58:36 -06:00
Steve Plimpton	c3ff8812b3	added XOR operator to variable command	2016-10-03 17:57:33 -06:00
Axel Kohlmeyer	03766dbda7	apply bugfix for MEAM provided by Wolfgang Verestek on lammps-users this closes lammps/#188	2016-10-03 16:28:59 -04:00
Axel Kohlmeyer	6e719f2d94	remove trailing whitespace	2016-10-03 07:07:28 -04:00
Axel Kohlmeyer	45d2cc2895	permission update for ch2lmp tool folder	2016-10-03 07:03:42 -04:00
Axel Kohlmeyer	690f91300b	rebuild charmm2lammps example output files with updated tools	2016-10-03 06:58:51 -04:00
Axel Kohlmeyer	3b94627dfe	properly handle -nohints flag, make -cmap flag take version as option. step version number	2016-10-03 06:52:30 -04:00
Axel Kohlmeyer	c2e11dffa2	import updated charmm2lammps.pl script from Rober Latour	2016-10-02 20:33:20 -04:00
Axel Kohlmeyer	1985db4fb1	correct designation of meam supporting USER-OMP and meam/spline not	2016-10-01 23:05:05 -04:00
Axel Kohlmeyer	a3e05a2bac	permission cleanup	2016-10-01 06:34:45 -04:00
Axel Kohlmeyer	035279de87	correct logic bug in bufix for fix tmd (cherry picked from commit `267c1ec957`)	2016-10-01 06:26:52 -04:00
Axel Kohlmeyer	d45e333f7c	restrict choice of weight factors and guarantee that weights are >= 0.001	2016-09-30 11:11:32 -04:00