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
141 changed files with 29679 additions and 8880 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="28 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
-28 Sep 2016 version :c,h4
+5 Oct 2016 version :c,h4

 Version info: :h4

--- a/doc/src/PDF/colvars-refman-lammps.pdf
+++ b/doc/src/PDF/colvars-refman-lammps.pdf
--- 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/fix_colvars.txt
+++ b/doc/src/fix_colvars.txt
@ -39,8 +39,8 @@ metadynamics, Steered Molecular Dynamics (SMD) and Umbrella Sampling
 (US) via a flexible harmonic restraint bias. The colvars library is
 hosted at "http://colvars.github.io/"_http://colvars.github.io/

-This documentation describes only the fix colvars command itself and 
-LAMMPS specific parts of the code.  The full documentation of the 
+This documentation describes only the fix colvars command itself and
+LAMMPS specific parts of the code.  The full documentation of the
 colvars library is available as "this supplementary PDF document"_PDF/colvars-refman-lammps.pdf

 A detailed discussion of the implementation of the portable collective
@ -122,7 +122,7 @@ not a limitation of functionality.

 [Default:]

-The default options are input = NULL, output = out, seed = 1966, unwrap yes, 
+The default options are input = NULL, output = out, seed = 1966, unwrap yes,
 and tstat = NULL.

 :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/examples/cmap/charmm22.cmap
+++ b/examples/cmap/charmm22.cmap
@ -1,4 +1,5 @@
-# Title: charmm correction map
+# DATE: 2016-09-26 CONTRIBUTOR: Robert Latour, latourr@clemson.edu  CITATION: TBA
+# Title: charmm22/charmm27 dihedral correction map

 #  alanine map, type 1

--- a/examples/prd/README
+++ b/examples/prd/README
@ -1,5 +1,6 @@
 Run this example as:

-mpirun -np 4 lmp_linux -partition 4x1 -in in.prd
+mpirun -np 4 lmp_g++ -partition 4x1 -in in.prd
+mpirun -np 8 lmp_g++ -partition 4x2 -in in.prd

 You should be able to use any number of replicas >= 3.
--- a/examples/prd/log.15Feb16.prd.g++.4
+++ b/examples/prd/log.15Feb16.prd.g++.4
@ -1,22 +0,0 @@
-LAMMPS (15 Feb 2016)
-Running on 4 partitions of processors
-Step CPU Clock Event Correlated Coincident Replica
-100 0.000 0 0 0 0 0
-100 0.539 0 0 0 0 0
-100 0.694 0 0 0 0 0
-100 0.850 0 0 0 0 0
-100 0.928 0 0 0 0 0
-200 1.555 400 1 0 4 1
-300 1.924 500 2 1 1 1
-500 3.495 772 3 0 2 2
-800 5.446 1328 4 0 2 2
-1000 7.037 1636 5 0 1 3
-1000 7.989 1636 5 0 1 3
-1200 8.838 1908 6 0 1 2
-1300 9.212 2008 7 1 1 2
-1900 12.291 4024 8 0 1 3
-1900 12.832 4024 8 0 1 3
-1900 13.157 4024 8 0 1 3
-1900 13.320 4024 8 0 1 3
-2100 14.090 4220 9 0 2 3
-Loop time of 14.0941 on 4 procs for 2000 steps with 511 atoms
--- a/examples/prd/log.29Sep16.prd.g++.4
+++ b/examples/prd/log.29Sep16.prd.g++.4
@ -0,0 +1,22 @@
+LAMMPS (28 Sep 2016)
+Running on 4 partitions of processors
+Step CPU Clock Event Correlated Coincident Replica
+100 0.000 0 0 0 0 0
+100 0.521 0 0 0 0 0
+100 0.670 0 0 0 0 0
+100 0.822 0 0 0 0 0
+100 0.896 0 0 0 0 0
+200 1.555 400 1 0 4 1
+300 1.918 500 2 1 1 1
+500 3.476 772 3 0 2 2
+800 5.379 1328 4 0 2 2
+1000 6.914 1636 5 0 1 3
+1000 7.859 1636 5 0 1 3
+1200 8.658 1908 6 0 1 2
+1300 9.018 2008 7 1 1 2
+1900 12.005 4024 8 0 1 3
+1900 12.539 4024 8 0 1 3
+1900 12.861 4024 8 0 1 3
+1900 13.027 4024 8 0 1 3
+2100 13.798 4220 9 0 2 3
+Loop time of 13.8021 on 4 procs for 2000 steps with 511 atoms
--- a/examples/prd/log.29Sep16.prd.g++.8
+++ b/examples/prd/log.29Sep16.prd.g++.8
@ -0,0 +1,22 @@
+LAMMPS (28 Sep 2016)
+Running on 4 partitions of processors
+Step CPU Clock Event Correlated Coincident Replica
+100 0.000 0 0 0 0 0
+100 0.289 0 0 0 0 0
+100 0.373 0 0 0 0 0
+100 0.458 0 0 0 0 0
+100 0.500 0 0 0 0 0
+200 0.800 400 1 0 4 1
+300 0.999 500 2 1 1 1
+500 1.834 772 3 0 2 2
+800 2.864 1328 4 0 2 2
+1000 3.706 1636 5 0 1 3
+1000 4.211 1636 5 0 1 3
+1200 4.655 1908 6 0 1 2
+1300 4.854 2008 7 1 1 2
+1900 6.503 4024 8 0 1 3
+1900 6.792 4024 8 0 1 3
+1900 6.968 4024 8 0 1 3
+1900 7.055 4024 8 0 1 3
+2100 7.463 4220 9 0 2 3
+Loop time of 7.46531 on 8 procs for 2000 steps with 511 atoms
--- a/examples/prd/screen.29Sep16.prd.g++.4.0
+++ b/examples/prd/screen.29Sep16.prd.g++.4.0
--- a/examples/prd/screen.29Sep16.prd.g++.4.1
+++ b/examples/prd/screen.29Sep16.prd.g++.4.1
--- a/examples/prd/screen.29Sep16.prd.g++.4.2
+++ b/examples/prd/screen.29Sep16.prd.g++.4.2
--- a/examples/prd/screen.29Sep16.prd.g++.4.3
+++ b/examples/prd/screen.29Sep16.prd.g++.4.3
--- a/examples/prd/screen.29Sep16.prd.g++.8.0
+++ b/examples/prd/screen.29Sep16.prd.g++.8.0
--- a/examples/prd/screen.29Sep16.prd.g++.8.1
+++ b/examples/prd/screen.29Sep16.prd.g++.8.1
--- a/examples/prd/screen.29Sep16.prd.g++.8.2
+++ b/examples/prd/screen.29Sep16.prd.g++.8.2
--- a/examples/prd/screen.29Sep16.prd.g++.8.3
+++ b/examples/prd/screen.29Sep16.prd.g++.8.3
--- a/examples/tad/README
+++ b/examples/tad/README
@ -1,5 +1,5 @@
 Run this example as:

-mpirun -np 3 lmp_linux -partition 3x1 -in in.tad
+mpirun -np 3 lmp_g++ -partition 3x1 -in in.tad

 You should be able to use any number of replicas >= 3.
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@ -150,7 +150,7 @@ colvar::colvar(std::string const &conf)
    feature_states[f_cv_linear]->enabled = lin;
  }

-  // Colvar is homogeneous iff:
+  // Colvar is homogeneous if:
  // - it is linear (hence not scripted)
  // - all cvcs have coefficient 1 or -1
  // i.e. sum or difference of cvcs
@ -375,28 +375,16 @@ colvar::colvar(std::string const &conf)

  {
    bool temp;
-    if (get_keyval(conf, "outputSystemForce", temp, false)) {
-      cvm::error("Colvar option outputSystemForce is deprecated.\n"
-        "Please use outputTotalForce, or outputSystemForce within an ABF bias.");
+    if (get_keyval(conf, "outputSystemForce", temp, false, colvarparse::parse_silent)) {
+      cvm::error("Option outputSystemForce is deprecated: only outputTotalForce is supported instead.\n"
+                 "The two are NOT identical: see http://colvars.github.io/totalforce.html.\n", INPUT_ERROR);
      return;
    }
  }

-  {
-    bool b_output_total_force;
-    get_keyval(conf, "outputTotalForce", b_output_total_force, false);
-    if (b_output_total_force) {
-      enable(f_cv_output_total_force);
-    }
-  }
-
-  {
-    bool b_output_applied_force;
-    get_keyval(conf, "outputAppliedForce", b_output_applied_force, false);
-    if (b_output_applied_force) {
-      enable(f_cv_output_applied_force);
-    }
-  }
+  get_keyval_feature(this, conf, "outputTotalForce", f_cv_output_total_force, false);
+  get_keyval_feature(this, conf, "outputAppliedForce", f_cv_output_applied_force, false);
+  get_keyval_feature(this, conf, "subtractAppliedForce", f_cv_subtract_applied_force, false);

  // Start in active state by default
  enable(f_cv_active);
@ -409,6 +397,8 @@ colvar::colvar(std::string const &conf)
  fj.type(value());
  ft.type(value());
  ft_reported.type(value());
+  f_old.type(value());
+  f_old.reset();

  if (cvm::b_analysis)
    parse_analysis(conf);
@ -519,6 +509,8 @@ int colvar::init_components(std::string const &conf)
    "number", "coordNum");
  error_code |= init_components_type<selfcoordnum>(conf, "self-coordination "
    "number", "selfCoordNum");
+  error_code |= init_components_type<groupcoordnum>(conf, "group-coordination "
+    "number", "groupCoord");
  error_code |= init_components_type<angle>(conf, "angle", "angle");
  error_code |= init_components_type<dipole_angle>(conf, "dipole angle", "dipoleAngle");
  error_code |= init_components_type<dihedral>(conf, "dihedral", "dihedral");
@ -1104,6 +1096,14 @@ int colvar::calc_colvar_properties()

  } else {

+    if (is_enabled(f_cv_subtract_applied_force)) {
+      // correct the total force only if it has been measured
+      // TODO add a specific test instead of relying on sq norm
+      if (ft.norm2() > 0.0) {
+        ft -= f_old;
+      }
+    }
+
    x_reported = x;
    ft_reported = ft;
  }
@ -1210,6 +1210,10 @@ cvm::real colvar::update_forces_energy()
    x_old = x;
  }

+  if (is_enabled(f_cv_subtract_applied_force)) {
+    f_old = f;
+  }
+
  if (cvm::debug())
    cvm::log("Done updating colvar \""+this->name+"\".\n");
  return (potential_energy + kinetic_energy);
@ -1640,15 +1644,9 @@ std::ostream & colvar::write_traj(std::ostream &os)
  }

  if (is_enabled(f_cv_output_applied_force)) {
-    if (is_enabled(f_cv_extended_Lagrangian)) {
-      os << " "
-         << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
-         << fr;
-    } else {
-      os << " "
-         << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
-         << f;
-    }
+    os << " "
+       << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
+       << applied_force();
  }

  return os;
--- a/lib/colvars/colvar.h
+++ b/lib/colvars/colvar.h
@ -175,6 +175,9 @@ public:
  /// (if defined) contribute to it
  colvarvalue f;

+  /// Applied force at the previous step (to be subtracted from total force if needed)
+  colvarvalue f_old;
+
  /// \brief Total force, as derived from the atomic trajectory;
  /// should equal the system force plus \link f \endlink
  colvarvalue ft;
@ -272,10 +275,13 @@ public:
  /// \brief Calculate the quantities associated to the colvar (but not to the CVCs)
  int calc_colvar_properties();

-  /// Get the current biasing force
-  inline colvarvalue bias_force() const
+  /// Get the current applied force
+  inline colvarvalue const applied_force() const
  {
-    return fb;
+    if (is_enabled(f_cv_extended_Lagrangian)) {
+      return fr;
+    }
+    return f;
  }

  /// Set the total biasing force to zero
@ -482,6 +488,7 @@ public:
  class dihedral;
  class coordnum;
  class selfcoordnum;
+  class groupcoordnum;
  class h_bond;
  class rmsd;
  class orientation_angle;
--- a/lib/colvars/colvarbias_histogram.cpp
+++ b/lib/colvars/colvarbias_histogram.cpp
@ -67,7 +67,7 @@ int colvarbias_histogram::init(std::string const &conf)

  if (colvar_array_size > 0) {
    weights.assign(colvar_array_size, 1.0);
-    get_keyval(conf, "weights", weights, weights, colvarparse::parse_silent);
+    get_keyval(conf, "weights", weights, weights);
  }

  for (i = 0; i < colvars.size(); i++) {
@ -79,7 +79,7 @@ int colvarbias_histogram::init(std::string const &conf)

  {
    std::string grid_conf;
-    if (key_lookup(conf, "grid", grid_conf)) {
+    if (key_lookup(conf, "histogramGrid", grid_conf)) {
      grid->parse_params(grid_conf);
    }
  }
--- a/lib/colvars/colvarbias_meta.cpp
+++ b/lib/colvars/colvarbias_meta.cpp
@ -92,7 +92,11 @@ int colvarbias_meta::init(std::string const &conf)
    get_keyval(conf, "keepHills", keep_hills, false);
    if (! get_keyval(conf, "writeFreeEnergyFile", dump_fes, true))
      get_keyval(conf, "dumpFreeEnergyFile", dump_fes, true, colvarparse::parse_silent);
-    get_keyval(conf, "saveFreeEnergyFile", dump_fes_save, false);
+    if (get_keyval(conf, "saveFreeEnergyFile", dump_fes_save, false, colvarparse::parse_silent)) {
+      cvm::log("Option \"saveFreeEnergyFile\" is deprecated, "
+               "please use \"keepFreeEnergyFiles\" instead.");
+    }
+    get_keyval(conf, "keepFreeEnergyFiles", dump_fes_save, dump_fes_save);

    hills_energy           = new colvar_grid_scalar(colvars);
    hills_energy_gradients = new colvar_grid_gradient(colvars);
--- a/lib/colvars/colvarbias_restraint.cpp
+++ b/lib/colvars/colvarbias_restraint.cpp
@ -612,3 +612,250 @@ cvm::real colvarbias_restraint_linear::restraint_convert_k(cvm::real k,



+colvarbias_restraint_histogram::colvarbias_restraint_histogram(char const *key)
+  : colvarbias(key)
+{
+  lower_boundary = 0.0;
+  upper_boundary = 0.0;
+  width = 0.0;
+  gaussian_width = 0.0;
+}
+
+
+int colvarbias_restraint_histogram::init(std::string const &conf)
+{
+  colvarbias::init(conf);
+
+  get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary);
+  get_keyval(conf, "upperBoundary", upper_boundary, upper_boundary);
+  get_keyval(conf, "width", width, width);
+
+  if (width <= 0.0) {
+    cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR);
+  }
+
+  get_keyval(conf, "gaussianWidth", gaussian_width, 2.0 * width, colvarparse::parse_silent);
+  get_keyval(conf, "gaussianSigma", gaussian_width, 2.0 * width);
+
+  if (lower_boundary >= upper_boundary) {
+    cvm::error("Error: the upper boundary, "+
+               cvm::to_str(upper_boundary)+
+               ", is not higher than the lower boundary, "+
+               cvm::to_str(lower_boundary)+".\n",
+               INPUT_ERROR);
+  }
+
+  cvm::real const nbins = (upper_boundary - lower_boundary) / width;
+  int const nbins_round = (int)(nbins);
+
+  if (std::fabs(nbins - cvm::real(nbins_round)) > 1.0E-10) {
+    cvm::log("Warning: grid interval ("+
+             cvm::to_str(lower_boundary, cvm::cv_width, cvm::cv_prec)+" - "+
+             cvm::to_str(upper_boundary, cvm::cv_width, cvm::cv_prec)+
+             ") is not commensurate to its bin width ("+
+             cvm::to_str(width, cvm::cv_width, cvm::cv_prec)+").\n");
+  }
+
+  p.resize(nbins_round);
+  ref_p.resize(nbins_round);
+  p_diff.resize(nbins_round);
+
+  bool const inline_ref_p =
+    get_keyval(conf, "refHistogram", ref_p.data_array(), ref_p.data_array());
+  std::string ref_p_file;
+  get_keyval(conf, "refHistogramFile", ref_p_file, std::string(""));
+  if (ref_p_file.size()) {
+    if (inline_ref_p) {
+      cvm::error("Error: cannot specify both refHistogram and refHistogramFile at the same time.\n",
+                 INPUT_ERROR);
+    } else {
+      std::ifstream is(ref_p_file.c_str());
+      std::string data_s = "";
+      std::string line;
+      while (getline_nocomments(is, line)) {
+        data_s.append(line+"\n");
+      }
+      if (data_s.size() == 0) {
+        cvm::error("Error: file \""+ref_p_file+"\" empty or unreadable.\n", FILE_ERROR);
+      }
+      is.close();
+      cvm::vector1d<cvm::real> data;
+      if (data.from_simple_string(data_s) != 0) {
+        cvm::error("Error: could not read histogram from file \""+ref_p_file+"\".\n");
+      }
+      if (data.size() == 2*ref_p.size()) {
+        // file contains both x and p(x)
+        size_t i;
+        for (i = 0; i < ref_p.size(); i++) {
+          ref_p[i] = data[2*i+1];
+        }
+      } else if (data.size() == ref_p.size()) {
+        ref_p = data;
+      } else {
+        cvm::error("Error: file \""+ref_p_file+"\" contains a histogram of different length.\n",
+                   INPUT_ERROR);
+      }
+    }
+  }
+  cvm::real const ref_integral = ref_p.sum() * width;
+  if (std::fabs(ref_integral - 1.0) > 1.0e-03) {
+    cvm::log("Reference distribution not normalized, normalizing to unity.\n");
+    ref_p /= ref_integral;
+  }
+
+  get_keyval(conf, "writeHistogram", b_write_histogram, false);
+  get_keyval(conf, "forceConstant", force_k, 1.0);
+
+  return COLVARS_OK;
+}
+
+
+colvarbias_restraint_histogram::~colvarbias_restraint_histogram()
+{
+  p.resize(0);
+  ref_p.resize(0);
+  p_diff.resize(0);
+}
+
+
+int colvarbias_restraint_histogram::update()
+{
+  if (cvm::debug())
+    cvm::log("Updating the histogram restraint bias \""+this->name+"\".\n");
+
+  size_t vector_size = 0;
+  size_t icv;
+  for (icv = 0; icv < colvars.size(); icv++) {
+    vector_size += colvars[icv]->value().size();
+  }
+
+  cvm::real const norm = 1.0/(std::sqrt(2.0*PI)*gaussian_width*vector_size);
+
+  // calculate the histogram
+  p.reset();
+  for (icv = 0; icv < colvars.size(); icv++) {
+    colvarvalue const &cv = colvars[icv]->value();
+    if (cv.type() == colvarvalue::type_scalar) {
+      cvm::real const cv_value = cv.real_value;
+      size_t igrid;
+      for (igrid = 0; igrid < p.size(); igrid++) {
+        cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
+        p[igrid] += norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
+                                    (2.0 * gaussian_width * gaussian_width));
+      }
+    } else if (cv.type() == colvarvalue::type_vector) {
+      size_t idim;
+      for (idim = 0; idim < cv.vector1d_value.size(); idim++) {
+        cvm::real const cv_value = cv.vector1d_value[idim];
+        size_t igrid;
+        for (igrid = 0; igrid < p.size(); igrid++) {
+          cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
+          p[igrid] += norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
+                                      (2.0 * gaussian_width * gaussian_width));
+        }
+      }
+    } else {
+      // TODO
+    }
+  }
+
+  cvm::real const force_k_cv = force_k * vector_size;
+
+  // calculate the difference between current and reference
+  p_diff = p - ref_p;
+  bias_energy = 0.5 * force_k_cv * p_diff * p_diff;
+
+  // calculate the forces
+  for (icv = 0; icv < colvars.size(); icv++) {
+    colvarvalue const &cv = colvars[icv]->value();
+    colvarvalue &cv_force = colvar_forces[icv];
+    cv_force.type(cv);
+    cv_force.reset();
+
+    if (cv.type() == colvarvalue::type_scalar) {
+      cvm::real const cv_value = cv.real_value;
+      cvm::real &force = cv_force.real_value;
+      size_t igrid;
+      for (igrid = 0; igrid < p.size(); igrid++) {
+        cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
+        force += force_k_cv * p_diff[igrid] *
+          norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
+                          (2.0 * gaussian_width * gaussian_width)) *
+          (-1.0 * (x_grid - cv_value) / (gaussian_width * gaussian_width));
+      }
+    } else if (cv.type() == colvarvalue::type_vector) {
+      size_t idim;
+      for (idim = 0; idim < cv.vector1d_value.size(); idim++) {
+        cvm::real const cv_value = cv.vector1d_value[idim];
+        cvm::real &force = cv_force.vector1d_value[idim];
+        size_t igrid;
+        for (igrid = 0; igrid < p.size(); igrid++) {
+          cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
+          force += force_k_cv * p_diff[igrid] *
+            norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
+                            (2.0 * gaussian_width * gaussian_width)) *
+            (-1.0 * (x_grid - cv_value) / (gaussian_width * gaussian_width));
+        }
+      }
+    } else {
+      // TODO
+    }
+  }
+
+  return COLVARS_OK;
+}
+
+
+std::ostream & colvarbias_restraint_histogram::write_restart(std::ostream &os)
+{
+  if (b_write_histogram) {
+    std::string file_name(cvm::output_prefix+"."+this->name+".hist.dat");
+    std::ofstream os(file_name.c_str());
+    os << "# " << cvm::wrap_string(colvars[0]->name, cvm::cv_width)
+       << "  " << "p(" << cvm::wrap_string(colvars[0]->name, cvm::cv_width-3)
+       << ")\n";
+    size_t igrid;
+    for (igrid = 0; igrid < p.size(); igrid++) {
+      cvm::real const x_grid = (lower_boundary + (igrid+1)*width);
+      os << "  "
+         << std::setprecision(cvm::cv_prec)
+         << std::setw(cvm::cv_width)
+         << x_grid
+         << "  "
+         << std::setprecision(cvm::cv_prec)
+         << std::setw(cvm::cv_width)
+         << p[igrid] << "\n";
+    }
+    os.close();
+  }
+  return os;
+}
+
+
+std::istream & colvarbias_restraint_histogram::read_restart(std::istream &is)
+{
+  return is;
+}
+
+
+std::ostream & colvarbias_restraint_histogram::write_traj_label(std::ostream &os)
+{
+  os << " ";
+  if (b_output_energy) {
+    os << " E_"
+       << cvm::wrap_string(this->name, cvm::en_width-2);
+  }
+  return os;
+}
+
+
+std::ostream & colvarbias_restraint_histogram::write_traj(std::ostream &os)
+{
+  os << " ";
+  if (b_output_energy) {
+    os << " "
+       << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
+       << bias_energy;
+  }
+  return os;
+}
--- a/lib/colvars/colvarbias_restraint.h
+++ b/lib/colvars/colvarbias_restraint.h
@ -168,4 +168,52 @@ protected:
 };


