Merge pull request #3 from openmm/master

Sync with official repo

tests/TestVariableLangevinIntegrator.h

@@ -310,6 +310,24 @@ void testArgonBox() {
     ke /= 1000;
     double expected = 1.5 * numParticles * BOLTZ * temp;
     ASSERT_USUALLY_EQUAL_TOL(expected, ke, 0.01);
+
+    // Compute the mean step size.
+    
+    double meanStep = 0;
+    for (int i = 0; i < 100; i++) {
+        integrator.step(1);
+        meanStep += integrator.getStepSize();
+    }
+    meanStep /= 100;
+    double maxStep = meanStep/2;
+    
+    // Now set a limit on the step size and see if it is obeyed.
+    
+    integrator.setMaximumStepSize(maxStep);
+    for (int i = 0; i < 100; i++) {
+        integrator.step(1);
+        ASSERT(integrator.getStepSize() <= maxStep*1.000001);
+    }
 }
 
 void runPlatformTests();

tests/TestVariableVerletIntegrator.h

@@ -289,6 +289,24 @@ void testArgonBox() {
         double energy = state.getKineticEnergy() + state.getPotentialEnergy();
         ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
     }
+
+    // Compute the mean step size.
+    
+    double meanStep = 0;
+    for (int i = 0; i < 100; i++) {
+        integrator.step(1);
+        meanStep += integrator.getStepSize();
+    }
+    meanStep /= 100;
+    double maxStep = meanStep/2;
+    
+    // Now set a limit on the step size and see if it is obeyed.
+    
+    integrator.setMaximumStepSize(maxStep);
+    for (int i = 0; i < 100; i++) {
+        integrator.step(1);
+        ASSERT(integrator.getStepSize() <= maxStep*1.000001);
+    }
 }
 
 void runPlatformTests();

wrappers/Doxyfile.in

@@ -582,7 +582,7 @@ INPUT_ENCODING         = UTF-8
 # *.c *.cc *.cxx *.cpp *.c++ *.java *.ii *.ixx *.ipp *.i++ *.inl *.h *.hh *.hxx 
 # *.hpp *.h++ *.idl *.odl *.cs *.php *.php3 *.inc *.m *.mm *.py *.f90
 
-FILE_PATTERNS          = 
+FILE_PATTERNS          = *.h
 
 # The RECURSIVE tag can be used to turn specify whether or not subdirectories 
 # should be searched for input files as well. Possible values are YES and NO. 

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/charmmparameterset.py

@@ -116,6 +116,7 @@ class CharmmParameterSet(object):
         self.nbfix_types = dict()
         self.nbthole_types = dict()
         self.parametersets = []
+        self.nbxmod = 5
 
         # Load all of the files
         tops, pars, strs = [], [], []
@@ -228,8 +229,17 @@ class CharmmParameterSet(object):
         nonbonded_types = dict() # Holder
         parameterset = None
         read_first_nonbonded = False
+        previous = ''
         for line in f:
-            line = line.strip()
+            line = previous+line.strip()
+            previous = ''
+            if line.endswith('-'):
+                # This will be continued on the next line.
+                previous = line[:-1]
+                continue
+            if line.startswith('!'):
+                # This is a comment.
+                continue
             if not line:
                 # This is a blank line
                 continue
@@ -258,6 +268,12 @@ class CharmmParameterSet(object):
             if line.startswith('NONBONDED'):
                 read_first_nonbonded = False
                 section = 'NONBONDED'
+                fields = line.upper().split()
+                if 'NBXMOD' in fields:
+                    nbxmod = int(fields[fields.index('NBXMOD')+1])
+                    if nbxmod not in list(range(-5, 6)):
+                        raise CharmmFileError('Unsupported value for NBXMOD: %d' % nbxmod)
+                    self.nbxmod = nbxmod
                 continue
             if line.startswith('NBFIX'):
                 section = 'NBFIX'

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
@@ -919,7 +934,7 @@ class CharmmPsfFile(object):
                 atom1=lpsite[1]
                 atom2=lpsite[2]
                 atom3=lpsite[3]
-                if atom3 > 0: 
+                if atom3 >= 0: 
                     if lpsite[4] > 0 : # relative lonepair type
                         r = lpsite[4] /10.0 # in nanometer
                         xweights = [-1.0, 0.0, 1.0]
@@ -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:
@@ -1308,26 +1323,28 @@ class CharmmPsfFile(object):
         # Add 1-4 interactions
         excluded_atom_pairs = set() # save these pairs so we don't zero them out
         sigma_scale = 2**(-1/6)
-        for tor in self.dihedral_parameter_list:
-            # First check to see if atoms 1 and 4 are already excluded because
-            # they are 1-2 or 1-3 pairs (would happen in 6-member rings or
-            # fewer). Then check that they're not already added as exclusions
-            if tor.atom1 in tor.atom4.bond_partners: continue
-            if tor.atom1 in tor.atom4.angle_partners: continue
-            key = min((tor.atom1.idx, tor.atom4.idx),
-                      (tor.atom4.idx, tor.atom1.idx))
-            if key in excluded_atom_pairs: continue # multiterm...
-            charge_prod = (tor.atom1.charge * tor.atom4.charge)
-            epsilon = (sqrt(abs(tor.atom1.type.epsilon_14) * ene_conv *
-                            abs(tor.atom4.type.epsilon_14) * ene_conv))
-            sigma = (tor.atom1.type.rmin_14 + tor.atom4.type.rmin_14) * (
-                     length_conv * sigma_scale)
-            force.addException(tor.atom1.idx, tor.atom4.idx,
-                               charge_prod, sigma, epsilon)
-            excluded_atom_pairs.add(
-                    min((tor.atom1.idx, tor.atom4.idx),
-                        (tor.atom4.idx, tor.atom1.idx))
-            )
+        nbxmod = abs(params.nbxmod)
+        if nbxmod == 5:
+            for tor in self.dihedral_parameter_list:
+                # First check to see if atoms 1 and 4 are already excluded because
+                # they are 1-2 or 1-3 pairs (would happen in 6-member rings or
+                # fewer). Then check that they're not already added as exclusions
+                if tor.atom1 in tor.atom4.bond_partners: continue
+                if tor.atom1 in tor.atom4.angle_partners: continue
+                key = min((tor.atom1.idx, tor.atom4.idx),
+                          (tor.atom4.idx, tor.atom1.idx))
+                if key in excluded_atom_pairs: continue # multiterm...
+                charge_prod = (tor.atom1.charge * tor.atom4.charge)
+                epsilon = (sqrt(abs(tor.atom1.type.epsilon_14) * ene_conv *
+                                abs(tor.atom4.type.epsilon_14) * ene_conv))
+                sigma = (tor.atom1.type.rmin_14 + tor.atom4.type.rmin_14) * (
+                         length_conv * sigma_scale)
+                force.addException(tor.atom1.idx, tor.atom4.idx,
+                                   charge_prod, sigma, epsilon)
+                excluded_atom_pairs.add(
+                        min((tor.atom1.idx, tor.atom4.idx),
+                            (tor.atom4.idx, tor.atom1.idx))
+                )
 
         # Add excluded atoms
         # Drude and lonepairs will be excluded based on their parent atoms
