Merge branch 'master' as of patch 6Jan17 into USER-DPD_kokkos

2017-01-09 13:12:02 -05:00
parent 49d644ec3a 51fa33a407
commit 26d50853e9
196 changed files with 4623 additions and 2885 deletions
--- a/doc/src/JPG/tutorial_additional_changes.png
+++ b/doc/src/JPG/tutorial_additional_changes.png
--- a/doc/src/JPG/tutorial_automated_checks.png
+++ b/doc/src/JPG/tutorial_automated_checks.png
--- a/doc/src/JPG/tutorial_automated_checks_passed.png
+++ b/doc/src/JPG/tutorial_automated_checks_passed.png
--- a/doc/src/JPG/tutorial_branch.png
+++ b/doc/src/JPG/tutorial_branch.png
--- a/doc/src/JPG/tutorial_changes_others.png
+++ b/doc/src/JPG/tutorial_changes_others.png
--- a/doc/src/JPG/tutorial_create_new_pull_request1.png
+++ b/doc/src/JPG/tutorial_create_new_pull_request1.png
--- a/doc/src/JPG/tutorial_create_new_pull_request2.png
+++ b/doc/src/JPG/tutorial_create_new_pull_request2.png
--- a/doc/src/JPG/tutorial_edits_maintainers.png
+++ b/doc/src/JPG/tutorial_edits_maintainers.png
--- a/doc/src/JPG/tutorial_fork.png
+++ b/doc/src/JPG/tutorial_fork.png
--- a/doc/src/JPG/tutorial_https_block.png
+++ b/doc/src/JPG/tutorial_https_block.png
--- a/doc/src/JPG/tutorial_new_pull_request.png
+++ b/doc/src/JPG/tutorial_new_pull_request.png
--- a/doc/src/JPG/tutorial_reverse_pull_request.png
+++ b/doc/src/JPG/tutorial_reverse_pull_request.png
--- a/doc/src/JPG/tutorial_reverse_pull_request2.png
+++ b/doc/src/JPG/tutorial_reverse_pull_request2.png
--- a/doc/src/JPG/tutorial_reverse_pull_request3.png
+++ b/doc/src/JPG/tutorial_reverse_pull_request3.png
--- a/doc/src/JPG/tutorial_reverse_pull_request4.png
+++ b/doc/src/JPG/tutorial_reverse_pull_request4.png
--- a/doc/src/JPG/tutorial_reverse_pull_request5.png
+++ b/doc/src/JPG/tutorial_reverse_pull_request5.png
--- a/doc/src/JPG/tutorial_reverse_pull_request6.png
+++ b/doc/src/JPG/tutorial_reverse_pull_request6.png
--- a/doc/src/JPG/tutorial_steve_assignee.png
+++ b/doc/src/JPG/tutorial_steve_assignee.png
--- 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="21 Dec 2016 version">
+<META NAME="docnumber" CONTENT="6 Jan 2017 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
-21 Dec 2016 version :c,h4
+6 Jan 2017 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
@ -581,7 +581,7 @@ USER-INTEL, k = KOKKOS, o = USER-OMP, t = OPT.
 "indent"_fix_indent.html,
 "langevin (k)"_fix_langevin.html,
 "lineforce"_fix_lineforce.html,
-"momentum"_fix_momentum.html,
+"momentum (k)"_fix_momentum.html,
 "move"_fix_move.html,
 "msst"_fix_msst.html,
 "neb"_fix_neb.html,
@ -702,6 +702,7 @@ package"_Section_start.html#start_3.
 "manifoldforce"_fix_manifoldforce.html,
 "meso/stationary"_fix_meso_stationary.html,
 "nve/manifold/rattle"_fix_nve_manifold_rattle.html,
 "nvk"_fix_nvk.html,
 "nvt/manifold/rattle"_fix_nvt_manifold_rattle.html,
 "nph/eff"_fix_nh_eff.html,
 "npt/eff"_fix_nh_eff.html,
--- a/doc/src/fix_colvars.txt
+++ b/doc/src/fix_colvars.txt
@ -31,21 +31,19 @@ fix abf all colvars colvars.inp tstat 1 :pre
 [Description:]
-This fix interfaces LAMMPS to a "collective variables" or "colvars"
+This fix interfaces LAMMPS to the collective variables "Colvars"
-module library which allows to calculate potentials of mean force
+library, which allows to calculate potentials of mean force
 (PMFs) for any set of colvars, using different sampling methods:
 currently implemented are the Adaptive Biasing Force (ABF) method,
 metadynamics, Steered Molecular Dynamics (SMD) and Umbrella Sampling
-(US) via a flexible harmonic restraint bias. The colvars library is
+(US) via a flexible harmonic restraint bias.
 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
 colvars library is available as "this supplementary PDF document"_PDF/colvars-refman-lammps.pdf
-A detailed discussion of the implementation of the portable collective
+The Colvars library is developed at "https://github.com/colvars/colvars"_https://github.com/colvars/colvars
-variable library is in "(Fiorin)"_#Fiorin. Additional information can
+A detailed discussion of its implementation is in "(Fiorin)"_#Fiorin.
 be found in "(Henin)"_#Henin.
 There are some example scripts for using this package with LAMMPS in the
 examples/USER/colvars directory.
@ -129,8 +127,3 @@ and tstat = NULL.
 :link(Fiorin)
 [(Fiorin)] Fiorin , Klein, Henin, Mol. Phys., DOI:10.1080/00268976.2013.813594
 :link(Henin)
 [(Henin)] Henin, Fiorin, Chipot, Klein, J. Chem. Theory Comput., 6,
 35-47 (2010)
--- a/doc/src/fix_momentum.txt
+++ b/doc/src/fix_momentum.txt
@ -7,6 +7,7 @@
 :line
 fix momentum command :h3
 fix momentum/kk command :h3
 [Syntax:]
@ -55,6 +56,29 @@ of atoms by rescaling the velocities after the momentum was removed.
 Note that the "velocity"_velocity.html command can be used to create
 initial velocities with zero aggregate linear and/or angular momentum.
 :line
 Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
 functionally the same as the corresponding style without the suffix.
 They have been optimized to run faster, depending on your available
 hardware, as discussed in "Section 5"_Section_accelerate.html
 of the manual.  The accelerated styles take the same arguments and
 should produce the same results, except for round-off and precision
 issues.
 These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
 USER-OMP and OPT packages, respectively.  They are only enabled if
 LAMMPS was built with those packages.  See the "Making
 LAMMPS"_Section_start.html#start_3 section for more info.
 You can specify the accelerated styles explicitly in your input script
 by including their suffix, or you can use the "-suffix command-line
 switch"_Section_start.html#start_7 when you invoke LAMMPS, or you can
 use the "suffix"_suffix.html command in your input script.
 See "Section 5"_Section_accelerate.html of the manual for
 more instructions on how to use the accelerated styles effectively.
 [Restart, fix_modify, output, run start/stop, minimize info:]
 No information about this fix is written to "binary restart
--- a/doc/src/fix_nvk.txt
+++ b/doc/src/fix_nvk.txt
@ -0,0 +1,71 @@
 "LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
 :link(lws,http://lammps.sandia.gov)
 :link(ld,Manual.html)
 :link(lc,Section_commands.html#comm)
 :line
 fix nvk command :h3
 [Syntax:]
 fix ID group-ID nvk :pre
 ID, group-ID are documented in "fix"_fix.html command
 nvk = style name of this fix command :ul
 [Examples:]
 fix 1 all nvk :pre
 [Description:]
 Perform constant kinetic energy integration using the Gaussian
 thermostat to update position and velocity for atoms in the group each
 timestep.  V is volume; K is kinetic energy. This creates a system
 trajectory consistent with the isokinetic ensemble.
 The equations of motion used are those of Minary et al in
 "(Minary)"_#nvk-Minary, a variant of those initially given by Zhang in 
 "(Zhang)"_#nvk-Zhang.
 The kinetic energy will be held constant at its value given when fix
 nvk is initiated. If a different kinetic energy is desired, the
 "velocity"_velocity.html command should be used to change the kinetic
 energy prior to this fix.
 :line
 [Restart, fix_modify, output, run start/stop, minimize info:]
 No information about this fix is written to "binary restart
 files"_restart.html.  None of the "fix_modify"_fix_modify.html options
 are relevant to this fix.  No global or per-atom quantities are stored
 by this fix for access by various "output
 commands"_Section_howto.html#howto_15.  No parameter of this fix can
 be used with the {start/stop} keywords of the "run"_run.html command.
 This fix is not invoked during "energy minimization"_minimize.html.
 [Restrictions:]
 The Gaussian thermostat only works when it is applied to all atoms in
 the simulation box. Therefore, the group must be set to all.
 This fix has not yet been implemented to work with the RESPA integrator.
 This fix is part of the USER-MISC package.  It is only enabled if LAMMPS
 was built with that package.  See the "Making
 LAMMPS"_Section_start.html#start_3 section for more info.
 [Related commands:] none
 [Default:] none
 :line
 :link(nvk-Minary)
 [(Minary)] Minary, Martyna, and Tuckerman, J Chem Phys, 18, 2510 (2003).
 :link(nvk-Zhang)
 [(Zhang)] Zhang, J Chem Phys, 106, 6102 (1997).
--- a/doc/src/fixes.txt
+++ b/doc/src/fixes.txt
@ -90,6 +90,7 @@ Fixes :h1
   fix_nve_noforce
   fix_nve_sphere
   fix_nve_tri
   fix_nvk
   fix_nvt_asphere
   fix_nvt_body
   fix_nvt_manifold_rattle
--- a/doc/src/lammps.book
+++ b/doc/src/lammps.book
@ -214,6 +214,7 @@ fix_nve_manifold_rattle.html
 fix_nve_noforce.html
 fix_nve_sphere.html
 fix_nve_tri.html
 fix_nvk.html
 fix_nvt_asphere.html
 fix_nvt_body.html
 fix_nvt_manifold_rattle.html
--- a/doc/src/tutorial_github.txt
+++ b/doc/src/tutorial_github.txt
@ -11,10 +11,22 @@ LAMMPS GitHub tutorial :h3
 :line
-This document briefly describes how to use GitHub to merge changes you
+This document describes the process of how to use GitHub to integrate
-make into LAMMPS, using GitHub. It assumes that you are familiar with
+changes or additions you have made to LAMMPS into the official LAMMPS
-git. You may want to have a look at the "Git
+distribution.  It uses the process of updating this very tutorial as
-book"_http://git-scm.com/book/ to reacquaint yourself.
+an example to describe the individual steps and options.  You need to
 be familiar with git and you may want to have a look at the
 "Git book"_http://git-scm.com/book/ to reacquaint yourself with some
 of the more advanced git features used below.
 As of fall 2016, submitting contributions to LAMMPS via pull requests
 on GitHub is the preferred option for integrating contributed features
 or improvements to LAMMPS, as it significantly reduces the amount of
 work required by the LAMMPS developers. Consequently, creating a pull
 request will increase your chances to have your contribution included
 and will reduce the time until the integration is complete. For more
 information on the requirements to have your code included into LAMMPS
 please see "Section 10.15"_Section_modify.html#mod_15
 :line
@ -30,106 +42,121 @@ username or e-mail address and password.
 [Forking the repository]
-To get changes into LAMMPS, you need to first fork the repository. At
+To get changes into LAMMPS, you need to first fork the `lammps/lammps`
-the time of writing, LAMMPS-ICMS is the preferred fork. Go to "LAMMPS
+repository on GitHub. At the time of writing, {master} is the preferred
-on GitHub"_https://github.com/lammps/lammps and make sure branch is
+target branch. Thus go to "LAMMPS on GitHub"_https://github.com/lammps/lammps
-set to "lammps-icms", see the figure below.
+and make sure branch is set to "master", as shown in the figure below.
 :c,image(JPG/tutorial_branch.png)
-Now, click on fork in the top right corner:
+If it is not, use the button to change it to {master}. Once it is, use the
 fork button to create a fork.
 :c,image(JPG/tutorial_fork.png)
-This will create your own fork of the LAMMPS repository. You can make
+
-changes in this fork and later file {pull requests} to allow the
+This will create a fork (which is essentially a copy, but uses less
-upstream repository to merge changes from your own fork into the one
+resources) of the LAMMPS repository under your own GitHub account. You
-we just forked from. At the same time, you can set things up, so you
+can make changes in this fork and later file {pull requests} to allow
-can include changes from upstream into your repository.
+the upstream repository to merge changes from your own fork into the one
 we just forked from (or others that were forked from the same repository).
 At the same time, you can set things up, so you can include changes from
 upstream into your repository and thus keep it in sync with the ongoing
 LAMMPS development.
 :line
 [Adding changes to your own fork]
-Before adding changes, it is better to first create a new branch that
+Additions to the upstream version of LAMMPS are handled using {feature
-will contain these changes, a so-called feature branch.
+branches}.  For every new feature, a so-called feature branch is
 created, which contains only those modification relevant to one specific
 feature. For example, adding a single fix would consist of creating a
 branch with only the fix header and source file and nothing else.  It is
 explained in more detail here: "feature branch
 workflow"_https://www.atlassian.com/git/tutorials/comparing-workflows/feature-branch-workflow.
 [Feature branches]
-Since LAMMPS is such a big project and most user contributions come in
+First of all, create a clone of your version on github on your local
-small portions, the most ideal workflow for LAMMPS is the so-called
+machine via HTTPS:
 "Feature branch" workflow. It is explained in great detail here:
 "feature branch
 workflow"_https://www.atlassian.com/git/tutorials/comparing-workflows/feature-branch-workflow.
-The idea is that every new feature for LAMMPS gets its own
+  $ git clone https://github.com/<your user name>/lammps.git <some name> :pre
 branch. This way, it is fairly painless to incorporate new features
 into the upstream repository. I will explain briefly here how to do
 it. In this feature branch, I will add a USER-package.
-I assume that git is installed on the local machine and you know how
+or, if you have set up your GitHub account for using SSH keys, via SSH:
 to use a command line.
-First of all, you need to clone your own fork of LAMMPS:
+  $ git clone git@github.com:<your user name>/lammps.git :pre
-  $ git clone https://github.com/<your user name>/lammps.git :pre
+You can find the proper URL by clicking the "Clone or download"-button:
 You can find the proper url to the right of the "HTTPS" block, see figure.
 :c,image(JPG/tutorial_https_block.png)
 The above command copies ("clones") the git repository to your local
-machine. You can use this local clone to make changes and test them
+machine to a directory with the name you chose. If none is given, it will
-without interfering with the repository on github. First, however, it
+default to "lammps". Typical names are "mylammps" or something similar.
 is recommended to make a new branch for a particular feature you would
 like added to LAMMPS. In this example, I will try adding a new
 USER-package called USER-MANIFOLD.
-To create a new branch, run the following git command in your repository:
+You can use this local clone to make changes and
 test them without interfering with the repository on Github.
-$ git checkout -b add-user-manifold :pre
+To pull changes from upstream into this copy, you can go to the directory
 and use git pull:
-The name of this new branch is "add-user-manifold" in my case. Just
+  $ cd mylammps
-name it after something that resembles the feature you want added to
+  $ git checkout master
-LAMMPS.
+  $ git pull https://github.com/lammps/lammps :pre
-Now that you've changed branches, you can edit the files as you see
+You can also add this URL as a remote:
-fit, add new files, and commit as much as you would like. Just
+
-remember that if halfway you decide to add another, unrelated feature,
+  $ git remote add lammps_upstream https://www.github.com/lammps/lammps :pre
-you should switch branches!
+
 At this point, you typically make a feature branch from the updated master
 branch for the feature you want to work on. This tutorial contains the
 workflow that updated this tutorial, and hence we will call the branch
 "github-tutorial-update":
 $ git checkout -b github-tutorial-update master :pre
 Now that we have changed branches, we can make our changes to our local
 repository. Just remember that if you want to start working on another,
 unrelated feature, you should switch branches!
 [After changes are made]
 After everything is done, add the files to the branch and commit them:
-  $ git add src/USER-MANIFOLD examples/USER/manifold/
+ $ git add doc/src/tutorial_github.txt
-  $ git add doc/fix_nv\{t,e\}_manifold_rattle.txt
+ $ git add doc/src/JPG/tutorial*.png :pre
  $ git add doc/fix_manifoldforce.txt doc/user_manifolds.txt :pre
-After the files are added, the change should be comitted:
+IMPORTANT NOTE: Do not use {git commit -a} (or {git add -A}).  The -a
 flag (or -A flag) will automatically include _all_ modified or new files
 and that is rarely the behavior you want.  It can easily lead to
 accidentally adding unrelated and unwanted changes into the repository.
 Instead it is preferable to explicitly use {git add}, {git rm}, {git mv}
 for adding, removing, renaming individual files, respectively, and then
 {git commit} to finalize the commit.  Carefully check all pending
 changes with {git status} before committing them.  If you find doing
 this on the command line too tedious, consider using a GUI, for example
 the one included in git distributions written in Tk, i.e. use {git gui}
 (on some Linux distributions it may be required to install an additional
 package to use it).
-  $ git commit -m 'Added user-manifold package' :pre
+After adding all files, the change set can be committed with some
 useful message that explains the change.
-The "-m" switch is used to add a message to the commit. Use this to
+  $ git commit -m 'Finally updated the github tutorial' :pre
 indicate what type of change was commited.
 [Wisdom by Axel]
 {"Do not use "git commit -a".  the -a flag will automatically include
 *all* modified or new files.  mercurial does that and it find it
 hugely annoying and often leading to accidental commits of files you
 don't want.  use git add, git rm, git mv for adding, removing,
 renaming and then git commit to finalize the commit.  personally, i
 find it very convenient to use the bundled gui for commits, i.e. git
 gui.  typically, i will do git add and other operations, but then
 verify and review them with git gui.  git gui also allows to do
 line-by-line unstaging and other convenient operations."}
 After the commit, the changes can be pushed to the same branch on GitHub:
 $ git push :pre
 Git will ask you for your user name and password on GitHub if you have
-not configured anything. If you correctly type your user name and
+not configured anything. If your local branch is not present on Github yet,
-password, the change should be added to your fork on GitHub.
+it will ask you to add it by running
  $ git push --set-upstream origin github-tutorial-update :pre
 If you correctly type your user name and
 password, the feature branch should be added to your fork on GitHub.
 If you want to make really sure you push to the right repository
 (which is good practice), you can provide it explicitly:
@ -140,16 +167,20 @@ or using an explicit URL:
 $ git push git@github.com:Pakketeretet2/lammps.git :pre
-After that, you can file a new pull request based on this
+:line
 branch. GitHub will now look like this:
-:c,image(JPG/tutorial_pull_request_feature_branch1.png)
+[Filing a pull request]
 Up to this point in the tutorial, all changes were to {your} clones of
 LAMMPS.  Eventually, however, you want this feature to be included into
 the official LAMMPS version.  To do this, you will want to file a pull
 request by clicking on the "New pull request" button:
 :c,image(JPG/tutorial_new_pull_request.png)
 Make sure that the current branch is set to the correct one, which, in
-this case, is "add-user-manifold". Now click "New pull request". If
+this case, is "github-tutorial-update". If done correctly, the only
-done correctly, the only changes you will see are those that were made
+changes you will see are those that were made on this branch.
 on this branch, so in my case, I will see nothing related to
 $\mathrm{pair\_dzugatov}.$
 This will open up a new window that lists changes made to the
 repository. If you are just adding new files, there is not much to do,
@ -158,36 +189,159 @@ changes in existing files. If all changes can automatically be merged,
 green text at the top will say so and you can click the "Create pull
 request" button, see image.
-:c,image(JPG/tutorial_pull_request2.png)
+:c,image(JPG/tutorial_create_new_pull_request1.png)
-After this you have to specify a short title and a comment with
+Before creating the pull request, make sure the short title is accurate
-details about your pull request. I guess here you write what your
+and add a comment with details about your pull request.  Here you write
-modifications do and why they should be incorporated upstream. After
+what your modifications do and why they should be incorporated upstream.
 that, click the "Create pull request" button, see image below.
-:c,image(JPG/tutorial_pull_request3.png)
+Note the checkbox that says "Allow edits from maintainers".
 This is checked by default checkbox (although in my version of Firefox, only the checkmark is visible):
-Now just write some nice comments, click "Comment", and that is it. It
+:c,image(JPG/tutorial_edits_maintainers.png)
 is now up to the maintainer(s) of the upstream repository to
 incorporate the changes into the repository and to close the pull
 request.
-:c,image(JPG/tutorial_pull_request4.png)
+If it is checked, maintainers can immediately add their own edits to the
 pull request.  This helps the inclusion of your branch significantly, as
 simple/trivial changes can be added directly to your pull request branch
 by the LAMMPS maintainers.  The alternative would be that they make
 changes on their own version of the branch and file a reverse pull
 request to you.  Just leave this box checked unless you have a very good
 reason not to.
 Now just write some nice comments and click on "Create pull request".
 :c,image(JPG/tutorial_create_new_pull_request2.png)
 :line
 [After filing a pull request]
 NOTE: When you submit a pull request (or ask for a pull request) for the
 first time, you will receive an invitation to become a LAMMPS project
 collaborator. Please accept this invite as being a collaborator will
 simplify certain administrative tasks and will probably speed up the
 merging of your feature, too.
 You will notice that after filing the pull request, some checks are
 performed automatically:
 :c,image(JPG/tutorial_automated_checks.png)
 If all is fine, you will see this:
 :c,image(JPG/tutorial_automated_checks_passed.png)
 If any of the checks are failing, your pull request will not be
 processed, as your changes may break compilation for certain
 configurations or may not merge cleanly. It is your responsibility
 to remove the reason(s) for the failed test(s). If you need help
 with this, please contact the LAMMPS developers by adding a comment
 explaining your problems with resolving the failed tests.
 A few further interesting things (can) happen to pull requests before
 they are included.
 [Additional changes]
-Before the pull request is accepted, any additional changes you push
+First of all, any additional changes you push into your branch in your
-into your repository will automatically become part of the pull
+repository will automatically become part of the pull request:
-request.
+
 :c,image(JPG/tutorial_additional_changes.png)
 This means you can add changes that should be part of the feature after
 filing the pull request, which is useful in case you have forgotten
 them, or if a developer has requested that something needs to be changed
 before the feature can be accepted into the official LAMMPS version.
 After each push, the automated checks are run again.
 [Assignees]
 There is an assignee label for pull requests. If the request has not
 been reviewed by any developer yet, it is not assigned to anyone. After
 revision, a developer can choose to assign it to either a) you, b) a
 LAMMPS developer (including him/herself) or c) Steve Plimpton (sjplimp).
 Case a) happens if changes are required on your part :ulb,l
 Case b) means that at the moment, it is being tested and reviewed by a
 LAMMPS developer with the expectation that some changes would be required.
 After the review, the developer can choose to implement changes directly
 or suggest them to you. :l
 Case c) means that the pull request has been assigned to the lead
 developer Steve Plimpton and means it is considered ready for merging. :ule,l
 In this case, Axel assigned the tutorial to Steve:
 :c,image(JPG/tutorial_steve_assignee.png)
 [Edits from LAMMPS maintainers]
 If you allowed edits from maintainers (the default), any LAMMPS
 maintainer can add changes to your pull request.  In this case, both
 Axel and Richard made changes to the tutorial:
 :c,image(JPG/tutorial_changes_others.png)
 [Reverse pull requests]
 Sometimes, however, you might not feel comfortable having other people
 push changes into your own branch, or maybe the maintainers are not sure
 their idea was the right one.  In such a case, they can make changes,
 reassign you as the assignee, and file a "reverse pull request", i.e.
 file a pull request in your GitHub repository to include changes in the
 branch, that you have submitted as a pull request yourself.  In that
 case, you can choose to merge their changes back into your branch,
 possibly make additional changes or corrections and proceed from there.
 It looks something like this:
 :c,image(JPG/tutorial_reverse_pull_request.png)
 For some reason, the highlighted button didn't work in my case, but I
 can go to my own repository and merge the pull request from there:
 :c,image(JPG/tutorial_reverse_pull_request2.png)
 Be sure to check the changes to see if you agree with them by clicking
 on the tab button:
 :c,image(JPG/tutorial_reverse_pull_request3.png)
 In this case, most of it is changes in the markup and a short rewrite of
 Axel's explanation of the "git gui" and "git add" commands.
 :c,image(JPG/tutorial_reverse_pull_request4.png)
 Because the changes are OK with us, we are going to merge by clicking on
 "Merge pull request".  After a merge it looks like this:
 :c,image(JPG/tutorial_reverse_pull_request5.png)
 Now, since in the meantime our local text for the tutorial also changed,
 we need to pull Axel's change back into our branch, and merge them:
 $ git add tutorial_github.txt
 $ git add JPG/tutorial_reverse_pull_request*.png
 $ git commit -m "Updated text and images on reverse pull requests"
 $ git pull :pre
 In this case, the merge was painless because git could auto-merge:
 :c,image(JPG/tutorial_reverse_pull_request6.png)
 With Axel's changes merged in and some final text updates, our feature
 branch is now perfect as far as we are concerned, so we are going to
 commit and push again:
 $ git add tutorial_github.txt
 $ git add JPG/tutorial_reverse_pull_request6.png
 $ git commit -m "Merged Axel's suggestions and updated text"
 $ git push git@github.com:Pakketeretet2/lammps :pre
 :line
 [After a merge]
-When everything is fine the feature branch is merged into the LAMMPS
+When everything is fine, the feature branch is merged into the master branch.
 repositories:
 :c,image(JPG/tutorial_merged.png)
@ -198,17 +352,29 @@ It is in principle safe to delete them from your own fork. This helps
 keep it a bit more tidy. Note that you first have to switch to another
 branch!
