Merge branch 'master' into gayberne

tests/TestCustomTorsionForce.h

@@ -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-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -186,6 +186,42 @@ void testIllegalVariable() {
     ASSERT(threwException);
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions.
+    
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
+    VerletIntegrator integrator(0.01);
+    CustomTorsionForce* torsions = new CustomTorsionForce("k*(1+cos(n*theta-theta0))");
+    torsions->addPerTorsionParameter("theta0");
+    torsions->addPerTorsionParameter("n");
+    torsions->addGlobalParameter("k", 1.1);
+    vector<double> parameters(2);
+    parameters[0] = M_PI/3;
+    parameters[1] = 2;
+    torsions->addTorsion(0, 1, 2, 3, parameters);
+    torsions->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(torsions);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(4);
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    positions[3] = Vec3(1, 0, 2);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    double torque = -2*1.1*std::sin(2*M_PI/3);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, -torque, 0), forces[3], TOL);
+    ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
+    ASSERT_EQUAL_TOL(1.1*(1+std::cos(2*M_PI/3)), state.getPotentialEnergy(), TOL);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -194,6 +230,7 @@ int main(int argc, char* argv[]) {
         testTorsions();
         testRange();
         testIllegalVariable();
+        testPeriodic();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestHarmonicAngleForce.h

@@ -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-2015 Stanford University and the Authors.s      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.s      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -93,12 +93,40 @@ void testAngles() {
     }
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions.
+    
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 1.5, 0), Vec3(0, 0, 3));
+    VerletIntegrator integrator(0.01);
+    HarmonicAngleForce* angles = new HarmonicAngleForce();
+    angles->addAngle(0, 1, 2, PI_M/3, 1.1);
+    system.addForce(angles);
+    angles->setUsesPeriodicBoundaryConditions(true);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    double torque = 1.1*PI_M/6;
+    ASSERT_EQUAL_VEC(Vec3(2*torque, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, -torque, 0), forces[2], TOL);
+    ASSERT_EQUAL_TOL(0.5*1.1*(PI_M/6)*(PI_M/6), state.getPotentialEnergy(), TOL);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
     try {
         initializeTests(argc, argv);
         testAngles();
+        testPeriodic();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestHarmonicBondForce.h

@@ -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-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -82,12 +82,37 @@ void testBonds() {
     }
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions.
+    
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
+    VerletIntegrator integrator(0.01);
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    bonds->addBond(0, 1, 1.2, 0.8);
+    bonds->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(bonds);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 2, 0);
+    positions[1] = Vec3(0, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    ASSERT_EQUAL_VEC(Vec3(0, -0.8*0.2, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, 0.8*0.2, 0), forces[1], TOL);
+    ASSERT_EQUAL_TOL(0.5*0.8*0.2*0.2, state.getPotentialEnergy(), TOL);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
     try {
         initializeTests(argc, argv);
         testBonds();
+        testPeriodic();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestMonteCarloAnisotropicBarostat.h

@@ -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-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman, Lee-Ping Wang                                      *
  * Contributors:                                                              *
  *                                                                            *
@@ -75,7 +75,7 @@ void testIdealGas() {
     // Test it for three different temperatures.
     
     for (int i = 0; i < 3; i++) {
-        barostat->setTemperature(temp[i]);
+        barostat->setDefaultTemperature(temp[i]);
         LangevinIntegrator integrator(temp[i], 0.1, 0.01);
         Context context(system, integrator, platform);
         context.setPositions(positions);
@@ -135,7 +135,7 @@ void testIdealGasAxis(int axis) {
     // Test it for three different temperatures.
     
     for (int i = 0; i < 3; i++) {
-        barostat->setTemperature(temp[i]);
+        barostat->setDefaultTemperature(temp[i]);
         LangevinIntegrator integrator(temp[i], 0.1, 0.01);
         Context context(system, integrator, platform);
         context.setPositions(positions);
@@ -371,7 +371,7 @@ void testEinsteinCrystal() {
             // Create the barostat.
             MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pres3[p], pres3[p], pres3[p]), temp, (a==0||a==3), (a==1||a==3), (a==2||a==3), frequency);
             system.addForce(barostat);
-            barostat->setTemperature(temp);
+            barostat->setDefaultTemperature(temp);
             LangevinIntegrator integrator(temp, 0.1, 0.01);
             Context context(system, integrator, platform);
             context.setPositions(positions);
@@ -417,7 +417,7 @@ void testEinsteinCrystal() {
         // Create the barostat.
         MonteCarloBarostat* barostat = new MonteCarloBarostat(pres3[p], temp, frequency);
         system.addForce(barostat);
-        barostat->setTemperature(temp);
+        barostat->setDefaultTemperature(temp);
         LangevinIntegrator integrator(temp, 0.1, 0.001);
         Context context(system, integrator, platform);
         context.setPositions(positions);

tests/TestMonteCarloBarostat.h

@@ -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-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -111,7 +111,7 @@ void testIdealGas() {
     // Test it for three different temperatures.
 
     for (int i = 0; i < 3; i++) {
-        barostat->setTemperature(temp[i]);
+        barostat->setDefaultTemperature(temp[i]);
         LangevinIntegrator integrator(temp[i], 0.1, 0.01);
         Context context(system, integrator, platform);
         context.setPositions(positions);

tests/TestPeriodicTorsionForce.h

@@ -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-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -88,12 +88,43 @@ void testPeriodicTorsions() {
     }
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions.
+    
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
+    VerletIntegrator integrator(0.01);
+    PeriodicTorsionForce* torsions = new PeriodicTorsionForce();
+    torsions->addTorsion(0, 1, 2, 3, 2, PI_M/3, 1.1);
+    torsions->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(torsions);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(4);
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    positions[3] = Vec3(1, 0, 2);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    double torque = -2*1.1*std::sin(2*PI_M/3);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, -torque, 0), forces[3], TOL);
+    ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
+    ASSERT_EQUAL_TOL(1.1*(1+std::cos(2*PI_M/3)), state.getPotentialEnergy(), TOL);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
     try {
         initializeTests(argc, argv);
         testPeriodicTorsions();
+        testPeriodic();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestRBTorsionForce.h

@@ -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-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -107,12 +107,53 @@ void testRBTorsions() {
     }
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions.
+    
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
+    VerletIntegrator integrator(0.01);
+    RBTorsionForce* torsions = new RBTorsionForce();
+    torsions->addTorsion(0, 1, 2, 3, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6);
+    torsions->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(torsions);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(4);
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    positions[3] = Vec3(0, 0, 2);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    double psi = 0.5*PI_M;
+    double torque = 0.0;
+    for (int i = 1; i < 6; ++i) {
+        double c = 0.1*(i+1);
+        torque += -c*i*std::pow(std::cos(psi), i-1)*std::sin(psi);
+    }
+    ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, -torque, 0), forces[3], TOL);
+    ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
+    double energy = 0.0;
+    for (int i = 0; i < 6; ++i) {
+        double c = 0.1*(i+1);
+        energy += c*std::pow(std::cos(psi), i);
+    }
+    ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
     try {
         initializeTests(argc, argv);
         testRBTorsions();
+        testPeriodic();
         runPlatformTests();
     }
     catch(const exception& e) {

wrappers/generateWrappers.py

@@ -397,6 +397,7 @@ extern OPENMM_EXPORT void %(name)s_insert(%(name)s* set, %(type)s value);""" % v
 /* These methods need to be handled specially, since their C++ APIs cannot be directly translated to C.
    Unlike the C++ versions, the return value is allocated on the heap, and you must delete it yourself. */
 extern OPENMM_EXPORT OpenMM_State* OpenMM_Context_getState(const OpenMM_Context* target, int types, int enforcePeriodicBox);
+extern OPENMM_EXPORT OpenMM_State* OpenMM_Context_getState_2(const OpenMM_Context* target, int types, int enforcePeriodicBox, int groups);
 extern OPENMM_EXPORT OpenMM_StringArray* OpenMM_Platform_loadPluginsFromDirectory(const char* directory);
 extern OPENMM_EXPORT OpenMM_StringArray* OpenMM_Platform_getPluginLoadFailures();
 extern OPENMM_EXPORT char* OpenMM_XmlSerializer_serializeSystem(const OpenMM_System* system);
@@ -801,6 +802,10 @@ OPENMM_EXPORT OpenMM_State* OpenMM_Context_getState(const OpenMM_Context* target
     State result = reinterpret_cast<const Context*>(target)->getState(types, enforcePeriodicBox);
     return reinterpret_cast<OpenMM_State*>(new State(result));
 }
+OPENMM_EXPORT OpenMM_State* OpenMM_Context_getState_2(const OpenMM_Context* target, int types, int enforcePeriodicBox, int groups) {
+    State result = reinterpret_cast<const Context*>(target)->getState(types, enforcePeriodicBox, groups);
+    return reinterpret_cast<OpenMM_State*>(new State(result));
+}
 OPENMM_EXPORT OpenMM_StringArray* OpenMM_Platform_loadPluginsFromDirectory(const char* directory) {
     vector<string> result = Platform::loadPluginsFromDirectory(string(directory));
     return reinterpret_cast<OpenMM_StringArray*>(new vector<string>(result));
@@ -1314,6 +1319,14 @@ MODULE OpenMM
             integer*4 enforcePeriodicBox
             type(OpenMM_State) result
         end subroutine
+        subroutine OpenMM_Context_getState_2(target, types, enforcePeriodicBox, groups, result)
+            use OpenMM_Types; implicit none
+            type (OpenMM_Context) target
+            integer*4 types
+            integer*4 enforcePeriodicBox
+            integer*4 groups
+            type(OpenMM_State) result
+        end subroutine
         subroutine OpenMM_Platform_loadPluginsFromDirectory(directory, result)
             use OpenMM_Types; implicit none
             character(*) directory
@@ -1991,6 +2004,12 @@ OPENMM_EXPORT void openmm_context_getstate_(const OpenMM_Context*& target, int c
 OPENMM_EXPORT void OPENMM_CONTEXT_GETSTATE(const OpenMM_Context*& target, int const& types, int const& enforcePeriodicBox, OpenMM_State*& result) {
     result = OpenMM_Context_getState(target, types, enforcePeriodicBox);
 }
+OPENMM_EXPORT void openmm_context_getstate_2_(const OpenMM_Context*& target, int const& types, int const& enforcePeriodicBox, int const& groups, OpenMM_State*& result) {
+    result = OpenMM_Context_getState_2(target, types, enforcePeriodicBox, groups);
+}
+OPENMM_EXPORT void OPENMM_CONTEXT_GETSTATE_2(const OpenMM_Context*& target, int const& types, int const& enforcePeriodicBox, int const& groups, OpenMM_State*& result) {
+    result = OpenMM_Context_getState_2(target, types, enforcePeriodicBox, groups);
+}
 OPENMM_EXPORT void openmm_platform_loadpluginsfromdirectory_(const char* directory, OpenMM_StringArray*& result, int length) {
     result = OpenMM_Platform_loadPluginsFromDirectory(makeString(directory, length).c_str());
 }

wrappers/python/CMakeLists.txt

@@ -108,17 +108,10 @@ else(SWIG_EXECUTABLE)
     set(SWIG_VERSION "0.0.0" CACHE STRING "Swig version" FORCE)
 endif(SWIG_EXECUTABLE)
 # Enforce swig version
-if (PYTHON_VERSION_MAJOR EQUAL 3)
-    string(COMPARE LESS "${SWIG_VERSION}" "2.0.4" SWIG_VERSION_ERROR)
-    if(SWIG_VERSION_ERROR)
-        message("Swig version must be 2.0.4 or greater for Python 3! (You have ${SWIG_VERSION})")
-    endif(SWIG_VERSION_ERROR)
-else(PYTHON_VERSION_MAJOR EQUAL 3)
-    string(COMPARE LESS "${SWIG_VERSION}" "1.3.39" SWIG_VERSION_ERROR)
-    if(SWIG_VERSION_ERROR)
-        message("Swig version must be 1.3.39 or greater for Python 2! (You have ${SWIG_VERSION})")
-    endif(SWIG_VERSION_ERROR)
-endif(PYTHON_VERSION_MAJOR EQUAL 3)
+string(COMPARE LESS "${SWIG_VERSION}" "3.0.5" SWIG_VERSION_ERROR)
+if(SWIG_VERSION_ERROR)
+    message("Swig version must be 3.0.5 or greater! (You have ${SWIG_VERSION})")
+endif(SWIG_VERSION_ERROR)
 
 find_package(Doxygen REQUIRED)
 mark_as_advanced(CLEAR DOXYGEN_EXECUTABLE)

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

@@ -8,11 +8,10 @@ Structures at Stanford, funded under the NIH Roadmap for Medical Research,
 grant U54 GM072970. See https://simtk.org.  This code was originally part of
 the ParmEd program and was ported for use with OpenMM.
 
-Copyright (c) 2014-2015 the Authors
+Copyright (c) 2014-2016 the Authors
 
 Author: Jason M. Swails
 Contributors:
-Date: August 19, 2014
 
 Permission is hereby granted, free of charge, to any person obtaining a
 copy of this software and associated documentation files (the "Software"),
@@ -111,7 +110,7 @@ def _strip_optunit(thing, unit):
 
 # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
-_resre = re.compile(r'(\d+)([a-zA-Z]*)')
+_resre = re.compile(r'(-?\d+)([a-zA-Z]*)')
 
 class CharmmPsfFile(object):
     """A chemical structure instantiated from CHARMM files.

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

@@ -120,88 +120,101 @@ class ForceField(object):
         self._forces = []
         self._scripts = []
         self._templateGenerators = []
-        for file in files:
-            self.loadFile(file)
+        self.loadFile(files)
 
-    def loadFile(self, file):
+    def loadFile(self, files):
         """Load an XML file and add the definitions from it to this ForceField.
 
         Parameters
         ----------
-        file : string or file
-            An XML file containing force field definitions.  It may be either an
-            absolute file path, a path relative to the current working
+        files : string or file or tuple
+            An XML file or tuple of XML files containing force field definitions.
+            Each entry may be either an absolute file path, a path relative to the current working
             directory, a path relative to this module's data subdirectory (for
             built in force fields), or an open file-like object with a read()
             method from which the forcefield XML data can be loaded.
         """
-        try:
-            # this handles either filenames or open file-like objects
-            tree = etree.parse(file)
-        except IOError:
-            tree = etree.parse(os.path.join(os.path.dirname(__file__), 'data', file))
-        except Exception as e:
-            # Fail with an error message about which file could not be read.
-            # TODO: Also handle case where fallback to 'data' directory encounters problems,
-            # but this is much less worrisome because we control those files.
-            msg  = str(e) + '\n'
-            if hasattr(file, 'name'):
-                filename = file.name
-            else:
-                filename = str(file)
-            msg += "ForceField.loadFile() encountered an error reading file '%s'\n" % filename
-            raise Exception(msg)
 
-        root = tree.getroot()
+        if not isinstance(files, tuple):
+            files = (files,)
+
+        trees = []
+
+        for file in files:
+            try:
+                # this handles either filenames or open file-like objects
+                tree = etree.parse(file)
+            except IOError:
+                tree = etree.parse(os.path.join(os.path.dirname(__file__), 'data', file))
+            except Exception as e:
+                # Fail with an error message about which file could not be read.
+                # TODO: Also handle case where fallback to 'data' directory encounters problems,
+                # but this is much less worrisome because we control those files.
+                msg  = str(e) + '\n'
+                if hasattr(file, 'name'):
+                    filename = file.name
+                else:
+                    filename = str(file)
+                msg += "ForceField.loadFile() encountered an error reading file '%s'\n" % filename
+                raise Exception(msg)
+
+            trees.append(tree)
+
 
         # Load the atom types.
 
-        if tree.getroot().find('AtomTypes') is not None:
-            for type in tree.getroot().find('AtomTypes').findall('Type'):
-                self.registerAtomType(type.attrib)
+        for tree in trees:
+            if tree.getroot().find('AtomTypes') is not None:
+                for type in tree.getroot().find('AtomTypes').findall('Type'):
+                    self.registerAtomType(type.attrib)
 
         # Load the residue templates.
 
-        if tree.getroot().find('Residues') is not None:
-            for residue in root.find('Residues').findall('Residue'):
-                resName = residue.attrib['name']
-                template = ForceField._TemplateData(resName)
-                atomIndices = {}
-                for atom in residue.findall('Atom'):
-                    params = {}
-                    for key in atom.attrib:
-                        if key not in ('name', 'type'):
-                            params[key] = _convertParameterToNumber(atom.attrib[key])
-                    atomName = atom.attrib['name']
-                    if atomName in atomIndices:
-                        raise ValueError('Residue '+resName+' contains multiple atoms named '+atomName)
-                    atomIndices[atomName] = len(template.atoms)
-                    typeName = atom.attrib['type']
-                    template.atoms.append(ForceField._TemplateAtomData(atomName, typeName, self._atomTypes[typeName].element, params))
-                for site in residue.findall('VirtualSite'):
-                    template.virtualSites.append(ForceField._VirtualSiteData(site, atomIndices))
-                for bond in residue.findall('Bond'):
-                    if 'atomName1' in bond.attrib:
-                        template.addBondByName(bond.attrib['atomName1'], bond.attrib['atomName2'])
-                    else:
-                        template.addBond(int(bond.attrib['from']), int(bond.attrib['to']))
-                for bond in residue.findall('ExternalBond'):
-                    if 'atomName' in bond.attrib:
-                        template.addExternalBondByName(bond.attrib['atomName'])
-                    else:
-                        template.addExternalBond(int(bond.attrib['from']))
-                self.registerResidueTemplate(template)
+        for tree in trees:
+            if tree.getroot().find('Residues') is not None:
+                for residue in tree.getroot().find('Residues').findall('Residue'):
+                    resName = residue.attrib['name']
+                    template = ForceField._TemplateData(resName)
+                    if 'overload' in residue.attrib:
+                        template.overloadLevel = int(residue.attrib['overload'])
+                    atomIndices = {}
+                    for atom in residue.findall('Atom'):
+                        params = {}
+                        for key in atom.attrib:
+                            if key not in ('name', 'type'):
+                                params[key] = _convertParameterToNumber(atom.attrib[key])
+                        atomName = atom.attrib['name']
+                        if atomName in atomIndices:
+                            raise ValueError('Residue '+resName+' contains multiple atoms named '+atomName)
+                        atomIndices[atomName] = len(template.atoms)
+                        typeName = atom.attrib['type']
+                        template.atoms.append(ForceField._TemplateAtomData(atomName, typeName, self._atomTypes[typeName].element, params))
+                    for site in residue.findall('VirtualSite'):
+                        template.virtualSites.append(ForceField._VirtualSiteData(site, atomIndices))
+                    for bond in residue.findall('Bond'):
+                        if 'atomName1' in bond.attrib:
+                            template.addBondByName(bond.attrib['atomName1'], bond.attrib['atomName2'])
+                        else:
+                            template.addBond(int(bond.attrib['from']), int(bond.attrib['to']))
+                    for bond in residue.findall('ExternalBond'):
+                        if 'atomName' in bond.attrib:
+                            template.addExternalBondByName(bond.attrib['atomName'])
+                        else:
+                            template.addExternalBond(int(bond.attrib['from']))
+                    self.registerResidueTemplate(template)
 
         # Load force definitions
 
-        for child in root:
-            if child.tag in parsers:
-                parsers[child.tag](child, self)
+        for tree in trees:
+            for child in tree.getroot():
+                if child.tag in parsers:
+                    parsers[child.tag](child, self)
 
         # Load scripts
 
-        for node in tree.getroot().findall('Script'):
-            self.registerScript(node.text)
+        for tree in trees:
+            for node in tree.getroot().findall('Script'):
+                self.registerScript(node.text)
 
     def getGenerators(self):
         """Get the list of all registered generators."""
@@ -237,7 +250,25 @@ class ForceField(object):
         self._templates[template.name] = template
         signature = _createResidueSignature([atom.element for atom in template.atoms])
         if signature in self._templateSignatures:
-            self._templateSignatures[signature].append(template)
+            registered = False
+            for regtemplate in self._templateSignatures[signature]:
+                if regtemplate.name == template.name:
+                    if regtemplate.overloadLevel > template.overloadLevel:
+                        # ok to break - this is done every time a template is
+                        # registered so there can only be one already existing
+                        # with same name at a time
+                        registered = True
+                        break
+                    elif regtemplate.overloadLevel < template.overloadLevel:
+                        self._templateSignatures[signature].remove(regtemplate)
+                        self._templateSignatures[signature].append(template)
+                        registered = True
+                    else:
+                        raise Exception('Residue template %s with the same overloadLevel %d already exists.' %
+                                         (template.name, template.overloadLevel)
+                                         )
+            if not registered:
+                self._templateSignatures[signature].append(template)
         else:
             self._templateSignatures[signature] = [template]
 
@@ -379,6 +410,7 @@ class ForceField(object):
             self.virtualSites = []
             self.bonds = []
             self.externalBonds = []
+            self.overloadLevel = 0
 
         def getAtomIndexByName(self, atom_name):
             """Look up an atom index by atom name, providing a helpful error message if not found."""
@@ -566,11 +598,16 @@ class ForceField(object):
         matches = None
         signature = _createResidueSignature([atom.element for atom in res.atoms()])
         if signature in self._templateSignatures:
+            allMatches = []
             for t in self._templateSignatures[signature]:
-                matches = _matchResidue(res, t, bondedToAtom)
-                if matches is not None:
-                    template = t
-                    break
+                match = _matchResidue(res, t, bondedToAtom)
+                if match is not None:
+                    allMatches.append((t, match))
+            if len(allMatches) == 1:
+                template = allMatches[0][0]
+                matches = allMatches[0][1]
+            elif len(allMatches) > 1:
+                raise Exception('Multiple matching templates found for residue %d (%s).' % (res.index+1, res.name))
         return [template, matches]
 
     def _buildBondedToAtomList(self, topology):
@@ -703,7 +740,7 @@ class ForceField(object):
         return [templates, unique_unmatched_residues]
 
     def createSystem(self, topology, nonbondedMethod=NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
-                     constraints=None, rigidWater=True, removeCMMotion=True, hydrogenMass=None, **args):
+                     constraints=None, rigidWater=True, removeCMMotion=True, hydrogenMass=None, residueTemplates=dict(), **args):
         """Construct an OpenMM System representing a Topology with this force field.
 
         Parameters
@@ -727,6 +764,13 @@ class ForceField(object):
             The mass to use for hydrogen atoms bound to heavy atoms.  Any mass
             added to a hydrogen is subtracted from the heavy atom to keep
             their total mass the same.
+        residueTemplates : dict=dict()
+           Key: Topology Residue object
+           Value: string, name of _TemplateData residue template object to use for
+                  (Key) residue
+           This allows user to specify which template to apply to particular Residues
+           in the event that multiple matching templates are available (e.g Fe2+ and Fe3+
+           templates in the ForceField for a monoatomic iron ion in the topology).
         args
              Arbitrary additional keyword arguments may also be specified.
              This allows extra parameters to be specified that are specific to
@@ -765,8 +809,15 @@ class ForceField(object):
 
         for chain in topology.chains():
             for res in chain.residues():
-                # Attempt to match one of the existing templates.
-                [template, matches] = self._getResidueTemplateMatches(res, bondedToAtom)
+                if res in residueTemplates:
+                    tname = residueTemplates[res]
+                    template = self._templates[tname]
+                    matches = _matchResidue(res, template, bondedToAtom)
+                    if matches is None:
+                        raise Exception('User-supplied template %s does not match the residue %d (%s)' % (tname, res.index+1, res.name))
+                else:
+                    # Attempt to match one of the existing templates.
+                    [template, matches] = self._getResidueTemplateMatches(res, bondedToAtom)
                 if matches is None:
                     # No existing templates match.  Try any registered residue template generators.
                     for generator in self._templateGenerators:
@@ -856,13 +907,13 @@ class ForceField(object):
         uniquePropers = set()
         for angle in data.angles:
             for atom in bondedToAtom[angle[0]]:
-                if atom != angle[1]:
+                if atom not in angle:
                     if atom < angle[2]:
                         uniquePropers.add((atom, angle[0], angle[1], angle[2]))
                     else:
                         uniquePropers.add((angle[2], angle[1], angle[0], atom))
             for atom in bondedToAtom[angle[2]]:
-                if atom != angle[1]:
+                if atom not in angle:
                     if atom > angle[0]:
                         uniquePropers.add((angle[0], angle[1], angle[2], atom))
                     else:
@@ -1183,8 +1234,12 @@ class HarmonicBondGenerator(object):
 
     @staticmethod
     def parseElement(element, ff):
-        generator = HarmonicBondGenerator(ff)
-        ff.registerGenerator(generator)
+        existing = [f for f in ff._forces if isinstance(f, HarmonicBondGenerator)]
+        if len(existing) == 0:
+            generator = HarmonicBondGenerator(ff)
+            ff.registerGenerator(generator)
+        else:
+            generator = existing[0]
         for bond in element.findall('Bond'):
             generator.registerBond(bond.attrib)
 
@@ -1236,8 +1291,12 @@ class HarmonicAngleGenerator(object):
 
     @staticmethod
     def parseElement(element, ff):
-        generator = HarmonicAngleGenerator(ff)
-        ff.registerGenerator(generator)
+        existing = [f for f in ff._forces if isinstance(f, HarmonicAngleGenerator)]
+        if len(existing) == 0:
+            generator = HarmonicAngleGenerator(ff)
+            ff.registerGenerator(generator)
+        else:
+            generator = existing[0]
         for angle in element.findall('Angle'):
             generator.registerAngle(angle.attrib)
 
@@ -1320,8 +1379,12 @@ class PeriodicTorsionGenerator(object):
 
     @staticmethod
     def parseElement(element, ff):
-        generator = PeriodicTorsionGenerator(ff)
-        ff.registerGenerator(generator)
+        existing = [f for f in ff._forces if isinstance(f, PeriodicTorsionGenerator)]
+        if len(existing) == 0:
+            generator = PeriodicTorsionGenerator(ff)
+            ff.registerGenerator(generator)
+        else:
+            generator = existing[0]
         for torsion in element.findall('Proper'):
             generator.registerProperTorsion(torsion.attrib)
         for torsion in element.findall('Improper'):
@@ -1419,8 +1482,12 @@ class RBTorsionGenerator(object):
 
     @staticmethod
     def parseElement(element, ff):
-        generator = RBTorsionGenerator(ff)
-        ff.registerGenerator(generator)
+        existing = [f for f in ff._forces if isinstance(f, RBTorsionGenerator)]
+        if len(existing) == 0:
+            generator = RBTorsionGenerator(ff)
+            ff.registerGenerator(generator)
+        else:
+            generator = existing[0]
         for torsion in element.findall('Proper'):
             types = ff._findAtomTypes(torsion.attrib, 4)
             if None not in types:
@@ -1523,8 +1590,12 @@ class CMAPTorsionGenerator(object):
 
     @staticmethod
     def parseElement(element, ff):
-        generator = CMAPTorsionGenerator(ff)
-        ff.registerGenerator(generator)
+        existing = [f for f in ff._forces if isinstance(f, CMAPTorsionGenerator)]
+        if len(existing) == 0:
+            generator = CMAPTorsionGenerator(ff)
+            ff.registerGenerator(generator)
+        else:
+            generator = existing[0]
         for map in element.findall('Map'):
             values = [float(x) for x in map.text.split()]
             size = sqrt(len(values))

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

@@ -35,7 +35,7 @@ __author__ = "Peter Eastman"
 __version__ = "1.0"
 
 from simtk.openmm.app import Topology, PDBFile, ForceField
-from simtk.openmm.app.forcefield import HAngles, AllBonds, _createResidueSignature, _matchResidue, DrudeGenerator
+from simtk.openmm.app.forcefield import HAngles, AllBonds, CutoffNonPeriodic, _createResidueSignature, _matchResidue, DrudeGenerator
 from simtk.openmm.app.topology import Residue
 from simtk.openmm.vec3 import Vec3
 from simtk.openmm import System, Context, NonbondedForce, CustomNonbondedForce, HarmonicBondForce, HarmonicAngleForce, VerletIntegrator, LocalEnergyMinimizer
@@ -857,7 +857,7 @@ class Modeller(object):
         if forcefield is not None:
             # Use the ForceField the user specified.
 
-            system = forcefield.createSystem(newTopology, rigidWater=False)
+            system = forcefield.createSystem(newTopology, rigidWater=False, nonbondedMethod=CutoffNonPeriodic)
             atoms = list(newTopology.atoms())
             for i in range(system.getNumParticles()):
                 if atoms[i].element != elem.hydrogen:
@@ -869,6 +869,8 @@ class Modeller(object):
 
             system = System()
             nonbonded = CustomNonbondedForce('100/((r/0.1)^4+1)')
+            nonbonded.setNonbondedMethod(CustomNonbondedForce.CutoffNonPeriodic);
+            nonbonded.setCutoffDistance(1*nanometer)
             bonds = HarmonicBondForce()
             angles = HarmonicAngleForce()
             system.addForce(nonbonded)

wrappers/python/simtk/openmm/app/pdbfile.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-2015 Stanford University and the Authors.
+Portions copyright (c) 2012-2016 Stanford University and the Authors.
 Authors: Peter Eastman
 Contributors:
 
@@ -58,6 +58,9 @@ class PDBFile(object):
 
     _residueNameReplacements = {}
     _atomNameReplacements = {}
+    _standardResidues = ['ALA', 'ASN', 'CYS', 'GLU', 'HIS', 'LEU', 'MET', 'PRO', 'THR', 'TYR',
+                         'ARG', 'ASP', 'GLN', 'GLY', 'ILE', 'LYS', 'PHE', 'SER', 'TRP', 'VAL',
+                         'A', 'G', 'C', 'U', 'I', 'DA', 'DG', 'DC', 'DT', 'DI', 'HOH']
 
     def __init__(self, file, extraParticleIdentifier='EP'):
         """Load a PDB file.
@@ -75,10 +78,6 @@ class PDBFile(object):
         metalElements = ['Al','As','Ba','Ca','Cd','Ce','Co','Cs','Cu','Dy','Fe','Gd','Hg','Ho','In','Ir','K','Li','Mg',
         'Mn','Mo','Na','Ni','Pb','Pd','Pt','Rb','Rh','Sm','Sr','Te','Tl','V','W','Yb','Zn']
         
-        standardResidues = ['ALA', 'ASN', 'CYS', 'GLU', 'HIS', 'LEU', 'MET', 'PRO', 'THR', 'TYR',
-                            'ARG', 'ASP', 'GLN', 'GLY', 'ILE', 'LYS', 'PHE', 'SER', 'TRP', 'VAL',
-                            'A', 'G', 'C', 'U', 'I', 'DA', 'DG', 'DC', 'DT', 'DI', 'HOH']
-
         top = Topology()
         ## The Topology read from the PDB file
         self.topology = top
@@ -173,9 +172,9 @@ class PDBFile(object):
                     if atomByNumber[i].element is not None and atomByNumber[j].element is not None:
                         if atomByNumber[i].element.symbol not in metalElements and atomByNumber[j].element.symbol not in metalElements:
                             connectBonds.append((atomByNumber[i], atomByNumber[j])) 
-                        elif atomByNumber[i].element.symbol in metalElements and atomByNumber[j].residue.name not in standardResidues:
+                        elif atomByNumber[i].element.symbol in metalElements and atomByNumber[j].residue.name not in PDBFile._standardResidues:
                             connectBonds.append((atomByNumber[i], atomByNumber[j])) 
-                        elif atomByNumber[j].element.symbol in metalElements and atomByNumber[i].residue.name not in standardResidues:
+                        elif atomByNumber[j].element.symbol in metalElements and atomByNumber[i].residue.name not in PDBFile._standardResidues:
                             connectBonds.append((atomByNumber[i], atomByNumber[j]))     
                     else:
                         connectBonds.append((atomByNumber[i], atomByNumber[j]))         
@@ -331,6 +330,8 @@ class PDBFile(object):
             raise ValueError('Particle position is NaN')
         if any(math.isinf(norm(pos)) for pos in positions):
             raise ValueError('Particle position is infinite')
+        nonHeterogens = PDBFile._standardResidues[:]
+        nonHeterogens.remove('HOH')
         atomIndex = 1
         posIndex = 0
         if modelIndex is not None:
@@ -350,6 +351,10 @@ class PDBFile(object):
                     resId = res.id
                 else:
                     resId = "%4d" % ((resIndex+1)%10000)
+                if res.name in nonHeterogens:
+                    recordName = "ATOM  "
+                else:
+                    recordName = "HETATM"
                 for atom in res.atoms():
                     if atom.element is not None:
                         symbol = atom.element.symbol
@@ -362,8 +367,8 @@ class PDBFile(object):
                     else:
                         atomName = atom.name
                     coords = positions[posIndex]
-                    line = "ATOM  %5d %-4s %3s %s%4s    %s%s%s  1.00  0.00          %2s  " % (
-                        atomIndex%100000, atomName, resName, chainName, resId, _format_83(coords[0]),
+                    line = "%s%5d %-4s %3s %s%4s    %s%s%s  1.00  0.00          %2s  " % (
+                        recordName, atomIndex%100000, atomName, resName, chainName, resId, _format_83(coords[0]),
                         _format_83(coords[1]), _format_83(coords[2]), symbol)
                     assert len(line) == 80, 'Fixed width overflow detected'
                     print(line, file=file)
@@ -388,12 +393,9 @@ class PDBFile(object):
         """
         # Identify bonds that should be listed as CONECT records.
 
-        standardResidues = ['ALA', 'ASN', 'CYS', 'GLU', 'HIS', 'LEU', 'MET', 'PRO', 'THR', 'TYR',
-                            'ARG', 'ASP', 'GLN', 'GLY', 'ILE', 'LYS', 'PHE', 'SER', 'TRP', 'VAL',
-                            'A', 'G', 'C', 'U', 'I', 'DA', 'DG', 'DC', 'DT', 'DI', 'HOH']
         conectBonds = []
         for atom1, atom2 in topology.bonds():
-            if atom1.residue.name not in standardResidues or atom2.residue.name not in standardResidues:
+            if atom1.residue.name not in PDBFile._standardResidues or atom2.residue.name not in PDBFile._standardResidues:
                 conectBonds.append((atom1, atom2))
             elif atom1.name == 'SG' and atom2.name == 'SG' and atom1.residue.name == 'CYS' and atom2.residue.name == 'CYS':
                 conectBonds.append((atom1, atom2))

wrappers/python/tests/TestForceField.py

@@ -448,6 +448,183 @@ class TestForceField(unittest.TestCase):
         self.assertEqual(templates[1].name, 'ALA')
         self.assertEqual(templates[2].name, 'CALA')
 
+    def test_Wildcard(self):
+        """Test that PeriodicTorsionForces using wildcard ('') for atom types / classes in the ffxml are correctly registered"""
+
+        # Use wildcards in types
+        xml = """
+<ForceField>
+ <AtomTypes>
+  <Type name="C" class="C" element="C" mass="12.010000"/>
+  <Type name="O" class="O" element="O" mass="16.000000"/>
+ </AtomTypes>
+ <PeriodicTorsionForce>
+  <Proper type1="" type2="C" type3="C" type4="" periodicity1="2" phase1="3.141593" k1="15.167000"/>
+  <Improper type1="C" type2="" type3="" type4="O" periodicity1="2" phase1="3.141593" k1="43.932000"/>
+ </PeriodicTorsionForce>
+</ForceField>"""
+
+        ff = ForceField(StringIO(xml))
+
+        self.assertEqual(len(ff._forces[0].proper), 1)
+        self.assertEqual(len(ff._forces[0].improper), 1)
+
+       # Use wildcards in classes
+        xml = """
+<ForceField>
+ <AtomTypes>
+  <Type name="C" class="C" element="C" mass="12.010000"/>
+  <Type name="O" class="O" element="O" mass="16.000000"/>
+ </AtomTypes>
+ <PeriodicTorsionForce>
+  <Proper class1="" class2="C" class3="C" class4="" periodicity1="2" phase1="3.141593" k1="15.167000"/>
+  <Improper class1="C" class2="" class3="" class4="O" periodicity1="2" phase1="3.141593" k1="43.932000"/>
+ </PeriodicTorsionForce>
+</ForceField>"""
+
+        ff = ForceField(StringIO(xml))
+
+        self.assertEqual(len(ff._forces[0].proper), 1)
+        self.assertEqual(len(ff._forces[0].improper), 1)
+
+    def test_ScalingFactorCombining(self):
+        """ Tests that FFs can be combined if their scaling factors are very close """
+        forcefield = ForceField('amber99sb.xml', os.path.join('systems', 'test_amber_ff.xml'))
+        # This would raise an exception if it didn't work
+
+    def test_MultipleFilesandForceTags(self):
+        """Test that the order of listing of multiple ffxmls does not matter.
+           Tests that one generator per force type is created and that the ffxml
+           defining atom types does not have to be listed first"""
+
+        ffxml = """<ForceField>
+ <Residues>
+  <Residue name="ACE-Test">
+   <Atom name="HH31" type="710"/>
+   <Atom name="CH3" type="711"/>
+   <Atom name="HH32" type="710"/>
+   <Atom name="HH33" type="710"/>
+   <Atom name="C" type="712"/>
+   <Atom name="O" type="713"/>
+   <Bond from="0" to="1"/>
+   <Bond from="1" to="2"/>
+   <Bond from="1" to="3"/>
+   <Bond from="1" to="4"/>
+   <Bond from="4" to="5"/>
+   <ExternalBond from="4"/>
+  </Residue>
+ </Residues>
+ <PeriodicTorsionForce>
+  <Proper class1="C" class2="C" class3="C" class4="C" periodicity1="2" phase1="3.14159265359" k1="10.46"/>
+  <Improper class1="C" class2="C" class3="C" class4="C" periodicity1="2" phase1="3.14159265359" k1="43.932"/>
+ </PeriodicTorsionForce>
+</ForceField>"""
+
+        ff1 = ForceField(StringIO(ffxml), 'amber99sbildn.xml')
+        ff2 = ForceField('amber99sbildn.xml', StringIO(ffxml))
+
+        self.assertEqual(len(ff1._forces), 4)
+        self.assertEqual(len(ff2._forces), 4)
+
+        pertorsion1 = ff1._forces[0]
+        pertorsion2 = ff2._forces[2]
+
+        self.assertEqual(len(pertorsion1.proper), 110)
+        self.assertEqual(len(pertorsion1.improper), 42)
+        self.assertEqual(len(pertorsion2.proper), 110)
+        self.assertEqual(len(pertorsion2.improper), 42)
+
+    def test_ResidueTemplateUserChoice(self):
+        """Test createSystem does not allow multiple matching templates, unless
+           user has specified which template to use via residueTemplates arg"""
+        ffxml = """<ForceField>
+ <AtomTypes>
+  <Type name="Fe2+" class="Fe2+" element="Fe" mass="55.85"/>
+  <Type name="Fe3+" class="Fe3+" element="Fe" mass="55.85"/>
+ </AtomTypes>
+ <Residues>
+  <Residue name="FE2">
+   <Atom name="FE2" type="Fe2+" charge="2.0"/>
+  </Residue>
+  <Residue name="FE">
+   <Atom name="FE" type="Fe3+" charge="3.0"/>
+  </Residue>
+ </Residues>
+ <NonbondedForce coulomb14scale="0.833333333333" lj14scale="0.5">
+  <UseAttributeFromResidue name="charge"/>
+  <Atom type="Fe2+" sigma="0.227535532613" epsilon="0.0150312292"/>
+  <Atom type="Fe3+" sigma="0.192790482606" epsilon="0.00046095128"/>
+ </NonbondedForce>
+</ForceField>"""
+
+        pdb_string = "ATOM      1 FE    FE A   1      20.956  27.448 -29.067  1.00  0.00          Fe"
+        ff = ForceField(StringIO(ffxml))
+        pdb = PDBFile(StringIO(pdb_string))
+
+        self.assertRaises(Exception, lambda: ff.createSystem(pdb.topology))
+        sys = ff.createSystem(pdb.topology, residueTemplates={list(pdb.topology.residues())[0] : 'FE2'})
+        # confirm charge
+        self.assertEqual(sys.getForce(0).getParticleParameters(0)[0]._value, 2.0)
+        sys = ff.createSystem(pdb.topology, residueTemplates={list(pdb.topology.residues())[0] : 'FE'})
+        # confirm charge
+        self.assertEqual(sys.getForce(0).getParticleParameters(0)[0]._value, 3.0)
+
+    def test_ResidueOverloading(self):
+        """Test residue overloading via overload tag in the XML"""
+
+        ffxml1 = """<ForceField>
+ <AtomTypes>
+  <Type name="Fe2+_tip3p_HFE" class="Fe2+_tip3p_HFE" element="Fe" mass="55.85"/>
+ </AtomTypes>
+ <Residues>
+  <Residue name="FE2">
+   <Atom name="FE2" type="Fe2+_tip3p_HFE" charge="2.0"/>
+  </Residue>
+ </Residues>
+ <NonbondedForce coulomb14scale="0.833333333333" lj14scale="0.5">
+  <UseAttributeFromResidue name="charge"/>
+  <Atom type="Fe2+_tip3p_HFE" sigma="0.227535532613" epsilon="0.0150312292"/>
+ </NonbondedForce>
+</ForceField>"""
+
+        ffxml2 = """<ForceField>
+ <AtomTypes>
+  <Type name="Fe2+_tip3p_standard" class="Fe2+_tip3p_standard" element="Fe" mass="55.85"/>
+ </AtomTypes>
+ <Residues>
+  <Residue name="FE2">
+   <Atom name="FE2" type="Fe2+_tip3p_standard" charge="2.0"/>
+  </Residue>
+ </Residues>
+ <NonbondedForce coulomb14scale="0.833333333333" lj14scale="0.5">
+  <UseAttributeFromResidue name="charge"/>
+  <Atom type="Fe2+_tip3p_standard" sigma="0.241077193129" epsilon="0.03940482832"/>
+ </NonbondedForce>
+</ForceField>"""
+
+        ffxml3 = """<ForceField>
+ <AtomTypes>
+  <Type name="Fe2+_tip3p_standard" class="Fe2+_tip3p_standard" element="Fe" mass="55.85"/>
+ </AtomTypes>
+ <Residues>
+  <Residue name="FE2" overload="1">
+   <Atom name="FE2" type="Fe2+_tip3p_standard" charge="2.0"/>
+  </Residue>
+ </Residues>
+ <NonbondedForce coulomb14scale="0.833333333333" lj14scale="0.5">
+  <UseAttributeFromResidue name="charge"/>
+  <Atom type="Fe2+_tip3p_standard" sigma="0.241077193129" epsilon="0.03940482832"/>
+ </NonbondedForce>
+</ForceField>"""
+
+        pdb_string = "ATOM      1 FE    FE A   1      20.956  27.448 -29.067  1.00  0.00          Fe"
+        pdb = PDBFile(StringIO(pdb_string))
+
+        self.assertRaises(Exception, lambda: ForceField(StringIO(ffxml1), StringIO(ffxml2)))
+        ff = ForceField(StringIO(ffxml1), StringIO(ffxml3))
+        self.assertEqual(ff._templates['FE2'].atoms[0].type, 'Fe2+_tip3p_standard')
+        ff.createSystem(pdb.topology)
+
 class AmoebaTestForceField(unittest.TestCase):
     """Test the ForceField.createSystem() method with the AMOEBA forcefield."""
 
@@ -535,49 +712,5 @@ class AmoebaTestForceField(unittest.TestCase):
             diff = norm(f1-f2)
             self.assertTrue(diff < 0.1 or diff/norm(f1) < 1e-3)
 
-    def test_Wildcard(self):
-        """Test that PeriodicTorsionForces using wildcard ('') for atom types / classes in the ffxml are correctly registered"""
-
-        # Use wildcards in types
-        xml = """
-<ForceField>
- <AtomTypes>
-  <Type name="C" class="C" element="C" mass="12.010000"/>
-  <Type name="O" class="O" element="O" mass="16.000000"/>
- </AtomTypes>
- <PeriodicTorsionForce>
-  <Proper type1="" type2="C" type3="C" type4="" periodicity1="2" phase1="3.141593" k1="15.167000"/>
-  <Improper type1="C" type2="" type3="" type4="O" periodicity1="2" phase1="3.141593" k1="43.932000"/>
- </PeriodicTorsionForce>
-</ForceField>"""
-
-        ff = ForceField(StringIO(xml))
-
-        self.assertEqual(len(ff._forces[0].proper), 1)
-        self.assertEqual(len(ff._forces[0].improper), 1)
-
-       # Use wildcards in classes
-        xml = """
-<ForceField>
- <AtomTypes>
-  <Type name="C" class="C" element="C" mass="12.010000"/>
-  <Type name="O" class="O" element="O" mass="16.000000"/>
- </AtomTypes>
- <PeriodicTorsionForce>
-  <Proper class1="" class2="C" class3="C" class4="" periodicity1="2" phase1="3.141593" k1="15.167000"/>
-  <Improper class1="C" class2="" class3="" class4="O" periodicity1="2" phase1="3.141593" k1="43.932000"/>
- </PeriodicTorsionForce>
-</ForceField>"""
-
-        ff = ForceField(StringIO(xml))
-
-        self.assertEqual(len(ff._forces[0].proper), 1)
-        self.assertEqual(len(ff._forces[0].improper), 1)
-
-    def test_ScalingFactorCombining(self):
-        """ Tests that FFs can be combined if their scaling factors are very close """
-        forcefield = ForceField('amber99sb.xml', os.path.join('systems', 'test_amber_ff.xml'))
-        # This would raise an exception if it didn't work
-
 if __name__ == '__main__':
     unittest.main()