+/// Restrain the 1D histogram of a set of variables (or of a multidimensional one)
+// TODO this could be reimplemented more cleanly as a derived class of both restraint and histogram
+class colvarbias_restraint_histogram : public colvarbias {
+
+public:
+
+  colvarbias_restraint_histogram(char const *key);
+  int init(std::string const &conf);
+  ~colvarbias_restraint_histogram();
+
+  virtual int update();
+
+  virtual std::istream & read_restart(std::istream &is);
+  virtual std::ostream & write_restart(std::ostream &os);
+  virtual std::ostream & write_traj_label(std::ostream &os);
+  virtual std::ostream & write_traj(std::ostream &os);
+
+protected:
+
+  /// Probability density
+  cvm::vector1d<cvm::real> p;
+
+  /// Reference probability density
+  cvm::vector1d<cvm::real> ref_p;
+
+  /// Difference between probability density and reference
+  cvm::vector1d<cvm::real> p_diff;
+
+  /// Lower boundary of the grid
+  cvm::real lower_boundary;
+
+  /// Upper boundary of the grid
+  cvm::real upper_boundary;
+
+  /// Resolution of the grid
+  cvm::real width;
+
+  /// Width of the Gaussians
+  cvm::real gaussian_width;
+
+  /// Restraint force constant
+  cvm::real force_k;
+
+  /// Write the histogram to a file
+  bool b_write_histogram;
+};
+
+
 #endif
--- a/lib/colvars/colvarcomp.cpp
+++ b/lib/colvars/colvarcomp.cpp
@ -54,6 +54,21 @@ colvar::cvc::cvc(std::string const &conf)
 }


+int colvar::cvc::init_total_force_params(std::string const &conf)
+{
+  if (get_keyval_feature(this, conf, "oneSiteSystemForce",
+                         f_cvc_one_site_total_force, is_enabled(f_cvc_one_site_total_force))) {
+    cvm::log("Warning: keyword \"oneSiteSystemForce\" is deprecated: "
+             "please use \"oneSiteTotalForce\" instead.\n");
+  }
+  if (get_keyval_feature(this, conf, "oneSiteTotalForce",
+                         f_cvc_one_site_total_force, is_enabled(f_cvc_one_site_total_force))) {
+    cvm::log("Computing total force on group 1 only");
+  }
+  return COLVARS_OK;
+}
+
+
 cvm::atom_group *colvar::cvc::parse_group(std::string const &conf,
                                          char const *group_key,
                                          bool optional)
@ -77,8 +92,6 @@ cvm::atom_group *colvar::cvc::parse_group(std::string const &conf,

      if (is_enabled(f_cvc_scalable)) {
        cvm::log("Will enable scalable calculation for group \""+group->key+"\".\n");
-      } else {
-        cvm::log("Scalable calculation is not available for group \""+group->key+"\" with the current configuration.\n");
      }
    }

--- a/lib/colvars/colvarcomp.h
+++ b/lib/colvars/colvarcomp.h
@ -108,6 +108,9 @@ public:
                   char const *group_key,
                   bool optional = false);

+  /// \brief Parse options pertaining to total force calculation
+  virtual int init_total_force_params(std::string const &conf);
+
  /// \brief After construction, set data related to dependency handling
  int setup();

@ -306,9 +309,6 @@ protected:
  cvm::rvector     dist_v;
  /// Use absolute positions, ignoring PBCs when present
  bool b_no_PBC;
-  /// Compute total force on first site only to avoid unwanted
-  /// coupling to other colvars (see e.g. Ciccotti et al., 2005)
-  bool b_1site_force;
 public:
  distance(std::string const &conf);
  distance();
@ -388,9 +388,6 @@ protected:
  cvm::atom_group  *ref2;
  /// Use absolute positions, ignoring PBCs when present
  bool b_no_PBC;
-  /// Compute total force on one site only to avoid unwanted
-  /// coupling to other colvars (see e.g. Ciccotti et al., 2005)
-  bool b_1site_force;
  /// Vector on which the distance vector is projected
  cvm::rvector axis;
  /// Norm of the axis
@ -854,6 +851,62 @@ public:
                                  colvarvalue const &x2) const;
 };

+
+/// \brief Colvar component: coordination number between two groups
+/// (colvarvalue::type_scalar type, range [0:N1*N2])
+class colvar::groupcoordnum
+  : public colvar::distance
+{
+protected:
+  /// \brief "Cutoff" for isotropic calculation (default)
+  cvm::real     r0;
+  /// \brief "Cutoff vector" for anisotropic calculation
+  cvm::rvector  r0_vec;
+  /// \brief Wheter dist/r0 or \vec{dist}*\vec{1/r0_vec} should ne be
+  /// used
+  bool b_anisotropic;
+  /// Integer exponent of the function numerator
+  int en;
+  /// Integer exponent of the function denominator
+  int ed;
+public:
+  /// Constructor
+  groupcoordnum(std::string const &conf);
+  groupcoordnum();
+  virtual inline ~groupcoordnum() {}
+  virtual void calc_value();
+  virtual void calc_gradients();
+  virtual void apply_force(colvarvalue const &force);
+  template<bool b_gradients>
+  /// \brief Calculate a coordination number through the function
+  /// (1-x**n)/(1-x**m), x = |A1-A2|/r0 \param r0 "cutoff" for the
+  /// coordination number \param exp_num \i n exponent \param exp_den
+  /// \i m exponent \param A1 atom \param A2 atom
+  static cvm::real switching_function(cvm::real const &r0,
+                                      int const &exp_num, int const &exp_den,
+                                      cvm::atom &A1, cvm::atom &A2);
+
+  /*
+  template<bool b_gradients>
+  /// \brief Calculate a coordination number through the function
+  /// (1-x**n)/(1-x**m), x = |(A1-A2)*(r0_vec)^-|1 \param r0_vec
+  /// vector of different cutoffs in the three directions \param
+  /// exp_num \i n exponent \param exp_den \i m exponent \param A1
+  /// atom \param A2 atom
+  static cvm::real switching_function(cvm::rvector const &r0_vec,
+                                      int const &exp_num, int const &exp_den,
+                                      cvm::atom &A1, cvm::atom &A2);
+
+  virtual cvm::real dist2(colvarvalue const &x1,
+                          colvarvalue const &x2) const;
+  virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
+                                  colvarvalue const &x2) const;
+  virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
+                                  colvarvalue const &x2) const;
+  */
+};
+
+
 /// \brief Colvar component: hydrogen bond, defined as the product of
 /// a colvar::coordnum and 1/2*(1-cos((180-ang)/ang_tol))
 /// (colvarvalue::type_scalar type, range [0:1])
--- a/lib/colvars/colvarcomp_angles.cpp
+++ b/lib/colvars/colvarcomp_angles.cpp
@ -18,9 +18,9 @@ colvar::angle::angle(std::string const &conf)
  group1 = parse_group(conf, "group1");
  group2 = parse_group(conf, "group2");
  group3 = parse_group(conf, "group3");
-  if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {
-    cvm::log("Computing total force on group 1 only");
-  }
+
+  init_total_force_params(conf);
+
  x.type(colvarvalue::type_scalar);
 }

@ -33,7 +33,6 @@ colvar::angle::angle(cvm::atom const &a1,
  provide(f_cvc_inv_gradient);
  provide(f_cvc_Jacobian);
  provide(f_cvc_com_based);
-  b_1site_force = false;

  group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
  group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
@ -94,7 +93,7 @@ void colvar::angle::calc_force_invgrads()
  // centered on group2, which means group2 is kept fixed
  // when propagating changes in the angle)

-  if (b_1site_force) {
+  if (is_enabled(f_cvc_one_site_total_force)) {
    group1->read_total_forces();
    cvm::real norm_fact = 1.0 / dxdr1.norm2();
    ft.real_value = norm_fact * dxdr1 * group1->total_force();
@ -140,9 +139,8 @@ colvar::dipole_angle::dipole_angle(std::string const &conf)
  group2 = parse_group(conf, "group2");
  group3 = parse_group(conf, "group3");

-  if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {
-    cvm::log("Computing total force on group 1 only");
-  }
+  init_total_force_params(conf);
+
  x.type(colvarvalue::type_scalar);
 }

@ -152,7 +150,6 @@ colvar::dipole_angle::dipole_angle(cvm::atom const &a1,
                      cvm::atom const &a3)
 {
  function_type = "dipole_angle";
-  b_1site_force = false;

  group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
  group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
@ -250,14 +247,13 @@ colvar::dihedral::dihedral(std::string const &conf)
  provide(f_cvc_Jacobian);
  provide(f_cvc_com_based);

-  if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {
-    cvm::log("Computing total force on group 1 only");
-  }
  group1 = parse_group(conf, "group1");
  group2 = parse_group(conf, "group2");
  group3 = parse_group(conf, "group3");
  group4 = parse_group(conf, "group4");

+  init_total_force_params(conf);
+
  x.type(colvarvalue::type_scalar);
 }

@ -422,7 +418,7 @@ void colvar::dihedral::calc_force_invgrads()
  cvm::real const fact4 = d34 * std::sqrt(1.0 - dot4 * dot4);

  group1->read_total_forces();
-  if ( b_1site_force ) {
+  if (is_enabled(f_cvc_one_site_total_force)) {
    // This is only measuring the force on group 1
    ft.real_value = PI/180.0 * fact1 * (cross1 * group1->total_force());
  } else {
--- a/lib/colvars/colvarcomp_coordnums.cpp
+++ b/lib/colvars/colvarcomp_coordnums.cpp
@ -338,3 +338,151 @@ void colvar::selfcoordnum::apply_force(colvarvalue const &force)
  }
 }

+
+// groupcoordnum member functions
+colvar::groupcoordnum::groupcoordnum(std::string const &conf)
+  : distance(conf), b_anisotropic(false)
+{
+  function_type = "groupcoordnum";
+  x.type(colvarvalue::type_scalar);
+
+  // group1 and group2 are already initialized by distance()
+  if (group1->b_dummy || group2->b_dummy)
+    cvm::fatal_error("Error: neither group can be a dummy atom\n");
+
+  bool const b_scale = get_keyval(conf, "cutoff", r0,
+                                  cvm::real(4.0 * cvm::unit_angstrom()));
+
+  if (get_keyval(conf, "cutoff3", r0_vec,
+                 cvm::rvector(4.0, 4.0, 4.0), parse_silent)) {
+
+    if (b_scale)
+      cvm::fatal_error("Error: cannot specify \"scale\" and "
+                       "\"scale3\" at the same time.\n");
+    b_anisotropic = true;
+    // remove meaningless negative signs
+    if (r0_vec.x < 0.0) r0_vec.x *= -1.0;
+    if (r0_vec.y < 0.0) r0_vec.y *= -1.0;
+    if (r0_vec.z < 0.0) r0_vec.z *= -1.0;
+  }
+
+  get_keyval(conf, "expNumer", en, int(6) );
+  get_keyval(conf, "expDenom", ed, int(12));
+
+  if ( (en%2) || (ed%2) ) {
+    cvm::fatal_error("Error: odd exponents provided, can only use even ones.\n");
+  }
+
+}
+
+
+colvar::groupcoordnum::groupcoordnum()
+  : b_anisotropic(false)
+{
+  function_type = "groupcoordnum";
+  x.type(colvarvalue::type_scalar);
+}
+
+
+template<bool calculate_gradients>
+cvm::real colvar::groupcoordnum::switching_function(cvm::real const &r0,
+                                                    int const &en,
+                                                    int const &ed,
+                                                    cvm::atom &A1,
+                                                    cvm::atom &A2)
+{
+  cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
+  cvm::real const l2 = diff.norm2()/(r0*r0);
+
+  // Assume en and ed are even integers, and avoid sqrt in the following
+  int const en2 = en/2;
+  int const ed2 = ed/2;
+
+  cvm::real const xn = std::pow(l2, en2);
+  cvm::real const xd = std::pow(l2, ed2);
+  cvm::real const func = (1.0-xn)/(1.0-xd);
+
+  if (calculate_gradients) {
+    cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
+    cvm::rvector const dl2dx = (2.0/(r0*r0))*diff;
+    A1.grad += (-1.0)*dFdl2*dl2dx;
+    A2.grad +=        dFdl2*dl2dx;
+  }
+
+  return func;
+}
+
+
+#if 0  // AMG: I don't think there's any reason to support anisotropic,
+       //      and I don't have those flags below in calc_value, but
+       //      if I need them, I'll also need to uncomment this method
+template<bool calculate_gradients>
+cvm::real colvar::groupcoordnum::switching_function(cvm::rvector const &r0_vec,
+                                                    int const &en,
+                                                    int const &ed,
+                                                    cvm::atom &A1,
+                                                    cvm::atom &A2)
+{
+  cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
+  cvm::rvector const scal_diff(diff.x/r0_vec.x, diff.y/r0_vec.y, diff.z/r0_vec.z);
+  cvm::real const l2 = scal_diff.norm2();
+
+  // Assume en and ed are even integers, and avoid sqrt in the following
+  int const en2 = en/2;
+  int const ed2 = ed/2;
+
+  cvm::real const xn = std::pow(l2, en2);
+  cvm::real const xd = std::pow(l2, ed2);
+  cvm::real const func = (1.0-xn)/(1.0-xd);
+
+  if (calculate_gradients) {
+    cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
+    cvm::rvector const dl2dx((2.0/(r0_vec.x*r0_vec.x))*diff.x,
+                             (2.0/(r0_vec.y*r0_vec.y))*diff.y,
+                             (2.0/(r0_vec.z*r0_vec.z))*diff.z);
+    A1.grad += (-1.0)*dFdl2*dl2dx;
+    A2.grad +=        dFdl2*dl2dx;
+  }
+  return func;
+}
+#endif
+
+
+
+void colvar::groupcoordnum::calc_value()
+{
+
+  // create fake atoms to hold the com coordinates
+  cvm::atom group1_com_atom;
+  cvm::atom group2_com_atom;
+  group1_com_atom.pos = group1->center_of_mass();
+  group2_com_atom.pos = group2->center_of_mass();
+
+  x.real_value = coordnum::switching_function<false>(r0, en, ed,
+                                                     group1_com_atom, group2_com_atom);
+
+}
+
+
+void colvar::groupcoordnum::calc_gradients()
+{
+  cvm::atom group1_com_atom;
+  cvm::atom group2_com_atom;
+  group1_com_atom.pos = group1->center_of_mass();
+  group2_com_atom.pos = group2->center_of_mass();
+
+  coordnum::switching_function<true>(r0, en, ed, group1_com_atom, group2_com_atom);
+  group1->set_weighted_gradient(group1_com_atom.grad);
+  group2->set_weighted_gradient(group2_com_atom.grad);
+
+}
+
+
+void colvar::groupcoordnum::apply_force(colvarvalue const &force)
+{
+  if (!group1->noforce)
+    group1->apply_colvar_force(force.real_value);
+
+  if (!group2->noforce)
+    group2->apply_colvar_force(force.real_value);
+}
--- a/lib/colvars/colvarcomp_distances.cpp
+++ b/lib/colvars/colvarcomp_distances.cpp
@ -18,14 +18,14 @@ colvar::distance::distance(std::string const &conf)
  provide(f_cvc_Jacobian);
  provide(f_cvc_com_based);

+  group1 = parse_group(conf, "group1");
+  group2 = parse_group(conf, "group2");
+
  if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) {
    cvm::log("Computing distance using absolute positions (not minimal-image)");
  }
-  if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {
-    cvm::log("Computing total force on group 1 only");
-  }
-  group1 = parse_group(conf, "group1");
-  group2 = parse_group(conf, "group2");
+
+  init_total_force_params(conf);

  x.type(colvarvalue::type_scalar);
 }