-$ git checkout lammps-icms
+$ git checkout master
-$ git pull lammps-icms
+$ git pull master
-$ git branch -d add-user-manifold :pre
+$ git branch -d github-tutorial-update :pre
 If you do not pull first, it is not really a problem but git will warn
 you at the next statement that you are deleting a local branch that
 was not yet fully merged into HEAD. This is because git does not yet
-know your branch just got merged into lammps-icms upstream. If you
+know your branch just got merged into LAMMPS upstream. If you
 first delete and then pull, everything should still be fine.
 Finally, if you delete the branch locally, you might want to push this
 to your remote(s) as well:
-$ git push origin :add-user-manifold :pre
+$ git push origin :github-tutorial-update :pre
 [Recent changes in the workflow]
 Some changes to the workflow are not captured in this tutorial.  For
 example, in addition to the master branch, to which all new features
 should be submitted, there is now also an "unstable" and a "stable"
 branch; these have the same content as "master", but are only updated
 after a patch release or stable release was made.
 Furthermore, the naming of the patches now follow the pattern
 "patch_<Day><Month><Year>" to simplify comparisons between releases.
 Finally, all patches and submissions are subject to automatic testing
 and code checks to make sure they at the very least compile.
--- a/lib/colvars/Makefile.colvars
+++ b/lib/colvars/Makefile.colvars
@ -0,0 +1,119 @@
 # library build -*- makefile -*- for colvars module
 # which file will be copied to Makefile.lammps
 EXTRAMAKE = Makefile.lammps.empty
 # ------ SETTINGS ------
 CXX =		g++
 CXXFLAGS =	-O2 -g -Wall -fPIC -funroll-loops # -DCOLVARS_DEBUG
 ARCHIVE =	ar
 ARCHFLAG =	-rscv
 SHELL =		/bin/sh
 # ------ DEFINITIONS ------
 SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias_alb.cpp colvarbias.cpp  \
 colvarbias_histogram.cpp colvarbias_meta.cpp colvarbias_restraint.cpp      \
 colvarcomp_angles.cpp colvarcomp_coordnums.cpp colvarcomp.cpp              \
 colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp   \
 colvardeps.cpp colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp  \
 colvarscript.cpp colvartypes.cpp colvarvalue.cpp
 LIB = libcolvars.a
 OBJ = $(SRC:.cpp=.o)
 EXE = #colvars_standalone
 # ------ MAKE PROCEDURE ------
 default: $(LIB) $(EXE) Makefile.lammps
 Makefile.lammps:
 	@cp $(EXTRAMAKE) Makefile.lammps
 $(LIB):	$(OBJ)
 	$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
 colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB)
 	$(CXX) -o $@ $(CXXFLAGS) $^
 # ------ MAKE FLAGS ------
 .SUFFIXES:
 .SUFFIXES: .cpp .o
 .PHONY: default clean
 # ------ COMPILE RULES ------
 .cpp.o:
 	$(CXX) $(CXXFLAGS) -c $<
 # ------ DEPENDENCIES ------
 #
 colvaratoms.o: colvaratoms.cpp colvarmodule.h colvartypes.h colvarproxy.h \
 colvarvalue.h colvarparse.h colvardeps.h colvaratoms.h
 colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h colvartypes.h \
 colvarproxy.h colvarvalue.h colvar.h colvarparse.h colvardeps.h \
 colvarbias_abf.h colvarbias.h colvargrid.h
 colvarbias_alb.o: colvarbias_alb.cpp colvarmodule.h colvartypes.h \
 colvarproxy.h colvarvalue.h colvarbias_alb.h colvar.h colvarparse.h \
 colvardeps.h colvarbias_restraint.h colvarbias.h
 colvarbias.o: colvarbias.cpp colvarmodule.h colvartypes.h colvarproxy.h \
 colvarvalue.h colvarbias.h colvar.h colvarparse.h colvardeps.h
 colvarbias_histogram.o: colvarbias_histogram.cpp colvarmodule.h \
 colvartypes.h colvarproxy.h colvarvalue.h colvar.h colvarparse.h \
 colvardeps.h colvarbias_histogram.h colvarbias.h colvargrid.h
 colvarbias_meta.o: colvarbias_meta.cpp colvar.h colvarmodule.h \
 colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
 colvarbias_meta.h colvarbias.h colvargrid.h
 colvarbias_restraint.o: colvarbias_restraint.cpp colvarmodule.h \
 colvartypes.h colvarproxy.h colvarvalue.h colvarbias_restraint.h \
 colvarbias.h colvar.h colvarparse.h colvardeps.h
 colvarcomp_angles.o: colvarcomp_angles.cpp colvarmodule.h colvartypes.h \
 colvarproxy.h colvarvalue.h colvar.h colvarparse.h colvardeps.h \
 colvarcomp.h colvaratoms.h
 colvarcomp_coordnums.o: colvarcomp_coordnums.cpp colvarmodule.h \
 colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
 colvaratoms.h colvar.h colvarcomp.h
 colvarcomp.o: colvarcomp.cpp colvarmodule.h colvartypes.h colvarproxy.h \
 colvarvalue.h colvar.h colvarparse.h colvardeps.h colvarcomp.h \
 colvaratoms.h
 colvarcomp_distances.o: colvarcomp_distances.cpp colvarmodule.h \
 colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
 colvar.h colvarcomp.h colvaratoms.h
 colvarcomp_protein.o: colvarcomp_protein.cpp colvarmodule.h colvartypes.h \
 colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h colvar.h \
 colvarcomp.h colvaratoms.h
 colvarcomp_rotations.o: colvarcomp_rotations.cpp colvarmodule.h \
 colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
 colvar.h colvarcomp.h colvaratoms.h
 colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \
 colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
 colvaratoms.h colvarscript.h colvarbias.h
 colvardeps.o: colvardeps.cpp colvardeps.h colvarmodule.h colvartypes.h \
 colvarproxy.h colvarvalue.h colvarparse.h
 colvargrid.o: colvargrid.cpp colvarmodule.h colvartypes.h colvarproxy.h \
 colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
 colvaratoms.h colvargrid.h
 colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h \
 colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h colvar.h \
 colvarbias.h colvarbias_abf.h colvargrid.h colvarbias_alb.h \
 colvarbias_restraint.h colvarbias_histogram.h colvarbias_meta.h \
 colvarscript.h
 colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \
 colvarvalue.h colvarparse.h
 colvarscript.o: colvarscript.cpp colvarscript.h colvarmodule.h \
 colvartypes.h colvarproxy.h colvarvalue.h colvarbias.h colvar.h \
 colvarparse.h colvardeps.h
 colvartypes.o: colvartypes.cpp colvarmodule.h colvartypes.h colvarproxy.h \
 colvarvalue.h colvarparse.h
 colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \
 colvarvalue.h
 # ------ CLEAN ------
 clean:
 	-rm *.o *~ $(LIB)
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include "colvarmodule.h"
 #include "colvarvalue.h"
 #include "colvarparse.h"
@ -163,38 +170,38 @@ colvar::colvar(std::string const &conf)
    }
    feature_states[f_cv_homogeneous]->enabled = homogeneous;
  }
-  // Colvar is deemed periodic iff:
+
  // Colvar is deemed periodic if:
  // - it is homogeneous
  // - all cvcs are periodic
  // - all cvcs have the same period
-
+  if (cvcs[0]->b_periodic) { // TODO make this a CVC feature
-  b_periodic = cvcs[0]->b_periodic && is_enabled(f_cv_homogeneous);
+    bool b_periodic = true;
-  period = cvcs[0]->period;
+    period = cvcs[0]->period;
-  for (i = 1; i < cvcs.size(); i++) {
+    for (i = 1; i < cvcs.size(); i++) {
-    if (!cvcs[i]->b_periodic || cvcs[i]->period != period) {
+      if (!cvcs[i]->b_periodic || cvcs[i]->period != period) {
-      b_periodic = false;
+        b_periodic = false;
-      period = 0.0;
+        period = 0.0;
        cvm::log("Warning: although one component is periodic, this colvar will "
                 "not be treated as periodic, either because the exponent is not "
                 "1, or because components of different periodicity are defined.  "
                 "Make sure that you know what you are doing!");
      }
    }
    feature_states[f_cv_periodic]->enabled = b_periodic;
  }
  feature_states[f_cv_periodic]->enabled = b_periodic;
  // check that cvcs are compatible
  for (i = 0; i < cvcs.size(); i++) {
    if ((cvcs[i])->b_periodic && !b_periodic) {
        cvm::log("Warning: although this component is periodic, the colvar will "
                  "not be treated as periodic, either because the exponent is not "
                  "1, or because multiple components are present. Make sure that "
                  "you know what you are doing!");
    }
    // components may have different types only for scripted functions
    if (!is_enabled(f_cv_scripted) && (colvarvalue::check_types(cvcs[i]->value(),
-                                                           cvcs[0]->value())) ) {
+                                                                cvcs[0]->value())) ) {
      cvm::error("ERROR: you are definining this collective variable "
                 "by using components of different types. "
                 "You must use the same type in order to "
-                 " sum them together.\n", INPUT_ERROR);
+                 "sum them together.\n", INPUT_ERROR);
      return;
    }
  }
@ -207,16 +214,15 @@ colvar::colvar(std::string const &conf)
  // at this point, the colvar's type is defined
  f.type(value());
  f_accumulated.type(value());
-  fb.type(value());
+
  reset_bias_force();
  get_keyval(conf, "width", width, 1.0);
  if (width <= 0.0) {
    cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR);
    return;
  }
  // NOTE: not porting wall stuff to new deps, as this will change to a separate bias
  // the grid functions will wait a little as well
  lower_boundary.type(value());
  lower_wall.type(value());
@ -308,6 +314,9 @@ colvar::colvar(std::string const &conf)
      enable(f_cv_extended_Lagrangian);
      provide(f_cv_Langevin);
      // The extended mass will apply forces
      enable(f_cv_gradient);
      xr.type(value());
      vr.type(value());
      fr.type(value());
@ -400,6 +409,9 @@ colvar::colvar(std::string const &conf)
  f_old.type(value());
  f_old.reset();
  x_restart.type(value());
  after_restart = false;
  if (cvm::b_analysis)
    parse_analysis(conf);
@ -761,9 +773,13 @@ int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
    return error_code;
  }
  if (cvm::step_relative() > 0) {
    // Total force depends on Jacobian derivative from previous timestep
    error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
  }
  // atom coordinates are updated by the next line
  error_code |= calc_cvc_values(first_cvc, num_cvcs);
  error_code |= calc_cvc_gradients(first_cvc, num_cvcs);
  error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
  error_code |= calc_cvc_Jacobians(first_cvc, num_cvcs);
  if (cvm::debug())
@ -780,9 +796,12 @@ int colvar::collect_cvc_data()
  int error_code = COLVARS_OK;
  if (cvm::step_relative() > 0) {
    // Total force depends on Jacobian derivative from previous timestep
    error_code |= collect_cvc_total_forces();
  }
  error_code |= collect_cvc_values();
  error_code |= collect_cvc_gradients();
  error_code |= collect_cvc_total_forces();
  error_code |= collect_cvc_Jacobians();
  error_code |= calc_colvar_properties();
@ -872,6 +891,22 @@ int colvar::collect_cvc_values()
    cvm::log("Colvar \""+this->name+"\" has value "+
              cvm::to_str(x, cvm::cv_width, cvm::cv_prec)+".\n");
  if (after_restart) {
    after_restart = false;
    if (cvm::proxy->simulation_running()) {
      cvm::real const jump2 = dist2(x, x_restart) / (width*width);
      if (jump2 > 0.25) {
        cvm::error("Error: the calculated value of colvar \""+name+
                   "\":\n"+cvm::to_str(x)+"\n differs greatly from the value "
                   "last read from the state file:\n"+cvm::to_str(x_restart)+
                   "\nPossible causes are changes in configuration, "
                   "wrong state file, or how PBC wrapping is handled.\n",
                   INPUT_ERROR);
        return INPUT_ERROR;
      }
    }
  }
  return COLVARS_OK;
 }
@ -979,22 +1014,17 @@ int colvar::calc_cvc_total_force(int first_cvc, size_t num_cvcs)
    if (cvm::debug())
      cvm::log("Calculating total force of colvar \""+this->name+"\".\n");
-    // if (!tasks[task_extended_lagrangian] && (cvm::step_relative() > 0)) {
+    cvm::increase_depth();
   // Disabled check to allow for explicit total force calculation
    // even with extended Lagrangian