@@ -1353,21 +1370,24 @@ class CharmmPsfFile(object):
                             force.addException(excludeterm[j], excludeterm[i], 0.0, 0.1, 0.0)
         # Exclude all bonds and angles, as well as the lonepair/Drude attached onto them
         for atom in self.atom_list:
-            for atom2 in atom.bond_partners:
-                if atom2.idx > atom.idx:
-                    for excludeatom in [atom.idx]+parent_exclude_list[atom.idx]:
-                        for excludeatom2 in [atom2.idx]+parent_exclude_list[atom2.idx]:
-                            force.addException(excludeatom, excludeatom2, 0.0, 0.1, 0.0)
-            for atom2 in atom.angle_partners:
-                if atom2.idx > atom.idx:
-                    for excludeatom in [atom.idx]+parent_exclude_list[atom.idx]:
-                        for excludeatom2 in [atom2.idx]+parent_exclude_list[atom2.idx]:
-                            force.addException(excludeatom, excludeatom2, 0.0, 0.1, 0.0)
-            for atom2 in atom.dihedral_partners:
-                if atom2.idx <= atom.idx: continue
-                if ((atom.idx, atom2.idx) in excluded_atom_pairs):
-                    continue
-                force.addException(atom.idx, atom2.idx, 0.0, 0.1, 0.0)
+            if nbxmod > 1:
+                for atom2 in atom.bond_partners:
+                    if atom2.idx > atom.idx:
+                        for excludeatom in [atom.idx]+parent_exclude_list[atom.idx]:
+                            for excludeatom2 in [atom2.idx]+parent_exclude_list[atom2.idx]:
+                                force.addException(excludeatom, excludeatom2, 0.0, 0.1, 0.0)
+            if nbxmod > 2:
+                for atom2 in atom.angle_partners:
+                    if atom2.idx > atom.idx:
+                        for excludeatom in [atom.idx]+parent_exclude_list[atom.idx]:
+                            for excludeatom2 in [atom2.idx]+parent_exclude_list[atom2.idx]:
+                                force.addException(excludeatom, excludeatom2, 0.0, 0.1, 0.0)
+            if nbxmod > 3:
+                for atom2 in atom.dihedral_partners:
+                    if atom2.idx <= atom.idx: continue
+                    if ((atom.idx, atom2.idx) in excluded_atom_pairs):
+                        continue
+                    force.addException(atom.idx, atom2.idx, 0.0, 0.1, 0.0)
         system.addForce(force)
 
         # Add Drude particles (Drude force)

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/data/pdbNames.xml

@@ -239,19 +239,19 @@
   <Atom name="HN2" alt1="2H" alt2="HT2" alt3="H2"/>
  </Residue>
  <Residue name="UNK" type="Protein"/>
- <Residue name="A" alt1="ADE" type="Nucleic">
+ <Residue name="A" alt1="ADE" alt2="A3" alt3="A5" type="Nucleic">
   <Atom name="H61" alt1="1H6"/>
   <Atom name="H62" alt1="2H6"/>
  </Residue>
- <Residue name="G" alt1="GUA" type="Nucleic">
+ <Residue name="G" alt1="GUA" alt2="G3" alt3="G5" type="Nucleic">
   <Atom name="H21" alt1="1H2"/>
   <Atom name="H22" alt1="2H2"/>
  </Residue>
- <Residue name="C" alt1="CYT" type="Nucleic">
+ <Residue name="C" alt1="CYT" alt2="C3" alt3="C5" type="Nucleic">
   <Atom name="H41" alt1="1H4"/>
   <Atom name="H42" alt1="2H4"/>
  </Residue>
- <Residue name="U" alt1="URA" type="Nucleic"/>
+ <Residue name="U" alt1="URA" alt2="U3" alt3="U5" type="Nucleic"/>
  <Residue name="DA" alt1="DAD" alt2="DA3" alt3="DA5" type="Nucleic">
   <Atom name="H61" alt1="1H6"/>
   <Atom name="H62" alt1="2H6"/>

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
@@ -634,6 +635,29 @@ class ForceField(object):
             atom = self.getAtomIndexByName(atom_name)
             self.addExternalBond(atom)
 
+        def areParametersIdentical(self, template2, matchingAtoms, matchingAtoms2):
+            """Get whether this template and another one both assign identical atom types and parameters to all atoms.
+            
+            Parameters
+            ----------
+            template2: _TemplateData
+                the template to compare this one to
+            matchingAtoms: list
+                the indices of atoms in this template that atoms of the residue are matched to
+            matchingAtoms2: list
+                the indices of atoms in template2 that atoms of the residue are matched to
+            """
+            atoms1 = [self.atoms[m] for m in matchingAtoms]
+            atoms2 = [template2.atoms[m] for m in matchingAtoms2]
+            if any(a1.type != a2.type or a1.parameters != a2.parameters for a1,a2 in zip(atoms1, atoms2)):
+                return False
+            # Properly comparing virtual sites really needs a much more complicated analysis.  This simple check
+            # could easily fail for templates containing vsites, even if they're actually identical.  Since we
+            # currently have no force fields that include both patches and vsites, I'm not going to worry about it now.
+            if self.virtualSites != template2.virtualSites:
+                return False
+            return True
+
     class _TemplateAtomData(object):
         """Inner class used to encapsulate data about an atom in a residue template definition."""
         def __init__(self, name, type, element, parameters={}):
