Merge remote-tracking branch 'origin/master' into charmm

tools/ci/requirements-conda.txt

+numpy
+cython
+pip
+nose
+

tools/conda-recipe/README.md

+This is a recipe for building the current development package into a conda binary.
+
+The installation on travis-ci is done by building the conda package, installing
+it, running the tests, and then if successful pushing the package to binstar
+(and the docs to AWS S3). The binstar auth token is an encrypted environment
+variable generated using:
+
+binstar auth -n openmm-travis -o omnia --max-age 22896000 -c --scopes api:write
+
+and then saved in the environment variable BINSTAR_TOKEN.
+
+You can set up travis to store an encrypted token via
+
+gem install travis travis encrypt BINSTAR_TOKEN=xx
+
+where xx is the token output by binstar. The final command should print a line (containing 'secure') for inclusion in your .travis.yml file.

tools/conda-recipe/build.sh

+#!/bin/bash
+
+CMAKE_FLAGS="-DCMAKE_INSTALL_PREFIX=$PREFIX"
+
+if [[ "$OSTYPE" == "linux-gnu" ]]; then
+    # setting the rpath so that libOpenMMPME.so finds the right libfftw3
+    CMAKE_FLAGS+=" -DCMAKE_INSTALL_RPATH=.."
+    CMAKE_FLAGS+=" -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++"
+    #CMAKE_FLAGS+=" -DOPENCL_LIBRARY=/opt/AMDAPP/lib/x86_64/libOpenCL.so" # TEST
+
+elif [[ "$OSTYPE" == "darwin"* ]]; then
+    export MACOSX_DEPLOYMENT_TARGET="10.7"
+    CMAKE_FLAGS+=" -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++"
+fi
+
+# Set location for FFTW3 on both linux and mac
+CMAKE_FLAGS+=" -DFFTW_INCLUDES=$PREFIX/include"
+if [[ "$OSTYPE" == "linux-gnu" ]]; then
+    CMAKE_FLAGS+=" -DFFTW_LIBRARY=$PREFIX/lib/libfftw3f.so"
+    CMAKE_FLAGS+=" -DFFTW_THREADS_LIBRARY=$PREFIX/lib/libfftw3f_threads.so"
+elif [[ "$OSTYPE" == "darwin"* ]]; then
+    CMAKE_FLAGS+=" -DFFTW_LIBRARY=$PREFIX/lib/libfftw3f.dylib"
+    CMAKE_FLAGS+=" -DFFTW_THREADS_LIBRARY=$PREFIX/lib/libfftw3f_threads.dylib"
+fi
+
+# Copy source to current directory.
+cp -r $RECIPE_DIR/../.. .
+
+# Build in subdirectory.
+mkdir build
+cd build
+cmake .. $CMAKE_FLAGS
+make -j4
+make install
+
+# Run C tests.
+# Exclude OpenCL tests because @peastman suspects mesa on travis implementation is broken.
+# @jchodera and @pgrinaway suspect travis is working, but AMD OpenCL tests are actually failing due to a bug.
+ctest -j2 -V -E "[A-Za-z]+OpenCL[A-Za-z]+"
+
+# Install Python wrappers.
+export OPENMM_INCLUDE_PATH=$PREFIX/include
+export OPENMM_LIB_PATH=$PREFIX/lib
+cd python
+$PYTHON setup.py install
+cd ..
+
+# Remove one random file
+#rm $PREFIX/bin/TestReferenceHarmonicBondForce
+
+# Copy all tests to bin directory so they will be distributed with install package.
+cp `find . -name "Test*" -type f -maxdepth 1` $PREFIX/bin
+

tools/conda-recipe/meta.yaml

+package:
+  name: openmm
+  version: !!str dev
+
+requirements:
+  build:
+    - cmake
+    - python
+    - fftw3f
+
+  run:
+    - python
+    - fftw3f
+    - numpy
+
+about:
+  home: https://simtk.org/home/openmm
+  license: GPL

tools/conda-recipe/plugin-dir.patch