-    if (cvm::step_relative() > 0) {
+    for (i = first_cvc, cvc_count = 0;
-      cvm::increase_depth();
+        (i < cvcs.size()) && (cvc_count < cvc_max_count);
-      // get from the cvcs the total forces from the PREVIOUS step
+        i++) {
-      for (i = first_cvc, cvc_count = 0;
+      if (!cvcs[i]->is_enabled()) continue;
-          (i < cvcs.size()) && (cvc_count < cvc_max_count);
+      cvc_count++;
-          i++) {
+      (cvcs[i])->calc_force_invgrads();
        if (!cvcs[i]->is_enabled()) continue;
        cvc_count++;
        (cvcs[i])->calc_force_invgrads();
      }
      cvm::decrease_depth();
    }
    cvm::decrease_depth();
    if (cvm::debug())
      cvm::log("Done calculating total force of colvar \""+this->name+"\".\n");
@ -1013,6 +1043,11 @@ int colvar::collect_cvc_total_forces()
      // get from the cvcs the total forces from the PREVIOUS step
      for (size_t i = 0; i < cvcs.size();  i++) {
        if (!cvcs[i]->is_enabled()) continue;
            if (cvm::debug())
            cvm::log("Colvar component no. "+cvm::to_str(i+1)+
                " within colvar \""+this->name+"\" has total force "+
                cvm::to_str((cvcs[i])->total_force(),
                cvm::cv_width, cvm::cv_prec)+".\n");
        // linear combination is assumed
        ft += (cvcs[i])->total_force() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
      }
@ -1056,6 +1091,11 @@ int colvar::collect_cvc_Jacobians()
    fj.reset();
    for (size_t i = 0; i < cvcs.size(); i++) {
      if (!cvcs[i]->is_enabled()) continue;
        if (cvm::debug())
          cvm::log("Colvar component no. "+cvm::to_str(i+1)+
            " within colvar \""+this->name+"\" has Jacobian derivative"+
            cvm::to_str((cvcs[i])->Jacobian_derivative(),
            cvm::cv_width, cvm::cv_prec)+".\n");
      // linear combination is assumed
      fj += (cvcs[i])->Jacobian_derivative() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
    }
@ -1085,7 +1125,7 @@ int colvar::calc_colvar_properties()
    // TODO: put it in the restart information
    if (cvm::step_relative() == 0) {
      xr = x;
-      vr = 0.0; // (already 0; added for clarity)
+      vr.reset(); // (already 0; added for clarity)
    }
    // report the restraint center as "value"
@ -1128,15 +1168,16 @@ cvm::real colvar::update_forces_energy()
  if (is_enabled(f_cv_Jacobian)) {
    // the instantaneous Jacobian force was not included in the reported total force;
    // instead, it is subtracted from the applied force (silent Jacobian correction)
    // This requires the Jacobian term for the *current* timestep
    if (is_enabled(f_cv_hide_Jacobian))
      f -= fj;
  }
-  if (is_enabled(f_cv_lower_wall) || is_enabled(f_cv_upper_wall)) {
+  // Wall force
  colvarvalue fw(x);
  fw.reset();
-    // Wall force
+  if (is_enabled(f_cv_lower_wall) || is_enabled(f_cv_upper_wall)) {
    colvarvalue fw(x);
    fw.reset();
    if (cvm::debug())
      cvm::log("Calculating wall forces for colvar \""+this->name+"\".\n");
@ -1144,12 +1185,11 @@ cvm::real colvar::update_forces_energy()
    // For a periodic colvar, both walls may be applicable at the same time
    // in which case we pick the closer one
    if ( (!is_enabled(f_cv_upper_wall)) ||
-         (this->dist2(x_reported, lower_wall) < this->dist2(x_reported, upper_wall)) ) {
+         (this->dist2(x, lower_wall) < this->dist2(x, upper_wall)) ) {
-      cvm::real const grad = this->dist2_lgrad(x_reported, lower_wall);
+      cvm::real const grad = this->dist2_lgrad(x, lower_wall);
      if (grad < 0.0) {
        fw = -0.5 * lower_wall_k * grad;
        f += fw;
        if (cvm::debug())
          cvm::log("Applying a lower wall force("+
                    cvm::to_str(fw)+") to \""+this->name+"\".\n");
@ -1157,10 +1197,9 @@ cvm::real colvar::update_forces_energy()
    } else {
-      cvm::real const grad = this->dist2_lgrad(x_reported, upper_wall);
+      cvm::real const grad = this->dist2_lgrad(x, upper_wall);
      if (grad > 0.0) {
        fw = -0.5 * upper_wall_k * grad;
        f += fw;
        if (cvm::debug())
          cvm::log("Applying an upper wall force("+
                    cvm::to_str(fw)+") to \""+this->name+"\".\n");
@ -1168,17 +1207,26 @@ cvm::real colvar::update_forces_energy()
    }
  }
  // At this point f is the force f from external biases that will be applied to the
  // extended variable if there is one
  if (is_enabled(f_cv_extended_Lagrangian)) {
    if (cvm::debug()) {
      cvm::log("Updating extended-Lagrangian degrees of freedom.\n");
    }
    cvm::real dt = cvm::dt();
-    cvm::real f_ext;
+    colvarvalue f_ext(fr.type());
    f_ext.reset();
    // the total force is applied to the fictitious mass, while the
-    // atoms only feel the harmonic force
+    // atoms only feel the harmonic force + wall force
    // fr: bias force on extended variable (without harmonic spring), for output in trajectory
    // f_ext: total force on extended variable (including harmonic spring)
-    // f: - initially, external biasing force (including wall forces)
+    // f: - initially, external biasing force
-    //    - after this code block, colvar force to be applied to atomic coordinates, ie. spring force
+    //    - after this code block, colvar force to be applied to atomic coordinates
    //      ie. spring force + wall force
    fr    = f;
    f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
    f     =     (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x);
@ -1200,15 +1248,24 @@ cvm::real colvar::update_forces_energy()
    potential_energy = 0.5 * ext_force_k * this->dist2(xr, x);
    // leap to v_(i+1/2)
    if (is_enabled(f_cv_Langevin)) {
-      vr -= dt * ext_gamma * vr.real_value;
+      vr -= dt * ext_gamma * vr;
-      vr += dt * ext_sigma * cvm::rand_gaussian() / ext_mass;
+      colvarvalue rnd(x);
      rnd.set_random();
      vr += dt * ext_sigma * rnd / ext_mass;
    }
    vr  += (0.5 * dt) * f_ext / ext_mass;
    xr  += dt * vr;
    xr.apply_constraints();
-    if (this->b_periodic) this->wrap(xr);
+    if (this->is_enabled(f_cv_periodic)) this->wrap(xr);
  }
  // TODO remove the wall force
  f += fw;
  // Now adding the force on the actual colvar (for those biases who
  // bypass the extended Lagrangian mass)
  f += fb_actual;
  // Store force to be applied, possibly summed over several timesteps
  f_accumulated += f;
  if (is_enabled(f_cv_fdiff_velocity)) {
@ -1425,14 +1482,15 @@ std::istream & colvar::read_restart(std::istream &is)
    }
  }
-  if ( !(get_keyval(conf, "x", x,
+  if ( !(get_keyval(conf, "x", x, x, colvarparse::parse_silent)) ) {
                    colvarvalue(x.type()), colvarparse::parse_silent)) ) {
    cvm::log("Error: restart file does not contain "
             "the value of the colvar \""+
             name+"\" .\n");
  } else {
    cvm::log("Restarting collective variable \""+name+"\" from value: "+
             cvm::to_str(x)+"\n");
    x_restart = x;
    after_restart = true;
  }
  if (is_enabled(f_cv_extended_Lagrangian)) {
--- a/lib/colvars/colvar.h
+++ b/lib/colvars/colvar.h
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #ifndef COLVAR_H
 #define COLVAR_H
@ -170,6 +177,9 @@ public:
  /// the biases are updated
  colvarvalue fb;
  /// \brief Bias force to the actual value (only useful with extended Lagrangian)
  colvarvalue fb_actual;
  /// \brief Total \em applied force; fr (if extended_lagrangian
  /// is defined), fb (if biases are applied) and the walls' forces
  /// (if defined) contribute to it
@ -183,13 +193,9 @@ public:
  colvarvalue ft;
-  /// Period, if it is a constant
+  /// Period, if this variable is periodic
  cvm::real period;
  /// \brief Same as above, but also takes into account components
  /// with a variable period, such as distanceZ
  bool b_periodic;
  /// \brief Expand the boundaries of multiples of width, to keep the
  /// value always within range
@ -290,6 +296,9 @@ public:
  /// Add to the total force from biases
  void add_bias_force(colvarvalue const &force);
  /// Apply a force to the actual value (only meaningful with extended Lagrangian)
  void add_bias_force_actual_value(colvarvalue const &force);
  /// \brief Collect all forces on this colvar, integrate internal
  /// equations of motion of internal degrees of freedom; see also
  /// colvar::communicate_forces()
@ -386,6 +395,12 @@ protected:
  /// Previous value (to calculate velocities during analysis)
  colvarvalue            x_old;
  /// Value read from the most recent state file (if any)
  colvarvalue            x_restart;
  /// True if a state file was just read
  bool                   after_restart;
  /// Time series of values and velocities used in correlation
  /// functions
  std::list< std::list<colvarvalue> > acf_x_history, acf_v_history;
@ -577,9 +592,20 @@ inline void colvar::add_bias_force(colvarvalue const &force)
 }
 inline void colvar::add_bias_force_actual_value(colvarvalue const &force)
 {
  if (cvm::debug()) {
    cvm::log("Adding biasing force "+cvm::to_str(force)+" to colvar \""+name+"\".\n");
  }
  fb_actual += force;
 }
 inline void colvar::reset_bias_force() {
  fb.type(value());
  fb.reset();
  fb_actual.type(value());
  fb_actual.reset();
 }
 #endif
--- a/lib/colvars/colvaratoms.cpp
+++ b/lib/colvars/colvaratoms.cpp
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include "colvarmodule.h"
 #include "colvarparse.h"
 #include "colvaratoms.h"
@ -171,7 +178,10 @@ int cvm::atom_group::remove_atom(cvm::atom_iter ai)
 int cvm::atom_group::init()
 {
-  if (!key.size()) key = "atoms";
+  if (!key.size()) key = "unnamed";
  description = "atom group " + key;
  // These will be overwritten by parse(), if initializing from a config string
  atoms.clear();
  // TODO: check with proxy whether atom forces etc are available
@ -179,6 +189,7 @@ int cvm::atom_group::init()
  index = -1;
  b_dummy = false;
  b_center = false;
  b_rotate = false;
  b_user_defined_fit = false;
@ -440,6 +451,7 @@ int cvm::atom_group::parse(std::string const &conf)
  if (b_print_atom_ids) {
    cvm::log("Internal definition of the atom group:\n");
    cvm::log(print_atom_ids());
  }
  cvm::decrease_depth();
--- a/lib/colvars/colvaratoms.h
+++ b/lib/colvars/colvaratoms.h
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #ifndef COLVARATOMS_H
 #define COLVARATOMS_H
@ -64,7 +71,7 @@ public:
  /// \brief Initialize an atom for collective variable calculation
  /// and get its internal identifier \param atom_number Atom index in
-  /// the system topology (starting from 1)
+  /// the system topology (1-based)
  atom(int atom_number);
  /// \brief Initialize an atom for collective variable calculation
@ -453,6 +460,8 @@ public:
  /// are not used, either because they were not defined (e.g because
  /// the colvar has not a scalar value) or the biases require to
  /// micromanage the force.
  /// This function will be phased out eventually, in favor of
  /// apply_colvar_force() once that is implemented for non-scalar values
  void apply_force(cvm::rvector const &force);
 };
--- a/lib/colvars/colvarbias.cpp
+++ b/lib/colvars/colvarbias.cpp
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include "colvarmodule.h"
 #include "colvarvalue.h"
 #include "colvarbias.h"
@ -29,6 +36,7 @@ colvarbias::colvarbias(char const *key)
  has_data = false;
  b_output_energy = false;
  reset();
  state_file_step = 0;
  // Start in active state by default
  enable(f_cvb_active);
@ -141,20 +149,25 @@ int colvarbias::clear()
 int colvarbias::add_colvar(std::string const &cv_name)
 {
  if (colvar *cv = cvm::colvar_by_name(cv_name)) {
-    // Removed this as nor all biases apply forces eg histogram
+
    // cv->enable(colvar::task_gradients);
    if (cvm::debug()) {
      cvm::log("Applying this bias to collective variable \""+
               cv->name+"\".\n");
    }
    colvars.push_back(cv);
    colvar_forces.push_back(colvarvalue());
    colvar_forces.back().type(cv->value()); // make sure each force is initialized to zero
    colvar_forces.back().is_derivative(); // colvar constraints are not applied to the force
    colvar_forces.back().reset();
    cv->biases.push_back(this); // add back-reference to this bias to colvar
    if (is_enabled(f_cvb_apply_force)) {
      cv->enable(f_cv_gradient);
    }
    // Add dependency link.
    // All biases need at least the value of each colvar
    // although possibly not at all timesteps
@ -165,6 +178,7 @@ int colvarbias::add_colvar(std::string const &cv_name)
               cv_name+"\".\n", INPUT_ERROR);
    return INPUT_ERROR;
  }
  return COLVARS_OK;
 }
@ -190,16 +204,19 @@ void colvarbias::communicate_forces()
 }
-void colvarbias::change_configuration(std::string const &conf)
+int colvarbias::change_configuration(std::string const &conf)
 {
-  cvm::error("Error: change_configuration() not implemented.\n");
+  cvm::error("Error: change_configuration() not implemented.\n",
             COLVARS_NOT_IMPLEMENTED);
  return COLVARS_NOT_IMPLEMENTED;
 }
 cvm::real colvarbias::energy_difference(std::string const &conf)
 {
-  cvm::error("Error: energy_difference() not implemented.\n");
+  cvm::error("Error: energy_difference() not implemented.\n",
-  return 0.;
+             COLVARS_NOT_IMPLEMENTED);
  return 0.0;
 }
@ -225,6 +242,118 @@ int colvarbias::replica_share()
  return COLVARS_NOT_IMPLEMENTED;
 }
 std::string const colvarbias::get_state_params() const
 {
  std::ostringstream os;
  os << "step " << cvm::step_absolute() << "\n"
     << "name " << this->name << "\n";
  return os.str();
 }
 int colvarbias::set_state_params(std::string const &conf)
 {
  std::string new_name = "";
  if (colvarparse::get_keyval(conf, "name", new_name,
                              std::string(""), colvarparse::parse_silent) &&
      (new_name != this->name)) {
    cvm::error("Error: in the state file, the "
               "\""+bias_type+"\" block has a different name, \""+new_name+
               "\": different system?\n", INPUT_ERROR);
  }
  if (name.size() == 0) {
    cvm::error("Error: \""+bias_type+"\" block within the restart file "
               "has no identifiers.\n", INPUT_ERROR);
  }
  colvarparse::get_keyval(conf, "step", state_file_step,
                          cvm::step_absolute(), colvarparse::parse_silent);
  return COLVARS_OK;
 }
 std::ostream & colvarbias::write_state(std::ostream &os)
 {
  if (cvm::debug()) {
    cvm::log("Writing state file for bias \""+name+"\"\n");
  }
  os.setf(std::ios::scientific, std::ios::floatfield);
  os.precision(cvm::cv_prec);
  os << bias_type << " {\n"
     << "  configuration {\n";
  std::istringstream is(get_state_params());
  std::string line;
  while (std::getline(is, line)) {
    os << "    " << line << "\n";
  }
  os << "  }\n";
  write_state_data(os);
  os << "}\n\n";
  return os;
 }
 std::istream & colvarbias::read_state(std::istream &is)
 {
  size_t const start_pos = is.tellg();
  std::string key, brace, conf;
  if ( !(is >> key)   || !(key == bias_type) ||
       !(is >> brace) || !(brace == "{") ||
       !(is >> colvarparse::read_block("configuration", conf)) ||
       (set_state_params(conf) != COLVARS_OK) ) {
    cvm::error("Error: in reading state configuration for \""+bias_type+"\" bias \""+
               this->name+"\" at position "+
               cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR);
    is.clear();
    is.seekg(start_pos, std::ios::beg);
    is.setstate(std::ios::failbit);
    return is;
  }
  if (!read_state_data(is)) {
    cvm::error("Error: in reading state data for \""+bias_type+"\" bias \""+
               this->name+"\" at position "+
               cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR);
    is.clear();
    is.seekg(start_pos, std::ios::beg);
    is.setstate(std::ios::failbit);
  }
  is >> brace;
  if (brace != "}") {
    cvm::error("Error: corrupt restart information for \""+bias_type+"\" bias \""+
               this->name+"\": no matching brace at position "+
               cvm::to_str(is.tellg())+" in stream.\n");
    is.setstate(std::ios::failbit);
  }
  return is;
 }
 std::istream & colvarbias::read_state_data_key(std::istream &is, char const *key)
 {
  size_t const start_pos = is.tellg();
  std::string key_in;
  if ( !(is >> key_in) ||
       !(key_in == to_lower_cppstr(std::string(key))) ) {
    cvm::error("Error: in reading restart configuration for "+
               bias_type+" bias \""+this->name+"\" at position "+
               cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR);
    is.clear();
    is.seekg(start_pos, std::ios::beg);
    is.setstate(std::ios::failbit);
    return is;
  }
  return is;
 }
 std::ostream & colvarbias::write_traj_label(std::ostream &os)
 {
  os << " ";
--- a/lib/colvars/colvarbias.h
+++ b/lib/colvars/colvarbias.h
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #ifndef COLVARBIAS_H
 #define COLVARBIAS_H
@ -9,7 +16,8 @@
 /// \brief Collective variable bias, base class
-class colvarbias : public colvarparse, public colvardeps {
+class colvarbias
  : public virtual colvarparse, public virtual colvardeps {
 public:
  /// Name of this bias
@ -24,6 +32,12 @@ public:
  /// Add a new collective variable to this bias
  int add_colvar(std::string const &cv_name);
  /// Add a new collective variable to this bias
  size_t number_of_colvars() const
  {
    return colvars.size();
  }
  /// Retrieve colvar values and calculate their biasing forces
  /// Return bias energy
  virtual int update();
@ -31,7 +45,7 @@ public:
  // TODO: move update_bias here (share with metadynamics)
  /// Load new configuration - force constant and/or centers only
-  virtual void change_configuration(std::string const &conf);
+  virtual int change_configuration(std::string const &conf);
  /// Calculate change in energy from using alternate configuration
  virtual cvm::real energy_difference(std::string const &conf);
@ -49,7 +63,7 @@ public:
  virtual void analyze() {}
  /// Send forces to the collective variables
-  void communicate_forces();
+  virtual void communicate_forces();
  /// \brief Constructor
  colvarbias(char const *key);
@ -60,13 +74,11 @@ public:
  /// \brief Set to zero all mutable data
  virtual int reset();
-protected:
+private:
  /// Default constructor
  colvarbias();
 private:
  /// Copy constructor
  colvarbias(colvarbias &);
@ -78,28 +90,59 @@ public:
  /// Destructor
  virtual ~colvarbias();
-  /// Read the bias configuration from a restart file
+  /// Write the values of specific mutable properties to a string
-  virtual std::istream & read_restart(std::istream &is) = 0;
+  virtual std::string const get_state_params() const;
-  /// Write the bias configuration to a restart file
+  /// Read the values of specific mutable properties from a string
-  virtual std::ostream & write_restart(std::ostream &os) = 0;
+  virtual int set_state_params(std::string const &state_conf);
  /// Write all mutable data not already written by get_state_params()
  virtual std::ostream & write_state_data(std::ostream &os)
  {
    return os;
  }
  /// Read all mutable data not already set by set_state_params()
  virtual std::istream & read_state_data(std::istream &is)
  {
    return is;
  }
  /// Read a keyword from the state data (typically a header)
  std::istream & read_state_data_key(std::istream &is, char const *key);
  /// Write the bias configuration to a restart file or other stream
  virtual std::ostream & write_state(std::ostream &os);
  /// Read the bias configuration from a restart file or other stream
  virtual std::istream & read_state(std::istream &is);
  /// Write a label to the trajectory file (comment line)
  virtual std::ostream & write_traj_label(std::ostream &os);
  /// (Re)initialize the output files (does not write them yet)
  virtual int setup_output() { return COLVARS_OK; }
  /// Output quantities such as the bias energy to the trajectory file
  virtual std::ostream & write_traj(std::ostream &os);
-  /// Write output files (if defined, e.g. in analysis mode)
+  /// (Re)initialize the output files (does not write them yet)
  virtual int setup_output()
  {
    return COLVARS_OK;
  }
  /// Write any output files that this bias may have (e.g. PMF files)
  virtual int write_output_files()
  {
    return COLVARS_OK;
  }
-  inline cvm::real get_energy() {
+  /// If this bias is communicating with other replicas through files, send it to them
  virtual int write_state_to_replicas()
  {
    return COLVARS_OK;
  }
  inline cvm::real get_energy()
  {
    return bias_energy;
  }
@ -107,9 +150,11 @@ public:
  static std::vector<feature *> cvb_features;
  /// \brief Implementation of the feature list accessor for colvarbias
-  virtual std::vector<feature *> &features() {
+  virtual std::vector<feature *> &features()
  {
    return cvb_features;
  }
 protected:
  /// \brief Pointers to collective variables to which the bias is
@ -130,6 +175,9 @@ protected:
  /// (for history-dependent biases)
  bool                     has_data;
  /// \brief Step number read from the last state file
  size_t                   state_file_step;
 };
 #endif
--- a/lib/colvars/colvarbias_abf.cpp
+++ b/lib/colvars/colvarbias_abf.cpp
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include "colvarmodule.h"
 #include "colvar.h"
 #include "colvarbias_abf.h"
@ -23,6 +30,8 @@ int colvarbias_abf::init(std::string const &conf)
 {
  colvarbias::init(conf);
  provide(f_cvb_scalar_variables);
  enable(f_cvb_scalar_variables);
  provide(f_cvb_history_dependent);
  // TODO relax this in case of VMD plugin
@ -590,16 +599,10 @@ void colvarbias_abf::read_gradients_samples()
 }
-std::ostream & colvarbias_abf::write_restart(std::ostream& os)
+std::ostream & colvarbias_abf::write_state_data(std::ostream& os)
 {
  std::ios::fmtflags flags(os.flags());
  os << "abf {\n"
     << "  configuration {\n"
     << "    name " << this->name << "\n";
  os << "  }\n";
  os.setf(std::ios::fmtflags(0), std::ios::floatfield); // default floating-point format
  os << "\nsamples\n";
  samples->write_raw(os, 8);
@ -617,117 +620,47 @@ std::ostream & colvarbias_abf::write_restart(std::ostream& os)
    z_gradients->write_raw(os, 8);
  }
  os << "}\n\n";
  os.flags(flags);
  return os;
 }
-std::istream & colvarbias_abf::read_restart(std::istream& is)
+std::istream & colvarbias_abf::read_state_data(std::istream& is)
 {
  if ( input_prefix.size() > 0 ) {
    cvm::error("ERROR: cannot provide both inputPrefix and a colvars state file.\n", INPUT_ERROR);
  }
-  size_t const start_pos = is.tellg();
+  if (! read_state_data_key(is, "samples")) {
  cvm::log("Restarting ABF bias \""+
            this->name+"\".\n");
  std::string key, brace, conf;
  if ( !(is >> key)   || !(key == "abf") ||
       !(is >> brace) || !(brace == "{") ||
       !(is >> colvarparse::read_block("configuration", conf)) ) {
    cvm::log("Error: in reading restart configuration for ABF bias \""+
              this->name+"\" at position "+
              cvm::to_str(is.tellg())+" in stream.\n");
    is.clear();
    is.seekg(start_pos, std::ios::beg);
    is.setstate(std::ios::failbit);
    return is;
  }
  std::string name = "";
  if ( (colvarparse::get_keyval(conf, "name", name, std::string(""), colvarparse::parse_silent)) &&
         (name != this->name) )
    cvm::error("Error: in the restart file, the "
                      "\"abf\" block has wrong name(" + name + ")\n");
  if ( name == "" ) {
    cvm::error("Error: \"abf\" block in the restart file has no name.\n");
  }
  if ( !(is >> key)   || !(key == "samples")) {
    cvm::log("Error: in reading restart configuration for ABF bias \""+
              this->name+"\" at position "+
              cvm::to_str(is.tellg())+" in stream.\n");
    is.clear();
    is.seekg(start_pos, std::ios::beg);
    is.setstate(std::ios::failbit);
    return is;
  }
  if (! samples->read_raw(is)) {
    is.clear();
    is.seekg(start_pos, std::ios::beg);
    is.setstate(std::ios::failbit);
    return is;
  }
-  if ( !(is >> key)   || !(key == "gradient")) {
+  if (! read_state_data_key(is, "gradient")) {
    cvm::log("Error: in reading restart configuration for ABF bias \""+
              this->name+"\" at position "+
              cvm::to_str(is.tellg())+" in stream.\n");
    is.clear();
    is.seekg(start_pos, std::ios::beg);
    is.setstate(std::ios::failbit);
    return is;
  }
  if (! gradients->read_raw(is)) {
    is.clear();
    is.seekg(start_pos, std::ios::beg);
    is.setstate(std::ios::failbit);
    return is;
  }
  if (z_gradients) {
-    if ( !(is >> key) || !(key == "z_samples")) {
+
-      cvm::log("Error: in reading restart configuration for ABF bias \""+
+    if (! read_state_data_key(is, "z_samples")) {
                this->name+"\" at position "+
                cvm::to_str(is.tellg())+" in stream.\n");
      is.clear();
      is.seekg(start_pos, std::ios::beg);
      is.setstate(std::ios::failbit);
      return is;
    }
    if (! z_samples->read_raw(is)) {
      is.clear();
      is.seekg(start_pos, std::ios::beg);
      is.setstate(std::ios::failbit);
      return is;
    }
-    if ( !(is >> key)   || !(key == "z_gradient")) {
+    if (! read_state_data_key(is, "z_gradient")) {
      cvm::log("Error: in reading restart configuration for ABF bias \""+
                this->name+"\" at position "+
                cvm::to_str(is.tellg())+" in stream.\n");
      is.clear();
      is.seekg(start_pos, std::ios::beg);
      is.setstate(std::ios::failbit);
      return is;
    }
    if (! z_gradients->read_raw(is)) {
      is.clear();
      is.seekg(start_pos, std::ios::beg);
      is.setstate(std::ios::failbit);
      return is;
    }
  }
-  is >> brace;
+
  if (brace != "}") {
    cvm::error("Error: corrupt restart information for ABF bias \""+
                      this->name+"\": no matching brace at position "+
                      cvm::to_str(is.tellg())+" in the restart file.\n");
    is.setstate(std::ios::failbit);
  }
  return is;
 }
--- a/lib/colvars/colvarbias_abf.h
+++ b/lib/colvars/colvarbias_abf.h
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #ifndef COLVARBIAS_ABF_H
 #define COLVARBIAS_ABF_H
@ -107,8 +114,8 @@ private:
  /// Read human-readable FE gradients and sample count (if not using restart)
  void		  read_gradients_samples();
-  std::istream& read_restart(std::istream&);
+  std::istream& read_state_data(std::istream&);
-  std::ostream& write_restart(std::ostream&);
+  std::ostream& write_state_data(std::ostream&);
 };
 #endif
--- a/lib/colvars/colvarbias_alb.cpp
+++ b/lib/colvars/colvarbias_alb.cpp
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
@ -33,6 +40,9 @@ int colvarbias_alb::init(std::string const &conf)
 {
  colvarbias::init(conf);
  provide(f_cvb_scalar_variables);
  enable(f_cvb_scalar_variables);
  provide(f_cvb_history_dependent);
  size_t i;
@ -239,37 +249,8 @@ int colvarbias_alb::update()
 }
-std::istream & colvarbias_alb::read_restart(std::istream &is)
+int colvarbias_alb::set_state_params(std::string const &conf)
 {
  size_t const start_pos = is.tellg();
  cvm::log("Restarting adaptive linear bias \""+
            this->name+"\".\n");
  std::string key, brace, conf;
  if ( !(is >> key)   || !(key == "ALB") ||
       !(is >> brace) || !(brace == "{") ||
       !(is >> colvarparse::read_block("configuration", conf)) ) {
    cvm::log("Error: in reading restart configuration for restraint bias \""+
              this->name+"\" at position "+
              cvm::to_str(is.tellg())+" in stream.\n");
    is.clear();
    is.seekg(start_pos, std::ios::beg);
    is.setstate(std::ios::failbit);
    return is;
  }
  std::string name = "";
  if ( (colvarparse::get_keyval(conf, "name", name, std::string(""), colvarparse::parse_silent)) &&
       (name != this->name) )
    cvm::fatal_error("Error: in the restart file, the "
                      "\"ALB\" block has a wrong name\n");
  if (name.size() == 0) {
    cvm::fatal_error("Error: \"ALB\" block in the restart file "
                      "has no identifiers.\n");
  }
  if (!get_keyval(conf, "setCoupling", set_coupling))
    cvm::fatal_error("Error: current setCoupling  is missing from the restart.\n");
@ -299,23 +280,13 @@ std::istream & colvarbias_alb::read_restart(std::istream &is)
  if (!get_keyval(conf, "b_equilibration", b_equilibration))
    cvm::fatal_error("Error: current updateCalls is missing from the restart.\n");
-  is >> brace;
+  return COLVARS_OK;
  if (brace != "}") {
    cvm::fatal_error("Error: corrupt restart information for adaptive linear bias \""+
                      this->name+"\": no matching brace at position "+
                      cvm::to_str(is.tellg())+" in the restart file.\n");
    is.setstate(std::ios::failbit);
  }
  return is;
 }
-std::ostream & colvarbias_alb::write_restart(std::ostream &os)
+std::string const colvarbias_alb::get_state_params() const
 {
-  os << "ALB {\n"
+  std::ostringstream os;
     << "  configuration {\n"
     << "    name " << this->name << "\n";
  os << "    setCoupling ";
  size_t i;
  for (i = 0; i < colvars.size(); i++) {
@ -358,10 +329,7 @@ std::ostream & colvarbias_alb::write_restart(std::ostream &os)
  else
    os << "    b_equilibration no\n";
-  os << "  }\n"
+  return os.str();
     << "}\n\n";
  return os;
 }
--- a/lib/colvars/colvarbias_alb.h
+++ b/lib/colvars/colvarbias_alb.h
@ -1,10 +1,18 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #ifndef COLVARBIAS_ALB_H
 #define COLVARBIAS_ALB_H
 #include "colvar.h"
-#include "colvarbias_restraint.h"
+#include "colvarbias.h"
 class colvarbias_alb : public colvarbias {
@ -15,8 +23,8 @@ public:
  virtual int init(std::string const &conf);
  virtual int update();
-  virtual std::istream & read_restart(std::istream &is);
+  virtual std::string const get_state_params() const;
-  virtual std::ostream & write_restart(std::ostream &os);
+  virtual int set_state_params(std::string const &conf);
  virtual std::ostream & write_traj_label(std::ostream &os);
  virtual std::ostream & write_traj(std::ostream &os);
--- a/lib/colvars/colvarbias_histogram.cpp
+++ b/lib/colvars/colvarbias_histogram.cpp
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include "colvarmodule.h"
 #include "colvar.h"
 #include "colvarbias_histogram.h"
@ -17,6 +24,9 @@ int colvarbias_histogram::init(std::string const &conf)
 {
  colvarbias::init(conf);
  provide(f_cvb_scalar_variables);
  enable(f_cvb_scalar_variables);
  provide(f_cvb_history_dependent);
  enable(f_cvb_history_dependent);
@ -196,78 +206,25 @@ int colvarbias_histogram::write_output_files()
 }
-std::istream & colvarbias_histogram::read_restart(std::istream& is)
+std::istream & colvarbias_histogram::read_state_data(std::istream& is)
 {
-  size_t const start_pos = is.tellg();
+  if (! read_state_data_key(is, "grid")) {
  cvm::log("Restarting collective variable histogram \""+
            this->name+"\".\n");
  std::string key, brace, conf;
  if ( !(is >> key)   || !(key == "histogram") ||
       !(is >> brace) || !(brace == "{") ||
       !(is >> colvarparse::read_block("configuration", conf)) ) {
    cvm::log("Error: in reading restart configuration for histogram \""+
              this->name+"\" at position "+
              cvm::to_str(is.tellg())+" in stream.\n");
    is.clear();
    is.seekg(start_pos, std::ios::beg);
    is.setstate(std::ios::failbit);
    return is;
  }
  int id = -1;
  std::string name = "";
  if ( (colvarparse::get_keyval(conf, "name", name, std::string(""), colvarparse::parse_silent)) &&
         (name != this->name) )
    cvm::error("Error: in the restart file, the "
                      "\"histogram\" block has a wrong name: different system?\n");
  if ( (id == -1) && (name == "") ) {
    cvm::error("Error: \"histogram\" block in the restart file "
                      "has no name.\n");
  }
  if ( !(is >> key)   || !(key == "grid")) {
    cvm::error("Error: in reading restart configuration for histogram \""+
              this->name+"\" at position "+
              cvm::to_str(is.tellg())+" in stream.\n");
    is.clear();
    is.seekg(start_pos, std::ios::beg);
    is.setstate(std::ios::failbit);
    return is;
  }
  if (! grid->read_raw(is)) {
    is.clear();
    is.seekg(start_pos, std::ios::beg);
    is.setstate(std::ios::failbit);
    return is;
  }
  is >> brace;
  if (brace != "}") {
    cvm::error("Error: corrupt restart information for ABF bias \""+
                this->name+"\": no matching brace at position "+
                cvm::to_str(is.tellg())+" in the restart file.\n");
    is.setstate(std::ios::failbit);
  }
  return is;
 }
-std::ostream & colvarbias_histogram::write_restart(std::ostream& os)
+
 std::ostream & colvarbias_histogram::write_state_data(std::ostream& os)
 {
  std::ios::fmtflags flags(os.flags());
  os.setf(std::ios::fmtflags(0), std::ios::floatfield);
  os << "histogram {\n"
     << "  configuration {\n"
     << "    name " << this->name << "\n";
  os << "  }\n";
  os << "grid\n";
  grid->write_raw(os, 8);
  os << "}\n\n";
  os.flags(flags);
  return os;
 }
--- a/lib/colvars/colvarbias_histogram.h
+++ b/lib/colvars/colvarbias_histogram.h
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #ifndef COLVARBIAS_HISTOGRAM_H
 #define COLVARBIAS_HISTOGRAM_H
@ -35,8 +42,8 @@ protected:
  /// If colvar_array_size is larger than 1, weigh each one by this number before accumulating the histogram
  std::vector<cvm::real> weights;
-  virtual std::istream& read_restart(std::istream&);
+  virtual std::istream & read_state_data(std::istream &is);
-  virtual std::ostream& write_restart(std::ostream&);
+  virtual std::ostream & write_state_data(std::ostream &os);
 };
 #endif
--- a/lib/colvars/colvarbias_meta.cpp
+++ b/lib/colvars/colvarbias_meta.cpp
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include <iostream>
 #include <sstream>
 #include <fstream>
@ -26,10 +33,9 @@
 colvarbias_meta::colvarbias_meta(char const *key)
-  : colvarbias(key),
+  : colvarbias(key)
    new_hills_begin(hills.end()),
    state_file_step(0)
 {
  new_hills_begin = hills.end();
 }
@ -164,11 +170,32 @@ int colvarbias_meta::init(std::string const &conf)
  target_dist = NULL;
  get_keyval(conf, "ebMeta", ebmeta, false);
  if(ebmeta){
    if (use_grids && expand_grids) {
      cvm::fatal_error("Error: expandBoundaries is not supported with "
                       "ebMeta please allocate wide enough boundaries for "
                       "each colvar ahead of time and set targetdistfile "
                       "accordingly. \n");
    }
    target_dist = new colvar_grid_scalar();
    target_dist->init_from_colvars(colvars);
    get_keyval(conf, "targetdistfile", target_dist_file);
    std::ifstream targetdiststream(target_dist_file.c_str());
    target_dist->read_multicol(targetdiststream);
    cvm::real min_val = target_dist->minimum_value();
    if(min_val<0){
      cvm::error("Error: Target distribution of ebMeta "
                 "has negative values!.\n", INPUT_ERROR);
    }
    cvm::real min_pos_val = target_dist->minimum_pos_value();
    if(min_pos_val<=0){
      cvm::error("Error: Target distribution of ebMeta has negative "
                 "or zero minimum positive value!.\n", INPUT_ERROR);
    }
    if(min_val==0){
      cvm::log("WARNING: Target distribution has zero values.\n");
      cvm::log("Zeros will be converted to the minimum positive value.\n");
      target_dist->remove_zeros(min_pos_val);
    }
    // normalize target distribution and multiply by effective volume = exp(differential entropy)
    target_dist->multiply_constant(1.0/target_dist->integral());
    cvm::real volume = std::exp(target_dist->entropy());
@ -404,8 +431,10 @@ int colvarbias_meta::update()
    if (ebmeta) {
       hills_scale *= 1.0/target_dist->value(target_dist->get_colvars_index());
-       if(cvm::step_absolute() <= ebmeta_equil_steps) {
+       if(cvm::step_absolute() <= long(ebmeta_equil_steps)) {
-         cvm::real const hills_lambda=(cvm::real(ebmeta_equil_steps - cvm::step_absolute()))/(cvm::real(ebmeta_equil_steps));
+         cvm::real const hills_lambda =
           (cvm::real(long(ebmeta_equil_steps) - cvm::step_absolute())) /
           (cvm::real(ebmeta_equil_steps));
         hills_scale = hills_lambda + (1-hills_lambda)*hills_scale;
       }
    }
@ -972,7 +1001,7 @@ void colvarbias_meta::read_replica_files()
               (replicas[ir])->replica_state_file+"\".\n");
      std::ifstream is((replicas[ir])->replica_state_file.c_str());
-      if (! (replicas[ir])->read_restart(is)) {
+      if (! (replicas[ir])->read_state(is)) {
        cvm::log("Reading from file \""+(replicas[ir])->replica_state_file+
                 "\" failed or incomplete: will try again in "+
                 cvm::to_str(replica_update_freq)+" steps.\n");
@ -1068,63 +1097,24 @@ void colvarbias_meta::read_replica_files()
 }
-// **********************************************************************
+int colvarbias_meta::set_state_params(std::string const &state_conf)
 // input functions
 // **********************************************************************
 std::istream & colvarbias_meta::read_restart(std::istream& is)
 {
-  size_t const start_pos = is.tellg();
+  std::string new_replica = "";
-
+  if (colvarparse::get_keyval(state_conf, "replicaID", new_replica,
-  if (comm == single_replica) {
+                              std::string(""), colvarparse::parse_silent) &&
-    // if using a multiple replicas scheme, output messages
+      (new_replica != this->replica_id)) {
-    // are printed before and after calling this function
+    cvm::error("Error: in the state file, the "
-    cvm::log("Restarting metadynamics bias \""+this->name+"\""+
+               "\"metadynamics\" block has a different replicaID: different system?\n",
-             ".\n");
+               INPUT_ERROR);
-  }
+    return INPUT_ERROR;
  std::string key, brace, conf;
  if ( !(is >> key)   || !(key == "metadynamics") ||
       !(is >> brace) || !(brace == "{") ||
       !(is >> colvarparse::read_block("configuration", conf)) ) {
    if (comm == single_replica)
      cvm::log("Error: in reading restart configuration for metadynamics bias \""+
               this->name+"\""+
               ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
               (replica_state_file_in_sync ? ("at position "+
                                              cvm::to_str(start_pos)+
                                              " in the state file") : "")+".\n");
    is.clear();
    is.seekg(start_pos, std::ios::beg);
    is.setstate(std::ios::failbit);
    return is;
  }
-  std::string name = "";
+  return COLVARS_OK;
-  if ( colvarparse::get_keyval(conf, "name", name,
+}
                               std::string(""), colvarparse::parse_silent) &&
       (name != this->name) )
    cvm::fatal_error("Error: in the restart file, the "
                     "\"metadynamics\" block has a different name: different system?\n");
  if (name.size() == 0) {
    cvm::fatal_error("Error: \"metadynamics\" block within the restart file "
                     "has no identifiers.\n");
  }
  if (comm != single_replica) {
    std::string replica = "";
    if (colvarparse::get_keyval(conf, "replicaID", replica,
                                std::string(""), colvarparse::parse_silent) &&
        (replica != this->replica_id))
      cvm::fatal_error("Error: in the restart file, the "
                       "\"metadynamics\" block has a different replicaID: different system?\n");
    colvarparse::get_keyval(conf, "step", state_file_step,
                            cvm::step_absolute(), colvarparse::parse_silent);
  }
 std::istream & colvarbias_meta::read_state_data(std::istream& is)
 {
  bool grids_from_restart_file = use_grids;
  if (use_grids) {
@ -1155,6 +1145,7 @@ std::istream & colvarbias_meta::read_restart(std::istream& is)
    }
    size_t const hills_energy_pos = is.tellg();
    std::string key;
    if (!(is >> key)) {
      if (hills_energy_backup != NULL) {
        delete hills_energy;
@ -1316,17 +1307,6 @@ std::istream & colvarbias_meta::read_restart(std::istream& is)
    }
  }
  is >> brace;
  if (brace != "}") {
    cvm::log("Incomplete restart information for metadynamics bias \""+
             this->name+"\""+
             ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
             ": no closing brace at position "+
             cvm::to_str(is.tellg())+" in the file.\n");
    is.setstate(std::ios::failbit);
    return is;
  }
  if (cvm::debug())
    cvm::log("colvarbias_meta::read_restart() done\n");
@ -1424,13 +1404,6 @@ std::istream & colvarbias_meta::read_hill(std::istream &is)
 }
 // **********************************************************************
 // output functions
 // **********************************************************************
 int colvarbias_meta::setup_output()
 {
@ -1537,16 +1510,17 @@ int colvarbias_meta::setup_output()
 }
-std::ostream & colvarbias_meta::write_restart(std::ostream& os)
+std::string const colvarbias_meta::get_state_params() const
 {
-  os << "metadynamics {\n"
+  std::ostringstream os;
     << "  configuration {\n"
     << "    step " << cvm::step_absolute() << "\n"
     << "    name " << this->name << "\n";
  if (this->comm != single_replica)
-    os << "    replicaID " << this->replica_id << "\n";
+    os << "replicaID " << this->replica_id << "\n";
-  os << "  }\n\n";
+  return (colvarbias::get_state_params() + os.str());
 }
 std::ostream & colvarbias_meta::write_state_data(std::ostream& os)
 {
  if (use_grids) {
    // this is a very good time to project hills, if you haven't done
@ -1578,8 +1552,12 @@ std::ostream & colvarbias_meta::write_restart(std::ostream& os)
    }
  }
-  os << "}\n\n";
+  return os;
 }
 int colvarbias_meta::write_state_to_replicas()
 {
  if (comm != single_replica) {
    write_replica_state_file();
    // schedule to reread the state files of the other replicas (they
@ -1588,12 +1566,16 @@ std::ostream & colvarbias_meta::write_restart(std::ostream& os)
      (replicas[ir])->replica_state_file_in_sync = false;
    }
  }
  return COLVARS_OK;
 }
 int colvarbias_meta::write_output_files()
 {
  if (dump_fes) {
    write_pmf();
  }
-
+  return COLVARS_OK;
  return os;
 }
@ -1665,53 +1647,67 @@ void colvarbias_meta::write_pmf()
-void colvarbias_meta::write_replica_state_file()
+int colvarbias_meta::write_replica_state_file()
 {
-  // write down also the restart for the other replicas: TODO: this
+  if (cvm::debug()) {
-  // is duplicated code, that could be cleaned up later
+    cvm::log("Writing replica state file for bias \""+name+"\"\n");
  }
  // write down also the restart for the other replicas
  cvm::backup_file(replica_state_file.c_str());
-  cvm::ofstream rep_state_os(replica_state_file.c_str());
+  std::ostream *rep_state_os = cvm::proxy->output_stream(replica_state_file);
-  if (!rep_state_os.is_open())
+  if (rep_state_os == NULL) {
-    cvm::fatal_error("Error: in opening file \""+
+    cvm::error("Error: in opening file \""+
-                     replica_state_file+"\" for writing.\n");
+               replica_state_file+"\" for writing.\n", FILE_ERROR);
-
+    return FILE_ERROR;
  rep_state_os.setf(std::ios::scientific, std::ios::floatfield);
  rep_state_os << "\n"
               << "metadynamics {\n"
               << "  configuration {\n"
               << "    name " << this->name << "\n"
               << "    step " << cvm::step_absolute() << "\n";
  if (this->comm != single_replica) {
    rep_state_os << "    replicaID " << this->replica_id << "\n";
  }
  rep_state_os << "  }\n\n";
  rep_state_os << "  hills_energy\n";
  rep_state_os << std::setprecision(cvm::cv_prec)
               << std::setw(cvm::cv_width);
  hills_energy->write_restart(rep_state_os);
  rep_state_os << "  hills_energy_gradients\n";
  rep_state_os << std::setprecision(cvm::cv_prec)
               << std::setw(cvm::cv_width);
  hills_energy_gradients->write_restart(rep_state_os);
  if ( (!use_grids) || keep_hills ) {
    // write all hills currently in memory
    for (std::list<hill>::const_iterator h = this->hills.begin();
         h != this->hills.end();
         h++) {
      rep_state_os << *h;
    }
  } else {
    // write just those that are near the grid boundaries
    for (std::list<hill>::const_iterator h = this->hills_off_grid.begin();
         h != this->hills_off_grid.end();
         h++) {
      rep_state_os << *h;
    }
  }
-  rep_state_os << "}\n\n";
+  rep_state_os->setf(std::ios::scientific, std::ios::floatfield);
-  rep_state_os.close();
+
  if (!write_state(*rep_state_os)) {
    cvm::error("Error: in writing to file \""+
               replica_state_file+"\".\n", FILE_ERROR);
    cvm::proxy->close_output_stream(replica_state_file);
    return FILE_ERROR;
  }
  cvm::proxy->close_output_stream(replica_state_file);
  // rep_state_os.setf(std::ios::scientific, std::ios::floatfield);
  // rep_state_os << "\n"
  //              << "metadynamics {\n"
  //              << "  configuration {\n"
  //              << "    name " << this->name << "\n"
  //              << "    step " << cvm::step_absolute() << "\n";
  // if (this->comm != single_replica) {
  //   rep_state_os << "    replicaID " << this->replica_id << "\n";
  // }
  // rep_state_os << "  }\n\n";
  // rep_state_os << "  hills_energy\n";
  // rep_state_os << std::setprecision(cvm::cv_prec)
  //              << std::setw(cvm::cv_width);
  // hills_energy->write_restart(rep_state_os);
  // rep_state_os << "  hills_energy_gradients\n";
  // rep_state_os << std::setprecision(cvm::cv_prec)
  //              << std::setw(cvm::cv_width);
  // hills_energy_gradients->write_restart(rep_state_os);
  // if ( (!use_grids) || keep_hills ) {
  //   // write all hills currently in memory
  //   for (std::list<hill>::const_iterator h = this->hills.begin();
  //        h != this->hills.end();
  //        h++) {
  //     rep_state_os << *h;
  //   }
  // } else {
  //   // write just those that are near the grid boundaries
  //   for (std::list<hill>::const_iterator h = this->hills_off_grid.begin();
  //        h != this->hills_off_grid.end();
  //        h++) {
  //     rep_state_os << *h;
  //   }
  // }
  // rep_state_os << "}\n\n";
  // rep_state_os.close();
  // reopen the hills file
  replica_hills_os.close();
@ -1721,8 +1717,11 @@ void colvarbias_meta::write_replica_state_file()
    cvm::fatal_error("Error: in opening file \""+
                     replica_hills_file+"\" for writing.\n");
  replica_hills_os.setf(std::ios::scientific, std::ios::floatfield);
  return COLVARS_OK;
 }
 std::string colvarbias_meta::hill::output_traj()
 {
  std::ostringstream os;
--- a/lib/colvars/colvarbias_meta.h
+++ b/lib/colvars/colvarbias_meta.h
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #ifndef COLVARBIAS_META_H
 #define COLVARBIAS_META_H
@ -31,10 +38,16 @@ public:
  virtual int init(std::string const &conf);
  virtual ~colvarbias_meta();
  virtual int update();
-  virtual std::istream & read_restart(std::istream &is);
+
-  virtual std::ostream & write_restart(std::ostream &os);
+  virtual std::string const get_state_params() const;
  virtual int set_state_params(std::string const &state_conf);
  virtual std::ostream & write_state_data(std::ostream &os);
  virtual std::istream & read_state_data(std::istream &os);
  virtual int setup_output();
  virtual int write_output_files();
  virtual void write_pmf();
  virtual int write_state_to_replicas();
  class hill;
  typedef std::list<hill>::iterator hill_iter;
@ -77,13 +90,6 @@ protected:
  /// Read a hill from a file
  std::istream & read_hill(std::istream &is);
  /// \brief step present in a state file
  ///
  /// When using grids and reading state files containing them
  /// (multiple replicas), this is used to check whether a hill is
  /// newer or older than the grids
  size_t                   state_file_step;
  /// \brief Add a new hill; if a .hills trajectory is written,
  /// write it there; if there is more than one replica, communicate
  /// it to the others
@ -187,7 +193,7 @@ protected:
  virtual void read_replica_files();
  /// \brief Write data to other replicas
-  virtual void write_replica_state_file();
+  virtual int write_replica_state_file();
  /// \brief Additional, "mirror" metadynamics biases, to collect info
  /// from the other replicas
--- a/lib/colvars/colvarbias_restraint.cpp
+++ b/lib/colvars/colvarbias_restraint.cpp
--- a/lib/colvars/colvarbias_restraint.h
+++ b/lib/colvars/colvarbias_restraint.h
@ -1,36 +1,43 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #ifndef COLVARBIAS_RESTRAINT_H
 #define COLVARBIAS_RESTRAINT_H
 #include "colvarbias.h"
-/// \brief Bias restraint, optionally moving towards a target
+/// \brief Most general definition of a colvar restraint:
 /// see derived classes for specific types
 /// (implementation of \link colvarbias \endlink)
-class colvarbias_restraint : public colvarbias {
+class colvarbias_restraint
  : public virtual colvarbias
 {
 public:
  /// Retrieve colvar values and calculate their biasing forces
  virtual int update();
  // TODO the following can be supplanted by a new call to init()
  /// Load new configuration - force constant and/or centers only
-  virtual void change_configuration(std::string const &conf);
+  virtual int change_configuration(std::string const &conf) { return COLVARS_NOT_IMPLEMENTED; }
  /// Calculate change in energy from using alternate configuration
-  virtual cvm::real energy_difference(std::string const &conf);
+  virtual cvm::real energy_difference(std::string const &conf) { return 0.0; }
-  /// Read the bias configuration from a restart file
+  virtual std::string const get_state_params() const;
-  virtual std::istream & read_restart(std::istream &is);
+  virtual int set_state_params(std::string const &conf);
  // virtual std::ostream & write_state_data(std::ostream &os);
  // virtual std::istream & read_state_data(std::istream &os);
  virtual std::ostream & write_state(std::ostream &os);
  virtual std::istream & read_state(std::istream &is);
  /// Write the bias configuration to a restart file
  virtual std::ostream & write_restart(std::ostream &os);
  /// Write a label to the trajectory file (comment line)
  virtual std::ostream & write_traj_label(std::ostream &os);
  /// Output quantities such as the bias energy to the trajectory file
  virtual std::ostream & write_traj(std::ostream &os);
  /// \brief Constructor
@ -42,26 +49,110 @@ public:
 protected:
-  /// \brief Potential function
+  /// \brief Potential function for the i-th colvar
-  virtual cvm::real restraint_potential(cvm::real k, colvar const *x,
+  virtual cvm::real restraint_potential(size_t i) const = 0;
                                        colvarvalue const &xcenter) const = 0;
-  /// \brief Force function
+  /// \brief Force function for the i-th colvar
-  virtual colvarvalue restraint_force(cvm::real k, colvar const *x,
+  virtual colvarvalue const restraint_force(size_t i) const = 0;
                                      colvarvalue const &xcenter) const = 0;
-  ///\brief Unit scaling
+  /// \brief Derivative of the potential function with respect to the force constant
-  virtual cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const = 0;
+  virtual cvm::real d_restraint_potential_dk(size_t i) const = 0;
 };
 /// Definition and parsing of the restraint centers
 class colvarbias_restraint_centers
  : public virtual colvarbias_restraint
 {
 public:
  colvarbias_restraint_centers(char const *key);
  virtual int init(std::string const &conf);
  virtual int change_configuration(std::string const &conf);
 protected:
  /// \brief Restraint centers
  std::vector<colvarvalue> colvar_centers;
-  /// \brief Restraint centers without wrapping or constraints applied
+  /// \brief Restraint centers outside the domain of the colvars (no wrapping or constraints applied)
  std::vector<colvarvalue> colvar_centers_raw;
 };
 /// Definition and parsing of the force constant
 class colvarbias_restraint_k
  : public virtual colvarbias_restraint
 {
 public:
  colvarbias_restraint_k(char const *key);
  virtual int init(std::string const &conf);
  virtual int change_configuration(std::string const &conf);
 protected:
  /// \brief Restraint force constant
  cvm::real force_k;
 };
 /// Options to change the restraint configuration over time (shared between centers and k moving)
 class colvarbias_restraint_moving
  : public virtual colvarparse {
 public:
  colvarbias_restraint_moving(char const *key);
  // Note: despite the diamond inheritance, most of this function gets only executed once
  virtual int init(std::string const &conf);
  virtual int update() { return COLVARS_OK; }
  virtual int change_configuration(std::string const &conf) { return COLVARS_NOT_IMPLEMENTED; }
  virtual std::string const get_state_params() const;
  virtual int set_state_params(std::string const &conf);
 protected:
  /// \brief Moving target?
  bool b_chg_centers;
  /// \brief Changing force constant?
  bool b_chg_force_k;
  /// \brief Number of stages over which to perform the change
  /// If zero, perform a continuous change
  int target_nstages;
  /// \brief Number of current stage of the perturbation
  int stage;
    /// \brief Lambda-schedule for custom varying force constant
  std::vector<cvm::real> lambda_schedule;
  /// \brief Number of steps required to reach the target force constant
  /// or restraint centers
  long target_nsteps;
 };
 /// Options to change the restraint centers over time
 class colvarbias_restraint_centers_moving
  : public virtual colvarbias_restraint_centers,
    public virtual colvarbias_restraint_moving
 {
 public:
  colvarbias_restraint_centers_moving(char const *key);
  virtual int init(std::string const &conf);
  virtual int update();
  virtual int change_configuration(std::string const &conf) { return COLVARS_NOT_IMPLEMENTED; }
  virtual std::string const get_state_params() const;
  virtual int set_state_params(std::string const &conf);
  virtual std::ostream & write_traj_label(std::ostream &os);
  virtual std::ostream & write_traj(std::ostream &os);
 protected:
  /// \brief New restraint centers
  std::vector<colvarvalue> target_centers;
@ -78,11 +169,29 @@ protected:
  /// \brief Accumulated work
  cvm::real acc_work;
-  /// \brief Restraint force constant
+  /// Update the accumulated work
-  cvm::real force_k;
+  int update_acc_work();
 };
-  /// \brief Changing force constant?
+
-  bool b_chg_force_k;
+/// Options to change the restraint force constant over time
 class colvarbias_restraint_k_moving
  : public virtual colvarbias_restraint_k,
    public virtual colvarbias_restraint_moving
 {
 public:
  colvarbias_restraint_k_moving(char const *key);
  virtual int init(std::string const &conf);
  virtual int update();
  virtual int change_configuration(std::string const &conf) { return COLVARS_NOT_IMPLEMENTED; }
  virtual std::string const get_state_params() const;
  virtual int set_state_params(std::string const &conf);
  virtual std::ostream & write_traj_label(std::ostream &os);
  virtual std::ostream & write_traj(std::ostream &os);
 protected:
  /// \brief Restraint force constant (target value)
  cvm::real target_force_k;
@ -90,9 +199,6 @@ protected:
  /// \brief Restraint force constant (starting value)
  cvm::real starting_force_k;
  /// \brief Lambda-schedule for custom varying force constant
  std::vector<cvm::real> lambda_schedule;
  /// \brief Exponent for varying the force constant
  cvm::real force_k_exp;
@ -100,71 +206,91 @@ protected:
  /// (in TI, would be the accumulating FE derivative)
  cvm::real restraint_FE;
  /// \brief Equilibration steps for restraint FE calculation through TI
  cvm::real target_equil_steps;
  /// \brief Number of stages over which to perform the change
  /// If zero, perform a continuous change
  int target_nstages;
  /// \brief Number of current stage of the perturbation
  int stage;
  /// \brief Number of steps required to reach the target force constant
  /// or restraint centers
  long target_nsteps;
 };
 /// \brief Harmonic bias restraint
 /// (implementation of \link colvarbias_restraint \endlink)
-class colvarbias_restraint_harmonic : public colvarbias_restraint {
+class colvarbias_restraint_harmonic
-
+  : public colvarbias_restraint_centers_moving,
    public colvarbias_restraint_k_moving
 {
 public:
  colvarbias_restraint_harmonic(char const *key);
  virtual int init(std::string const &conf);
-  // no additional members, destructor not needed
+  virtual int update();
  virtual std::string const get_state_params() const;
  virtual int set_state_params(std::string const &conf);
  virtual std::ostream & write_traj_label(std::ostream &os);
  virtual std::ostream & write_traj(std::ostream &os);
  virtual int change_configuration(std::string const &conf);
  virtual cvm::real energy_difference(std::string const &conf);
 protected:
-  /// \brief Potential function
+  virtual cvm::real restraint_potential(size_t i) const;
-  virtual cvm::real restraint_potential(cvm::real k, colvar const *x,
+  virtual colvarvalue const restraint_force(size_t i) const;
-                                        colvarvalue const &xcenter) const;
+  virtual cvm::real d_restraint_potential_dk(size_t i) const;
 };
  /// \brief Force function
  virtual colvarvalue restraint_force(cvm::real k, colvar const *x,
                                      colvarvalue const &xcenter) const;
-  ///\brief Unit scaling
+/// \brief Wall restraint
-  virtual cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const;
+/// (implementation of \link colvarbias_restraint \endlink)
 class colvarbias_restraint_harmonic_walls
  : public colvarbias_restraint_k_moving
 {
 public:
  colvarbias_restraint_harmonic_walls(char const *key);
  virtual int init(std::string const &conf);
  virtual int update();
  virtual void communicate_forces();
  virtual std::string const get_state_params() const;
  virtual int set_state_params(std::string const &conf);
  virtual std::ostream & write_traj_label(std::ostream &os);
  virtual std::ostream & write_traj(std::ostream &os);
 protected:
  /// \brief Location of the lower walls
  std::vector<colvarvalue> lower_walls;
  /// \brief Location of the upper walls
  std::vector<colvarvalue> upper_walls;
  virtual cvm::real colvar_distance(size_t i) const;
  virtual cvm::real restraint_potential(size_t i) const;
  virtual colvarvalue const restraint_force(size_t i) const;
  virtual cvm::real d_restraint_potential_dk(size_t i) const;
 };
 /// \brief Linear bias restraint
 /// (implementation of \link colvarbias_restraint \endlink)
-class colvarbias_restraint_linear : public colvarbias_restraint {
+class colvarbias_restraint_linear
  : public colvarbias_restraint_centers_moving,
    public colvarbias_restraint_k_moving
 {
 public:
  colvarbias_restraint_linear(char const *key);
  virtual int init(std::string const &conf);
-  // no additional members, destructor not needed
+  virtual int update();
  virtual int change_configuration(std::string const &conf);
  virtual cvm::real energy_difference(std::string const &conf);
  virtual std::string const get_state_params() const;
  virtual int set_state_params(std::string const &conf);
  virtual std::ostream & write_traj_label(std::ostream &os);
  virtual std::ostream & write_traj(std::ostream &os);
 protected:
-  /// \brief Potential function
+  virtual cvm::real restraint_potential(size_t i) const;
-  virtual cvm::real restraint_potential(cvm::real k, colvar const *x,
+  virtual colvarvalue const restraint_force(size_t i) const;
-                                        colvarvalue const &xcenter) const;
+  virtual cvm::real d_restraint_potential_dk(size_t i) const;
  /// \brief Force function
  virtual colvarvalue restraint_force(cvm::real k, colvar const *x,
                                      colvarvalue const &xcenter) const;
  ///\brief Unit scaling
  virtual cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const;
 };
--- a/lib/colvars/colvarcomp.cpp
+++ b/lib/colvars/colvarcomp.cpp
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include "colvarmodule.h"
 #include "colvarvalue.h"
 #include "colvar.h"
@ -154,15 +161,17 @@ void colvar::cvc::read_data()
 void colvar::cvc::calc_force_invgrads()
 {
-  cvm::fatal_error("Error: calculation of inverse gradients is not implemented "
+  cvm::error("Error: calculation of inverse gradients is not implemented "
-                    "for colvar components of type \""+function_type+"\".\n");
+             "for colvar components of type \""+function_type+"\".\n",
             COLVARS_NOT_IMPLEMENTED);
 }
 void colvar::cvc::calc_Jacobian_derivative()
 {
-  cvm::fatal_error("Error: calculation of inverse gradients is not implemented "
+  cvm::error("Error: calculation of inverse gradients is not implemented "
-                    "for colvar components of type \""+function_type+"\".\n");
+             "for colvar components of type \""+function_type+"\".\n",
             COLVARS_NOT_IMPLEMENTED);
 }
@ -281,6 +290,33 @@ void colvar::cvc::debug_gradients(cvm::atom_group *group)
 }
 cvm::real colvar::cvc::dist2(colvarvalue const &x1,
                             colvarvalue const &x2) const
 {
  return x1.dist2(x2);
 }
 colvarvalue colvar::cvc::dist2_lgrad(colvarvalue const &x1,
                                     colvarvalue const &x2) const
 {
  return x1.dist2_grad(x2);
 }
 colvarvalue colvar::cvc::dist2_rgrad(colvarvalue const &x1,
                                     colvarvalue const &x2) const
 {
  return x2.dist2_grad(x1);
 }
 void colvar::cvc::wrap(colvarvalue &x) const
 {
  return;
 }
 // Static members
 std::vector<colvardeps::feature *> colvar::cvc::cvc_features;
--- a/lib/colvars/colvarcomp.h
+++ b/lib/colvars/colvarcomp.h
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #ifndef COLVARCOMP_H
 #define COLVARCOMP_H
@ -105,8 +112,8 @@ public:
  /// options from the provided configuration string
  /// Returns reference to new group
  cvm::atom_group *parse_group(std::string const &conf,
-                   char const *group_key,
+                               char const *group_key,
-                   bool optional = false);
+                               bool optional = false);
  /// \brief Parse options pertaining to total force calculation
  virtual int init_total_force_params(std::string const &conf);
@ -130,7 +137,7 @@ public:
  }
  /// \brief Obtain data needed for the calculation for the backend
-  void read_data();
+  virtual void read_data();
  /// \brief Calculate the variable
  virtual void calc_value() = 0;
@ -151,17 +158,14 @@ public:
  /// \brief Return the previously calculated value
-  virtual colvarvalue const & value() const;
+  colvarvalue const & value() const;
  // /// \brief Return const pointer to the previously calculated value
  // virtual const colvarvalue *p_value() const;
  /// \brief Return the previously calculated total force
-  virtual colvarvalue const & total_force() const;
+  colvarvalue const & total_force() const;
  /// \brief Return the previously calculated divergence of the
  /// inverse atomic gradients
-  virtual colvarvalue const & Jacobian_derivative() const;
+  colvarvalue const & Jacobian_derivative() const;
  /// \brief Apply the collective variable force, by communicating the
  /// atomic forces to the simulation program (\b Note: the \link ft
@ -247,52 +251,24 @@ protected:
 };
 inline colvarvalue const & colvar::cvc::value() const
 {
  return x;
 }
 // inline const colvarvalue * colvar::cvc::p_value() const
 // {
 //   return &x;
 // }
 inline colvarvalue const & colvar::cvc::total_force() const
 {
  return ft;
 }
 inline colvarvalue const & colvar::cvc::Jacobian_derivative() const
 {
  return jd;
 }
 inline cvm::real colvar::cvc::dist2(colvarvalue const &x1,
                                    colvarvalue const &x2) const
 {
  return x1.dist2(x2);
 }
 inline colvarvalue colvar::cvc::dist2_lgrad(colvarvalue const &x1,
                                            colvarvalue const &x2) const
 {
  return x1.dist2_grad(x2);
 }
 inline colvarvalue colvar::cvc::dist2_rgrad(colvarvalue const &x1,
                                            colvarvalue const &x2) const
 {
  return x2.dist2_grad(x1);
 }
 inline void colvar::cvc::wrap(colvarvalue &x) const
 {
  return;
 }
 /// \brief Colvar component: distance between the centers of mass of
 /// two groups (colvarvalue::type_scalar type, range [0:*))
@ -312,7 +288,7 @@ protected:
 public:
  distance(std::string const &conf);
  distance();
-  virtual inline ~distance() {}
+  virtual ~distance() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void calc_force_invgrads();
@ -327,6 +303,7 @@ public:
 };
 // \brief Colvar component: distance vector between centers of mass
 // of two groups (\link colvarvalue::type_3vector \endlink type,
 // range (-*:*)x(-*:*)x(-*:*))
@ -336,7 +313,7 @@ class colvar::distance_vec
 public:
  distance_vec(std::string const &conf);
  distance_vec();
-  virtual inline ~distance_vec() {}
+  virtual ~distance_vec() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force(colvarvalue const &force);
@ -352,6 +329,7 @@ public:
 };
 /// \brief Colvar component: distance unit vector (direction) between
 /// centers of mass of two groups (colvarvalue::type_unit3vector type,
 /// range [-1:1]x[-1:1]x[-1:1])
@ -361,19 +339,14 @@ class colvar::distance_dir
 public:
  distance_dir(std::string const &conf);
  distance_dir();
-  virtual inline ~distance_dir() {}
+  virtual ~distance_dir() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force(colvarvalue const &force);
  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: projection of the distance vector along
 /// an axis(colvarvalue::type_scalar type, range (-*:*))
 class colvar::distance_z
@ -399,7 +372,7 @@ protected:
 public:
  distance_z(std::string const &conf);
  distance_z();
-  virtual inline ~distance_z() {}
+  virtual ~distance_z() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void calc_force_invgrads();
@ -416,6 +389,7 @@ public:
 };
 /// \brief Colvar component: projection of the distance vector on a
 /// plane (colvarvalue::type_scalar type, range [0:*))
 class colvar::distance_xy
@ -429,7 +403,7 @@ protected:
 public:
  distance_xy(std::string const &conf);
  distance_xy();
-  virtual inline ~distance_xy() {}
+  virtual ~distance_xy() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void calc_force_invgrads();
@ -444,6 +418,7 @@ public:
 };
 /// \brief Colvar component: average distance between two groups of atoms, weighted as the sixth power,
 /// as in NMR refinements(colvarvalue::type_scalar type, range (0:*))
 class colvar::distance_inv
@ -455,7 +430,7 @@ protected:
 public:
  distance_inv(std::string const &conf);
  distance_inv();
-  virtual inline ~distance_inv() {}
+  virtual ~distance_inv() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force(colvarvalue const &force);
@ -468,6 +443,7 @@ public:
 };
 /// \brief Colvar component: N1xN2 vector of pairwise distances
 /// (colvarvalue::type_vector type, range (0:*) for each component)
 class colvar::distance_pairs
@ -483,16 +459,10 @@ protected:
 public:
  distance_pairs(std::string const &conf);
  distance_pairs();
-  virtual inline ~distance_pairs() {}
+  virtual ~distance_pairs() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force(colvarvalue const &force);
  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;
 };
@ -509,7 +479,7 @@ public:
  /// Constructor
  gyration(std::string const &conf);
  gyration();
-  virtual inline ~gyration() {}
+  virtual ~gyration() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void calc_force_invgrads();
@ -524,6 +494,7 @@ public:
 };
 /// \brief Colvar component: moment of inertia of an atom group
 /// (colvarvalue::type_scalar type, range [0:*))
 class colvar::inertia
@ -533,7 +504,7 @@ public:
  /// Constructor
  inertia(std::string const &conf);
  inertia();
-  virtual inline ~inertia() {}
+  virtual ~inertia() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force(colvarvalue const &force);
@ -546,6 +517,7 @@ public:
 };
 /// \brief Colvar component: moment of inertia of an atom group
 /// around a user-defined axis (colvarvalue::type_scalar type, range [0:*))
 class colvar::inertia_z
@ -558,7 +530,7 @@ public:
  /// Constructor
  inertia_z(std::string const &conf);
  inertia_z();
-  virtual inline ~inertia_z() {}
+  virtual ~inertia_z() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force(colvarvalue const &force);
@ -571,6 +543,7 @@ public:
 };
 /// \brief Colvar component: projection of 3N coordinates onto an
 /// eigenvector(colvarvalue::type_scalar type, range (-*:*))
 class colvar::eigenvector
@ -597,7 +570,7 @@ public:
  /// Constructor
  eigenvector(std::string const &conf);
-  virtual inline ~eigenvector() {}
+  virtual ~eigenvector() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void calc_force_invgrads();
@ -645,7 +618,7 @@ public:
  /// \brief Initialize the three groups after three atoms
  angle(cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3);
  angle();
-  virtual inline ~angle() {}
+  virtual ~angle() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void calc_force_invgrads();
@ -659,6 +632,8 @@ public:
                                  colvarvalue const &x2) const;
 };
 /// \brief Colvar component: angle between the dipole of a molecule and an axis
 /// formed by two groups of atoms(colvarvalue::type_scalar type, range [0:PI])
 class colvar::dipole_angle