@@ -868,7 +892,7 @@ class ForceField(object):
 
         Returns
         -------
-        template : _ForceFieldTemplate
+        template : _TemplateData
             The matching forcefield residue template, or None if no matches are found.
         matches : list
             a list specifying which atom of the template each atom of the residue
@@ -890,7 +914,13 @@ class ForceField(object):
                 template = allMatches[0][0]
                 matches = allMatches[0][1]
             elif len(allMatches) > 1:
-                raise Exception('Multiple matching templates found for residue %d (%s): %s.' % (res.index+1, res.name, ', '.join(match[0].name for match in allMatches)))
+                # We found multiple matches.  This is OK if and only if they assign identical types and parameters to all atoms.
+                t1, m1 = allMatches[0]
+                for t2, m2 in allMatches[1:]:
+                    if not t1.areParametersIdentical(t2, m1, m2):
+                        raise Exception('Multiple non-identical matching templates found for residue %d (%s): %s.' % (res.index+1, res.name, ', '.join(match[0].name for match in allMatches)))
+                template = allMatches[0][0]
+                matches = allMatches[0][1]
         return [template, matches]
 
     def _buildBondedToAtomList(self, topology):
@@ -2276,10 +2306,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):
@@ -2288,15 +2319,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):
@@ -2320,7 +2367,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):
@@ -2337,10 +2395,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):
@@ -2359,14 +2418,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'):
@@ -2437,12 +2510,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)
 
@@ -2509,9 +2594,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()
@@ -2746,7 +2834,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)
@@ -2806,7 +2897,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()
@@ -2889,7 +2983,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)
@@ -3027,7 +3124,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)
@@ -4363,7 +4463,7 @@ class AmoebaVdwGenerator(object):
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
 
         sigmaMap = {'ARITHMETIC':1, 'GEOMETRIC':1, 'CUBIC-MEAN':1}
-        epsilonMap = {'ARITHMETIC':1, 'GEOMETRIC':1, 'HARMONIC':1, 'HHG':1}
+        epsilonMap = {'ARITHMETIC':1, 'GEOMETRIC':1, 'HARMONIC':1, 'W-H':1, 'HHG':1}
 
         force = mm.AmoebaVdwForce()
         sys.addForce(force)

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

@@ -563,20 +563,27 @@ class GromacsTopFile(object):
                     if atomName in atomReplacements:
                         atomName = atomReplacements[atomName]
 
-                    # Try to guess the element.
-
-                    upper = atomName.upper()
-                    if upper.startswith('CL'):
-                        element = elem.chlorine
-                    elif upper.startswith('NA'):
-                        element = elem.sodium
-                    elif upper.startswith('MG'):
-                        element = elem.magnesium
+                    # Try to determine the element.
+
+                    atomicNumber = self._atomTypes[fields[1]][2]
+                    if atomicNumber is None:
+                        # Try to guess the element from the name.
+                        upper = atomName.upper()
+                        if upper.startswith('CL'):
+                            element = elem.chlorine
+                        elif upper.startswith('NA'):
+                            element = elem.sodium
+                        elif upper.startswith('MG'):
+                            element = elem.magnesium
+                        else:
+                            try:
+                                element = elem.get_by_symbol(atomName[0])
+                            except KeyError:
+                                element = None
+                    elif atomicNumber == '0':
+                        element = None
                     else:
-                        try:
-                            element = elem.get_by_symbol(atomName[0])
-                        except KeyError:
-                            element = None
+                        element = elem.Element.getByAtomicNumber(int(atomicNumber))
                     atoms.append(top.addAtom(atomName, element, r))
 
                 # Add bonds to the topology

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)
@@ -389,18 +389,20 @@ def _createEnergyTerms(force, solventDielectric, soluteDielectric, SA, cutoff, k
                                 "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))"+params, CustomGBForce.ParticlePairNoExclusions)
 
 
-_SCREEN_PARAMETERS = { # normal, GBn, GBn2
-        E.hydrogen : (0.85, 1.09085413633, 1.425952),
-        E.carbon : (0.72, 0.48435382330, 1.058554),
-        E.nitrogen : (0.79, 0.700147318409, 0.733599),
-        E.oxygen : (0.85, 1.06557401132, 1.061039),
-        E.fluorine : (0.88, 0.5, 0.5),
-        E.phosphorus : (0.86, 0.5, 0.5),
-        E.sulfur : (0.96, 0.602256336067, -0.703469),
-        None : (0.8, 0.5, 0.5) # default
+_SCREEN_PARAMETERS = { # normal, GBn, GBn2, GBn2 nucleic
+        E.hydrogen : (0.85, 1.09085413633, 1.425952, 1.696538),
+        E.carbon : (0.72, 0.48435382330, 1.058554, 1.268902),
+        E.nitrogen : (0.79, 0.700147318409, 0.733599, 1.4259728),
+        E.oxygen : (0.85, 1.06557401132, 1.061039, 0.1840098),
+        E.fluorine : (0.88, 0.5, 0.5, 0.5),
+        E.phosphorus : (0.86, 0.5, 0.5, 1.5450597),
+        E.sulfur : (0.96, 0.602256336067, -0.703469, 0.05),
+        None : (0.8, 0.5, 0.5, 0.5) # default
 }
 _SCREEN_PARAMETERS[E.deuterium] = _SCREEN_PARAMETERS[E.hydrogen]
 
+_NUCLEIC_ACID_RESIDUES = ['A', 'C', 'G', 'U', 'DA', 'DC', 'DG', 'DT']
+
 def _screen_parameter(atom):
     return _SCREEN_PARAMETERS.get(atom.element, _SCREEN_PARAMETERS[None])
 
@@ -541,7 +543,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 +609,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 +677,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)
@@ -873,7 +875,16 @@ class GBSAGBn2Force(GBSAGBnForce):
         E.oxygen:       [0.867814, 0.876635, 0.387882],
         E.sulfur:       [0.867814, 0.876635, 0.387882],
     }
