Merge master

.travis.yml

@@ -7,11 +7,14 @@ addons:
       - doxygen
       - python-numpy
       - python-scipy
+      - libfftw3-dev
 
 matrix:
   include:
     - sudo: required
-      env: OPENCL=true
+      env: ==CPU_OPENCL==
+           OPENCL=true
+           CUDA=false
            CC=gcc
            CXX=g++
            CMAKE_FLAGS="
@@ -28,48 +31,104 @@ matrix:
            -DOPENMM_BUILD_EXAMPLES=OFF"
       addons: {apt: {packages: []}}
 
+    - sudo: required
+      env: ==CUDA_COMPILE==
+           CUDA=true
+           OPENCL=false
+           CUDA_VERSION="7.0-28"
+           CMAKE_FLAGS="
+             -DOPENMM_BUILD_CUDA_TESTS=OFF
+             -DOPENMM_BUILD_OPENCL_TESTS=OFF
+             -DOPENMM_BUILD_PYTHON_WRAPPERS=OFF
+             -DOPENMM_BUILD_REFERENCE_TESTS=OFF
+             -DOPENMM_BUILD_SERIALIZATION_TESTS=OFF
+             -DOPENMM_BUILD_C_AND_FORTRAN_WRAPPERS=OFF
+             -DOPENMM_BUILD_EXAMPLES=OFF
+             -DOPENCL_LIBRARY=/usr/local/cuda-7.0/lib64/libOpenCL.so"
+      addons: {apt: {packages: []}}
+
+    - language: objective-c
+      os: osx
+      env: ==OSX==
+           OPENCL=false
+           CUDA=false
+           CMAKE_FLAGS="
+             -DOPENMM_BUILD_OPENCL_TESTS=OFF
+             -DSWIG_EXECUTABLE=/usr/local/Cellar/swig/3.0.2/bin/swig"
+      addons: {apt: {packages: []}}
+
     - sudo: false
       python: 2.7_with_system_site_packages
-      env: OPENCL=false
+      env: ==STATIC_LIB==
+           OPENCL=false
+           CUDA=false
            CC=clang
            CXX=clang++
            CMAKE_FLAGS="-DOPENMM_BUILD_STATIC_LIB=ON"
 
     - sudo: false
       python: 2.7_with_system_site_packages
-      env: OPENCL=false
+      env: ==PYTNON_2==
+           OPENCL=false
+           CUDA=false
            CC=clang
            CXX=clang++
-           CMAKE_FLAGS="-DOPENMM_BUILD_STATIC_LIB=OFF"
+           CMAKE_FLAGS=""
 
     - sudo: false
       python: 3.4
-      env: OPENCL=false
+      env: ==PYTHON_3==
+           OPENCL=false
+           CUDA=false
            CC=gcc
            CXX=g++
-           CMAKE_FLAGS="-DOPENMM_BUILD_STATIC_LIB=OFF"
+           CMAKE_FLAGS=""
 
 before_install:
+  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then
+      brew install doxygen swig fftw;
+      sudo easy_install pytest;
+    fi
   - if [[ "$OPENCL" == "true" ]]; then
-      sudo apt-get -yq update &>> ~/apt-get-update.log;
+      sudo apt-get -yq update > /dev/null 2>&1 ;
       sudo apt-get install -qq fglrx=2:8.960-0ubuntu1 opencl-headers;
     fi
-  - if [[ "$OPENCL" == "false" ]]; then
-      wget https://anaconda.org/anaconda/swig/3.0.2/download/linux-64/swig-3.0.2-0.tar.bz2;
+  # Install swig for Python wrappers. However, testing CUDA and OpenCL, we
+  # skip the Python wrapper for speed. We're not using anaconda python,
+  # but this is a fast way to get an apparently functional precompiled
+  # build of swig that's more modern than what's in apt.
+  - if [[ "$OPENCL" == "false" && "$CUDA" == "false" && "$TRAVIS_OS_NAME" == "linux" ]]; then
+      wget https://anaconda.org/omnia/swig/3.0.7/download/linux-64/swig-3.0.7-0.tar.bz2;
       mkdir $HOME/swig;