@ -691,7 +666,7 @@ public:
  /// \brief Initialize the three groups after three atoms
  dipole_angle (cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3);
  dipole_angle();
-  virtual inline ~dipole_angle() {}
+  virtual ~dipole_angle() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force (colvarvalue const &force);
@ -703,6 +678,8 @@ public:
                                   colvarvalue const &x2) const;
 };
 /// \brief Colvar component: dihedral between the centers of mass of
 /// four groups (colvarvalue::type_scalar type, range [-PI:PI])
 class colvar::dihedral
@ -732,7 +709,7 @@ public:
  /// \brief Initialize the four groups after four atoms
  dihedral(cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3, cvm::atom const &a4);
  dihedral();
-  virtual inline ~dihedral() {}
+  virtual ~dihedral() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void calc_force_invgrads();
@ -753,6 +730,7 @@ public:
 };
 /// \brief Colvar component: coordination number between two groups
 /// (colvarvalue::type_scalar type, range [0:N1*N2])
 class colvar::coordnum
@ -781,7 +759,7 @@ public:
  /// Constructor
  coordnum(std::string const &conf);
  coordnum();
-  virtual inline ~coordnum() {}
+  virtual ~coordnum() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force(colvarvalue const &force);
@ -812,6 +790,8 @@ public:
                                  colvarvalue const &x2) const;
 };
 /// \brief Colvar component: self-coordination number within a group
 /// (colvarvalue::type_scalar type, range [0:N*(N-1)/2])
 class colvar::selfcoordnum