-    _default_atom_params = [0.8, 4.85, 0.5]
+    _atom_params_nucleic = {
+        E.hydrogen:     [0.537050, 0.362861, 0.116704],
+        E.deuterium:    [0.537050, 0.362861, 0.116704],
+        E.carbon:       [0.331670, 0.196842, 0.093422],
+        E.nitrogen:     [0.686311, 0.463189, 0.138722],
+        E.oxygen:       [0.606344, 0.463006, 0.142262],
+        E.sulfur:       [0.606344, 0.463006, 0.142262],
+        E.phosphorus:   [0.418365, 0.290054, 0.1064245],
+    }
+    _default_atom_params = [1.0, 0.8, 4.851]
 
     @classmethod
     def getStandardParameters(cls, topology):
@@ -897,14 +908,23 @@ class GBSAGBn2Force(GBSAGBnForce):
             radii = numpy.empty([natoms,5], numpy.double)
             radii[:,0] = _mbondi3_radii(topology)/10
             for atom, rad in zip(topology.atoms(), radii):
-                rad[1] = _screen_parameter(atom)[2]
-                rad[2:] = cls._atom_params.get(atom.element, cls._default_atom_params)
+                if atom.residue.name in _NUCLEIC_ACID_RESIDUES:
+                    rad[1] = _screen_parameter(atom)[3]
+                    rad[2:] = cls._atom_params_nucleic.get(atom.element, cls._default_atom_params)
+                else:
+                    rad[1] = _screen_parameter(atom)[2]
+                    rad[2:] = cls._atom_params.get(atom.element, cls._default_atom_params)
         else:
             radii = [[r/10, 0, 0, 0, 0] for r in _mbondi3_radii(topology)]
             for atom, rad in zip(topology.atoms(), radii):
-                rad[1] = _screen_parameter(atom)[2]
-                for i, p in enumerate(cls._atom_params.get(atom.element, cls._default_atom_params)):
-                    rad[2+i] = p
+                if atom.residue.name in _NUCLEIC_ACID_RESIDUES:
+                    rad[1] = _screen_parameter(atom)[3]
+                    for i, p in enumerate(cls._atom_params_nucleic.get(atom.element, cls._default_atom_params)):
+                        rad[2+i] = p
+                else:
+                    rad[1] = _screen_parameter(atom)[2]
+                    for i, p in enumerate(cls._atom_params.get(atom.element, cls._default_atom_params)):
+                        rad[2+i] = p
         return radii
 
     def _addEnergyTerms(self):

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'])

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

@@ -258,13 +258,16 @@ class Modeller(object):
         self.topology = newTopology
         self.positions = newPositions
 
-    def _addIons(self, forcefield, replaceableMols, ionCutoff=0.05*nanometer, positiveIon='Na+', negativeIon='Cl-', ionicStrength=0*molar, neutralize=True):
+    def _addIons(self, forcefield, numWaters, replaceableMols, ionCutoff=0.05*nanometer, positiveIon='Na+', negativeIon='Cl-', ionicStrength=0*molar, neutralize=True):
         """Adds ions to the system by replacing certain molecules.
 
         Parameters
         ----------
         forcefield : ForceField
             the ForceField to use to determine the total charge of the system.
+        numWaters : int
+            the total number of water molecules in the simulation box, used to
+            calculate the number of ions / concentration to add.
         replaceableMols : dict
             the molecules to replace by ions, as a dictionary of residue:positions
         ionCutoff: distance=0.5*nanometer
@@ -326,51 +329,52 @@ class Modeller(object):
                 numPositive -= totalCharge
 
         if ionicStrength > 0 * molar:
-            numIons = (numReplaceableMols - numPositive - numNegative) * ionicStrength / (55.4 * molar)  # Pure water is about 55.4 molar (depending on temperature)
+            numIons = (numWaters - numPositive - numNegative) * ionicStrength / (55.4 * molar)  # Pure water is about 55.4 molar (depending on temperature)
             numPairs = int(floor(numIons + 0.5))
             numPositive += numPairs
             numNegative += numPairs
         totalIons = numPositive + numNegative
 
-        # Randomly select a set of waters
-        # while ensuring ions are not placed too close to each other.
-        modeller = Modeller(self.topology, self.positions)
-
-        replaceableList = list(replaceableMols.keys())
-        numAddedIons = 0
-        numTrials = 10  # Attempts to add ions N times before quitting
-        toReplace = []  # list of molecules to be replaced
-        while numAddedIons < totalIons:
-            pickedMol = random.choice(replaceableList)
-            replaceableList.remove(pickedMol)
-            # Check distance to other ions
-            for pos in ionPositions:
-                distance = norm(pos - replaceableMols[pickedMol])
-                if distance <= ionCutoff:
-                    numTrials -= 1
-                    break
-            else:
-                toReplace.append(pickedMol)
-                ionPositions.append(replaceableMols[pickedMol])
-                numAddedIons += 1
+        if totalIons > 0:
+            # Randomly select a set of waters
+            # while ensuring ions are not placed too close to each other.
+            modeller = Modeller(self.topology, self.positions)
+
+            replaceableList = list(replaceableMols.keys())
+            numAddedIons = 0
+            numTrials = 10  # Attempts to add ions N times before quitting
+            toReplace = []  # list of molecules to be replaced
+            while numAddedIons < totalIons:
+                pickedMol = random.choice(replaceableList)
+                replaceableList.remove(pickedMol)
+                # Check distance to other ions
+                for pos in ionPositions:
+                    distance = norm(pos - replaceableMols[pickedMol])
+                    if distance <= ionCutoff:
+                        numTrials -= 1
+                        break
+                else:
+                    toReplace.append(pickedMol)
+                    ionPositions.append(replaceableMols[pickedMol])
+                    numAddedIons += 1
 
-                n_trials = 10
+                    n_trials = 10
 
-            if n_trials == 0:
-                raise ValueError('Could not add more than {} ions to the system'.format(numAddedIons))
+                if n_trials == 0:
+                    raise ValueError('Could not add more than {} ions to the system'.format(numAddedIons))
 
-        # Replace waters/ions in the topology
-        modeller.delete(toReplace)
-        ionChain = modeller.topology.addChain()
-        for i, water in enumerate(toReplace):
-            element = (positiveElement if i < numPositive else negativeElement)
-            newResidue = modeller.topology.addResidue(element.symbol.upper(), ionChain)
-            modeller.topology.addAtom(element.symbol, element, newResidue)
-            modeller.positions.append(replaceableMols[water])
+            # Replace waters/ions in the topology
+            modeller.delete(toReplace)
+            ionChain = modeller.topology.addChain()
+            for i, water in enumerate(toReplace):
+                element = (positiveElement if i < numPositive else negativeElement)
+                newResidue = modeller.topology.addResidue(element.symbol.upper(), ionChain)
+                modeller.topology.addAtom(element.symbol, element, newResidue)
+                modeller.positions.append(replaceableMols[water])
 
-        # Update topology/positions
-        self.topology = modeller.topology
-        self.positions = modeller.positions
+            # Update topology/positions
+            self.topology = modeller.topology
+            self.positions = modeller.positions
 
     def addSolvent(self, forcefield, model='tip3p', boxSize=None, boxVectors=None, padding=None, numAdded=None, positiveIon='Na+', negativeIon='Cl-', ionicStrength=0*molar, neutralize=True):
         """Add solvent (both water and ions) to the model to fill a rectangular box.