-      tar -xjvf swig-3.0.2-0.tar.bz2 -C $HOME/swig;
+      tar -xjvf swig-3.0.7-0.tar.bz2 -C $HOME/swig;
       export PATH=$HOME/swig/bin:$PATH;
-      export SWIG_LIB=$HOME/swig/share/swig/3.0.2;
+      export SWIG_LIB=$HOME/swig/share/swig/3.0.7;
+    fi
+
+  - if [[ "$CUDA" == "true" ]]; then
+      wget "http://developer.download.nvidia.com/compute/cuda/repos/ubuntu1204/x86_64/cuda-repo-ubuntu1204_${CUDA_VERSION}_amd64.deb";
+      sudo dpkg -i cuda-repo-ubuntu1204_${CUDA_VERSION}_amd64.deb;
+      sudo apt-get update -qq;
+      export CUDA_APT=${CUDA_VERSION%-*};
+      export CUDA_APT=${CUDA_APT/./-};
+      sudo apt-get install -y cuda-drivers cuda-core-${CUDA_APT} cuda-cudart-dev-${CUDA_APT} cuda-cufft-dev-${CUDA_APT};
+      sudo apt-get clean;
+      export CUDA_HOME=/usr/local/cuda-${CUDA_VERSION%%-*};
+      export LD_LIBRARY_PATH=${CUDA_HOME}/lib64:${LD_LIBRARY_PATH};
+      export PATH=${CUDA_HOME}/bin:${PATH};
     fi
 
 script:
   - CTEST_STOP_TIME=$(python -c "from datetime import datetime, timedelta; import sys; sys.stdout.write((datetime.now() + timedelta(minutes=25)).strftime('%H:%M:%S'))")
   - cmake . $CMAKE_FLAGS -DCMAKE_INSTALL_PREFIX=$HOME/OpenMM
-
   - make -j2 install
   - if [[ "$OPENCL" == "true" ]]; then ./TestOpenCLDeviceQuery; fi
-  - if [[ "$OPENCL" == "false" ]]; then
-      make PythonInstall;
+  - if [[ "$OPENCL" == "false" && "$CUDA" == "false" ]]; then
+      if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then
+          sudo make PythonInstall;
+      else
+          make PythonInstall;
+      fi;
       python -m simtk.testInstallation;
       (cd python/tests && py.test -v *);
     fi

CMakeLists.txt

@@ -109,6 +109,7 @@ ELSE( CMAKE_SIZEOF_VOID_P EQUAL 8 )
   SET( LIB64  )
 ENDIF( CMAKE_SIZEOF_VOID_P EQUAL 8 )
 
+
 IF (APPLE AND (NOT PNACL))
     # Build 64 bit binaries compatible with OS X 10.7
     IF (NOT CMAKE_OSX_DEPLOYMENT_TARGET)
@@ -117,6 +118,14 @@ IF (APPLE AND (NOT PNACL))
     IF (NOT CMAKE_OSX_ARCHITECTURES)
         SET (CMAKE_OSX_ARCHITECTURES "x86_64" CACHE STRING "The processor architectures to build for" FORCE)
     ENDIF (NOT CMAKE_OSX_ARCHITECTURES)
+    IF (NOT CMAKE_OSX_SYSROOT)
+        EXECUTE_PROCESS(COMMAND "xcrun" "--show-sdk-path" OUTPUT_VARIABLE XCRUN_OSX_SYSROOT RESULT_VARIABLE XCRUN_OSX_SYSROOT_STATUS OUTPUT_STRIP_TRAILING_WHITESPACE)
+        IF (XCRUN_OSX_SYSROOT_STATUS EQUAL 0)
+            SET (CMAKE_OSX_SYSROOT "${XCRUN_OSX_SYSROOT}" CACHE STRING "SDK Path" FORCE)
+        ENDIF (XCRUN_OSX_SYSROOT_STATUS EQUAL 0)
+        UNSET(XCRUN_OSX_SYSROOT)
+        UNSET(XCRUN_OSX_SYSROOT_STATUS)
+    ENDIF (NOT CMAKE_OSX_SYSROOT)
 
     # Improve the linking behavior of Mac libraries
     SET (CMAKE_INSTALL_NAME_DIR "@rpath")