+diff --git wrappers/python/simtk/openmm/__init__.py wrappers/python/simtk/openmm/__init__.py
+index 7e47b11..0ef290a 100644
+--- wrappers/python/simtk/openmm/__init__.py
++++ wrappers/python/simtk/openmm/__init__.py
+@@ -13,6 +13,10 @@ import os, os.path
+ from simtk.openmm.openmm import *
+ from simtk.openmm.vec3 import Vec3
+ from simtk.openmm import version
+-if os.getenv('OPENMM_PLUGIN_DIR') is None and os.path.isdir(version.openmm_library_path):
++_plugin_dir = os.path.abspath(os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', '..', '..', '..', 'plugins'))
++if os.path.isdir(_plugin_dir):
++    pluginLoadedLibNames = Platform.loadPluginsFromDirectory(_plugin_dir)
++    del _plugin_dir
++elif os.getenv('OPENMM_PLUGIN_DIR') is None and os.path.isdir(version.openmm_library_path):
+     pluginLoadedLibNames = Platform.loadPluginsFromDirectory(os.path.join(version.openmm_library_path, 'plugins'))
+ else:
+     pluginLoadedLibNames = Platform.loadPluginsFromDirectory(Platform.getDefaultPluginsDirectory())
\ No newline at end of file

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

@@ -13,7 +13,7 @@ Copyright (c) 2014 the Authors
 
 Author: Jason Deckman
 Contributors: Jason M. Swails
-Date: April 19, 2014
+Date: June 6, 2014
 
 Permission is hereby granted, free of charge, to any person obtaining a
 copy of this software and associated documentation files (the "Software"),
@@ -64,7 +64,7 @@ class CharmmCrdFile(object):
         self.resname = []                  # Residue name
         self.resid = []                    # Residue ID
         self.attype = []                   # Atom type
-        self.positions = []*u.angstroms    # 3N atomic coordinates
+        self.positions = []                # 3N atomic coordinates
         self.title = []                    # .crd file title block
         self.segid = []                    # Segment ID
         self.weighting = []                # Atom weighting
@@ -105,7 +105,7 @@ class CharmmCrdFile(object):
                 self.resname.append(line[2])
                 self.attype.append(line[3])
                 pos = Vec3(float(line[4]), float(line[5]), float(line[6]))
-                self.positions.append(pos * u.angstroms)
+                self.positions.append(pos)
                 self.segid.append(line[7])
                 self.resid.append(int(line[8]))
                 self.weighting.append(float(line[9]))
@@ -120,6 +120,10 @@ class CharmmCrdFile(object):
         except (ValueError, IndexError), e:
             raise CharmmFileError('Error parsing CHARMM coordinate file')
 
+        # Apply units to the positions now. Do it this way to allow for
+        # (possible) numpy functionality in the future.
+        self.positions = u.Quantity(self.positions, u.angstroms)
+
 class CharmmRstFile(object):
     """
     Reads and parses data, velocities and coordinates from a CHARMM restart
@@ -209,8 +213,9 @@ class CharmmRstFile(object):
         self.velocities = [v * ONE_TIMESCALE for v in self.velocities]
 
         # Add units to positions and velocities
-        self.positions *= u.angstroms
-        self.velocities *= u.angstroms / u.picoseconds
+        self.positions = u.Quantity(self.positions, u.angstroms)
+        self.positionsold = u.Quantity(self.positionsold, u.angstroms)
+        self.velocities = u.Quantity(self.velocities, u.angstroms/u.picoseconds)
 
     def _scan(self, handle, str, r=0): # read lines in file until str is found
         scanning = True

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

@@ -158,10 +158,13 @@ class CharmmPsfFile(object):
         title = list()
         for i in range(ntitle):
             title.append(psf.readline().rstrip())
-        # Skip the blank line
-        psf.readline()
+        # Skip all blank lines. Most of the time there is only 1, but I've seen
+        # a file that has 2
+        line = psf.readline().strip()
+        while not line:
+            line = psf.readline().strip()
         # Next is the number of atoms
-        natom = conv(psf.readline().strip(), int, 'natom')
+        natom = conv(line, int, 'natom')
         # Parse all of the atoms
         residue_list = ResidueList()
         atom_list = AtomList()
@@ -1182,8 +1185,11 @@ class CharmmPsfFile(object):
                 if switchDistance >= nonbondedCutoff:
                     raise ValueError('switchDistance is too large compared '
                                      'to the cutoff!')
-                    force.setUseSwitchingFunction(True)
-                    force.setSwitchingDistance(switchDistance)
+                if abs(switchDistance) != switchDistance:
+                    # Detects negatives for both Quantity and float
+                    raise ValueError('switchDistance must be non-negative!')
+                force.setUseSwitchingFunction(True)
+                force.setSwitchingDistance(switchDistance)
 
         else: # periodic
             # Set up box vectors (from inpcrd if available, or fall back to
@@ -1223,8 +1229,11 @@ class CharmmPsfFile(object):
                 if switchDistance >= nonbondedCutoff:
                     raise ValueError('switchDistance is too large compared '
                                      'to the cutoff!')
-                    force.setUseSwitchingFunction(True)
-                    force.setSwitchingDistance(switchDistance)
+                if abs(switchDistance) != switchDistance:
+                    # Detects negatives for both Quantity and float
+                    raise ValueError('switchDistance must be non-negative!')
+                force.setUseSwitchingFunction(True)
+                force.setSwitchingDistance(switchDistance)
 
             if ewaldErrorTolerance is not None:
                 force.setEwaldErrorTolerance(ewaldErrorTolerance)

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

-#!/bin/env python
-
 """
 element.py: Used for managing elements.
 

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

@@ -6,7 +6,7 @@ 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) 2012-2013 Stanford University and the Authors.
+Portions copyright (c) 2012-2014 Stanford University and the Authors.
 Authors: Peter Eastman
 Contributors:
 
@@ -111,7 +111,7 @@ class GromacsTopFile(object):
                 if len(fields) < 2:
                     raise ValueError('Illegal line in .top file: '+line)
                 name = fields[1]
-                valueStart = stripped.find(name, len(command))+len(name)
+                valueStart = stripped.find(name, len(command))+len(name)+1
                 value = line[valueStart:].strip()
                 self._defines[name] = value
             elif command == '#ifdef':
@@ -745,6 +745,12 @@ class GromacsTopFile(object):
                     atom1params = nb.getParticleParameters(baseAtomIndex+atoms[0])
                     atom2params = nb.getParticleParameters(baseAtomIndex+atoms[1])
                     exceptions.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], atom1params[0]*atom2params[0]*fudgeQQ, params[0], params[1]))
+                for fields in moleculeType.exclusions:
+                    atoms = [int(x)-1 for x in fields]
+                    for atom in atoms[1:]:
+                        if atom > atoms[0]:
+                            exceptions.append((baseAtomIndex+atoms[0], baseAtomIndex+atom, 0, 0, 0))
+                    
 
         # Create nonbonded exceptions.
 

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

@@ -12,9 +12,9 @@ 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) 2012-2013 Stanford University and the Authors.
+Portions copyright (c) 2012-2014 Stanford University and the Authors.
 Authors: Randall J. Radmer, John D. Chodera, Peter Eastman