@@ -621,8 +625,11 @@ class Modeller(object):
                         if atom.element == _oxygen:
                             waterPos[residue] = newPositions[atom.index]
 
+        # Total number of waters in the box
+        numTotalWaters = len(waterPos)
+
         # Add ions to neutralize the system.
-        self._addIons(forcefield, waterPos, positiveIon=positiveIon, negativeIon=negativeIon, ionicStrength=ionicStrength, neutralize=neutralize)
+        self._addIons(forcefield, numTotalWaters, waterPos, positiveIon=positiveIon, negativeIon=negativeIon, ionicStrength=ionicStrength, neutralize=neutralize)
 
     class _ResidueData:
         """Inner class used to encapsulate data about the hydrogens for a residue."""
@@ -712,6 +719,8 @@ class Modeller(object):
         automatically, this function will only add hydrogens.  It will never remove ones that are already present in the
         model, regardless of the specified pH.
 
+        In all cases, the positions of existing atoms (including existing hydrogens) are not modified.
+
         Definitions for standard amino acids and nucleotides are built in.  You can call loadHydrogenDefinitions() to load
         additional definitions for other residue types.
 
@@ -930,14 +939,16 @@ class Modeller(object):
 
         # The hydrogens were added at random positions.  Now perform an energy minimization to fix them up.
 
+        addedH = set(newIndices)  # keep track of Hs added
+
         if forcefield is not None:
             # Use the ForceField the user specified.
 
             system = forcefield.createSystem(newTopology, rigidWater=False, nonbondedMethod=CutoffNonPeriodic)
             atoms = list(newTopology.atoms())
             for i in range(system.getNumParticles()):
-                if atoms[i].element != elem.hydrogen:
-                    # This is a heavy atom, so make it immobile.
+                if i not in addedH:
+                    # Existing atom, make it immobile.
                     system.setParticleMass(i, 0)
         else:
             # Create a System that restrains the distance of each hydrogen from its parent atom
@@ -955,7 +966,7 @@ class Modeller(object):
             bondedTo = []
             for atom in newTopology.atoms():
                 nonbonded.addParticle([])
-                if atom.element != elem.hydrogen:
+                if atom.index not in addedH:  # make immobile
                     system.addParticle(0.0)
                 else:
                     system.addParticle(1.0)