@ -38,7 +38,6 @@ colvar::distance::distance()
  provide(f_cvc_inv_gradient);
  provide(f_cvc_Jacobian);
  provide(f_cvc_com_based);
-  b_1site_force = false;
  b_no_PBC = false;
  x.type(colvarvalue::type_scalar);
 }
@ -67,7 +66,7 @@ void colvar::distance::calc_gradients()
 void colvar::distance::calc_force_invgrads()
 {
  group1->read_total_forces();
-  if ( b_1site_force ) {
+  if (is_enabled(f_cvc_one_site_total_force)) {
    ft.real_value = -1.0 * (group1->total_force() * dist_v.unit());
  } else {
    group2->read_total_forces();
@ -97,6 +96,7 @@ colvar::distance_vec::distance_vec(std::string const &conf)
  : distance(conf)
 {
  function_type = "distance_vec";
+  provide(f_cvc_com_based);
  x.type(colvarvalue::type_3vector);
 }

@ -105,6 +105,7 @@ colvar::distance_vec::distance_vec()
  : distance()
 {
  function_type = "distance_vec";
+  provide(f_cvc_com_based);
  x.type(colvarvalue::type_3vector);
 }

@ -185,9 +186,9 @@ colvar::distance_z::distance_z(std::string const &conf)
  if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) {
    cvm::log("Computing distance using absolute positions (not minimal-image)");
  }
-  if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {
-    cvm::log("Computing total force on group \"main\" only");
-  }
+
+  init_total_force_params(conf);
+
 }

 colvar::distance_z::distance_z()
@ -251,7 +252,7 @@ void colvar::distance_z::calc_force_invgrads()
 {
  main->read_total_forces();

-  if (fixed_axis && !b_1site_force) {
+  if (fixed_axis && !is_enabled(f_cvc_one_site_total_force)) {
    ref1->read_total_forces();
    ft.real_value = 0.5 * ((main->total_force() - ref1->total_force()) * axis);
  } else {
@ -351,7 +352,7 @@ void colvar::distance_xy::calc_force_invgrads()
 {
  main->read_total_forces();

-  if (fixed_axis && !b_1site_force) {
+  if (fixed_axis && !is_enabled(f_cvc_one_site_total_force)) {
    ref1->read_total_forces();
    ft.real_value = 0.5 / x.real_value * ((main->total_force() - ref1->total_force()) * dist_v_ortho);
  } else {
@ -382,6 +383,7 @@ colvar::distance_dir::distance_dir(std::string const &conf)
  : distance(conf)
 {
  function_type = "distance_dir";
+  provide(f_cvc_com_based);
  x.type(colvarvalue::type_unit3vector);
 }

@ -390,6 +392,7 @@ colvar::distance_dir::distance_dir()
  : distance()
 {
  function_type = "distance_dir";
+  provide(f_cvc_com_based);
  x.type(colvarvalue::type_unit3vector);
 }

@ -461,7 +464,6 @@ colvar::distance_inv::distance_inv()
 {
  function_type = "distance_inv";
  exponent = 6;
-  b_1site_force = false;
  x.type(colvarvalue::type_scalar);
 }

--- a/lib/colvars/colvardeps.cpp
+++ b/lib/colvars/colvardeps.cpp
@ -293,6 +293,9 @@ void colvardeps::init_cv_requires() {
    f_description(f_cv_output_total_force, "output total force");
    f_req_self(f_cv_output_total_force, f_cv_total_force);

+    f_description(f_cv_subtract_applied_force, "subtract applied force from total force");
+    f_req_self(f_cv_subtract_applied_force, f_cv_total_force);
+
    f_description(f_cv_lower_boundary, "lower boundary");
    f_req_self(f_cv_lower_boundary, f_cv_scalar);

@ -376,6 +379,11 @@ void colvardeps::init_cvc_requires() {

    f_description(f_cvc_com_based, "depends on group centers of mass");

+    // Compute total force on first site only to avoid unwanted
+    // coupling to other colvars (see e.g. Ciccotti et al., 2005)
+    f_description(f_cvc_one_site_total_force, "compute total collective force only from one group center");
+    f_req_self(f_cvc_one_site_total_force, f_cvc_com_based);
+
    f_description(f_cvc_scalable, "scalable calculation");
    f_req_self(f_cvc_scalable, f_cvc_scalable_com);

--- a/lib/colvars/colvardeps.h
+++ b/lib/colvars/colvardeps.h
@ -176,6 +176,8 @@ public:
    f_cv_total_force,
    /// \brief Calculate total force from atomic forces
    f_cv_total_force_calc,
+    /// \brief Subtract the applied force from the total force
+    f_cv_subtract_applied_force,
    /// \brief Estimate Jacobian derivative
    f_cv_Jacobian,
    /// \brief Do not report the Jacobian force as part of the total force
@ -236,6 +238,7 @@ public:
    /// \brief If enabled, calc_gradients() will call debug_gradients() for every group needed
    f_cvc_debug_gradient,
    f_cvc_Jacobian,
+    f_cvc_one_site_total_force,
    f_cvc_com_based,
    f_cvc_scalable,
    f_cvc_scalable_com,
--- a/lib/colvars/colvargrid.h
+++ b/lib/colvars/colvargrid.h
@ -382,8 +382,8 @@ public:
  inline int current_bin_scalar(int const i, int const iv) const
  {
    return value_to_bin_scalar(actual_value[i] ?
-        cv[i]->actual_value().vector1d_value[iv] :
-        cv[i]->value().vector1d_value[iv], i);
+                               cv[i]->actual_value().vector1d_value[iv] :
+                               cv[i]->value().vector1d_value[iv], i);
  }

  /// \brief Use the lower boundary and the width to report which bin
@ -395,8 +395,8 @@ public:

  /// \brief Same as the standard version, but uses another grid definition
  inline int value_to_bin_scalar(colvarvalue const &value,
-                                  colvarvalue const &new_offset,
-                                  cvm::real   const &new_width) const
+                                 colvarvalue const &new_offset,
+                                 cvm::real   const &new_width) const
  {
    return (int) std::floor( (value.real_value - new_offset.real_value) / new_width );
  }
@ -410,22 +410,22 @@ public:

  /// \brief Same as the standard version, but uses different parameters
  inline colvarvalue bin_to_value_scalar(int const &i_bin,
-                                          colvarvalue const &new_offset,
-                                          cvm::real const &new_width) const
+                                         colvarvalue const &new_offset,
+                                         cvm::real const &new_width) const
  {
    return new_offset.real_value + new_width * (0.5 + i_bin);
  }

  /// Set the value at the point with index ix
  inline void set_value(std::vector<int> const &ix,
-                         T const &t,
-                         size_t const &imult = 0)
+                        T const &t,
+                        size_t const &imult = 0)
  {
    data[this->address(ix)+imult] = t;
    has_data = true;
  }

- /// \brief Get the change from this to other_grid
+  /// \brief Get the change from this to other_grid
  /// and store the result in this.
  /// this_grid := other_grid - this_grid
  /// Grids must have the same dimensions.
@ -434,13 +434,13 @@ public:

    if (other_grid.multiplicity() != this->multiplicity()) {
      cvm::error("Error: trying to subtract two grids with "
-                        "different multiplicity.\n");
+                 "different multiplicity.\n");
      return;
    }

    if (other_grid.data.size() != this->data.size()) {
      cvm::error("Error: trying to subtract two grids with "
-                        "different size.\n");
+                 "different size.\n");
      return;
    }

@ -457,13 +457,13 @@ public:
  {
    if (other_grid.multiplicity() != this->multiplicity()) {
      cvm::error("Error: trying to copy two grids with "
-                        "different multiplicity.\n");
+                 "different multiplicity.\n");
      return;
    }

    if (other_grid.data.size() != this->data.size()) {
      cvm::error("Error: trying to copy two grids with "
-                        "different size.\n");
+                 "different size.\n");
      return;
    }

@ -493,7 +493,7 @@ public:
  /// \brief Get the binned value indexed by ix, or the first of them
  /// if the multiplicity is larger than 1
  inline T const & value(std::vector<int> const &ix,
-                          size_t const &imult = 0) const
+                         size_t const &imult = 0) const
  {
    return data[this->address(ix) + imult];
  }
@ -541,7 +541,7 @@ public:
  /// boundaries; a negative number is returned if the given point is
  /// off-grid
  inline cvm::real bin_distance_from_boundaries(std::vector<colvarvalue> const &values,
-                                                 bool skip_hard_boundaries = false)
+                                                bool skip_hard_boundaries = false)
  {
    cvm::real minimum = 1.0E+16;
    for (size_t i = 0; i < nd; i++) {
@ -574,7 +574,7 @@ public:
  {
    if (other_grid.multiplicity() != this->multiplicity()) {
      cvm::error("Error: trying to merge two grids with values of "
-                        "different multiplicity.\n");
+                 "different multiplicity.\n");
      return;
    }

@ -593,8 +593,8 @@ public:
      for (size_t i = 0; i < nd; i++) {
        oix[i] =
          value_to_bin_scalar(bin_to_value_scalar(ix[i], gb[i], gw[i]),
-                               ogb[i],
-                               ogw[i]);
+                              ogb[i],
+                              ogw[i]);
      }

      if (! other_grid.index_ok(oix)) {
@ -614,11 +614,11 @@ public:
  /// \brief Add data from another grid of the same type, AND
  /// identical definition (boundaries, widths)
  void add_grid(colvar_grid<T> const &other_grid,
-                 cvm::real scale_factor = 1.0)
+                cvm::real scale_factor = 1.0)
  {
    if (other_grid.multiplicity() != this->multiplicity()) {
      cvm::error("Error: trying to sum togetehr two grids with values of "
-                        "different multiplicity.\n");
+                 "different multiplicity.\n");
      return;
    }
    if (scale_factor != 1.0)
@ -636,7 +636,7 @@ public:
  /// \brief Return the value suitable for output purposes (so that it
  /// may be rescaled or manipulated without changing it permanently)
  virtual inline T value_output(std::vector<int> const &ix,
-                                 size_t const &imult = 0)
+                                size_t const &imult = 0)
  {
    return value(ix, imult);
  }
@ -645,9 +645,9 @@ public:
  /// into the internal representation (the two may be different,
  /// e.g. when using colvar_grid_count)
  virtual inline void value_input(std::vector<int> const &ix,
-                                   T const &t,
-                                   size_t const &imult = 0,
-                                   bool add = false)
+                                  T const &t,
+                                  size_t const &imult = 0,
+                                  bool add = false)
  {
    if ( add )
      data[address(ix) + imult] += t;
@ -737,7 +737,8 @@ public:
  }

  /// Read a grid definition from a config string
-  int parse_params(std::string const &conf)
+  int parse_params(std::string const &conf,
+                   colvarparse::Parse_Mode const parse_mode = colvarparse::parse_normal)
  {
    if (cvm::debug()) cvm::log("Reading grid configuration from string.\n");

@ -746,30 +747,33 @@ public:

    {
      size_t nd_in = 0;
+      // this is only used in state files
      colvarparse::get_keyval(conf, "n_colvars", nd_in, nd, colvarparse::parse_silent);
      if (nd_in != nd) {
        cvm::error("Error: trying to read data for a grid "
-                    "that contains a different number of colvars ("+
-                    cvm::to_str(nd_in)+") than the grid defined "
-                    "in the configuration file("+cvm::to_str(nd)+
-                    ").\n");
+                   "that contains a different number of colvars ("+
+                   cvm::to_str(nd_in)+") than the grid defined "
+                   "in the configuration file("+cvm::to_str(nd)+
+                   ").\n");
        return COLVARS_ERROR;
      }
    }

+    // underscore keywords are used in state file
    colvarparse::get_keyval(conf, "lower_boundaries",
-                             lower_boundaries, lower_boundaries, colvarparse::parse_silent);
+                            lower_boundaries, lower_boundaries, colvarparse::parse_silent);
    colvarparse::get_keyval(conf, "upper_boundaries",
-                             upper_boundaries, upper_boundaries, colvarparse::parse_silent);
+                            upper_boundaries, upper_boundaries, colvarparse::parse_silent);

-    // support also camel case
+    // camel case keywords are used in config file
    colvarparse::get_keyval(conf, "lowerBoundaries",
-                             lower_boundaries, lower_boundaries, colvarparse::parse_silent);
+                            lower_boundaries, lower_boundaries, parse_mode);
    colvarparse::get_keyval(conf, "upperBoundaries",
-                             upper_boundaries, upper_boundaries, colvarparse::parse_silent);
+                            upper_boundaries, upper_boundaries, parse_mode);

-    colvarparse::get_keyval(conf, "widths", widths, widths, colvarparse::parse_silent);
+    colvarparse::get_keyval(conf, "widths", widths, widths, parse_mode);

+    // only used in state file
    colvarparse::get_keyval(conf, "sizes", nx, nx, colvarparse::parse_silent);

    if (nd < lower_boundaries.size()) nd = lower_boundaries.size();
@ -808,13 +812,13 @@ public:
  {
    for (size_t i = 0; i < nd; i++) {
      if ( (std::sqrt(cv[i]->dist2(cv[i]->lower_boundary,
-                                     lower_boundaries[i])) > 1.0E-10) ||
+                                   lower_boundaries[i])) > 1.0E-10) ||
           (std::sqrt(cv[i]->dist2(cv[i]->upper_boundary,
-                                     upper_boundaries[i])) > 1.0E-10) ||
+                                   upper_boundaries[i])) > 1.0E-10) ||
           (std::sqrt(cv[i]->dist2(cv[i]->width,
-                                     widths[i])) > 1.0E-10) ) {
+                                   widths[i])) > 1.0E-10) ) {
        cvm::error("Error: restart information for a grid is "
-                    "inconsistent with that of its colvars.\n");
+                   "inconsistent with that of its colvars.\n");
        return;
      }
    }
@ -830,19 +834,19 @@ public:
      // matter: boundaries should be EXACTLY the same (otherwise,
      // map_grid() should be used)
      if ( (std::fabs(other_grid.lower_boundaries[i] -
-                       lower_boundaries[i]) > 1.0E-10) ||
+                      lower_boundaries[i]) > 1.0E-10) ||
           (std::fabs(other_grid.upper_boundaries[i] -
-                       upper_boundaries[i]) > 1.0E-10) ||
+                      upper_boundaries[i]) > 1.0E-10) ||
           (std::fabs(other_grid.widths[i] -
-                       widths[i]) > 1.0E-10) ||
+                      widths[i]) > 1.0E-10) ||
           (data.size() != other_grid.data.size()) ) {
-      cvm::error("Error: inconsistency between "
-                  "two grids that are supposed to be equal, "
-                  "aside from the data stored.\n");
-      return;
+        cvm::error("Error: inconsistency between "
+                   "two grids that are supposed to be equal, "
+                   "aside from the data stored.\n");
+        return;
+      }
    }
  }