@ -830,7 +810,7 @@ public:
  /// Constructor
  selfcoordnum(std::string const &conf);
  selfcoordnum();
-  virtual inline ~selfcoordnum() {}
+  virtual ~selfcoordnum() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force(colvarvalue const &force);
@ -852,6 +832,7 @@ public:
 };
 /// \brief Colvar component: coordination number between two groups
 /// (colvarvalue::type_scalar type, range [0:N1*N2])
 class colvar::groupcoordnum
@ -873,7 +854,7 @@ public:
  /// Constructor
  groupcoordnum(std::string const &conf);
  groupcoordnum();
-  virtual inline ~groupcoordnum() {}
+  virtual ~groupcoordnum() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force(colvarvalue const &force);
@ -896,6 +877,7 @@ public:
  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;
@ -903,10 +885,10 @@ public:
                                  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])
@ -941,6 +923,7 @@ public:
 };
 /// \brief Colvar component: alpha helix content of a contiguous
 /// segment of 5 or more residues, implemented as a sum of phi/psi
 /// dihedral angles and hydrogen bonds (colvarvalue::type_scalar type,
@ -969,7 +952,7 @@ public:
 //   alpha_dihedrals (std::string const &conf);
 //   alpha_dihedrals();
-//   virtual inline ~alpha_dihedrals() {}
+//   virtual ~alpha_dihedrals() {}
 //   virtual void calc_value();
 //   virtual void calc_gradients();
 //   virtual void apply_force (colvarvalue const &force);
@ -982,6 +965,7 @@ public:
 // };
 /// \brief Colvar component: alpha helix content of a contiguous
 /// segment of 5 or more residues, implemented as a sum of Ca-Ca-Ca
 /// angles and hydrogen bonds (colvarvalue::type_scalar type, range
@ -1022,6 +1006,8 @@ public:
                                  colvarvalue const &x2) const;
 };
 /// \brief Colvar component: dihedPC
 /// Projection of the config onto a dihedral principal component
 /// See e.g. Altis et al., J. Chem. Phys 126, 244111 (2007)
@ -1050,6 +1036,8 @@ public:
                                  colvarvalue const &x2) const;
 };
 /// \brief Colvar component: orientation in space of an atom group,
 /// with respect to a set of reference coordinates
 /// (colvarvalue::type_quaternion type, range
@ -1078,7 +1066,7 @@ public:
  orientation(std::string const &conf);
  orientation();
-  virtual inline ~orientation() {}
+  virtual ~orientation() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force(colvarvalue const &force);
@ -1091,6 +1079,7 @@ public:
 };
 /// \brief Colvar component: angle of rotation with respect to a set
 /// of reference coordinates (colvarvalue::type_scalar type, range
 /// [0:PI))
@ -1101,7 +1090,7 @@ public:
  orientation_angle(std::string const &conf);
  orientation_angle();
-  virtual inline ~orientation_angle() {}
+  virtual ~orientation_angle() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force(colvarvalue const &force);
@ -1114,6 +1103,7 @@ public:
 };
 /// \brief Colvar component: cosine of the angle of rotation with respect to a set
 /// of reference coordinates (colvarvalue::type_scalar type, range
 /// [-1:1])
@ -1124,7 +1114,7 @@ public:
  orientation_proj(std::string const &conf);
  orientation_proj();
-  virtual inline ~orientation_proj() {}
+  virtual ~orientation_proj() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force(colvarvalue const &force);
@ -1137,6 +1127,7 @@ public:
 };
 /// \brief Colvar component: projection of the orientation vector onto
 /// a predefined axis (colvarvalue::type_scalar type, range [-1:1])
 class colvar::tilt
@ -1150,7 +1141,7 @@ public:
  tilt(std::string const &conf);
  tilt();
-  virtual inline ~tilt() {}
+  virtual ~tilt() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force(colvarvalue const &force);
@ -1177,7 +1168,7 @@ public:
  spin_angle(std::string const &conf);
  spin_angle();
-  virtual inline ~spin_angle() {}
+  virtual ~spin_angle() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force(colvarvalue const &force);
@ -1215,7 +1206,7 @@ public:
  /// Constructor
  rmsd(std::string const &conf);
-  virtual inline ~rmsd() {}
+  virtual ~rmsd() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void calc_force_invgrads();
@ -1230,6 +1221,7 @@ public:
 };
 // \brief Colvar component: flat vector of Cartesian coordinates
 // Mostly useful to compute scripted colvar values
 class colvar::cartesian
@ -1243,7 +1235,7 @@ protected:
 public:
  cartesian(std::string const &conf);
  cartesian();
-  virtual inline ~cartesian() {}
+  virtual ~cartesian() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force(colvarvalue const &force);
@ -1260,255 +1252,26 @@ public:
 #define simple_scalar_dist_functions(TYPE)                              \
                                                                        \
-  inline cvm::real colvar::TYPE::dist2(colvarvalue const &x1,           \
+                                                                        \
-                                       colvarvalue const &x2) const     \
+  cvm::real colvar::TYPE::dist2(colvarvalue const &x1,                  \
                                colvarvalue const &x2) const            \
  {                                                                     \
    return (x1.real_value - x2.real_value)*(x1.real_value - x2.real_value); \
  }                                                                     \
                                                                        \
-  inline colvarvalue colvar::TYPE::dist2_lgrad(colvarvalue const &x1,   \
+                                                                        \
-                                               colvarvalue const &x2) const \
+  colvarvalue colvar::TYPE::dist2_lgrad(colvarvalue const &x1,          \
                                        colvarvalue const &x2) const    \
  {                                                                     \
    return 2.0 * (x1.real_value - x2.real_value);                       \
  }                                                                     \
                                                                        \
-  inline colvarvalue colvar::TYPE::dist2_rgrad(colvarvalue const &x1,   \
+                                                                        \
-                                               colvarvalue const &x2) const \
+  colvarvalue colvar::TYPE::dist2_rgrad(colvarvalue const &x1,          \
                                        colvarvalue const &x2) const    \
  {                                                                     \
    return this->dist2_lgrad(x2, x1);                                   \
  }                                                                     \
                                                                        \
 simple_scalar_dist_functions(distance)
 // NOTE: distance_z has explicit functions, see below
 simple_scalar_dist_functions(distance_xy)
 simple_scalar_dist_functions(distance_inv)
 simple_scalar_dist_functions(angle)
 simple_scalar_dist_functions(dipole_angle)
 simple_scalar_dist_functions(coordnum)
 simple_scalar_dist_functions(selfcoordnum)
 simple_scalar_dist_functions(h_bond)
 simple_scalar_dist_functions(gyration)
 simple_scalar_dist_functions(inertia)
 simple_scalar_dist_functions(inertia_z)
 simple_scalar_dist_functions(rmsd)
 simple_scalar_dist_functions(orientation_angle)
 simple_scalar_dist_functions(orientation_proj)
 simple_scalar_dist_functions(tilt)
 simple_scalar_dist_functions(eigenvector)
 //  simple_scalar_dist_functions (alpha_dihedrals)
 simple_scalar_dist_functions(alpha_angles)
 simple_scalar_dist_functions(dihedPC)
 // metrics functions for cvc implementations with a periodicity
 inline cvm::real colvar::dihedral::dist2(colvarvalue const &x1,
                                           colvarvalue const &x2) const
 {
 cvm::real diff = x1.real_value - x2.real_value;
 diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
 return diff * diff;
 }
 inline colvarvalue colvar::dihedral::dist2_lgrad(colvarvalue const &x1,
                                                   colvarvalue const &x2) const
 {
 cvm::real diff = x1.real_value - x2.real_value;
 diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
 return 2.0 * diff;
 }
 inline colvarvalue colvar::dihedral::dist2_rgrad(colvarvalue const &x1,
                                                   colvarvalue const &x2) const
 {
 cvm::real diff = x1.real_value - x2.real_value;
 diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
 return (-2.0) * diff;
 }
 inline void colvar::dihedral::wrap(colvarvalue &x) const
 {
 if ((x.real_value - wrap_center) >= 180.0) {
 x.real_value -= 360.0;
 return;
 }
 if ((x.real_value - wrap_center) < -180.0) {
 x.real_value += 360.0;
 return;
 }
 return;
 }
 inline cvm::real colvar::spin_angle::dist2(colvarvalue const &x1,
                                             colvarvalue const &x2) const
 {
 cvm::real diff = x1.real_value - x2.real_value;
 diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
 return diff * diff;
 }
 inline colvarvalue colvar::spin_angle::dist2_lgrad(colvarvalue const &x1,
                                                     colvarvalue const &x2) const
 {
 cvm::real diff = x1.real_value - x2.real_value;
 diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
 return 2.0 * diff;
 }
 inline colvarvalue colvar::spin_angle::dist2_rgrad(colvarvalue const &x1,
                                                     colvarvalue const &x2) const
 {
 cvm::real diff = x1.real_value - x2.real_value;
 diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
 return (-2.0) * diff;
 }
 inline void colvar::spin_angle::wrap(colvarvalue &x) const
 {
 if ((x.real_value - wrap_center) >= 180.0) {
 x.real_value -= 360.0;
 return;
 }
 if ((x.real_value - wrap_center) < -180.0) {
 x.real_value += 360.0;
 return;
 }
 return;
 }
 // Projected distance
 // Differences should always be wrapped around 0 (ignoring wrap_center)
 inline cvm::real colvar::distance_z::dist2(colvarvalue const &x1,
                                             colvarvalue const &x2) const
 {
  cvm::real diff = x1.real_value - x2.real_value;
  if (period != 0.0) {
    cvm::real shift = std::floor(diff/period + 0.5);
    diff -= shift * period;
  }
  return diff * diff;
 }
 inline colvarvalue colvar::distance_z::dist2_lgrad(colvarvalue const &x1,
                                                   colvarvalue const &x2) const
 {
  cvm::real diff = x1.real_value - x2.real_value;
  if (period != 0.0) {
    cvm::real shift = std::floor(diff/period + 0.5);
    diff -= shift * period;
  }
  return 2.0 * diff;
 }
 inline colvarvalue colvar::distance_z::dist2_rgrad(colvarvalue const &x1,
                                                   colvarvalue const &x2) const
 {
  cvm::real diff = x1.real_value - x2.real_value;
  if (period != 0.0) {
    cvm::real shift = std::floor(diff/period + 0.5);
    diff -= shift * period;
  }
  return (-2.0) * diff;
 }
 inline void colvar::distance_z::wrap(colvarvalue &x) const
 {
  if (! this->b_periodic) {
    // don't wrap if the period has not been set
    return;
  }
  cvm::real shift = std::floor((x.real_value - wrap_center) / period + 0.5);
  x.real_value -= shift * period;
  return;
 }
 // distance between three dimensional vectors
 //
 // TODO apply PBC to distance_vec
 // Note: differences should be centered around (0, 0, 0)!
 inline cvm::real colvar::distance_vec::dist2(colvarvalue const &x1,
                                             colvarvalue const &x2) const
 {
  return cvm::position_dist2(x1.rvector_value, x2.rvector_value);
 }
 inline colvarvalue colvar::distance_vec::dist2_lgrad(colvarvalue const &x1,
                                                     colvarvalue const &x2) const
 {
  return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
 }
 inline colvarvalue colvar::distance_vec::dist2_rgrad(colvarvalue const &x1,
                                                     colvarvalue const &x2) const
 {
  return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
 }
 inline cvm::real colvar::distance_dir::dist2(colvarvalue const &x1,
                                             colvarvalue const &x2) const
 {
  return (x1.rvector_value - x2.rvector_value).norm2();
 }
 inline colvarvalue colvar::distance_dir::dist2_lgrad(colvarvalue const &x1,
                                                     colvarvalue const &x2) const
 {
  return colvarvalue((x1.rvector_value - x2.rvector_value), colvarvalue::type_unit3vector);
 }
 inline colvarvalue colvar::distance_dir::dist2_rgrad(colvarvalue const &x1,
                                                     colvarvalue const &x2) const
 {
  return colvarvalue((x2.rvector_value - x1.rvector_value), colvarvalue::type_unit3vector);
 }
 inline cvm::real colvar::distance_pairs::dist2(colvarvalue const &x1,
                                               colvarvalue const &x2) const
 {
  return (x1.vector1d_value - x2.vector1d_value).norm2();
 }
 inline colvarvalue colvar::distance_pairs::dist2_lgrad(colvarvalue const &x1,
                                                       colvarvalue const &x2) const
 {
  return colvarvalue((x1.vector1d_value - x2.vector1d_value), colvarvalue::type_vector);
 }
 inline colvarvalue colvar::distance_pairs::dist2_rgrad(colvarvalue const &x1,
                                                       colvarvalue const &x2) const
 {
  return colvarvalue((x2.vector1d_value - x1.vector1d_value), colvarvalue::type_vector);
 }
 // distance between quaternions
 inline cvm::real colvar::orientation::dist2(colvarvalue const &x1,
                                            colvarvalue const &x2) const
 {
  return x1.quaternion_value.dist2(x2);
 }
 inline colvarvalue colvar::orientation::dist2_lgrad(colvarvalue const &x1,
                                                    colvarvalue const &x2) const
 {
  return x1.quaternion_value.dist2_grad(x2);
 }
 inline colvarvalue colvar::orientation::dist2_rgrad(colvarvalue const &x1,
                                                    colvarvalue const &x2) const
 {
  return x2.quaternion_value.dist2_grad(x1);
 }
 #endif
--- a/lib/colvars/colvarcomp_angles.cpp
+++ b/lib/colvars/colvarcomp_angles.cpp
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include "colvarmodule.h"
 #include "colvar.h"
 #include "colvarcomp.h"
@ -7,6 +14,7 @@
 #include <cmath>
 colvar::angle::angle(std::string const &conf)
  : cvc(conf)
 {
@ -85,6 +93,7 @@ void colvar::angle::calc_gradients()
  group3->set_weighted_gradient(dxdr3);
 }
 void colvar::angle::calc_force_invgrads()
 {
  // This uses a force measurement on groups 1 and 3 only
@ -107,6 +116,7 @@ void colvar::angle::calc_force_invgrads()
  return;
 }
 void colvar::angle::calc_Jacobian_derivative()
 {
  // det(J) = (2 pi) r^2 * sin(theta)
@ -129,6 +139,8 @@ void colvar::angle::apply_force(colvarvalue const &force)
 }
 simple_scalar_dist_functions(angle)
 colvar::dipole_angle::dipole_angle(std::string const &conf)
@ -235,6 +247,8 @@ void colvar::dipole_angle::apply_force(colvarvalue const &force)
 }
 simple_scalar_dist_functions(dipole_angle)
 colvar::dihedral::dihedral(std::string const &conf)
@ -453,3 +467,46 @@ void colvar::dihedral::apply_force(colvarvalue const &force)
 }
 // metrics functions for cvc implementations with a periodicity
 cvm::real colvar::dihedral::dist2(colvarvalue const &x1,
                                  colvarvalue const &x2) const
 {
  cvm::real diff = x1.real_value - x2.real_value;
  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
  return diff * diff;
 }
 colvarvalue colvar::dihedral::dist2_lgrad(colvarvalue const &x1,
                                          colvarvalue const &x2) const
 {
  cvm::real diff = x1.real_value - x2.real_value;
  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
  return 2.0 * diff;
 }
 colvarvalue colvar::dihedral::dist2_rgrad(colvarvalue const &x1,
                                          colvarvalue const &x2) const
 {
  cvm::real diff = x1.real_value - x2.real_value;
  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
  return (-2.0) * diff;
 }
 void colvar::dihedral::wrap(colvarvalue &x) const
 {
  if ((x.real_value - wrap_center) >= 180.0) {
    x.real_value -= 360.0;
    return;
  }
  if ((x.real_value - wrap_center) < -180.0) {
    x.real_value += 360.0;
    return;
  }
  return;
 }
--- a/lib/colvars/colvarcomp_coordnums.cpp
+++ b/lib/colvars/colvarcomp_coordnums.cpp
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include <cmath>
 #include "colvarmodule.h"
@ -13,10 +20,10 @@
 template<bool calculate_gradients>
 cvm::real colvar::coordnum::switching_function(cvm::real const &r0,
-                                                int const &en,
+                                               int const &en,
-                                                int const &ed,
+                                               int const &ed,
-                                                cvm::atom &A1,
+                                               cvm::atom &A1,
-                                                cvm::atom &A2)
+                                               cvm::atom &A2)
 {
  cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
  cvm::real const l2 = diff.norm2()/(r0*r0);
@ -42,10 +49,10 @@ cvm::real colvar::coordnum::switching_function(cvm::real const &r0,
 template<bool calculate_gradients>
 cvm::real colvar::coordnum::switching_function(cvm::rvector const &r0_vec,
-                                                int const &en,
+                                               int const &en,
-                                                int const &ed,
+                                               int const &ed,
-                                                cvm::atom &A1,
+                                               cvm::atom &A1,
-                                                cvm::atom &A2)
+                                               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);
@ -190,6 +197,7 @@ void colvar::coordnum::calc_gradients()
  }
 }
 void colvar::coordnum::apply_force(colvarvalue const &force)
 {
  if (!group1->noforce)
@ -200,6 +208,9 @@ void colvar::coordnum::apply_force(colvarvalue const &force)
 }
 simple_scalar_dist_functions(coordnum)
 // h_bond member functions
@ -252,6 +263,7 @@ colvar::h_bond::h_bond(cvm::atom const &acceptor,
  atom_groups[0]->add_atom(donor);
 }
 colvar::h_bond::h_bond()
  : cvc()
 {
@ -284,6 +296,8 @@ void colvar::h_bond::apply_force(colvarvalue const &force)
 }
 simple_scalar_dist_functions(h_bond)
 colvar::selfcoordnum::selfcoordnum(std::string const &conf)
@ -339,6 +353,9 @@ void colvar::selfcoordnum::apply_force(colvarvalue const &force)
 }
 simple_scalar_dist_functions(selfcoordnum)
 // groupcoordnum member functions
 colvar::groupcoordnum::groupcoordnum(std::string const &conf)
  : distance(conf), b_anisotropic(false)
@ -448,7 +465,6 @@ cvm::real colvar::groupcoordnum::switching_function(cvm::rvector const &r0_vec,
 #endif
 void colvar::groupcoordnum::calc_value()
 {
@ -460,7 +476,6 @@ void colvar::groupcoordnum::calc_value()
  x.real_value = coordnum::switching_function<false>(r0, en, ed,
                                                     group1_com_atom, group2_com_atom);
 }
@ -486,3 +501,6 @@ void colvar::groupcoordnum::apply_force(colvarvalue const &force)
  if (!group2->noforce)
    group2->apply_colvar_force(force.real_value);
 }
 simple_scalar_dist_functions(groupcoordnum)
--- a/lib/colvars/colvarcomp_distances.cpp
+++ b/lib/colvars/colvarcomp_distances.cpp
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include <cmath>
 #include "colvarmodule.h"