@@ -1221,33 +1232,33 @@ class Modeller(object):
 
     def addMembrane(self, forcefield, lipidType='POPC', membraneCenterZ=0*nanometer, minimumPadding=1*nanometer, positiveIon='Na+', negativeIon='Cl-', ionicStrength=0*molar, neutralize=True):
         """Add a lipid membrane to the model.
-        
+
         This method actually adds both a membrane and a water box.  It is best to build them together,
         both to avoid adding waters inside the membrane and to ensure that lipid head groups are properly
         solvated.  For that reason, this method includes many of the same arguments as addSolvent().
-        
+
         The membrane is added in the XY plane, and the existing protein is assumed to already be oriented
         and positioned correctly.  When possible, it is recommended to start with a model from the
         Orientations of Proteins in Membranes (OPM) database at http://opm.phar.umich.edu.  Otherwise, it
         is up to you to select the protein position yourself.
-        
+
         The algorithm is based on the one described in Wolf et al., J. Comp. Chem. 31, pp. 2169-2174 (2010).
         It begins by tiling copies of a pre-equilibrated membrane patch to create a membrane of the desired
         size.  Next it scales down the protein by 50% along the X and Y axes.  Any lipid within a cutoff
         distance of the scaled protein is removed.  It also ensures that equal numbers of lipids are removed
         from each leaf of the membrane.  Finally, 1000 steps of molecular dynamics are performed to let
         the membrane relax while the protein is gradually scaled back up to its original size.
-        
+
         The size of the membrane and water box are determined by the minimumPadding argument.  All
         pre-existing atoms are guaranteed to be at least this far from any edge of the periodic box.  It
         is also possible for the periodic box to have more padding than requested.  In particular, it only
         adds whole copies of the pre-equilibrated membrane patch, so the box dimensions will always be
         integer multiples of the patch size.  That may lead to a larger membrane than what you requested.
-        
+
         This method has built in support for POPC and POPE lipids.  You can also build other types of
         membranes by providing a pre-equilibrated, solvated membrane patch that can be tiled in the XY
         plane to form the membrane.
-        
+
         Parameters
         ----------
         forcefield : ForceField
@@ -1282,8 +1293,8 @@ class Modeller(object):
             membraneCenterZ = membraneCenterZ.value_in_unit(nanometer)
         if is_quantity(minimumPadding):
             minimumPadding = minimumPadding.value_in_unit(nanometer)
-        
-        # Figure out how many copies of the membrane patch we need in each direction. 
+
+        # Figure out how many copies of the membrane patch we need in each direction.
 
         proteinPos = self.positions.value_in_unit(nanometer)
         proteinMinPos = Vec3(*[min((p[i] for p in proteinPos)) for i in range(3)])
@@ -1298,15 +1309,15 @@ class Modeller(object):
         patchCenterPos = (patchMinPos+patchMaxPos)/2
         nx = int(ceil((proteinSize[0]+2*minimumPadding)/patchSize[0]))
         ny = int(ceil((proteinSize[1]+2*minimumPadding)/patchSize[1]))
-       
+
         # Record the bonds for each residue.
-        
+
         resBonds = defaultdict(list)
         for bond in patch.topology.bonds():
             resBonds[bond[0].residue].append(bond)
-        
+
         # Identify which leaf of the membrane each lipid is in.
-        
+
         numLipidAtoms = 0
         resMeanZ = {}
         membraneMeanZ = 0.0
@@ -1322,17 +1333,17 @@ class Modeller(object):
                 resMeanZ[res] = sumZ/numResAtoms
         membraneMeanZ /= numLipidAtoms
         lipidLeaf = dict((res, 0 if resMeanZ[res] < membraneMeanZ else 1) for res in resMeanZ)
-        
+
         # Compute scaled positions for the protein.
-        
+
         scaledProteinPos = [None]*len(proteinPos)
         for i, p in enumerate(proteinPos):
             p = p-proteinCenterPos
             p = Vec3(0.5*p[0], 0.5*p[1], p[2])
             scaledProteinPos[i] = p+proteinCenterPos
-        
+
         # Create a new Topology for the membrane.
-        
+
         membraneTopology = Topology()
         membranePos = []
         boxSizeZ = patchSize[2]
@@ -1341,14 +1352,13 @@ class Modeller(object):
         else:
             boxSizeZ = max(boxSizeZ, proteinSize[2]+2*minimumPadding)
         membraneTopology.setUnitCellDimensions((nx*patchSize[0], ny*patchSize[1], boxSizeZ))
-        
+
         # Add membrane patches.  We exclude any water that is within a cutoff distance of either the actual or scaled
         # protein, and any lipid that is within a cutoff distance of the scaled protein.  We also keep track of how
         # many lipids have been excluded from each leaf of the membrane, so we can make sure exactly the same
         # number get removed from each leaf.
-        
+
         overlapCutoff = 0.22
-        chain = membraneTopology.addChain()
         addedWater = []
         addedLipids = []
         removedFromLeaf = [0, 0]
@@ -1395,28 +1405,32 @@ class Modeller(object):
         del cellLists
         del cells
         del proteinCells
-        
+
         # Add the lipids.
-        
+
         newAtoms = {}
         lipidChain = membraneTopology.addChain()
+        lipidResNum = 1  # renumber lipid residues to handle large patches
         for (nearest, residue, pos) in addedLipids:
             if skipFromLeaf[lipidLeaf[residue]] > 0:
                 # Remove the same number of residues from each leaf.
                 skipFromLeaf[lipidLeaf[residue]] -= 1
             else:
-                newResidue = membraneTopology.addResidue(residue.name, lipidChain, residue.id, residue.insertionCode)
+                newResidue = membraneTopology.addResidue(residue.name, lipidChain, lipidResNum, residue.insertionCode)
+                lipidResNum += 1
+
                 for atom in residue.atoms():
                     newAtom = membraneTopology.addAtom(atom.name, atom.element, newResidue, atom.id)
                     newAtoms[atom] = newAtom
                 membranePos += pos
                 for bond in resBonds[residue]:
                     membraneTopology.addBond(newAtoms[bond[0]], newAtoms[bond[1]], bond.type, bond.order)
+
         del lipidLeaf
         del addedLipids
-        
+
         # Add the solvent.
-        
+
         solventChain = membraneTopology.addChain()
         for (residue, pos) in addedWater:
             newResidue = membraneTopology.addResidue(residue.name, solventChain, residue.id, residue.insertionCode)
@@ -1432,7 +1446,7 @@ class Modeller(object):
         gc.collect()
 
         # Create a System for the lipids, then add in the protein as stationary particles.
-        
+
         system = forcefield.createSystem(membraneTopology, nonbondedMethod=CutoffPeriodic)
         proteinSystem = forcefield.createSystem(self.topology, nonbondedMethod=CutoffNonPeriodic)
         numMembraneParticles = system.getNumParticles()
@@ -1454,9 +1468,9 @@ class Modeller(object):
         del membranePos
         del scaledProteinCells
         gc.collect()
-        
+
         # Run a simulation while slowly scaling up the protein so the membrane can relax.
-        
+
         integrator = LangevinIntegrator(10.0, 50.0, 0.001)
         context = Context(system, integrator)
         context.setPositions(mergedPositions)
@@ -1479,24 +1493,24 @@ class Modeller(object):
                     mergedPositions[j+numMembraneParticles] = (weight1*proteinPos[j] + weight2*scaledProteinPos[j])
             context.setPositions(mergedPositions)
             integrator.step(20)
-        
+
         # Add the membrane to the protein.
-        
+
         modeller = Modeller(self.topology, self.positions)
         modeller.add(membraneTopology, context.getState(getPositions=True).getPositions()[:numMembraneParticles])
         modeller.topology.setPeriodicBoxVectors(membraneTopology.getPeriodicBoxVectors())
         del context
         del system
         del integrator
-        
+
         # Depending on the box size, we may need to add more water beyond what was included with the membrane patch.
-        
+
         needExtraWater = (boxSizeZ > patchSize[2])
         if needExtraWater:
             modeller.addSolvent(forcefield, neutralize=False)
-        
+
         # Record the positions of all waters that have been added.
-        
+
         waterPos = {}
         for chain in list(modeller.topology.chains())[-2:]:
             for residue in chain.residues():
@@ -1534,6 +1548,11 @@ class Modeller(object):
                         if atom.element == elem.oxygen:
                             waterPos[residue] = modeller.positions[atom.index]
 
+        # Total number of water molecules
+        # Use this number to avoid underestimating the concentration of ions
+        # in _addIons after we exclude waters close to lipids.
+        numTotalWaters = len(waterPos)
+
         # Calculate lipid Z boundaries
         lipidNames = {res.name for res in patch.topology.residues() if res.name != 'HOH'}
         lipidZMax = sys.float_info.min
@@ -1555,12 +1574,12 @@ class Modeller(object):
             if lowerZBoundary < waterZ.value_in_unit(nanometer) < upperZBoundary:
                 del waterPos[wRes]
 
-        self._addIons(forcefield, waterPos, positiveIon=positiveIon, negativeIon=negativeIon, ionicStrength=ionicStrength, neutralize=neutralize)
+        self._addIons(forcefield, numTotalWaters, waterPos, positiveIon=positiveIon, negativeIon=negativeIon, ionicStrength=ionicStrength, neutralize=neutralize)
 
 
 class _CellList(object):
     """This class organizes atom positions into cells, so the neighbors of a point can be quickly retrieved"""
-    
+
     def __init__(self, positions, maxCutoff, vectors, periodic):
         self.positions = positions[:]
         self.cells = {}

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

@@ -69,8 +69,8 @@ class PDBFile(object):
 
         Parameters
         ----------
-        file : string
-            the name of the file to load
+        file : string or file
+            the name of the file to load.  Alternatively you can pass an open file object.
         extraParticleIdentifier : string='EP'
             if this value appears in the element column for an ATOM record, the Atom's element will be set to None to mark it as an extra particle
         """