-Contributors: Christoph Klein, Michael R. Shirts
+Contributors: Christoph Klein, Michael R. Shirts, Jason Swails
 
 Permission is hereby granted, free of charge, to any person obtaining a
 copy of this software and associated documentation files (the "Software"),
@@ -76,6 +76,10 @@ POINTER_LABEL_LIST = POINTER_LABELS.replace(',', '').split()
 VELSCALE = 20.455 # velocity conversion factor to angstroms/picosecond
 TINY = 1.0e-8
 
+class NbfixPresent(Exception):
+    """ Exception raised when NBFIX is used for the Lennard-Jones terms """
+    pass
+
 class PrmtopLoader(object):
     """Parsed AMBER prmtop file.
 
@@ -297,7 +301,11 @@ class PrmtopLoader(object):
         return self.residuePointerDict[iAtom]
 
     def getNonbondTerms(self):
-        """Return list of all rVdw, epsilon pairs for each atom"""
+        """
+        Return list of all rVdw, epsilon pairs for each atom. If off-diagonal
+        elements of the Lennard-Jones A and B coefficient matrices are found,
+        NbfixPresent exception is raised
+        """
         try:
             return self._nonbondTerms
         except AttributeError:
@@ -305,9 +313,10 @@ class PrmtopLoader(object):
         self._nonbondTerms=[]
         lengthConversionFactor = units.angstrom.conversion_factor_to(units.nanometer)
         energyConversionFactor = units.kilocalorie_per_mole.conversion_factor_to(units.kilojoule_per_mole)
+        numTypes = self.getNumTypes()
+        atomTypeIndexes=self._getAtomTypeIndexes()
+        type_parameters = [(0, 0) for i in range(numTypes)]
         for iAtom in range(self.getNumAtoms()):
-            numTypes=self.getNumTypes()
-            atomTypeIndexes=self._getAtomTypeIndexes()
             index=(numTypes+1)*(atomTypeIndexes[iAtom]-1)
             nbIndex=int(self._raw_data['NONBONDED_PARM_INDEX'][index])-1
             if nbIndex<0:
@@ -320,9 +329,31 @@ class PrmtopLoader(object):
             except ZeroDivisionError:
                 rMin = 1.0
                 epsilon = 0.0
+            type_parameters[atomTypeIndexes[iAtom]-1] = (rMin/2.0, epsilon)
             rVdw = rMin/2.0*lengthConversionFactor
             epsilon = epsilon*energyConversionFactor
             self._nonbondTerms.append( (rVdw, epsilon) )
+        # Check if we have any off-diagonal modified LJ terms that would require
+        # an NBFIX-like solution
+        for i in range(numTypes):
+            for j in range(numTypes):
+                index = int(self._raw_data['NONBONDED_PARM_INDEX'][numTypes*i+j]) - 1
+                rij = type_parameters[i][0] + type_parameters[j][0]
+                wdij = sqrt(type_parameters[i][1] * type_parameters[j][1])
+                a = float(self._raw_data['LENNARD_JONES_ACOEF'][index])
+                b = float(self._raw_data['LENNARD_JONES_BCOEF'][index])
+                if a == 0 or b == 0:
+                    if a != 0 or b != 0 or (wdij != 0 and rij != 0):
+                        del self._nonbondTerms
+                        raise NbfixPresent('Off-diagonal Lennard-Jones elements'
+                                           ' found. Cannot determine LJ '
+                                           'parameters for individual atoms.')
+                elif (abs((a - (wdij * rij ** 12)) / a) > 1e-6 or
+                      abs((b - (2 * wdij * rij**6)) / b) > 1e-6):
+                    del self._nonbondTerms
+                    raise NbfixPresent('Off-diagonal Lennard-Jones elements '
+                                       'found. Cannot determine LJ parameters '
+                                       'for individual atoms.')
         return self._nonbondTerms
 
     def _getBonds(self, bondPointers):
@@ -429,27 +460,67 @@ class PrmtopLoader(object):
                           +self._raw_data["DIHEDRALS_WITHOUT_HYDROGEN"]
         returnList=[]
         charges=self.getCharges()
-        nonbondTerms = self.getNonbondTerms()
-        for ii in range(0,len(dihedralPointers),5):
-            if int(dihedralPointers[ii+2])>0 and int(dihedralPointers[ii+3])>0:
-                iAtom = int(dihedralPointers[ii])//3
-                lAtom = int(dihedralPointers[ii+3])//3
-                iidx = int(dihedralPointers[ii+4]) - 1
-                chargeProd = charges[iAtom]*charges[lAtom]
-                (rVdwI, epsilonI) = nonbondTerms[iAtom]
-                (rVdwL, epsilonL) = nonbondTerms[lAtom]
-                rMin = (rVdwI+rVdwL)
-                epsilon = sqrt(epsilonI*epsilonL)
-                try:
-                    iScee = float(self._raw_data["SCEE_SCALE_FACTOR"][iidx])
-                except KeyError:
-                    iScee = 1.2
-                try:
-                    iScnb = float(self._raw_data["SCNB_SCALE_FACTOR"][iidx])
-                except KeyError:
-                    iScnb = 2.0
-
-                returnList.append((iAtom, lAtom, chargeProd, rMin, epsilon, iScee, iScnb))
+        try:
+            nonbondTerms = self.getNonbondTerms()
+        except NbfixPresent:
+            # We need to do the unit conversions here, since getNonbondTerms
+            # never finished and it has unit conversions in there
+            length_conv = units.angstrom.conversion_factor_to(units.nanometers)
+            ene_conv = units.kilocalories_per_mole.conversion_factor_to(
+                                units.kilojoules_per_mole)
+            parm_acoef = [float(x) for x in self._raw_data['LENNARD_JONES_ACOEF']]
+            parm_bcoef = [float(x) for x in self._raw_data['LENNARD_JONES_BCOEF']]
+            nbidx = [int(x) for x in self._raw_data['NONBONDED_PARM_INDEX']]
+            numTypes = self.getNumTypes()
+            atomTypeIndexes=self._getAtomTypeIndexes()
+            for ii in range(0, len(dihedralPointers), 5):
+                if int(dihedralPointers[ii+2])>0 and int(dihedralPointers[ii+3])>0:
+                    iAtom = int(dihedralPointers[ii])//3
+                    lAtom = int(dihedralPointers[ii+3])//3
+                    iidx = int(dihedralPointers[ii+4]) - 1
+                    chargeProd = charges[iAtom]*charges[lAtom]
+                    typ1 = atomTypeIndexes[iAtom] - 1
+                    typ2 = atomTypeIndexes[lAtom] - 1
+                    idx = nbidx[numTypes*typ1+typ2] - 1
+                    a = parm_acoef[idx]
+                    b = parm_bcoef[idx]
+                    try:
+                        epsilon = b * b / (4 * a) * ene_conv
+                        rMin = (2 * a / b) ** (1/6.0) * length_conv
+                    except ZeroDivisionError:
+                        rMin = 1
+                        epsilon = 0
+                    try:
+                        iScee = float(self._raw_data['SCEE_SCALE_FACTOR'][iidx])
+                    except KeyError:
+                        iScee = 1.2
+                    try:
+                        iScnb = float(self._raw_data['SCNB_SCALE_FACTOR'][iidx])
+                    except KeyError:
+                        iScnb = 2.0
+                    returnList.append((iAtom, lAtom, chargeProd, rMin, epsilon, iScee, iScnb))
+        else:
+            # This block gets hit if NbfixPresent is _not_ caught
+            for ii in range(0,len(dihedralPointers),5):
+                if int(dihedralPointers[ii+2])>0 and int(dihedralPointers[ii+3])>0:
+                    iAtom = int(dihedralPointers[ii])//3
+                    lAtom = int(dihedralPointers[ii+3])//3
+                    iidx = int(dihedralPointers[ii+4]) - 1
+                    chargeProd = charges[iAtom]*charges[lAtom]
+                    (rVdwI, epsilonI) = nonbondTerms[iAtom]
+                    (rVdwL, epsilonL) = nonbondTerms[lAtom]
+                    rMin = (rVdwI+rVdwL)
+                    epsilon = sqrt(epsilonI*epsilonL)
+                    try:
+                        iScee = float(self._raw_data["SCEE_SCALE_FACTOR"][iidx])
+                    except KeyError:
+                        iScee = 1.2
+                    try:
+                        iScnb = float(self._raw_data["SCNB_SCALE_FACTOR"][iidx])
+                    except KeyError:
+                        iScnb = 2.0
+    
+                    returnList.append((iAtom, lAtom, chargeProd, rMin, epsilon, iScee, iScnb))
         return returnList
 
     def getExcludedAtoms(self):
@@ -782,9 +853,54 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
 
     # Add per-particle nonbonded parameters.
     sigmaScale = 2**(-1./6.) * 2.0
-    for (charge, (rVdw, epsilon)) in zip(prmtop.getCharges(), prmtop.getNonbondTerms()):
-        sigma = rVdw * sigmaScale
-        force.addParticle(charge, sigma, epsilon)
+    nbfix = False
+    try:
+        nonbondTerms = prmtop.getNonbondTerms()
+    except NbfixPresent:
+        nbfix = True
+        for charge in prmtop.getCharges():
+            force.addParticle(charge, 1.0, 0.0)
+        numTypes = prmtop.getNumTypes()
+        parm_acoef = [float(x) for x in prmtop._raw_data['LENNARD_JONES_ACOEF']]
+        parm_bcoef = [float(x) for x in prmtop._raw_data['LENNARD_JONES_BCOEF']]
+        nbidx = [int(x) for x in prmtop._raw_data['NONBONDED_PARM_INDEX']]
+        acoef = [0 for i in range(numTypes*numTypes)]
+        bcoef = acoef[:] # copy
+        ene_conv = units.kilocalories_per_mole.conversion_factor_to(units.kilojoules_per_mole)
+        length_conv = units.angstroms.conversion_factor_to(units.nanometers)
+        afac = sqrt(ene_conv) * length_conv**6
+        bfac = ene_conv * length_conv**6
+        for i in range(numTypes):
+            for j in range(numTypes):
+                idx = nbidx[numTypes*i+j] - 1
+                acoef[i+numTypes*j] = sqrt(parm_acoef[idx]) * afac
+                bcoef[i+numTypes*j] = parm_bcoef[idx] * bfac
+        cforce = mm.CustomNonbondedForce('(a/r6)^2-b/r6; r6=r^6;'
+                                         'a=acoef(type1, type2);'
+                                         'b=bcoef(type1, type2);')
+        cforce.addTabulatedFunction('acoef',
+                    mm.Discrete2DFunction(numTypes, numTypes, acoef))
+        cforce.addTabulatedFunction('bcoef',
+                    mm.Discrete2DFunction(numTypes, numTypes, bcoef))
+        cforce.addPerParticleParameter('type')
+        for atom in prmtop._getAtomTypeIndexes():
+            cforce.addParticle((atom-1,))
+        # Now set the various properties based on the NonbondedForce object
+        if nonbondedMethod in ('PME', 'Ewald', 'CutoffPeriodic'):
+            cforce.setNonbondedMethod(cforce.CutoffPeriodic)
+            cforce.setCutoffDistance(nonbondedCutoff)
+        elif nonbondedMethod == 'CutoffNonPeriodic':
+            cforce.setNonbondedMethod(cforce.CutoffNonPeriodic)
+            cforce.setCutoffDistance(nonbondedCutoff)
+        elif nonbondedMethod == 'NoCutoff':
+            cforce.setNonbondedMethod(cforce.NoCutoff)
+        else:
+            raise ValueError('Unrecognized cutoff option %s' % nonbondedMethod)
+    else:
+        for (charge, (rVdw, epsilon)) in zip(prmtop.getCharges(), nonbondTerms):
+            sigma = rVdw * sigmaScale
+            force.addParticle(charge, sigma, epsilon)
+
 
     # Add 1-4 Interactions
     excludedAtomPairs = set()
@@ -807,6 +923,13 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
             if min((iAtom, jAtom), (jAtom, iAtom)) in excludedAtomPairs: continue
             force.addException(iAtom, jAtom, excludeParams[0], excludeParams[1], excludeParams[2])
 
+    # Copy the exceptions as exclusions to the CustomNonbondedForce if we have
+    # NBFIX terms
+    if nbfix:
+        for i in range(force.getNumExceptions()):
+            ii, jj, chg, sig, eps = force.getExceptionParameters(i)
+            cforce.addExclusion(ii, jj)
+        system.addForce(cforce)
     system.addForce(force)
 
     # Add virtual sites for water.

wrappers/python/simtk/openmm/app/internal/charmm/topologyobjects.py

@@ -12,7 +12,7 @@ Copyright (c) 2014 the Authors
 
 Author: Jason M. Swails
 Contributors:
-Date: April 18, 2014
+Date: July 3, 2014
 
 Permission is hereby granted, free of charge, to any person obtaining a
 copy of this software and associated documentation files (the "Software"),
@@ -421,8 +421,10 @@ class ResidueList(list):
         if self._last_residue is None:
             res = self._last_residue = Residue(resname, resnum)
             list.append(self, res)
-        elif self._last_residue != (resname, resnum):
-            if self._last_residue.idx == resnum:
+        elif (self._last_residue != (resname, resnum) or
+              system != self._last_residue.system):
+            if (self._last_residue.idx == resnum and
+                system == self._last_residue.system):
                 lresname = self._last_residue.resname
                 warnings.warn('Residue %d split into separate residues %s '
                               'and %s' % (resnum, lresname, resname),
@@ -499,7 +501,7 @@ class Angle(object):
         """ See if a bond or an atom is in this angle """
         if isinstance(thing, Bond):
             return self.atom2 in thing and (self.atom1 in thing or
-                                            self.atom2 in thing)
+                                            self.atom3 in thing)
         # Otherwise assume it's an atom
         return self.atom1 is thing or self.atom2 is thing or self.atom3 is thing
 

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

@@ -952,14 +952,14 @@ class Modeller(object):
                                     # This is a virtual site.  Compute its position by the correct rule.
 
                                     if site.type == 'average2':
-                                        position = site.weights[0]*templateAtomPositions[index+site.atoms[0]] + site.weights[1]*templateAtomPositions[index+site.atoms[1]]
+                                        position = site.weights[0]*templateAtomPositions[site.atoms[0]] + site.weights[1]*templateAtomPositions[site.atoms[1]]
                                     elif site.type == 'average3':
-                                        position = site.weights[0]*templateAtomPositions[index+site.atoms[0]] + site.weights[1]*templateAtomPositions[index+site.atoms[1]] + site.weights[2]*templateAtomPositions[index+site.atoms[2]]
+                                        position = site.weights[0]*templateAtomPositions[site.atoms[0]] + site.weights[1]*templateAtomPositions[site.atoms[1]] + site.weights[2]*templateAtomPositions[site.atoms[2]]
                                     elif site.type == 'outOfPlane':
-                                        v1 = templateAtomPositions[index+site.atoms[1]] - templateAtomPositions[index+site.atoms[0]]
-                                        v2 = templateAtomPositions[index+site.atoms[2]] - templateAtomPositions[index+site.atoms[0]]
+                                        v1 = templateAtomPositions[site.atoms[1]] - templateAtomPositions[site.atoms[0]]
+                                        v2 = templateAtomPositions[site.atoms[2]] - templateAtomPositions[site.atoms[0]]
                                         cross = Vec3(v1[1]*v2[2]-v1[2]*v2[1], v1[2]*v2[0]-v1[0]*v2[2], v1[0]*v2[1]-v1[1]*v2[0])
-                                        position = templateAtomPositions[index+site.atoms[0]] + site.weights[0]*v1 + site.weights[1]*v2 + site.weights[2]*cross
+                                        position = templateAtomPositions[site.atoms[0]] + site.weights[0]*v1 + site.weights[1]*v2 + site.weights[2]*cross
                             if position is None and atom.type in drudeTypeMap:
                                 # This is a Drude particle.  Put it on top of its parent atom.
 

wrappers/python/tests/TestAmberPrmtopFile.py

@@ -5,17 +5,14 @@ from simtk.openmm import *
 from simtk.unit import *
 import simtk.openmm.app.element as elem
 
+prmtop1 = AmberPrmtopFile('systems/alanine-dipeptide-explicit.prmtop')
+prmtop2 = AmberPrmtopFile('systems/alanine-dipeptide-implicit.prmtop')
+prmtop3 = AmberPrmtopFile('systems/ff14ipq.parm7')
+inpcrd3 = AmberInpcrdFile('systems/ff14ipq.rst7')
+
 class TestAmberPrmtopFile(unittest.TestCase):
 
     """Test the AmberPrmtopFile.createSystem() method."""
- 
-    def setUp(self):
-        """Set up the tests by loading the input files."""
-
-        # alanine dipeptide with explicit water
-        self.prmtop1 = AmberPrmtopFile('systems/alanine-dipeptide-explicit.prmtop')
-        # alanine dipeptide with implicit water
-        self.prmtop2 = AmberPrmtopFile('systems/alanine-dipeptide-implicit.prmtop')
 
     def test_NonbondedMethod(self):
         """Test all five options for the nonbondedMethod parameter."""
@@ -25,7 +22,7 @@ class TestAmberPrmtopFile(unittest.TestCase):
                      CutoffPeriodic:NonbondedForce.CutoffPeriodic, 
                      Ewald:NonbondedForce.Ewald, PME: NonbondedForce.PME}
         for method in methodMap:
-            system = self.prmtop1.createSystem(nonbondedMethod=method)
+            system = prmtop1.createSystem(nonbondedMethod=method)
             forces = system.getForces()
             self.assertTrue(any(isinstance(f, NonbondedForce) and 
                                 f.getNonbondedMethod()==methodMap[method] 
@@ -35,9 +32,9 @@ class TestAmberPrmtopFile(unittest.TestCase):
         """Test to make sure the nonbondedCutoff parameter is passed correctly."""
 
         for method in [CutoffNonPeriodic, CutoffPeriodic, Ewald, PME]:
-            system = self.prmtop1.createSystem(nonbondedMethod=method, 
-                                               nonbondedCutoff=2*nanometer, 
-                                               constraints=HBonds)
+            system = prmtop1.createSystem(nonbondedMethod=method, 
+                                          nonbondedCutoff=2*nanometer, 
+                                          constraints=HBonds)
             cutoff_distance = 0.0*nanometer
             cutoff_check = 2.0*nanometer
             for force in system.getForces():
@@ -49,9 +46,9 @@ class TestAmberPrmtopFile(unittest.TestCase):
         """Test to make sure the ewaldErrorTolerance parameter is passed correctly."""
 
         for method in [Ewald, PME]:
-            system = self.prmtop1.createSystem(nonbondedMethod=method, 
-                                               ewaldErrorTolerance=1e-6, 
-                                               constraints=HBonds)
+            system = prmtop1.createSystem(nonbondedMethod=method, 
+                                          ewaldErrorTolerance=1e-6, 
+                                          constraints=HBonds)
             tolerance = 0
             tolerance_check = 1e-6
             for force in system.getForces():
@@ -63,18 +60,18 @@ class TestAmberPrmtopFile(unittest.TestCase):
         """Test both options (True and False) for the removeCMMotion parameter."""
 
         for b in [True, False]:
-            system = self.prmtop1.createSystem(removeCMMotion=b)
+            system = prmtop1.createSystem(removeCMMotion=b)
             forces = system.getForces()
             self.assertEqual(any(isinstance(f, CMMotionRemover) for f in forces), b)
 
     def test_RigidWaterAndConstraints(self):
         """Test all eight options for the constraints and rigidWater parameters."""
 
-        topology = self.prmtop1.topology
+        topology = prmtop1.topology
         for constraints_value in [None, HBonds, AllBonds, HAngles]:
             for rigidWater_value in [True, False]:
-                system = self.prmtop1.createSystem(constraints=constraints_value, 
-                                                   rigidWater=rigidWater_value)
+                system = prmtop1.createSystem(constraints=constraints_value, 
+                                              rigidWater=rigidWater_value)
                 validateConstraints(self, topology, system, 
                                     constraints_value, rigidWater_value)
 
@@ -84,7 +81,7 @@ class TestAmberPrmtopFile(unittest.TestCase):
 
         """
         for implicitSolvent_value in [HCT, OBC1, OBC2, GBn]:
-            system = self.prmtop2.createSystem(implicitSolvent=implicitSolvent_value)
+            system = prmtop2.createSystem(implicitSolvent=implicitSolvent_value)
             forces = system.getForces()
             if implicitSolvent_value in set([HCT, OBC1, GBn]):
                 force_type = CustomGBForce
@@ -99,7 +96,7 @@ class TestAmberPrmtopFile(unittest.TestCase):
                      CutoffNonPeriodic:NonbondedForce.CutoffNonPeriodic}
         for implicitSolvent_value in [HCT, OBC1, OBC2, GBn]:
             for method in methodMap:
-                system = self.prmtop2.createSystem(implicitSolvent=implicitSolvent_value, 
+                system = prmtop2.createSystem(implicitSolvent=implicitSolvent_value, 
                                     solventDielectric=50.0, soluteDielectric=0.9, nonbondedMethod=method)
                 found_matching_solvent_dielectric=False
                 found_matching_solute_dielectric=False
@@ -136,10 +133,10 @@ class TestAmberPrmtopFile(unittest.TestCase):
     def test_HydrogenMass(self):
         """Test that altering the mass of hydrogens works correctly."""
         
-        topology = self.prmtop1.topology
+        topology = prmtop1.topology
         hydrogenMass = 4*amu
-        system1 = self.prmtop1.createSystem()
-        system2 = self.prmtop1.createSystem(hydrogenMass=hydrogenMass)
+        system1 = prmtop1.createSystem()
+        system2 = prmtop1.createSystem(hydrogenMass=hydrogenMass)
         for atom in topology.atoms():
             if atom.element == elem.hydrogen:
                 self.assertNotEqual(hydrogenMass, system1.getParticleMass(atom.index))
@@ -148,5 +145,35 @@ class TestAmberPrmtopFile(unittest.TestCase):
         totalMass2 = sum([system2.getParticleMass(i) for i in range(system2.getNumParticles())]).value_in_unit(amu)
         self.assertAlmostEqual(totalMass1, totalMass2)
 
+    def test_NBFIX(self):
+        """Test that prmtop files with modified off-diagonal LJ elements are treated properly"""
+        system = prmtop3.createSystem(nonbondedMethod=PME,
+                                      nonbondedCutoff=8*angstroms)
+        # Check the forces
+        has_nonbond_force = has_custom_nonbond_force = False
+        nonbond_exceptions = custom_nonbond_exclusions = 0
+        for force in system.getForces():
+            if isinstance(force, NonbondedForce):
+                has_nonbond_force = True
+                nonbond_exceptions = force.getNumExceptions()
+            elif isinstance(force, CustomNonbondedForce):
+                has_custom_nonbond_force = True
+                custom_nonbond_exceptions = force.getNumExclusions()
+        self.assertTrue(has_nonbond_force)
+        self.assertTrue(has_custom_nonbond_force)
+        self.assertEqual(nonbond_exceptions, custom_nonbond_exceptions)
+        integrator = VerletIntegrator(1.0*femtoseconds)
+        # Use reference platform, since it should always be present and
+        # 'working', and the system is plenty small so this won't be too slow
+        sim = Simulation(prmtop3.topology, system, integrator, Platform.getPlatformByName('Reference'))
+        # Check that the energy is about what we expect it to be
+        sim.context.setPeriodicBoxVectors(*inpcrd3.boxVectors)
+        sim.context.setPositions(inpcrd3.positions)
+        ene = sim.context.getState(getEnergy=True, enforcePeriodicBox=True).getPotentialEnergy()
+        ene = ene.value_in_unit(kilocalories_per_mole)
+        # Make sure the energy is relatively close to the value we get with
+        # Amber using this force field.
+        self.assertAlmostEqual(-7042.3903307/ene, 1, places=3)
+
 if __name__ == '__main__':
     unittest.main()

wrappers/python/tests/TestBytes.py

@@ -8,7 +8,8 @@ class TestBytes(unittest.TestCase):
         system.addParticle(1.0)
         refPositions = [(0,0,0)]
 
-        context = mm.Context(system, mm.VerletIntegrator(0))
+        platform = mm.Platform.getPlatformByName('Reference')
+        context = mm.Context(system, mm.VerletIntegrator(0), platform)
         context.setPositions(refPositions)
         chk = context.createCheckpoint()
         # check that the return value of createCheckpoint is of type bytes (non-unicode)
@@ -22,13 +23,6 @@ class TestBytes(unittest.TestCase):
 
         assert newPositions == refPositions
 
-        # try encoding the checkpoint in utf-8. OpenMM should be able to handle this too
-        context.setPositions([(12345, 12345, 123451)])
-        context.loadCheckpoint(chk.decode('utf-8'))
-        newPositions = context.getState(getPositions=True).getPositions()._value
-
-        assert newPositions == refPositions
-
 
 if __name__ == '__main__':
     unittest.main()

wrappers/python/tests/TestForceField.py

@@ -4,6 +4,7 @@ from simtk.openmm.app import *
 from simtk.openmm import *
 from simtk.unit import *
 import simtk.openmm.app.element as elem
+import simtk.openmm.app.forcefield as forcefield
 
 class TestForceField(unittest.TestCase):
     """Test the ForceField.createSystem() method."""
@@ -18,7 +19,8 @@ class TestForceField(unittest.TestCase):
         self.forcefield1 = ForceField('amber99sb.xml', 'tip3p.xml')
         self.topology1 = self.pdb1.topology
         self.topology1.setUnitCellDimensions(Vec3(2, 2, 2))
-        
+
+        # alalnine dipeptide with implicit water
         self.pdb2 = PDBFile('systems/alanine-dipeptide-implicit.pdb')
         self.forcefield2 = ForceField('amber99sb.xml', 'amber99_obc.xml')
 
@@ -38,6 +40,18 @@ class TestForceField(unittest.TestCase):
                                 f.getNonbondedMethod()==methodMap[method] 
                                 for f in forces))
 