@ -91,6 +98,9 @@ void colvar::distance::apply_force(colvarvalue const &force)
 }
 simple_scalar_dist_functions(distance)
 colvar::distance_vec::distance_vec(std::string const &conf)
  : distance(conf)
@ -138,6 +148,27 @@ void colvar::distance_vec::apply_force(colvarvalue const &force)
 }
 cvm::real colvar::distance_vec::dist2(colvarvalue const &x1,
                                      colvarvalue const &x2) const
 {
  return cvm::position_dist2(x1.rvector_value, x2.rvector_value);
 }
 colvarvalue colvar::distance_vec::dist2_lgrad(colvarvalue const &x1,
                                              colvarvalue const &x2) const
 {
  return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
 }
 colvarvalue colvar::distance_vec::dist2_rgrad(colvarvalue const &x1,
                                              colvarvalue const &x2) const
 {
  return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
 }
 colvar::distance_z::distance_z(std::string const &conf)
  : cvc(conf)
@ -191,6 +222,7 @@ colvar::distance_z::distance_z(std::string const &conf)
 }
 colvar::distance_z::distance_z()
 {
  function_type = "distance_z";
@ -200,6 +232,7 @@ colvar::distance_z::distance_z()
  x.type(colvarvalue::type_scalar);
 }
 void colvar::distance_z::calc_value()
 {
  if (fixed_axis) {
@ -227,6 +260,7 @@ void colvar::distance_z::calc_value()
  this->wrap(x);
 }
 void colvar::distance_z::calc_gradients()
 {
  main->set_weighted_gradient( axis );
@ -248,6 +282,7 @@ void colvar::distance_z::calc_gradients()
  }
 }
 void colvar::distance_z::calc_force_invgrads()
 {
  main->read_total_forces();
@ -260,11 +295,13 @@ void colvar::distance_z::calc_force_invgrads()
  }
 }
 void colvar::distance_z::calc_Jacobian_derivative()
 {
  jd.real_value = 0.0;
 }
 void colvar::distance_z::apply_force(colvarvalue const &force)
 {
  if (!ref1->noforce)
@ -278,6 +315,56 @@ void colvar::distance_z::apply_force(colvarvalue const &force)
 }
 // Differences should always be wrapped around 0 (ignoring wrap_center)
 cvm::real colvar::distance_z::dist2(colvarvalue const &x1,
                                    colvarvalue const &x2) const
 {
  cvm::real diff = x1.real_value - x2.real_value;
  if (b_periodic) {
    cvm::real shift = std::floor(diff/period + 0.5);
    diff -= shift * period;
  }
  return diff * diff;
 }
 colvarvalue colvar::distance_z::dist2_lgrad(colvarvalue const &x1,
                                            colvarvalue const &x2) const
 {
  cvm::real diff = x1.real_value - x2.real_value;
  if (b_periodic) {
    cvm::real shift = std::floor(diff/period + 0.5);
    diff -= shift * period;
  }
  return 2.0 * diff;
 }
 colvarvalue colvar::distance_z::dist2_rgrad(colvarvalue const &x1,
                                            colvarvalue const &x2) const
 {
  cvm::real diff = x1.real_value - x2.real_value;
  if (b_periodic) {
    cvm::real shift = std::floor(diff/period + 0.5);
    diff -= shift * period;
  }
  return (-2.0) * diff;
 }
 void colvar::distance_z::wrap(colvarvalue &x) const
 {
  if (!b_periodic) {
    // don't wrap if the period has not been set
    return;
  }
  cvm::real shift = std::floor((x.real_value - wrap_center) / period + 0.5);
  x.real_value -= shift * period;
  return;
 }
 colvar::distance_xy::distance_xy(std::string const &conf)
  : distance_z(conf)
@ -289,6 +376,7 @@ colvar::distance_xy::distance_xy(std::string const &conf)
  x.type(colvarvalue::type_scalar);
 }
 colvar::distance_xy::distance_xy()
  : distance_z()
 {
@ -299,6 +387,7 @@ colvar::distance_xy::distance_xy()
  x.type(colvarvalue::type_scalar);
 }
 void colvar::distance_xy::calc_value()
 {
  if (b_no_PBC) {
@ -321,6 +410,7 @@ void colvar::distance_xy::calc_value()
  x.real_value = dist_v_ortho.norm();
 }
 void colvar::distance_xy::calc_gradients()
 {
  // Intermediate quantity (r_P3 / r_12 where P is the projection
@ -348,6 +438,7 @@ void colvar::distance_xy::calc_gradients()
  }
 }
 void colvar::distance_xy::calc_force_invgrads()
 {
  main->read_total_forces();
@ -360,11 +451,13 @@ void colvar::distance_xy::calc_force_invgrads()
  }
 }
 void colvar::distance_xy::calc_Jacobian_derivative()
 {
  jd.real_value = x.real_value ? (1.0 / x.real_value) : 0.0;
 }
 void colvar::distance_xy::apply_force(colvarvalue const &force)
 {
  if (!ref1->noforce)
@ -378,6 +471,9 @@ void colvar::distance_xy::apply_force(colvarvalue const &force)
 }
 simple_scalar_dist_functions(distance_xy)
 colvar::distance_dir::distance_dir(std::string const &conf)
  : distance(conf)
@ -403,7 +499,7 @@ void colvar::distance_dir::calc_value()
    dist_v = group2->center_of_mass() - group1->center_of_mass();
  } else {
    dist_v = cvm::position_distance(group1->center_of_mass(),
-                                     group2->center_of_mass());
+                                    group2->center_of_mass());
  }
  x.rvector_value = dist_v.unit();
 }
@ -460,6 +556,7 @@ colvar::distance_inv::distance_inv(std::string const &conf)
  x.type(colvarvalue::type_scalar);
 }
 colvar::distance_inv::distance_inv()
 {
  function_type = "distance_inv";
@ -467,6 +564,7 @@ colvar::distance_inv::distance_inv()
  x.type(colvarvalue::type_scalar);
 }
 void colvar::distance_inv::calc_value()
 {
  x.real_value = 0.0;
@ -504,6 +602,7 @@ void colvar::distance_inv::calc_value()
  x.real_value = std::pow(x.real_value, -1.0/(cvm::real(exponent)));
 }
 void colvar::distance_inv::calc_gradients()
 {
  cvm::real const dxdsum = (-1.0/(cvm::real(exponent))) * std::pow(x.real_value, exponent+1) / cvm::real(group1->size() * group2->size());
@ -515,6 +614,7 @@ void colvar::distance_inv::calc_gradients()
  }
 }
 void colvar::distance_inv::apply_force(colvarvalue const &force)
 {
  if (!group1->noforce)
@ -525,6 +625,9 @@ void colvar::distance_inv::apply_force(colvarvalue const &force)
 }
 simple_scalar_dist_functions(distance_inv)
 colvar::distance_pairs::distance_pairs(std::string const &conf)
  : cvc(conf)
@ -579,11 +682,13 @@ void colvar::distance_pairs::calc_value()
  }
 }
 void colvar::distance_pairs::calc_gradients()
 {
  // will be calculated on the fly in apply_force()
 }
 void colvar::distance_pairs::apply_force(colvarvalue const &force)
 {
  if (b_no_PBC) {
@ -608,6 +713,7 @@ void colvar::distance_pairs::apply_force(colvarvalue const &force)
 }
 colvar::gyration::gyration(std::string const &conf)
  : cvc(conf)
 {
@ -621,6 +727,7 @@ colvar::gyration::gyration(std::string const &conf)
  } else {
    atoms->b_center = true;
    atoms->ref_pos.assign(1, cvm::atom_pos(0.0, 0.0, 0.0));
    atoms->fit_gradients.assign(atoms->size(), cvm::rvector(0.0, 0.0, 0.0));
  }
  x.type(colvarvalue::type_scalar);
@ -681,6 +788,9 @@ void colvar::gyration::apply_force(colvarvalue const &force)
 }
 simple_scalar_dist_functions(gyration)
 colvar::inertia::inertia(std::string const &conf)
  : gyration(conf)
@ -721,6 +831,10 @@ void colvar::inertia::apply_force(colvarvalue const &force)
 }
 simple_scalar_dist_functions(inertia_z)
 colvar::inertia_z::inertia_z(std::string const &conf)
  : inertia(conf)
 {
@ -771,6 +885,10 @@ void colvar::inertia_z::apply_force(colvarvalue const &force)
 }
 simple_scalar_dist_functions(inertia)
 colvar::rmsd::rmsd(std::string const &conf)
  : cvc(conf)
@ -970,6 +1088,8 @@ void colvar::rmsd::calc_Jacobian_derivative()
 }
 simple_scalar_dist_functions(rmsd)
 colvar::eigenvector::eigenvector(std::string const &conf)
@ -1254,6 +1374,10 @@ void colvar::eigenvector::calc_Jacobian_derivative()
 }
 simple_scalar_dist_functions(eigenvector)
 colvar::cartesian::cartesian(std::string const &conf)
  : cvc(conf)
 {
--- a/lib/colvars/colvarcomp_protein.cpp
+++ b/lib/colvars/colvarcomp_protein.cpp
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include <cmath>
 #include "colvarmodule.h"
@ -119,6 +126,7 @@ colvar::alpha_angles::alpha_angles()
  x.type(colvarvalue::type_scalar);
 }
 colvar::alpha_angles::~alpha_angles()
 {
  while (theta.size() != 0) {
@ -131,6 +139,7 @@ colvar::alpha_angles::~alpha_angles()
  }
 }
 void colvar::alpha_angles::calc_value()
 {
  x.real_value = 0.0;
@ -222,6 +231,10 @@ void colvar::alpha_angles::apply_force(colvarvalue const &force)
 }
 simple_scalar_dist_functions(alpha_angles)
 //////////////////////////////////////////////////////////////////////
 // dihedral principal component
 //////////////////////////////////////////////////////////////////////
@ -337,14 +350,22 @@ colvar::dihedPC::dihedPC(std::string const &conf)
  for (size_t i = 0; i < residues.size()-1; i++) {
    // Psi
    theta.push_back(new colvar::dihedral(cvm::atom(r[i  ], "N", sid),
-                                           cvm::atom(r[i  ], "CA", sid),
+                                         cvm::atom(r[i  ], "CA", sid),
-                                           cvm::atom(r[i  ], "C", sid),
+                                         cvm::atom(r[i  ], "C", sid),
-                                           cvm::atom(r[i+1], "N", sid)));
+                                         cvm::atom(r[i+1], "N", sid)));
    atom_groups.push_back(theta.back()->atom_groups[0]);
    atom_groups.push_back(theta.back()->atom_groups[1]);
    atom_groups.push_back(theta.back()->atom_groups[2]);
    atom_groups.push_back(theta.back()->atom_groups[3]);
    // Phi (next res)
    theta.push_back(new colvar::dihedral(cvm::atom(r[i  ], "C", sid),
-                                           cvm::atom(r[i+1], "N", sid),
+                                         cvm::atom(r[i+1], "N", sid),
-                                           cvm::atom(r[i+1], "CA", sid),
+                                         cvm::atom(r[i+1], "CA", sid),
-                                           cvm::atom(r[i+1], "C", sid)));
+                                         cvm::atom(r[i+1], "C", sid)));
    atom_groups.push_back(theta.back()->atom_groups[0]);
    atom_groups.push_back(theta.back()->atom_groups[1]);
    atom_groups.push_back(theta.back()->atom_groups[2]);
    atom_groups.push_back(theta.back()->atom_groups[3]);
  }
  if (cvm::debug())
@ -400,3 +421,6 @@ void colvar::dihedPC::apply_force(colvarvalue const &force)
                           coeffs[2*i+1] * dsindt) * force);
  }
 }
 simple_scalar_dist_functions(dihedPC)
--- a/lib/colvars/colvarcomp_rotations.cpp
+++ b/lib/colvars/colvarcomp_rotations.cpp
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include <cmath>
 #include "colvarmodule.h"
@ -123,6 +130,27 @@ void colvar::orientation::apply_force(colvarvalue const &force)
 }
 cvm::real colvar::orientation::dist2(colvarvalue const &x1,
                                     colvarvalue const &x2) const
 {
  return x1.quaternion_value.dist2(x2);
 }
 colvarvalue colvar::orientation::dist2_lgrad(colvarvalue const &x1,
                                             colvarvalue const &x2) const
 {
  return x1.quaternion_value.dist2_grad(x2);
 }
 colvarvalue colvar::orientation::dist2_rgrad(colvarvalue const &x1,
                                             colvarvalue const &x2) const
 {
  return x2.quaternion_value.dist2_grad(x1);
 }
 colvar::orientation_angle::orientation_angle(std::string const &conf)
  : orientation(conf)
@ -176,6 +204,9 @@ void colvar::orientation_angle::apply_force(colvarvalue const &force)
 }
 simple_scalar_dist_functions(orientation_angle)
 colvar::orientation_proj::orientation_proj(std::string const &conf)
  : orientation(conf)
@ -220,6 +251,9 @@ void colvar::orientation_proj::apply_force(colvarvalue const &force)
 }
 simple_scalar_dist_functions(orientation_proj)
 colvar::tilt::tilt(std::string const &conf)
  : orientation(conf)
@ -278,6 +312,9 @@ void colvar::tilt::apply_force(colvarvalue const &force)
 }
 simple_scalar_dist_functions(tilt)
 colvar::spin_angle::spin_angle(std::string const &conf)
  : orientation(conf)
@ -339,3 +376,46 @@ void colvar::spin_angle::apply_force(colvarvalue const &force)
    atoms->apply_colvar_force(fw);
  }
 }
 cvm::real colvar::spin_angle::dist2(colvarvalue const &x1,
                                    colvarvalue const &x2) const
 {
  cvm::real diff = x1.real_value - x2.real_value;
  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
  return diff * diff;
 }
 colvarvalue colvar::spin_angle::dist2_lgrad(colvarvalue const &x1,
                                            colvarvalue const &x2) const
 {
  cvm::real diff = x1.real_value - x2.real_value;
  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
  return 2.0 * diff;
 }
 colvarvalue colvar::spin_angle::dist2_rgrad(colvarvalue const &x1,
                                            colvarvalue const &x2) const
 {
  cvm::real diff = x1.real_value - x2.real_value;
  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
  return (-2.0) * diff;
 }
 void colvar::spin_angle::wrap(colvarvalue &x) const
 {
  if ((x.real_value - wrap_center) >= 180.0) {
    x.real_value -= 360.0;
    return;
  }
  if ((x.real_value - wrap_center) < -180.0) {
    x.real_value += 360.0;
    return;
  }
  return;
 }
--- a/lib/colvars/colvardeps.cpp
+++ b/lib/colvars/colvardeps.cpp
@ -1,3 +1,13 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include "colvardeps.h"
@ -219,6 +229,9 @@ void colvardeps::init_cvb_requires() {
  f_description(f_cvb_history_dependent, "history-dependent");
  f_description(f_cvb_scalar_variables, "require scalar variables");
  f_req_children(f_cvb_scalar_variables, f_cv_scalar);
  // Initialize feature_states for each instance
  feature_states.reserve(f_cvb_ntot);
  for (i = 0; i < f_cvb_ntot; i++) {
@ -229,6 +242,9 @@ void colvardeps::init_cvb_requires() {
  // some biases are not history-dependent
  feature_states[f_cvb_history_dependent]->available = false;
  // by default, biases should work with vector variables, too
  feature_states[f_cvb_scalar_variables]->available = false;
 }
@ -321,6 +337,7 @@ void colvardeps::init_cv_requires() {
    // The features below are set programmatically
    f_description(f_cv_scripted, "scripted");
    f_description(f_cv_periodic, "periodic");
    f_req_self(f_cv_periodic, f_cv_homogeneous);
    f_description(f_cv_scalar, "scalar");
    f_description(f_cv_linear, "linear");
    f_description(f_cv_homogeneous, "homogeneous");
@ -407,6 +424,11 @@ void colvardeps::init_cvc_requires() {
  // Each cvc specifies what other features are available
  feature_states[f_cvc_active]->available = true;
  feature_states[f_cvc_gradient]->available = true;
  // Features that are implemented by default if their requirements are
  feature_states[f_cvc_one_site_total_force]->available = true;
  // Features That are implemented only for certain simulation engine configurations
  feature_states[f_cvc_scalable_com]->available = (cvm::proxy->scalable_group_coms() == COLVARS_OK);
  feature_states[f_cvc_scalable]->available = feature_states[f_cvc_scalable_com]->available;
 }
--- a/lib/colvars/colvardeps.h
+++ b/lib/colvars/colvardeps.h
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #ifndef COLVARDEPS_H
 #define COLVARDEPS_H
@ -157,6 +164,7 @@ public:
    f_cvb_apply_force, // will apply forces
    f_cvb_get_total_force, // requires total forces
    f_cvb_history_dependent, // depends on simulation history
    f_cvb_scalar_variables, // requires scalar colvars
    f_cvb_ntot
  };
--- a/lib/colvars/colvargrid.cpp
+++ b/lib/colvars/colvargrid.cpp
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include "colvarmodule.h"
 #include "colvarvalue.h"
 #include "colvarparse.h"
@ -73,6 +80,21 @@ cvm::real colvar_grid_scalar::minimum_value() const
  return min;
 }
 cvm::real colvar_grid_scalar::minimum_pos_value() const
 {
  cvm::real minpos = data[0];
  size_t i;
  for (i = 0; i < nt; i++) {
    if(data[i] > 0) {
      minpos = data[i];
      break;
    }
  }
  for (i = 0; i < nt; i++) {
    if (data[i] > 0 && data[i] < minpos) minpos = data[i];
  }
  return minpos;
 }
 cvm::real colvar_grid_scalar::integral() const
 {
--- a/lib/colvars/colvargrid.h
+++ b/lib/colvars/colvargrid.h
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #ifndef COLVARGRID_H
 #define COLVARGRID_H
@ -378,6 +385,13 @@ public:
    return value_to_bin_scalar(actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
  }
  /// \brief Report the bin corresponding to the current value of variable i
  /// and assign first or last bin if out of boundaries
  inline int current_bin_scalar_bound(int const i) const
  {
    return value_to_bin_scalar_bound(actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
  }
  /// \brief Report the bin corresponding to the current value of item iv in variable i
  inline int current_bin_scalar(int const i, int const iv) const
  {
@ -393,6 +407,16 @@ public:
    return (int) std::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
  }
  /// \brief Use the lower boundary and the width to report which bin
  /// the provided value is in and assign first or last bin if out of boundaries
  inline int value_to_bin_scalar_bound(colvarvalue const &value, const int i) const
  {
    int bin_index = std::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
    if (bin_index < 0) bin_index=0;
    if (bin_index >=int(nx[i])) bin_index=int(nx[i])-1;
    return (int) bin_index;
  }
  /// \brief Same as the standard version, but uses another grid definition
  inline int value_to_bin_scalar(colvarvalue const &value,
                                 colvarvalue const &new_offset,
@ -514,6 +538,13 @@ public:
      data[i] *= a;
  }
  /// \brief Assign all zero elements a scalar constant (fast loop)
  inline void remove_zeros(cvm::real const &a)
  {
    for (size_t i = 0; i < nt; i++)
      if(data[i]==0) data[i] = a;
  }
  /// \brief Get the bin indices corresponding to the provided values of
  /// the colvars
@ -537,6 +568,17 @@ public:
    return index;
  }
  /// \brief Get the bin indices corresponding to the provided values of
  /// the colvars and assign first or last bin if out of boundaries
  inline std::vector<int> const get_colvars_index_bound() const
  {
    std::vector<int> index = new_index();
    for (size_t i = 0; i < nd; i++) {
      index[i] = current_bin_scalar_bound(i);
    }
    return index;
  }
  /// \brief Get the minimal distance (in number of bins) from the
  /// boundaries; a negative number is returned if the given point is
  /// off-grid
@ -1071,7 +1113,7 @@ public:
  {
    // write the header
    os << "object 1 class gridpositions counts";
-    int icv;
+    size_t icv;
    for (icv = 0; icv < number_of_colvars(); icv++) {
      os << " " << number_of_points(icv);
    }
@ -1166,45 +1208,49 @@ 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,
+  inline cvm::real log_gradient_finite_diff(const std::vector<int> &ix0,
-                                                   int n = 0)
+                                            int n = 0)
  {
-    cvm::real A0, A1;
+    int A0, A1, A2;
-    std::vector<int> ix;
+    std::vector<int> ix = ix0;
    // factor for mesh width, 2.0 for central finite difference
    // but only 1.0 on edges for non-PBC coordinates
    cvm::real factor;
    if (periodic[n]) {
      factor = 2.;
      ix = ix0;
      ix[n]--; wrap(ix);
      A0 = data[address(ix)];
      ix = ix0;
      ix[n]++; wrap(ix);
      A1 = data[address(ix)];
-    } else {
+      if (A0 * A1 == 0) {
-      factor = 0.;
+        return 0.; // can't handle empty bins
-      ix = ix0;
+      } else {
-      if (ix[n] > 0) { // not left edge
+        return (std::log((cvm::real)A1) - std::log((cvm::real)A0))
-        ix[n]--;
+          / (widths[n] * 2.);
        factor += 1.;
      }
    } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge
      ix[n]--;
      A0 = data[address(ix)];
      ix = ix0;
-      if (ix[n]+1 < nx[n]) { // not right edge
+      ix[n]++;
        ix[n]++;
        factor += 1.;
      }
      A1 = data[address(ix)];
-    }
+      if (A0 * A1 == 0) {
-    if (A0 == 0 || A1 == 0) {
+        return 0.; // can't handle empty bins
-      // can't handle empty bins
+      } else {
-      return 0.;
+        return (std::log((cvm::real)A1) - std::log((cvm::real)A0))
          / (widths[n] * 2.);
      }
    } else {
-      return (std::log((cvm::real)A1) - std::log((cvm::real)A0))
+      // edge: use 2nd order derivative
-        / (widths[n] * factor);
+      int increment = (ix[n] == 0 ? 1 : -1);
      // move right from left edge, or the other way around
      A0 = data[address(ix)];
      ix[n] += increment; A1 = data[address(ix)];
      ix[n] += increment; A2 = data[address(ix)];
      if (A0 * A1 * A2 == 0) {
        return 0.; // can't handle empty bins
      } else {
        return (-1.5 * std::log((cvm::real)A0) + 2. * std::log((cvm::real)A1)
          - 0.5 * std::log((cvm::real)A2)) * increment / widths[n];
      }
    }
  }
 };
@ -1322,6 +1368,9 @@ public:
  /// \brief Return the lowest value
  cvm::real minimum_value() const;
  /// \brief Return the lowest positive value
  cvm::real minimum_pos_value() const;
  /// \brief Calculates the integral of the map (uses widths if they are defined)
  cvm::real integral() const;
--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include <sstream>
 #include <string.h>
@ -296,6 +303,9 @@ int colvarmodule::parse_biases(std::string const &conf)
  /// initialize harmonic restraints
  parse_biases_type<colvarbias_restraint_harmonic>(conf, "harmonic", n_rest_biases);
  /// initialize harmonic walls restraints
  parse_biases_type<colvarbias_restraint_harmonic_walls>(conf, "harmonicWalls", n_rest_biases);
  /// initialize histograms
  parse_biases_type<colvarbias_histogram>(conf, "histogram", n_histo_biases);
@ -562,7 +572,6 @@ int colvarmodule::calc_colvars()
    colvars_smp_items.reserve(colvars.size());
    // set up a vector containing all components
    size_t num_colvar_items = 0;
    cvm::increase_depth();
    for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
@ -576,8 +585,6 @@ int colvarmodule::calc_colvars()
        colvars_smp.push_back(*cvi);
        colvars_smp_items.push_back(icvc);
      }
      num_colvar_items += num_items;
    }
    cvm::decrease_depth();
@ -641,7 +648,7 @@ int colvarmodule::calc_biases()
    for (bi = biases.begin(); bi != biases.end(); bi++) {
      error_code |= (*bi)->update();
      if (cvm::get_error()) {
-        return COLVARS_ERROR;
+        return error_code;
      }
    }
    cvm::decrease_depth();
@ -1007,7 +1014,7 @@ std::istream & colvarmodule::read_restart(std::istream &is)
  for (std::vector<colvarbias *>::iterator bi = biases.begin();
       bi != biases.end();
       bi++) {
-    if (!((*bi)->read_restart(is))) {
+    if (!((*bi)->read_state(is))) {
      cvm::error("Error: in reading restart configuration for bias \""+
                 (*bi)->name+"\".\n",
                 INPUT_ERROR);
@ -1070,15 +1077,15 @@ continue the previous simulation.\n\n");
    cvm::log(cvm::line_marker);
    // update this ahead of time in this special case
-    output_prefix = proxy->output_prefix();
+    output_prefix = proxy->input_prefix();
    cvm::log("All output files will now be saved with the prefix \""+output_prefix+".tmp.*\".\n");
    output_prefix = output_prefix+".tmp";
    write_output_files();
    cvm::log(cvm::line_marker);
    cvm::log("Please review the important warning above. After that, you may rename:\n\
 \""+output_prefix+".tmp.colvars.state\"\n\
 to:\n\
-\""+output_prefix+".colvars.state\"\n");
+\""+ proxy->input_prefix()+".colvars.state\"\n");
    output_prefix = output_prefix+".tmp";
    write_output_files();
    cvm::error("Exiting with error until issue is addressed.\n", FATAL_ERROR);
  }
@ -1120,6 +1127,7 @@ int colvarmodule::write_output_files()
       bi != biases.end();
       bi++) {
    (*bi)->write_output_files();
    (*bi)->write_state_to_replicas();
  }
  cvm::decrease_depth();