-}


  /// \brief Read grid entry in restart file
@ -853,7 +857,7 @@ public:
    if ((is >> key) && (key == std::string("grid_parameters"))) {
      is.seekg(start_pos, std::ios::beg);
      is >> colvarparse::read_block("grid_parameters", conf);
-      parse_params(conf);
+      parse_params(conf, colvarparse::parse_silent);
    } else {
      cvm::log("Grid parameters are missing in the restart file, using those from the configuration.\n");
      is.seekg(start_pos, std::ios::beg);
@ -871,11 +875,11 @@ public:
  }


-/// \brief Write the grid data without labels, as they are
-/// represented in memory
-/// \param buf_size Number of values per line
+  /// \brief Write the grid data without labels, as they are
+  /// represented in memory
+  /// \param buf_size Number of values per line
  std::ostream & write_raw(std::ostream &os,
-                            size_t const buf_size = 3)
+                           size_t const buf_size = 3)
  {
    std::streamsize const w = os.width();
    std::streamsize const p = os.precision();
@ -935,10 +939,10 @@ public:
    os << std::setw(2) << "# " << nd << "\n";
    for (size_t i = 0; i < nd; i++) {
      os << "# "
-        << std::setw(10) << lower_boundaries[i]
-        << std::setw(10) << widths[i]
-        << std::setw(10) << nx[i] << "  "
-        << periodic[i] << "\n";
+         << std::setw(10) << lower_boundaries[i]
+         << std::setw(10) << widths[i]
+         << std::setw(10) << nx[i] << "  "
+         << periodic[i] << "\n";
    }


@ -951,14 +955,14 @@ public:

      for (size_t i = 0; i < nd; i++) {
        os << " "
-          << std::setw(w) << std::setprecision(p)
-          << bin_to_value_scalar(ix[i], i);
+           << std::setw(w) << std::setprecision(p)
+           << bin_to_value_scalar(ix[i], i);
      }
      os << " ";
      for (size_t imult = 0; imult < mult; imult++) {
        os << " "
-          << std::setw(w) << std::setprecision(p)
-          << value_output(ix, imult);
+           << std::setw(w) << std::setprecision(p)
+           << value_output(ix, imult);
      }
      os << "\n";
    }
@ -986,7 +990,7 @@ public:

    if ( !(is >> hash) || (hash != "#") ) {
      cvm::error("Error reading grid at position "+
-                  cvm::to_str(is.tellg())+" in stream(read \"" + hash + "\")\n");
+                 cvm::to_str(is.tellg())+" in stream(read \"" + hash + "\")\n");
      return is;
    }

@ -1008,7 +1012,7 @@ public:
    for (size_t i = 0; i < nd; i++ ) {
      if ( !(is >> hash) || (hash != "#") ) {
        cvm::error("Error reading grid at position "+
-                    cvm::to_str(is.tellg())+" in stream(read \"" + hash + "\")\n");
+                   cvm::to_str(is.tellg())+" in stream(read \"" + hash + "\")\n");
        return is;
      }

@ -1016,10 +1020,10 @@ public:


      if ( (std::fabs(lower - lower_boundaries[i].real_value) > 1.0e-10) ||
-          (std::fabs(width - widths[i] ) > 1.0e-10) ||
-          (nx_read[i] != nx[i]) ) {
+           (std::fabs(width - widths[i] ) > 1.0e-10) ||
+           (nx_read[i] != nx[i]) ) {
        cvm::log("Warning: reading from different grid definition (colvar "
-                  + cvm::to_str(i+1) + "); remapping data on new grid.\n");
+                 + cvm::to_str(i+1) + "); remapping data on new grid.\n");
        remap = true;
      }
    }
@ -1063,7 +1067,6 @@ public:

  /// \brief Write the grid data without labels, as they are
  /// represented in memory
-  /// \param buf_size Number of values per line
  std::ostream & write_opendx(std::ostream &os)
  {
    // write the header
@ -1122,11 +1125,11 @@ public:

  /// Constructor
  colvar_grid_count(std::vector<int> const &nx_i,
-                     size_t const           &def_count = 0);
+                    size_t const           &def_count = 0);

  /// Constructor from a vector of colvars
  colvar_grid_count(std::vector<colvar *>  &colvars,
-                     size_t const           &def_count = 0);
+                    size_t const           &def_count = 0);

  /// Increment the counter at given position
  inline void incr_count(std::vector<int> const &ix)
@ -1136,7 +1139,7 @@ public:

  /// \brief Get the binned count indexed by ix from the newly read data
  inline size_t const & new_count(std::vector<int> const &ix,
-                                   size_t const &imult = 0)
+                                  size_t const &imult = 0)
  {
    return new_data[address(ix) + imult];
  }
@ -1145,9 +1148,9 @@ public:
  /// into the internal representation (it may have been rescaled or
  /// manipulated)
  virtual inline void value_input(std::vector<int> const &ix,
-                                   size_t const &t,
-                                   size_t const &imult = 0,
-                                   bool add = false)
+                                  size_t const &t,
+                                  size_t const &imult = 0,
+                                  bool add = false)
  {
    if (add) {
      data[address(ix)] += t;
@ -1164,7 +1167,7 @@ public:
  /// \brief Return the log-gradient from finite differences
  /// on the *same* grid for dimension n
  inline const cvm::real log_gradient_finite_diff( const std::vector<int> &ix0,
-    int n = 0)
+                                                   int n = 0)
  {
    cvm::real A0, A1;
    std::vector<int> ix;
@ -1377,7 +1380,7 @@ public:
  /// \brief Return the value of the function at ix divided by its
  /// number of samples (if the count grid is defined)
  virtual inline cvm::real value_output(std::vector<int> const &ix,
-                                         size_t const &imult = 0)
+                                        size_t const &imult = 0)
  {
    if (samples)
      return (samples->value(ix) > 0) ?
@ -1391,9 +1394,9 @@ public:
  /// into the internal representation (it may have been rescaled or
  /// manipulated)
  virtual inline void value_input(std::vector<int> const &ix,
-                                   cvm::real const &new_value,
-                                   size_t const &imult = 0,
-                                   bool add = false)
+                                  cvm::real const &new_value,
+                                  size_t const &imult = 0,
+                                  bool add = false)
  {
    if (add) {
      if (samples)
--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@ -293,6 +293,9 @@ int colvarmodule::parse_biases(std::string const &conf)
  /// initialize histograms
  parse_biases_type<colvarbias_histogram>(conf, "histogram", n_histo_biases);

+  /// initialize histogram restraints
+  parse_biases_type<colvarbias_restraint_histogram>(conf, "histogramRestraint", n_rest_biases);
+
  /// initialize linear restraints
  parse_biases_type<colvarbias_restraint_linear>(conf, "linear", n_rest_biases);

--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@ -4,7 +4,7 @@
 #define COLVARMODULE_H

 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2016-09-14"
+#define COLVARS_VERSION "2016-09-30"
 #endif

 #ifndef COLVARS_DEBUG
--- a/lib/colvars/colvarscript.cpp
+++ b/lib/colvars/colvarscript.cpp
@ -243,11 +243,17 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
  }

  if (subcmd == "getappliedforce") {
-    result = (cv->bias_force()).to_simple_string();
+    result = (cv->applied_force()).to_simple_string();
    return COLVARS_OK;
  }

  if (subcmd == "getsystemforce") {
+    // TODO warning here
+    result = (cv->total_force()).to_simple_string();
+    return COLVARS_OK;
+  }
+
+  if (subcmd == "gettotalforce") {
    result = (cv->total_force()).to_simple_string();
    return COLVARS_OK;
  }
--- a/lib/colvars/colvartypes.h
+++ b/lib/colvars/colvartypes.h
@ -57,6 +57,12 @@ public:
    }
  }

+  /// Return a reference to the data
+  inline std::vector<T> &data_array()
+  {
+    return data;
+  }
+
  inline ~vector1d()
  {
    data.clear();
@ -203,6 +209,16 @@ public:
    return std::sqrt(this->norm2());
  }