openmmapi/include/openmm/Context.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -165,6 +165,11 @@ public:
      * @param randomSeed     the random number seed to use when selecting velocities
      */
     void setVelocitiesToTemperature(double temperature, int randomSeed=osrngseed());
+    /**
+     * Get all adjustable parameters that have been defined by Force objects in the System, along
+     * with their current values.
+     */
+    const std::map<std::string, double>& getParameters() const;
     /**
      * Get the value of an adjustable parameter defined by a Force object in the System.
      * 

openmmapi/src/Context.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -114,21 +114,15 @@ State Context::getState(int types, bool enforcePeriodicBox, int groups) const {
                 center *= 1.0/molecules[i].size();
 
                 // Find the displacement to move it into the first periodic box.
-
-                int xcell = (int) floor(center[0]/periodicBoxSize[0][0]);
-                int ycell = (int) floor(center[1]/periodicBoxSize[1][1]);
-                int zcell = (int) floor(center[2]/periodicBoxSize[2][2]);
-                double dx = xcell*periodicBoxSize[0][0];
-                double dy = ycell*periodicBoxSize[1][1];
-                double dz = zcell*periodicBoxSize[2][2];
+                Vec3 diff;
+                diff -= periodicBoxSize[0]*static_cast<int>(center[0]/periodicBoxSize[0][0]);
+                diff -= periodicBoxSize[1]*static_cast<int>(center[1]/periodicBoxSize[1][1]);
+                diff -= periodicBoxSize[2]*static_cast<int>(center[2]/periodicBoxSize[2][2]);
 
                 // Translate all the particles in the molecule.
-                
                 for (int j = 0; j < (int) molecules[i].size(); j++) {
                     Vec3& pos = positions[molecules[i][j]];
-                    pos[0] -= dx;
-                    pos[1] -= dy;
-                    pos[2] -= dz;
+                    pos -= diff;
                 }
             }
         }
@@ -207,6 +201,10 @@ void Context::setVelocitiesToTemperature(double temperature, int randomSeed) {
     impl->applyVelocityConstraints(1e-5);
 }
 
+const map<string, double>& Context::getParameters() const {
+    return impl->getParameters();
+}
+
 double Context::getParameter(const string& name) const {
     return impl->getParameter(name);
 }

openmmapi/src/CustomNonbondedForce.cpp

@@ -263,9 +263,9 @@ void CustomNonbondedForce::setFunctionParameters(int index, const std::string& n
 
 int CustomNonbondedForce::addInteractionGroup(const std::set<int>& set1, const std::set<int>& set2) {
     for (set<int>::iterator it = set1.begin(); it != set1.end(); ++it)
-        ASSERT_VALID_INDEX(*it, particles);
+        ASSERT(*it >= 0);
     for (set<int>::iterator it = set2.begin(); it != set2.end(); ++it)
-        ASSERT_VALID_INDEX(*it, particles);
+        ASSERT(*it >= 0);
     interactionGroups.push_back(InteractionGroupInfo(set1, set2));
     return interactionGroups.size()-1;
 }

openmmapi/src/CustomNonbondedForceImpl.cpp

@@ -111,6 +111,26 @@ void CustomNonbondedForceImpl::initialize(ContextImpl& context) {
         if (cutoff > 0.5*boxVectors[0][0] || cutoff > 0.5*boxVectors[1][1] || cutoff > 0.5*boxVectors[2][2])
             throw OpenMMException("CustomNonbondedForce: The cutoff distance cannot be greater than half the periodic box size.");
     }
+    // Check that all interaction groups only specify particles that have been defined.
+    for (int group = 0; group < owner.getNumInteractionGroups(); group++) {
+        set<int> set1, set2;
+        owner.getInteractionGroupParameters(group, set1, set2);
+        for (set<int>::iterator it = set1.begin(); it != set1.end(); ++it)
+            if ((*it < 0) || (*it >= owner.getNumParticles())) {
+                stringstream msg;
+                msg << "CustomNonbondedForce: Interaction group " << group << " set1 contains a particle index (" << *it << ") "
+                    << "not present in system (" << owner.getNumParticles() << " particles).";
+                throw OpenMMException(msg.str());
+            }
+        for (set<int>::iterator it = set2.begin(); it != set2.end(); ++it)
+            if ((*it < 0) || (*it >= owner.getNumParticles())) {
+                stringstream msg;
+                msg << "CustomNonbondedForce: Interaction group " << group << " set2 contains a particle index (" << *it << ") "
+                    << "not present in system (" << owner.getNumParticles() << " particles).";
+                throw OpenMMException(msg.str());
+            }
+    }
+
     kernel.getAs<CalcCustomNonbondedForceKernel>().initialize(context.getSystem(), owner);
 }
 

plugins/rpmd/openmmapi/src/RPMDIntegrator.cpp

@@ -128,21 +128,15 @@ State RPMDIntegrator::getState(int copy, int types, bool enforcePeriodicBox, int
             center *= 1.0/molecules[i].size();
 
             // Find the displacement to move it into the first periodic box.
-
-            int xcell = (int) floor(center[0]/periodicBoxSize[0][0]);
-            int ycell = (int) floor(center[1]/periodicBoxSize[1][1]);
-            int zcell = (int) floor(center[2]/periodicBoxSize[2][2]);
-            double dx = xcell*periodicBoxSize[0][0];
-            double dy = ycell*periodicBoxSize[1][1];
-            double dz = zcell*periodicBoxSize[2][2];
+            Vec3 diff;
+            diff -= periodicBoxSize[0]*static_cast<int>(center[0]/periodicBoxSize[0][0]);
+            diff -= periodicBoxSize[1]*static_cast<int>(center[1]/periodicBoxSize[1][1]);
+            diff -= periodicBoxSize[2]*static_cast<int>(center[2]/periodicBoxSize[2][2]);
 
             // Translate all the particles in the molecule.
-
             for (int j = 0; j < (int) molecules[i].size(); j++) {
                 Vec3& pos = positions[molecules[i][j]];
-                pos[0] -= dx;
-                pos[1] -= dy;
-                pos[2] -= dz;
+                pos -= diff;
             }
         }
 
@@ -170,7 +164,7 @@ double RPMDIntegrator::computeKineticEnergy() {
 
 void RPMDIntegrator::step(int steps) {
     if (context == NULL)
-        throw OpenMMException("This Integrator is not bound to a context!");    
+        throw OpenMMException("This Integrator is not bound to a context!");
     if (!hasSetPosition) {
         // Initialize the positions from the context.
 

tests/TestEnforcePeriodicBox.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * 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) 2010-2015 Stanford University and the Authors.      *
+ * Authors: Robert McGibbon                                                   *
+ * 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/Platform.h"
+#include "openmm/VerletIntegrator.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+
+using namespace OpenMM;
+using namespace std;
+
+void testTruncatedOctahedron() {
+    const int numMolecules = 5;
+    const int numParticles = numMolecules*2;
+    const float cutoff = 2.0;
+    Vec3 a(6.7929, 0, 0);
+    Vec3 b(-2.264163559406279, 6.404455775962287, 0);
+    Vec3 c(-2.264163559406279, -3.2019384603140684, 5.54658849047036);
+        
+    System system;
+    system.setDefaultPeriodicBoxVectors(a, b, c);
+    NonbondedForce* force = new NonbondedForce();
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    
+    force->setCutoffDistance(cutoff);
+    force->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    
+    for (int i = 0; i < numMolecules; i++) {
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        force->addParticle(-1, 0.2, 0.2);
+        force->addParticle(1, 0.2, 0.2);
+        positions[2*i] = a*genrand_real2(sfmt) + b*genrand_real2(sfmt) + c*genrand_real2(sfmt);
+        positions[2*i+1] = positions[2*i] + Vec3(1.0, 0.0, 0.0);
+        system.addConstraint(2*i, 2*i+1, 1.0);
+    }
+    system.addForce(force);
+    
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, Platform::getPlatformByName("Reference"));
+    context.setPositions(positions);
+    State initialState = context.getState(State::Positions | State::Energy, true);
+    double initialEnergy = initialState.getPotentialEnergy();
+
+    context.setState(initialState);
+    State finalState = context.getState(State::Positions | State::Energy, true);
+    double finalEnergy = finalState.getPotentialEnergy();
+
+    ASSERT_EQUAL_TOL(initialEnergy, finalEnergy, 1e-4);
+}
+
+int main(int argc, char* argv[]) {
+    try {
+        testTruncatedOctahedron();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

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

@@ -256,14 +256,14 @@ class Modeller(object):
         self.topology = newTopology
         self.positions = newPositions
 
-    def addSolvent(self, forcefield, model='tip3p', boxSize=None, boxVectors=None, padding=None, numAdded=None, positiveIon='Na+', negativeIon='Cl-', ionicStrength=0*molar):
+    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.
 
         The algorithm works as follows:
 
         1. Water molecules are added to fill the box.
         2. Water molecules are removed if their distance to any solute atom is less than the sum of their van der Waals radii.
-        3. If the solute is charged, enough positive or negative ions are added to neutralize it.  Each ion is added by
+        3. If the solute is charged and neutralize=True, enough positive or negative ions are added to neutralize it.  Each ion is added by
            randomly selecting a water molecule and replacing it with the ion.
         4. Ion pairs are added to give the requested total ionic strength.
 
@@ -299,6 +299,8 @@ class Modeller(object):
         ionicStrength : concentration=0*molar
             the total concentration of ions (both positive and negative) to add.  This
             does not include ions that are added to neutralize the system.
+        neutralize : bool=True
+            whether to add ions to neutralize the system
         """
         if len([x for x in (boxSize, boxVectors, padding, numAdded) if x is not None]) > 1:
             raise ValueError('At most one of the following arguments may be specified: boxSize, boxVectors, padding, numAdded')
@@ -504,9 +506,6 @@ class Modeller(object):
 
         # Add ions to neutralize the system.
 
-        totalCharge = int(floor(0.5+sum((nonbonded.getParticleParameters(i)[0].value_in_unit(elementary_charge) for i in range(system.getNumParticles())))))
-        if abs(totalCharge) > len(addedWaters):
-            raise Exception('Cannot neutralize the system because the charge is greater than the number of available positions for ions')
         def addIon(element):
             # Replace a water by an ion.
             index = random.randint(0, len(addedWaters)-1)
@@ -514,8 +513,12 @@ class Modeller(object):
             newTopology.addAtom(element.symbol, element, newResidue)
             newPositions.append(addedWaters[index][1]*nanometer)
             del addedWaters[index]
-        for i in range(abs(totalCharge)):
-            addIon(positiveElement if totalCharge < 0 else negativeElement)
+        if neutralize:
+            totalCharge = int(floor(0.5+sum((nonbonded.getParticleParameters(i)[0].value_in_unit(elementary_charge) for i in range(system.getNumParticles())))))
+            if abs(totalCharge) > len(addedWaters):
+                raise Exception('Cannot neutralize the system because the charge is greater than the number of available positions for ions')
+            for i in range(abs(totalCharge)):
+                addIon(positiveElement if totalCharge < 0 else negativeElement)
 
         # Add ions based on the desired ionic strength.
 

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -308,6 +308,7 @@ UNITS = {
 ("AmoebaWcaDispersionForce",              "getShctd")                                      :  ( None, ()),
 
 ("Context", "getParameter") : (None, ()),
+("Context", "getParameters") : (None, ()),
 ("Context", "getMolecules") : (None, ()),
 ("CMAPTorsionForce", "getMapParameters") : (None, (None, 'unit.kilojoule_per_mole')),
 ("CMAPTorsionForce", "getTorsionParameters") : (None, ()),

wrappers/python/src/swig_doxygen/swig_lib/python/extend.i

+%inline %{
+    typedef int bitmask32t;
+%}
+
+%typemap(in) bitmask32t %{
+    $1 = 0;
+#if PY_VERSION_HEX >= 0x03000000
+    if (PyLong_Check($input)) {
+        unsigned long u = PyLong_AsUnsignedLongMask($input);
+#else
+    if (PyInt_Check($input)) {
+        unsigned long u = PyInt_AsUnsignedLongMask($input);
+#endif
+        // 64-bit Windows has 32-bit longs, but other platforms have
+        // 64-bit longs
+        $1 = u & 0xffffffff;
+    } else {
+        PyErr_SetString(PyExc_ValueError, "in method $symname, argument $argnum could not be converted to type $type");
+        SWIG_fail;
+    }
+%}
+
+
+
 %extend OpenMM::Context {
   PyObject *_getStateAsLists(int getPositions,
                             int getVelocities,
@@ -5,7 +29,7 @@
                             int getEnergy,
                             int getParameters,
                             int enforcePeriodic,
-                            int groups) {
+                            bitmask32t groups) {
     State state;
     PyThreadState* _savePythonThreadState = PyEval_SaveThread();
     int types = 0;
@@ -27,14 +51,9 @@
 
 
   %pythoncode %{
-    def getState(self,
-                 getPositions=False,
-                 getVelocities=False,
-                 getForces=False,
-                 getEnergy=False,
-                 getParameters=False,
-                 enforcePeriodicBox=False,
-                 groups=-1):
+    def getState(self, getPositions=False, getVelocities=False,
+                 getForces=False, getEnergy=False, getParameters=False,
+                 enforcePeriodicBox=False, groups=-1):
         """Get a State object recording the current state information stored in this context.
 
         Parameters
@@ -189,19 +208,7 @@ Parameters:
                  getParameters=False,
                  enforcePeriodicBox=False,
                  groups=-1):
-        """
-        getState(self,
-                 copy,
-                 getPositions = False,
-                 getVelocities = False,
-                 getForces = False,
-                 getEnergy = False,
-                 getParameters = False,
-                 enforcePeriodicBox = False,
-                 groups = -1)
-              -> State
-
-        Get a State object recording the current state information about one copy of the system.
+        """Get a State object recording the current state information about one copy of the system.
 
         Parameters
         ----------

wrappers/python/tests/TestForceGroups.py

@@ -20,7 +20,7 @@ class TestForceGroups(unittest.TestCase):
         self.context = context
 
     def test1(self):
-        n = 31  # Should be 32, but github issue #1198
+        n = 32
         for (i,j) in itertools.combinations(range(n), 2):
             groups = 1<<i | 1<<j
             e_0 = self.context.getState(getEnergy=True, groups=groups).getPotentialEnergy()._value
@@ -34,6 +34,12 @@ class TestForceGroups(unittest.TestCase):
             # groups must be an int or set
             self.context.getState(getEnergy=True, groups=(1, 2))
 
+    def test3(self):
+        e_0 = self.context.getState(getEnergy=True, groups=-1).getPotentialEnergy()._value
+        e_ref = sum(range(32))
+        self.assertEqual(e_0, e_ref)
+
+
 if __name__ == '__main__':
     unittest.main()
 

wrappers/python/tests/TestModeller.py

@@ -397,44 +397,46 @@ class TestModeller(unittest.TestCase):
 
     def test_addSolventNegativeSolvent(self):
         """ Test the addSolvent() method; test adding ions to a negatively charged solvent. """
-    
+
         topology_start = self.pdb.topology
         topology_start.setUnitCellDimensions(Vec3(3.5, 3.5, 3.5)*nanometers)
-               
-        # set up modeller with no solvent
-        modeller = Modeller(topology_start, self.positions)
-        modeller.deleteWater()
 
-        # add 5 Cl- ions to the original topology
-        topology_toAdd = Topology()
-        newChain = topology_toAdd.addChain()
-        for i in range(5):
-            topology_toAdd.addResidue('CL',  newChain)
-        residues = [residue for residue in topology_toAdd.residues()]
-        for i in range(5):
-            topology_toAdd.addAtom('Cl',Element.getBySymbol('Cl'), residues[i]) 
-        positions_toAdd = [Vec3(1.0,1.2,1.5), Vec3(1.7,1.0,1.4), Vec3(1.5,2.0,1.0),
-                           Vec3(2.0,2.0,2.0), Vec3(2.0,1.5,1.0)]*nanometers
-        modeller.add(topology_toAdd, positions_toAdd)
-        modeller.addSolvent(self.forcefield, ionicStrength=1.0*molar)
-        topology_after = modeller.getTopology()
+        for neutralize in (True, False):
+            # set up modeller with no solvent
+            modeller = Modeller(topology_start, self.positions)
+            modeller.deleteWater()
 
-        water_count = 0
-        sodium_count = 0
-        chlorine_count = 0
-        for residue in topology_after.residues():
-            if residue.name=='HOH':
-                water_count += 1
-            elif residue.name=='NA':
-                sodium_count += 1
-            elif residue.name=='CL':
-                chlorine_count += 1
+            # add 5 Cl- ions to the original topology
+            topology_toAdd = Topology()
+            newChain = topology_toAdd.addChain()
+            for i in range(5):
+                topology_toAdd.addResidue('CL',  newChain)
+            residues = [residue for residue in topology_toAdd.residues()]
+            for i in range(5):
+                topology_toAdd.addAtom('Cl',Element.getBySymbol('Cl'), residues[i]) 
+            positions_toAdd = [Vec3(1.0,1.2,1.5), Vec3(1.7,1.0,1.4), Vec3(1.5,2.0,1.0),
+                               Vec3(2.0,2.0,2.0), Vec3(2.0,1.5,1.0)]*nanometers
+            modeller.add(topology_toAdd, positions_toAdd)
+            modeller.addSolvent(self.forcefield, ionicStrength=1.0*molar, neutralize=neutralize)
+            topology_after = modeller.getTopology()
 
-        total_water_ions = water_count+sodium_count+chlorine_count
-        expected_ion_fraction = 1.0*molar/(55.4*molar)
-        expected_ions = math.floor((total_water_ions-10)*expected_ion_fraction/2+0.5)+5
-        self.assertEqual(sodium_count, expected_ions)
-        self.assertEqual(chlorine_count, expected_ions)
+            water_count = 0
+            sodium_count = 0
+            chlorine_count = 0
+            for residue in topology_after.residues():
+                if residue.name=='HOH':
+                    water_count += 1
+                elif residue.name=='NA':
+                    sodium_count += 1
+                elif residue.name=='CL':
+                    chlorine_count += 1
+
+            total_water_ions = water_count+sodium_count+chlorine_count
+            expected_ion_fraction = 1.0*molar/(55.4*molar)
+            expected_chlorine = math.floor((total_water_ions-10)*expected_ion_fraction/2+0.5)+5
+            expected_sodium = expected_chlorine if neutralize else expected_chlorine-5
+            self.assertEqual(sodium_count, expected_sodium)
+            self.assertEqual(chlorine_count, expected_chlorine)
 
     def test_addSolventPositiveSolvent(self):
         """ Test the addSolvent() method; test adding ions to a positively charged solvent. """
@@ -442,42 +444,44 @@ class TestModeller(unittest.TestCase):
         topology_start = self.pdb.topology
         topology_start.setUnitCellDimensions(Vec3(3.5, 3.5, 3.5)*nanometers)
 
-        # set up modeller with no solvent
-        modeller = Modeller(topology_start, self.positions)
-        modeller.deleteWater()
+        for neutralize in (True, False):
+            # set up modeller with no solvent
+            modeller = Modeller(topology_start, self.positions)
+            modeller.deleteWater()
 
-        # add 5 Na+ ions to the original topology
-        topology_toAdd = Topology()
-        newChain = topology_toAdd.addChain()
-        for i in range(5):
-            topology_toAdd.addResidue('NA', newChain)
-        residues = [residue for residue in topology_toAdd.residues()]
-        for i in range(5):
-             topology_toAdd.addAtom('Na',Element.getBySymbol('Na'), residues[i]) 
-        positions_toAdd = [Vec3(1.0,1.2,1.5), Vec3(1.7,1.0,1.4), Vec3(1.5,2.0,1.0),
-                           Vec3(2.0,2.0,2.0), Vec3(2.0,1.5,1.0)]*nanometers
-
-        # positions_toAdd doesn't need to change
-        modeller.add(topology_toAdd, positions_toAdd)
-        modeller.addSolvent(self.forcefield, ionicStrength=1.0*molar)
-        topology_after = modeller.getTopology()
+            # add 5 Na+ ions to the original topology
+            topology_toAdd = Topology()
+            newChain = topology_toAdd.addChain()
+            for i in range(5):
+                topology_toAdd.addResidue('NA', newChain)
+            residues = [residue for residue in topology_toAdd.residues()]
+            for i in range(5):
+                 topology_toAdd.addAtom('Na',Element.getBySymbol('Na'), residues[i]) 
+            positions_toAdd = [Vec3(1.0,1.2,1.5), Vec3(1.7,1.0,1.4), Vec3(1.5,2.0,1.0),
+                               Vec3(2.0,2.0,2.0), Vec3(2.0,1.5,1.0)]*nanometers
+
+            # positions_toAdd doesn't need to change
+            modeller.add(topology_toAdd, positions_toAdd)
+            modeller.addSolvent(self.forcefield, ionicStrength=1.0*molar, neutralize=neutralize)
+            topology_after = modeller.getTopology()
 
-        water_count = 0
-        sodium_count = 0
-        chlorine_count = 0
-        for residue in topology_after.residues():
-            if residue.name=='HOH':
-                water_count += 1
-            elif residue.name=='NA':
-                sodium_count += 1
-            elif residue.name=='CL':
-                chlorine_count += 1
+            water_count = 0
+            sodium_count = 0
+            chlorine_count = 0
+            for residue in topology_after.residues():
+                if residue.name=='HOH':
+                    water_count += 1
+                elif residue.name=='NA':
+                    sodium_count += 1
+                elif residue.name=='CL':
+                    chlorine_count += 1
 
-        total_water_ions = water_count+sodium_count+chlorine_count
-        expected_ion_fraction = 1.0*molar/(55.4*molar)
-        expected_ions = math.floor((total_water_ions-10)*expected_ion_fraction/2+0.5)+5
-        self.assertEqual(sodium_count, expected_ions)
-        self.assertEqual(chlorine_count, expected_ions)
+            total_water_ions = water_count+sodium_count+chlorine_count
+            expected_ion_fraction = 1.0*molar/(55.4*molar)
+            expected_sodium = math.floor((total_water_ions-10)*expected_ion_fraction/2+0.5)+5
+            expected_chlorine = expected_sodium if neutralize else expected_sodium-5
+            self.assertEqual(sodium_count, expected_sodium)
+            self.assertEqual(chlorine_count, expected_chlorine)
 
     def test_addSolventIons(self):
         """ Test the addSolvent() method with all possible choices for positive and negative ions. """