@@ -352,9 +352,11 @@ class PDBFile(object):
                     resName = res.name
                 if keepIds and len(res.id) < 5:
                     resId = res.id
-                    resIC = res.insertionCode
                 else:
                     resId = "%4d" % ((resIndex+1)%10000)
+                if len(res.insertionCode) == 1:
+                    resIC = res.insertionCode
+                else:
                     resIC = " "
                 if res.name in nonHeterogens:
                     recordName = "ATOM  "

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

@@ -397,7 +397,7 @@ class PDBxFile(object):
             for (resIndex, res) in enumerate(residues):
                 if keepIds:
                     resId = res.id
-                    resIC = (res.insertionCode if len(res.insertionCode) > 0 else '.')
+                    resIC = (res.insertionCode if res.insertionCode.strip() else '.')
                 else:
                     resId = resIndex + 1
                     resIC = '.'

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

+from __future__ import print_function
+
+"""
+simulatedtempering.py: Implements simulated tempering
+
+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) 2015 Stanford University and the Authors.
+Authors: Peter Eastman
+Contributors:
+
+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.
+"""
+__author__ = "Peter Eastman"
+__version__ = "1.0"
+
+import simtk.unit as unit
+import math
+import random
+from sys import stdout
+
+try:
+    import bz2
+    have_bz2 = True
+except: have_bz2 = False
+
+try:
+    import gzip
+    have_gzip = True
+except: have_gzip = False
+
+try:
+    import numpy
+    have_numpy = True
+except: have_numpy = False
+
+class SimulatedTempering(object):
+    """SimulatedTempering implements the simulated tempering algorithm for accelerated sampling.
+    
+    It runs a simulation while allowing the temperature to vary.  At high temperatures, it can more easily cross
+    energy barriers to explore a wider area of conformation space.  At low temperatures, it can thoroughly
+    explore each local region.  For details, see Marinari, E. and Parisi, G., Europhys. Lett. 19(6). pp. 451-458 (1992).
+    
+    The set of temperatures to sample can be specified in two ways.  First, you can explicitly provide a list
+    of temperatures by using the "temperatures" argument.  Alternatively, you can specify the minimum and
+    maximum temperatures, and the total number of temperatures to use.  The temperatures are chosen spaced
+    exponentially between the two extremes.  For example,
+    
+    st = SimulatedTempering(simulation, numTemperatures=15, minTemperature=300*kelvin, maxTemperature=450*kelvin)
+    
+    After creating the SimulatedTempering object, call step() on it to run the simulation.
+    
+    Transitions between temperatures are performed at regular intervals, as specified by the "tempChangeInterval"
+    argument.  For each transition, a new temperature is selected using the independence sampling method, as
+    described in Chodera, J. and Shirts, M., J. Chem. Phys. 135, 194110 (2011).
+    
+    Simulated tempering requires a "weight factor" for each temperature.  Ideally, these should be chosen so
+    the simulation spends equal time at every temperature.  You can specify the list of weights to use with the
+    optional "weights" argument.  If this is omitted, weights are selected automatically using the Wang-Landau
+    algorithm as described in Wang, F. and Landau, D. P., Phys. Rev. Lett. 86(10), pp. 2050-2053 (2001).
+    
+    To properly analyze the results of the simulation, it is important to know the temperature and weight factors
+    at every point in time.  The SimulatedTempering object functions as a reporter, writing this information
+    to a file or stdout at regular intervals (which should match the interval at which you save frames from the
+    simulation).  You can specify the output file and reporting interval with the "reportFile" and "reportInterval"
+    arguments.
+    """
+
+    def __init__(self, simulation, temperatures=None, numTemperatures=None, minTemperature=None, maxTemperature=None, weights=None, tempChangeInterval=25, reportInterval=1000, reportFile=stdout):
+        """Create a new SimulatedTempering.
+        
+        Parameters
+        ----------
+        simulation: Simulation
+            The Simulation defining the System, Context, and Integrator to use
+        temperatures: list
+            The list of temperatures to use for tempering, in increasing order
+        numTemperatures: int
+            The number of temperatures to use for tempering.  If temperatures is not None, this is ignored.
+        minTemperature: temperature
+            The minimum temperature to use for tempering.  If temperatures is not None, this is ignored.
+        maxTemperature: temperature
+            The maximum temperature to use for tempering.  If temperatures is not None, this is ignored.
+        weights: list
+            The weight factor for each temperature.  If none, weights are selected automatically.
+        tempChangeInterval: int
+            The interval (in time steps) at which to attempt transitions between temperatures
+        reportInterval: int
+            The interval (in time steps) at which to write information to the report file
+        reportFile: string or file
+            The file to write reporting information to, specified as a file name or file object
+        """        
+        self.simulation = simulation
+        if temperatures is None:
+            if unit.is_quantity(minTemperature):
+                minTemperature = minTemperature.value_in_unit(unit.kelvin)
+            if unit.is_quantity(maxTemperature):
+                maxTemperature = maxTemperature.value_in_unit(unit.kelvin)
+            self.temperatures = [minTemperature*((float(maxTemperature)/minTemperature)**(i/float(numTemperatures-1))) for i in range(numTemperatures)]*unit.kelvin
+        else:
+            numTemperatures = len(temperatures)
+            self.temperatures = [(t.value_in_unit(unit.kelvin) if unit.is_quantity(t) else t)*unit.kelvin for t in temperatures]
+            if any(self.temperatures[i] >= self.temperatures[i+1] for i in range(numTemperatures-1)):
+                raise ValueError('The temperatures must be in strictly increasing order')
+        self.tempChangeInterval = tempChangeInterval
+        self.reportInterval = reportInterval
+        self.inverseTemperatures = [1.0/(unit.MOLAR_GAS_CONSTANT_R*t) for t in self.temperatures]
+
+        # If necessary, open the file we will write reports to.
+
+        self._openedFile = isinstance(reportFile, str)
+        if self._openedFile:
+            # Detect the desired compression scheme from the filename extension
+            # and open all files unbuffered
+            if reportFile.endswith('.gz'):
+                if not have_gzip:
+                    raise RuntimeError("Cannot write .gz file because Python could not import gzip library")
+                self._out = gzip.GzipFile(fileobj=open(reportFile, 'wb', 0))
+            elif reportFile.endswith('.bz2'):
+                if not have_bz2:
+                    raise RuntimeError("Cannot write .bz2 file because Python could not import bz2 library")
+                self._out = bz2.BZ2File(reportFile, 'w', 0)
+            else:
+                self._out = open(reportFile, 'w', 1)
+        else:
+            self._out = reportFile
+        
+        # Initialize the weights.
+        
+        if weights is None:
+            self._weights = [0.0]*numTemperatures
+            self._updateWeights = True
+            self._weightUpdateFactor = 1.0
+            self._histogram = [0]*numTemperatures
+            self._hasMadeTransition = False
+        else:
+            self._weights = weights
+            self._updateWeights = False
+
+        # Select the initial temperature.
+        
+        self.currentTemperature = 0
+        self.simulation.integrator.setTemperature(self.temperatures[self.currentTemperature])
+        
+        # Add a reporter to the simulation which will handle the updates and reports.
+        
+        class STReporter(object):
+            def __init__(self, st):
+                self.st = st
+
+            def describeNextReport(self, simulation):
+                st = self.st
+                steps1 = st.tempChangeInterval - simulation.currentStep%st.tempChangeInterval
+                steps2 = st.reportInterval - simulation.currentStep%st.reportInterval
+                steps = min(steps1, steps2)
+                isUpdateAttempt = (steps1 == steps)
+                return (steps, False, isUpdateAttempt, False, isUpdateAttempt)
+
+            def report(self, simulation, state):
+                st = self.st
+                if simulation.currentStep%st.tempChangeInterval == 0:
+                    st._attemptTemperatureChange(state)
+                if simulation.currentStep%st.reportInterval == 0:
+                    st._writeReport()
+        
+        simulation.reporters.append(STReporter(self))
+        
+        # Write out the header line.
+        
+        headers = ['Steps', 'Temperature (K)']
+        for t in self.temperatures:
+            headers.append('%gK Weight' % t.value_in_unit(unit.kelvin))
+        print('#"%s"' % ('"\t"').join(headers), file=self._out)
+
+    def __del__(self):
+        if self._openedFile:
+            self._out.close()
+    
+    @property
+    def weights(self):
+        return [x-self._weights[0] for x in self._weights]
+
+    def step(self, steps):
+        """Advance the simulation by integrating a specified number of time steps."""
+        self.simulation.step(steps)
+    
+    def _attemptTemperatureChange(self, state):
+        """Attempt to move to a different temperature."""
+        
+        # Compute the probability for each temperature.  This is done in log space to avoid overflow.
+
+        logProbability = [(self._weights[i]-self.inverseTemperatures[i]*state.getPotentialEnergy()) for i in range(len(self._weights))]
+        maxLogProb = max(logProbability)
+        offset = maxLogProb + math.log(sum(math.exp(x-maxLogProb) for x in logProbability))
+        probability = [math.exp(x-offset) for x in logProbability]
+        r = random.random()
+        for j in range(len(probability)):
+            if r < probability[j]:
+                if j != self.currentTemperature:
+                    # Rescale the velocities.
+                    
+                    scale = math.sqrt(self.temperatures[j]/self.temperatures[self.currentTemperature])
+                    if have_numpy:
+                        velocities = scale*state.getVelocities(asNumpy=True).value_in_unit(unit.nanometers/unit.picoseconds)
+                    else:
+                        velocities = [v*scale for v in state.getVelocities().value_in_unit(unit.nanometers/unit.picoseconds)]
+                    self.simulation.context.setVelocities(velocities)
+
+                    # Select this temperature.
+                    
+                    self._hasMadeTransition = True
+                    self.currentTemperature = j
+                    self.simulation.integrator.setTemperature(self.temperatures[j])
+                if self._updateWeights:
+                    # Update the weight factors.
+                    
+                    self._weights[j] -= self._weightUpdateFactor
+                    self._histogram[j] += 1
+                    minCounts = min(self._histogram)
+                    if minCounts > 20 and minCounts >= 0.2*sum(self._histogram)/len(self._histogram):
+                        # Reduce the weight update factor and reset the histogram.
+                        
+                        self._weightUpdateFactor *= 0.5
+                        self._histogram = [0]*len(self.temperatures)
+                        self._weights = [x-self._weights[0] for x in self._weights]
+                    elif not self._hasMadeTransition and probability[self.currentTemperature] > 0.99:
+                        # Rapidly increase the weight update factor at the start of the simulation to find
+                        # a reasonable starting value.
+                        
+                        self._weightUpdateFactor *= 2.0
+                        self._histogram = [0]*len(self.temperatures)
+                return
+            r -= probability[j]
+
+    def _writeReport(self):
+        """Write out a line to the report."""
+        temperature = self.temperatures[self.currentTemperature].value_in_unit(unit.kelvin)
+        values = [temperature]+self.weights
+        print(('%d\t' % self.simulation.currentStep) + '\t'.join('%g' % v for v in values), file=self._out)