@ -1212,20 +1220,30 @@ std::ostream & colvarmodule::write_restart(std::ostream &os)
     << "  version " << std::string(COLVARS_VERSION) << "\n"
     << "}\n\n";
  int error_code = COLVARS_OK;
  cvm::increase_depth();
  for (std::vector<colvar *>::iterator cvi = colvars.begin();
       cvi != colvars.end();
       cvi++) {
    (*cvi)->write_restart(os);
    error_code |= (*cvi)->write_output_files();
  }
  for (std::vector<colvarbias *>::iterator bi = biases.begin();
       bi != biases.end();
       bi++) {
-    (*bi)->write_restart(os);
+    (*bi)->write_state(os);
    error_code |= (*bi)->write_state_to_replicas();
    error_code |= (*bi)->write_output_files();
  }
  cvm::decrease_depth();
  if (error_code != COLVARS_OK) {
    // TODO make this function return an int instead
    os.setstate(std::ios::failbit);
  }
  return os;
 }
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@ -1,10 +1,17 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #ifndef COLVARMODULE_H
 #define COLVARMODULE_H
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2016-10-21"
+#define COLVARS_VERSION "2016-12-27"
 #endif
 #ifndef COLVARS_DEBUG
@ -198,7 +205,7 @@ public:
  }
  /// \brief How many objects are configured yet?
-  inline size_t const size() const
+  inline size_t size() const
  {
    return colvars.size() + biases.size();
  }
--- a/lib/colvars/colvarparse.cpp
+++ b/lib/colvars/colvarparse.cpp
@ -1,3 +1,11 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include <sstream>
--- a/lib/colvars/colvarparse.h
+++ b/lib/colvars/colvarparse.h
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #ifndef COLVARPARSE_H
 #define COLVARPARSE_H
--- a/lib/colvars/colvarproxy.h
+++ b/lib/colvars/colvarproxy.h
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #ifndef COLVARPROXY_H
 #define COLVARPROXY_H
@ -24,11 +31,18 @@ public:
  /// Pointer to the main object
  colvarmodule *colvars;
-  /// Default constructor
+  /// Constructor
-  inline colvarproxy() : script(NULL), b_smp_active(true) {}
+  colvarproxy()
  {
    colvars = NULL;
    b_simulation_running = true;
    b_smp_active = true;
    script = NULL;
  }
-  /// Default destructor
+  /// Destructor
-  virtual ~colvarproxy() {}
+  virtual ~colvarproxy()
  {}
  /// (Re)initialize required member data after construction
  virtual int setup()
@ -104,6 +118,19 @@ public:
    return 0;
  }
 protected:
  /// Whether a simulation is running (and try to prevent irrecovarable errors)
  bool b_simulation_running;
 public:
  /// Whether a simulation is running (and try to prevent irrecovarable errors)
  virtual bool simulation_running() const
  {
    return b_simulation_running;
  }
 protected:
  /// \brief Currently opened output files: by default, these are ofstream objects.
--- a/lib/colvars/colvarscript.cpp
+++ b/lib/colvars/colvarscript.cpp
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include <cstdlib>
 #include <stdlib.h>
 #include <string.h>
--- a/lib/colvars/colvarscript.h
+++ b/lib/colvars/colvarscript.h
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #ifndef COLVARSCRIPT_H
 #define COLVARSCRIPT_H
--- a/lib/colvars/colvartypes.cpp
+++ b/lib/colvars/colvartypes.cpp
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include <stdlib.h>
 #include <string.h>
--- a/lib/colvars/colvartypes.h
+++ b/lib/colvars/colvartypes.h
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #ifndef COLVARTYPES_H
 #define COLVARTYPES_H
--- a/lib/colvars/colvarvalue.cpp
+++ b/lib/colvars/colvarvalue.cpp
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #include <vector>
 #include <sstream>
 #include <iostream>
@ -70,45 +77,30 @@ void colvarvalue::set_elem(int const icv, colvarvalue const &x)
 }
-colvarvalue colvarvalue::inverse() const
+void colvarvalue::set_random()
 {
-  switch (value_type) {
+  size_t ic;
  switch (this->type()) {
  case colvarvalue::type_scalar:
-    return colvarvalue(1.0/real_value);
+    this->real_value = cvm::rand_gaussian();
    break;
  case colvarvalue::type_3vector:
  case colvarvalue::type_unit3vector:
  case colvarvalue::type_unit3vectorderiv:
-    return colvarvalue(cvm::rvector(1.0/rvector_value.x,
+    this->rvector_value.x = cvm::rand_gaussian();
-                                    1.0/rvector_value.y,
+    this->rvector_value.y = cvm::rand_gaussian();
-                                    1.0/rvector_value.z));
+    this->rvector_value.z = cvm::rand_gaussian();
    break;
  case colvarvalue::type_quaternion:
  case colvarvalue::type_quaternionderiv:
-    return colvarvalue(cvm::quaternion(1.0/quaternion_value.q0,
+    this->quaternion_value.q0 = cvm::rand_gaussian();
-                                       1.0/quaternion_value.q1,
+    this->quaternion_value.q1 = cvm::rand_gaussian();
-                                       1.0/quaternion_value.q2,
+    this->quaternion_value.q2 = cvm::rand_gaussian();
-                                       1.0/quaternion_value.q3));
+    this->quaternion_value.q3 = cvm::rand_gaussian();
    break;
  case colvarvalue::type_vector:
-    {
+    for (ic = 0; ic < this->vector1d_value.size(); ic++) {
-      cvm::vector1d<cvm::real> result(vector1d_value);
+      this->vector1d_value[ic] = cvm::rand_gaussian();
      if (elem_types.size() > 0) {
        // if we have information about non-scalar types, use it
        size_t i;
        for (i = 0; i < elem_types.size(); i++) {
          result.sliceassign(elem_indices[i], elem_indices[i]+elem_sizes[i],
                             cvm::vector1d<cvm::real>((this->get_elem(i)).inverse()));
        }
      } else {
        size_t i;
        for (i = 0; i < result.size(); i++) {
          if (result[i] != 0.0) {
            result = 1.0/result[i];
          }
        }
      }
      return colvarvalue(result, type_vector);
    }
    break;
  case colvarvalue::type_notset:
@ -116,7 +108,6 @@ colvarvalue colvarvalue::inverse() const
    undef_op();
    break;
  }
  return colvarvalue();
 }
@ -288,7 +279,7 @@ colvarvalue colvarvalue::dist2_grad(colvarvalue const &x2) const
                                       (-1.0) * sin_t * v2.z +
                                       cos_t/sin_t * (v1.z - cos_t*v2.z)
                                       ),
-                          colvarvalue::type_unit3vector );
+                          colvarvalue::type_unit3vectorderiv );
    }
  case colvarvalue::type_quaternion:
  case colvarvalue::type_quaternionderiv:
--- a/lib/colvars/colvarvalue.h
+++ b/lib/colvars/colvarvalue.h
@ -1,5 +1,12 @@
 // -*- c++ -*-
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
 // Please update all Colvars source files before making any changes.
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 #ifndef COLVARVALUE_H
 #define COLVARVALUE_H
@ -128,30 +135,23 @@ public:
  {}
  /// \brief Copy constructor from rvector base type (Note: this sets
-  /// automatically a type \link type_3vector \endlink , if you want a
+  /// by default a type \link type_3vector \endlink , if you want a
  /// \link type_unit3vector \endlink you must set it explicitly)
-  inline colvarvalue(cvm::rvector const &v)
+  inline colvarvalue(cvm::rvector const &v, Type vti = type_3vector)
    : value_type(type_3vector), rvector_value(v)
  {}
  /// \brief Copy constructor from rvector base type (additional
  /// argument to make possible to choose a \link type_unit3vector
  /// \endlink
  inline colvarvalue(cvm::rvector const &v, Type const &vti)
    : value_type(vti), rvector_value(v)
  {}
  /// \brief Copy constructor from quaternion base type
-  inline colvarvalue(cvm::quaternion const &q)
+  inline colvarvalue(cvm::quaternion const &q, Type vti = type_quaternion)
-    : value_type(type_quaternion), quaternion_value(q)
+    : value_type(vti), quaternion_value(q)
  {}
  /// Copy constructor from vector1d base type
  colvarvalue(cvm::vector1d<cvm::real> const &v, Type vti = type_vector);
  /// Copy constructor from another \link colvarvalue \endlink
  colvarvalue(colvarvalue const &x);
  /// Copy constructor from vector1d base type
  colvarvalue(cvm::vector1d<cvm::real> const &v, Type const &vti);
  /// Set to the null value for the data type currently defined
  void reset();
@ -211,10 +211,6 @@ public:
    return std::sqrt(this->norm2());
  }
  /// \brief Return the value whose scalar product with this value is
  /// 1
  inline colvarvalue inverse() const;
  /// Square distance between this \link colvarvalue \endlink and another
  cvm::real dist2(colvarvalue const &x2) const;
@ -297,6 +293,9 @@ public:
  /// Set elements of the vector from a single colvarvalue
  void set_elem(int const i_begin, int const i_end, colvarvalue const &x);
  /// Make each element a random number in N(0,1)
  void set_random();
  /// Get a single colvarvalue out of elements of the vector
  colvarvalue const get_elem(int const icv) const;
@ -533,7 +532,7 @@ inline colvarvalue::colvarvalue(colvarvalue const &x)
  }
 }
-inline colvarvalue::colvarvalue(cvm::vector1d<cvm::real> const &v, Type const &vti)
+inline colvarvalue::colvarvalue(cvm::vector1d<cvm::real> const &v, Type vti)
 {
  if ((vti != type_vector) && (v.size() != num_dimensions(vti))) {
    cvm::error("Error: trying to initialize a variable of type \""+type_desc(vti)+
@ -617,11 +616,22 @@ inline int colvarvalue::check_types_assign(colvarvalue::Type const &vt1,
  }
  if (vt1 != type_notset) {
-    if (vt1 != vt2) {
+    if (((vt1 == type_unit3vector) &&
-      cvm::error("Trying to assign a colvar value with type \""+
+         (vt2 == type_unit3vectorderiv)) ||
-                 type_desc(vt2)+"\" to one with type \""+
+        ((vt2 == type_unit3vector) &&
-                 type_desc(vt1)+"\".\n");
+         (vt1 == type_unit3vectorderiv)) ||
-      return COLVARS_ERROR;
+        ((vt1 == type_quaternion) &&
         (vt2 == type_quaternionderiv)) ||
        ((vt2 == type_quaternion) &&
         (vt1 == type_quaternionderiv))) {
      return COLVARS_OK;
    } else {
      if (vt1 != vt2) {
        cvm::error("Trying to assign a colvar value with type \""+
                   type_desc(vt2)+"\" to one with type \""+
                   type_desc(vt1)+"\".\n");
        return COLVARS_ERROR;
      }
    }
  }
  return COLVARS_OK;
--- a/src/.gitignore
+++ b/src/.gitignore
@ -24,6 +24,7 @@
 /kokkos.cpp
 /kokkos.h
 /kokkos_type.h
 /kokkos_few.h
 /manifold*.cpp
 /manifold*.h
--- a/src/KOKKOS/Install.sh
+++ b/src/KOKKOS/Install.sh
@ -97,6 +97,8 @@ action fix_reaxc_species_kokkos.cpp fix_reaxc_species.cpp
 action fix_reaxc_species_kokkos.h fix_reaxc_species.h
 action fix_setforce_kokkos.cpp
 action fix_setforce_kokkos.h
 action fix_momentum_kokkos.cpp
 action fix_momentum_kokkos.h
 action fix_wall_reflect_kokkos.cpp
 action fix_wall_reflect_kokkos.h
 action fix_dpd_energy_kokkos.cpp fix_dpd_energy.cpp
@ -108,6 +110,7 @@ action improper_harmonic_kokkos.h improper_harmonic.h
 action kokkos.cpp
 action kokkos.h
 action kokkos_type.h
 action kokkos_few.h
 action memory_kokkos.h
 action modify_kokkos.cpp
 action modify_kokkos.h
--- a/src/KOKKOS/atom_vec_angle_kokkos.cpp
+++ b/src/KOKKOS/atom_vec_angle_kokkos.cpp
@ -1593,7 +1593,7 @@ int AtomVecAngleKokkos::unpack_restart(double *buf)
  double **extra = atom->extra;
  if (atom->nextra_store) {
-    int size = static_cast<int> (ubuf(buf[m++]).i) - m;
+    int size = static_cast<int> (buf[0]) - m;
    for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++];
  }
--- a/src/KOKKOS/atom_vec_atomic_kokkos.cpp
+++ b/src/KOKKOS/atom_vec_atomic_kokkos.cpp
@ -1232,7 +1232,7 @@ int AtomVecAtomicKokkos::unpack_restart(double *buf)
  double **extra = atom->extra;
  if (atom->nextra_store) {
-    int size = static_cast<int> (ubuf(buf[m++]).i) - m;
+    int size = static_cast<int> (buf[0]) - m;
    for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++];
  }
--- a/src/KOKKOS/atom_vec_bond_kokkos.cpp
+++ b/src/KOKKOS/atom_vec_bond_kokkos.cpp
@ -1460,7 +1460,7 @@ int AtomVecBondKokkos::unpack_restart(double *buf)
  double **extra = atom->extra;
  if (atom->nextra_store) {
-    int size = static_cast<int> (ubuf(buf[m++]).i) - m;
+    int size = static_cast<int> (buf[0]) - m;
    for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++];
  }
--- a/src/KOKKOS/atom_vec_charge_kokkos.cpp
+++ b/src/KOKKOS/atom_vec_charge_kokkos.cpp
@ -1314,7 +1314,7 @@ int AtomVecChargeKokkos::unpack_restart(double *buf)
  double **extra = atom->extra;
  if (atom->nextra_store) {
-    int size = static_cast<int> (ubuf(buf[m++]).i) - m;
+    int size = static_cast<int> (buf[0]) - m;
    for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++];
  }
--- a/src/KOKKOS/atom_vec_full_kokkos.cpp
+++ b/src/KOKKOS/atom_vec_full_kokkos.cpp
@ -1936,7 +1936,7 @@ int AtomVecFullKokkos::unpack_restart(double *buf)
  double **extra = atom->extra;
  if (atom->nextra_store) {
-    int size = static_cast<int> (ubuf(buf[m++]).i) - m;
+    int size = static_cast<int> (buf[0]) - m;
    for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++];
  }
--- a/src/KOKKOS/atom_vec_molecular_kokkos.cpp
+++ b/src/KOKKOS/atom_vec_molecular_kokkos.cpp
@ -1850,7 +1850,7 @@ int AtomVecMolecularKokkos::unpack_restart(double *buf)
  double **extra = atom->extra;
  if (atom->nextra_store) {
-    int size = static_cast<int> (ubuf(buf[m++]).i) - m;
+    int size = static_cast<int> (buf[0]) - m;
    for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++];
  }
--- a/src/KOKKOS/comm_tiled_kokkos.cpp
+++ b/src/KOKKOS/comm_tiled_kokkos.cpp
@ -262,5 +262,5 @@ void CommTiledKokkos::forward_comm_array(int nsize, double **array)
 int CommTiledKokkos::exchange_variable(int n, double *inbuf, double *&outbuf)
 {
-  CommTiled::exchange_variable(n,inbuf,outbuf);
+  return CommTiled::exchange_variable(n,inbuf,outbuf);
 }
--- a/src/KOKKOS/compute_temp_kokkos.h
+++ b/src/KOKKOS/compute_temp_kokkos.h
@ -45,14 +45,13 @@ namespace LAMMPS_NS {
    }
    KOKKOS_INLINE_FUNCTION
-    volatile s_CTEMP& operator+=(const volatile s_CTEMP &rhs) volatile {
+    void operator+=(const volatile s_CTEMP &rhs) volatile {
      t0 += rhs.t0;
      t1 += rhs.t1;
      t2 += rhs.t2;
      t3 += rhs.t3;
      t4 += rhs.t4;
      t5 += rhs.t5;
      return *this;
    }
  };
  typedef s_CTEMP CTEMP;
