Merge branch 'master' into feature/cuda-docker

CONTRIBUTING.md

+## How to Contribute to OpenMM Development
+
+We welcome anyone who wants to contribute to the project, whether by adding a feature,
+fixing a bug, or improving documentation.  The process is quite simple.
+
+First, it is always best to begin by opening an issue on Github that describes the change you
+want to make.  This gives everyone a chance to discuss it before you put in a lot of work.
+For bug fixes, we will confirm that the behavior is actually a bug and that the proposed fix
+is correct.  For new features, we will decide whether the proposed feature is something we
+want and discuss possible designs for it.
+
+Once everyone is in agreement, the next step is to
+[create a pull request](https://help.github.com/en/articles/about-pull-requests) with the code changes.
+For larger features, feel free to create the pull request even before the implementaton is
+finished so as to get early feedback on the code.  When doing this, put the letters "WIP" at
+the start of the title of the pull request to indicate it is still a work in progress.
+
+For new features, consult the [New Feature Checklist](https://github.com/openmm/openmm/wiki/Checklist-when-adding-a-new-feature),
+which lists various items that need to be included before the feature can be merged (documentation,
+tests, serialization, support for all APIs, etc.).  Not every item is necessarily applicable to
+every new feature, but usually at least some of them are.
+
+The core developers will review the pull request and may suggest changes.  Simply push the
+changes to the branch that is being pulled from, and they will automatically be added to the
+pull request.  In addition, the full test suite is automatically run on every pull request,
+and rerun every time a change is added.  Once the tests are passing and everyone is satisfied
+with the code, the pull request will be merged.  Congratulations on a succesful contribution!
\ No newline at end of file

README.md

-[![Build Status](https://travis-ci.org/pandegroup/openmm.svg?branch=master)](https://travis-ci.org/pandegroup/openmm?branch=master)
+[![Build Status](https://travis-ci.org/openmm/openmm.svg?branch=master)](https://travis-ci.org/openmm/openmm?branch=master)
 [![Anaconda Cloud Badge](https://anaconda.org/omnia/openmm/badges/downloads.svg)](https://anaconda.org/omnia/openmm)
 
 ## OpenMM: A High Performance Molecular Dynamics Library

SUPPORT.md

+## How to Get Support for OpenMM
+
+There are two main venues for getting support for OpenMM: the [discussion forum](https://simtk.org/forums/viewforum.php?f=161)
+and the [Github repository](https://github.com/openmm/openmm).  There is some overlap
+between the two, but generally speaking the forum is for user oriented issues while the
+repository is for developer oriented issues.  If you have a question about how to use OpenMM
+(including writing programs that access it through its public API), post on the forum.  If
+you want to suggest a change to the code, or if you think you have found a bug,
+open an issue on Github.  The core developers monitor both, so don't worry if you aren't
+sure which one is most appropriate for your question.  We will see it either way.
+
+You also may want to consult the [documentation](http://docs.openmm.org/).  It is quite
+thorough, and you may be able to find the answer to your question.
\ No newline at end of file

devtools/forcefield-scripts/processCharmmForceField.py

@@ -377,7 +377,7 @@ for values in sorted(uniqueCmaps, key=lambda x: uniqueCmaps[x]):
 for map in cmaps:
     print('   <Torsion map="%d" class1="%s" class2="%s" class3="%s" class4="%s" class5="%s"/>' % (uniqueCmaps[map.values], map.classes[0], map.classes[1], map.classes[2], map.classes[3], map.classes[4]))
 print(' </CMAPTorsionForce>')
-print(' <NonbondedForce coulomb14scale="1.0" lj14scale="1.0">')
+print(' <NonbondedForce coulomb14scale="1.0" lj14scale="1.0" useDispersionCorrection="False">')
 for type in atomTypes:
     print('  <Atom type="%d" charge="%g" sigma="1.0" epsilon="1.0"/>' % (type.type, type.charge))
 print(' </NonbondedForce>')
@@ -464,6 +464,7 @@ print("""]
 
 customNonbondedForce = mm.CustomNonbondedForce('4*eps*((sig/r)^12-(sig/r)^6); eps=epsilon(type1, type2); sig=sigma(type1, type2)')
 customNonbondedForce.setNonbondedMethod(min(nonbondedForce.getNonbondedMethod(), 2))
+customNonbondedForce.setUseLongRangeCorrection(False)
 customNonbondedForce.addTabulatedFunction('epsilon', mm.Discrete2DFunction(numAtomClasses, numAtomClasses, epsilon))
 customNonbondedForce.addTabulatedFunction('sigma', mm.Discrete2DFunction(numAtomClasses, numAtomClasses, sigma))
 customNonbondedForce.addPerParticleParameter('type')

docs-source/usersguide/application.rst

@@ -69,7 +69,7 @@ Anaconda or Miniconda.
 
     conda install -c omnia -c conda-forge openmm
 
-This installs a version of OpenMM that is compiled to work with CUDA 9.2.
+This installs a version of OpenMM that is compiled to work with CUDA 10.1.
 Alternatively you can request a version that is compiled for a specific CUDA
 version with the command
 ::
@@ -78,7 +78,7 @@ version with the command
 
 where :code:`cuda92` should be replaced with the particular CUDA version
 installed on your computer.  Supported values are :code:`cuda75`, :code:`cuda80`,
-:code:`cuda90`, :code:`cuda91`, :code:`cuda92`, and :code:`cuda100`.  Because
+:code:`cuda90`, :code:`cuda91`, :code:`cuda92`, :code:`cuda100`, and :code:`cuda101`.  Because
 different CUDA releases are not binary compatible with each other, OpenMM can
 only work with the particular CUDA version it was compiled with.
 
@@ -1101,16 +1101,6 @@ algorithm\ :cite:`Tuckerman1992`.  This allows some forces in the system to be e
 frequently than others.  For details on how to use it, consult the API
 documentation.
 
-aMD Integrator
---------------
-
-There are three different integrator types that implement variations of the
-aMD\ :cite:`Hamelberg2007` accelerated sampling algorithm: :class:`AMDIntegrator`,
-:class:`AMDForceGroupIntegrator`, and :class:`DualAMDIntegrator`.  They
-perform integration on a modified potential energy surface to allow much faster
-sampling of conformations.  For details on how to use them, consult the API
-documentation.
-
 Compound Integrator
 -------------------
 
@@ -1371,6 +1361,55 @@ a checkpoint file every 5,000 steps, for example:
 
 Note that the checkpoint reporter will overwrite the last checkpoint file.
 
+
+Enhanced Sampling Methods
+=========================
+
+In many situations, the goal of a simulation is to sample the range of configurations
+accessible to a system.  It does not matter whether the simulation represents a
+single, physically realistic trajectory, only whether it produces a correct distribution
+of states.  In this case, a variety of methods can be used to sample configuration
+space much more quickly and efficiently than a single physical trajectory would.
+These are known as enhanced sampling methods.  OpenMM offers several that you
+can choose from.  They are briefly described here.  Consult the API documentation
+for more detailed descriptions and example code.
+
+Simulated Tempering
+-------------------
+
+Simulated tempering\ :cite:`Marinari1992` involves making frequent changes to the
+temperature of a simulation.  At high temperatures, it can quickly cross energy barriers
+and explore a wide range of configurations.  At lower temperatures, it more thoroughly
+explores each local region of configuration space.  This is a powerful method to
+speed up sampling when you do not know in advance what motions you want to sample.
+Simply specify the range of temperatures to simulate and the algorithm handles
+everything for you mostly automatically.
+
+Metadynamics
+------------
+
+Metadynamics\ :cite:`Barducci2008` is used when you do know in advance what
+motions you want to sample.  You specify one or more collective variables, and the
+algorithm adds a biasing potential to make the simulation explore a wide range of
+values for those variables.  It does this by periodically adding "bumps" to the biasing
+potential at the current values of the collective variables.  This encourages the simulation
+to move away from regions it has already explored and sample a wide range of values.
+At the end of the simulation, the biasing potential can be used to calculate the
+free energy of the system as a function of the collective variables.
+
+Accelerated Molecular Dynamics (aMD)
+------------------------------------
+
+aMD\ :cite:`Hamelberg2007` is another method that can be used when you do not know in
+advance what motions you want to accelerate.  It alters the potential energy surface
+by adding a "boost" potential whenever the potential energy is below a threshold.
+This makes local minima shallower and allows more frequent transitions between them.
+The boost can be applied to the total potential energy, to just a subset of interactions
+(typically the dihedral torsions), or both.  There are separate integrator classes
+for each of these options: :class:`AMDIntegrator`, :class:`AMDForceGroupIntegrator`,
+and :class:`DualAMDIntegrator`.
+
+
 .. _model-building-and-editing:
 
 Model Building and Editing

docs-source/usersguide/references.bib

@@ -29,6 +29,18 @@
    type = {Journal Article}
 }
 
+@article{Barducci2008,
+  title = {Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method},
+  author = {Barducci, Alessandro and Bussi, Giovanni and Parrinello, Michele},
+  journal = {Physical Review Letters},
+  volume = {100},
+  issue = {2},
+  pages = {020603},
+  year = {2008},
+  publisher = {American Physical Society},
+  doi = {10.1103/PhysRevLett.100.020603},
+}
+
 @article{Berendsen1987
    author = {Berendsen, H. J. C. and Grigera, J. R. and Straatsma, T. P.},
    title = {The missing term in effective pair potentials},
@@ -290,6 +302,18 @@
    type = {Journal Article}
 }
 
+@article{Marinari1992,
+    doi = {10.1209/0295-5075/19/6/002},
+    year = 1992,
+    publisher = {{IOP} Publishing},
+    volume = {19},
+    number = {6},
+    pages = {451--458},
+    author = {E Marinari and G Parisi},
+    title = {Simulated Tempering: A New Monte Carlo Scheme},
+    journal = {Europhysics Letters ({EPL})},
+}
+
 @article{Markland2008
    author = {Markland, Thomas E. and Manolopoulos, David E.},
    title = {An efficient ring polymer contraction scheme for imaginary time path integral simulations},

olla/include/openmm/Platform.h

@@ -213,19 +213,24 @@ public:
      */
     static void loadPluginLibrary(const std::string& file);
     /**
-     * Load multiple dynamic libraries (DLLs) which contain OpenMM plugins from a single directory.
-     * This method loops over every file contained in the specified directory and calls loadPluginLibrary()
+     * Load multiple dynamic libraries (DLLs) which contain OpenMM plugins from one or more directories.
+     * Multiple fully-qualified paths can be joined together with ':' on unix-like systems
+     * (or ';' on windows-like systems); each will be searched for plugins, in-order. For example,
+     * '/foo/plugins:/bar/plugins' will search both `/foo/plugins` and `/bar/plugins`. If an
+     * identically-named plugin is encountered twice it will be loaded at both points; be careful!!!
+     *
+     * This method loops over every file contained in the specified directories and calls loadPluginLibrary()
      * for each one.  If an error occurs while trying to load a particular file, that file is simply
      * ignored. You can retrieve a list of all such errors by calling getPluginLoadFailures().
      *
-     * @param directory    the path to the directory containing libraries to load
+     * @param directory    a ':' (unix) or ';' (windows) deliminated list of paths containing libraries to load
      * @return the names of all files which were successfully loaded as libraries
      */
     static std::vector<std::string> loadPluginsFromDirectory(const std::string& directory);
     /**
      * Get the default directory from which to load plugins.  If the environment variable
      * OPENMM_PLUGIN_DIR is set, this returns its value.  Otherwise, it returns a platform
-     * specific default location.
+     * specific default location. 
      *
      * @return the path to the default plugin directory
      */

olla/src/Platform.cpp

@@ -37,7 +37,6 @@
 #include "openmm/internal/ContextImpl.h"
 #ifdef WIN32
 #include <windows.h>
-#include <sstream>
 #else
 #ifndef __PNACL__
     #include <dlfcn.h>
@@ -45,6 +44,7 @@
 #include <dirent.h>
 #include <cstdlib>
 #endif
+#include <sstream>
 #include <set>
 #include <algorithm>
 
@@ -261,31 +261,42 @@ void Platform::loadPluginLibrary(const string& file) {
 vector<string> Platform::loadPluginsFromDirectory(const string& directory) {
     vector<string> files;
     char dirSeparator;
+    char pathSeparator;
+    stringstream sdirectory(directory);
 #ifdef WIN32
     dirSeparator = '\\';
+    pathSeparator = ';';
     WIN32_FIND_DATA fileInfo;
-    string filePattern(directory + dirSeparator + "*.dll");
-    HANDLE findHandle = FindFirstFile(filePattern.c_str(), &fileInfo);
-    if (findHandle != INVALID_HANDLE_VALUE) {
-        do {
-            if (fileInfo.cFileName[0] != '.')
-                files.push_back(string(fileInfo.cFileName));
-        } while (FindNextFile(findHandle, &fileInfo));
-        FindClose(findHandle);
+
+    for (string path; std::getline(sdirectory, path, pathSeparator);) {
+        string filePattern(path + dirSeparator + "*.dll");
+        HANDLE findHandle = FindFirstFile(filePattern.c_str(), &fileInfo);
+        if (findHandle != INVALID_HANDLE_VALUE) {
+            do {
+                if (fileInfo.cFileName[0] != '.')
+                    files.push_back(path+dirSeparator+string(fileInfo.cFileName));
+            } while (FindNextFile(findHandle, &fileInfo));
+            FindClose(findHandle);
+        }
     }
     vector<HMODULE> plugins;
 #else
-    dirSeparator = '/';
     DIR* dir;
+    dirSeparator = '/';
+    pathSeparator = ':';
     struct dirent *entry;
-    dir = opendir(directory.c_str());
-    if (dir != NULL) {
-        while ((entry = readdir(dir)) != NULL) {
-            if (entry->d_name[0] != '.')
-                files.push_back(string(entry->d_name));
+
+    for (string path; std::getline(sdirectory, path, pathSeparator);) {
+        dir = opendir(path.c_str());
+        if (dir != NULL) {
+            while ((entry = readdir(dir)) != NULL) {
+                if (entry->d_name[0] != '.')
+                    files.push_back(path+dirSeparator+string(entry->d_name));
+            }
+            closedir(dir);
         }
-        closedir(dir);
     }
+
     vector<void*> plugins;
 #endif
     vector<string> loadedLibraries;
@@ -294,7 +305,7 @@ vector<string> Platform::loadPluginsFromDirectory(const string& directory) {
 
     for (unsigned int i = 0; i < files.size(); ++i) {
         try {
-            plugins.push_back(loadOneLibrary(directory+dirSeparator+files[i]));
+            plugins.push_back(loadOneLibrary(files[i]));
             loadedLibraries.push_back(files[i]);
         } catch (OpenMMException& ex) {
 	    pluginLoadFailures.push_back(ex.what());

platforms/opencl/include/OpenCLContext.h

@@ -32,6 +32,18 @@
 #include <string>
 #define __CL_ENABLE_EXCEPTIONS
 #define CL_USE_DEPRECATED_OPENCL_1_1_APIS
+#ifndef CL_DEVICE_SIMD_PER_COMPUTE_UNIT_AMD
+  #define CL_DEVICE_SIMD_PER_COMPUTE_UNIT_AMD 0x4040
+#endif
+#ifndef CL_DEVICE_SIMD_WIDTH_AMD
+  #define CL_DEVICE_SIMD_WIDTH_AMD 0x4041
+#endif
+#ifndef CL_DEVICE_SIMD_INSTRUCTION_WIDTH_AMD
+  #define CL_DEVICE_SIMD_INSTRUCTION_WIDTH_AMD 0x4042
+#endif
+#ifndef CL_DEVICE_WAVEFRONT_WIDTH_AMD
+  #define CL_DEVICE_WAVEFRONT_WIDTH_AMD 0x4043
+#endif
 #ifdef _MSC_VER
     // Prevent Windows from defining macros that interfere with other code.
     #define NOMINMAX

platforms/opencl/src/OpenCLContext.cpp

@@ -131,7 +131,6 @@ OpenCLContext::OpenCLContext(const System& system, int platformIndex, int device
                     // This attribute does not ensure that all queries are supported by the runtime (it may be an older runtime,
                     // or the CPU device) so still have to check for errors.
                     try {
-#ifdef CL_DEVICE_SIMD_WIDTH_AMD
                         processingElementsPerComputeUnit =
                             // AMD GPUs either have a single VLIW SIMD or multiple scalar SIMDs.
                             // The SIMD width is the number of threads the SIMD executes per cycle.
@@ -145,7 +144,6 @@ OpenCLContext::OpenCLContext(const System& system, int platformIndex, int device
                         // Just in case any of the queries return 0.
                         if (processingElementsPerComputeUnit <= 0)
                             processingElementsPerComputeUnit = 1;
-#endif
                     }
                     catch (cl::Error err) {
                         // Runtime does not support the queries so use default.
@@ -221,7 +219,6 @@ OpenCLContext::OpenCLContext(const System& system, int platformIndex, int device
                     // This attribute does not ensure that all queries are supported by the runtime so still have to
                     // check for errors.
                     try {
-#ifdef CL_DEVICE_SIMD_PER_COMPUTE_UNIT_AMD
                         // Must catch cl:Error as will fail if runtime does not support queries.
 
                         cl_uint simdPerComputeUnit = device.getInfo<CL_DEVICE_SIMD_PER_COMPUTE_UNIT_AMD>();
@@ -230,12 +227,15 @@ OpenCLContext::OpenCLContext(const System& system, int platformIndex, int device
                         // If the GPU has multiple SIMDs per compute unit then it is uses the scalar instruction
                         // set instead of the VLIW instruction set. It therefore needs more thread blocks per
                         // compute unit to hide memory latency.
-                        if (simdPerComputeUnit > 1)
-                            numThreadBlocksPerComputeUnit = 4 * simdPerComputeUnit;
+                        if (simdPerComputeUnit > 1) {
+                            if (simdWidth == 32)
+                                numThreadBlocksPerComputeUnit = 6*simdPerComputeUnit; // Navi seems to like more thread blocks than older GPUs
+                            else
+                                numThreadBlocksPerComputeUnit = 4*simdPerComputeUnit;
+                        }
 
                         // If the queries are supported then must be newer than SDK 2.4.
                         amdPostSdk2_4 = true;
-#endif
                     }
                     catch (cl::Error err) {
                         // Runtime does not support the query so is unlikely to be the newer scalar GPU.
@@ -254,7 +254,7 @@ OpenCLContext::OpenCLContext(const System& system, int platformIndex, int device
         if (supportsDoublePrecision)
             compilationDefines["SUPPORTS_DOUBLE_PRECISION"] = "";
         if (simdWidth >= 32)
-            compilationDefines["SYNC_WARPS"] = "";
+            compilationDefines["SYNC_WARPS"] = "mem_fence(CLK_LOCAL_MEM_FENCE)";
         else
             compilationDefines["SYNC_WARPS"] = "barrier(CLK_LOCAL_MEM_FENCE)";
         vector<cl::Device> contextDevices;
@@ -729,7 +729,7 @@ void OpenCLContext::clearAutoclearBuffers() {
         executeKernel(clearTwoBuffersKernel, max(autoclearBufferSizes[base], autoclearBufferSizes[base+1]), 128);
     }
     else if (total-base == 1) {
-        clearBuffer(*autoclearBuffers[base], autoclearBufferSizes[base]);
+        clearBuffer(*autoclearBuffers[base], autoclearBufferSizes[base]*4);
     }
 }
 

wrappers/python/simtk/openmm/app/__init__.py

@@ -6,7 +6,7 @@ from __future__ import absolute_import
 __docformat__ = "epytext en"
 
 __author__ = "Peter Eastman"
-__copyright__ = "Copyright 2016, Stanford University and Peter Eastman"
+__copyright__ = "Copyright 2016-2019, Stanford University and Peter Eastman"
 __credits__ = []
 __license__ = "MIT"
 __maintainer__ = "Peter Eastman"
@@ -32,6 +32,8 @@ from .checkpointreporter import CheckpointReporter
 from .charmmcrdfiles import CharmmCrdFile, CharmmRstFile
 from .charmmparameterset import CharmmParameterSet
 from .charmmpsffile import CharmmPsfFile, CharmmPSFWarning
+from .simulatedtempering import SimulatedTempering
+from .metadynamics import Metadynamics, BiasVariable
 
 # Enumerated values
 

wrappers/python/simtk/openmm/app/charmmcrdfiles.py

@@ -47,7 +47,7 @@ ONE_TIMESCALE = 1 / TIMESCALE
 class CharmmCrdFile(object):
     """
     Reads and parses a CHARMM coordinate file (.crd) into its components,
-    namely the coordinates, CHARMM atom types, resid, resname, etc.
+    namely the coordinates, CHARMM atom types, resname, etc.
 
     Attributes
     ----------
@@ -69,7 +69,6 @@ class CharmmCrdFile(object):
         self.atomno = []                   # Atom number
         self.resno = []                    # Residue number
         self.resname = []                  # Residue name
-        self.resid = []                    # Residue ID
         self.attype = []                   # Atom type
         self.positions = []                # 3N atomic coordinates
         self.title = []                    # .crd file title block
@@ -105,7 +104,7 @@ class CharmmCrdFile(object):
             try:
                 self.natom = int(line.strip().split()[0])
     
-                for row in range(self.natom):
+                for _ in range(self.natom):
                     line = crdfile.readline().strip().split()
                     self.atomno.append(int(line[0]))
                     self.resno.append(int(line[1]))
@@ -114,7 +113,6 @@ class CharmmCrdFile(object):
                     pos = Vec3(float(line[4]), float(line[5]), float(line[6]))
                     self.positions.append(pos)
                     self.segid.append(line[7])
-                    self.resid.append(int(line[8]))
                     self.weighting.append(float(line[9]))
     
                 if self.natom != len(self.positions):
@@ -124,7 +122,7 @@ class CharmmCrdFile(object):
                                            len(self.positions))
                     )
     
-            except (ValueError, IndexError) as e:
+            except (ValueError, IndexError):
                 raise CharmmFileError('Error parsing CHARMM coordinate file')
 
         # Apply units to the positions now. Do it this way to allow for

wrappers/python/simtk/openmm/app/charmmpsffile.py

@@ -166,7 +166,7 @@ class CharmmPsfFile(object):
     GB_FORCE_GROUP = 6
 
     @_catchindexerror
-    def __init__(self, psf_name):
+    def __init__(self, psf_name, periodicBoxVectors=None, unitCellDimensions=None):
         """Opens and parses a PSF file, then instantiates a CharmmPsfFile
         instance from the data.
 
@@ -174,6 +174,11 @@ class CharmmPsfFile(object):
         ----------
         psf_name : str
             Name of the PSF file (it must exist)
+        periodicBoxVectors : tuple of Vec3
+            the vectors defining the periodic box
+        unitCellDimensions : Vec3
+            the dimensions of the crystallographic unit cell.  For
+            non-rectangular unit cells, specify periodicBoxVectors instead.
 
         Raises
         ------
@@ -358,7 +363,10 @@ class CharmmPsfFile(object):
         set_molecules(atom_list)
         molecule_list = [atom.marked for atom in atom_list]
         if len(holder) == len(atom_list):
-            if molecule_list != holder:
+            if len(molecule_list) != len(holder):
+                # The MOLNT section is only used for fluctuating charge models,
+                # which are currently not supported anyway.
+                # Therefore, we only check the lengths of the lists now rather than their contents.
                 warnings.warn('Detected PSF molecule section that is WRONG. '
                               'Resetting molecularity.', CharmmPSFWarning)
         # We have a CHARMM PSF file; now do NUMLP/NUMLPH sections
@@ -449,7 +457,14 @@ class CharmmPsfFile(object):
         self.group_list = group_list
         self.title = title
         self.flags = psf_flags
-        self.box_vectors = None
+        if unitCellDimensions is not None:
+            if periodicBoxVectors is not None:
+                raise ValueError("specify either periodicBoxVectors or unitCellDimensions, but not both")
+            if u.is_quantity(unitCellDimensions):
+                unitCellDimensions = unitCellDimensions.value_in_unit(u.nanometers)
+            self.box_vectors = (Vec3(unitCellDimensions[0], 0, 0), Vec3(0, unitCellDimensions[1], 0), Vec3(0, 0, unitCellDimensions[2]))*u.nanometers
+        else:
+            self.box_vectors = periodicBoxVectors
 
     @staticmethod
     def _convert(string, type, message):
@@ -710,7 +725,7 @@ class CharmmPsfFile(object):
                 last_residue = None
             if resid != last_residue:
                 last_residue = resid
-                residue = topology.addResidue(atom.residue.resname, chain, resid)
+                residue = topology.addResidue(atom.residue.resname, chain, str(atom.residue.idx), atom.residue.inscode)
             if atom.type is not None:
                 # This is the most reliable way of determining the element
                 atomic_num = atom.type.atomic_number
@@ -1082,6 +1097,7 @@ class CharmmPsfFile(object):
         # Add nonbonded terms now
         if verbose: print('Adding nonbonded interactions...')
         force = mm.NonbondedForce()
+        force.setUseDispersionCorrection(False)
         force.setForceGroup(self.NONBONDED_FORCE_GROUP)
         if not hasbox: # non-periodic
             if nonbondedMethod is ff.NoCutoff:
@@ -1221,7 +1237,6 @@ class CharmmPsfFile(object):
             if (nonbondedMethod in (ff.PME, ff.LJPME, ff.Ewald, ff.CutoffPeriodic)):
                 cforce.setNonbondedMethod(cforce.CutoffPeriodic)
                 cforce.setCutoffDistance(nonbondedCutoff)
-                cforce.setUseLongRangeCorrection(True)
             elif nonbondedMethod is ff.NoCutoff:
                 cforce.setNonbondedMethod(cforce.NoCutoff)
             elif nonbondedMethod is ff.CutoffNonPeriodic:

wrappers/python/simtk/openmm/app/data/charmm36.xml

@@ -102920,7 +102920,7 @@
     <Torsion map="6" type1="C" type2="NH1" type3="CT" type4="C" type5="NH1"/>
     <Torsion map="7" type1="C" type2="N" type3="CT2" type4="C" type5="N"/>
   </CMAPTorsionForce>
-  <NonbondedForce coulomb14scale="1.0" lj14scale="1.0">
+  <NonbondedForce coulomb14scale="1.0" lj14scale="1.0" useDispersionCorrection="False">
     <UseAttributeFromResidue name="charge"/>
     <Atom epsilon="0.0" sigma="1.0" type="H"/>
     <Atom epsilon="0.0" sigma="1.0" type="HC"/>
@@ -103335,7 +103335,7 @@
     <Atom epsilon="0.0" sigma="1.0" type="NE"/>
     <Atom epsilon="0.0" sigma="1.0" type="DUM"/>
   </NonbondedForce>
-  <LennardJonesForce lj14scale="1.0">
+  <LennardJonesForce lj14scale="1.0" useDispersionCorrection="False">
     <Atom epsilon="0.192464" sigma="0.04000135244450124" type="H"/>
     <Atom epsilon="0.192464" sigma="0.04000135244450124" type="HC"/>
     <Atom epsilon="0.092048" sigma="0.2351972615890496" type="HA"/>

wrappers/python/simtk/openmm/app/data/charmm_polar_2013.xml

@@ -7200,7 +7200,7 @@
    <Torsion map="2" class1="CD2O1A" class2="ND3A3" class3="CD31C" class4="CD2O1A" class5="ND2A2"/>
    <Torsion map="2" class1="CD2O1A" class2="ND3A3" class3="CD31C" class4="CD2O1A" class5="ND3A3"/>
  </CMAPTorsionForce>
- <NonbondedForce coulomb14scale="1.0" lj14scale="1.0">
+ <NonbondedForce coulomb14scale="1.0" lj14scale="1.0" useDispersionCorrection="False">
   <Atom type="0" charge="1.71639" sigma="1.0" epsilon="1.0"/>
   <Atom type="1" charge="-1.11466" sigma="1.0" epsilon="1.0"/>
   <Atom type="2" charge="0.55733" sigma="1.0" epsilon="1.0"/>
@@ -11370,6 +11370,7 @@ sigma14 = [
 customNonbondedForce = mm.CustomNonbondedForce('4*eps*((sig/r)^12-(sig/r)^6); eps=epsilon(type1, type2); sig=sigma(type1, type2)')
 customNonbondedForce.setNonbondedMethod(min(nonbondedForce.getNonbondedMethod(), 2))
+customNonbondedForce.setUseLongRangeCorrection(False)
 customNonbondedForce.setCutoffDistance(nonbondedForce.getCutoffDistance())
 customNonbondedForce.addTabulatedFunction('epsilon', mm.Discrete2DFunction(numAtomClasses, numAtomClasses, epsilon))
 customNonbondedForce.addTabulatedFunction('sigma', mm.Discrete2DFunction(numAtomClasses, numAtomClasses, sigma))

wrappers/python/simtk/openmm/app/dcdreporter.py

@@ -101,7 +101,10 @@ class DCDReporter(object):
         """
 
         if self._dcd is None:
-            self._dcd = DCDFile(self._out, simulation.topology, simulation.integrator.getStepSize(), 0, self._reportInterval, self._append)
+            self._dcd = DCDFile(
+                self._out, simulation.topology, simulation.integrator.getStepSize(),
+                simulation.currentStep, self._reportInterval, self._append
+            )
         self._dcd.writeModel(state.getPositions(), periodicBoxVectors=state.getPeriodicBoxVectors())
 
     def __del__(self):

wrappers/python/simtk/openmm/app/forcefield.py

@@ -37,6 +37,7 @@ import os
 import itertools
 import xml.etree.ElementTree as etree
 import math
+import warnings
 from math import sqrt, cos
 from copy import deepcopy
 from collections import defaultdict
@@ -2273,10 +2274,11 @@ class NonbondedGenerator(object):
 
     SCALETOL = 1e-5
 
-    def __init__(self, forcefield, coulomb14scale, lj14scale):
+    def __init__(self, forcefield, coulomb14scale, lj14scale, useDispersionCorrection):
         self.ff = forcefield
         self.coulomb14scale = coulomb14scale
         self.lj14scale = lj14scale
+        self.useDispersionCorrection = useDispersionCorrection
         self.params = ForceField._AtomTypeParameters(forcefield, 'NonbondedForce', 'Atom', ('charge', 'sigma', 'epsilon'))
 
     def registerAtom(self, parameters):
@@ -2285,15 +2287,31 @@ class NonbondedGenerator(object):
     @staticmethod
     def parseElement(element, ff):
         existing = [f for f in ff._forces if isinstance(f, NonbondedGenerator)]
+        if 'useDispersionCorrection' in element.attrib:
+            useDispersionCorrection = bool(eval(element.attrib.get('useDispersionCorrection')))
+        else:
+            useDispersionCorrection = None
         if len(existing) == 0:
-            generator = NonbondedGenerator(ff, float(element.attrib['coulomb14scale']), float(element.attrib['lj14scale']))
+            generator = NonbondedGenerator(
+                ff,
+                float(element.attrib['coulomb14scale']),
+                float(element.attrib['lj14scale']),
+                useDispersionCorrection
+            )
             ff.registerGenerator(generator)
         else:
             # Multiple <NonbondedForce> tags were found, probably in different files.  Simply add more types to the existing one.
             generator = existing[0]
-            if abs(generator.coulomb14scale - float(element.attrib['coulomb14scale'])) > NonbondedGenerator.SCALETOL or \
-                    abs(generator.lj14scale - float(element.attrib['lj14scale'])) > NonbondedGenerator.SCALETOL:
+            if (abs(generator.coulomb14scale - float(element.attrib['coulomb14scale'])) > NonbondedGenerator.SCALETOL
+                or abs(generator.lj14scale - float(element.attrib['lj14scale'])) > NonbondedGenerator.SCALETOL
+            ):
                 raise ValueError('Found multiple NonbondedForce tags with different 1-4 scales')
+            if (
+                    generator.useDispersionCorrection is not None
+                    and useDispersionCorrection is not None
+                    and generator.useDispersionCorrection != useDispersionCorrection
+            ):
+                raise ValueError('Found multiple NonbondedForce tags with different useDispersionCorrection settings.')
         generator.params.parseDefinitions(element)
 
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
@@ -2317,7 +2335,18 @@ class NonbondedGenerator(object):
         if 'ewaldErrorTolerance' in args:
             force.setEwaldErrorTolerance(args['ewaldErrorTolerance'])
         if 'useDispersionCorrection' in args:
-            force.setUseDispersionCorrection(bool(args['useDispersionCorrection']))
+            lrc_keyword = bool(args['useDispersionCorrection'])
+            if self.useDispersionCorrection is not None and lrc_keyword != self.useDispersionCorrection:
+                warnings.warn(
+                    "Conflicting settings for useDispersionCorrection in createSystem() and forcefield file. "
+                    "Using the one specified in createSystem()."
+                )
+            force.setUseDispersionCorrection(lrc_keyword)
+        elif self.useDispersionCorrection is not None:
+            force.setUseDispersionCorrection(self.useDispersionCorrection)
+        else:
+            # by default
+            force.setUseDispersionCorrection(True)
         sys.addForce(force)
 
     def postprocessSystem(self, sys, data, args):
@@ -2334,10 +2363,11 @@ parsers["NonbondedForce"] = NonbondedGenerator.parseElement
 class LennardJonesGenerator(object):
     """A NBFix generator to construct the L-J force with NBFIX implemented as a lookup table"""
 
-    def __init__(self, forcefield, lj14scale):
+    def __init__(self, forcefield, lj14scale, useDispersionCorrection):
         self.ff = forcefield
         self.nbfixTypes = {}
         self.lj14scale = lj14scale
+        self.useDispersionCorrection = useDispersionCorrection
         self.ljTypes = ForceField._AtomTypeParameters(forcefield, 'LennardJonesForce', 'Atom', ('sigma', 'epsilon'))
 
     def registerNBFIX(self, parameters):
@@ -2356,14 +2386,28 @@ class LennardJonesGenerator(object):
     @staticmethod
     def parseElement(element, ff):
         existing = [f for f in ff._forces if isinstance(f, LennardJonesGenerator)]
+        if 'useDispersionCorrection' in element.attrib:
+            useDispersionCorrection = bool(eval(element.attrib.get('useDispersionCorrection')))
+        else:
+            useDispersionCorrection = None
         if len(existing) == 0:
-            generator = LennardJonesGenerator(ff, float(element.attrib['lj14scale']))
+            generator = LennardJonesGenerator(
+                ff,
+                float(element.attrib['lj14scale']),
+                useDispersionCorrection=useDispersionCorrection
+            )
             ff.registerGenerator(generator)
         else:
             # Multiple <LennardJonesForce> tags were found, probably in different files
             generator = existing[0]
             if abs(generator.lj14scale - float(element.attrib['lj14scale'])) > NonbondedGenerator.SCALETOL:
                 raise ValueError('Found multiple LennardJonesForce tags with different 1-4 scales')
+            if (
+                    generator.useDispersionCorrection is not None
+                    and useDispersionCorrection is not None
+                    and generator.useDispersionCorrection != useDispersionCorrection
+            ):
+                raise ValueError('Found multiple LennardJonesForce tags with different useDispersionCorrection settings.')
         for LJ in element.findall('Atom'):
             generator.registerLennardJones(LJ.attrib)
         for Nbfix in element.findall('NBFixPair'):
@@ -2434,12 +2478,24 @@ class LennardJonesGenerator(object):
         if args['switchDistance'] is not None:
             self.force.setUseSwitchingFunction(True)
             self.force.setSwitchingDistance(args['switchDistance'])
+        if 'useDispersionCorrection' in args:
+            lrc_keyword = bool(args['useDispersionCorrection'])
+            if self.useDispersionCorrection is not None and lrc_keyword != self.useDispersionCorrection:
+                warnings.warn(
+                    "Conflicting settings for useDispersionCorrection in createSystem() and forcefield file. "
+                    "Using the one specified in createSystem()."
+                )
+            self.force.setUseLongRangeCorrection(lrc_keyword)
+        elif self.useDispersionCorrection is not None:
+            self.force.setUseLongRangeCorrection(self.useDispersionCorrection)
+        else:
+            # by default
+            self.force.setUseLongRangeCorrection(True)
 
         # Add the particles
 
         for atom in data.atoms:
             self.force.addParticle((typeToMergedType[data.atomType[atom]],))
-        self.force.setUseLongRangeCorrection(True)
         self.force.setCutoffDistance(nonbondedCutoff)
         sys.addForce(self.force)
 
@@ -2506,9 +2562,12 @@ class GBSAOBCGenerator(object):
         generator.params.parseDefinitions(element)
 
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
-        methodMap = {NoCutoff:mm.NonbondedForce.NoCutoff,
-                     CutoffNonPeriodic:mm.NonbondedForce.CutoffNonPeriodic,
-                     CutoffPeriodic:mm.NonbondedForce.CutoffPeriodic}
+        methodMap = {NoCutoff:mm.GBSAOBCForce.NoCutoff,
+                     CutoffNonPeriodic:mm.GBSAOBCForce.CutoffNonPeriodic,
+                     CutoffPeriodic:mm.GBSAOBCForce.CutoffPeriodic,
+                     Ewald:mm.GBSAOBCForce.CutoffPeriodic,
+                     PME:mm.GBSAOBCForce.CutoffPeriodic,
+                     LJPME:mm.GBSAOBCForce.CutoffPeriodic}
         if nonbondedMethod not in methodMap:
             raise ValueError('Illegal nonbonded method for GBSAOBCForce')
         force = mm.GBSAOBCForce()
@@ -2743,7 +2802,10 @@ class CustomNonbondedGenerator(object):
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
         methodMap = {NoCutoff:mm.CustomNonbondedForce.NoCutoff,
                      CutoffNonPeriodic:mm.CustomNonbondedForce.CutoffNonPeriodic,
-                     CutoffPeriodic:mm.CustomNonbondedForce.CutoffPeriodic}
+                     CutoffPeriodic:mm.CustomNonbondedForce.CutoffPeriodic,
+                     Ewald:mm.CustomNonbondedForce.CutoffPeriodic,
+                     PME:mm.CustomNonbondedForce.CutoffPeriodic,
+                     LJPME:mm.CustomNonbondedForce.CutoffPeriodic}
         if nonbondedMethod not in methodMap:
             raise ValueError('Illegal nonbonded method for CustomNonbondedForce')
         force = mm.CustomNonbondedForce(self.energy)
@@ -2803,7 +2865,10 @@ class CustomGBGenerator(object):
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
         methodMap = {NoCutoff:mm.CustomGBForce.NoCutoff,
                      CutoffNonPeriodic:mm.CustomGBForce.CutoffNonPeriodic,
-                     CutoffPeriodic:mm.CustomGBForce.CutoffPeriodic}
+                     CutoffPeriodic:mm.CustomGBForce.CutoffPeriodic,
+                     Ewald:mm.CustomGBForce.CutoffPeriodic,
+                     PME:mm.CustomGBForce.CutoffPeriodic,
+                     LJPME:mm.CustomGBForce.CutoffPeriodic}
         if nonbondedMethod not in methodMap:
             raise ValueError('Illegal nonbonded method for CustomGBForce')
         force = mm.CustomGBForce()
@@ -2886,7 +2951,10 @@ class CustomHbondGenerator(object):
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
         methodMap = {NoCutoff:mm.CustomHbondForce.NoCutoff,
                      CutoffNonPeriodic:mm.CustomHbondForce.CutoffNonPeriodic,
-                     CutoffPeriodic:mm.CustomHbondForce.CutoffPeriodic}
+                     CutoffPeriodic:mm.CustomHbondForce.CutoffPeriodic,
+                     Ewald:mm.CustomHbondForce.CutoffPeriodic,
+                     PME:mm.CustomHbondForce.CutoffPeriodic,
+                     LJPME:mm.CustomHbondForce.CutoffPeriodic}
         if nonbondedMethod not in methodMap:
             raise ValueError('Illegal nonbonded method for CustomNonbondedForce')
         force = mm.CustomHbondForce(self.energy)
@@ -3024,7 +3092,10 @@ class CustomManyParticleGenerator(object):
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
         methodMap = {NoCutoff:mm.CustomManyParticleForce.NoCutoff,
                      CutoffNonPeriodic:mm.CustomManyParticleForce.CutoffNonPeriodic,
-                     CutoffPeriodic:mm.CustomManyParticleForce.CutoffPeriodic}
+                     CutoffPeriodic:mm.CustomManyParticleForce.CutoffPeriodic,
+                     Ewald:mm.CustomManyParticleForce.CutoffPeriodic,
+                     PME:mm.CustomManyParticleForce.CutoffPeriodic,
+                     LJPME:mm.CustomManyParticleForce.CutoffPeriodic}
         if nonbondedMethod not in methodMap:
             raise ValueError('Illegal nonbonded method for CustomManyParticleForce')
         force = mm.CustomManyParticleForce(self.particlesPerSet, self.energy)

wrappers/python/simtk/openmm/app/internal/amber_file_parser.py

@@ -989,6 +989,8 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
                 gb = mm.GBSAOBCForce()
                 gb.setSoluteDielectric(soluteDielectric)
                 gb.setSolventDielectric(solventDielectric)
+                if gbsaModel is None:
+                    gb.setSurfaceAreaEnergy(0)
         elif gbmodel == 'GBn':
             gb = customgb.GBSAGBnForce(solventDielectric, soluteDielectric, gbsaModel, cutoff, implicitSolventKappa)
         elif gbmodel == 'GBn2':

wrappers/python/simtk/openmm/app/internal/customgbforces.py

@@ -250,7 +250,7 @@ def _bondi_radii(topology):
         E.silicon:      2.1,
         E.phosphorus:   1.85,
         E.sulfur:       1.8,
-        E.chlorine:     1.5,
+        E.chlorine:     1.7,
     }
 
     if _have_numpy:
@@ -272,7 +272,7 @@ def _mbondi_radii(topology, all_bonds = None):
         E.silicon:      2.1,
         E.phosphorus:   1.85,
         E.sulfur:       1.8,
-        E.chlorine:     1.5,
+        E.chlorine:     1.7,
     }
     if _have_numpy:
         radii = numpy.empty(topology.getNumAtoms(), numpy.double)
@@ -310,7 +310,7 @@ def _mbondi2_radii(topology, all_bonds = None):
         E.silicon:      2.1,
         E.phosphorus:   1.85,
         E.sulfur:       1.8,
-        E.chlorine:     1.5,
+        E.chlorine:     1.7,
     }
     if _have_numpy:
         radii = numpy.empty(topology.getNumAtoms(), numpy.double)
@@ -541,7 +541,7 @@ class GBSAHCTForce(CustomAmberGBForceBase):
         self.addPerParticleParameter("charge")
         self.addPerParticleParameter("or") # Offset radius
         self.addPerParticleParameter("sr") # Scaled offset radius
-        self.addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
+        self.addComputedValue("I", "select(step(r+sr2-or1), 0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r), 0);"
                                    "U=r+sr2;"
                                    "L=max(or1, D);"
                                    "D=abs(r-sr2)",
@@ -607,7 +607,7 @@ class GBSAOBC1Force(CustomAmberGBForceBase):
         self.addPerParticleParameter("charge")
         self.addPerParticleParameter("or") # Offset radius
         self.addPerParticleParameter("sr") # Scaled offset radius
-        self.addComputedValue("I",  "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
+        self.addComputedValue("I",  "select(step(r+sr2-or1), 0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r), 0);"
                                     "U=r+sr2;"
                                     "L=max(or1, D);"
                                     "D=abs(r-sr2)", CustomGBForce.ParticlePairNoExclusions)
@@ -675,7 +675,7 @@ class GBSAOBC2Force(GBSAOBC1Force):
         self.addPerParticleParameter("charge")
         self.addPerParticleParameter("or") # Offset radius
         self.addPerParticleParameter("sr") # Scaled offset radius
-        self.addComputedValue("I",  "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
+        self.addComputedValue("I",  "select(step(r+sr2-or1), 0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r), 0);"
                                     "U=r+sr2;"
                                     "L=max(or1, D);"
                                     "D=abs(r-sr2)", CustomGBForce.ParticlePairNoExclusions)

wrappers/python/simtk/openmm/app/metadynamics.py

+"""
+metadynamics.py: Well-tempered metadynamics
+
+This is part of the OpenMM molecular simulation toolkit originating from
+Simbios, the NIH National Center for Physics-Based Simulation of
+Biological Structures at Stanford, funded under the NIH Roadmap for
+Medical Research, grant U54 GM072970. See https://simtk.org.
+
+Portions copyright (c) 2018-2019 Stanford University and the Authors.
+Authors: Peter Eastman
+
+Permission is hereby granted, free of charge, to any person obtaining a
+copy of this software and associated documentation files (the "Software"),
+to deal in the Software without restriction, including without limitation
+the rights to use, copy, modify, merge, publish, distribute, sublicense,
+and/or sell copies of the Software, and to permit persons to whom the
+Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
+DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
+OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
+USE OR OTHER DEALINGS IN THE SOFTWARE.
+"""
+import simtk.openmm as mm
+import simtk.unit as unit
+from collections import namedtuple
+from functools import reduce
+import os
+import re
+try:
+    import numpy as np
+except:
+    pass
+
+
+class Metadynamics(object):
+    """Performs metadynamics.
+
+    This class implements well-tempered metadynamics, as described in Barducci et al.,
+    "Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method"
+    (https://doi.org/10.1103/PhysRevLett.100.020603).  You specify from one to three
+    collective variables whose sampling should be accelerated.  A biasing force that
+    depends on the collective variables is added to the simulation.  Initially the bias
+    is zero.  As the simulation runs, Gaussian bumps are periodically added to the bias
+    at the current location of the simulation.  This pushes the simulation away from areas
+    it has already explored, encouraging it to sample other regions.  At the end of the
+    simulation, the bias function can be used to calculate the system's free energy as a
+    function of the collective variables.
+
+    To use the class you create a Metadynamics object, passing to it the System you want
+    to simulate and a list of BiasVariable objects defining the collective variables.
+    It creates a biasing force and adds it to the System.  You then run the simulation
+    as usual, but call step() on the Metadynamics object instead of on the Simulation.
+
+    You can optionally specify a directory on disk where the current bias function should
+    periodically be written.  In addition, it loads biases from any other files in the
+    same directory and includes them in the simulation.  It loads files when the
+    Metqdynamics object is first created, and also checks for any new files every time it
+    updates its own bias on disk.
+
+    This serves two important functions.  First, it lets you stop a metadynamics run and
+    resume it later.  When you begin the new simulation, it will load the biases computed
+    in the earlier simulation and continue adding to them.  Second, it provides an easy
+    way to parallelize metadynamics sampling across many computers.  Just point all of
+    them to a shared directory on disk.  Each process will save its biases to that
+    directory, and also load in and apply the biases added by other processes.
+    """
+
+    def __init__(self, system, variables, temperature, biasFactor, height, frequency, saveFrequency=None, biasDir=None):
+        """Create a Metadynamics object.
+
+        Parameters
+        ----------
+        system: System
+            the System to simulate.  A CustomCVForce implementing the bias is created and
+            added to the System.
+        variables: list of BiasVariables
+            the collective variables to sample
+        temperature: temperature
+            the temperature at which the simulation is being run.  This is used in computing
+            the free energy.
+        biasFactor: float
+            used in scaling the height of the Gaussians added to the bias.  The collective
+            variables are sampled as if the effective temperature of the simulation were
+            temperature*biasFactor.
+        height: energy
+            the initial height of the Gaussians to add
+        frequency: int
+            the interval in time steps at which Gaussians should be added to the bias potential
+        saveFrequency: int (optional)
+            the interval in time steps at which to write out the current biases to disk.  At
+            the same time it writes biases, it also checks for updated biases written by other
+            processes and loads them in.  This must be a multiple of frequency.
+        biasDir: str (optional)
+            the directory to which biases should be written, and from which biases written by
+            other processes should be loaded
+        """
+        if not unit.is_quantity(temperature):
+            temperature = temperature*unit.kelvin
+        if not unit.is_quantity(height):
+            height = height*unit.kilojoules_per_mole
+        if biasFactor < 1.0:
+            raise ValueError('biasFactor must be >= 1')
+        if (saveFrequency is None and biasDir is not None) or (saveFrequency is not None and biasDir is None):
+            raise ValueError('Must specify both saveFrequency and biasDir')
+        if saveFrequency is not None and (saveFrequency < frequency or saveFrequency%frequency != 0):
+            raise ValueError('saveFrequency must be a multiple of frequency')
+        self.variables = variables
+        self.temperature = temperature
+        self.biasFactor = biasFactor
+        self.height = height
+        self.frequency = frequency
+        self.biasDir = biasDir
+        self.saveFrequency = saveFrequency
+        self._id = np.random.randint(0x7FFFFFFF)
+        self._saveIndex = 0
+        self._selfBias = np.zeros(tuple(v.gridWidth for v in variables))
+        self._totalBias = np.zeros(tuple(v.gridWidth for v in variables))
+        self._loadedBiases = {}
+        self._deltaT = temperature*(biasFactor-1)
+        varNames = ['cv%d' % i for i in range(len(variables))]
+        self._force = mm.CustomCVForce('table(%s)' % ', '.join(varNames))
+        for name, var in zip(varNames, variables):
+            self._force.addCollectiveVariable(name, var.force)
+        widths = [v.gridWidth for v in variables]
+        mins = [v.minValue for v in variables]
+        maxs = [v.maxValue for v in variables]
+        if len(variables) == 1:
+            self._table = mm.Continuous1DFunction(self._totalBias.flatten(), mins[0], maxs[0])
+        elif len(variables) == 2:
+            self._table = mm.Continuous2DFunction(widths[0], widths[1], self._totalBias.flatten(), mins[0], maxs[0], mins[1], maxs[1])
+        elif len(variables) == 3:
+            self._table = mm.Continuous3DFunction(widths[0], widths[1], widths[2], self._totalBias.flatten(), mins[0], maxs[0], mins[1], maxs[1], mins[2], maxs[2])
+        else:
+            raise ValueError('Metadynamics requires 1, 2, or 3 collective variables')
+        self._force.addTabulatedFunction('table', self._table)
+        self._force.setForceGroup(31)
+        system.addForce(self._force)
+        self._syncWithDisk()
+
+    def step(self, simulation, steps):
+        """Advance the simulation by integrating a specified number of time steps.
+
+        Parameters
+        ----------
+        simulation: Simulation
+            the Simulation to advance
+        steps: int
+            the number of time steps to integrate
+        """
+        stepsToGo = steps
+        while stepsToGo > 0:
+            nextSteps = stepsToGo
+            if simulation.currentStep % self.frequency == 0:
+                nextSteps = min(nextSteps, self.frequency)
+            else:
+                nextSteps = min(nextSteps, simulation.currentStep % self.frequency)
+            simulation.step(nextSteps)
+            if simulation.currentStep % self.frequency == 0:
+                position = self._force.getCollectiveVariableValues(simulation.context)
+                energy = simulation.context.getState(getEnergy=True, groups={31}).getPotentialEnergy()
+                height = self.height*np.exp(-energy/(unit.MOLAR_GAS_CONSTANT_R*self._deltaT))
+                self._addGaussian(position, height, simulation.context)
+            if self.saveFrequency is not None and simulation.currentStep % self.saveFrequency == 0:
+                self._syncWithDisk()
+            stepsToGo -= nextSteps
+
+    def getFreeEnergy(self):
+        """Get the free energy of the system as a function of the collective variables.
+
+        The result is returned as a N-dimensional NumPy array, where N is the number of collective
+        variables.  The values are in kJ/mole.  The i'th position along an axis corresponds to
+        minValue + i*(maxValue-minValue)/gridWidth.
+        """
+        return -((self.temperature+self._deltaT)/self._deltaT)*self._totalBias*unit.kilojoules_per_mole
+
+    def getCollectiveVariables(self, simulation):
+        """Get the current values of all collective variables in a Simulation."""
+        return self._force.getCollectiveVariableValues(simulation.context)
+
+    def _addGaussian(self, position, height, context):
+        """Add a Gaussian to the bias function."""
+        # Compute a Gaussian along each axis.
+
+        axisGaussians = []
+        for i,v in enumerate(self.variables):
+            x = (position[i]-v.minValue) / (v.maxValue-v.minValue)
+            if v.periodic:
+                x = x % 1.0
+            dist = np.abs(np.linspace(0, 1.0, num=v.gridWidth) - x)
+            if v.periodic:
+                dist = np.min(np.array([dist, np.abs(dist-1)]), axis=0)
+            axisGaussians.append(np.exp(-dist*dist*v.gridWidth/v.biasWidth))
+
+        # Compute their outer product.
+
+        if len(self.variables) == 1:
+            gaussian = axisGaussians[0]
+        else:
+            gaussian = reduce(np.multiply.outer, reversed(axisGaussians))
+
+        # Add it to the bias.
+
+        height = height.value_in_unit(unit.kilojoules_per_mole)
+        self._selfBias += height*gaussian
+        self._totalBias += height*gaussian
+        widths = [v.gridWidth for v in self.variables]
+        mins = [v.minValue for v in self.variables]
+        maxs = [v.maxValue for v in self.variables]
+        if len(self.variables) == 1:
+            self._table.setFunctionParameters(self._totalBias.flatten(), mins[0], maxs[0])
+        elif len(self.variables) == 2:
+            self._table.setFunctionParameters(widths[0], widths[1], self._totalBias.flatten(), mins[0], maxs[0], mins[1], maxs[1])
+        elif len(self.variables) == 3:
+            self._table.setFunctionParameters(widths[0], widths[1], widths[2], self._totalBias.flatten(), mins[0], maxs[0], mins[1], maxs[1], mins[2], maxs[2])
+        self._force.updateParametersInContext(context)
+
+    def _syncWithDisk(self):
+        """Save biases to disk, and check for updated files created by other processes."""
+        if self.biasDir is None:
+            return
+
+        # Use a safe save to write out the biases to disk, then delete the older file.
+
+        oldName = os.path.join(self.biasDir, 'bias_%d_%d.npy' % (self._id, self._saveIndex))
+        self._saveIndex += 1
+        tempName = os.path.join(self.biasDir, 'temp_%d_%d.npy' % (self._id, self._saveIndex))
+        fileName = os.path.join(self.biasDir, 'bias_%d_%d.npy' % (self._id, self._saveIndex))
+        np.save(tempName, self._selfBias)
+        os.rename(tempName, fileName)
+        if os.path.exists(oldName):
+            os.remove(oldName)
+
+        # Check for any files updated by other processes.
+
+        fileLoaded = False
+        pattern = re.compile('bias_(.*)_(.*)\.npy')
+        for filename in os.listdir(self.biasDir):
+            match = pattern.match(filename)
+            if match is not None:
+                matchId = int(match.group(1))
+                matchIndex = int(match.group(2))
+                if matchId != self._id and (matchId not in self._loadedBiases or matchIndex > self._loadedBiases[matchId].index):
+                    try:
+                        data = np.load(os.path.join(self.biasDir, filename))
+                        self._loadedBiases[matchId] = _LoadedBias(matchId, matchIndex, data)
+                        fileLoaded = True
+                    except IOError:
+                        # There's a tiny chance the file could get deleted by another process between when
+                        # we check the directory and when we try to load it.  If so, just ignore the error
+                        # and keep using whatever version of that process' biases we last loaded.
+                        pass
+
+        # If we loaded any files, recompute the total bias from all processes.
+
+        if fileLoaded:
+            self._totalBias = self._selfBias
+            for bias in self._loadedBiases.values():
+                self._totalBias += bias.bias
+
+
+class BiasVariable(object):
+    """A collective variable that can be used to bias a simulation with metadynamics."""
+
+    def __init__(self, force, minValue, maxValue, biasWidth, periodic=False, gridWidth=None):
+        """Create a BiasVariable.
+
+        Parameters
+        ----------
+        force: Force
+            the Force object whose potential energy defines the collective variable
+        minValue: float
+            the minimum value the collective variable can take.  If it should ever go below this,
+            the bias force will be set to 0.
+        maxValue: float
+            the maximum value the collective variable can take.  If it should ever go above this,
+            the bias force will be set to 0.
+        biasWidth: float
+            the width (standard deviation) of the Gaussians added to the bias during metadynamics
+        periodic: bool
+            whether this is a periodic variable, such that minValue and maxValue are physical equivalent
+        gridWidth: int
+            the number of grid points to use when tabulating the bias function.  If this is omitted,
+            a reasonable value is chosen automatically.
+        """
+        self.force = force
+        self.minValue = minValue
+        self.maxValue = maxValue
+        self.biasWidth = biasWidth
+        self.periodic = periodic
+        if gridWidth is None:
+            self.gridWidth = int(np.ceil(5*(maxValue-minValue)/biasWidth))
+        else:
+            self.gridWidth = gridWidth
+
+
+_LoadedBias = namedtuple('LoadedBias', ['id', 'index', 'bias'])