+    def test_DispersionCorrection(self):
+        """Test to make sure the nonbondedCutoff parameter is passed correctly."""
+
+        for useDispersionCorrection in [True, False]:
+            system = self.forcefield1.createSystem(self.pdb1.topology,
+                                                   nonbondedCutoff=2*nanometer, 
+                                                   useDispersionCorrection=useDispersionCorrection)
+
+            for force in system.getForces():
+                if isinstance(force, NonbondedForce):
+                    self.assertEqual(useDispersionCorrection, force.getUseDispersionCorrection())
+
     def test_Cutoff(self):
         """Test to make sure the nonbondedCutoff parameter is passed correctly."""
 
@@ -119,6 +133,113 @@ class TestForceField(unittest.TestCase):
         totalMass1 = sum([system1.getParticleMass(i) for i in range(system1.getNumParticles())]).value_in_unit(amu)
         totalMass2 = sum([system2.getParticleMass(i) for i in range(system2.getNumParticles())]).value_in_unit(amu)
         self.assertAlmostEqual(totalMass1, totalMass2)
+    
+    def test_Forces(self):
+        """Compute forces and compare them to ones generated with a previous version of OpenMM to ensure they haven't changed."""
+        
+        pdb = PDBFile('systems/lysozyme-implicit.pdb')
+        system = self.forcefield2.createSystem(pdb.topology)
+        integrator = VerletIntegrator(0.001)
+        context = Context(system, integrator)
+        context.setPositions(pdb.positions)
+        state1 = context.getState(getForces=True)
+        state2 = XmlSerializer.deserialize(open('systems/lysozyme-implicit-forces.xml').read())
+        numDifferences = 0
+        for f1, f2, in zip(state1.getForces().value_in_unit(kilojoules_per_mole/nanometer), state2.getForces().value_in_unit(kilojoules_per_mole/nanometer)):
+            diff = norm(f1-f2)
+            if diff > 0.1 and diff/norm(f1) > 1e-3:
+                numDifferences += 1
+        self.assertTrue(numDifferences < system.getNumParticles()/20) # Tolerate occasional differences from numerical error
+    
+    def test_ProgrammaticForceField(self):
+        """Test building a ForceField programmatically."""
+        
+        # Build the ForceField for TIP3P programmatically.
+        ff = ForceField()
+        ff.registerAtomType({'name':'tip3p-O', 'class':'OW', 'mass':15.99943*daltons, 'element':elem.oxygen})
+        ff.registerAtomType({'name':'tip3p-H', 'class':'HW', 'mass':1.007947*daltons, 'element':elem.hydrogen})
+        residue = ForceField._TemplateData('HOH')
+        residue.atoms.append(ForceField._TemplateAtomData('O', 'tip3p-O', elem.oxygen))
+        residue.atoms.append(ForceField._TemplateAtomData('H1', 'tip3p-H', elem.hydrogen))
+        residue.atoms.append(ForceField._TemplateAtomData('H2', 'tip3p-H', elem.hydrogen))
+        residue.addBond(0, 1)
+        residue.addBond(0, 2)
+        ff.registerResidueTemplate(residue)
+        bonds = forcefield.HarmonicBondGenerator(ff)
+        bonds.registerBond({'class1':'OW', 'class2':'HW', 'length':0.09572*nanometers, 'k':462750.4*kilojoules_per_mole/nanometer})
+        ff.registerGenerator(bonds)
+        angles = forcefield.HarmonicAngleGenerator(ff)
+        angles.registerAngle({'class1':'HW', 'class2':'OW', 'class3':'HW', 'angle':1.82421813418*radians, 'k':836.8*kilojoules_per_mole/radian})
+        ff.registerGenerator(angles)
+        nonbonded = forcefield.NonbondedGenerator(ff, 0.833333, 0.5)
+        nonbonded.registerAtom({'type':'tip3p-O', 'charge':-0.834, 'sigma':0.31507524065751241*nanometers, 'epsilon':0.635968*kilojoules_per_mole})
+        nonbonded.registerAtom({'type':'tip3p-H', 'charge':0.417, 'sigma':1*nanometers, 'epsilon':0*kilojoules_per_mole})
+        ff.registerGenerator(nonbonded)
+        
+        # Build a water box.
+        modeller = Modeller(Topology(), [])
+        modeller.addSolvent(ff, boxSize=Vec3(3, 3, 3)*nanometers)
+        
+        # Create a system using the programmatic force field as well as one from an XML file.
+        system1 = ff.createSystem(modeller.topology)
+        ff2 = ForceField('tip3p.xml')
+        system2 = ff2.createSystem(modeller.topology)
+        self.assertEqual(XmlSerializer.serialize(system1), XmlSerializer.serialize(system2))
+
+class AmoebaTestForceField(unittest.TestCase):
+    """Test the ForceField.createSystem() method with the AMOEBA forcefield."""
+ 
+    def setUp(self):
+        """Set up the tests by loading the input pdb files and force field 
+        xml files.
+
+        """
+
+        self.pdb1 = PDBFile('systems/amoeba-ion-in-water.pdb')
+        self.forcefield1 = ForceField('amoeba2009.xml')
+        self.topology1 = self.pdb1.topology
+
+
+    def test_NonbondedMethod(self):
+        """Test all five options for the nonbondedMethod parameter."""
+
+        methodMap = {NoCutoff:AmoebaMultipoleForce.NoCutoff,
+                     PME:AmoebaMultipoleForce.PME}
+
+        for method in methodMap:
+            system = self.forcefield1.createSystem(self.pdb1.topology,
+                                                  nonbondedMethod=method)
+            forces = system.getForces()
+            self.assertTrue(any(isinstance(f, AmoebaMultipoleForce) and
+                                f.getNonbondedMethod()==methodMap[method]
+                                for f in forces))
+    def test_Cutoff(self):
+        """Test to make sure the nonbondedCutoff parameter is passed correctly."""
+
+        cutoff_distance = 0.7*nanometer
+        for method in [NoCutoff, PME]:
+            system = self.forcefield1.createSystem(self.pdb1.topology,
+                                                   nonbondedMethod=method,
+                                                   nonbondedCutoff=cutoff_distance,
+                                                   constraints=None)
+
+            for force in system.getForces():
+                if isinstance(force, AmoebaVdwForce):
+                    self.assertEqual(force.getCutoff(), cutoff_distance)
+                if isinstance(force, AmoebaMultipoleForce):
+                    self.assertEqual(force.getCutoffDistance(), cutoff_distance)
+
+    def test_DispersionCorrection(self):
+        """Test to make sure the nonbondedCutoff parameter is passed correctly."""
+
+        for useDispersionCorrection in [True, False]:
+            system = self.forcefield1.createSystem(self.pdb1.topology,
+                                                   nonbondedMethod=PME,
+                                                   useDispersionCorrection=useDispersionCorrection)
+
+            for force in system.getForces():
+                if isinstance(force, AmoebaVdwForce):
+                    self.assertEqual(useDispersionCorrection, force.getUseDispersionCorrection())
 
 
 if __name__ == '__main__':