--- a/src/KOKKOS/domain_kokkos.h
+++ b/src/KOKKOS/domain_kokkos.h
@ -16,6 +16,7 @@
 #include "domain.h"
 #include "kokkos_type.h"
 #include "kokkos_few.h"
 namespace LAMMPS_NS {
@ -35,6 +36,10 @@ class DomainKokkos : public Domain {
  void image_flip(int, int, int);
  void x2lamda(int);
  void lamda2x(int);
  // these lines bring in the x2lamda signatures from Domain
  // that are not overloaded here
  using Domain::x2lamda;
  using Domain::lamda2x;
  int closest_image(const int, int) const;
@ -50,6 +55,10 @@ class DomainKokkos : public Domain {
  KOKKOS_INLINE_FUNCTION
  void operator()(TagDomain_x2lamda, const int&) const;
  static KOKKOS_INLINE_FUNCTION
  Few<double,3> unmap(Few<double,3> prd, Few<double,6> h, int triclinic,
      Few<double,3> x, imageint image);
 private:
  double lo[3],hi[3],period[3];
  int n_flip, m_flip, p_flip;
@ -57,6 +66,26 @@ class DomainKokkos : public Domain {
  ArrayTypes<LMPDeviceType>::t_imageint_1d image;
 };
 KOKKOS_INLINE_FUNCTION
 Few<double,3> DomainKokkos::unmap(Few<double,3> prd, Few<double,6> h,
    int triclinic, Few<double,3> x, imageint image)
 {
  int xbox = (image & IMGMASK) - IMGMAX;
  int ybox = (image >> IMGBITS & IMGMASK) - IMGMAX;
  int zbox = (image >> IMG2BITS) - IMGMAX;
  Few<double,3> y;
  if (triclinic == 0) {
    y[0] = x[0] + xbox*prd[0];
    y[1] = x[1] + ybox*prd[1];
    y[2] = x[2] + zbox*prd[2];
  } else {
    y[0] = x[0] + h[0]*xbox + h[5]*ybox + h[4]*zbox;
    y[1] = x[1] + h[1]*ybox + h[3]*zbox;
    y[2] = x[2] + h[2]*zbox;
  }
  return y;
 }
 }
 #endif
--- a/src/KOKKOS/fix_momentum_kokkos.cpp
+++ b/src/KOKKOS/fix_momentum_kokkos.cpp
@ -0,0 +1,204 @@
 /* ----------------------------------------------------------------------
   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 <stdlib.h>
 #include <string.h>
 #include "fix_momentum_kokkos.h"
 #include "atom_kokkos.h"
 #include "atom_masks.h"
 #include "domain.h"
 #include "domain_kokkos.h"
 #include "group.h"
 #include "error.h"
 #include "force.h"
 #include "kokkos_few.h"
 using namespace LAMMPS_NS;
 using namespace FixConst;
 /* ----------------------------------------------------------------------
   Contributing author: Dan Ibanez (SNL)
 ------------------------------------------------------------------------- */
 /* ---------------------------------------------------------------------- */
 template<class DeviceType>
 FixMomentumKokkos<DeviceType>::FixMomentumKokkos(LAMMPS *lmp, int narg, char **arg) :
  FixMomentum(lmp, narg, arg)
 {
  kokkosable = 1;
  atomKK = (AtomKokkos *) atom;
  execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
  datamask_read = EMPTY_MASK;
  datamask_modify = EMPTY_MASK;
 }
 /* ---------------------------------------------------------------------- */
 template<class DeviceType>
 static double get_kinetic_energy(
    AtomKokkos* atomKK,
    MPI_Comm world,
    int groupbit,
    int nlocal,
    typename ArrayTypes<DeviceType>::t_v_array_randomread v,
    typename ArrayTypes<DeviceType>::t_int_1d_randomread mask)
 {
  using AT = ArrayTypes<DeviceType>;
  auto execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
  double ke=0.0;
  if (atomKK->rmass) {
    atomKK->sync(execution_space, RMASS_MASK);
    typename AT::t_float_1d_randomread rmass = atomKK->k_rmass.view<DeviceType>();
    Kokkos::parallel_reduce(nlocal, LAMMPS_LAMBDA(int i, double& update) {
      if (mask(i) & groupbit)
        update += rmass(i) *
          (v(i,0)*v(i,0) + v(i,1)*v(i,1) + v(i,2)*v(i,2));
    }, ke);
  } else {
    // D.I. : why is there no MASS_MASK ?
    atomKK->sync(execution_space, TYPE_MASK);
    typename AT::t_int_1d_randomread type = atomKK->k_type.view<DeviceType>();
    typename AT::t_float_1d_randomread mass = atomKK->k_mass.view<DeviceType>();
    Kokkos::parallel_reduce(nlocal, LAMMPS_LAMBDA(int i, double& update) {
      if (mask(i) & groupbit)
        update += mass(type(i)) *
          (v(i,0)*v(i,0) + v(i,1)*v(i,1) + v(i,2)*v(i,2));
    }, ke);
  }
  double ke_total;
  MPI_Allreduce(&ke,&ke_total,1,MPI_DOUBLE,MPI_SUM,world);
  return ke_total;
 }
 template<class DeviceType>
 void FixMomentumKokkos<DeviceType>::end_of_step()
 {
  atomKK->sync(execution_space, V_MASK | MASK_MASK);
  typename AT::t_v_array v = atomKK->k_v.view<DeviceType>();
  typename AT::t_int_1d_randomread mask = atomKK->k_mask.view<DeviceType>();
  const int nlocal = atom->nlocal;
  double ekin_old,ekin_new;
  ekin_old = ekin_new = 0.0;
  if (dynamic)
    masstotal = group->mass(igroup); // change once Group is ported to Kokkos
  // do nothing if group is empty, i.e. mass is zero;
  if (masstotal == 0.0) return;
  // compute kinetic energy before momentum removal, if needed
  if (rescale) {
    ekin_old = get_kinetic_energy<DeviceType>(atomKK, world, groupbit, nlocal, v, mask);
  }
  auto groupbit2 = groupbit;
  if (linear) {
    /* this is needed because Group is not Kokkos-aware ! */
    atomKK->sync(ExecutionSpaceFromDevice<LMPHostType>::space,
        V_MASK | MASK_MASK | TYPE_MASK | RMASS_MASK);
    Few<double, 3> vcm;
    group->vcm(igroup,masstotal,&vcm[0]);
    // adjust velocities by vcm to zero linear momentum
    // only adjust a component if flag is set
    auto xflag2 = xflag;
    auto yflag2 = yflag;
    auto zflag2 = zflag;
    Kokkos::parallel_for(nlocal, LAMMPS_LAMBDA(int i) {
      if (mask(i) & groupbit2) {
        if (xflag2) v(i,0) -= vcm[0];
        if (yflag2) v(i,1) -= vcm[1];
        if (zflag2) v(i,2) -= vcm[2];
      }
    });
    atomKK->modified(execution_space, V_MASK);
  }
  if (angular) {
    Few<double, 3> xcm, angmom, omega;
    double inertia[3][3];
    /* syncs for each Kokkos-unaware Group method */
    atomKK->sync(ExecutionSpaceFromDevice<LMPHostType>::space,
        X_MASK | MASK_MASK | TYPE_MASK | IMAGE_MASK | RMASS_MASK);
    group->xcm(igroup,masstotal,&xcm[0]);
    atomKK->sync(ExecutionSpaceFromDevice<LMPHostType>::space,
        X_MASK | V_MASK | MASK_MASK | TYPE_MASK | IMAGE_MASK | RMASS_MASK);
    group->angmom(igroup,&xcm[0],&angmom[0]);
    atomKK->sync(ExecutionSpaceFromDevice<LMPHostType>::space,
        X_MASK | MASK_MASK | TYPE_MASK | IMAGE_MASK | RMASS_MASK);
    group->inertia(igroup,&xcm[0],inertia);
    group->omega(&angmom[0],inertia,&omega[0]);
    // adjust velocities to zero omega
    // vnew_i = v_i - w x r_i
    // must use unwrapped coords to compute r_i correctly
    atomKK->sync(execution_space, X_MASK | IMAGE_MASK);
    typename AT::t_x_array_randomread x = atomKK->k_x.view<DeviceType>();
    typename AT::t_imageint_1d_randomread image = atomKK->k_image.view<DeviceType>();
    int nlocal = atom->nlocal;
    auto prd = Few<double,3>(domain->prd);
    auto h = Few<double,6>(domain->h);
    auto triclinic = domain->triclinic;
    Kokkos::parallel_for(nlocal, LAMMPS_LAMBDA(int i) {
      if (mask[i] & groupbit2) {
        Few<double,3> x_i;
        x_i[0] = x(i,0);
        x_i[1] = x(i,1);
        x_i[2] = x(i,2);
        auto unwrap = DomainKokkos::unmap(prd,h,triclinic,x_i,image(i));
        auto dx = unwrap[0] - xcm[0];
        auto dy = unwrap[1] - xcm[1];
        auto dz = unwrap[2] - xcm[2];
        v(i,0) -= omega[1]*dz - omega[2]*dy;
        v(i,1) -= omega[2]*dx - omega[0]*dz;
        v(i,2) -= omega[0]*dy - omega[1]*dx;
      }
    });
    atomKK->modified(execution_space, V_MASK);
  }
  // compute kinetic energy after momentum removal, if needed
  if (rescale) {
    ekin_new = get_kinetic_energy<DeviceType>(atomKK, world, groupbit, nlocal, v, mask);
    double factor = 1.0;
    if (ekin_new != 0.0) factor = sqrt(ekin_old/ekin_new);
    Kokkos::parallel_for(nlocal, LAMMPS_LAMBDA(int i) {
      if (mask(i) & groupbit2) {
        v(i,0) *= factor;
        v(i,1) *= factor;
        v(i,2) *= factor;
      }
    });
    atomKK->modified(execution_space, V_MASK);
  }
 }
 namespace LAMMPS_NS {
 template class FixMomentumKokkos<LMPDeviceType>;
 #ifdef KOKKOS_HAVE_CUDA
 template class FixMomentumKokkos<LMPHostType>;
 #endif
 }
--- a/src/KOKKOS/fix_momentum_kokkos.h
+++ b/src/KOKKOS/fix_momentum_kokkos.h
@ -0,0 +1,46 @@
 /* -*- 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.
 ------------------------------------------------------------------------- */
 #ifdef FIX_CLASS
 FixStyle(momentum/kk,FixMomentumKokkos<LMPDeviceType>)
 FixStyle(momentum/kk/device,FixMomentumKokkos<LMPDeviceType>)
 FixStyle(momentum/kk/host,FixMomentumKokkos<LMPHostType>)
 #else
 #ifndef LMP_FIX_MOMENTUM_KOKKOS_H
 #define LMP_FIX_MOMENTUM_KOKKOS_H
 #include "fix_momentum.h"
 #include "kokkos_type.h"
 namespace LAMMPS_NS {
 template<class DeviceType>
 class FixMomentumKokkos : public FixMomentum {
 public:
  typedef ArrayTypes<DeviceType> AT;
  FixMomentumKokkos(class LAMMPS *, int, char **);
  void end_of_step();
 };
 }
 #endif
 #endif
 /* ERROR/WARNING messages:
 */
--- a/src/KOKKOS/fix_setforce_kokkos.h
+++ b/src/KOKKOS/fix_setforce_kokkos.h
@ -43,11 +43,10 @@ struct s_double_3 {
  }
  KOKKOS_INLINE_FUNCTION
-  volatile s_double_3& operator+=(const volatile s_double_3 &rhs) volatile {
+  void operator+=(const volatile s_double_3 &rhs) volatile {
    d0 += rhs.d0;
    d1 += rhs.d1;
    d2 += rhs.d2;
    return *this;
  }
 };
 typedef s_double_3 double_3;
--- a/src/KOKKOS/kokkos_few.h
+++ b/src/KOKKOS/kokkos_few.h
@ -0,0 +1,68 @@
 #ifndef KOKKOS_FEW_H
 #define KOKKOS_FEW_H
 #include <Kokkos_Core.hpp>
 template <typename T, std::size_t n>
 class Few {
  alignas(T) char array_[n * sizeof(T)];
 public:
  enum { size = n };
  Few(std::initializer_list<T> l) {
    std::size_t i = 0;
    for (auto it = l.begin(); it != l.end(); ++it) {
      new (data() + (i++)) T(*it);
    }
  }
  KOKKOS_INLINE_FUNCTION Few(T const a[]) {
    for (std::size_t i = 0; i < n; ++i) new (data() + i) T(a[i]);
  }
  KOKKOS_INLINE_FUNCTION Few() {
    for (std::size_t i = 0; i < n; ++i) new (data() + i) T();
  }
  KOKKOS_INLINE_FUNCTION ~Few() {
    for (std::size_t i = 0; i < n; ++i) (data()[i]).~T();
  }
  KOKKOS_INLINE_FUNCTION Few(Few<T, n> const& rhs) {
    for (std::size_t i = 0; i < n; ++i) new (data() + i) T(rhs[i]);
  }
  KOKKOS_INLINE_FUNCTION Few(Few<T, n> const volatile& rhs) {
    for (std::size_t i = 0; i < n; ++i) new (data() + i) T(rhs[i]);
  }
  KOKKOS_INLINE_FUNCTION void operator=(Few<T, n> const& rhs) volatile {
    for (std::size_t i = 0; i < n; ++i) data()[i] = rhs[i];
  }
  KOKKOS_INLINE_FUNCTION void operator=(Few<T, n> const& rhs) {
    for (std::size_t i = 0; i < n; ++i) data()[i] = rhs[i];
  }
  KOKKOS_INLINE_FUNCTION void operator=(Few<T, n> const volatile& rhs) {
    for (std::size_t i = 0; i < n; ++i) data()[i] = rhs[i];
  }
  KOKKOS_INLINE_FUNCTION T* data() {
    return reinterpret_cast<T*>(array_);
  }
  KOKKOS_INLINE_FUNCTION T const* data() const {
    return reinterpret_cast<T const*>(array_);
  }
  KOKKOS_INLINE_FUNCTION T volatile* data() volatile {
    return reinterpret_cast<T volatile*>(array_);
  }
  KOKKOS_INLINE_FUNCTION T const volatile* data() const volatile {
    return reinterpret_cast<T const volatile*>(array_);
  }
  KOKKOS_INLINE_FUNCTION T& operator[](std::size_t i) {
    return data()[i];
  }
  KOKKOS_INLINE_FUNCTION T const& operator[](std::size_t i) const {
    return data()[i];
  }
  KOKKOS_INLINE_FUNCTION T volatile& operator[](std::size_t i) volatile {
    return data()[i];
  }
  KOKKOS_INLINE_FUNCTION T const volatile& operator[](std::size_t i) const volatile {
    return data()[i];
  }
 };
 #endif
--- a/src/KOKKOS/modify_kokkos.cpp
+++ b/src/KOKKOS/modify_kokkos.cpp
@ -81,7 +81,7 @@ void ModifyKokkos::setup_pre_exchange()
      atomKK->sync(fix[list_min_pre_exchange[i]]->execution_space,
                   fix[list_min_pre_exchange[i]]->datamask_read);
      if (!fix[list_min_pre_exchange[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
-      fix[list_min_pre_exchange[i]]->min_setup_pre_exchange();
+      fix[list_min_pre_exchange[i]]->setup_pre_exchange();
      lmp->kokkos->auto_sync = 0;
      atomKK->modified(fix[list_min_pre_exchange[i]]->execution_space,
                       fix[list_min_pre_exchange[i]]->datamask_modify);
@ -110,7 +110,7 @@ void ModifyKokkos::setup_pre_neighbor()
      atomKK->sync(fix[list_min_pre_neighbor[i]]->execution_space,
                   fix[list_min_pre_neighbor[i]]->datamask_read);
      if (!fix[list_min_pre_neighbor[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
-      fix[list_min_pre_neighbor[i]]->min_setup_pre_neighbor();
+      fix[list_min_pre_neighbor[i]]->setup_pre_neighbor();
      lmp->kokkos->auto_sync = 0;
      atomKK->modified(fix[list_min_pre_neighbor[i]]->execution_space,
                       fix[list_min_pre_neighbor[i]]->datamask_modify);
@ -139,13 +139,42 @@ void ModifyKokkos::setup_pre_force(int vflag)
      atomKK->sync(fix[list_min_pre_force[i]]->execution_space,
                   fix[list_min_pre_force[i]]->datamask_read);
      if (!fix[list_min_pre_force[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
-      fix[list_min_pre_force[i]]->min_setup_pre_force(vflag);
+      fix[list_min_pre_force[i]]->setup_pre_force(vflag);
      lmp->kokkos->auto_sync = 0;
      atomKK->modified(fix[list_min_pre_force[i]]->execution_space,
                       fix[list_min_pre_force[i]]->datamask_modify);
    }
 }
 /* ----------------------------------------------------------------------
   setup pre_reverse call, only for fixes that define pre_reverse
   called from Verlet, RESPA, Min
 ------------------------------------------------------------------------- */
 void ModifyKokkos::setup_pre_reverse(int eflag, int vflag)
 {
  if (update->whichflag == 1)
    for (int i = 0; i < n_pre_reverse; i++) {
      atomKK->sync(fix[list_pre_reverse[i]]->execution_space,
                   fix[list_pre_reverse[i]]->datamask_read);
      if (!fix[list_pre_reverse[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
      fix[list_pre_reverse[i]]->setup_pre_reverse(eflag,vflag);
      lmp->kokkos->auto_sync = 0;
      atomKK->modified(fix[list_pre_reverse[i]]->execution_space,
                       fix[list_pre_reverse[i]]->datamask_modify);
    }
  else if (update->whichflag == 2)
    for (int i = 0; i < n_min_pre_reverse; i++) {
      atomKK->sync(fix[list_min_pre_reverse[i]]->execution_space,
                   fix[list_min_pre_reverse[i]]->datamask_read);
      if (!fix[list_min_pre_reverse[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
      fix[list_min_pre_reverse[i]]->setup_pre_reverse(eflag,vflag);
      lmp->kokkos->auto_sync = 0;
      atomKK->modified(fix[list_min_pre_reverse[i]]->execution_space,
                       fix[list_min_pre_reverse[i]]->datamask_modify);
    }
 }
 /* ----------------------------------------------------------------------
   1st half of integrate call, only for relevant fixes
 ------------------------------------------------------------------------- */
@ -231,6 +260,23 @@ void ModifyKokkos::pre_force(int vflag)
  }
 }
 /* ----------------------------------------------------------------------
   pre_reverse call, only for relevant fixes
 ------------------------------------------------------------------------- */
 void ModifyKokkos::pre_reverse(int eflag, int vflag)
 {
  for (int i = 0; i < n_pre_reverse; i++) {
    atomKK->sync(fix[list_pre_reverse[i]]->execution_space,
                 fix[list_pre_reverse[i]]->datamask_read);
    if (!fix[list_pre_reverse[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
    fix[list_pre_reverse[i]]->pre_reverse(eflag,vflag);
    lmp->kokkos->auto_sync = 0;
    atomKK->modified(fix[list_pre_reverse[i]]->execution_space,
                     fix[list_pre_reverse[i]]->datamask_modify);
  }
 }
 /* ----------------------------------------------------------------------
   post_force call, only for relevant fixes
 ------------------------------------------------------------------------- */
@ -476,6 +522,23 @@ void ModifyKokkos::min_pre_force(int vflag)
  }
 }
 /* ----------------------------------------------------------------------
   minimizer pre-reverse call, only for relevant fixes
 ------------------------------------------------------------------------- */
 void ModifyKokkos::min_pre_reverse(int eflag, int vflag)
 {
  for (int i = 0; i < n_min_pre_reverse; i++) {
    atomKK->sync(fix[list_min_pre_reverse[i]]->execution_space,
                 fix[list_min_pre_reverse[i]]->datamask_read);
    if (!fix[list_min_pre_reverse[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
    fix[list_min_pre_reverse[i]]->min_pre_reverse(eflag,vflag);
    lmp->kokkos->auto_sync = 0;
    atomKK->modified(fix[list_min_pre_reverse[i]]->execution_space,
                     fix[list_min_pre_reverse[i]]->datamask_modify);
  }
 }
 /* ----------------------------------------------------------------------
   minimizer force adjustment call, only for relevant fixes
 ------------------------------------------------------------------------- */
@ -658,4 +721,3 @@ int ModifyKokkos::min_reset_ref()
  }
  return itmpall;
 }
--- a/src/KOKKOS/modify_kokkos.h
+++ b/src/KOKKOS/modify_kokkos.h
@ -26,12 +26,14 @@ class ModifyKokkos : public Modify {
  void setup_pre_exchange();
  void setup_pre_neighbor();
  void setup_pre_force(int);
  void setup_pre_reverse(int, int);
  void initial_integrate(int);
  void post_integrate();
  void pre_decide();
  void pre_exchange();
  void pre_neighbor();
  void pre_force(int);
  void pre_reverse(int,int);
  void post_force(int);
  void final_integrate();
  void end_of_step();
@ -48,6 +50,7 @@ class ModifyKokkos : public Modify {
  void min_pre_exchange();
  void min_pre_neighbor();
  void min_pre_force(int);
  void min_pre_reverse(int,int);
  void min_post_force(int);
  double min_energy(double *);
--- a/src/KOKKOS/pair_buck_coul_cut_kokkos.h
+++ b/src/KOKKOS/pair_buck_coul_cut_kokkos.h
@ -75,7 +75,7 @@ class PairBuckCoulCutKokkos : public PairBuckCoulCut {
  Kokkos::DualView<params_buck_coul**,Kokkos::LayoutRight,DeviceType> k_params;
  typename Kokkos::DualView<params_buck_coul**,
    Kokkos::LayoutRight,DeviceType>::t_dev_const_um params;
-  // hardwired to space for 15 atom types
+  // hardwired to space for 12 atom types
  params_buck_coul m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
--- a/src/KOKKOS/pair_buck_coul_long_kokkos.h
+++ b/src/KOKKOS/pair_buck_coul_long_kokkos.h
@ -76,7 +76,7 @@ class PairBuckCoulLongKokkos : public PairBuckCoulLong {
  Kokkos::DualView<params_buck_coul**,Kokkos::LayoutRight,DeviceType> k_params;
  typename Kokkos::DualView<params_buck_coul**,
    Kokkos::LayoutRight,DeviceType>::t_dev_const_um params;
-  // hardwired to space for 15 atom types
+  // hardwired to space for 12 atom types
  params_buck_coul m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
--- a/src/KOKKOS/pair_buck_kokkos.h
+++ b/src/KOKKOS/pair_buck_kokkos.h
@ -68,7 +68,7 @@ class PairBuckKokkos : public PairBuck {
  Kokkos::DualView<params_buck**,Kokkos::LayoutRight,DeviceType> k_params;
  typename Kokkos::DualView<params_buck**,Kokkos::LayoutRight,DeviceType>::t_dev_const_um params;
-  params_buck m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];  // hardwired to space for 15 atom types
+  params_buck m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];  // hardwired to space for 12 atom types
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  typename ArrayTypes<DeviceType>::t_x_array_randomread x;
  typename ArrayTypes<DeviceType>::t_x_array c_x;
--- a/src/KOKKOS/pair_coul_cut_kokkos.h
+++ b/src/KOKKOS/pair_coul_cut_kokkos.h
@ -81,7 +81,7 @@ class PairCoulCutKokkos : public PairCoulCut {
  Kokkos::DualView<params_coul**,Kokkos::LayoutRight,DeviceType> k_params;
  typename Kokkos::DualView<params_coul**,
    Kokkos::LayoutRight,DeviceType>::t_dev_const_um params;
-  // hardwired to space for 15 atom types
+  // hardwired to space for 12 atom types
  params_coul m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
--- a/src/KOKKOS/pair_coul_debye_kokkos.h
+++ b/src/KOKKOS/pair_coul_debye_kokkos.h
@ -79,7 +79,7 @@ class PairCoulDebyeKokkos : public PairCoulDebye {
  Kokkos::DualView<params_coul**,Kokkos::LayoutRight,DeviceType> k_params;
  typename Kokkos::DualView<params_coul**,
    Kokkos::LayoutRight,DeviceType>::t_dev_const_um params;
-  // hardwired to space for 15 atom types
+  // hardwired to space for 12 atom types
  params_coul m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
--- a/src/KOKKOS/pair_coul_long_kokkos.h
+++ b/src/KOKKOS/pair_coul_long_kokkos.h
@ -80,7 +80,7 @@ class PairCoulLongKokkos : public PairCoulLong {
  Kokkos::DualView<params_coul**,Kokkos::LayoutRight,DeviceType> k_params;
  typename Kokkos::DualView<params_coul**,
    Kokkos::LayoutRight,DeviceType>::t_dev_const_um params;
-  // hardwired to space for 15 atom types
+  // hardwired to space for 12 atom types
  params_coul m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
--- a/src/KOKKOS/pair_lj_charmm_coul_charmm_implicit_kokkos.cpp
+++ b/src/KOKKOS/pair_lj_charmm_coul_charmm_implicit_kokkos.cpp
@ -490,14 +490,11 @@ double PairLJCharmmCoulCharmmImplicitKokkos<DeviceType>::init_one(int i, int j)
    m_cut_coulsq[j][i] = m_cut_coulsq[i][j] = cut_coulsqm;
  }
-  k_cutsq.h_view(i,j) = cutone*cutone;
+  k_cutsq.h_view(i,j) = k_cutsq.h_view(j,i) = cutone*cutone;
  k_cutsq.h_view(j,i) = k_cutsq.h_view(i,j);
  k_cutsq.template modify<LMPHostType>();
-  k_cut_ljsq.h_view(i,j) = cut_ljsqm;
+  k_cut_ljsq.h_view(i,j) = k_cut_ljsq.h_view(j,i) = cut_ljsqm;
  k_cut_ljsq.h_view(j,i) = k_cut_ljsq.h_view(i,j);
  k_cut_ljsq.template modify<LMPHostType>();
-  k_cut_coulsq.h_view(i,j) = cut_coulsqm;
+  k_cut_coulsq.h_view(i,j) = k_cut_coulsq.h_view(j,i) = cut_coulsqm;
  k_cut_coulsq.h_view(j,i) = k_cut_coulsq.h_view(i,j);
  k_cut_coulsq.template modify<LMPHostType>();
  k_params.template modify<LMPHostType>();
--- a/src/KOKKOS/pair_lj_charmm_coul_charmm_implicit_kokkos.h
+++ b/src/KOKKOS/pair_lj_charmm_coul_charmm_implicit_kokkos.h
@ -76,7 +76,7 @@ class PairLJCharmmCoulCharmmImplicitKokkos : public PairLJCharmmCoulCharmmImplic
  Kokkos::DualView<params_lj_coul**,Kokkos::LayoutRight,DeviceType> k_params;
  typename Kokkos::DualView<params_lj_coul**,
    Kokkos::LayoutRight,DeviceType>::t_dev_const_um params;
-  // hardwired to space for 15 atom types
+  // hardwired to space for 12 atom types
  params_lj_coul m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
--- a/src/KOKKOS/pair_lj_charmm_coul_charmm_kokkos.cpp
+++ b/src/KOKKOS/pair_lj_charmm_coul_charmm_kokkos.cpp
@ -491,14 +491,11 @@ double PairLJCharmmCoulCharmmKokkos<DeviceType>::init_one(int i, int j)
    m_cut_coulsq[j][i] = m_cut_coulsq[i][j] = cut_coulsqm;
  }
-  k_cutsq.h_view(i,j) = cutone*cutone;
+  k_cutsq.h_view(i,j) = k_cutsq.h_view(j,i) = cutone*cutone;
  k_cutsq.h_view(j,i) = k_cutsq.h_view(i,j);
  k_cutsq.template modify<LMPHostType>();
-  k_cut_ljsq.h_view(i,j) = cut_ljsqm;
+  k_cut_ljsq.h_view(i,j) = k_cut_ljsq.h_view(j,i) = cut_ljsqm;
  k_cut_ljsq.h_view(j,i) = k_cut_ljsq.h_view(i,j);
  k_cut_ljsq.template modify<LMPHostType>();
-  k_cut_coulsq.h_view(i,j) = cut_coulsqm;
+  k_cut_coulsq.h_view(i,j) = k_cut_coulsq.h_view(j,i) = cut_coulsqm;
  k_cut_coulsq.h_view(j,i) = k_cut_coulsq.h_view(i,j);
  k_cut_coulsq.template modify<LMPHostType>();
  k_params.template modify<LMPHostType>();
--- a/src/KOKKOS/pair_lj_charmm_coul_charmm_kokkos.h
+++ b/src/KOKKOS/pair_lj_charmm_coul_charmm_kokkos.h
@ -76,7 +76,7 @@ class PairLJCharmmCoulCharmmKokkos : public PairLJCharmmCoulCharmm {
  Kokkos::DualView<params_lj_coul**,Kokkos::LayoutRight,DeviceType> k_params;
  typename Kokkos::DualView<params_lj_coul**,
    Kokkos::LayoutRight,DeviceType>::t_dev_const_um params;
-  // hardwired to space for 15 atom types
+  // hardwired to space for 12 atom types
  params_lj_coul m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
--- a/src/KOKKOS/pair_lj_charmm_coul_long_kokkos.cpp
+++ b/src/KOKKOS/pair_lj_charmm_coul_long_kokkos.cpp
@ -521,14 +521,11 @@ double PairLJCharmmCoulLongKokkos<DeviceType>::init_one(int i, int j)
    m_cut_coulsq[j][i] = m_cut_coulsq[i][j] = cut_coulsqm;
  }
-  k_cutsq.h_view(i,j) = cutone*cutone;
+  k_cutsq.h_view(i,j) = k_cutsq.h_view(j,i) = cutone*cutone;
  k_cutsq.h_view(j,i) = k_cutsq.h_view(i,j);
  k_cutsq.template modify<LMPHostType>();
-  k_cut_ljsq.h_view(i,j) = cut_ljsqm;
+  k_cut_ljsq.h_view(i,j) = k_cut_ljsq.h_view(j,i) = cut_ljsqm;
  k_cut_ljsq.h_view(j,i) = k_cut_ljsq.h_view(i,j);
  k_cut_ljsq.template modify<LMPHostType>();
-  k_cut_coulsq.h_view(i,j) = cut_coulsqm;
+  k_cut_coulsq.h_view(i,j) = k_cut_coulsq.h_view(j,i) = cut_coulsqm;
  k_cut_coulsq.h_view(j,i) = k_cut_coulsq.h_view(i,j);
  k_cut_coulsq.template modify<LMPHostType>();
  k_params.template modify<LMPHostType>();
--- a/src/KOKKOS/pair_lj_charmm_coul_long_kokkos.h
+++ b/src/KOKKOS/pair_lj_charmm_coul_long_kokkos.h
@ -75,7 +75,7 @@ class PairLJCharmmCoulLongKokkos : public PairLJCharmmCoulLong {
  Kokkos::DualView<params_lj_coul**,Kokkos::LayoutRight,DeviceType> k_params;
  typename Kokkos::DualView<params_lj_coul**,
    Kokkos::LayoutRight,DeviceType>::t_dev_const_um params;
-  // hardwired to space for 15 atom types
+  // hardwired to space for 12 atom types
  params_lj_coul m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
--- a/src/KOKKOS/pair_lj_class2_coul_cut_kokkos.cpp
+++ b/src/KOKKOS/pair_lj_class2_coul_cut_kokkos.cpp
@ -326,11 +326,12 @@ double PairLJClass2CoulCutKokkos<DeviceType>::init_one(int i, int j)
    m_cut_ljsq[j][i] = m_cut_ljsq[i][j] = cut_ljsqm;
    m_cut_coulsq[j][i] = m_cut_coulsq[i][j] = cut_coulsqm;
  }
-  k_cutsq.h_view(i,j) = cutone*cutone;
+
  k_cutsq.h_view(i,j) = k_cutsq.h_view(j,i) = cutone*cutone;
  k_cutsq.template modify<LMPHostType>();
-  k_cut_ljsq.h_view(i,j) = cut_ljsqm;
+  k_cut_ljsq.h_view(i,j) = k_cut_ljsq.h_view(j,i) = cut_ljsqm;
  k_cut_ljsq.template modify<LMPHostType>();
-  k_cut_coulsq.h_view(i,j) = cut_coulsqm;
+  k_cut_coulsq.h_view(i,j) = k_cut_coulsq.h_view(j,i) = cut_coulsqm;
  k_cut_coulsq.template modify<LMPHostType>();
  k_params.template modify<LMPHostType>();
--- a/src/KOKKOS/pair_lj_class2_coul_cut_kokkos.h
+++ b/src/KOKKOS/pair_lj_class2_coul_cut_kokkos.h
@ -75,7 +75,7 @@ class PairLJClass2CoulCutKokkos : public PairLJClass2CoulCut {
  Kokkos::DualView<params_lj_coul**,Kokkos::LayoutRight,DeviceType> k_params;
  typename Kokkos::DualView<params_lj_coul**,
    Kokkos::LayoutRight,DeviceType>::t_dev_const_um params;
-  // hardwired to space for 15 atom types
+  // hardwired to space for 12 atom types
  params_lj_coul m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
--- a/src/KOKKOS/pair_lj_class2_coul_long_kokkos.cpp
+++ b/src/KOKKOS/pair_lj_class2_coul_long_kokkos.cpp
@ -480,11 +480,11 @@ double PairLJClass2CoulLongKokkos<DeviceType>::init_one(int i, int j)
    m_cut_coulsq[j][i] = m_cut_coulsq[i][j] = cut_coulsqm;
  }
-  k_cutsq.h_view(i,j) = cutone*cutone;
+  k_cutsq.h_view(i,j) = k_cutsq.h_view(j,i) = cutone*cutone;
  k_cutsq.template modify<LMPHostType>();
-  k_cut_ljsq.h_view(i,j) = cut_ljsqm;
+  k_cut_ljsq.h_view(i,j) = k_cut_ljsq.h_view(j,i) = cut_ljsqm;
  k_cut_ljsq.template modify<LMPHostType>();
-  k_cut_coulsq.h_view(i,j) = cut_coulsqm;
+  k_cut_coulsq.h_view(i,j) = k_cut_coulsq.h_view(j,i) = cut_coulsqm;
  k_cut_coulsq.template modify<LMPHostType>();
  k_params.template modify<LMPHostType>();
--- a/src/KOKKOS/pair_lj_class2_coul_long_kokkos.h
+++ b/src/KOKKOS/pair_lj_class2_coul_long_kokkos.h
@ -76,7 +76,7 @@ class PairLJClass2CoulLongKokkos : public PairLJClass2CoulLong {
  Kokkos::DualView<params_lj_coul**,Kokkos::LayoutRight,DeviceType> k_params;
  typename Kokkos::DualView<params_lj_coul**,
    Kokkos::LayoutRight,DeviceType>::t_dev_const_um params;
-  // hardwired to space for 15 atom types
+  // hardwired to space for 12 atom types
  params_lj_coul m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
--- a/src/KOKKOS/pair_lj_class2_kokkos.cpp
+++ b/src/KOKKOS/pair_lj_class2_kokkos.cpp
@ -258,7 +258,7 @@ double PairLJClass2Kokkos<DeviceType>::init_one(int i, int j)
    m_params[i][j] = m_params[j][i] = k_params.h_view(i,j);
    m_cutsq[j][i] = m_cutsq[i][j] = cutone*cutone;
  }
-  k_cutsq.h_view(i,j) = cutone*cutone;
+  k_cutsq.h_view(i,j) = k_cutsq.h_view(j,i) = cutone*cutone;
  k_cutsq.template modify<LMPHostType>();
  k_params.template modify<LMPHostType>();
--- a/Show More
+++ b/Show More