+  inline cvm::real sum() const
+  {
+    cvm::real result = 0.0;
+    size_t i;
+    for (i = 0; i < this->size(); i++) {
+      result += (*this)[i];
+    }
+    return result;
+  }
+
  /// Slicing
  inline vector1d<T> const slice(size_t const i1, size_t const i2) const
  {
@ -295,11 +311,23 @@ public:
  {
    std::stringstream stream(s);
    size_t i = 0;
-    while ((stream >> (*this)[i]) && (i < this->size())) {
-      i++;
-    }
-    if (i < this->size()) {
-      return COLVARS_ERROR;
+    if (this->size()) {
+      while ((stream >> (*this)[i]) && (i < this->size())) {
+        i++;
+      }
+      if (i < this->size()) {
+        return COLVARS_ERROR;
+      }
+    } else {
+      T input;
+      while (stream >> input) {
+        if ((i % 100) == 0) {
+          data.reserve(data.size()+100);
+        }
+        data.resize(data.size()+1);
+        data[i] = input;
+        i++;
+      }
    }
    return COLVARS_OK;
  }
@ -434,6 +462,12 @@ public:
    this->clear();
  }

+  /// Return a reference to the data
+  inline std::vector<T> &data_array()
+  {
+    return data;
+  }
+
  inline row & operator [] (size_t const i)
  {
    return rows[i];
--- a/lib/kokkos/config/configure_compton_cpu.sh
+++ b/lib/kokkos/config/configure_compton_cpu.sh
--- a/lib/kokkos/config/configure_compton_mic.sh
+++ b/lib/kokkos/config/configure_compton_mic.sh
--- a/lib/kokkos/config/configure_kokkos.sh
+++ b/lib/kokkos/config/configure_kokkos.sh
--- a/lib/kokkos/config/configure_kokkos_nvidia.sh
+++ b/lib/kokkos/config/configure_kokkos_nvidia.sh
--- a/lib/kokkos/config/configure_shannon.sh
+++ b/lib/kokkos/config/configure_shannon.sh
--- 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/potentials/charmm22.cmap
+++ b/potentials/charmm22.cmap
--- a/potentials/charmm36.cmap
+++ b/potentials/charmm36.cmap
--- 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/.gitignore
+++ b/src/.gitignore
@ -303,6 +303,8 @@
 /fix_bond_create.h
 /fix_bond_swap.cpp
 /fix_bond_swap.h
+/fix_cmap.cpp
+/fix_cmap.h
 /fix_deposit.cpp
 /fix_deposit.h
 /fix_efield.cpp
--- a/src/KOKKOS/comm_tiled_kokkos.cpp
+++ b/src/KOKKOS/comm_tiled_kokkos.cpp
@ -0,0 +1,266 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include <string.h>
+#include "comm_tiled_kokkos.h"
+#include "comm_brick.h"
+#include "atom_kokkos.h"
+#include "atom_vec.h"
+#include "domain.h"
+#include "force.h"
+#include "pair.h"
+#include "neighbor.h"
+#include "modify.h"
+#include "fix.h"
+#include "compute.h"
+#include "output.h"
+#include "dump.h"
+#include "memory.h"
+#include "error.h"
+#include "atom_masks.h"
+
+using namespace LAMMPS_NS;
+
+#define BUFFACTOR 1.5
+#define BUFFACTOR 1.5
+#define BUFMIN 1000
+#define BUFEXTRA 1000
+#define EPSILON 1.0e-6
+
+#define DELTA_PROCS 16
+
+enum{SINGLE,MULTI};               // same as in Comm
+enum{LAYOUT_UNIFORM,LAYOUT_NONUNIFORM,LAYOUT_TILED};    // several files
+
+/* ---------------------------------------------------------------------- */
+
+CommTiledKokkos::CommTiledKokkos(LAMMPS *lmp) : CommTiled(lmp)
+{
+
+}
+
+/* ---------------------------------------------------------------------- */
+//IMPORTANT: we *MUST* pass "*oldcomm" to the Comm initializer here, as
+//           the code below *requires* that the (implicit) copy constructor
+//           for Comm is run and thus creating a shallow copy of "oldcomm".
+//           The call to Comm::copy_arrays() then converts the shallow copy
+//           into a deep copy of the class with the new layout.
+
+CommTiledKokkos::CommTiledKokkos(LAMMPS *lmp, Comm *oldcomm) : CommTiled(lmp,oldcomm)
+{
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+CommTiledKokkos::~CommTiledKokkos()
+{
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+
+/* ----------------------------------------------------------------------
+   forward communication of atom coords every timestep
+   other per-atom attributes may also be sent via pack/unpack routines
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::forward_comm(int dummy)
+{
+  if (comm_x_only) {
+    atomKK->sync(Host,X_MASK);
+    atomKK->modified(Host,X_MASK);
+  } else if (ghost_velocity) {
+    atomKK->sync(Host,X_MASK | V_MASK);
+    atomKK->modified(Host,X_MASK | V_MASK);
+  } else {
+    atomKK->sync(Host,ALL_MASK);
+    atomKK->modified(Host,ALL_MASK);
+  }
+
+  CommTiled::forward_comm(dummy);
+}
+
+/* ----------------------------------------------------------------------
+   reverse communication of forces on atoms every timestep
+   other per-atom attributes may also be sent via pack/unpack routines
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::reverse_comm()
+{
+  if (comm_f_only)
+    atomKK->sync(Host,F_MASK);
+  else
+    atomKK->sync(Host,ALL_MASK);
+  CommTiled::reverse_comm();
+  if (comm_f_only)
+    atomKK->modified(Host,F_MASK);
+  else
+    atomKK->modified(Host,ALL_MASK);
+  atomKK->sync(Device,ALL_MASK);
+}
+
+/* ----------------------------------------------------------------------
+   exchange: move atoms to correct processors
+   atoms exchanged with procs that touch sub-box in each of 3 dims
+   send out atoms that have left my box, receive ones entering my box
+   atoms will be lost if not inside a touching proc's box
+     can happen if atom moves outside of non-periodic bounary
+     or if atom moves more than one proc away
+   this routine called before every reneighboring
+   for triclinic, atoms must be in lamda coords (0-1) before exchange is called
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::exchange()
+{
+  atomKK->sync(Host,ALL_MASK);
+  CommTiled::exchange();
+  atomKK->modified(Host,ALL_MASK);
+}
+
+/* ----------------------------------------------------------------------
+   borders: list nearby atoms to send to neighboring procs at every timestep
+   one list is created per swap/proc that will be made
+   as list is made, actually do communication
+   this does equivalent of a forward_comm(), so don't need to explicitly
+     call forward_comm() on reneighboring timestep
+   this routine is called before every reneighboring
+   for triclinic, atoms must be in lamda coords (0-1) before borders is called
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::borders()
+{
+  atomKK->sync(Host,ALL_MASK);
+  CommTiled::borders();
+  atomKK->modified(Host,ALL_MASK);
+}
+
+/* ----------------------------------------------------------------------
+   forward communication invoked by a Pair
+   nsize used only to set recv buffer limit
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::forward_comm_pair(Pair *pair)
+{
+  CommTiled::forward_comm_pair(pair);
+}
+
+/* ----------------------------------------------------------------------
+   reverse communication invoked by a Pair
+   nsize used only to set recv buffer limit
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::reverse_comm_pair(Pair *pair)
+{
+  CommTiled::reverse_comm_pair(pair);
+}
+
+/* ----------------------------------------------------------------------
+   forward communication invoked by a Fix
+   size/nsize used only to set recv buffer limit
+   size = 0 (default) -> use comm_forward from Fix
+   size > 0 -> Fix passes max size per atom
+   the latter is only useful if Fix does several comm modes,
+     some are smaller than max stored in its comm_forward
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::forward_comm_fix(Fix *fix, int size)
+{
+  CommTiled::forward_comm_fix(fix,size);
+}
+
+/* ----------------------------------------------------------------------
+   reverse communication invoked by a Fix
+   size/nsize used only to set recv buffer limit
+   size = 0 (default) -> use comm_forward from Fix
+   size > 0 -> Fix passes max size per atom
+   the latter is only useful if Fix does several comm modes,
+     some are smaller than max stored in its comm_forward
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::reverse_comm_fix(Fix *fix, int size)
+{
+  CommTiled::reverse_comm_fix(fix,size);
+}
+
+/* ----------------------------------------------------------------------
+   reverse communication invoked by a Fix with variable size data
+   query fix for all pack sizes to insure buf_send is big enough
+   handshake sizes before irregular comm to insure buf_recv is big enough
+   NOTE: how to setup one big buf recv with correct offsets ??
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::reverse_comm_fix_variable(Fix *fix)
+{
+  CommTiled::reverse_comm_fix_variable(fix);
+}
+
+/* ----------------------------------------------------------------------
+   forward communication invoked by a Compute
+   nsize used only to set recv buffer limit
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::forward_comm_compute(Compute *compute)
+{
+  CommTiled::forward_comm_compute(compute);
+}
+
+/* ----------------------------------------------------------------------
+   reverse communication invoked by a Compute
+   nsize used only to set recv buffer limit
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::reverse_comm_compute(Compute *compute)
+{
+  CommTiled::reverse_comm_compute(compute);
+}
+
+/* ----------------------------------------------------------------------
+   forward communication invoked by a Dump
+   nsize used only to set recv buffer limit
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::forward_comm_dump(Dump *dump)
+{
+  CommTiled::forward_comm_dump(dump);
+}
+
+/* ----------------------------------------------------------------------
+   reverse communication invoked by a Dump
+   nsize used only to set recv buffer limit
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::reverse_comm_dump(Dump *dump)
+{
+  CommTiled::reverse_comm_dump(dump);
+}
+
+/* ----------------------------------------------------------------------
+   forward communication of Nsize values in per-atom array
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::forward_comm_array(int nsize, double **array)
+{
+  CommTiled::forward_comm_array(nsize,array);
+}
+
+/* ----------------------------------------------------------------------
+   exchange info provided with all 6 stencil neighbors
+   NOTE: this method is currently not used
+------------------------------------------------------------------------- */
+
+int CommTiledKokkos::exchange_variable(int n, double *inbuf, double *&outbuf)
+{
+  CommTiled::exchange_variable(n,inbuf,outbuf);
+}
--- a/src/KOKKOS/comm_tiled_kokkos.h
+++ b/src/KOKKOS/comm_tiled_kokkos.h
@ -0,0 +1,59 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef LMP_COMM_TILED_KOKKOS_H
+#define LMP_COMM_TILED_KOKKOS_H
+
+#include "comm_tiled.h"
+#include "kokkos_type.h"
+
+namespace LAMMPS_NS {
+
+class CommTiledKokkos : public CommTiled {
+ public:
+  CommTiledKokkos(class LAMMPS *);
+  CommTiledKokkos(class LAMMPS *, class Comm *);
+  virtual ~CommTiledKokkos();
+
+  void forward_comm(int dummy = 0);    // forward comm of atom coords
+  void reverse_comm();                 // reverse comm of forces
+  void exchange();                     // move atoms to new procs
+  void borders();                      // setup list of atoms to comm
+
+  void forward_comm_pair(class Pair *);    // forward comm from a Pair
+  void reverse_comm_pair(class Pair *);    // reverse comm from a Pair
+  void forward_comm_fix(class Fix *, int size=0);
+                                                   // forward comm from a Fix
+  void reverse_comm_fix(class Fix *, int size=0);
+                                                   // reverse comm from a Fix
+  void reverse_comm_fix_variable(class Fix *);
+                                     // variable size reverse comm from a Fix
+  void forward_comm_compute(class Compute *);  // forward from a Compute
+  void reverse_comm_compute(class Compute *);  // reverse from a Compute
+  void forward_comm_dump(class Dump *);    // forward comm from a Dump
+  void reverse_comm_dump(class Dump *);    // reverse comm from a Dump
+
+  void forward_comm_array(int, double **);          // forward comm of array
+  int exchange_variable(int, double *, double *&);  // exchange on neigh stencil
+
+ private:
+
+};
+
+}
+
+#endif
+
+/* ERROR/WARNING messages:
+
+*/
--- a/src/KOKKOS/neighbor_kokkos.cpp
+++ b/src/KOKKOS/neighbor_kokkos.cpp
@ -280,15 +280,19 @@ void NeighborKokkos::choose_build(int index, NeighRequest *rq)
  if (rq->kokkos_host != 0) {
    PairPtrHost pb = NULL;
    if (rq->ghost) {
-      if (rq->full) pb = &NeighborKokkos::full_bin_kokkos<LMPHostType,0,1>;
-      else if (rq->half) &NeighborKokkos::full_bin_kokkos<LMPHostType,1,1>;
-      pair_build_host[index] = pb;
+      if (rq->full) {
+        if (rq->full_cluster) pb = &NeighborKokkos::full_bin_cluster_kokkos<LMPHostType>;
+        else pb = &NeighborKokkos::full_bin_kokkos<LMPHostType,0,1>;
+      }
+      else if (rq->half) pb = &NeighborKokkos::full_bin_kokkos<LMPHostType,1,1>;
    } else {
-      if (rq->full) pb = &NeighborKokkos::full_bin_kokkos<LMPHostType,0,0>;
+      if (rq->full) {
+        if (rq->full_cluster) pb = &NeighborKokkos::full_bin_cluster_kokkos<LMPHostType>;
+        else pb = &NeighborKokkos::full_bin_kokkos<LMPHostType,0,0>;
+      }
      else if (rq->half) pb = &NeighborKokkos::full_bin_kokkos<LMPHostType,1,0>;
-      pair_build_host[index] = pb;
    }
-    return;
+    pair_build_host[index] = pb;
  }
  if (rq->kokkos_device != 0) {
    PairPtrDevice pb = NULL;
--- 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
@ -13,11 +13,11 @@

 /* ----------------------------------------------------------------------
   Implementation of the CHARMM CMAP; adds an extra energy term for the
-   peptide backbone dihedrals.  The tools/ch2lmp.pl conversion script, which
-   generates an extra section in the LAMMPS data file, is needed in order to
-   generate the info used by this fix style.
+   peptide backbone dihedrals.  The tools/ch2lmp/charmm2lammps.pl
+   conversion script, which generates an extra section in the LAMMPS data
+   file, is needed in order to generate the info used by this fix style.

-   Contributing authors: 
+   Contributing authors:
   Xiaohu Hu, CMB/ORNL (hux2@ornl.gov)
   David Hyde-Volpe, Tigran Abramyan, and Robert A. Latour (Clemson University)
   Chris Lorenz (Kings College-London)
@ -27,11 +27,11 @@
   - MacKerell et al., J. Comput. Chem. 25(2004):1400-1415.
 -------------------------------------------------------------------------*/

-#include "mpi.h"
-#include "math.h"
-#include "stdlib.h"
-#include "string.h"
-#include "stdio.h"
+#include <mpi.h>
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
 #include "fix_cmap.h"
 #include "atom.h"
 #include "atom_vec.h"
@ -80,7 +80,7 @@ FixCMAP::FixCMAP(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)

  MPI_Comm_rank(world,&me);
  MPI_Comm_size(world,&nprocs);
- 
+
  // allocate memory for CMAP data

  memory->create(g_axis,CMAPDIM,"cmap:g_axis");
@ -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
@ -123,8 +122,8 @@ FixCMAP::~FixCMAP()
 {
  // unregister callbacks to this fix from Atom class

-  atom->delete_callback(id,0); 
-  atom->delete_callback(id,1); 
+  atom->delete_callback(id,0);
+  atom->delete_callback(id,1);

  memory->destroy(g_axis);
  memory->destroy(cmapgrid);
@ -173,7 +172,7 @@ void FixCMAP::init()
  }

  // pre-compute the derivatives of the maps
-  
+
  for (i = 0; i < 6; i++)
    set_map_derivatives(cmapgrid[i],d1cmapgrid[i],d2cmapgrid[i],d12cmapgrid[i]);

@ -194,7 +193,7 @@ void FixCMAP::setup(int vflag)
    ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
    post_force_respa(vflag,nlevels_respa-1,0);
    ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
-  }  
+  }
 }

 /* --------------------------------------------------------------------- */
@ -241,8 +240,8 @@ void FixCMAP::pre_neighbor()
      atom3 = atom->map(crossterm_atom3[i][m]);
      atom4 = atom->map(crossterm_atom4[i][m]);
      atom5 = atom->map(crossterm_atom5[i][m]);
-      
-      if (atom1 == -1 || atom2 == -1 || atom3 == -1 || 
+
+      if (atom1 == -1 || atom2 == -1 || atom3 == -1 ||
          atom4 == -1 || atom5 == -1) {
        char str[128];
        sprintf(str,"CMAP atoms "
@ -260,7 +259,7 @@ void FixCMAP::pre_neighbor()
      atom4 = domain->closest_image(i,atom4);
      atom5 = domain->closest_image(i,atom5);

-      if (i <= atom1 && i <= atom2 && i <= atom3 && 
+      if (i <= atom1 && i <= atom2 && i <= atom3 &&
          i <= atom4 && i <= atom5) {
        if (ncrosstermlist == maxcrossterm) {
          maxcrossterm += LISTDELTA;
@ -274,7 +273,7 @@ void FixCMAP::pre_neighbor()
        crosstermlist[ncrosstermlist][5] = crossterm_type[i][m];
        ncrosstermlist++;
      }
-    } 
+    }
  }
 }

@ -304,7 +303,7 @@ void FixCMAP::post_force(int vflag)
  double dpr32r43,dpr45r43,r43,vb12x,vb12y,vb12z,vb54x,vb54y,vb54z;
  // cross-term dihedral angles
  double phi,psi,phi1,psi1;
-  double f1[3],f2[3],f3[3],f4[3],f5[3],vcmap[6]; 
+  double f1[3],f2[3],f3[3],f4[3],f5[3],vcmap[6];
  double gs[4],d1gs[4],d2gs[4],d12gs[4];
  double engfraction;
  // vectors needed for the gradient/force calculation
@ -316,12 +315,12 @@ void FixCMAP::post_force(int vflag)
  // Definition of cross-term dihedrals

  //         phi dihedral
-  //   |--------------------| 
+  //   |--------------------|
  //   a1-----a2-----a3-----a4-----a5    cross-term atoms
  //   C      N      CA     C      N     cross-term atom types
  //          |--------------------|
  //               psi dihedral
- 
+
  double **x = atom->x;
  double **f = atom->f;
  int nlocal = atom->nlocal;
@ -342,7 +341,7 @@ void FixCMAP::post_force(int vflag)
    if (type == 0) continue;

    // calculate bond vectors for both dihedrals
-      
+
    // phi
    // vb21 = r2 - r1

@ -369,50 +368,50 @@ void FixCMAP::post_force(int vflag)
      vb43x = -1.0*vb34x;
      vb43y = -1.0*vb34y;
      vb43z = -1.0*vb34z;
-    
-      vb45x = x[i4][0] - x[i5][0]; 
-      vb45y = x[i4][1] - x[i5][1]; 
-      vb45z = x[i4][2] - x[i5][2]; 
+
+      vb45x = x[i4][0] - x[i5][0];
+      vb45y = x[i4][1] - x[i5][1];
+      vb45z = x[i4][2] - x[i5][2];
      vb54x = -1.0*vb45x;
      vb54y = -1.0*vb45y;
      vb54z = -1.0*vb45z;
- 
+
      // calculate normal vectors for planes that define the dihedral angles
-      
+
      a1x = vb12y*vb23z - vb12z*vb23y;
      a1y = vb12z*vb23x - vb12x*vb23z;
      a1z = vb12x*vb23y - vb12y*vb23x;
-      
+
      b1x = vb43y*vb23z - vb43z*vb23y;
      b1y = vb43z*vb23x - vb43x*vb23z;
      b1z = vb43x*vb23y - vb43y*vb23x;
-      
+
      a2x = vb23y*vb34z - vb23z*vb34y;
      a2y = vb23z*vb34x - vb23x*vb34z;
      a2z = vb23x*vb34y - vb23y*vb34x;
-      
+
      b2x = vb45y*vb43z - vb45z*vb43y;
      b2y = vb45z*vb43x - vb45x*vb43z;
      b2z = vb45x*vb43y - vb45y*vb43x;
-      
+
      // calculate terms used later in calculations
-      
+
      r32 = sqrt(vb32x*vb32x + vb32y*vb32y + vb32z*vb32z);
      a1sq = a1x*a1x + a1y*a1y + a1z*a1z;
      b1sq = b1x*b1x + b1y*b1y + b1z*b1z;
-      
+
      r43 = sqrt(vb43x*vb43x + vb43y*vb43y + vb43z*vb43z);
      a2sq = a2x*a2x + a2y*a2y + a2z*a2z;
      b2sq = b2x*b2x + b2y*b2y + b2z*b2z;
-      //if (a1sq<0.0001 || b1sq<0.0001 || a2sq<0.0001 || b2sq<0.0001) 
+      //if (a1sq<0.0001 || b1sq<0.0001 || a2sq<0.0001 || b2sq<0.0001)
      //  printf("a1sq b1sq a2sq b2sq: %f %f %f %f \n",a1sq,b1sq,a2sq,b2sq);
      if (a1sq<0.0001 || b1sq<0.0001 || a2sq<0.0001 || b2sq<0.0001) continue;
      dpr21r32 = vb21x*vb32x + vb21y*vb32y + vb21z*vb32z;
      dpr34r32 = vb34x*vb32x + vb34y*vb32y + vb34z*vb32z;
      dpr32r43 = vb32x*vb43x + vb32y*vb43y + vb32z*vb43z;
      dpr45r43 = vb45x*vb43x + vb45y*vb43y + vb45z*vb43z;
-      
-      // calculate the backbone dihedral angles as VMD and GROMACS 
+
+      // calculate the backbone dihedral angles as VMD and GROMACS

      phi = dihedral_angle_atan2(vb21x,vb21y,vb21z,a1x,a1y,a1z,b1x,b1y,b1z,r32);
      psi = dihedral_angle_atan2(vb32x,vb32y,vb32z,a2x,a2y,a2z,b2x,b2y,b2z,r43);
@ -425,7 +424,7 @@ void FixCMAP::post_force(int vflag)
      psi1 = psi;
      if (psi1 < 0.0) psi1 += 360.0;

-      // find the neighbor grid point index 
+      // find the neighbor grid point index

      li1 = int(((phi1+CMAPXMIN2)/CMAPDX)+((CMAPDIM*1.0)/2.0));
      li2 = int(((psi1+CMAPXMIN2)/CMAPDX)+((CMAPDIM*1.0)/2.0));
@ -436,87 +435,87 @@ void FixCMAP::post_force(int vflag)
      mli4 = li4 % CMAPDIM;
      mli31 = (li3+1) % CMAPDIM;
      mli41 = (li4+1)  %CMAPDIM;
-      mli1 = li1 % CMAPDIM; 
+      mli1 = li1 % CMAPDIM;
      mli2 = li2 % CMAPDIM;
      mli11 = (li1+1) % CMAPDIM;
      mli21 = (li2+1)  %CMAPDIM;
      t1 = type-1;
      if (t1 < 0 || t1 > 5) error->all(FLERR,"Invalid CMAP crossterm_type");
-      
-      // determine the values and derivatives for the grid square points 
+
+      // determine the values and derivatives for the grid square points

      gs[0] = cmapgrid[t1][mli3][mli4];
      gs[1] = cmapgrid[t1][mli31][mli4];
      gs[2] = cmapgrid[t1][mli31][mli41];
      gs[3] = cmapgrid[t1][mli3][mli41];
-      d1gs[0] = d1cmapgrid[t1][mli1][mli2];  
+      d1gs[0] = d1cmapgrid[t1][mli1][mli2];
      d1gs[1] = d1cmapgrid[t1][mli11][mli2];
      d1gs[2] = d1cmapgrid[t1][mli11][mli21];
      d1gs[3] = d1cmapgrid[t1][mli1][mli21];
-      
-      d2gs[0] = d2cmapgrid[t1][mli1][mli2];  
+
+      d2gs[0] = d2cmapgrid[t1][mli1][mli2];
      d2gs[1] = d2cmapgrid[t1][mli11][mli2];
      d2gs[2] = d2cmapgrid[t1][mli11][mli21];
      d2gs[3] = d2cmapgrid[t1][mli1][mli21];
-      
-      d12gs[0] = d12cmapgrid[t1][mli1][mli2];  
+
+      d12gs[0] = d12cmapgrid[t1][mli1][mli2];
      d12gs[1] = d12cmapgrid[t1][mli11][mli2];
      d12gs[2] = d12cmapgrid[t1][mli11][mli21];
      d12gs[3] = d12cmapgrid[t1][mli1][mli21];
-      
-      // calculate the cmap energy and the gradient (dE/dphi,dE/dpsi) 
-      
-      bc_interpol(phi,psi,li3,li4,gs,d1gs,d2gs,d12gs);
-      
-      // sum up cmap energy contributions
-      
-      engfraction = 0.2 * E;
-      if (i1 < nlocal) ecmap += engfraction; 
-      if (i2 < nlocal) ecmap += engfraction; 
-      if (i3 < nlocal) ecmap += engfraction; 
-      if (i4 < nlocal) ecmap += engfraction; 
-      if (i5 < nlocal) ecmap += engfraction; 

-      // calculate the derivatives dphi/dr_i 
-      
+      // calculate the cmap energy and the gradient (dE/dphi,dE/dpsi)
+
+      bc_interpol(phi,psi,li3,li4,gs,d1gs,d2gs,d12gs);
+
+      // sum up cmap energy contributions
+
+      engfraction = 0.2 * E;
+      if (i1 < nlocal) ecmap += engfraction;
+      if (i2 < nlocal) ecmap += engfraction;
+      if (i3 < nlocal) ecmap += engfraction;
+      if (i4 < nlocal) ecmap += engfraction;
+      if (i5 < nlocal) ecmap += engfraction;
+
+      // calculate the derivatives dphi/dr_i
+
      dphidr1x = 1.0*r32/a1sq*a1x;
      dphidr1y = 1.0*r32/a1sq*a1y;
      dphidr1z = 1.0*r32/a1sq*a1z;
-      
-      dphidr2x = -1.0*r32/a1sq*a1x - dpr21r32/a1sq/r32*a1x + 
+
+      dphidr2x = -1.0*r32/a1sq*a1x - dpr21r32/a1sq/r32*a1x +
        dpr34r32/b1sq/r32*b1x;
-      dphidr2y = -1.0*r32/a1sq*a1y - dpr21r32/a1sq/r32*a1y + 
+      dphidr2y = -1.0*r32/a1sq*a1y - dpr21r32/a1sq/r32*a1y +
        dpr34r32/b1sq/r32*b1y;
-      dphidr2z = -1.0*r32/a1sq*a1z - dpr21r32/a1sq/r32*a1z + 
+      dphidr2z = -1.0*r32/a1sq*a1z - dpr21r32/a1sq/r32*a1z +
        dpr34r32/b1sq/r32*b1z;

      dphidr3x = dpr34r32/b1sq/r32*b1x - dpr21r32/a1sq/r32*a1x - r32/b1sq*b1x;
      dphidr3y = dpr34r32/b1sq/r32*b1y - dpr21r32/a1sq/r32*a1y - r32/b1sq*b1y;
      dphidr3z = dpr34r32/b1sq/r32*b1z - dpr21r32/a1sq/r32*a1z - r32/b1sq*b1z;
-    
+
      dphidr4x = r32/b1sq*b1x;
      dphidr4y = r32/b1sq*b1y;
      dphidr4z = r32/b1sq*b1z;

-      // calculate the derivatives dpsi/dr_i 
-      
+      // calculate the derivatives dpsi/dr_i
+
      dpsidr1x = 1.0*r43/a2sq*a2x;
      dpsidr1y = 1.0*r43/a2sq*a2y;
      dpsidr1z = 1.0*r43/a2sq*a2z;
-      
+
      dpsidr2x = r43/a2sq*a2x + dpr32r43/a2sq/r43*a2x - dpr45r43/b2sq/r43*b2x;
      dpsidr2y = r43/a2sq*a2y + dpr32r43/a2sq/r43*a2y - dpr45r43/b2sq/r43*b2y;
      dpsidr2z = r43/a2sq*a2z + dpr32r43/a2sq/r43*a2z - dpr45r43/b2sq/r43*b2z;
-      
+
      dpsidr3x = dpr45r43/b2sq/r43*b2x - dpr32r43/a2sq/r43*a2x - r43/b2sq*b2x;
      dpsidr3y = dpr45r43/b2sq/r43*b2y - dpr32r43/a2sq/r43*a2y - r43/b2sq*b2y;
      dpsidr3z = dpr45r43/b2sq/r43*b2z - dpr32r43/a2sq/r43*a2z - r43/b2sq*b2z;
-      
+
      dpsidr4x = r43/b2sq*b2x;
      dpsidr4y = r43/b2sq*b2y;
      dpsidr4z = r43/b2sq*b2z;

-      // calculate forces on cross-term atoms: F = -(dE/dPhi)*(dPhi/dr) 
+      // calculate forces on cross-term atoms: F = -(dE/dPhi)*(dPhi/dr)

      f1[0] = dEdPhi*dphidr1x;
      f1[1] = dEdPhi*dphidr1y;
@ -533,9 +532,9 @@ void FixCMAP::post_force(int vflag)
      f5[0] = -dEdPsi*dpsidr4x;
      f5[1] = -dEdPsi*dpsidr4y;
      f5[2] = -dEdPsi*dpsidr4z;
-      
+
      // apply force to each of the 5 atoms
-      
+
      if (i1 < nlocal) {
        f[i1][0] += f1[0];
        f[i1][1] += f1[1];
@ -561,7 +560,7 @@ void FixCMAP::post_force(int vflag)
        f[i5][1] += f5[1];
        f[i5][2] += f5[2];
      }
-      
+
      // tally energy and/or virial

      if (evflag) {
@ -622,15 +621,18 @@ double FixCMAP::compute_scalar()

 void FixCMAP::read_grid_map(char *cmapfile)
 {
-  char line[MAXLINE];
-  char *chunk;
+  char linebuf[MAXLINE];
+  char *chunk,*line;
  int i1, i2, i3, i4, i5, i6, j1, j2, j3, j4, j5, j6, counter;
-  
-  FILE *fp = fopen(cmapfile,"r");
-  if (fp == NULL) {
-    char str[128];
-    sprintf(str,"Cannot open fix cmap file %s",cmapfile);
-    error->one(FLERR,str);
+
+  FILE *fp = NULL;
+  if (comm->me == 0) {
+    fp = force->open_potential(cmapfile);
+    if (fp == NULL) {
+      char str[128];
+      sprintf(str,"Cannot open fix cmap file %s",cmapfile);
+      error->one(FLERR,str);
+    }
  }

  for (int ix1 = 0; ix1 < 6; ix1++)
@ -642,8 +644,24 @@ void FixCMAP::read_grid_map(char *cmapfile)
  i1 = i2 = i3 = i4 = i5 = i6 = 0;
  j1 = j2 = j3 = j4 = j5 = j6 = 0;

-  while (fgets(line,MAXLINE,fp) != NULL) {
-    if (line == "" || line[0] == '#') { ;; }
+  int done = 0;
+
+  while (!done) {
+    // only read on rank 0 and broadcast to all other ranks
+    if (comm->me == 0)
+      done = (fgets(linebuf,MAXLINE,fp) == NULL);
+
+    MPI_Bcast(&done,1,MPI_INT,0,world);
+    if (done) continue;
+
+    MPI_Bcast(linebuf,MAXLINE,MPI_CHAR,0,world);
+
+    // remove leading whitespace
+    line = linebuf;
+    while (line && (*line == ' ' || *line == '\t' || *line == '\r')) ++line;
+
+    // skip if empty line or comment
+    if (!line || *line =='\n' || *line == '\0' || *line == '#') continue;

    // read in the cmap grid point values
    // NOTE: The order to read the 6 grid maps is HARD-CODED, thus errors
@ -655,101 +673,99 @@ void FixCMAP::read_grid_map(char *cmapfile)
    // 3. Proline map
    // 4. Two adjacent prolines map
    // 5. Glycine map
-    // 6. Glycine before proline map	
-    
-    else {
-      chunk = strtok(line, " \r\n");
-      while (chunk != NULL) {
+    // 6. Glycine before proline map

-        // alanine map 
-        
-        if (counter < CMAPDIM*CMAPDIM) {
-          cmapgrid[0][i1][j1] = atof(chunk);
-          chunk = strtok(NULL, " \r\n");
-          j1++;
-          if (j1 == CMAPDIM) {
-            j1 = 0;
-            i1++;
-          }
-          counter++;
+    chunk = strtok(line, " \r\n");
+    while (chunk != NULL) {
+
+      // alanine map
+
+      if (counter < CMAPDIM*CMAPDIM) {
+        cmapgrid[0][i1][j1] = atof(chunk);
+        chunk = strtok(NULL, " \r\n");
+        j1++;
+        if (j1 == CMAPDIM) {
+          j1 = 0;
+          i1++;
        }
-        
-        // alanine-proline map 
-        
-        else if (counter >= CMAPDIM*CMAPDIM &&
-                 counter < 2*CMAPDIM*CMAPDIM) {
-          cmapgrid[1][i2][j2]= atof(chunk);
-          chunk = strtok(NULL, " \r\n");
-          j2++;
-          if (j2 == CMAPDIM) {
-            j2 = 0;
-            i2++;
-          }
-          counter++;
-        }
-        
-        // proline map 
-        
-        else if (counter >= 2*CMAPDIM*CMAPDIM &&
-                 counter < 3*CMAPDIM*CMAPDIM) {
-          cmapgrid[2][i3][j3] = atof(chunk);
-          chunk = strtok(NULL, " \r\n");
-          j3++;
-          if (j3 == CMAPDIM) {
-            j3 = 0;
-            i3++;
-          }
-          counter++;
-        }
-        
-        // 2 adjacent prolines map 
-        
-        else if (counter >= 3*CMAPDIM*CMAPDIM &&
-                 counter < 4*CMAPDIM*CMAPDIM) {
-          cmapgrid[3][i4][j4] = atof(chunk);
-          chunk = strtok(NULL, " \r\n");
-          j4++;
-          if (j4 == CMAPDIM) {
-            j4 = 0;
-            i4++;
-          }
-          counter++;
-        }
-        
-        // glycine map 
-        
-        else if (counter >= 4*CMAPDIM*CMAPDIM &&
-                 counter < 5*CMAPDIM*CMAPDIM) {
-          cmapgrid[4][i5][j5] = atof(chunk);
-          chunk = strtok(NULL, " \r\n");
-          j5++;
-          if (j5 == CMAPDIM) {
-            j5 = 0;
-            i5++;
-          }
-          counter++;
-        }
-        
-        // glycine-proline map 
-        
-        else if (counter >= 5*CMAPDIM*CMAPDIM &&
-                 counter < 6*CMAPDIM*CMAPDIM) {
-          cmapgrid[5][i6][j6] = atof(chunk);
-          chunk = strtok(NULL, " \r\n");
-          j6++;
-          if (j6 == CMAPDIM) {
-            j6 = 0;
-            i6++;
-          }
-          counter++;
-        }
-        
-        else break;
+        counter++;
      }
+
+      // alanine-proline map
+
+      else if (counter >= CMAPDIM*CMAPDIM &&
+               counter < 2*CMAPDIM*CMAPDIM) {
+        cmapgrid[1][i2][j2]= atof(chunk);
+        chunk = strtok(NULL, " \r\n");
+        j2++;
+        if (j2 == CMAPDIM) {
+          j2 = 0;
+          i2++;
+        }
+        counter++;
+      }
+
+      // proline map
+
+      else if (counter >= 2*CMAPDIM*CMAPDIM &&
+               counter < 3*CMAPDIM*CMAPDIM) {
+        cmapgrid[2][i3][j3] = atof(chunk);
+        chunk = strtok(NULL, " \r\n");
+        j3++;
+        if (j3 == CMAPDIM) {
+          j3 = 0;
+          i3++;
+        }
+        counter++;
+      }
+
+      // 2 adjacent prolines map
+
+      else if (counter >= 3*CMAPDIM*CMAPDIM &&
+               counter < 4*CMAPDIM*CMAPDIM) {
+        cmapgrid[3][i4][j4] = atof(chunk);
+        chunk = strtok(NULL, " \r\n");
+        j4++;
+        if (j4 == CMAPDIM) {
+          j4 = 0;
+          i4++;
+        }
+        counter++;
+      }
+
+      // glycine map
+
+      else if (counter >= 4*CMAPDIM*CMAPDIM &&
+               counter < 5*CMAPDIM*CMAPDIM) {
+        cmapgrid[4][i5][j5] = atof(chunk);
+        chunk = strtok(NULL, " \r\n");
+        j5++;
+        if (j5 == CMAPDIM) {
+          j5 = 0;
+          i5++;
+        }
+        counter++;
+      }
+
+      // glycine-proline map
+
+      else if (counter >= 5*CMAPDIM*CMAPDIM &&
+               counter < 6*CMAPDIM*CMAPDIM) {
+        cmapgrid[5][i6][j6] = atof(chunk);
+        chunk = strtok(NULL, " \r\n");
+        j6++;
+        if (j6 == CMAPDIM) {
+          j6 = 0;
+          i6++;
+        }
+        counter++;
+      }
+
+      else break;
    }
  }
-  
-  fclose(fp);
+
+  if (comm->me == 0) fclose(fp);
 }

 /* ---------------------------------------------------------------------- */
@ -783,18 +799,18 @@ void FixCMAP::spline(double *y, double *ddy, int n)
 /* ---------------------------------------------------------------------- */

 void FixCMAP::spl_interpolate(double x, double *y, double *ddy, double &yo,
-		double &dyo)
+                              double &dyo)
 {
  // perform a 1D cubic spline interpolation

  int ix;
  double a,b,a1,b1,a2,b2;
-  
+
  ix = int((x-CMAPXMIN)/CMAPDX-(1./2.));

  a = (CMAPXMIN+(ix*1.0)*CMAPDX-x)/CMAPDX;
  b = (x-CMAPXMIN-(((ix-1)*1.0)*CMAPDX))/CMAPDX;
-  
+
  a1 = a*a*a-a;
  b1 = b*b*b-b;

@ -807,7 +823,7 @@ void FixCMAP::spl_interpolate(double x, double *y, double *ddy, double &yo,
 /* ---------------------------------------------------------------------- */

 void FixCMAP::set_map_derivatives(double **map, double **d1yo, double **d2yo,
-                                   double **d12yo)
+                                  double **d12yo)
 {
  // precompute the gradient and cross-derivatives of the map grid points.
  // use the bicubic spline to calculate the derivatives
@ -833,7 +849,7 @@ void FixCMAP::set_map_derivatives(double **map, double **d1yo, double **d2yo,
  memory->create(tddmap,CMAPDIM*2,CMAPDIM*2,"cmap:tddmap");

  // periodically expand the original map
-  // use the expanded map for bicubic spline interpolation, 
+  // use the expanded map for bicubic spline interpolation,
  //   which is used to obtain the derivatives
  // actual interpolation is done with bicubic interpolation

@ -844,7 +860,7 @@ void FixCMAP::set_map_derivatives(double **map, double **d1yo, double **d2yo,
      tmap[i][j] = map[ii][jj];
    }
  }
-  
+
  for (i = 0; i < CMAPDIM*2; i++)
    spline(tmap[i], tddmap[i], CMAPDIM*2);

@ -869,7 +885,7 @@ void FixCMAP::set_map_derivatives(double **map, double **d1yo, double **d2yo,
        tmp_y[k] = tyyk;
        tmp_dy[k] = tdyk;
      }
-      
+
      spline(tmp_y,tmp_ddy,CMAPDIM+xm+xm);
      ix = int((phi-CMAPXMIN)/CMAPDX);
      a = (CMAPXMIN+((ix+1)*1.0)*CMAPDX-phi)/CMAPDX;
@ -899,7 +915,7 @@ void FixCMAP::set_map_derivatives(double **map, double **d1yo, double **d2yo,
      d12yo[i%p][j%p] = d12y;
    }
  }
-  
+
  memory->destroy(tmp_y);
  memory->destroy(tmp_dy);
  memory->destroy(tmp_ddy);
@ -911,7 +927,7 @@ void FixCMAP::set_map_derivatives(double **map, double **d1yo, double **d2yo,

 double FixCMAP::dihedral_angle_atan2(double fx, double fy, double fz,
                                      double ax, double ay, double az,
-                                      double bx, double by, double bz, 
+                                      double bx, double by, double bz,
                                      double absg)
 {
  // calculate the dihedral angle
@ -953,19 +969,19 @@ void FixCMAP::bc_coeff(double *gs, double *d1gs, double *d2gs, double *d12gs)
      2,-2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0,
      -6, 6,-6, 6,-3,-3, 3, 3,-4, 4, 2,-2,-2,-2,-1,-1,
-      4,-4, 4,-4, 2, 2,-2,-2, 2,-2,-2, 2, 1, 1, 1, 1 
+      4,-4, 4,-4, 2, 2,-2,-2, 2,-2,-2, 2, 1, 1, 1, 1
    };

-  int i, j, k, l, in;
+  int i, j, k, in;
  double xx, x[16];
-  
+
  for (i = 0; i < 4; i++) {
    x[i] = gs[i];
    x[i+4] = d1gs[i]*CMAPDX;
    x[i+8] = d2gs[i]*CMAPDX;
    x[i+12] = d12gs[i]*CMAPDX*CMAPDX;
  }
-  
+
  in = 0;
  for (i = 0; i < 4; i++) {
    for (j = 0; j < 4; j++) {
@ -1006,7 +1022,7 @@ void FixCMAP::bc_interpol(double x1, double x2, int low1, int low2, double *gs,
    dEdPhi = u*dEdPhi + (3.0*cij[3][i]*t+2.0*cij[2][i])*t+cij[1][i];
    dEdPsi = t*dEdPsi + (3.0*cij[i][3]*u+2.0*cij[i][2])*u+cij[i][1];
  }
-  
+
  dEdPhi *= (180.0/MY_PI/CMAPDX);
  dEdPsi *= (180.0/MY_PI/CMAPDX);
 }
@ -1023,7 +1039,7 @@ void FixCMAP::read_data_header(char *line)
    sscanf(line,BIGINT_FORMAT,&ncmap);
  } else error->all(FLERR,"Invalid read data header line for fix cmap");

-  // didn't set in constructor b/c this fix could be defined 
+  // didn't set in constructor b/c this fix could be defined
  // before newton command

  newton_bond = force->newton_bond;
@ -1117,7 +1133,7 @@ void FixCMAP::read_data_section(char *keyword, int n, char *buf,
      crossterm_atom5[m][num_crossterm[m]] = atom5;
      num_crossterm[m]++;
    }
-    
+
    if ((m = atom->map(atom5)) >= 0) {
      if (num_crossterm[m] == CMAPMAX)
        error->one(FLERR,"Too many CMAP crossterms for one atom");
@ -1226,7 +1242,7 @@ void FixCMAP::write_data_section(int mth, FILE *fp,
                                  int n, double **buf, int index)
 {
  for (int i = 0; i < n; i++)
-    fprintf(fp,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT 
+    fprintf(fp,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT
            " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT "\n",
            index+i,(int) ubuf(buf[i][0]).i,(tagint) ubuf(buf[i][1]).i,
            (tagint) ubuf(buf[i][2]).i,(tagint) ubuf(buf[i][3]).i,
@ -1261,8 +1277,8 @@ void FixCMAP::restart(char *buf)
  ncmap = *((bigint *) buf);
 }

-/* ---------------------------------------------------------------------- 
-   pack values in local atom-based arrays for restart file 
+/* ----------------------------------------------------------------------
+   pack values in local atom-based arrays for restart file
 ------------------------------------------------------------------------- */

 int FixCMAP::pack_restart(int i, double *buf)
@ -1275,14 +1291,14 @@ int FixCMAP::pack_restart(int i, double *buf)
    buf[n++] = ubuf(crossterm_atom3[i][m]).d;
    buf[n++] = ubuf(crossterm_atom4[i][m]).d;
    buf[n++] = ubuf(crossterm_atom5[i][m]).d;
-  }  
+  }
  buf[0] = n;

-  return n; 
+  return n;
 }

-/* ---------------------------------------------------------------------- 
-   unpack values from atom->extra array to restart the fix 
+/* ----------------------------------------------------------------------
+   unpack values from atom->extra array to restart the fix
 ------------------------------------------------------------------------- */

 void FixCMAP::unpack_restart(int nlocal, int nth)
@ -1290,9 +1306,9 @@ void FixCMAP::unpack_restart(int nlocal, int nth)
  double **extra = atom->extra;

  // skip to Nth set of extra values
-   
+
   int n = 0;
-   for (int i = 0; i < nth; i++) n += static_cast<int> (extra[nlocal][n]); 
+   for (int i = 0; i < nth; i++) n += static_cast<int> (extra[nlocal][n]);

   int count = static_cast<int> (extra[nlocal][n++]);
   num_crossterm[nlocal] = (count-1)/6;
@ -1307,17 +1323,17 @@ void FixCMAP::unpack_restart(int nlocal, int nth)
   }
 }

-/* ---------------------------------------------------------------------- 
-   maxsize of any atom's restart data 
+/* ----------------------------------------------------------------------
+   maxsize of any atom's restart data
 ------------------------------------------------------------------------- */

 int FixCMAP::maxsize_restart()
 {
-  return 1 + CMAPMAX*6; 
+  return 1 + CMAPMAX*6;
 }

-/* ---------------------------------------------------------------------- 
-   size of atom nlocal's restart data 
+/* ----------------------------------------------------------------------
+   size of atom nlocal's restart data
 ------------------------------------------------------------------------- */

 int FixCMAP::size_restart(int nlocal)
@ -1330,7 +1346,7 @@ int FixCMAP::size_restart(int nlocal)
 ------------------------------------------------------------------------- */

 void FixCMAP::grow_arrays(int nmax)
-{ 
+{
  num_crossterm = memory->grow(num_crossterm,nmax,"cmap:num_crossterm");
  crossterm_type = memory->grow(crossterm_type,nmax,CMAPMAX,
                                "cmap:crossterm_type");
@ -1382,7 +1398,7 @@ void FixCMAP::set_arrays(int i)
 /* ----------------------------------------------------------------------
   pack values in local atom-based array for exchange with another proc
 ------------------------------------------------------------------------- */
-  
+
 int FixCMAP::pack_exchange(int i, double *buf)
 {
  int n = 0;
--- a/src/MOLECULE/fix_cmap.h
+++ b/src/MOLECULE/fix_cmap.h
@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, Sandia National Laboratories
   Steve Plimpton, sjplimp@sandia.gov
--- a/src/Make.py
+++ b/src/Make.py
@ -26,7 +26,8 @@ libclasses = ("atc","awpmd","colvars","cuda","gpu","h5md",
 buildclasses = ("intel","kokkos")
 makeclasses = ("cc","flags","mpi","fft","jpg","png")

-setargs = ("gzip","#gzip","ffmpeg","#ffmpeg","smallbig","bigbig","smallsmall","exceptions","#exceptions")
+setargs = ("gzip","#gzip","ffmpeg","#ffmpeg","smallbig","bigbig",
+           "smallsmall","exceptions","#exceptions")
 actionargs = ("lib-all","file","clean","exe")

 gpubuildflag = 0
@ -85,7 +86,9 @@ def switch2str(switches,switch_order):
 def compile_check(compiler,ccflags,warn):
  open("tmpauto.cpp",'w').write("int main(int, char **) {}\n")
  tmp = "%s %s -c tmpauto.cpp" % (compiler,ccflags)
-  txt = subprocess.check_output(tmp,stderr=subprocess.STDOUT,shell=True).decode()
+  try: txt = subprocess.check_output(tmp,stderr=subprocess.STDOUT,
+                                     shell=True).decode()
+  except subprocess.CalledProcessError as e: txt = e.output
  flag = 1
  if txt or not os.path.isfile("tmpauto.o"):
    flag = 0
@ -104,7 +107,9 @@ def compile_check(compiler,ccflags,warn):
 def link_check(linker,linkflags,libs,warn):
  open("tmpauto.cpp",'w').write("int main(int, char **) {}\n")
  tmp = "%s %s -o tmpauto tmpauto.cpp %s" % (linker,linkflags,libs)
-  txt = subprocess.check_output(tmp,stderr=subprocess.STDOUT,shell=True).decode()
+  try: txt = subprocess.check_output(tmp,stderr=subprocess.STDOUT,
+                                     shell=True).decode()
+  except subprocess.CalledProcessError as e: txt = e.output
  flag = 1
  if txt or not os.path.isfile("tmpauto"):
    flag = 0
@ -530,7 +535,8 @@ class Actions(object):
    
    if caller == "file" or "file" not in self.alist:
      # make certain that 'MAKE/MINE' folder exists.
-      subprocess.check_output("mkdir -p %s/MAKE/MINE" % dir.src,stderr=subprocess.STDOUT,shell=True)
+      subprocess.check_output("mkdir -p %s/MAKE/MINE" % dir.src,
+                              stderr=subprocess.STDOUT,shell=True)
      make.write("%s/MAKE/MINE/Makefile.auto" % dir.src,1)
      print("Created src/MAKE/MINE/Makefile.auto")
      
@ -588,7 +594,7 @@ class Actions(object):
    else:
      print(tmp)
      try: subprocess.check_output(tmp,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)

    if not os.path.isfile("%s/lmp_auto" % dir.src):
      error('Unsuccessful "make auto"')
@ -857,9 +863,12 @@ class Packages(object):
    if self.plist: print("Installing packages ...")
    for one in self.plist:
      if one == "orig": continue
-      subprocess.check_output("cd %s; make %s" % (dir.src,one),stderr=subprocess.STDOUT,shell=True)
+      subprocess.check_output("cd %s; make %s" % (dir.src,one),
+                              stderr=subprocess.STDOUT,shell=True)
    if self.plist and verbose:
-      txt = subprocess.check_output("cd %s; make ps" % dir.src,stderr=subprocess.STDOUT,shell=True).decode()
+      txt = subprocess.check_output("cd %s; make ps" % dir.src,
+                                    stderr=subprocess.STDOUT,
+                                    shell=True).decode()
      print("Package status after installation:")
      print(txt)
      
@ -869,12 +878,16 @@ class Packages(object):
  def uninstall(self):
    if not self.plist or self.plist[-1] != "orig": return
    print("Restoring packages to original state ...")
-    subprocess.check_output("cd %s; make no-all" % dir.src,stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; make no-all" % dir.src,
+                            stderr=subprocess.STDOUT,shell=True)
    for one in self.all:
      if self.original[one]:
-        subprocess.check_output("cd %s; make yes-%s" % (dir.src,one),stderr=subprocess.STDOUT,shell=True)
+        subprocess.check_output("cd %s; make yes-%s" % (dir.src,one),
+                                stderr=subprocess.STDOUT,shell=True)
    if verbose:
-      txt = subprocess.check_output("cd %s; make ps" % dir.src,stderr=subprocess.STDOUT,shell=True).decode()
+      txt = subprocess.check_output("cd %s; make ps" % dir.src,
+                                    stderr=subprocess.STDOUT,
+                                    shell=True).decode()
      print("Restored package status:")
      print(txt)
      
@ -968,7 +981,8 @@ class Settings(object):
  def help(self):
    return """
 -s set1 set2 ...
-  possible settings = gzip #gzip ffmpeg #ffmpeg smallbig bigbig smallsmall exceptions #exceptions
+  possible settings = gzip #gzip ffmpeg #ffmpeg
+                      smallbig bigbig smallsmall exceptions #exceptions
  alter LAMMPS ifdef settings in Makefile.auto
    only happens if new Makefile.auto is created by use of "file" action
  gzip and #gzip turn on/off LAMMPS_GZIP setting
@ -1057,7 +1071,8 @@ class ATC(object):
      make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps)
    make.write("%s/Makefile.auto" % libdir)
    
-    subprocess.check_output("cd %s; make -f Makefile.auto clean" % libdir,stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; make -f Makefile.auto clean" %
+                            libdir,stderr=subprocess.STDOUT,shell=True)
    if jmake: txt = "cd %s; make -j %d -f Makefile.auto" % (libdir,jmake.n)
    else: txt = "cd %s; make -f Makefile.auto" % libdir
    
@ -1066,7 +1081,7 @@ class ATC(object):
    if verbose: subprocess.call(txt,shell=True)
    else:
      try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)

    if not os.path.isfile("%s/libatc.a" % libdir) or \
          not os.path.isfile("%s/Makefile.lammps" % libdir):
@ -1108,7 +1123,8 @@ class AWPMD(object):
      make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps)
    make.write("%s/Makefile.auto" % libdir)

-    subprocess.check_output("cd %s; make -f Makefile.auto clean" % libdir,stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; make -f Makefile.auto clean" %
+                            libdir,stderr=subprocess.STDOUT,shell=True)
    if jmake: txt = "cd %s; make -j %d -f Makefile.auto" % (libdir,jmake.n)
    else: txt = "cd %s; make -f Makefile.auto" % libdir

@ -1117,7 +1133,7 @@ class AWPMD(object):
    if verbose: subprocess.call(txt,shell=True)
    else:
      try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)
   
    if not os.path.isfile("%s/libawpmd.a" % libdir) or \
          not os.path.isfile("%s/Makefile.lammps" % libdir):
@ -1159,7 +1175,8 @@ class COLVARS(object):
      make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps)
    make.write("%s/Makefile.auto" % libdir)

-    subprocess.check_output("cd %s; make -f Makefile.auto clean" % libdir,stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; make -f Makefile.auto clean" %
+                            libdir,stderr=subprocess.STDOUT,shell=True)
    if jmake: txt = "cd %s; make -j %d -f Makefile.auto" % (libdir,jmake.n)
    else: txt = "cd %s; make -f Makefile.auto" % libdir

@ -1168,7 +1185,7 @@ class COLVARS(object):
    if verbose: subprocess.call(txt,shell=True)
    else:
      try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)

    if not os.path.isfile("%s/libcolvars.a" % libdir) or \
          not os.path.isfile("%s/Makefile.lammps" % libdir):
@ -1211,7 +1228,8 @@ class CUDA(object):
          
  def build(self): 
    libdir = dir.lib + "/cuda"
-    subprocess.check_output("cd %s; make clean" % libdir,stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; make clean" % libdir,
+                            stderr=subprocess.STDOUT,shell=True)
    if self.mode == "double": n = 2
    elif self.mode == "mixed": n = 3
    elif self.mode == "single": n = 1
@ -1225,7 +1243,7 @@ class CUDA(object):
    if verbose: subprocess.call(txt,shell=True)
    else:
      try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)

    if not os.path.isfile("%s/liblammpscuda.a" % libdir) or \
          not os.path.isfile("%s/Makefile.lammps" % libdir):
@ -1306,7 +1324,8 @@ class GPU(object):
    make = "make"
    if "shannon" == platform.node(): make = "srun make"

-    subprocess.check_output("cd %s; %s -f Makefile.auto clean" % (libdir,make),stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; %s -f Makefile.auto clean" %
+                            (libdir,make),stderr=subprocess.STDOUT,shell=True)
    if jmake: txt = "cd %s; %s -j %d -f Makefile.auto" % (libdir,make,jmake.n)
    else: txt = "cd %s; %s -f Makefile.auto" % (libdir,make)

@ -1315,7 +1334,7 @@ class GPU(object):
    if verbose: subprocess.call(txt,shell=True)
    else:
      try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)

    if not os.path.isfile("%s/libgpu.a" % libdir) or \
          not os.path.isfile("%s/Makefile.lammps" % libdir):
@ -1357,7 +1376,8 @@ class H5MD(object):
      make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps)
    make.write("%s/Makefile.auto" % libdir)

-    subprocess.check_output("cd %s; make clean" % libdir,stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; make clean" % libdir,
+                            stderr=subprocess.STDOUT,shell=True)
    txt = "cd %s; make" % libdir

    # if verbose, print output as build proceeds, else only print if fails
@ -1365,7 +1385,7 @@ class H5MD(object):
    if verbose: subprocess.call(txt,shell=True)
    else:
      try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)

    if not os.path.isfile("%s/libch5md.a" % libdir) or \
          not os.path.isfile("%s/Makefile.lammps" % libdir):
@ -1407,7 +1427,8 @@ class MEAM(object):
      make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps)
    make.write("%s/Makefile.auto" % libdir)

-    subprocess.check_output("cd %s; make -f Makefile.auto clean" % libdir,stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; make -f Makefile.auto clean" %
+                            libdir,stderr=subprocess.STDOUT,shell=True)
    # do not use -j for MEAM build, parallel build does not work
    txt = "cd %s; make -f Makefile.auto" % libdir

@ -1416,7 +1437,7 @@ class MEAM(object):
    if verbose: subprocess.call(txt,shell=True)
    else:
      try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)

    if not os.path.isfile("%s/libmeam.a" % libdir) or \
          not os.path.isfile("%s/Makefile.lammps" % libdir):
@ -1458,7 +1479,8 @@ class POEMS(object):
      make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps)
    make.write("%s/Makefile.auto" % libdir)

-    subprocess.check_output("cd %s; make -f Makefile.auto clean" % libdir,stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; make -f Makefile.auto clean" % libdir,
+                            stderr=subprocess.STDOUT,shell=True)
    if jmake: txt = "cd %s; make -j %d -f Makefile.auto" % (libdir,jmake.n)
    else: txt = "cd %s; make -f Makefile.auto" % libdir

@ -1467,7 +1489,7 @@ class POEMS(object):
    if verbose: subprocess.call(txt,shell=True)
    else:
      try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)

    if not os.path.isfile("%s/libpoems.a" % libdir) or \
          not os.path.isfile("%s/Makefile.lammps" % libdir):
@ -1504,9 +1526,10 @@ class PYTHON(object):
    libdir = dir.lib + "/python"
    if self.lammpsflag:
      subprocess.check_output("cd %s; cp Makefile.lammps.%s Makefile.lammps" %
-                         (libdir,self.lammps))
+                              (libdir,self.lammps))
    if not os.path.isfile("%s/Makefile.lammps.%s" % (libdir,self.lammps)):
-      error("Unsuccessful creation of lib/python/Makefile.lammps.%s file" % self.lammps)
+      error("Unsuccessful creation of lib/python/Makefile.lammps.%s file" %
+            self.lammps)
    else: print("Created lib/python/Makefile.lammps file")

 # QMMM lib
@ -1544,7 +1567,8 @@ class QMMM(object):
      make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps)
    make.write("%s/Makefile.auto" % libdir)

-    subprocess.check_output("cd %s; make -f Makefile.auto clean" % libdir,stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; make -f Makefile.auto clean" %
+                            libdir,stderr=subprocess.STDOUT,shell=True)
    if jmake: txt = "cd %s; make -j %d -f Makefile.auto" % (libdir,jmake.n)
    else: txt = "cd %s; make -f Makefile.auto" % libdir

@ -1553,7 +1577,7 @@ class QMMM(object):
    if verbose: subprocess.call(txt,shell=True)
    else:
      try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)
   
    if not os.path.isfile("%s/libqmmm.a" % libdir) or \
          not os.path.isfile("%s/Makefile.lammps" % libdir):
@ -1595,7 +1619,8 @@ class REAX(object):
      make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps)
    make.write("%s/Makefile.auto" % libdir)

-    commands.getoutput("cd %s; make -f Makefile.auto clean" % libdir)
+    cmd = "cd %s; make -f Makefile.auto clean" % libdir
+    subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
    if jmake: txt = "cd %s; make -j %d -f Makefile.auto" % (libdir,jmake.n)
    else: txt = "cd %s; make -f Makefile.auto" % libdir
    
@ -1604,7 +1629,7 @@ class REAX(object):
    if verbose: subprocess.call(txt,shell=True)
    else:
      try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)

    if not os.path.isfile("%s/libreax.a" % libdir) or \
          not os.path.isfile("%s/Makefile.lammps" % libdir):
@ -1641,7 +1666,8 @@ class VORONOI(object):
    if not self.install: return
    libdir = dir.lib + "/voronoi"
    cmd = "cd %s; python install.py %s" % (libdir,self.install)
-    txt = subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True).decode()
+    txt = subprocess.check_output(cmd,stderr=subprocess.STDOUT,
+                                  shell=True).decode()
    if verbose: print(txt)
    print("Created lib/voronoi library")

--- a/src/Makefile
+++ b/src/Makefile
@ -92,7 +92,7 @@ help:
 	@echo 'make package-update (pu) replace src files with updated package files'
 	@echo 'make package-overwrite   replace package files with src files'
 	@echo 'make package-diff (pd)   diff src files against package files'
-	@echo 'make package-purge       purge obsolete copies of package sources'
+	@echo 'make purge               purge obsolete copies of source files'
 	@echo ''
 	@echo 'make machine             build LAMMPS for machine'
 	@echo 'make mode=lib machine    build LAMMPS as static lib for machine'
@ -314,7 +314,7 @@ package-diff pd:
 	@echo ''
 	@for p in $(PACKUSERUC); do $(SHELL) Package.sh $$p diff; done

-package-purge: Purge.list
+purge: Purge.list
 	@echo 'Purging obsolete and auto-generated source files'
 	@for f in `grep -v '#' Purge.list` ;		\
 	    do test -f $$f && rm $$f && echo $$f || : ;		\
--- a/src/REPLICA/fix_neb.cpp
+++ b/src/REPLICA/fix_neb.cpp
@ -142,7 +142,7 @@ void FixNEB::init()
  if (count > MAXSMALLINT) error->all(FLERR,"Too many active NEB atoms");
  nebatoms = count;

-  // comm style for inter-replica exchange of coords
+  // comm mode for inter-replica exchange of coords

  if (nreplica == nprocs_universe &&
      nebatoms == atom->natoms && atom->sortfreq == 0) 
@ -392,7 +392,7 @@ void FixNEB::inter_replica_comm()
  // -----------------------------------------------------

  // single proc per replica
-  // all atoms are NEB atoms and no atom sorting is enabled
+  // all atoms are NEB atoms and no atom sorting
  // direct comm of x -> xprev and x -> xnext

  if (cmode == SINGLE_PROC_DIRECT) {
@ -414,7 +414,7 @@ void FixNEB::inter_replica_comm()
  // single proc per replica
  // but only some atoms are NEB atoms or atom sorting is enabled
  // send atom IDs and coords of only NEB atoms to prev/next proc
-  // recv proc uses atom->map() to match received coords to owned atoms
+  // recv procs use atom->map() to match received coords to owned atoms

  if (cmode == SINGLE_PROC_MAP) {
    m = 0;
--- a/src/REPLICA/prd.cpp
+++ b/src/REPLICA/prd.cpp
@ -52,6 +52,8 @@

 using namespace LAMMPS_NS;

+enum{SINGLE_PROC_DIRECT,SINGLE_PROC_MAP,MULTI_PROC};
+
 /* ---------------------------------------------------------------------- */

 PRD::PRD(LAMMPS *lmp) : Pointers(lmp) {}
@ -114,30 +116,35 @@ void PRD::command(int narg, char **arg)
  int color = me;
  MPI_Comm_split(universe->uworld,color,0,&comm_replica);

-  // equal_size_replicas = 1 if all replicas have same # of procs
-  // no longer used
+  // comm mode for inter-replica exchange of coords

-  //flag = 0;
-  //if (nreplica*nprocs == nprocs_universe) flag = 1;
-  //MPI_Allreduce(&flag,&equal_size_replicas,1,MPI_INT,MPI_MIN,
-  //              universe->uworld);
+  if (nreplica == nprocs_universe && atom->sortfreq == 0) 
+    cmode = SINGLE_PROC_DIRECT;
+  else if (nreplica == nprocs_universe) cmode = SINGLE_PROC_MAP;
+  else cmode = MULTI_PROC;

-  // workspace for inter-replica communication via gathers
+  // workspace for inter-replica communication

  natoms = atom->natoms;

-  displacements = NULL;
  tagall = NULL;
  xall = NULL;
  imageall = NULL;

-  if (nreplica != nprocs_universe) {
-    displacements = new int[nprocs];
+  if (cmode != SINGLE_PROC_DIRECT) {
    memory->create(tagall,natoms,"prd:tagall");
    memory->create(xall,natoms,3,"prd:xall");
    memory->create(imageall,natoms,"prd:imageall");
  }

+  counts = NULL;
+  displacements = NULL;
+
+  if (cmode == MULTI_PROC) {
+    memory->create(counts,nprocs,"prd:counts");
+    memory->create(displacements,nprocs,"prd:displacements");
+  }
+
  // random_select = same RNG for each replica, for multiple event selection
  // random_clock = same RNG for each replica, for clock updates
  // random_dephase = unique RNG for each replica, for dephasing
@ -238,7 +245,7 @@ void PRD::command(int narg, char **arg)
  if (domain->box_change)
    error->all(FLERR,"Cannot use PRD with a changing box");

-  // cannot use PRD with time-dependent fixes or regions or atom sorting
+  // cannot use PRD with time-dependent fixes or regions

  for (int i = 0; i < modify->nfix; i++)
    if (modify->fix[i]->time_depend)
@ -248,9 +255,6 @@ void PRD::command(int narg, char **arg)
    if (domain->regions[i]->dynamic_check())
      error->all(FLERR,"Cannot use PRD with a time-dependent region defined");

-  if (atom->sortfreq > 0)
-    error->all(FLERR,"Cannot use PRD with atom_modify sort enabled");
-
  // perform PRD simulation

  if (me_universe == 0 && universe->uscreen)
@ -433,12 +437,14 @@ void PRD::command(int narg, char **arg)
      fprintf(universe->uscreen,
              "Loop time of %g on %d procs for %d steps with " BIGINT_FORMAT
              " atoms\n",
-              timer->get_wall(Timer::TOTAL),nprocs_universe,nsteps,atom->natoms);
+              timer->get_wall(Timer::TOTAL),nprocs_universe,
+              nsteps,atom->natoms);
    if (universe->ulogfile)
      fprintf(universe->ulogfile,
              "Loop time of %g on %d procs for %d steps with " BIGINT_FORMAT
              " atoms\n",
-              timer->get_wall(Timer::TOTAL),nprocs_universe,nsteps,atom->natoms);
+              timer->get_wall(Timer::TOTAL),nprocs_universe,
+              nsteps,atom->natoms);
  }

  if (me == 0) {
@ -461,10 +467,11 @@ void PRD::command(int narg, char **arg)

  // clean up

-  delete [] displacements;
  memory->destroy(tagall);
  memory->destroy(xall);
  memory->destroy(imageall);
+  memory->destroy(counts);
+  memory->destroy(displacements);

  delete [] id_compute;
  MPI_Comm_free(&comm_replica);
@ -776,57 +783,96 @@ void PRD::log_event()

 void PRD::replicate(int ireplica)
 {
-  int nreplica = universe->nworlds;
-  int nprocs_universe = universe->nprocs;
  int i,m;

-  if (nreplica == nprocs_universe) {
-    MPI_Bcast(atom->image,atom->nlocal,MPI_INT,ireplica,comm_replica);
+  // -----------------------------------------------------
+  // 3 cases: two for single proc per replica
+  //          one for multiple procs per replica
+  // -----------------------------------------------------
+
+  // single proc per replica, no atom sorting
+  // direct bcast of image and x
+
+  if (cmode == SINGLE_PROC_DIRECT) {
    MPI_Bcast(atom->x[0],3*atom->nlocal,MPI_DOUBLE,ireplica,comm_replica);
+    MPI_Bcast(atom->image,atom->nlocal,MPI_INT,ireplica,comm_replica);
+    return;
+  }

-  } else {
-    int *counts = new int[nprocs];
-
-    if (universe->iworld == ireplica) {
-      MPI_Gather(&atom->nlocal,1,MPI_INT,counts,1,MPI_INT,0,world);
-      displacements[0] = 0;
-      for (i = 0; i < nprocs-1; i++)
-        displacements[i+1] = displacements[i] + counts[i];
-      MPI_Gatherv(atom->tag,atom->nlocal,MPI_LMP_TAGINT,
-                  tagall,counts,displacements,MPI_LMP_TAGINT,0,world);
-      MPI_Gatherv(atom->image,atom->nlocal,MPI_INT,
-                        imageall,counts,displacements,MPI_INT,0,world);
-      for (i = 0; i < nprocs; i++) counts[i] *= 3;
-      for (i = 0; i < nprocs-1; i++)
-        displacements[i+1] = displacements[i] + counts[i];
-      MPI_Gatherv(atom->x[0],3*atom->nlocal,MPI_DOUBLE,
-                        xall[0],counts,displacements,MPI_DOUBLE,0,world);
-    }
-
-    if (me == 0) {
-      MPI_Bcast(tagall,natoms,MPI_INT,ireplica,comm_replica);
-      MPI_Bcast(imageall,natoms,MPI_INT,ireplica,comm_replica);
-      MPI_Bcast(xall[0],3*natoms,MPI_DOUBLE,ireplica,comm_replica);
-    }
-
-    MPI_Bcast(tagall,natoms,MPI_INT,0,world);
-    MPI_Bcast(imageall,natoms,MPI_INT,0,world);
-    MPI_Bcast(xall[0],3*natoms,MPI_DOUBLE,0,world);
+  // single proc per replica, atom sorting is enabled
+  // bcast atom IDs, x, image via tagall, xall, imageall
+  // recv procs use atom->map() to match received info to owned atoms

+  if (cmode == SINGLE_PROC_MAP) {
    double **x = atom->x;
+    tagint *tag = atom->tag;
+    imageint *image = atom->image;
    int nlocal = atom->nlocal;

-    for (i = 0; i < natoms; i++) {
-      m = atom->map(tagall[i]);
-      if (m >= 0 && m < nlocal) {
-        x[m][0] = xall[i][0];
-        x[m][1] = xall[i][1];
-        x[m][2] = xall[i][2];
-        atom->image[m] = imageall[i];
-      }
+    if (universe->iworld == ireplica) {
+      memcpy(tagall,tag,nlocal*sizeof(tagint));
+      memcpy(xall[0],x[0],3*nlocal*sizeof(double));
+      memcpy(imageall,image,nlocal*sizeof(imageint));
    }

-    delete [] counts;
+    MPI_Bcast(tagall,natoms,MPI_INT,ireplica,comm_replica);
+    MPI_Bcast(xall[0],3*natoms,MPI_DOUBLE,ireplica,comm_replica);
+    MPI_Bcast(imageall,natoms,MPI_INT,ireplica,comm_replica);
+
+    for (i = 0; i < nlocal; i++) {
+      m = atom->map(tagall[i]);
+      x[m][0] = xall[i][0];
+      x[m][1] = xall[i][1];
+      x[m][2] = xall[i][2];
+      atom->image[m] = imageall[i];
+    }
+
+    return;
+  }
+
+  // multiple procs per replica
+  // MPI_Gather all atom IDs, x, image to root proc of ireplica
+  // bcast to root of other replicas
+  // bcast within each replica
+  // each proc extracts info for atoms it owns via atom->map()
+  // NOTE: assumes imagint and tagint are always the same size
+
+  if (universe->iworld == ireplica) {
+    MPI_Gather(&atom->nlocal,1,MPI_INT,counts,1,MPI_INT,0,world);
+    displacements[0] = 0;
+    for (i = 0; i < nprocs-1; i++)
+      displacements[i+1] = displacements[i] + counts[i];
+    MPI_Gatherv(atom->tag,atom->nlocal,MPI_LMP_TAGINT,
+                tagall,counts,displacements,MPI_LMP_TAGINT,0,world);
+    MPI_Gatherv(atom->image,atom->nlocal,MPI_LMP_TAGINT,
+                imageall,counts,displacements,MPI_LMP_TAGINT,0,world);
+    for (i = 0; i < nprocs; i++) counts[i] *= 3;
+    for (i = 0; i < nprocs-1; i++)
+      displacements[i+1] = displacements[i] + counts[i];
+    MPI_Gatherv(atom->x[0],3*atom->nlocal,MPI_DOUBLE,
+                xall[0],counts,displacements,MPI_DOUBLE,0,world);
+  }
+  
+  if (me == 0) {
+    MPI_Bcast(tagall,natoms,MPI_INT,ireplica,comm_replica);
+    MPI_Bcast(imageall,natoms,MPI_INT,ireplica,comm_replica);
+    MPI_Bcast(xall[0],3*natoms,MPI_DOUBLE,ireplica,comm_replica);
+  }
+  
+  MPI_Bcast(tagall,natoms,MPI_INT,0,world);
+  MPI_Bcast(imageall,natoms,MPI_INT,0,world);
+  MPI_Bcast(xall[0],3*natoms,MPI_DOUBLE,0,world);
+  
+  double **x = atom->x;
+  int nlocal = atom->nlocal;
+  
+  for (i = 0; i < natoms; i++) {
+    m = atom->map(tagall[i]);
+    if (m < 0 || m >= nlocal) continue;
+    x[m][0] = xall[i][0];
+    x[m][1] = xall[i][1];
+    x[m][2] = xall[i][2];
+    atom->image[m] = imageall[i];
  }
 }

--- a/src/REPLICA/prd.h
+++ b/src/REPLICA/prd.h
@ -34,7 +34,7 @@ class PRD : protected Pointers {
  int me,nprocs;
  int t_event,n_dephase,t_dephase,t_corr;
  double etol,ftol,temp_dephase;
-  int maxiter,maxeval,temp_flag,stepmode;
+  int maxiter,maxeval,temp_flag,stepmode,cmode;
  char *loop_setting,*dist_setting;

  int equal_size_replicas,natoms;
@ -46,9 +46,10 @@ class PRD : protected Pointers {
  double time_start;

  MPI_Comm comm_replica;
+  int *counts,*displacements;
  tagint *tagall;
-  int *displacements,*imageall;
  double **xall;
+  imageint *imageall;

  int ncoincident;

--- a/src/REPLICA/tad.cpp
+++ b/src/REPLICA/tad.cpp
@ -43,8 +43,6 @@
 #include "fix_store.h"
 #include "force.h"
 #include "pair.h"
-#include "random_park.h"
-#include "random_mars.h"
 #include "output.h"
 #include "dump.h"
 #include "finish.h"
@ -249,7 +247,7 @@ void TAD::command(int narg, char **arg)
  // need this line if quench() does only setup_minimal()
  // update->minimize->setup();

-  // This should work with if uncommented, but does not
+  // this should work with if statement uncommented, but does not
  // if (universe->iworld == 0) {

  fix_event->store_state_quench();
@ -399,12 +397,14 @@ void TAD::command(int narg, char **arg)
      fprintf(universe->uscreen,
              "Loop time of %g on %d procs for %d steps with " BIGINT_FORMAT
              " atoms\n",
-              timer->get_wall(Timer::TOTAL),nprocs_universe,nsteps,atom->natoms);
+              timer->get_wall(Timer::TOTAL),nprocs_universe,
+              nsteps,atom->natoms);
    if (universe->ulogfile)
      fprintf(universe->ulogfile,
              "Loop time of %g on %d procs for %d steps with " BIGINT_FORMAT
              " atoms\n",
-              timer->get_wall(Timer::TOTAL),nprocs_universe,nsteps,atom->natoms);
+              timer->get_wall(Timer::TOTAL),nprocs_universe,
+              nsteps,atom->natoms);
  }

  if ((me_universe == 0) && ulogfile_neb) fclose(ulogfile_neb);
@ -874,7 +874,7 @@ void TAD::revert_state()
 }

 /* ----------------------------------------------------------------------
-   Initialize list of possible events
+   initialize list of possible events
 ------------------------------------------------------------------------- */

 void TAD::initialize_event_list() {
@ -890,7 +890,7 @@ void TAD::initialize_event_list() {
 }

 /* ----------------------------------------------------------------------
-   Delete list of possible events
+   delete list of possible events
 ------------------------------------------------------------------------- */

 void TAD::delete_event_list() {
--- 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/USER-REAXC/pair_reax_c.cpp
+++ b/src/USER-REAXC/pair_reax_c.cpp
@ -348,7 +348,7 @@ void PairReaxC::init_style( )

  int iqeq;
  for (iqeq = 0; iqeq < modify->nfix; iqeq++)
-    if (strcmp(modify->fix[iqeq]->style,"qeq/reax") == 0) break;
+    if (strstr(modify->fix[iqeq]->style,"qeq/reax")) break;
  if (iqeq == modify->nfix && qeqflag == 1)
    error->all(FLERR,"Pair reax/c requires use of fix qeq/reax");

--- 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.cpp
+++ b/src/comm_tiled.cpp
@ -45,9 +45,6 @@ enum{LAYOUT_UNIFORM,LAYOUT_NONUNIFORM,LAYOUT_TILED};    // several files

 CommTiled::CommTiled(LAMMPS *lmp) : Comm(lmp)
 {
-  if (lmp->kokkos)
-    error->all(FLERR,"KOKKOS package does not yet support comm_style tiled");
-
  style = 1;
  layout = LAYOUT_UNIFORM;
  pbc_flag = NULL;
@ -63,9 +60,6 @@ CommTiled::CommTiled(LAMMPS *lmp) : Comm(lmp)

 CommTiled::CommTiled(LAMMPS *lmp, Comm *oldcomm) : Comm(*oldcomm)
 {
-  if (lmp->kokkos)
-    error->all(FLERR,"KOKKOS package does not yet support comm_style tiled");
-
  style = 1;
  layout = oldcomm->layout;
  Comm::copy_arrays(oldcomm);
--- a/src/comm_tiled.h
+++ b/src/comm_tiled.h
@ -26,26 +26,26 @@ class CommTiled : public Comm {

  void init();
  void setup();                        // setup comm pattern
-  void forward_comm(int dummy = 0);    // forward comm of atom coords
-  void reverse_comm();                 // reverse comm of forces
-  void exchange();                     // move atoms to new procs
-  void borders();                      // setup list of atoms to comm
+  virtual void forward_comm(int dummy = 0);    // forward comm of atom coords
+  virtual void reverse_comm();                 // reverse comm of forces
+  virtual void exchange();                     // move atoms to new procs
+  virtual void borders();                      // setup list of atoms to comm

-  void forward_comm_pair(class Pair *);    // forward comm from a Pair
-  void reverse_comm_pair(class Pair *);    // reverse comm from a Pair
+  virtual void forward_comm_pair(class Pair *);    // forward comm from a Pair
+  virtual void reverse_comm_pair(class Pair *);    // reverse comm from a Pair
  virtual void forward_comm_fix(class Fix *, int size=0);
                                                   // forward comm from a Fix
  virtual void reverse_comm_fix(class Fix *, int size=0);
                                                   // reverse comm from a Fix
  virtual void reverse_comm_fix_variable(class Fix *);
                                     // variable size reverse comm from a Fix
-  void forward_comm_compute(class Compute *);  // forward from a Compute
-  void reverse_comm_compute(class Compute *);  // reverse from a Compute
-  void forward_comm_dump(class Dump *);    // forward comm from a Dump
-  void reverse_comm_dump(class Dump *);    // reverse comm from a Dump
+  virtual void forward_comm_compute(class Compute *);  // forward from a Compute
+  virtual void reverse_comm_compute(class Compute *);  // reverse from a Compute
+  virtual void forward_comm_dump(class Dump *);    // forward comm from a Dump
+  virtual void reverse_comm_dump(class Dump *);    // reverse comm from a Dump

-  void forward_comm_array(int, double **);          // forward comm of array
-  int exchange_variable(int, double *, double *&);  // exchange on neigh stencil
+  virtual void forward_comm_array(int, double **);          // forward comm of array
+  virtual int exchange_variable(int, double *, double *&);  // exchange on neigh stencil

  void coord2proc_setup();
  int coord2proc(double *, int &, int &, int &);
@ -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_addforce.cpp
+++ b/src/fix_addforce.cpp
@ -203,13 +203,10 @@ void FixAddForce::init()
      update->whichflag == 2 && estyle == NONE)
    error->all(FLERR,"Must use variable energy with fix addforce");

-  int max_respa = 0;
-
-  if (strstr(update->integrate_style,"respa"))
-    max_respa = ((Respa *) update->integrate)->nlevels-1;
-
-  if (respa_level >= 0)
-    ilevel_respa = MIN(respa_level,max_respa);
+  if (strstr(update->integrate_style,"respa")) {
+    ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+    if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+  }
 }

 /* ---------------------------------------------------------------------- */
--- a/Show More
+++ b/Show More