Merge branch 'master' into update_ug

plugins/cpupme/src/CpuPmeKernels.cpp

@@ -37,6 +37,7 @@
 #include "openmm/internal/hardware.h"
 #include "openmm/internal/vectorize.h"
 #include <cmath>
+#include <algorithm>
 #include <cstring>
 #include <sstream>
 #include <cstdlib>
@@ -395,7 +396,7 @@ void CpuCalcPmeReciprocalForceKernel::initialize(int xsize, int ysize, int zsize
     
     // Initialize the b-spline moduli.
 
-    int maxSize = max(max(gridx, gridy), gridz);
+    int maxSize = std::max(std::max(gridx, gridy), gridz);
     vector<double> data(PME_ORDER);
     vector<double> ddata(PME_ORDER);
     vector<double> bsplinesData(maxSize);
@@ -521,7 +522,7 @@ void CpuCalcPmeReciprocalForceKernel::runWorkerThread(ThreadPool& threads, int i
     int gridStart = 4*((index*gridSize)/numThreads);
     int gridEnd = 4*(((index+1)*gridSize)/numThreads);
     int complexSize = gridx*gridy*(gridz/2+1);
-    int complexStart = max(1, ((index*complexSize)/numThreads));
+    int complexStart = std::max(1, ((index*complexSize)/numThreads));
     int complexEnd = (((index+1)*complexSize)/numThreads);
     spreadCharge(posq, tempGrid[index], gridx, gridy, gridz, numParticles, periodicBoxVectors, recipBoxVectors, atomicCounter);
     threads.syncThreads();

plugins/rpmd/openmmapi/src/RPMDIntegrator.cpp

@@ -129,9 +129,9 @@ State RPMDIntegrator::getState(int copy, int types, bool enforcePeriodicBox, int
 
             // Find the displacement to move it into the first periodic box.
             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]);
+            diff += periodicBoxSize[2]*floor(center[2]/periodicBoxSize[2][2]);
+            diff += periodicBoxSize[1]*floor((center[1]-diff[1])/periodicBoxSize[1][1]);
+            diff += periodicBoxSize[0]*floor((center[0]-diff[0])/periodicBoxSize[0][0]);
 
             // Translate all the particles in the molecule.
             for (int j = 0; j < (int) molecules[i].size(); j++) {

serialization/include/openmm/serialization/CustomCentroidBondForceProxy.h

+#ifndef OPENMM_CUSTOMCENTROIDBONDFORCE_PROXY_H_
+#define OPENMM_CUSTOMCENTROIDBONDFORCE_PROXY_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   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) 2015 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/internal/windowsExport.h"
+#include "openmm/serialization/SerializationProxy.h"
+
+namespace OpenMM {
+
+/**
+ * This is a proxy for serializing CustomCentroidBondForce objects.
+ */
+
+class OPENMM_EXPORT CustomCentroidBondForceProxy : public SerializationProxy {
+public:
+    CustomCentroidBondForceProxy();
+    void serialize(const void* object, SerializationNode& node) const;
+    void* deserialize(const SerializationNode& node) const;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CUSTOMCENTROIDBONDFORCE_PROXY_H_*/

serialization/src/CustomCentroidBondForceProxy.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: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/serialization/CustomCentroidBondForceProxy.h"
+#include "openmm/serialization/SerializationNode.h"
+#include "openmm/Force.h"
+#include "openmm/CustomCentroidBondForce.h"
+#include <sstream>
+
+using namespace OpenMM;
+using namespace std;
+
+CustomCentroidBondForceProxy::CustomCentroidBondForceProxy() : SerializationProxy("CustomCentroidBondForce") {
+}
+
+void CustomCentroidBondForceProxy::serialize(const void* object, SerializationNode& node) const {
+    node.setIntProperty("version", 1);
+    const CustomCentroidBondForce& force = *reinterpret_cast<const CustomCentroidBondForce*>(object);
+    node.setIntProperty("forceGroup", force.getForceGroup());
+    node.setIntProperty("groups", force.getNumGroupsPerBond());
+    node.setStringProperty("energy", force.getEnergyFunction());
+    SerializationNode& perBondParams = node.createChildNode("PerBondParameters");
+    for (int i = 0; i < force.getNumPerBondParameters(); i++) {
+        perBondParams.createChildNode("Parameter").setStringProperty("name", force.getPerBondParameterName(i));
+    }
+    SerializationNode& globalParams = node.createChildNode("GlobalParameters");
+    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
+        globalParams.createChildNode("Parameter").setStringProperty("name", force.getGlobalParameterName(i)).setDoubleProperty("default", force.getGlobalParameterDefaultValue(i));
+    }
+    SerializationNode& groups = node.createChildNode("Groups");
+    for (int i = 0; i < force.getNumGroups(); i++) {
+        vector<int> particles;
+        vector<double> weights;
+        force.getGroupParameters(i, particles, weights);
+        SerializationNode& group = groups.createChildNode("Group");
+        for (int j = 0; j < (int) particles.size(); j++) {
+            SerializationNode& node = group.createChildNode("Particle");
+            node.setIntProperty("p", particles[j]);
+            if (j < weights.size())
+                node.setDoubleProperty("weight", weights[j]);
+        }
+    }
+    SerializationNode& bonds = node.createChildNode("Bonds");
+    for (int i = 0; i < force.getNumBonds(); i++) {
+        vector<int> groups;
+        vector<double> params;
+        force.getBondParameters(i, groups, params);
+        SerializationNode& node = bonds.createChildNode("Bond");
+        for (int j = 0; j < (int) groups.size(); j++) {
+            stringstream key;
+            key << "g";
+            key << j+1;
+            node.setIntProperty(key.str(), groups[j]);
+        }
+        for (int j = 0; j < (int) params.size(); j++) {
+            stringstream key;
+            key << "param";
+            key << j+1;
+            node.setDoubleProperty(key.str(), params[j]);
+        }
+    }
+    SerializationNode& functions = node.createChildNode("Functions");
+    for (int i = 0; i < force.getNumTabulatedFunctions(); i++)
+        functions.createChildNode("Function", &force.getTabulatedFunction(i)).setStringProperty("name", force.getTabulatedFunctionName(i));
+}
+
+void* CustomCentroidBondForceProxy::deserialize(const SerializationNode& node) const {
+    if (node.getIntProperty("version") != 1)
+        throw OpenMMException("Unsupported version number");
+    CustomCentroidBondForce* force = NULL;
+    try {
+        CustomCentroidBondForce* force = new CustomCentroidBondForce(node.getIntProperty("groups"), node.getStringProperty("energy"));
+        force->setForceGroup(node.getIntProperty("forceGroup", 0));
+        const SerializationNode& perBondParams = node.getChildNode("PerBondParameters");
+        for (int i = 0; i < (int) perBondParams.getChildren().size(); i++) {
+            const SerializationNode& parameter = perBondParams.getChildren()[i];
+            force->addPerBondParameter(parameter.getStringProperty("name"));
+        }
+        const SerializationNode& globalParams = node.getChildNode("GlobalParameters");
+        for (int i = 0; i < (int) globalParams.getChildren().size(); i++) {
+            const SerializationNode& parameter = globalParams.getChildren()[i];
+            force->addGlobalParameter(parameter.getStringProperty("name"), parameter.getDoubleProperty("default"));
+        }
+        const SerializationNode& groups = node.getChildNode("Groups");
+        for (int i = 0; i < (int) groups.getChildren().size(); i++) {
+            const SerializationNode& group = groups.getChildren()[i];
+            vector<int> particles;
+            vector<double> weights;
+            for (int j = 0; j < (int) group.getChildren().size(); j++) {
+                particles.push_back(group.getChildren()[j].getIntProperty("p"));
+                if (group.getChildren()[j].hasProperty("weight"))
+                    weights.push_back(group.getChildren()[j].getDoubleProperty("weight"));
+            }
+            force->addGroup(particles, weights);
+        }
+        const SerializationNode& bonds = node.getChildNode("Bonds");
+        vector<int> bondGroups(force->getNumGroupsPerBond());
+        vector<double> params(force->getNumPerBondParameters());
+        for (int i = 0; i < (int) bonds.getChildren().size(); i++) {
+            const SerializationNode& bond = bonds.getChildren()[i];
+            for (int j = 0; j < (int) bondGroups.size(); j++) {
+                stringstream key;
+                key << "g";
+                key << j+1;
+                bondGroups[j] = bond.getIntProperty(key.str());
+            }
+            for (int j = 0; j < (int) params.size(); j++) {
+                stringstream key;
+                key << "param";
+                key << j+1;
+                params[j] = bond.getDoubleProperty(key.str());
+            }
+            force->addBond(bondGroups, params);
+        }
+        const SerializationNode& functions = node.getChildNode("Functions");
+        for (int i = 0; i < (int) functions.getChildren().size(); i++) {
+            const SerializationNode& function = functions.getChildren()[i];
+            if (function.hasProperty("type")) {
+                force->addTabulatedFunction(function.getStringProperty("name"), function.decodeObject<TabulatedFunction>());
+            }
+            else {
+                // This is an old file created before TabulatedFunction existed.
+
+                const SerializationNode& valuesNode = function.getChildNode("Values");
+                vector<double> values;
+                for (int j = 0; j < (int) valuesNode.getChildren().size(); j++)
+                    values.push_back(valuesNode.getChildren()[j].getDoubleProperty("v"));
+                force->addTabulatedFunction(function.getStringProperty("name"), new Continuous1DFunction(values, function.getDoubleProperty("min"), function.getDoubleProperty("max")));
+            }
+        }
+        return force;
+    }
+    catch (...) {
+        if (force != NULL)
+            delete force;
+        throw;
+    }
+}

serialization/src/SerializationProxyRegistration.cpp

@@ -36,6 +36,7 @@
 #include "openmm/CompoundIntegrator.h"
 #include "openmm/CustomAngleForce.h"
 #include "openmm/CustomBondForce.h"
+#include "openmm/CustomCentroidBondForce.h"
 #include "openmm/CustomCompoundBondForce.h"
 #include "openmm/CustomExternalForce.h"
 #include "openmm/CustomGBForce.h"
@@ -68,6 +69,7 @@
 #include "openmm/serialization/CompoundIntegratorProxy.h"
 #include "openmm/serialization/CustomAngleForceProxy.h"
 #include "openmm/serialization/CustomBondForceProxy.h"
+#include "openmm/serialization/CustomCentroidBondForceProxy.h"
 #include "openmm/serialization/CustomCompoundBondForceProxy.h"
 #include "openmm/serialization/CustomExternalForceProxy.h"
 #include "openmm/serialization/CustomGBForceProxy.h"
@@ -118,6 +120,7 @@ extern "C" void registerSerializationProxies() {
     SerializationProxy::registerProxy(typeid(Continuous3DFunction), new Continuous3DFunctionProxy());
     SerializationProxy::registerProxy(typeid(CustomAngleForce), new CustomAngleForceProxy());
     SerializationProxy::registerProxy(typeid(CustomBondForce), new CustomBondForceProxy());
+    SerializationProxy::registerProxy(typeid(CustomCentroidBondForce), new CustomCentroidBondForceProxy());
     SerializationProxy::registerProxy(typeid(CustomCompoundBondForce), new CustomCompoundBondForceProxy());
     SerializationProxy::registerProxy(typeid(CustomExternalForce), new CustomExternalForceProxy());
     SerializationProxy::registerProxy(typeid(CustomGBForce), new CustomGBForceProxy());

serialization/tests/TestSerializeCustomCentroidBondForce.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: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/CustomCentroidBondForce.h"
+#include "openmm/serialization/XmlSerializer.h"
+#include <iostream>
+#include <sstream>
+
+using namespace OpenMM;
+using namespace std;
+
+void testSerialization() {
+    // Create a Force.
+
+    CustomCentroidBondForce force(3, "5*sin(distance(g1,g2))^2+y*z");
+    force.setForceGroup(3);
+    force.addGlobalParameter("x", 1.3);
+    force.addGlobalParameter("y", 2.221);
+    force.addPerBondParameter("z");
+    for (int i = 0; i < 3; i++) {
+        vector<int> particles;
+        vector<double> weights;
+        for (int j = 0; j < i+1; j++) {
+            particles.push_back(i+j);
+            if (i < 2)
+                weights.push_back(1.0/(i+1));
+        }
+        force.addGroup(particles, weights);
+    }
+    vector<int> groups(3);
+    vector<double> params(1);
+    groups[0] = 0;
+    groups[1] = 1;
+    groups[2] = 2;
+    params[0] = 1.0;
+    force.addBond(groups, params);
+    groups[0] = 2;
+    groups[1] = 3;
+    groups[2] = 4;
+    params[0] = -3.3;
+    force.addBond(groups, params);
+    groups[0] = 3;
+    groups[1] = 5;
+    groups[2] = 1;
+    params[0] = 2.1;
+    force.addBond(groups, params);
+    vector<double> values(10);
+    for (int i = 0; i < 10; i++)
+        values[i] = sin((double) i);
+    force.addTabulatedFunction("f", new Continuous1DFunction(values, 0.5, 1.5));
+
+    // Serialize and then deserialize it.
+
+    stringstream buffer;
+    XmlSerializer::serialize<CustomCentroidBondForce>(&force, "Force", buffer);
+    CustomCentroidBondForce* copy = XmlSerializer::deserialize<CustomCentroidBondForce>(buffer);
+
+    // Compare the two forces to see if they are identical.
+
+    CustomCentroidBondForce& force2 = *copy;
+    ASSERT_EQUAL(force.getForceGroup(), force2.getForceGroup());
+    ASSERT_EQUAL(force.getNumGroupsPerBond(), force2.getNumGroupsPerBond());
+    ASSERT_EQUAL(force.getEnergyFunction(), force2.getEnergyFunction());
+    ASSERT_EQUAL(force.getNumPerBondParameters(), force2.getNumPerBondParameters());
+    for (int i = 0; i < force.getNumPerBondParameters(); i++)
+        ASSERT_EQUAL(force.getPerBondParameterName(i), force2.getPerBondParameterName(i));
+    ASSERT_EQUAL(force.getNumGlobalParameters(), force2.getNumGlobalParameters());
+    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
+        ASSERT_EQUAL(force.getGlobalParameterName(i), force2.getGlobalParameterName(i));
+        ASSERT_EQUAL(force.getGlobalParameterDefaultValue(i), force2.getGlobalParameterDefaultValue(i));
+    }
+    ASSERT_EQUAL(force.getNumGroups(), force2.getNumGroups());
+    for (int i = 0; i < force.getNumGroups(); i++) {
+        vector<int> particles1, particles2;
+        vector<double> weights1, weights2;
+        force.getGroupParameters(i, particles1, weights1);
+        force2.getGroupParameters(i, particles2, weights2);
+        ASSERT_EQUAL(weights1.size(), weights2.size());
+        for (int j = 0; j < (int) weights1.size(); j++)
+            ASSERT_EQUAL(weights1[j], weights2[j]);
+        ASSERT_EQUAL(particles1.size(), particles2.size());
+        for (int j = 0; j < (int) particles1.size(); j++)
+            ASSERT_EQUAL(particles1[j], particles2[j]);
+    }
+    ASSERT_EQUAL(force.getNumBonds(), force2.getNumBonds());
+    for (int i = 0; i < force.getNumBonds(); i++) {
+        vector<int> groups1, groups2;
+        vector<double> params1, params2;
+        force.getBondParameters(i, groups1, params1);
+        force2.getBondParameters(i, groups2, params2);
+        ASSERT_EQUAL(params1.size(), params2.size());
+        for (int j = 0; j < (int) params1.size(); j++)
+            ASSERT_EQUAL(params1[j], params2[j]);
+        ASSERT_EQUAL(groups1.size(), groups2.size());
+        for (int j = 0; j < (int) groups1.size(); j++)
+            ASSERT_EQUAL(groups1[j], groups2[j]);
+    }
+    ASSERT_EQUAL(force.getNumTabulatedFunctions(), force2.getNumTabulatedFunctions());
+    for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
+        double min1, min2, max1, max2;
+        vector<double> val1, val2;
+        dynamic_cast<Continuous1DFunction&>(force.getTabulatedFunction(i)).getFunctionParameters(val1, min1, max1);
+        dynamic_cast<Continuous1DFunction&>(force2.getTabulatedFunction(i)).getFunctionParameters(val2, min2, max2);
+        ASSERT_EQUAL(force.getTabulatedFunctionName(i), force2.getTabulatedFunctionName(i));
+        ASSERT_EQUAL(min1, min2);
+        ASSERT_EQUAL(max1, max2);
+        ASSERT_EQUAL(val1.size(), val2.size());
+        for (int j = 0; j < (int) val1.size(); j++)
+            ASSERT_EQUAL(val1[j], val2[j]);
+    }
+}
+
+int main() {
+    try {
+        testSerialization();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

tests/TestEnforcePeriodicBox.cpp

@@ -41,7 +41,7 @@ using namespace OpenMM;
 using namespace std;
 
 void testTruncatedOctahedron() {
-    const int numMolecules = 5;
+    const int numMolecules = 50;
     const int numParticles = numMolecules*2;
     const float cutoff = 2.0;
     Vec3 a(6.7929, 0, 0);
@@ -63,7 +63,7 @@ void testTruncatedOctahedron() {
         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] = a*(5*genrand_real2(sfmt)-2) + b*(5*genrand_real2(sfmt)-2) + c*(5*genrand_real2(sfmt)-2);
         positions[2*i+1] = positions[2*i] + Vec3(1.0, 0.0, 0.0);
         system.addConstraint(2*i, 2*i+1, 1.0);
     }
@@ -73,6 +73,15 @@ void testTruncatedOctahedron() {
     Context context(system, integrator, Platform::getPlatformByName("Reference"));
     context.setPositions(positions);
     State initialState = context.getState(State::Positions | State::Energy, true);
+    for (int i = 0; i < numMolecules; i++) {
+        Vec3 center = (initialState.getPositions()[2*i]+initialState.getPositions()[2*i+1])*0.5;
+        ASSERT(center[0] >= 0.0);
+        ASSERT(center[1] >= 0.0);
+        ASSERT(center[2] >= 0.0);
+        ASSERT(center[0] <= a[0]);
+        ASSERT(center[1] <= b[1]);
+        ASSERT(center[2] <= c[2]);
+    }
     double initialEnergy = initialState.getPotentialEnergy();
 
     context.setState(initialState);

wrappers/python/CMakeLists.txt

@@ -20,7 +20,6 @@ set(STAGING_OUTPUT_FILES "") # Will contain all required package files
 file(GLOB STAGING_INPUT_FILES RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}"
     "${CMAKE_CURRENT_SOURCE_DIR}/MANIFEST.in"
     "${CMAKE_CURRENT_SOURCE_DIR}/README.txt"
-    "${CMAKE_CURRENT_SOURCE_DIR}/filterPythonFiles.py"
 )
 configure_file(${CMAKE_CURRENT_SOURCE_DIR}/setup.py ${OPENMM_PYTHON_STAGING_DIR}/setup.py)
 
@@ -51,6 +50,7 @@ foreach(SUBDIR ${SUBDIRS})
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*README.txt"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.py"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.i"
+        "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.ini"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.sh"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.xml"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.pdb"
@@ -185,7 +185,7 @@ add_custom_command(
     COMMENT "Creating OpenMM Python swig input files..."
 )
 
-INSTALL_FILES(/include FILES "${SWIG_OPENMM_DIR}/OpenMMSwigHeaders.i")
+INSTALL_FILES(/include/swig FILES "${SWIG_OPENMM_DIR}/OpenMMSwigHeaders.i" "${SWIG_OPENMM_DIR}/swig_lib/python/typemaps.i")
 
 #~ swig -python -c++ \
 #~      -outdir $PYTHON_PACKAGE_DIR \

wrappers/python/filterPythonFiles.py

-
-from __future__ import print_function
-import re
-import sys
-import textwrap
-from inspect import cleandoc
-from numpydoc.docscrape import NumpyDocString
-
-
-# Doxygen does a bad job of generating documentation based on docstrings.
-# This script is run as a filter
-# on each file, and converts the docstrings into Doxygen style comments
-# so we get better documentation.
-
-input = open(sys.argv[1])
-
-while True:
-    line = input.readline()
-    if len(line) == 0:
-        break
-    stripped = line.lstrip()
-    if stripped.startswith('def') or stripped.startswith('class'):
-        prefix = line[:len(line)-len(stripped)]
-        split = stripped.split()
-        if split[0] == 'class' and split[1][0].islower():
-            # Classes that start with a lowercase letter were defined by SWIG.  We want to hide them.
-            print("%s## @private" % prefix, end='')
-        if split[1][0] == '_' and split[1][1] != '_':
-            # Names starting with a single _ are assumed to be private.
-            print("%s## @private" % prefix, end='')
-
-        # We're at the start of a class or function definition.  Find all lines that contain the declaration.
-
-        declaration = line
-        while len(line) > 0 and line.find(':') == -1:
-            line = input.readline()
-            declaration += line
-
-        # Now look for a docstring.
-        docstringlines = []
-        line = input.readline()
-        l = line.strip()
-        if l.startswith('"""') or l.startswith("'''"):
-            if l.count('"""') == 2 or l.count("'''") == 2:
-                docstringlines.append(l[3:-3])
-            else:
-                docstringlines.append(l[3:])
-                line = input.readline()
-                l = line.strip()
-                while l.find('"""') == -1 and l.find("'''") == -1:
-                    docstringlines.append(line.rstrip())
-                    line = input.readline()
-                    l = line.strip()
-                if line.rstrip().endswith('"""') or line.rstrip().endswith("'''"):
-                    # last line
-                    docstringlines.append(line.rstrip()[:-3])
-
-        docstring = NumpyDocString(cleandoc('\n'.join(docstringlines)))
-        # print(docstringlines)
-
-        # Print out the docstring in Doxygen syntax, followed by the declaration.
-        for line in docstring['Summary']:
-            print('{prefix}## {line}'.format(prefix=prefix, line=line))
-        if len(docstring['Extended Summary']) > 0:
-            print('{prefix}##'.format(prefix=prefix))
-            for line in docstring['Extended Summary']:
-                print('{prefix}## {line}'.format(prefix=prefix, line=line.strip()))
-        print('{prefix}##'.format(prefix=prefix))
-        for name, type, descr in docstring['Parameters']:
-            print('{prefix}## @param {name} ({type}) {descr}'.format(prefix=prefix, type=type, name=name, descr=' '.join(descr)))
-        for name, type, descr in docstring['Returns']:
-            if type == '':
-                type = name
-            print('{prefix}## @return ({type}) {descr}'.format(prefix=prefix, type=type, name=name, descr=' '.join(descr)))
-
-        print(declaration, end='')
-        if len(docstringlines) == 0:
-            print(line, end='')
-    else:
-        print(line, end='')

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

@@ -38,7 +38,7 @@ from simtk.openmm.app import PDBFile
 from simtk.openmm.app.internal import amber_file_parser
 from . import forcefield as ff
 from . import element as elem
-import simtk.unit as unit
+import simtk.unit as u
 import simtk.openmm as mm
 from simtk.openmm.app.internal.unitcell import computePeriodicBoxVectors
 
@@ -69,6 +69,15 @@ class GBn2(object):
         return 'GBn2'
 GBn2 = GBn2()
 
+def _strip_optunit(thing, unit):
+    """
+    Strips optional units, converting to specified unit type. If no unit
+    present, it just returns the number
+    """
+    if u.is_quantity(thing):
+        return thing.value_in_unit(unit)
+    return thing
+
 class AmberPrmtopFile(object):
     """AmberPrmtopFile parses an AMBER prmtop file and constructs a Topology and (optionally) an OpenMM System from it."""
 
@@ -145,12 +154,13 @@ class AmberPrmtopFile(object):
             box = prmtop.getBoxBetaAndDimensions()
             top.setPeriodicBoxVectors(computePeriodicBoxVectors(*(box[1:4] + box[0:1]*3)))
 
-    def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
+    def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*u.nanometer,
                      constraints=None, rigidWater=True, implicitSolvent=None,
-                     implicitSolventSaltConc=0.0*(unit.moles/unit.liter),
-                     implicitSolventKappa=None, temperature=298.15*unit.kelvin,
+                     implicitSolventSaltConc=0.0*(u.moles/u.liter),
+                     implicitSolventKappa=None, temperature=298.15*u.kelvin,
                      soluteDielectric=1.0, solventDielectric=78.5,
-                     removeCMMotion=True, hydrogenMass=None, ewaldErrorTolerance=0.0005):
+                     removeCMMotion=True, hydrogenMass=None, ewaldErrorTolerance=0.0005,
+                     switchDistance=0.0*u.nanometer):
         """Construct an OpenMM System representing the topology described by this
         prmtop file.
 
@@ -192,6 +202,11 @@ class AmberPrmtopFile(object):
             total mass the same.
         ewaldErrorTolerance : float=0.0005
             The error tolerance to use if nonbondedMethod is Ewald or PME.
+        switchDistance : float=0*nanometers
+            The distance at which the potential energy switching function is
+            turned on for Lennard-Jones interactions. If the switchDistance is 0
+            or evaluates to boolean False, no switching function will be used.
+            Values greater than nonbondedCutoff or less than 0 raise ValueError
 
         Returns
         -------
@@ -237,10 +252,10 @@ class AmberPrmtopFile(object):
             raise ValueError('Illegal value for implicit solvent model')
         # If implicitSolventKappa is None, compute it from the salt concentration
         if implicitSolvent is not None and implicitSolventKappa is None:
-            if unit.is_quantity(implicitSolventSaltConc):
-                implicitSolventSaltConc = implicitSolventSaltConc.value_in_unit(unit.moles/unit.liter)
-            if unit.is_quantity(temperature):
-                temperature = temperature.value_in_unit(unit.kelvin)
+            if u.is_quantity(implicitSolventSaltConc):
+                implicitSolventSaltConc = implicitSolventSaltConc.value_in_unit(u.moles/u.liter)
+            if u.is_quantity(temperature):
+                temperature = temperature.value_in_unit(u.kelvin)
             # The constant is 1 / sqrt( epsilon_0 * kB / (2 * NA * q^2 * 1000) )
             # where NA is avogadro's number, epsilon_0 is the permittivity of
             # free space, q is the elementary charge (this number matches
@@ -271,4 +286,17 @@ class AmberPrmtopFile(object):
                 force.setEwaldErrorTolerance(ewaldErrorTolerance)
         if removeCMMotion:
             sys.addForce(mm.CMMotionRemover())
+
+        if switchDistance and nonbondedMethod is not ff.NoCutoff:
+            # make sure it's legal
+            if (_strip_optunit(switchDistance, u.nanometer) >=
+                    _strip_optunit(nonbondedCutoff, u.nanometer)):
+                raise ValueError('switchDistance is too large compared '
+                                 'to the cutoff!')
+            if _strip_optunit(switchDistance, u.nanometer) < 0:
+                # Detects negatives for both Quantity and float
+                raise ValueError('switchDistance must be non-negative!')
+            force.setUseSwitchingFunction(True)
+            force.setSwitchingDistance(switchDistance)
+
         return sys

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

@@ -32,9 +32,7 @@ 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.
 """
-from __future__ import division
-from __future__ import absolute_import
-from __future__ import print_function
+from __future__ import division, absolute_import, print_function
 
 from functools import wraps
 from math import pi, cos, sin, sqrt
@@ -47,7 +45,7 @@ import simtk.unit as u
 from simtk.openmm.app import (forcefield as ff, Topology, element)
 from simtk.openmm.app.amberprmtopfile import HCT, OBC1, OBC2, GBn, GBn2
 from simtk.openmm.app.internal.customgbforces import (GBSAHCTForce,
-                GBSAOBC1Force, GBSAOBC2Force, GBSAGBnForce, GBSAGBn2Force, convertParameters)
+                GBSAOBC1Force, GBSAOBC2Force, GBSAGBnForce, GBSAGBn2Force)
 from simtk.openmm.app.internal.unitcell import computePeriodicBoxVectors
 # CHARMM imports
 from simtk.openmm.app.internal.charmm.topologyobjects import (
@@ -666,88 +664,6 @@ class CharmmPsfFile(object):
 
         return topology
 
-    def _get_gb_params(self, gb_model=HCT):
-        """ Gets the GB parameters. Need this method to special-case GB neck """
-        screen = [0 for atom in self.atom_list]
-        if gb_model is GBn:
-            radii = _bondi_radii(self.atom_list)
-            screen = [0.5 for atom in self.atom_list]
-            for i, atom in enumerate(self.atom_list):
-                if atom.type.atomic_number == 6:
-                    screen[i] = 0.48435382330
-                elif atom.type.atomic_number == 1:
-                    screen[i] = 1.09085413633
-                elif atom.type.atomic_number == 7:
-                    screen[i] = 0.700147318409
-                elif atom.type.atomic_number == 8:
-                    screen[i] = 1.06557401132
-                elif atom.type.atomic_number == 16:
-                    screen[i] = 0.602256336067
-        elif gb_model is GBn2:
-            radii = _mbondi3_radii(self.atom_list)
-            # Add non-optimized values as defaults
-            alpha = [1.0 for i in self.atom_list]
-            beta = [0.8 for i in self.atom_list]
-            gamma = [4.85 for i in self.atom_list]
-            screen = [0.5 for i in self.atom_list]
-            for i, atom in enumerate(self.atom_list):
-                if atom.type.atomic_number == 6:
-                    screen[i] = 1.058554
-                    alpha[i] = 0.733756
-                    beta[i] = 0.506378
-                    gamma[i] = 0.205844
-                elif atom.type.atomic_number == 1:
-                    screen[i] = 1.425952
-                    alpha[i] = 0.788440
-                    beta[i] = 0.798699
-                    gamma[i] = 0.437334
-                elif atom.type.atomic_number == 7:
-                    screen[i] = 0.733599
-                    alpha[i] = 0.503364
-                    beta[i] = 0.316828
-                    gamma[i] = 0.192915
-                elif atom.type.atomic_number == 8:
-                    screen[i] = 1.061039
-                    alpha[i] = 0.867814
-                    beta[i] = 0.876635
-                    gamma[i] = 0.387882
-                elif atom.type.atomic_number == 16:
-                    screen[i] = -0.703469
-                    alpha[i] = 0.867814
-                    beta[i] = 0.876635
-                    gamma[i] = 0.387882
-        else:
-            # Set the default screening parameters
-            for i, atom in enumerate(self.atom_list):
-                if atom.type.atomic_number == 1:
-                    screen[i] = 0.85
-                elif atom.type.atomic_number == 6:
-                    screen[i] = 0.72
-                elif atom.type.atomic_number == 7:
-                    screen[i] = 0.79
-                elif atom.type.atomic_number == 8:
-                    screen[i] = 0.85
-                elif atom.type.atomic_number == 9:
-                    screen[i] = 0.88
-                elif atom.type.atomic_number == 15:
-                    screen[i] = 0.86
-                elif atom.type.atomic_number == 16:
-                    screen[i] = 0.96
-                else:
-                    screen[i] = 0.8
-            # Determine which radii set we need
-            if gb_model is OBC1 or gb_model is OBC2:
-                radii = _mbondi2_radii(self.atom_list)
-            elif gb_model is HCT:
-                radii = _mbondi_radii(self.atom_list)
-
-        length_conv = u.angstrom.conversion_factor_to(u.nanometer)
-        radii = [x * length_conv for x in radii]
-
-        if gb_model is GBn2:
-            return zip(radii, screen, alpha, beta, gamma)
-        return zip(radii, screen)
-
     def createSystem(self, params, nonbondedMethod=ff.NoCutoff,
                      nonbondedCutoff=1.0*u.nanometer,
                      switchDistance=0.0*u.nanometer,
@@ -1249,9 +1165,6 @@ class CharmmPsfFile(object):
         # Add GB model if we're doing one
         if implicitSolvent is not None:
             if verbose: print('Adding GB parameters...')
-            gb_parms = self._get_gb_params(implicitSolvent)
-            gb_parms = convertParameters(gb_parms, str(implicitSolvent))
-
             # If implicitSolventKappa is None, compute it from salt
             # concentration
             if implicitSolventKappa is None:
@@ -1290,6 +1203,7 @@ class CharmmPsfFile(object):
             elif implicitSolvent is GBn2:
                 gb = GBSAGBn2Force(solventDielectric, soluteDielectric, None,
                                    cutoff, kappa=implicitSolventKappa)
+            gb_parms = gb.getStandardParameters(self.topology)
             for atom, gb_parm in zip(self.atom_list, gb_parms):
                 gb.addParticle([atom.charge] + list(gb_parm))
             # Set cutoff method

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

 <ForceField>
-  <AmoebaGeneralizedKirkwoodForce solventDielectric="78.3" soluteDielectric="1.0" dielectricOffset="0.009" includeCavityTerm="1" probeRadius="0.14" surfaceAreaFactor="-170.351730663">
-   <GeneralizedKirkwood type="1" charge="-0.22620" shct="0.79"  /> 
-   <GeneralizedKirkwood type="2" charge="-0.15245" shct="0.72"  /> 
-   <GeneralizedKirkwood type="2" charge="-0.26188" shct="0.72"  /> 
-   <GeneralizedKirkwood type="2" charge="-0.00060" shct="0.72"  /> 
-   <GeneralizedKirkwood type="3" charge="0.84374" shct="0.72"  /> 
-   <GeneralizedKirkwood type="4" charge="0.12752" shct="0.85"  /> 
-   <GeneralizedKirkwood type="5" charge="-0.73597" shct="0.85"  /> 
-   <GeneralizedKirkwood type="6" charge="0.07168" shct="0.85"  /> 
-   <GeneralizedKirkwood type="7" charge="-0.14985" shct="0.79"  /> 
-   <GeneralizedKirkwood type="7" charge="-0.07700" shct="0.79"  /> 
-   <GeneralizedKirkwood type="7" charge="-0.14168" shct="0.79"  /> 
-   <GeneralizedKirkwood type="8" charge="-0.21238" shct="0.72"  /> 
-   <GeneralizedKirkwood type="8" charge="-0.35660" shct="0.72"  /> 
-   <GeneralizedKirkwood type="8" charge="0.01822" shct="0.72"  /> 
-   <GeneralizedKirkwood type="9" charge="0.85068" shct="0.72"  /> 
-   <GeneralizedKirkwood type="9" charge="0.86830" shct="0.72"  /> 
-   <GeneralizedKirkwood type="9" charge="0.86623" shct="0.72"  /> 
-   <GeneralizedKirkwood type="10" charge="0.12992" shct="0.85"  /> 
-   <GeneralizedKirkwood type="10" charge="0.13014" shct="0.85"  /> 
-   <GeneralizedKirkwood type="10" charge="0.13192" shct="0.85"  /> 
-   <GeneralizedKirkwood type="11" charge="-0.77770" shct="0.85"  /> 
-   <GeneralizedKirkwood type="11" charge="-0.78568" shct="0.85"  /> 
-   <GeneralizedKirkwood type="11" charge="-0.77214" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08544" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08646" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.09030" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.08972" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.09885" shct="0.85"  /> 
-   <GeneralizedKirkwood type="12" charge="0.09273" shct="0.85"  /> 
-   <GeneralizedKirkwood type="13" charge="-0.15440" shct="0.72"  /> 
-   <GeneralizedKirkwood type="14" charge="0.07484" shct="0.85"  /> 
-   <GeneralizedKirkwood type="15" charge="0.01482" shct="0.72"  /> 
-   <GeneralizedKirkwood type="16" charge="0.06618" shct="0.85"  /> 
-   <GeneralizedKirkwood type="17" charge="-0.17773" shct="0.72"  /> 
-   <GeneralizedKirkwood type="18" charge="0.05743" shct="0.85"  /> 
-   <GeneralizedKirkwood type="19" charge="-0.05880" shct="0.72"  /> 
-   <GeneralizedKirkwood type="21" charge="0.02425" shct="0.72"  /> 
-   <GeneralizedKirkwood type="23" charge="-0.11850" shct="0.72"  /> 
-   <GeneralizedKirkwood type="20" charge="0.04696" shct="0.85"  /> 
-   <GeneralizedKirkwood type="22" charge="0.00841" shct="0.85"  /> 
-   <GeneralizedKirkwood type="24" charge="0.03989" shct="0.85"  /> 
-   <GeneralizedKirkwood type="25" charge="0.02186" shct="0.72"  /> 
-   <GeneralizedKirkwood type="26" charge="0.06618" shct="0.85"  /> 
-   <GeneralizedKirkwood type="27" charge="-0.17773" shct="0.72"  /> 
-   <GeneralizedKirkwood type="28" charge="0.05743" shct="0.85"  /> 
-   <GeneralizedKirkwood type="29" charge="-0.13782" shct="0.72"  /> 
-   <GeneralizedKirkwood type="30" charge="0.05838" shct="0.85"  /> 
-   <GeneralizedKirkwood type="31" charge="-0.16689" shct="0.72"  /> 
-   <GeneralizedKirkwood type="32" charge="0.05849" shct="0.85"  /> 
-   <GeneralizedKirkwood type="33" charge="-0.31424" shct="0.72"  /> 
-   <GeneralizedKirkwood type="33" charge="-0.44882" shct="0.72"  /> 
-   <GeneralizedKirkwood type="33" charge="-0.01057" shct="0.72"  /> 
-   <GeneralizedKirkwood type="34" charge="0.18740" shct="0.72"  /> 
-   <GeneralizedKirkwood type="35" charge="0.01982" shct="0.85"  /> 
-   <GeneralizedKirkwood type="36" charge="-0.41295" shct="0.85"  /> 
-   <GeneralizedKirkwood type="37" charge="0.27409" shct="0.85"  /> 
-   <GeneralizedKirkwood type="38" charge="0.17854" shct="0.72"  /> 
-   <GeneralizedKirkwood type="39" charge="0.02134" shct="0.85"  /> 
-   <GeneralizedKirkwood type="40" charge="-0.15307" shct="0.72"  /> 
-   <GeneralizedKirkwood type="41" charge="0.07192" shct="0.85"  /> 
-   <GeneralizedKirkwood type="42" charge="-0.39938" shct="0.85"  /> 
-   <GeneralizedKirkwood type="43" charge="0.22557" shct="0.85"  /> 
-   <GeneralizedKirkwood type="44" charge="-0.24710" shct="0.72"  /> 
-   <GeneralizedKirkwood type="44" charge="-0.37021" shct="0.72"  /> 
-   <GeneralizedKirkwood type="44" charge="-0.00633" shct="0.72"  /> 
-   <GeneralizedKirkwood type="45" charge="-0.15207" shct="0.72"  /> 
-   <GeneralizedKirkwood type="45" charge="-0.15207" shct="0.72"  /> 
-   <GeneralizedKirkwood type="46" charge="0.07591" shct="0.85"  /> 
-   <GeneralizedKirkwood type="47" charge="-0.04712" shct="0.96"  /> 
-   <GeneralizedKirkwood type="48" charge="0.11129" shct="0.85"  /> 
-   <GeneralizedKirkwood type="49" charge="0.06417" shct="0.96"  /> 
-   <GeneralizedKirkwood type="50" charge="0.02101" shct="0.72"  /> 
-   <GeneralizedKirkwood type="51" charge="0.00952" shct="0.85"  /> 
-   <GeneralizedKirkwood type="52" charge="-0.97001" shct="0.96"  /> 
-   <GeneralizedKirkwood type="53" charge="-0.11893" shct="0.79"  /> 
-   <GeneralizedKirkwood type="54" charge="-0.38385" shct="0.72"  /> 
-   <GeneralizedKirkwood type="54" charge="-0.44047" shct="0.72"  /> 
-   <GeneralizedKirkwood type="55" charge="0.98874" shct="0.72"  /> 
-   <GeneralizedKirkwood type="56" charge="-0.82816" shct="0.85"  /> 
-   <GeneralizedKirkwood type="57" charge="0.10715" shct="0.85"  /> 
-   <GeneralizedKirkwood type="58" charge="-0.09495" shct="0.72"  /> 
-   <GeneralizedKirkwood type="58" charge="-0.12680" shct="0.72"  /> 
-   <GeneralizedKirkwood type="59" charge="0.09275" shct="0.85"  /> 
-   <GeneralizedKirkwood type="59" charge="0.11510" shct="0.85"  /> 
-   <GeneralizedKirkwood type="60" charge="-0.11624" shct="0.72"  /> 
-   <GeneralizedKirkwood type="60" charge="-0.14420" shct="0.72"  /> 
-   <GeneralizedKirkwood type="61" charge="0.09217" shct="0.85"  /> 
-   <GeneralizedKirkwood type="61" charge="0.12121" shct="0.85"  /> 
-   <GeneralizedKirkwood type="62" charge="-0.02016" shct="0.72"  /> 
-   <GeneralizedKirkwood type="63" charge="0.04828" shct="0.85"  /> 
-   <GeneralizedKirkwood type="64" charge="-0.01895" shct="0.72"  /> 
-   <GeneralizedKirkwood type="65" charge="0.06839" shct="0.85"  /> 
-   <GeneralizedKirkwood type="66" charge="-0.06275" shct="0.72"  /> 
-   <GeneralizedKirkwood type="67" charge="0.00890" shct="0.72"  /> 
-   <GeneralizedKirkwood type="68" charge="0.02818" shct="0.85"  /> 
-   <GeneralizedKirkwood type="69" charge="-0.05402" shct="0.72"  /> 
-   <GeneralizedKirkwood type="70" charge="0.03565" shct="0.85"  /> 
-   <GeneralizedKirkwood type="71" charge="-0.05826" shct="0.72"  /> 
-   <GeneralizedKirkwood type="72" charge="0.03588" shct="0.85"  /> 
-   <GeneralizedKirkwood type="73" charge="-0.08206" shct="0.72"  /> 
-   <GeneralizedKirkwood type="74" charge="0.09740" shct="0.85"  /> 
-   <GeneralizedKirkwood type="75" charge="-0.04596" shct="0.72"  /> 
-   <GeneralizedKirkwood type="76" charge="0.05037" shct="0.72"  /> 
-   <GeneralizedKirkwood type="77" charge="0.01064" shct="0.85"  /> 
-   <GeneralizedKirkwood type="78" charge="-0.12462" shct="0.72"  /> 
-   <GeneralizedKirkwood type="79" charge="0.01162" shct="0.85"  /> 
-   <GeneralizedKirkwood type="80" charge="0.34249" shct="0.72"  /> 
-   <GeneralizedKirkwood type="81" charge="-0.46444" shct="0.85"  /> 
-   <GeneralizedKirkwood type="82" charge="0.22927" shct="0.85"  /> 
-   <GeneralizedKirkwood type="83" charge="-0.14811" shct="0.72"  /> 
-   <GeneralizedKirkwood type="84" charge="0.08044" shct="0.85"  /> 
-   <GeneralizedKirkwood type="85" charge="-0.21500" shct="0.72"  /> 
-   <GeneralizedKirkwood type="86" charge="-0.04459" shct="0.72"  /> 
-   <GeneralizedKirkwood type="87" charge="0.01849" shct="0.85"  /> 
-   <GeneralizedKirkwood type="88" charge="-0.15802" shct="0.72"  /> 
-   <GeneralizedKirkwood type="89" charge="-0.03234" shct="0.85"  /> 
-   <GeneralizedKirkwood type="90" charge="0.62054" shct="0.72"  /> 
-   <GeneralizedKirkwood type="91" charge="-0.91150" shct="0.85"  /> 
-   <GeneralizedKirkwood type="92" charge="-0.16392" shct="0.72"  /> 
-   <GeneralizedKirkwood type="93" charge="0.09307" shct="0.85"  /> 
-   <GeneralizedKirkwood type="94" charge="-0.10518" shct="0.72"  /> 
-   <GeneralizedKirkwood type="95" charge="0.06363" shct="0.72"  /> 
-   <GeneralizedKirkwood type="96" charge="0.02858" shct="0.85"  /> 
-   <GeneralizedKirkwood type="97" charge="0.00670" shct="0.72"  /> 
-   <GeneralizedKirkwood type="98" charge="-0.04014" shct="0.79"  /> 
-   <GeneralizedKirkwood type="99" charge="0.08943" shct="0.85"  /> 
-   <GeneralizedKirkwood type="100" charge="0.27213" shct="0.72"  /> 
-   <GeneralizedKirkwood type="101" charge="-0.21794" shct="0.72"  /> 
-   <GeneralizedKirkwood type="102" charge="0.00801" shct="0.85"  /> 
-   <GeneralizedKirkwood type="103" charge="-0.11761" shct="0.72"  /> 
-   <GeneralizedKirkwood type="104" charge="0.00440" shct="0.85"  /> 
-   <GeneralizedKirkwood type="105" charge="0.10201" shct="0.72"  /> 
-   <GeneralizedKirkwood type="106" charge="0.00775" shct="0.85"  /> 
-   <GeneralizedKirkwood type="107" charge="-0.05764" shct="0.72"  /> 
-   <GeneralizedKirkwood type="108" charge="0.00377" shct="0.85"  /> 
-   <GeneralizedKirkwood type="109" charge="-0.04139" shct="0.72"  /> 
-   <GeneralizedKirkwood type="110" charge="0.12549" shct="0.85"  /> 
-   <GeneralizedKirkwood type="111" charge="0.11571" shct="0.72"  /> 
-   <GeneralizedKirkwood type="112" charge="-0.06925" shct="0.79"  /> 
-   <GeneralizedKirkwood type="113" charge="0.21415" shct="0.85"  /> 
-   <GeneralizedKirkwood type="114" charge="0.01888" shct="0.72"  /> 
-   <GeneralizedKirkwood type="115" charge="0.06585" shct="0.85"  /> 
-   <GeneralizedKirkwood type="116" charge="0.36202" shct="0.72"  /> 
-   <GeneralizedKirkwood type="117" charge="0.08951" shct="0.85"  /> 
-   <GeneralizedKirkwood type="118" charge="-0.07427" shct="0.79"  /> 
-   <GeneralizedKirkwood type="119" charge="0.13793" shct="0.85"  /> 
-   <GeneralizedKirkwood type="120" charge="-0.07944" shct="0.72"  /> 
-   <GeneralizedKirkwood type="121" charge="0.08230" shct="0.85"  /> 
-   <GeneralizedKirkwood type="122" charge="0.02915" shct="0.72"  /> 
-   <GeneralizedKirkwood type="123" charge="-0.20791" shct="0.79"  /> 
-   <GeneralizedKirkwood type="124" charge="0.22642" shct="0.85"  /> 
-   <GeneralizedKirkwood type="125" charge="0.02645" shct="0.72"  /> 
-   <GeneralizedKirkwood type="126" charge="0.05039" shct="0.85"  /> 
-   <GeneralizedKirkwood type="127" charge="0.27710" shct="0.72"  /> 
-   <GeneralizedKirkwood type="128" charge="0.07110" shct="0.85"  /> 
-   <GeneralizedKirkwood type="129" charge="-0.48774" shct="0.79"  /> 
-   <GeneralizedKirkwood type="130" charge="-0.07191" shct="0.72"  /> 
-   <GeneralizedKirkwood type="131" charge="0.09317" shct="0.85"  /> 
-   <GeneralizedKirkwood type="132" charge="0.21765" shct="0.72"  /> 
-   <GeneralizedKirkwood type="133" charge="-0.58085" shct="0.79"  /> 
-   <GeneralizedKirkwood type="134" charge="-0.05319" shct="0.72"  /> 
-   <GeneralizedKirkwood type="135" charge="0.03309" shct="0.85"  /> 
-   <GeneralizedKirkwood type="136" charge="0.31901" shct="0.72"  /> 
-   <GeneralizedKirkwood type="137" charge="0.03395" shct="0.85"  /> 
-   <GeneralizedKirkwood type="138" charge="-0.11867" shct="0.79"  /> 
-   <GeneralizedKirkwood type="139" charge="0.10470" shct="0.85"  /> 
-   <GeneralizedKirkwood type="140" charge="-0.33040" shct="0.72"  /> 
-   <GeneralizedKirkwood type="141" charge="0.04893" shct="0.85"  /> 
-   <GeneralizedKirkwood type="142" charge="1.01644" shct="0.72"  /> 
-   <GeneralizedKirkwood type="143" charge="-0.85689" shct="0.85"  /> 
-   <GeneralizedKirkwood type="144" charge="-0.09907" shct="0.72"  /> 
-   <GeneralizedKirkwood type="144" charge="-0.20674" shct="0.72"  /> 
-   <GeneralizedKirkwood type="145" charge="0.09104" shct="0.85"  /> 
-   <GeneralizedKirkwood type="146" charge="0.80267" shct="0.72"  /> 
-   <GeneralizedKirkwood type="147" charge="-0.61184" shct="0.85"  /> 
-   <GeneralizedKirkwood type="148" charge="-0.44451" shct="0.85"  /> 
-   <GeneralizedKirkwood type="149" charge="0.24079" shct="0.85"  /> 
-   <GeneralizedKirkwood type="150" charge="-0.15072" shct="0.72"  /> 
-   <GeneralizedKirkwood type="150" charge="-0.14858" shct="0.72"  /> 
-   <GeneralizedKirkwood type="151" charge="0.09679" shct="0.85"  /> 
-   <GeneralizedKirkwood type="152" charge="0.76960" shct="0.72"  /> 
-   <GeneralizedKirkwood type="153" charge="-0.72950" shct="0.85"  /> 
-   <GeneralizedKirkwood type="154" charge="-0.38020" shct="0.79"  /> 
-   <GeneralizedKirkwood type="155" charge="0.18368" shct="0.85"  /> 
-   <GeneralizedKirkwood type="156" charge="0.00410" shct="0.72"  /> 
-   <GeneralizedKirkwood type="157" charge="0.04082" shct="0.85"  /> 
-   <GeneralizedKirkwood type="158" charge="-0.36014" shct="0.72"  /> 
-   <GeneralizedKirkwood type="159" charge="-0.00356" shct="0.85"  /> 
-   <GeneralizedKirkwood type="160" charge="1.14596" shct="0.72"  /> 
-   <GeneralizedKirkwood type="161" charge="-0.89716" shct="0.85"  /> 
-   <GeneralizedKirkwood type="162" charge="-0.07714" shct="0.72"  /> 
-   <GeneralizedKirkwood type="163" charge="0.09378" shct="0.85"  /> 
-   <GeneralizedKirkwood type="164" charge="-0.12108" shct="0.72"  /> 
-   <GeneralizedKirkwood type="165" charge="0.05947" shct="0.85"  /> 
-   <GeneralizedKirkwood type="166" charge="-0.02949" shct="0.72"  /> 
-   <GeneralizedKirkwood type="167" charge="0.06373" shct="0.85"  /> 
-   <GeneralizedKirkwood type="168" charge="-0.07285" shct="0.72"  /> 
-   <GeneralizedKirkwood type="169" charge="0.00251" shct="0.85"  /> 
-   <GeneralizedKirkwood type="170" charge="0.06969" shct="0.96"  /> 
-   <GeneralizedKirkwood type="171" charge="-0.15553" shct="0.72"  /> 
-   <GeneralizedKirkwood type="172" charge="0.04194" shct="0.85"  /> 
-   <GeneralizedKirkwood type="173" charge="0.01417" shct="0.72"  /> 
-   <GeneralizedKirkwood type="174" charge="0.06168" shct="0.85"  /> 
-   <GeneralizedKirkwood type="175" charge="-0.23796" shct="0.72"  /> 
-   <GeneralizedKirkwood type="175" charge="-0.14658" shct="0.72"  /> 
-   <GeneralizedKirkwood type="176" charge="0.07914" shct="0.85"  /> 
-   <GeneralizedKirkwood type="177" charge="-0.12517" shct="0.72"  /> 
-   <GeneralizedKirkwood type="178" charge="0.07947" shct="0.85"  /> 
-   <GeneralizedKirkwood type="179" charge="-0.00637" shct="0.72"  /> 
-   <GeneralizedKirkwood type="180" charge="0.10424" shct="0.85"  /> 
-   <GeneralizedKirkwood type="181" charge="0.10679" shct="0.79"  /> 
-   <GeneralizedKirkwood type="182" charge="0.19274" shct="0.85"  /> 
-   <GeneralizedKirkwood type="183" charge="-0.11624" shct="0.72"  /> 
-   <GeneralizedKirkwood type="184" charge="0.05672" shct="0.85"  /> 
-   <GeneralizedKirkwood type="185" charge="0.11242" shct="0.72"  /> 
-   <GeneralizedKirkwood type="186" charge="0.02390" shct="0.85"  /> 
-   <GeneralizedKirkwood type="187" charge="-0.30941" shct="0.79"  /> 
-   <GeneralizedKirkwood type="188" charge="0.08213" shct="0.85"  /> 
-   <GeneralizedKirkwood type="189" charge="-0.05810" shct="0.72"  /> 
-   <GeneralizedKirkwood type="190" charge="0.07332" shct="0.85"  /> 
-   <GeneralizedKirkwood type="191" charge="-0.12460" shct="0.72"  /> 
-   <GeneralizedKirkwood type="192" charge="0.10013" shct="0.85"  /> 
-   <GeneralizedKirkwood type="193" charge="-0.00504" shct="0.72"  /> 
-   <GeneralizedKirkwood type="194" charge="0.06386" shct="0.85"  /> 
-   <GeneralizedKirkwood type="195" charge="-0.26789" shct="0.79"  /> 
-   <GeneralizedKirkwood type="196" charge="0.14895" shct="0.85"  /> 
-   <GeneralizedKirkwood type="197" charge="0.90680" shct="0.72"  /> 
-   <GeneralizedKirkwood type="198" charge="-0.30169" shct="0.79"  /> 
-   <GeneralizedKirkwood type="199" charge="0.14969" shct="0.85"  /> 
-   <GeneralizedKirkwood type="200" charge="0.01417" shct="0.72"  /> 
-   <GeneralizedKirkwood type="201" charge="0.06168" shct="0.85"  /> 
-   <GeneralizedKirkwood type="202" charge="-0.12127" shct="0.72"  /> 
-   <GeneralizedKirkwood type="203" charge="0.08337" shct="0.85"  /> 
-   <GeneralizedKirkwood type="204" charge="-0.00637" shct="0.72"  /> 
-   <GeneralizedKirkwood type="205" charge="0.10424" shct="0.85"  /> 
-   <GeneralizedKirkwood type="206" charge="0.10679" shct="0.79"  /> 
-   <GeneralizedKirkwood type="207" charge="0.19274" shct="0.85"  /> 
-   <GeneralizedKirkwood type="208" charge="-0.19495" shct="0.72"  /> 
-   <GeneralizedKirkwood type="209" charge="0.07797" shct="0.85"  /> 
-   <GeneralizedKirkwood type="210" charge="0.70656" shct="0.72"  /> 
-   <GeneralizedKirkwood type="211" charge="-0.74615" shct="0.85"  /> 
-   <GeneralizedKirkwood type="212" charge="-0.25887" shct="0.79"  /> 
-   <GeneralizedKirkwood type="213" charge="0.12975" shct="0.85"  /> 
-   <GeneralizedKirkwood type="214" charge="-0.27516" shct="0.79"  /> 
-   <GeneralizedKirkwood type="215" charge="0.12524" shct="0.85"  /> 
-   <GeneralizedKirkwood type="216" charge="-0.00902" shct="0.72"  /> 
-   <GeneralizedKirkwood type="217" charge="0.05319" shct="0.85"  /> 
-   <GeneralizedKirkwood type="218" charge="0.11164" shct="0.79"  /> 
-   <GeneralizedKirkwood type="218" charge="0.11164" shct="0.79"  /> 
-   <GeneralizedKirkwood type="218" charge="0.11164" shct="0.79"  /> 
-   <GeneralizedKirkwood type="218" charge="0.11164" shct="0.79"  /> 
-   <GeneralizedKirkwood type="219" charge="0.21240" shct="0.85"  /> 
-   <GeneralizedKirkwood type="219" charge="0.21240" shct="0.85"  /> 
-   <GeneralizedKirkwood type="219" charge="0.21240" shct="0.85"  /> 
-   <GeneralizedKirkwood type="219" charge="0.21240" shct="0.85"  /> 
-   <GeneralizedKirkwood type="220" charge="1.02669" shct="0.72"  /> 
-   <GeneralizedKirkwood type="221" charge="-0.90443" shct="0.85"  /> 
-   <GeneralizedKirkwood type="222" charge="0.07219" shct="0.79"  /> 
-   <GeneralizedKirkwood type="223" charge="0.23262" shct="0.85"  /> 
-   <GeneralizedKirkwood type="224" charge="-0.16428" shct="0.72"  /> 
-   <GeneralizedKirkwood type="225" charge="0.88092" shct="0.72"  /> 
-   <GeneralizedKirkwood type="226" charge="-0.74794" shct="0.85"  /> 
-   <GeneralizedKirkwood type="227" charge="0.12611" shct="0.85"  /> 
-   <GeneralizedKirkwood type="228" charge="-0.04567" shct="0.72"  /> 
-   <GeneralizedKirkwood type="229" charge="0.10559" shct="0.85"  /> 
-   <GeneralizedKirkwood type="230" charge="-0.12218" shct="0.79"  /> 
-   <GeneralizedKirkwood type="231" charge="0.10618" shct="0.72"  /> 
-   <GeneralizedKirkwood type="232" charge="0.03294" shct="0.72"  /> 
-   <GeneralizedKirkwood type="233" charge="-0.29610" shct="0.79"  /> 
-   <GeneralizedKirkwood type="234" charge="0.12091" shct="0.72"  /> 
-   <GeneralizedKirkwood type="235" charge="-0.28556" shct="0.79"  /> 
-   <GeneralizedKirkwood type="236" charge="0.09440" shct="0.72"  /> 
-   <GeneralizedKirkwood type="237" charge="-0.30504" shct="0.79"  /> 
-   <GeneralizedKirkwood type="238" charge="0.11701" shct="0.72"  /> 
-   <GeneralizedKirkwood type="239" charge="0.10987" shct="0.85"  /> 
-   <GeneralizedKirkwood type="240" charge="-0.28195" shct="0.79"  /> 
-   <GeneralizedKirkwood type="241" charge="0.23408" shct="0.85"  /> 
-   <GeneralizedKirkwood type="242" charge="0.23408" shct="0.85"  /> 
-   <GeneralizedKirkwood type="243" charge="0.12346" shct="0.85"  /> 
-   <GeneralizedKirkwood type="244" charge="-0.12896" shct="0.79"  /> 
-   <GeneralizedKirkwood type="245" charge="0.25555" shct="0.72"  /> 
-   <GeneralizedKirkwood type="246" charge="-0.29612" shct="0.79"  /> 
-   <GeneralizedKirkwood type="247" charge="0.10206" shct="0.72"  /> 
-   <GeneralizedKirkwood type="248" charge="-0.08445" shct="0.72"  /> 
-   <GeneralizedKirkwood type="249" charge="0.09191" shct="0.72"  /> 
-   <GeneralizedKirkwood type="250" charge="-0.42963" shct="0.85"  /> 
-   <GeneralizedKirkwood type="251" charge="-0.28291" shct="0.79"  /> 
-   <GeneralizedKirkwood type="252" charge="0.22172" shct="0.85"  /> 
-   <GeneralizedKirkwood type="253" charge="0.22172" shct="0.85"  /> 
-   <GeneralizedKirkwood type="254" charge="0.09969" shct="0.85"  /> 
-   <GeneralizedKirkwood type="255" charge="0.11152" shct="0.85"  /> 
-   <GeneralizedKirkwood type="256" charge="-0.12123" shct="0.79"  /> 
-   <GeneralizedKirkwood type="257" charge="0.11988" shct="0.72"  /> 
-   <GeneralizedKirkwood type="258" charge="0.05939" shct="0.72"  /> 
-   <GeneralizedKirkwood type="259" charge="-0.25427" shct="0.79"  /> 
-   <GeneralizedKirkwood type="260" charge="0.05900" shct="0.72"  /> 
-   <GeneralizedKirkwood type="261" charge="-0.32183" shct="0.79"  /> 
-   <GeneralizedKirkwood type="262" charge="0.16650" shct="0.72"  /> 
-   <GeneralizedKirkwood type="263" charge="-0.20681" shct="0.79"  /> 
-   <GeneralizedKirkwood type="264" charge="0.22596" shct="0.72"  /> 
-   <GeneralizedKirkwood type="265" charge="0.25827" shct="0.85"  /> 
-   <GeneralizedKirkwood type="266" charge="-0.30237" shct="0.79"  /> 
-   <GeneralizedKirkwood type="267" charge="0.22815" shct="0.85"  /> 
-   <GeneralizedKirkwood type="268" charge="0.22815" shct="0.85"  /> 
-   <GeneralizedKirkwood type="269" charge="-0.37776" shct="0.85"  /> 
-   <GeneralizedKirkwood type="270" charge="0.12107" shct="0.85"  /> 
-   <GeneralizedKirkwood type="271" charge="-0.15031" shct="0.79"  /> 
-   <GeneralizedKirkwood type="272" charge="0.25814" shct="0.72"  /> 
-   <GeneralizedKirkwood type="273" charge="-0.18494" shct="0.79"  /> 
-   <GeneralizedKirkwood type="274" charge="0.20468" shct="0.72"  /> 
-   <GeneralizedKirkwood type="275" charge="-0.10776" shct="0.72"  /> 
-   <GeneralizedKirkwood type="276" charge="0.02014" shct="0.72"  /> 
-   <GeneralizedKirkwood type="277" charge="-0.42656" shct="0.85"  /> 
-   <GeneralizedKirkwood type="278" charge="0.27344" shct="0.85"  /> 
-   <GeneralizedKirkwood type="279" charge="-0.40427" shct="0.85"  /> 
-   <GeneralizedKirkwood type="280" charge="-0.09928" shct="0.72"  /> 
-   <GeneralizedKirkwood type="281" charge="0.12786" shct="0.85"  /> 
-   <GeneralizedKirkwood type="282" charge="0.11524" shct="0.85"  /> 
-   <GeneralizedKirkwood type="283" charge="-0.11970" shct="0.79"  /> 
-   <GeneralizedKirkwood type="284" charge="0.26026" shct="0.72"  /> 
-   <GeneralizedKirkwood type="285" charge="-0.18804" shct="0.79"  /> 
-   <GeneralizedKirkwood type="286" charge="0.20160" shct="0.72"  /> 
-   <GeneralizedKirkwood type="287" charge="-0.05028" shct="0.72"  /> 
-   <GeneralizedKirkwood type="288" charge="0.07568" shct="0.72"  /> 
-   <GeneralizedKirkwood type="289" charge="-0.42057" shct="0.85"  /> 
-   <GeneralizedKirkwood type="290" charge="0.27376" shct="0.85"  /> 
-   <GeneralizedKirkwood type="291" charge="-0.38720" shct="0.85"  /> 
-   <GeneralizedKirkwood type="292" charge="0.12060" shct="0.85"  /> 
-   <GeneralizedKirkwood type="293" charge="0.11599" shct="0.85"  /> 
-   <GeneralizedKirkwood type="294" charge="-0.52467" shct="0.85"  /> 
-   <GeneralizedKirkwood type="295" charge="0.05049" shct="0.72"  /> 
-   <GeneralizedKirkwood type="296" charge="0.07556" shct="0.85"  /> 
-   <GeneralizedKirkwood type="297" charge="0.07556" shct="0.85"  /> 
-   <GeneralizedKirkwood type="298" charge="0.10423" shct="0.72"  /> 
-   <GeneralizedKirkwood type="299" charge="0.13028" shct="0.85"  /> 
-   <GeneralizedKirkwood type="300" charge="-0.40197" shct="0.85"  /> 
-   <GeneralizedKirkwood type="301" charge="0.06316" shct="0.72"  /> 
-   <GeneralizedKirkwood type="302" charge="0.18443" shct="0.85"  /> 
-   <GeneralizedKirkwood type="303" charge="0.10143" shct="0.72"  /> 
-   <GeneralizedKirkwood type="304" charge="0.11227" shct="0.85"  /> 
-   <GeneralizedKirkwood type="305" charge="0.06213" shct="0.72"  /> 
-   <GeneralizedKirkwood type="306" charge="0.11711" shct="0.85"  /> 
-   <GeneralizedKirkwood type="307" charge="-0.47891" shct="0.85"  /> 
-   <GeneralizedKirkwood type="308" charge="0.32469" shct="0.85"  /> 
-   <GeneralizedKirkwood type="309" charge="-0.53467" shct="0.85"  /> 
-   <GeneralizedKirkwood type="310" charge="-0.52467" shct="0.85"  /> 
-   <GeneralizedKirkwood type="311" charge="0.05049" shct="0.72"  /> 
-   <GeneralizedKirkwood type="312" charge="0.07556" shct="0.85"  /> 
-   <GeneralizedKirkwood type="313" charge="0.07556" shct="0.85"  /> 
-   <GeneralizedKirkwood type="314" charge="0.10423" shct="0.72"  /> 
-   <GeneralizedKirkwood type="315" charge="0.13028" shct="0.85"  /> 
-   <GeneralizedKirkwood type="316" charge="-0.40197" shct="0.85"  /> 
-   <GeneralizedKirkwood type="317" charge="0.06316" shct="0.72"  /> 
-   <GeneralizedKirkwood type="318" charge="0.18443" shct="0.85"  /> 
-   <GeneralizedKirkwood type="319" charge="0.10143" shct="0.72"  /> 
-   <GeneralizedKirkwood type="320" charge="0.11227" shct="0.85"  /> 
-   <GeneralizedKirkwood type="321" charge="0.06213" shct="0.72"  /> 
-   <GeneralizedKirkwood type="322" charge="0.11711" shct="0.85"  /> 
-   <GeneralizedKirkwood type="323" charge="-0.47891" shct="0.85"  /> 
-   <GeneralizedKirkwood type="324" charge="0.32469" shct="0.85"  /> 
-   <GeneralizedKirkwood type="325" charge="-0.53467" shct="0.85"  /> 
-   <GeneralizedKirkwood type="326" charge="-0.52467" shct="0.85"  /> 
-   <GeneralizedKirkwood type="327" charge="0.00960" shct="0.72"  /> 
-   <GeneralizedKirkwood type="328" charge="0.07556" shct="0.85"  /> 
-   <GeneralizedKirkwood type="329" charge="0.07556" shct="0.85"  /> 
-   <GeneralizedKirkwood type="330" charge="0.15826" shct="0.72"  /> 
-   <GeneralizedKirkwood type="331" charge="0.12800" shct="0.85"  /> 
-   <GeneralizedKirkwood type="332" charge="-0.38842" shct="0.85"  /> 
-   <GeneralizedKirkwood type="333" charge="0.06720" shct="0.72"  /> 
-   <GeneralizedKirkwood type="334" charge="0.12287" shct="0.85"  /> 
-   <GeneralizedKirkwood type="335" charge="0.07060" shct="0.72"  /> 
-   <GeneralizedKirkwood type="336" charge="0.10297" shct="0.85"  /> 
-   <GeneralizedKirkwood type="337" charge="-0.10232" shct="0.72"  /> 
-   <GeneralizedKirkwood type="338" charge="0.10225" shct="0.85"  /> 
-   <GeneralizedKirkwood type="339" charge="0.09833" shct="0.85"  /> 
-   <GeneralizedKirkwood type="340" charge="-0.53467" shct="0.85"  /> 
-   <GeneralizedKirkwood type="341" charge="-0.52467" shct="0.85"  /> 
-   <GeneralizedKirkwood type="342" charge="0.00960" shct="0.72"  /> 
-   <GeneralizedKirkwood type="343" charge="0.07556" shct="0.85"  /> 
-   <GeneralizedKirkwood type="344" charge="0.07556" shct="0.85"  /> 
-   <GeneralizedKirkwood type="345" charge="0.15826" shct="0.72"  /> 
-   <GeneralizedKirkwood type="346" charge="0.12800" shct="0.85"  /> 
-   <GeneralizedKirkwood type="347" charge="-0.38842" shct="0.85"  /> 
-   <GeneralizedKirkwood type="348" charge="0.06720" shct="0.72"  /> 
-   <GeneralizedKirkwood type="349" charge="0.12287" shct="0.85"  /> 
-   <GeneralizedKirkwood type="350" charge="0.07060" shct="0.72"  /> 
-   <GeneralizedKirkwood type="351" charge="0.10297" shct="0.85"  /> 
-   <GeneralizedKirkwood type="352" charge="-0.10232" shct="0.72"  /> 
-   <GeneralizedKirkwood type="353" charge="0.10225" shct="0.85"  /> 
-   <GeneralizedKirkwood type="354" charge="0.09833" shct="0.85"  /> 
-   <GeneralizedKirkwood type="355" charge="-0.53467" shct="0.85"  /> 
-   <GeneralizedKirkwood type="356" charge="0.90070" shct="0.86"  /> 
-   <GeneralizedKirkwood type="357" charge="-0.62196" shct="0.85"  /> 
-   <GeneralizedKirkwood type="358" charge="-0.53491" shct="0.85"  /> 
-   <GeneralizedKirkwood type="358" charge="-0.53491" shct="0.85"  /> 
-   <GeneralizedKirkwood type="359" charge="0.34636" shct="0.85"  /> 
-   <GeneralizedKirkwood type="359" charge="0.34636" shct="0.85"  /> 
-   <GeneralizedKirkwood type="359" charge="0.34636" shct="0.85"  /> 
-   <GeneralizedKirkwood type="360" charge="-0.58166" shct="0.85"  /> 
-   <GeneralizedKirkwood type="360" charge="-0.58166" shct="0.85"  /> 
-   <GeneralizedKirkwood type="361" charge="1.02369" shct="0.86"  /> 
-   <GeneralizedKirkwood type="362" charge="-0.87686" shct="0.85"  /> 
-   <GeneralizedKirkwood type="363" charge="-0.50819" shct="0.85"  /> 
-   <GeneralizedKirkwood type="363" charge="-0.50819" shct="0.85"  /> 
-   <GeneralizedKirkwood type="364" charge="0.29418" shct="0.85"  /> 
-   <GeneralizedKirkwood type="364" charge="0.29418" shct="0.85"  /> 
-   <GeneralizedKirkwood type="365" charge="-0.60712" shct="0.85"  /> 
-   <GeneralizedKirkwood type="365" charge="-0.60712" shct="0.85"  /> 
-   <GeneralizedKirkwood type="366" charge="1.02369" shct="0.86"  /> 
-   <GeneralizedKirkwood type="367" charge="-0.87686" shct="0.85"  /> 
-   <GeneralizedKirkwood type="368" charge="0.90070" shct="0.86"  /> 
-   <GeneralizedKirkwood type="369" charge="-0.62196" shct="0.85"  /> 
-   <GeneralizedKirkwood type="370" charge="-0.53491" shct="0.85"  /> 
-   <GeneralizedKirkwood type="370" charge="-0.53491" shct="0.85"  /> 
-   <GeneralizedKirkwood type="371" charge="0.34636" shct="0.85"  /> 
-   <GeneralizedKirkwood type="371" charge="0.34636" shct="0.85"  /> 
-   <GeneralizedKirkwood type="372" charge="-0.58166" shct="0.85"  /> 
-   <GeneralizedKirkwood type="372" charge="-0.58166" shct="0.85"  /> 
-   <GeneralizedKirkwood type="373" charge="1.02369" shct="0.86"  /> 
-   <GeneralizedKirkwood type="374" charge="-0.87686" shct="0.85"  /> 
-   <GeneralizedKirkwood type="375" charge="-0.50819" shct="0.85"  /> 
-   <GeneralizedKirkwood type="375" charge="-0.50819" shct="0.85"  /> 
-   <GeneralizedKirkwood type="376" charge="0.29418" shct="0.85"  /> 
-   <GeneralizedKirkwood type="376" charge="0.29418" shct="0.85"  /> 
-   <GeneralizedKirkwood type="377" charge="-0.60712" shct="0.85"  /> 
-   <GeneralizedKirkwood type="377" charge="-0.60712" shct="0.85"  /> 
-   <GeneralizedKirkwood type="378" charge="1.02369" shct="0.86"  /> 
-   <GeneralizedKirkwood type="379" charge="-0.87686" shct="0.85"  /> 
-   <GeneralizedKirkwood type="380" charge="-0.51966" shct="0.85"  /> 
-   <GeneralizedKirkwood type="381" charge="0.25983" shct="0.85"  /> 
-   <GeneralizedKirkwood type="382" charge="1.00000" shct="0.8"  /> 
-   <GeneralizedKirkwood type="383" charge="1.00000" shct="0.8"  /> 
-   <GeneralizedKirkwood type="384" charge="1.00000" shct="0.8"  /> 
-   <GeneralizedKirkwood type="385" charge="1.00000" shct="0.8"  /> 
-   <GeneralizedKirkwood type="386" charge="1.00000" shct="0.8"  /> 
-   <GeneralizedKirkwood type="387" charge="2.00000" shct="0.8"  /> 
-   <GeneralizedKirkwood type="388" charge="2.00000" shct="0.8"  /> 
-   <GeneralizedKirkwood type="389" charge="2.00000" shct="0.8"  /> 
-   <GeneralizedKirkwood type="390" charge="-1.00000" shct="0.8"  /> 
+ <Info>
+  <Source>amoebabio09.prm</Source>
+  <DateGenerated>2015-11-19</DateGenerated>
+  <Reference></Reference>
+ </Info>
+ <AmoebaGeneralizedKirkwoodForce solventDielectric="78.3" soluteDielectric="1.0" includeCavityTerm="1" probeRadius="0.14" surfaceAreaFactor="-170.351730663">
+  <GeneralizedKirkwood type="1" charge="-0.22620" shct="0.79"  /> 
+  <GeneralizedKirkwood type="2" charge="-0.15245" shct="0.72"  /> 
+  <GeneralizedKirkwood type="2" charge="-0.26188" shct="0.72"  /> 
+  <GeneralizedKirkwood type="2" charge="-0.06378" shct="0.72"  /> 
+  <GeneralizedKirkwood type="2" charge="-0.00060" shct="0.72"  /> 
+  <GeneralizedKirkwood type="2" charge="-0.15431" shct="0.72"  /> 
+  <GeneralizedKirkwood type="3" charge="0.84374" shct="0.72"  /> 
+  <GeneralizedKirkwood type="4" charge="0.12752" shct="0.85"  /> 
+  <GeneralizedKirkwood type="5" charge="-0.73597" shct="0.85"  /> 
+  <GeneralizedKirkwood type="6" charge="0.07168" shct="0.85"  /> 
+  <GeneralizedKirkwood type="7" charge="-0.14985" shct="0.79"  /> 
+  <GeneralizedKirkwood type="7" charge="-0.07700" shct="0.79"  /> 
+  <GeneralizedKirkwood type="7" charge="-0.14168" shct="0.79"  /> 
+  <GeneralizedKirkwood type="8" charge="-0.21238" shct="0.72"  /> 
+  <GeneralizedKirkwood type="8" charge="-0.35660" shct="0.72"  /> 
+  <GeneralizedKirkwood type="8" charge="-0.15850" shct="0.72"  /> 
+  <GeneralizedKirkwood type="8" charge="0.01822" shct="0.72"  /> 
+  <GeneralizedKirkwood type="8" charge="-0.13549" shct="0.72"  /> 
+  <GeneralizedKirkwood type="9" charge="0.85068" shct="0.72"  /> 
+  <GeneralizedKirkwood type="9" charge="0.86830" shct="0.72"  /> 
+  <GeneralizedKirkwood type="9" charge="0.86623" shct="0.72"  /> 
+  <GeneralizedKirkwood type="10" charge="0.12992" shct="0.85"  /> 
+  <GeneralizedKirkwood type="10" charge="0.13014" shct="0.85"  /> 
+  <GeneralizedKirkwood type="10" charge="0.13192" shct="0.85"  /> 
+  <GeneralizedKirkwood type="11" charge="-0.77770" shct="0.85"  /> 
+  <GeneralizedKirkwood type="11" charge="-0.78568" shct="0.85"  /> 
+  <GeneralizedKirkwood type="11" charge="-0.77214" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08921" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08544" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08646" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.09030" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.08972" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.09885" shct="0.85"  /> 
+  <GeneralizedKirkwood type="12" charge="0.09273" shct="0.85"  /> 
+  <GeneralizedKirkwood type="13" charge="-0.15440" shct="0.72"  /> 
+  <GeneralizedKirkwood type="14" charge="0.07484" shct="0.85"  /> 
+  <GeneralizedKirkwood type="15" charge="0.01482" shct="0.72"  /> 
+  <GeneralizedKirkwood type="16" charge="0.06618" shct="0.85"  /> 
+  <GeneralizedKirkwood type="17" charge="-0.17773" shct="0.72"  /> 
+  <GeneralizedKirkwood type="18" charge="0.05743" shct="0.85"  /> 
+  <GeneralizedKirkwood type="19" charge="-0.05880" shct="0.72"  /> 
+  <GeneralizedKirkwood type="21" charge="0.02425" shct="0.72"  /> 
+  <GeneralizedKirkwood type="23" charge="-0.11850" shct="0.72"  /> 
+  <GeneralizedKirkwood type="20" charge="0.04696" shct="0.85"  /> 
+  <GeneralizedKirkwood type="22" charge="0.00841" shct="0.85"  /> 
+  <GeneralizedKirkwood type="24" charge="0.03989" shct="0.85"  /> 
+  <GeneralizedKirkwood type="25" charge="0.02186" shct="0.72"  /> 
+  <GeneralizedKirkwood type="26" charge="0.06618" shct="0.85"  /> 
+  <GeneralizedKirkwood type="27" charge="-0.17773" shct="0.72"  /> 
+  <GeneralizedKirkwood type="28" charge="0.05743" shct="0.85"  /> 
+  <GeneralizedKirkwood type="29" charge="-0.13782" shct="0.72"  /> 
+  <GeneralizedKirkwood type="30" charge="0.05838" shct="0.85"  /> 
+  <GeneralizedKirkwood type="31" charge="-0.16689" shct="0.72"  /> 
+  <GeneralizedKirkwood type="32" charge="0.05849" shct="0.85"  /> 
+  <GeneralizedKirkwood type="33" charge="-0.31424" shct="0.72"  /> 
+  <GeneralizedKirkwood type="33" charge="-0.44882" shct="0.72"  /> 
+  <GeneralizedKirkwood type="33" charge="-0.01057" shct="0.72"  /> 
+  <GeneralizedKirkwood type="33" charge="-0.16428" shct="0.72"  /> 
+  <GeneralizedKirkwood type="34" charge="0.18740" shct="0.72"  /> 
+  <GeneralizedKirkwood type="35" charge="0.01982" shct="0.85"  /> 
+  <GeneralizedKirkwood type="36" charge="-0.41295" shct="0.85"  /> 
+  <GeneralizedKirkwood type="37" charge="0.27409" shct="0.85"  /> 
+  <GeneralizedKirkwood type="38" charge="0.17854" shct="0.72"  /> 
+  <GeneralizedKirkwood type="39" charge="0.02134" shct="0.85"  /> 
+  <GeneralizedKirkwood type="40" charge="-0.15307" shct="0.72"  /> 
+  <GeneralizedKirkwood type="41" charge="0.07192" shct="0.85"  /> 
+  <GeneralizedKirkwood type="42" charge="-0.39938" shct="0.85"  /> 
+  <GeneralizedKirkwood type="43" charge="0.22557" shct="0.85"  /> 
+  <GeneralizedKirkwood type="44" charge="-0.24710" shct="0.72"  /> 
+  <GeneralizedKirkwood type="44" charge="-0.37021" shct="0.72"  /> 
+  <GeneralizedKirkwood type="44" charge="-0.00633" shct="0.72"  /> 
+  <GeneralizedKirkwood type="44" charge="-0.16004" shct="0.72"  /> 
+  <GeneralizedKirkwood type="45" charge="-0.15207" shct="0.72"  /> 
+  <GeneralizedKirkwood type="45" charge="-0.15207" shct="0.72"  /> 
+  <GeneralizedKirkwood type="46" charge="0.07591" shct="0.85"  /> 
+  <GeneralizedKirkwood type="47" charge="-0.04712" shct="0.96"  /> 
+  <GeneralizedKirkwood type="48" charge="0.11129" shct="0.85"  /> 
+  <GeneralizedKirkwood type="49" charge="0.06417" shct="0.96"  /> 
+  <GeneralizedKirkwood type="50" charge="0.02101" shct="0.72"  /> 
+  <GeneralizedKirkwood type="51" charge="0.00952" shct="0.85"  /> 
+  <GeneralizedKirkwood type="52" charge="-0.97001" shct="0.96"  /> 
+  <GeneralizedKirkwood type="53" charge="-0.11893" shct="0.79"  /> 
+  <GeneralizedKirkwood type="54" charge="-0.38385" shct="0.72"  /> 
+  <GeneralizedKirkwood type="54" charge="-0.44047" shct="0.72"  /> 
+  <GeneralizedKirkwood type="55" charge="0.98874" shct="0.72"  /> 
+  <GeneralizedKirkwood type="56" charge="-0.82816" shct="0.85"  /> 
+  <GeneralizedKirkwood type="57" charge="0.10715" shct="0.85"  /> 
+  <GeneralizedKirkwood type="58" charge="-0.09495" shct="0.72"  /> 
+  <GeneralizedKirkwood type="58" charge="-0.12680" shct="0.72"  /> 
+  <GeneralizedKirkwood type="59" charge="0.09275" shct="0.85"  /> 
+  <GeneralizedKirkwood type="59" charge="0.11510" shct="0.85"  /> 
+  <GeneralizedKirkwood type="60" charge="-0.11624" shct="0.72"  /> 
+  <GeneralizedKirkwood type="60" charge="-0.14420" shct="0.72"  /> 
+  <GeneralizedKirkwood type="61" charge="0.09217" shct="0.85"  /> 
+  <GeneralizedKirkwood type="61" charge="0.12121" shct="0.85"  /> 
+  <GeneralizedKirkwood type="62" charge="-0.02016" shct="0.72"  /> 
+  <GeneralizedKirkwood type="63" charge="0.04828" shct="0.85"  /> 
+  <GeneralizedKirkwood type="64" charge="-0.01895" shct="0.72"  /> 
+  <GeneralizedKirkwood type="65" charge="0.06839" shct="0.85"  /> 
+  <GeneralizedKirkwood type="66" charge="-0.06275" shct="0.72"  /> 
+  <GeneralizedKirkwood type="67" charge="0.00890" shct="0.72"  /> 
+  <GeneralizedKirkwood type="68" charge="0.02818" shct="0.85"  /> 
+  <GeneralizedKirkwood type="69" charge="-0.05402" shct="0.72"  /> 
+  <GeneralizedKirkwood type="70" charge="0.03565" shct="0.85"  /> 
+  <GeneralizedKirkwood type="71" charge="-0.05826" shct="0.72"  /> 
+  <GeneralizedKirkwood type="72" charge="0.03588" shct="0.85"  /> 
+  <GeneralizedKirkwood type="73" charge="-0.08206" shct="0.72"  /> 
+  <GeneralizedKirkwood type="74" charge="0.09740" shct="0.85"  /> 
+  <GeneralizedKirkwood type="75" charge="-0.04596" shct="0.72"  /> 
+  <GeneralizedKirkwood type="76" charge="0.05037" shct="0.72"  /> 
+  <GeneralizedKirkwood type="77" charge="0.01064" shct="0.85"  /> 
+  <GeneralizedKirkwood type="78" charge="-0.12462" shct="0.72"  /> 
+  <GeneralizedKirkwood type="79" charge="0.01162" shct="0.85"  /> 
+  <GeneralizedKirkwood type="80" charge="0.34249" shct="0.72"  /> 
+  <GeneralizedKirkwood type="81" charge="-0.46444" shct="0.85"  /> 
+  <GeneralizedKirkwood type="82" charge="0.22927" shct="0.85"  /> 
+  <GeneralizedKirkwood type="83" charge="-0.14811" shct="0.72"  /> 
+  <GeneralizedKirkwood type="84" charge="0.08044" shct="0.85"  /> 
+  <GeneralizedKirkwood type="85" charge="-0.21500" shct="0.72"  /> 
+  <GeneralizedKirkwood type="86" charge="-0.04459" shct="0.72"  /> 
+  <GeneralizedKirkwood type="87" charge="0.01849" shct="0.85"  /> 
+  <GeneralizedKirkwood type="88" charge="-0.15802" shct="0.72"  /> 
+  <GeneralizedKirkwood type="89" charge="-0.03234" shct="0.85"  /> 
+  <GeneralizedKirkwood type="90" charge="0.62054" shct="0.72"  /> 
+  <GeneralizedKirkwood type="91" charge="-0.91150" shct="0.85"  /> 
+  <GeneralizedKirkwood type="92" charge="-0.16392" shct="0.72"  /> 
+  <GeneralizedKirkwood type="93" charge="0.09307" shct="0.85"  /> 
+  <GeneralizedKirkwood type="94" charge="-0.10518" shct="0.72"  /> 
+  <GeneralizedKirkwood type="95" charge="0.06363" shct="0.72"  /> 
+  <GeneralizedKirkwood type="96" charge="0.02858" shct="0.85"  /> 
+  <GeneralizedKirkwood type="97" charge="0.00670" shct="0.72"  /> 
+  <GeneralizedKirkwood type="98" charge="-0.04014" shct="0.79"  /> 
+  <GeneralizedKirkwood type="99" charge="0.08943" shct="0.85"  /> 
+  <GeneralizedKirkwood type="100" charge="0.27213" shct="0.72"  /> 
+  <GeneralizedKirkwood type="101" charge="-0.21794" shct="0.72"  /> 
+  <GeneralizedKirkwood type="102" charge="0.00801" shct="0.85"  /> 
+  <GeneralizedKirkwood type="103" charge="-0.11761" shct="0.72"  /> 
+  <GeneralizedKirkwood type="104" charge="0.00440" shct="0.85"  /> 
+  <GeneralizedKirkwood type="105" charge="0.10201" shct="0.72"  /> 
+  <GeneralizedKirkwood type="106" charge="0.00775" shct="0.85"  /> 
+  <GeneralizedKirkwood type="107" charge="-0.05764" shct="0.72"  /> 
+  <GeneralizedKirkwood type="108" charge="0.00377" shct="0.85"  /> 
+  <GeneralizedKirkwood type="109" charge="-0.04139" shct="0.72"  /> 
+  <GeneralizedKirkwood type="110" charge="0.12549" shct="0.85"  /> 
+  <GeneralizedKirkwood type="111" charge="0.11571" shct="0.72"  /> 
+  <GeneralizedKirkwood type="112" charge="-0.06925" shct="0.79"  /> 
+  <GeneralizedKirkwood type="113" charge="0.21415" shct="0.85"  /> 
+  <GeneralizedKirkwood type="114" charge="0.01888" shct="0.72"  /> 
+  <GeneralizedKirkwood type="115" charge="0.06585" shct="0.85"  /> 
+  <GeneralizedKirkwood type="116" charge="0.36202" shct="0.72"  /> 
+  <GeneralizedKirkwood type="117" charge="0.08951" shct="0.85"  /> 
+  <GeneralizedKirkwood type="118" charge="-0.07427" shct="0.79"  /> 
+  <GeneralizedKirkwood type="119" charge="0.13793" shct="0.85"  /> 
+  <GeneralizedKirkwood type="120" charge="-0.07944" shct="0.72"  /> 
+  <GeneralizedKirkwood type="121" charge="0.08230" shct="0.85"  /> 
+  <GeneralizedKirkwood type="122" charge="0.02915" shct="0.72"  /> 
+  <GeneralizedKirkwood type="123" charge="-0.20791" shct="0.79"  /> 
+  <GeneralizedKirkwood type="124" charge="0.22642" shct="0.85"  /> 
+  <GeneralizedKirkwood type="125" charge="0.02645" shct="0.72"  /> 
+  <GeneralizedKirkwood type="126" charge="0.05039" shct="0.85"  /> 
+  <GeneralizedKirkwood type="127" charge="0.27710" shct="0.72"  /> 
+  <GeneralizedKirkwood type="128" charge="0.07110" shct="0.85"  /> 
+  <GeneralizedKirkwood type="129" charge="-0.48774" shct="0.79"  /> 
+  <GeneralizedKirkwood type="130" charge="-0.07191" shct="0.72"  /> 
+  <GeneralizedKirkwood type="131" charge="0.09317" shct="0.85"  /> 
+  <GeneralizedKirkwood type="132" charge="0.21765" shct="0.72"  /> 
+  <GeneralizedKirkwood type="133" charge="-0.58085" shct="0.79"  /> 
+  <GeneralizedKirkwood type="134" charge="-0.05319" shct="0.72"  /> 
+  <GeneralizedKirkwood type="135" charge="0.03309" shct="0.85"  /> 
+  <GeneralizedKirkwood type="136" charge="0.31901" shct="0.72"  /> 
+  <GeneralizedKirkwood type="137" charge="0.03395" shct="0.85"  /> 
+  <GeneralizedKirkwood type="138" charge="-0.11867" shct="0.79"  /> 
+  <GeneralizedKirkwood type="139" charge="0.10470" shct="0.85"  /> 
+  <GeneralizedKirkwood type="140" charge="-0.33040" shct="0.72"  /> 
+  <GeneralizedKirkwood type="141" charge="0.04893" shct="0.85"  /> 
+  <GeneralizedKirkwood type="142" charge="1.01644" shct="0.72"  /> 
+  <GeneralizedKirkwood type="143" charge="-0.85689" shct="0.85"  /> 
+  <GeneralizedKirkwood type="144" charge="-0.09907" shct="0.72"  /> 
+  <GeneralizedKirkwood type="145" charge="0.09104" shct="0.85"  /> 
+  <GeneralizedKirkwood type="146" charge="0.80267" shct="0.72"  /> 
+  <GeneralizedKirkwood type="147" charge="-0.61184" shct="0.85"  /> 
+  <GeneralizedKirkwood type="148" charge="-0.44451" shct="0.85"  /> 
+  <GeneralizedKirkwood type="149" charge="0.24079" shct="0.85"  /> 
+  <GeneralizedKirkwood type="150" charge="-0.15072" shct="0.72"  /> 
+  <GeneralizedKirkwood type="151" charge="0.09679" shct="0.85"  /> 
+  <GeneralizedKirkwood type="152" charge="0.76960" shct="0.72"  /> 
+  <GeneralizedKirkwood type="153" charge="-0.72950" shct="0.85"  /> 
+  <GeneralizedKirkwood type="154" charge="-0.38020" shct="0.79"  /> 
+  <GeneralizedKirkwood type="155" charge="0.18368" shct="0.85"  /> 
+  <GeneralizedKirkwood type="156" charge="0.00410" shct="0.72"  /> 
+  <GeneralizedKirkwood type="157" charge="0.04082" shct="0.85"  /> 
+  <GeneralizedKirkwood type="158" charge="-0.36014" shct="0.72"  /> 
+  <GeneralizedKirkwood type="159" charge="-0.00356" shct="0.85"  /> 
+  <GeneralizedKirkwood type="160" charge="1.14596" shct="0.72"  /> 
+  <GeneralizedKirkwood type="161" charge="-0.89716" shct="0.85"  /> 
+  <GeneralizedKirkwood type="162" charge="-0.07714" shct="0.72"  /> 
+  <GeneralizedKirkwood type="163" charge="0.09378" shct="0.85"  /> 
+  <GeneralizedKirkwood type="164" charge="-0.20674" shct="0.72"  /> 
+  <GeneralizedKirkwood type="165" charge="0.09104" shct="0.85"  /> 
+  <GeneralizedKirkwood type="166" charge="0.80267" shct="0.72"  /> 
+  <GeneralizedKirkwood type="167" charge="-0.61184" shct="0.85"  /> 
+  <GeneralizedKirkwood type="168" charge="-0.44451" shct="0.85"  /> 
+  <GeneralizedKirkwood type="169" charge="0.24079" shct="0.85"  /> 
+  <GeneralizedKirkwood type="170" charge="-0.12108" shct="0.72"  /> 
+  <GeneralizedKirkwood type="171" charge="0.05947" shct="0.85"  /> 
+  <GeneralizedKirkwood type="172" charge="-0.14858" shct="0.72"  /> 
+  <GeneralizedKirkwood type="173" charge="0.09679" shct="0.85"  /> 
+  <GeneralizedKirkwood type="174" charge="0.76960" shct="0.72"  /> 
+  <GeneralizedKirkwood type="175" charge="-0.72950" shct="0.85"  /> 
+  <GeneralizedKirkwood type="176" charge="-0.38020" shct="0.79"  /> 
+  <GeneralizedKirkwood type="177" charge="0.18368" shct="0.85"  /> 
+  <GeneralizedKirkwood type="178" charge="-0.02949" shct="0.72"  /> 
+  <GeneralizedKirkwood type="179" charge="0.06373" shct="0.85"  /> 
+  <GeneralizedKirkwood type="180" charge="-0.07285" shct="0.72"  /> 
+  <GeneralizedKirkwood type="181" charge="0.00251" shct="0.85"  /> 
+  <GeneralizedKirkwood type="182" charge="0.06969" shct="0.96"  /> 
+  <GeneralizedKirkwood type="183" charge="-0.15553" shct="0.72"  /> 
+  <GeneralizedKirkwood type="184" charge="0.04194" shct="0.85"  /> 
+  <GeneralizedKirkwood type="185" charge="0.01417" shct="0.72"  /> 
+  <GeneralizedKirkwood type="186" charge="0.06168" shct="0.85"  /> 
+  <GeneralizedKirkwood type="187" charge="-0.14658" shct="0.72"  /> 
+  <GeneralizedKirkwood type="188" charge="0.07914" shct="0.85"  /> 
+  <GeneralizedKirkwood type="189" charge="-0.12517" shct="0.72"  /> 
+  <GeneralizedKirkwood type="190" charge="0.07947" shct="0.85"  /> 
+  <GeneralizedKirkwood type="191" charge="-0.00637" shct="0.72"  /> 
+  <GeneralizedKirkwood type="192" charge="0.10424" shct="0.85"  /> 
+  <GeneralizedKirkwood type="193" charge="0.10679" shct="0.79"  /> 
+  <GeneralizedKirkwood type="194" charge="0.19274" shct="0.85"  /> 
+  <GeneralizedKirkwood type="195" charge="0.01417" shct="0.72"  /> 
+  <GeneralizedKirkwood type="196" charge="0.06168" shct="0.85"  /> 
+  <GeneralizedKirkwood type="197" charge="-0.23796" shct="0.72"  /> 
+  <GeneralizedKirkwood type="198" charge="0.07914" shct="0.85"  /> 
+  <GeneralizedKirkwood type="199" charge="-0.11624" shct="0.72"  /> 
+  <GeneralizedKirkwood type="200" charge="0.05672" shct="0.85"  /> 
+  <GeneralizedKirkwood type="201" charge="0.11242" shct="0.72"  /> 
+  <GeneralizedKirkwood type="202" charge="0.02390" shct="0.85"  /> 
+  <GeneralizedKirkwood type="203" charge="-0.30941" shct="0.79"  /> 
+  <GeneralizedKirkwood type="204" charge="0.08213" shct="0.85"  /> 
+  <GeneralizedKirkwood type="205" charge="-0.05810" shct="0.72"  /> 
+  <GeneralizedKirkwood type="206" charge="0.07332" shct="0.85"  /> 
+  <GeneralizedKirkwood type="207" charge="-0.12460" shct="0.72"  /> 
+  <GeneralizedKirkwood type="208" charge="0.10013" shct="0.85"  /> 
+  <GeneralizedKirkwood type="209" charge="-0.00504" shct="0.72"  /> 
+  <GeneralizedKirkwood type="210" charge="0.06386" shct="0.85"  /> 
+  <GeneralizedKirkwood type="211" charge="-0.26789" shct="0.79"  /> 
+  <GeneralizedKirkwood type="212" charge="0.14895" shct="0.85"  /> 
+  <GeneralizedKirkwood type="213" charge="0.90680" shct="0.72"  /> 
+  <GeneralizedKirkwood type="214" charge="-0.30169" shct="0.79"  /> 
+  <GeneralizedKirkwood type="215" charge="0.14969" shct="0.85"  /> 
+  <GeneralizedKirkwood type="216" charge="0.01417" shct="0.72"  /> 
+  <GeneralizedKirkwood type="217" charge="0.06168" shct="0.85"  /> 
+  <GeneralizedKirkwood type="218" charge="-0.12127" shct="0.72"  /> 
+  <GeneralizedKirkwood type="219" charge="0.08337" shct="0.85"  /> 
+  <GeneralizedKirkwood type="220" charge="-0.00637" shct="0.72"  /> 
+  <GeneralizedKirkwood type="221" charge="0.10424" shct="0.85"  /> 
+  <GeneralizedKirkwood type="222" charge="0.10679" shct="0.79"  /> 
+  <GeneralizedKirkwood type="223" charge="0.19274" shct="0.85"  /> 
+  <GeneralizedKirkwood type="224" charge="-0.19495" shct="0.72"  /> 
+  <GeneralizedKirkwood type="225" charge="0.07797" shct="0.85"  /> 
+  <GeneralizedKirkwood type="226" charge="0.70656" shct="0.72"  /> 
+  <GeneralizedKirkwood type="227" charge="-0.74615" shct="0.85"  /> 
+  <GeneralizedKirkwood type="228" charge="-0.25887" shct="0.79"  /> 
+  <GeneralizedKirkwood type="229" charge="0.12975" shct="0.85"  /> 
+  <GeneralizedKirkwood type="230" charge="-0.27516" shct="0.79"  /> 
+  <GeneralizedKirkwood type="231" charge="0.12524" shct="0.85"  /> 
+  <GeneralizedKirkwood type="232" charge="-0.00902" shct="0.72"  /> 
+  <GeneralizedKirkwood type="233" charge="0.05319" shct="0.85"  /> 
+  <GeneralizedKirkwood type="234" charge="0.11164" shct="0.79"  /> 
+  <GeneralizedKirkwood type="234" charge="0.11164" shct="0.79"  /> 
+  <GeneralizedKirkwood type="234" charge="0.11164" shct="0.79"  /> 
+  <GeneralizedKirkwood type="234" charge="0.11164" shct="0.79"  /> 
+  <GeneralizedKirkwood type="235" charge="0.21240" shct="0.85"  /> 
+  <GeneralizedKirkwood type="235" charge="0.21240" shct="0.85"  /> 
+  <GeneralizedKirkwood type="235" charge="0.21240" shct="0.85"  /> 
+  <GeneralizedKirkwood type="235" charge="0.21240" shct="0.85"  /> 
+  <GeneralizedKirkwood type="236" charge="-0.29203" shct="0.79"  /> 
+  <GeneralizedKirkwood type="236" charge="-0.29203" shct="0.79"  /> 
+  <GeneralizedKirkwood type="236" charge="-0.29203" shct="0.79"  /> 
+  <GeneralizedKirkwood type="236" charge="-0.29203" shct="0.79"  /> 
+  <GeneralizedKirkwood type="237" charge="0.09729" shct="0.85"  /> 
+  <GeneralizedKirkwood type="238" charge="1.02669" shct="0.72"  /> 
+  <GeneralizedKirkwood type="239" charge="-0.90443" shct="0.85"  /> 
+  <GeneralizedKirkwood type="240" charge="0.90327" shct="0.72"  /> 
+  <GeneralizedKirkwood type="240" charge="0.90327" shct="0.72"  /> 
+  <GeneralizedKirkwood type="240" charge="0.90327" shct="0.72"  /> 
+  <GeneralizedKirkwood type="240" charge="0.90327" shct="0.72"  /> 
+  <GeneralizedKirkwood type="241" charge="-0.65062" shct="0.85"  /> 
+  <GeneralizedKirkwood type="241" charge="-0.65062" shct="0.85"  /> 
+  <GeneralizedKirkwood type="241" charge="-0.65062" shct="0.85"  /> 
+  <GeneralizedKirkwood type="241" charge="-0.65062" shct="0.85"  /> 
+  <GeneralizedKirkwood type="242" charge="-0.47859" shct="0.85"  /> 
+  <GeneralizedKirkwood type="243" charge="0.24567" shct="0.85"  /> 
+  <GeneralizedKirkwood type="244" charge="0.07219" shct="0.79"  /> 
+  <GeneralizedKirkwood type="245" charge="0.23262" shct="0.85"  /> 
+  <GeneralizedKirkwood type="246" charge="-0.16428" shct="0.72"  /> 
+  <GeneralizedKirkwood type="247" charge="0.93900" shct="0.72"  /> 
+  <GeneralizedKirkwood type="248" charge="-0.74794" shct="0.85"  /> 
+  <GeneralizedKirkwood type="249" charge="0.12611" shct="0.85"  /> 
+  <GeneralizedKirkwood type="250" charge="-0.04567" shct="0.72"  /> 
+  <GeneralizedKirkwood type="251" charge="0.10559" shct="0.85"  /> 
+  <GeneralizedKirkwood type="252" charge="-0.12218" shct="0.79"  /> 
+  <GeneralizedKirkwood type="253" charge="0.10618" shct="0.72"  /> 
+  <GeneralizedKirkwood type="254" charge="0.03294" shct="0.72"  /> 
+  <GeneralizedKirkwood type="255" charge="-0.29610" shct="0.79"  /> 
+  <GeneralizedKirkwood type="256" charge="0.12091" shct="0.72"  /> 
+  <GeneralizedKirkwood type="257" charge="-0.28556" shct="0.79"  /> 
+  <GeneralizedKirkwood type="258" charge="0.09440" shct="0.72"  /> 
+  <GeneralizedKirkwood type="259" charge="-0.30504" shct="0.79"  /> 
+  <GeneralizedKirkwood type="260" charge="0.11701" shct="0.72"  /> 
+  <GeneralizedKirkwood type="261" charge="0.10987" shct="0.85"  /> 
+  <GeneralizedKirkwood type="262" charge="-0.28195" shct="0.79"  /> 
+  <GeneralizedKirkwood type="263" charge="0.23408" shct="0.85"  /> 
+  <GeneralizedKirkwood type="264" charge="0.23408" shct="0.85"  /> 
+  <GeneralizedKirkwood type="265" charge="0.12346" shct="0.85"  /> 
+  <GeneralizedKirkwood type="266" charge="-0.12896" shct="0.79"  /> 
+  <GeneralizedKirkwood type="267" charge="0.25555" shct="0.72"  /> 
+  <GeneralizedKirkwood type="268" charge="-0.29612" shct="0.79"  /> 
+  <GeneralizedKirkwood type="269" charge="0.10206" shct="0.72"  /> 
+  <GeneralizedKirkwood type="270" charge="-0.08445" shct="0.72"  /> 
+  <GeneralizedKirkwood type="271" charge="0.09191" shct="0.72"  /> 
+  <GeneralizedKirkwood type="272" charge="-0.42963" shct="0.85"  /> 
+  <GeneralizedKirkwood type="273" charge="-0.28291" shct="0.79"  /> 
+  <GeneralizedKirkwood type="274" charge="0.22172" shct="0.85"  /> 
+  <GeneralizedKirkwood type="275" charge="0.22172" shct="0.85"  /> 
+  <GeneralizedKirkwood type="276" charge="0.09969" shct="0.85"  /> 
+  <GeneralizedKirkwood type="277" charge="0.11152" shct="0.85"  /> 
+  <GeneralizedKirkwood type="278" charge="-0.12123" shct="0.79"  /> 
+  <GeneralizedKirkwood type="279" charge="0.11988" shct="0.72"  /> 
+  <GeneralizedKirkwood type="280" charge="0.05939" shct="0.72"  /> 
+  <GeneralizedKirkwood type="281" charge="-0.25427" shct="0.79"  /> 
+  <GeneralizedKirkwood type="282" charge="0.05900" shct="0.72"  /> 
+  <GeneralizedKirkwood type="283" charge="-0.32183" shct="0.79"  /> 
+  <GeneralizedKirkwood type="284" charge="0.16650" shct="0.72"  /> 
+  <GeneralizedKirkwood type="285" charge="-0.20681" shct="0.79"  /> 
+  <GeneralizedKirkwood type="286" charge="0.22596" shct="0.72"  /> 
+  <GeneralizedKirkwood type="287" charge="0.25827" shct="0.85"  /> 
+  <GeneralizedKirkwood type="288" charge="-0.30237" shct="0.79"  /> 
+  <GeneralizedKirkwood type="289" charge="0.22815" shct="0.85"  /> 
+  <GeneralizedKirkwood type="290" charge="0.22815" shct="0.85"  /> 
+  <GeneralizedKirkwood type="291" charge="-0.37776" shct="0.85"  /> 
+  <GeneralizedKirkwood type="292" charge="0.12107" shct="0.85"  /> 
+  <GeneralizedKirkwood type="293" charge="-0.15031" shct="0.79"  /> 
+  <GeneralizedKirkwood type="294" charge="0.25814" shct="0.72"  /> 
+  <GeneralizedKirkwood type="295" charge="-0.18494" shct="0.79"  /> 
+  <GeneralizedKirkwood type="296" charge="0.20468" shct="0.72"  /> 
+  <GeneralizedKirkwood type="297" charge="-0.10776" shct="0.72"  /> 
+  <GeneralizedKirkwood type="298" charge="0.02014" shct="0.72"  /> 
+  <GeneralizedKirkwood type="299" charge="-0.42656" shct="0.85"  /> 
+  <GeneralizedKirkwood type="300" charge="0.27344" shct="0.85"  /> 
+  <GeneralizedKirkwood type="301" charge="-0.40427" shct="0.85"  /> 
+  <GeneralizedKirkwood type="302" charge="-0.09928" shct="0.72"  /> 
+  <GeneralizedKirkwood type="303" charge="0.12786" shct="0.85"  /> 
+  <GeneralizedKirkwood type="304" charge="0.11524" shct="0.85"  /> 
+  <GeneralizedKirkwood type="305" charge="-0.11970" shct="0.79"  /> 
+  <GeneralizedKirkwood type="306" charge="0.26026" shct="0.72"  /> 
+  <GeneralizedKirkwood type="307" charge="-0.18804" shct="0.79"  /> 
+  <GeneralizedKirkwood type="308" charge="0.20160" shct="0.72"  /> 
+  <GeneralizedKirkwood type="309" charge="-0.05028" shct="0.72"  /> 
+  <GeneralizedKirkwood type="310" charge="0.07568" shct="0.72"  /> 
+  <GeneralizedKirkwood type="311" charge="-0.42057" shct="0.85"  /> 
+  <GeneralizedKirkwood type="312" charge="0.27376" shct="0.85"  /> 
+  <GeneralizedKirkwood type="313" charge="-0.38720" shct="0.85"  /> 
+  <GeneralizedKirkwood type="314" charge="0.12060" shct="0.85"  /> 
+  <GeneralizedKirkwood type="315" charge="0.11599" shct="0.85"  /> 
+  <GeneralizedKirkwood type="316" charge="-0.52467" shct="0.85"  /> 
+  <GeneralizedKirkwood type="317" charge="0.05049" shct="0.72"  /> 
+  <GeneralizedKirkwood type="318" charge="0.07556" shct="0.85"  /> 
+  <GeneralizedKirkwood type="319" charge="0.07556" shct="0.85"  /> 
+  <GeneralizedKirkwood type="320" charge="0.10423" shct="0.72"  /> 
+  <GeneralizedKirkwood type="321" charge="0.13028" shct="0.85"  /> 
+  <GeneralizedKirkwood type="322" charge="-0.40197" shct="0.85"  /> 
+  <GeneralizedKirkwood type="323" charge="0.06316" shct="0.72"  /> 
+  <GeneralizedKirkwood type="324" charge="0.18443" shct="0.85"  /> 
+  <GeneralizedKirkwood type="325" charge="0.10143" shct="0.72"  /> 
+  <GeneralizedKirkwood type="326" charge="0.11227" shct="0.85"  /> 
+  <GeneralizedKirkwood type="327" charge="0.06213" shct="0.72"  /> 
+  <GeneralizedKirkwood type="328" charge="0.11711" shct="0.85"  /> 
+  <GeneralizedKirkwood type="329" charge="-0.47891" shct="0.85"  /> 
+  <GeneralizedKirkwood type="330" charge="0.32469" shct="0.85"  /> 
+  <GeneralizedKirkwood type="331" charge="-0.53467" shct="0.85"  /> 
+  <GeneralizedKirkwood type="332" charge="-0.52467" shct="0.85"  /> 
+  <GeneralizedKirkwood type="333" charge="0.05049" shct="0.72"  /> 
+  <GeneralizedKirkwood type="334" charge="0.07556" shct="0.85"  /> 
+  <GeneralizedKirkwood type="335" charge="0.07556" shct="0.85"  /> 
+  <GeneralizedKirkwood type="336" charge="0.10423" shct="0.72"  /> 
+  <GeneralizedKirkwood type="337" charge="0.13028" shct="0.85"  /> 
+  <GeneralizedKirkwood type="338" charge="-0.40197" shct="0.85"  /> 
+  <GeneralizedKirkwood type="339" charge="0.06316" shct="0.72"  /> 
+  <GeneralizedKirkwood type="340" charge="0.18443" shct="0.85"  /> 
+  <GeneralizedKirkwood type="341" charge="0.10143" shct="0.72"  /> 
+  <GeneralizedKirkwood type="342" charge="0.11227" shct="0.85"  /> 
+  <GeneralizedKirkwood type="343" charge="0.06213" shct="0.72"  /> 
+  <GeneralizedKirkwood type="344" charge="0.11711" shct="0.85"  /> 
+  <GeneralizedKirkwood type="345" charge="-0.47891" shct="0.85"  /> 
+  <GeneralizedKirkwood type="346" charge="0.32469" shct="0.85"  /> 
+  <GeneralizedKirkwood type="347" charge="-0.53467" shct="0.85"  /> 
+  <GeneralizedKirkwood type="348" charge="-0.52467" shct="0.85"  /> 
+  <GeneralizedKirkwood type="349" charge="0.00960" shct="0.72"  /> 
+  <GeneralizedKirkwood type="350" charge="0.07556" shct="0.85"  /> 
+  <GeneralizedKirkwood type="351" charge="0.07556" shct="0.85"  /> 
+  <GeneralizedKirkwood type="352" charge="0.15826" shct="0.72"  /> 
+  <GeneralizedKirkwood type="353" charge="0.12800" shct="0.85"  /> 
+  <GeneralizedKirkwood type="354" charge="-0.38842" shct="0.85"  /> 
+  <GeneralizedKirkwood type="355" charge="0.06720" shct="0.72"  /> 
+  <GeneralizedKirkwood type="356" charge="0.12287" shct="0.85"  /> 
+  <GeneralizedKirkwood type="357" charge="0.07060" shct="0.72"  /> 
+  <GeneralizedKirkwood type="358" charge="0.10297" shct="0.85"  /> 
+  <GeneralizedKirkwood type="359" charge="-0.10232" shct="0.72"  /> 
+  <GeneralizedKirkwood type="360" charge="0.10225" shct="0.85"  /> 
+  <GeneralizedKirkwood type="361" charge="0.09833" shct="0.85"  /> 
+  <GeneralizedKirkwood type="362" charge="-0.53467" shct="0.85"  /> 
+  <GeneralizedKirkwood type="363" charge="-0.52467" shct="0.85"  /> 
+  <GeneralizedKirkwood type="364" charge="0.00960" shct="0.72"  /> 
+  <GeneralizedKirkwood type="365" charge="0.07556" shct="0.85"  /> 
+  <GeneralizedKirkwood type="366" charge="0.07556" shct="0.85"  /> 
+  <GeneralizedKirkwood type="367" charge="0.15826" shct="0.72"  /> 
+  <GeneralizedKirkwood type="368" charge="0.12800" shct="0.85"  /> 
+  <GeneralizedKirkwood type="369" charge="-0.38842" shct="0.85"  /> 
+  <GeneralizedKirkwood type="370" charge="0.06720" shct="0.72"  /> 
+  <GeneralizedKirkwood type="371" charge="0.12287" shct="0.85"  /> 
+  <GeneralizedKirkwood type="372" charge="0.07060" shct="0.72"  /> 
+  <GeneralizedKirkwood type="373" charge="0.10297" shct="0.85"  /> 
+  <GeneralizedKirkwood type="374" charge="-0.10232" shct="0.72"  /> 
+  <GeneralizedKirkwood type="375" charge="0.10225" shct="0.85"  /> 
+  <GeneralizedKirkwood type="376" charge="0.09833" shct="0.85"  /> 
+  <GeneralizedKirkwood type="377" charge="-0.53467" shct="0.85"  /> 
+  <GeneralizedKirkwood type="378" charge="0.90070" shct="0.86"  /> 
+  <GeneralizedKirkwood type="379" charge="-0.62196" shct="0.85"  /> 
+  <GeneralizedKirkwood type="380" charge="-0.53491" shct="0.85"  /> 
+  <GeneralizedKirkwood type="380" charge="-0.53491" shct="0.85"  /> 
+  <GeneralizedKirkwood type="381" charge="0.34636" shct="0.85"  /> 
+  <GeneralizedKirkwood type="381" charge="0.34636" shct="0.85"  /> 
+  <GeneralizedKirkwood type="381" charge="0.34636" shct="0.85"  /> 
+  <GeneralizedKirkwood type="382" charge="-0.58166" shct="0.85"  /> 
+  <GeneralizedKirkwood type="382" charge="-0.58166" shct="0.85"  /> 
+  <GeneralizedKirkwood type="383" charge="1.02369" shct="0.86"  /> 
+  <GeneralizedKirkwood type="384" charge="-0.87686" shct="0.85"  /> 
+  <GeneralizedKirkwood type="385" charge="-0.50819" shct="0.85"  /> 
+  <GeneralizedKirkwood type="385" charge="-0.50819" shct="0.85"  /> 
+  <GeneralizedKirkwood type="386" charge="0.29418" shct="0.85"  /> 
+  <GeneralizedKirkwood type="386" charge="0.29418" shct="0.85"  /> 
+  <GeneralizedKirkwood type="387" charge="-0.60712" shct="0.85"  /> 
+  <GeneralizedKirkwood type="387" charge="-0.60712" shct="0.85"  /> 
+  <GeneralizedKirkwood type="388" charge="1.02369" shct="0.86"  /> 
+  <GeneralizedKirkwood type="389" charge="-0.87686" shct="0.85"  /> 
+  <GeneralizedKirkwood type="390" charge="0.90070" shct="0.86"  /> 
+  <GeneralizedKirkwood type="391" charge="-0.62196" shct="0.85"  /> 
+  <GeneralizedKirkwood type="392" charge="-0.53491" shct="0.85"  /> 
+  <GeneralizedKirkwood type="392" charge="-0.53491" shct="0.85"  /> 
+  <GeneralizedKirkwood type="393" charge="0.34636" shct="0.85"  /> 
+  <GeneralizedKirkwood type="393" charge="0.34636" shct="0.85"  /> 
+  <GeneralizedKirkwood type="394" charge="-0.58166" shct="0.85"  /> 
+  <GeneralizedKirkwood type="394" charge="-0.58166" shct="0.85"  /> 
+  <GeneralizedKirkwood type="395" charge="1.02369" shct="0.86"  /> 
+  <GeneralizedKirkwood type="396" charge="-0.87686" shct="0.85"  /> 
+  <GeneralizedKirkwood type="397" charge="-0.50819" shct="0.85"  /> 
+  <GeneralizedKirkwood type="397" charge="-0.50819" shct="0.85"  /> 
+  <GeneralizedKirkwood type="398" charge="0.29418" shct="0.85"  /> 
+  <GeneralizedKirkwood type="398" charge="0.29418" shct="0.85"  /> 
+  <GeneralizedKirkwood type="399" charge="-0.60712" shct="0.85"  /> 
+  <GeneralizedKirkwood type="399" charge="-0.60712" shct="0.85"  /> 
+  <GeneralizedKirkwood type="400" charge="1.02369" shct="0.86"  /> 
+  <GeneralizedKirkwood type="401" charge="-0.87686" shct="0.85"  /> 
+  <GeneralizedKirkwood type="402" charge="-0.51966" shct="0.85"  /> 
+  <GeneralizedKirkwood type="403" charge="0.25983" shct="0.85"  /> 
+  <GeneralizedKirkwood type="404" charge="1.00000" shct="0.8"  /> 
+  <GeneralizedKirkwood type="405" charge="1.00000" shct="0.8"  /> 
+  <GeneralizedKirkwood type="406" charge="1.00000" shct="0.8"  /> 
+  <GeneralizedKirkwood type="407" charge="1.00000" shct="0.8"  /> 
+  <GeneralizedKirkwood type="408" charge="1.00000" shct="0.8"  /> 
+  <GeneralizedKirkwood type="409" charge="2.00000" shct="0.8"  /> 
+  <GeneralizedKirkwood type="410" charge="2.00000" shct="0.8"  /> 
+  <GeneralizedKirkwood type="411" charge="2.00000" shct="0.8"  /> 
+  <GeneralizedKirkwood type="412" charge="2.00000" shct="0.8"  /> 
+  <GeneralizedKirkwood type="413" charge="-1.00000" shct="0.8"  /> 
  </AmoebaGeneralizedKirkwoodForce>
-  <AmoebaWcaDispersionForce epso="0.46024" epsh="0.056484" rmino="0.17025" rminh="0.13275" awater="33.428" slevy="1.0"  dispoff="0.026" shctd="0.81" >
-   <WcaDispersion class="1" radius="0.1855" epsilon="0.46024" /> 
-   <WcaDispersion class="2" radius="0.191" epsilon="0.422584" /> 
-   <WcaDispersion class="3" radius="0.191" epsilon="0.443504" /> 
-   <WcaDispersion class="4" radius="0.1295" epsilon="0.092048" /> 
-   <WcaDispersion class="5" radius="0.165" epsilon="0.468608" /> 
-   <WcaDispersion class="6" radius="0.147" epsilon="0.108784" /> 
-   <WcaDispersion class="7" radius="0.1825" epsilon="0.422584" /> 
-   <WcaDispersion class="8" radius="0.191" epsilon="0.422584" /> 
-   <WcaDispersion class="9" radius="0.149" epsilon="0.100416" /> 
-   <WcaDispersion class="10" radius="0.17025" epsilon="0.46024" /> 
-   <WcaDispersion class="11" radius="0.13275" epsilon="0.056484" /> 
-   <WcaDispersion class="12" radius="0.1955" epsilon="1.61084" /> 
-   <WcaDispersion class="13" radius="0.15" epsilon="0.110876" /> 
-   <WcaDispersion class="15" radius="0.1855" epsilon="0.46024" /> 
-   <WcaDispersion class="16" radius="0.191" epsilon="0.422584" /> 
-   <WcaDispersion class="17" radius="0.19" epsilon="0.372376" /> 
-   <WcaDispersion class="18" radius="0.149" epsilon="0.108784" /> 
-   <WcaDispersion class="19" radius="0.17025" epsilon="0.46024" /> 
-   <WcaDispersion class="20" radius="0.13275" epsilon="0.056484" /> 
-   <WcaDispersion class="22" radius="0.19" epsilon="0.422584" /> 
-   <WcaDispersion class="23" radius="0.149" epsilon="0.108784" /> 
-   <WcaDispersion class="24" radius="0.1855" epsilon="0.46024" /> 
-   <WcaDispersion class="25" radius="0.1295" epsilon="0.092048" /> 
-   <WcaDispersion class="26" radius="0.19" epsilon="0.422584" /> 
-   <WcaDispersion class="27" radius="0.149" epsilon="0.108784" /> 
-   <WcaDispersion class="28" radius="0.19" epsilon="0.422584" /> 
-   <WcaDispersion class="29" radius="0.1855" epsilon="0.46024" /> 
-   <WcaDispersion class="30" radius="0.191" epsilon="0.443504" /> 
-   <WcaDispersion class="31" radius="0.1725" epsilon="0.468608" /> 
-   <WcaDispersion class="32" radius="0.165" epsilon="0.468608" /> 
-   <WcaDispersion class="33" radius="0.17025" epsilon="0.46024" /> 
-   <WcaDispersion class="34" radius="0.13275" epsilon="0.06276" /> 
-   <WcaDispersion class="35" radius="0.1855" epsilon="0.43932" /> 
-   <WcaDispersion class="36" radius="0.124" epsilon="0.075312" /> 
-   <WcaDispersion class="38" radius="0.121" epsilon="0.08368" /> 
-   <WcaDispersion class="39" radius="0.1825" epsilon="0.422584" /> 
-   <WcaDispersion class="40" radius="0.191" epsilon="0.422584" /> 
-   <WcaDispersion class="41" radius="0.1855" epsilon="0.43932" /> 
-   <WcaDispersion class="42" radius="0.124" epsilon="0.075312" /> 
-   <WcaDispersion class="43" radius="0.184" epsilon="0.43932" /> 
-   <WcaDispersion class="44" radius="0.189" epsilon="0.422584" /> 
-   <WcaDispersion class="45" radius="0.175" epsilon="0.422584" /> 
-   <WcaDispersion class="46" radius="0.189" epsilon="0.422584" /> 
-   <WcaDispersion class="47" radius="0.175" epsilon="0.43932" /> 
-   <WcaDispersion class="48" radius="0.189" epsilon="0.422584" /> 
-   <WcaDispersion class="49" radius="0.189" epsilon="0.422584" /> 
-   <WcaDispersion class="50" radius="0.13" epsilon="0.075312" /> 
-   <WcaDispersion class="51" radius="0.175" epsilon="0.422584" /> 
-   <WcaDispersion class="52" radius="0.1275" epsilon="0.079496" /> 
-   <WcaDispersion class="53" radius="0.189" epsilon="0.422584" /> 
-   <WcaDispersion class="54" radius="0.189" epsilon="0.422584" /> 
-   <WcaDispersion class="55" radius="0.165" epsilon="0.468608" /> 
-   <WcaDispersion class="56" radius="0.14" epsilon="0.100416" /> 
-   <WcaDispersion class="57" radius="0.14" epsilon="0.100416" /> 
-   <WcaDispersion class="58" radius="0.175" epsilon="0.43932" /> 
-   <WcaDispersion class="59" radius="0.115" epsilon="0.0422584" /> 
-   <WcaDispersion class="60" radius="0.191" epsilon="0.422584" /> 
-   <WcaDispersion class="61" radius="0.149" epsilon="0.100416" /> 
-   <WcaDispersion class="62" radius="0.17025" epsilon="0.468608" /> 
-   <WcaDispersion class="63" radius="0.191" epsilon="0.422584" /> 
-   <WcaDispersion class="64" radius="0.149" epsilon="0.100416" /> 
-   <WcaDispersion class="65" radius="0.149" epsilon="0.100416" /> 
-   <WcaDispersion class="66" radius="0.17025" epsilon="0.46024" /> 
-   <WcaDispersion class="67" radius="0.13275" epsilon="0.056484" /> 
-   <WcaDispersion class="68" radius="0.17025" epsilon="0.46024" /> 
-   <WcaDispersion class="69" radius="0.191" epsilon="0.422584" /> 
-   <WcaDispersion class="70" radius="0.2225" epsilon="1.63176" /> 
-   <WcaDispersion class="71" radius="0.168" epsilon="0.468608" /> 
-   <WcaDispersion class="72" radius="0.13275" epsilon="0.056484" /> 
-   <WcaDispersion class="73" radius="0.17025" epsilon="0.46024" /> 
-   <WcaDispersion class="74" radius="0.13275" epsilon="0.056484" /> 
-   <WcaDispersion class="75" radius="0.119" epsilon="0.33472" /> 
-   <WcaDispersion class="76" radius="0.151" epsilon="1.08784" /> 
-   <WcaDispersion class="77" radius="0.1855" epsilon="1.4644" /> 
-   <WcaDispersion class="78" radius="0.207" epsilon="1.84096" /> 
-   <WcaDispersion class="79" radius="0.2185" epsilon="2.21752" /> 
-   <WcaDispersion class="80" radius="0.094" epsilon="0.380744" /> 
-   <WcaDispersion class="81" radius="0.1275" epsilon="3.5564" /> 
-   <WcaDispersion class="82" radius="0.1565" epsilon="6.6944" /> 
-   <WcaDispersion class="83" radius="0.2065" epsilon="1.42256" /> 
-   <WcaDispersion class="14" radius="0.2175" epsilon="1.61084" /> 
-   <WcaDispersion class="21" radius="0.166" epsilon="0.468608" /> 
-   <WcaDispersion class="37" radius="0.135" epsilon="0.08368" /> 
+ <AmoebaWcaDispersionForce epso="0.46024" epsh="0.056484" rmino="0.17025" rminh="0.13275" awater="33.428" slevy="1.0"  dispoff="0.026" shctd="0.81" >
+  <WcaDispersion class="1" radius="0.1855" epsilon="0.46024" /> 
+  <WcaDispersion class="2" radius="0.191" epsilon="0.422584" /> 
+  <WcaDispersion class="3" radius="0.191" epsilon="0.443504" /> 
+  <WcaDispersion class="4" radius="0.1295" epsilon="0.092048" /> 
+  <WcaDispersion class="5" radius="0.165" epsilon="0.468608" /> 
+  <WcaDispersion class="6" radius="0.147" epsilon="0.108784" /> 
+  <WcaDispersion class="7" radius="0.1825" epsilon="0.422584" /> 
+  <WcaDispersion class="8" radius="0.191" epsilon="0.422584" /> 
+  <WcaDispersion class="9" radius="0.149" epsilon="0.100416" /> 
+  <WcaDispersion class="10" radius="0.17025" epsilon="0.46024" /> 
+  <WcaDispersion class="11" radius="0.13275" epsilon="0.056484" /> 
+  <WcaDispersion class="12" radius="0.1955" epsilon="1.61084" /> 
+  <WcaDispersion class="13" radius="0.15" epsilon="0.110876" /> 
+  <WcaDispersion class="14" radius="0.2175" epsilon="1.61084" /> 
+  <WcaDispersion class="15" radius="0.1855" epsilon="0.46024" /> 
+  <WcaDispersion class="16" radius="0.191" epsilon="0.422584" /> 
+  <WcaDispersion class="17" radius="0.19" epsilon="0.372376" /> 
+  <WcaDispersion class="18" radius="0.149" epsilon="0.108784" /> 
+  <WcaDispersion class="19" radius="0.17025" epsilon="0.46024" /> 
+  <WcaDispersion class="20" radius="0.13275" epsilon="0.056484" /> 
+  <WcaDispersion class="21" radius="0.166" epsilon="0.468608" /> 
+  <WcaDispersion class="22" radius="0.19" epsilon="0.422584" /> 
+  <WcaDispersion class="23" radius="0.149" epsilon="0.108784" /> 
+  <WcaDispersion class="24" radius="0.1855" epsilon="0.46024" /> 
+  <WcaDispersion class="25" radius="0.1295" epsilon="0.092048" /> 
+  <WcaDispersion class="26" radius="0.19" epsilon="0.422584" /> 
+  <WcaDispersion class="27" radius="0.149" epsilon="0.108784" /> 
+  <WcaDispersion class="28" radius="0.19" epsilon="0.422584" /> 
+  <WcaDispersion class="29" radius="0.1855" epsilon="0.46024" /> 
+  <WcaDispersion class="30" radius="0.191" epsilon="0.443504" /> 
+  <WcaDispersion class="31" radius="0.1725" epsilon="0.468608" /> 
+  <WcaDispersion class="32" radius="0.165" epsilon="0.468608" /> 
+  <WcaDispersion class="33" radius="0.17025" epsilon="0.46024" /> 
+  <WcaDispersion class="34" radius="0.13275" epsilon="0.06276" /> 
+  <WcaDispersion class="35" radius="0.1855" epsilon="0.43932" /> 
+  <WcaDispersion class="36" radius="0.124" epsilon="0.075312" /> 
+  <WcaDispersion class="37" radius="0.135" epsilon="0.08368" /> 
+  <WcaDispersion class="38" radius="0.121" epsilon="0.08368" /> 
+  <WcaDispersion class="39" radius="0.1825" epsilon="0.422584" /> 
+  <WcaDispersion class="40" radius="0.191" epsilon="0.422584" /> 
+  <WcaDispersion class="41" radius="0.1855" epsilon="0.43932" /> 
+  <WcaDispersion class="42" radius="0.124" epsilon="0.075312" /> 
+  <WcaDispersion class="43" radius="0.184" epsilon="0.43932" /> 
+  <WcaDispersion class="44" radius="0.189" epsilon="0.422584" /> 
+  <WcaDispersion class="45" radius="0.175" epsilon="0.422584" /> 
+  <WcaDispersion class="46" radius="0.189" epsilon="0.422584" /> 
+  <WcaDispersion class="47" radius="0.175" epsilon="0.43932" /> 
+  <WcaDispersion class="48" radius="0.189" epsilon="0.422584" /> 
+  <WcaDispersion class="49" radius="0.189" epsilon="0.422584" /> 
+  <WcaDispersion class="50" radius="0.13" epsilon="0.075312" /> 
+  <WcaDispersion class="51" radius="0.175" epsilon="0.422584" /> 
+  <WcaDispersion class="52" radius="0.1275" epsilon="0.079496" /> 
+  <WcaDispersion class="53" radius="0.189" epsilon="0.422584" /> 
+  <WcaDispersion class="54" radius="0.189" epsilon="0.422584" /> 
+  <WcaDispersion class="55" radius="0.165" epsilon="0.468608" /> 
+  <WcaDispersion class="56" radius="0.14" epsilon="0.100416" /> 
+  <WcaDispersion class="57" radius="0.14" epsilon="0.100416" /> 
+  <WcaDispersion class="58" radius="0.175" epsilon="0.43932" /> 
+  <WcaDispersion class="59" radius="0.115" epsilon="0.0422584" /> 
+  <WcaDispersion class="60" radius="0.191" epsilon="0.422584" /> 
+  <WcaDispersion class="61" radius="0.149" epsilon="0.100416" /> 
+  <WcaDispersion class="62" radius="0.17025" epsilon="0.468608" /> 
+  <WcaDispersion class="63" radius="0.191" epsilon="0.422584" /> 
+  <WcaDispersion class="64" radius="0.149" epsilon="0.100416" /> 
+  <WcaDispersion class="65" radius="0.149" epsilon="0.100416" /> 
+  <WcaDispersion class="66" radius="0.17025" epsilon="0.46024" /> 
+  <WcaDispersion class="67" radius="0.13275" epsilon="0.056484" /> 
+  <WcaDispersion class="68" radius="0.17025" epsilon="0.46024" /> 
+  <WcaDispersion class="69" radius="0.191" epsilon="0.422584" /> 
+  <WcaDispersion class="70" radius="0.2225" epsilon="1.63176" /> 
+  <WcaDispersion class="71" radius="0.168" epsilon="0.468608" /> 
+  <WcaDispersion class="72" radius="0.13275" epsilon="0.056484" /> 
+  <WcaDispersion class="73" radius="0.17025" epsilon="0.46024" /> 
+  <WcaDispersion class="74" radius="0.13275" epsilon="0.056484" /> 
+  <WcaDispersion class="75" radius="0.119" epsilon="0.33472" /> 
+  <WcaDispersion class="76" radius="0.151" epsilon="1.08784" /> 
+  <WcaDispersion class="77" radius="0.1855" epsilon="1.4644" /> 
+  <WcaDispersion class="78" radius="0.207" epsilon="1.84096" /> 
+  <WcaDispersion class="79" radius="0.2185" epsilon="2.21752" /> 
+  <WcaDispersion class="80" radius="0.094" epsilon="0.380744" /> 
+  <WcaDispersion class="81" radius="0.147" epsilon="1.2552" /> 
+  <WcaDispersion class="82" radius="0.1815" epsilon="1.4644" /> 
+  <WcaDispersion class="83" radius="0.134" epsilon="0.928848" /> 
+  <WcaDispersion class="84" radius="0.2065" epsilon="1.42256" /> 
  </AmoebaWcaDispersionForce>
 </ForceField>

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

@@ -220,6 +220,7 @@
   <Variant name="HID"/>
   <Variant name="HIE"/>
   <Variant name="HIP"/>
+  <Variant name="HIN"/>
   <H name="H" parent="N"/>
   <H name="H2" parent="N" terminal="N"/>
   <H name="H3" parent="N" terminal="N"/>

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

@@ -28,8 +28,8 @@ 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.
 """
-from __future__ import absolute_import
-from __future__ import print_function
+from __future__ import absolute_import, print_function
+
 __author__ = "Peter Eastman"
 __version__ = "1.0"
 
@@ -119,11 +119,12 @@ class ForceField(object):
         self._atomClasses = {'':set()}
         self._forces = []
         self._scripts = []
+        self._templateGenerators = []
         for file in files:
             self.loadFile(file)
 
     def loadFile(self, file):
-        """Load an XML file and add the definitions from it to this FieldField.
+        """Load an XML file and add the definitions from it to this ForceField.
 
         Parameters
         ----------
@@ -139,6 +140,18 @@ class ForceField(object):
             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()
 
         # Load the atom types.
@@ -163,21 +176,20 @@ class ForceField(object):
                     if atomName in atomIndices:
                         raise ValueError('Residue '+resName+' contains multiple atoms named '+atomName)
                     atomIndices[atomName] = len(template.atoms)
-                    template.atoms.append(ForceField._TemplateAtomData(atomName, atom.attrib['type'], self._atomTypes[atom.attrib['type']][2], params))
+                    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.addBond(atomIndices[bond.attrib['atomName1']], atomIndices[bond.attrib['atomName2']])
+                        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:
-                        b = atomIndices[bond.attrib['atomName']]
+                        template.addExternalBondByName(bond.attrib['atomName'])
                     else:
-                        b = int(bond.attrib['from'])
-                    template.externalBonds.append(b)
-                    template.atoms[b].externalBonds += 1
+                        template.addExternalBond(int(bond.attrib['from']))
                 self.registerResidueTemplate(template)
 
         # Load force definitions
@@ -211,7 +223,7 @@ class ForceField(object):
             element = parameters['element']
             if not isinstance(element, elem.Element):
                 element = elem.get_by_symbol(element)
-        self._atomTypes[name] = (atomClass, mass, element)
+        self._atomTypes[name] = ForceField._AtomType(name, atomClass, mass, element)
         if atomClass in self._atomClasses:
             typeSet = self._atomClasses[atomClass]
         else:
@@ -233,8 +245,66 @@ class ForceField(object):
         """Register a new script to be executed after building the System."""
         self._scripts.append(script)
 
+    def registerTemplateGenerator(self, generator):
+        """Register a residue template generator that can be used to parameterize residues that do not match existing forcefield templates.
+
+        This functionality can be used to add handlers to parameterize small molecules or unnatural/modified residues.
+
+        .. CAUTION:: This method is experimental, and its API is subject to change.
+
+        Parameters
+        ----------
+        generator : function
+            A function that will be called when a residue is encountered that does not match an existing forcefield template.
+
+        When a residue without a template is encountered, the `generator` function is called with:
+
+        ::
+           success = generator(forcefield, residue)
+        ```
+
+        where `forcefield` is the calling `ForceField` object and `residue` is a simtk.openmm.app.topology.Residue object.
+
+        `generator` must conform to the following API:
+        ::
+          Parameters
+           ----------
+           forcefield : simtk.openmm.app.ForceField
+               The ForceField object to which residue templates and/or parameters are to be added.
+           residue : simtk.openmm.app.Topology.Residue
+               The residue topology for which a template is to be generated.
+
+           Returns
+           -------
+           success : bool
+               If the generator is able to successfully parameterize the residue, `True` is returned.
+               If the generator cannot parameterize the residue, it should return `False` and not modify `forcefield`.
+
+           The generator should either register a residue template directly with `forcefield.registerResidueTemplate(template)`
+           or it should call `forcefield.loadFile(file)` to load residue definitions from an ffxml file.
+
+           It can also use the `ForceField` programmatic API to add additional atom types (via `forcefield.registerAtomType(parameters)`)
+           or additional parameters.
+
+        """
+        self._templateGenerators.append(generator)
+
     def _findAtomTypes(self, attrib, num):
-        """Parse the attributes on an XML tag to find the set of atom types for each atom it involves."""
+        """Parse the attributes on an XML tag to find the set of atom types for each atom it involves.
+
+        Parameters
+        ----------
+        attrib : dict of attributes
+            The dictionary of attributes for an XML parameter tag.
+        num : int
+            The number of atom specifiers (e.g. 'class1' through 'class4') to extract.
+
+        Returns
+        -------
+        types : list
+            A list of atom types that match.
+
+        """
         types = []
         for i in range(num):
             if num == 1:
@@ -250,11 +320,15 @@ class ForceField(object):
                     types.append(None) # Unknown atom class
                 else:
                     types.append(self._atomClasses[attrib[classAttrib]])
-            else:
-                if typeAttrib not in attrib or attrib[typeAttrib] not in self._atomTypes:
+            elif typeAttrib in attrib:
+                if attrib[typeAttrib] == '':
+                    types.append(self._atomClasses[''])
+                elif attrib[typeAttrib] not in self._atomTypes:
                     types.append(None) # Unknown atom type
                 else:
                     types.append([attrib[typeAttrib]])
+            else:
+                types.append(None) # Unknown atom type
         return types
 
     def _parseTorsion(self, attrib):
@@ -306,11 +380,39 @@ class ForceField(object):
             self.bonds = []
             self.externalBonds = []
 
+        def getAtomIndexByName(self, atom_name):
+            """Look up an atom index by atom name, providing a helpful error message if not found."""
+            for (index, atom) in enumerate(self.atoms):
+                if atom.name == atom_name:
+                    return index
+
+            # Provide a helpful error message if atom name not found.
+            msg =  "Atom name '%s' not found in residue template '%s'.\n" % (atom_name, self.name)
+            msg += "Possible atom names are: %s" % str(atomIndices.keys())
+            raise ValueError(msg)
+
         def addBond(self, atom1, atom2):
+            """Add a bond between two atoms in a template given their indices in the template."""
             self.bonds.append((atom1, atom2))
             self.atoms[atom1].bondedTo.append(atom2)
             self.atoms[atom2].bondedTo.append(atom1)
 
+        def addBondByName(self, atom1_name, atom2_name):
+            """Add a bond between two atoms in a template given their atom names."""
+            atom1 = self.getAtomIndexByName(atom1_name)
+            atom2 = self.getAtomIndexByName(atom2_name)
+            self.addBond(atom1, atom2)
+
+        def addExternalBond(self, atom_index):
+            """Designate that an atom in a residue template has an external bond, using atom index within template."""
+            self.externalBonds.append(atom_index)
+            self.atoms[atom_index].externalBonds += 1
+
+        def addExternalBondByName(self, atom_name):
+            """Designate that an atom in a residue template has an external bond, using atom name within template."""
+            atom = self.getAtomIndexByName(atom_name)
+            self.addExternalBond(atom)
+
     class _TemplateAtomData:
         """Inner class used to encapsulate data about an atom in a residue template definition."""
         def __init__(self, name, type, element, parameters={}):
@@ -362,6 +464,14 @@ class ForceField(object):
             else:
                 self.excludeWith = self.atoms[0]
 
+    class _AtomType:
+        """Inner class used to record atom types and associated properties."""
+        def __init__(self, name, atomClass, mass, element):
+            self.name = name
+            self.atomClass = atomClass
+            self.mass = mass
+            self.element = element
+
     class _AtomTypeParameters:
         """Inner class used to record parameter values for atom types."""
         def __init__(self, forcefield, forceName, atomTag, paramNames):
@@ -433,6 +543,165 @@ class ForceField(object):
                 raise ValueError('%s: No parameters defined for atom type %s' % (self.forceName, t))
 
 
+    def _getResidueTemplateMatches(self, res, bondedToAtom):
+        """Return the residue template matches, or None if none are found.
+
+        Parameters
+        ----------
+        res : Topology.Residue
+            The residue for which template matches are to be retrieved.
+        bondedToAtom : list of set of int
+            bondedToAtom[i] is the set of atoms bonded to atom index i
+
+        Returns
+        -------
+        template : _ForceFieldTemplate
+            The matching forcefield residue template, or None if no matches are found.
+        matches : list
+            a list specifying which atom of the template each atom of the residue
+            corresponds to, or None if it does not match the template
+
+        """
+        template = None
+        matches = None
+        signature = _createResidueSignature([atom.element for atom in res.atoms()])
+        if signature in self._templateSignatures:
+            for t in self._templateSignatures[signature]:
+                matches = _matchResidue(res, t, bondedToAtom)
+                if matches is not None:
+                    template = t
+                    break
+        return [template, matches]
+
+    def _buildBondedToAtomList(self, topology):
+        """Build a list of which atom indices are bonded to each atom.
+
+        Parameters
+        ----------
+        topology : Topology
+            The Topology whose bonds are to be indexed.
+
+        Returns
+        -------
+        bondedToAtom : list of set of int
+            bondedToAtom[index] is the set of atom indices bonded to atom `index`
+
+        """
+        bondedToAtom = []
+        for atom in topology.atoms():
+            bondedToAtom.append(set())
+        for (atom1, atom2) in topology.bonds():
+            bondedToAtom[atom1.index].add(atom2.index)
+            bondedToAtom[atom2.index].add(atom1.index)
+        return bondedToAtom
+
+    def getUnmatchedResidues(self, topology):
+        """Return a list of Residue objects from specified topology for which no forcefield templates are available.
+
+        .. CAUTION:: This method is experimental, and its API is subject to change.
+
+        Parameters
+        ----------
+        topology : Topology
+            The Topology whose residues are to be checked against the forcefield residue templates.
+
+        Returns
+        -------
+        unmatched_residues : list of Residue
+            List of Residue objects from `topology` for which no forcefield residue templates are available.
+            Note that multiple instances of the same residue appearing at different points in the topology may be returned.
+
+        This method may be of use in generating missing residue templates or diagnosing parameterization failures.
+        """
+        # Find the template matching each residue, compiling a list of residues for which no templates are available.
+        bondedToAtom = self._buildBondedToAtomList(topology)
+        unmatched_residues = list() # list of unmatched residues
+        for res in topology.residues():
+            # Attempt to match one of the existing templates.
+            [template, matches] = self._getResidueTemplateMatches(res, bondedToAtom)
+            if matches is None:
+                # No existing templates match.
+                unmatched_residues.append(res)
+
+        return unmatched_residues
+
+    def getMatchingTemplates(self, topology):
+        """Return a list of forcefield residue templates matching residues in the specified topology.
+
+        .. CAUTION:: This method is experimental, and its API is subject to change.
+
+        Parameters
+        ----------
+        topology : Topology
+            The Topology whose residues are to be checked against the forcefield residue templates.
+
+        Returns
+        -------
+        templates : list of _TemplateData
+            List of forcefield residue templates corresponding to residues in the topology.
+            templates[index] is template corresponding to residue `index` in topology.residues()
+
+        This method may be of use in debugging issues related to parameter assignment.
+        """
+        # Find the template matching each residue, compiling a list of residues for which no templates are available.
+        bondedToAtom = self._buildBondedToAtomList(topology)
+        templates = list() # list of templates matching the corresponding residues
+        for res in topology.residues():
+            # Attempt to match one of the existing templates.
+            [template, matches] = self._getResidueTemplateMatches(res, bondedToAtom)
+            # Raise an exception if we have found no templates that match.
+            if matches is None:
+                raise ValueError('No template found for residue %d (%s).  %s' % (res.index+1, res.name, _findMatchErrors(self, res)))
+            else:
+                templates.append(template)
+
+        return templates
+
+    def generateTemplatesForUnmatchedResidues(self, topology):
+        """Generate forcefield residue templates for residues in specified topology for which no forcefield templates are available.
+
+        .. CAUTION:: This method is experimental, and its API is subject to change.
+
+        Parameters
+        ----------
+        topology : Topology
+            The Topology whose residues are to be checked against the forcefield residue templates.
+
+        Returns
+        -------
+        templates : list of _TemplateData
+            List of forcefield residue templates corresponding to residues in `topology` for which no forcefield templates are currently available.
+            Atom types will be set to `None`, but template name, atom names, elements, and connectivity will be taken from corresponding Residue objects.
+        residues : list of Residue
+            List of Residue objects that were used to generate the templates.
+            `residues[index]` is the Residue that was used to generate the template `templates[index]`
+
+        """
+        # Get a non-unique list of unmatched residues.
+        unmatched_residues = self.getUnmatchedResidues(topology)
+        # Generate a unique list of unmatched residues by comparing fingerprints.
+        bondedToAtom = self._buildBondedToAtomList(topology)
+        unique_unmatched_residues = list() # list of unique unmatched Residue objects from topology
+        templates = list() # corresponding _TemplateData templates
+        signatures = set()
+        for residue in unmatched_residues:
+            signature = _createResidueSignature([ atom.element for atom in residue.atoms() ])
+            template = _createResidueTemplate(residue)
+            is_unique = True
+            if signature in signatures:
+                # Signature is the same as an existing residue; check connectivity.
+                for check_residue in unique_unmatched_residues:
+                    matches = _matchResidue(check_residue, template, bondedToAtom)
+                    if matches is not None:
+                        is_unique = False
+            if is_unique:
+                # Residue is unique.
+                unique_unmatched_residues.append(residue)
+                signatures.add(signature)
+                templates.append(template)
+
+        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):
         """Construct an OpenMM System representing a Topology with this force field.
@@ -492,20 +761,30 @@ class ForceField(object):
             data.atomBonds[bond.atom2].append(i)
 
         # Find the template matching each residue and assign atom types.
+        # If no templates are found, attempt to use residue template generators to create new templates (and potentially atom types/parameters).
 
         for chain in topology.chains():
             for res in chain.residues():
-                template = None
-                matches = None
-                signature = _createResidueSignature([atom.element for atom in res.atoms()])
-                if signature in self._templateSignatures:
-                    for t in self._templateSignatures[signature]:
-                        matches = _matchResidue(res, t, bondedToAtom)
-                        if matches is not None:
-                            template = t
-                            break
+                # 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:
+                        if generator(self, res):
+                            # This generator has registered a new residue template that should match.
+                            [template, matches] = self._getResidueTemplateMatches(res, bondedToAtom)
+                            if matches is None:
+                                # Something went wrong because the generated template does not match the residue signature.
+                                raise Exception('The residue handler %s indicated it had correctly parameterized residue %s, but the generated template did not match the residue signature.' % (generator.__class__.__name__, str(res)))
+                            else:
+                                # We successfully generated a residue template.  Break out of the for loop.
+                                break
+
+                # Raise an exception if we have found no templates that match.
                 if matches is None:
                     raise ValueError('No template found for residue %d (%s).  %s' % (res.index+1, res.name, _findMatchErrors(self, res)))
+
+                # Store parameters for the matched residue template.
                 matchAtoms = dict(zip(matches, res.atoms()))
                 for atom, match in zip(res.atoms(), matches):
                     data.atomType[atom] = template.atoms[match].type
@@ -518,11 +797,26 @@ class ForceField(object):
 
         sys = mm.System()
         for atom in topology.atoms():
-            sys.addParticle(self._atomTypes[data.atomType[atom]][1])
+            # Look up the atom type name, returning a helpful error message if it cannot be found.
+            if atom not in data.atomType:
+                raise Exception("Could not identify atom type for atom '%s'." % str(atom))
+            typename = data.atomType[atom]
+
+            # Look up the type name in the list of registered atom types, returning a helpful error message if it cannot be found.
+            if typename not in self._atomTypes:
+                msg  = "Could not find typename '%s' for atom '%s' in list of known atom types.\n" % (typename, str(atom))
+                msg += "Known atom types are: %s" % str(self._atomTypes.keys())
+                raise Exception(msg)
 
-        # Adjust masses.
+            # Add the particle to the OpenMM system.
+            mass = self._atomTypes[typename].mass
+            sys.addParticle(mass)
+
+        # Adjust hydrogen masses if requested.
 
         if hydrogenMass is not None:
+            if not unit.is_quantity(hydrogenMass):
+                hydrogenMass *= unit.dalton
             for atom1, atom2 in topology.bonds():
                 if atom1.element == elem.hydrogen:
                     (atom1, atom2) = (atom2, atom1)
@@ -536,7 +830,7 @@ class ForceField(object):
         boxVectors = topology.getPeriodicBoxVectors()
         if boxVectors is not None:
             sys.setDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2])
-        elif nonbondedMethod is not NoCutoff and nonbondedMethod is not CutoffNonPeriodic:
+        elif nonbondedMethod not in [NoCutoff, CutoffNonPeriodic]:
             raise ValueError('Requested periodic boundary conditions for a Topology that does not specify periodic box dimensions')
 
         # Make a list of all unique angles
@@ -650,7 +944,7 @@ class ForceField(object):
         if removeCMMotion:
             sys.addForce(mm.CMMotionRemover())
 
-        # Let generators do postprocessing
+        # Let force generators do postprocessing
 
         for force in self._forces:
             if 'postprocessSystem' in dir(force):
@@ -690,7 +984,6 @@ def _createResidueSignature(elements):
         s += element.symbol+str(count)
     return s
 
-
 def _matchResidue(res, template, bondedToAtom):
     """Determine whether a residue matches a template and return a list of corresponding atoms.
 
@@ -778,7 +1071,7 @@ def _findAtomMatches(atoms, template, bondedTo, externalBonds, matches, hasMatch
 def _findMatchErrors(forcefield, res):
     """Try to guess why a residue failed to match any template and return an error message."""
     residueCounts = _countResidueAtoms([atom.element for atom in res.atoms()])
-    numResidueAtoms = sum(residueCounts.itervalues())
+    numResidueAtoms = sum(residueCounts.values())
     numResidueHeavyAtoms = sum(residueCounts[element] for element in residueCounts if element not in (None, elem.hydrogen))
 
     # Loop over templates and see how closely each one might match.
@@ -797,7 +1090,7 @@ def _findMatchErrors(forcefield, res):
 
         # If there are too many missing atoms, discard this template.
 
-        numTemplateAtoms = sum(templateCounts.itervalues())
+        numTemplateAtoms = sum(templateCounts.values())
         numTemplateHeavyAtoms = sum(templateCounts[element] for element in templateCounts if element not in (None, elem.hydrogen))
         if numTemplateAtoms > numBestMatchAtoms:
             continue
@@ -835,6 +1128,34 @@ def _findMatchErrors(forcefield, res):
         return 'The set of atoms is similar to %s, but it is missing %d atoms.' % (bestMatchName, numBestMatchAtoms-numResidueAtoms)
     return 'This might mean your input topology is missing some atoms or bonds, or possibly that you are using the wrong force field.'
 
+def _createResidueTemplate(residue):
+    """Create a _TemplateData template from a Residue object.
+
+    Parameters
+    ----------
+    residue : Residue
+        The Residue from which the template is to be constructed.
+
+    Returns
+    -------
+    template : _TemplateData
+        The residue template, with atom types set to None.
+
+    This method may be useful in creating new residue templates for residues without templates defined by the ForceField.
+
+    """
+    template = ForceField._TemplateData(residue.name)
+    for atom in residue.atoms():
+        template.atoms.append(ForceField._TemplateAtomData(atom.name, None, atom.element))
+    for (atom1,atom2) in residue.internal_bonds():
+        template.addBondByName(atom1.name, atom2.name)
+    residue_atoms = [ atom for atom in residue.atoms() ]
+    for (atom1,atom2) in residue.external_bonds():
+        if atom1 in residue_atoms:
+            template.addExternalBondByName(atom1.name)
+        elif atom2 in residue_atoms:
+            template.addExternalBondByName(atom2.name)
+    return template
 
 # The following classes are generators that know how to create Force subclasses and add them to a System that is being
 # created.  Each generator class must define two methods: 1) a static method that takes an etree Element and a ForceField,
@@ -1041,14 +1362,16 @@ class PeriodicTorsionGenerator:
             type2 = data.atomType[data.atoms[torsion[1]]]
             type3 = data.atomType[data.atoms[torsion[2]]]
             type4 = data.atomType[data.atoms[torsion[3]]]
-            done = False
+            match = None
             for tordef in self.improper:
-                if done:
-                    break
                 types1 = tordef.types1
                 types2 = tordef.types2
                 types3 = tordef.types3
                 types4 = tordef.types4
+                hasWildcard = (wildcard in (types1, types2, types3, types4))
+                if match is not None and hasWildcard:
+                    # Prefer specific definitions over ones with wildcards
+                    continue
                 if type1 in types1:
                     for (t2, t3, t4) in itertools.permutations(((type2, 1), (type3, 2), (type4, 3))):
                         if t2[0] in types2 and t3[0] in types3 and t4[0] in types4:
@@ -1063,11 +1386,13 @@ class PeriodicTorsionGenerator:
                                 (a1, a2) = (a2, a1)
                             elif e1 != elem.carbon and (e2 == elem.carbon or e1.mass < e2.mass):
                                 (a1, a2) = (a2, a1)
-                            for i in range(len(tordef.phase)):
-                                if tordef.k[i] != 0:
-                                    force.addTorsion(a1, a2, torsion[0], torsion[t4[1]], tordef.periodicity[i], tordef.phase[i], tordef.k[i])
-                            done = True
+                            match = (a1, a2, torsion[0], torsion[t4[1]], tordef)
                             break
+            if match is not None:
+                (a1, a2, a3, a4, tordef) = match
+                for i in range(len(tordef.phase)):
+                    if tordef.k[i] != 0:
+                        force.addTorsion(a1, a2, a3, a4, tordef.periodicity[i], tordef.phase[i], tordef.k[i])
 
 parsers["PeriodicTorsionForce"] = PeriodicTorsionGenerator.parseElement
 
@@ -1138,18 +1463,20 @@ class RBTorsionGenerator:
             type2 = data.atomType[data.atoms[torsion[1]]]
             type3 = data.atomType[data.atoms[torsion[2]]]
             type4 = data.atomType[data.atoms[torsion[3]]]
-            done = False
+            match = None
             for tordef in self.improper:
-                if done:
-                    break
                 types1 = tordef.types1
                 types2 = tordef.types2
                 types3 = tordef.types3
                 types4 = tordef.types4
+                hasWildcard = (wildcard in (types1, types2, types3, types4))
+                if match is not None and hasWildcard:
+                    # Prefer specific definitions over ones with wildcards
+                    continue
                 if type1 in types1:
                     for (t2, t3, t4) in itertools.permutations(((type2, 1), (type3, 2), (type4, 3))):
                         if t2[0] in types2 and t3[0] in types3 and t4[0] in types4:
-                            if wildcard in (types1, types2, types3, types4):
+                            if hasWildcard:
                                 # Workaround to be more consistent with AMBER.  It uses wildcards to define most of its
                                 # impropers, which leaves the ordering ambiguous.  It then follows some bizarre rules
                                 # to pick the order.
@@ -1161,12 +1488,14 @@ class RBTorsionGenerator:
                                     (a1, a2) = (a2, a1)
                                 elif e1 != elem.carbon and (e2 == elem.carbon or e1.mass < e2.mass):
                                     (a1, a2) = (a2, a1)
-                                force.addTorsion(a1, a2, torsion[0], torsion[t4[1]], tordef.c[0], tordef.c[1], tordef.c[2], tordef.c[3], tordef.c[4], tordef.c[5])
+                                match = (a1, a2, torsion[0], torsion[t4[1]], tordef)
                             else:
                                 # There are no wildcards, so the order is unambiguous.
-                                force.addTorsion(torsion[0], torsion[t2[1]], torsion[t3[1]], torsion[t4[1]], tordef.c[0], tordef.c[1], tordef.c[2], tordef.c[3], tordef.c[4], tordef.c[5])
-                            done = True
+                                match = (torsion[0], torsion[t2[1]], torsion[t3[1]], torsion[t4[1]], tordef)
                             break
+            if match is not None:
+                (a1, a2, a3, a4, tordef) = match
+                force.addTorsion(a1, a2, a3, a4, tordef.c[0], tordef.c[1], tordef.c[2], tordef.c[3], tordef.c[4], tordef.c[5])
 
 parsers["RBTorsionForce"] = RBTorsionGenerator.parseElement
 
@@ -1565,18 +1894,20 @@ class CustomTorsionGenerator:
             type2 = data.atomType[data.atoms[torsion[1]]]
             type3 = data.atomType[data.atoms[torsion[2]]]
             type4 = data.atomType[data.atoms[torsion[3]]]
-            done = False
+            match = None
             for tordef in self.improper:
-                if done:
-                    break
                 types1 = tordef.types1
                 types2 = tordef.types2
                 types3 = tordef.types3
                 types4 = tordef.types4
+                hasWildcard = (wildcard in (types1, types2, types3, types4))
+                if match is not None and hasWildcard:
+                    # Prefer specific definitions over ones with wildcards
+                    continue
                 if type1 in types1:
                     for (t2, t3, t4) in itertools.permutations(((type2, 1), (type3, 2), (type4, 3))):
                         if t2[0] in types2 and t3[0] in types3 and t4[0] in types4:
-                            if wildcard in (types1, types2, types3, types4):
+                            if hasWildcard:
                                 # Workaround to be more consistent with AMBER.  It uses wildcards to define most of its
                                 # impropers, which leaves the ordering ambiguous.  It then follows some bizarre rules
                                 # to pick the order.
@@ -1588,12 +1919,14 @@ class CustomTorsionGenerator:
                                     (a1, a2) = (a2, a1)
                                 elif e1 != elem.carbon and (e2 == elem.carbon or e1.mass < e2.mass):
                                     (a1, a2) = (a2, a1)
-                                force.addTorsion(a1, a2, torsion[0], torsion[t4[1]], tordef.paramValues)
+                                match = (a1, a2, torsion[0], torsion[t4[1]], tordef)
                             else:
                                 # There are no wildcards, so the order is unambiguous.
-                                force.addTorsion(torsion[0], torsion[t2[1]], torsion[t3[1]], torsion[t4[1]], tordef.paramValues)
-                            done = True
+                                match = (torsion[0], torsion[t2[1]], torsion[t3[1]], torsion[t4[1]], tordef)
                             break
+            if match is not None:
+                (a1, a2, a3, a4, tordef) = match
+                force.addTorsion(a1, a2, a3, a4, tordef.paramValues)
 
 parsers["CustomTorsionForce"] = CustomTorsionGenerator.parseElement
 
@@ -2625,7 +2958,7 @@ class AmoebaPiTorsionGenerator:
             atom1 = bond.atom1
             atom2 = bond.atom2
 
-            if (len(data.atomBonds[atom1]) == 3 and len(data.atomBonds[atom1]) == 3):
+            if (len(data.atomBonds[atom1]) == 3 and len(data.atomBonds[atom2]) == 3):
 
                 type1 = data.atomType[data.atoms[atom1]]
                 type2 = data.atomType[data.atoms[atom2]]
@@ -2910,7 +3243,7 @@ class AmoebaTorsionTorsionGenerator:
 
                                 # match in reverse order
 
-                                if (type5 in types1 and type4 in types2 and type3 in types3 and type2 in types4 and type1 in types5):
+                                elif (type5 in types1 and type4 in types2 and type3 in types3 and type2 in types4 and type1 in types5):
                                     chiralAtomIndex = self.getChiralAtomIndex(data, sys, ib, ic, id)
                                     force.addTorsionTorsion(ie, id, ic, ib, ia, chiralAtomIndex, self.gridIndex[i])
 
@@ -3601,6 +3934,8 @@ class AmoebaMultipoleGenerator:
                 polarizationType = args['polarization']
                 if (polarizationType.lower() == 'direct'):
                     force.setPolarizationType(mm.AmoebaMultipoleForce.Direct)
+                elif (polarizationType.lower() == 'extrapolated'):
+                    force.setPolarizationType(mm.AmoebaMultipoleForce.Extrapolated)
                 else:
                     force.setPolarizationType(mm.AmoebaMultipoleForce.Mutual)
 
@@ -4393,3 +4728,5 @@ class DrudeGenerator:
                     drude.addScreenedPair(drude1, drude2, thole1+thole2)
 
 parsers["DrudeForce"] = DrudeGenerator.parseElement
+
+#=============================================================================================

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

@@ -259,7 +259,7 @@ class GromacsTopFile(object):
         """Process the [ defaults ] line."""
         fields = line.split()
         if len(fields) < 4:
-            raise ValueError('Too few fields in [ defaults ] line: '+line);
+            raise ValueError('Too few fields in [ defaults ] line: '+line)
         if fields[0] != '1':
             raise ValueError('Unsupported nonbonded type: '+fields[0])
         if fields[1] != '2':
@@ -272,7 +272,7 @@ class GromacsTopFile(object):
         """Process a line in the [ moleculetypes ] category."""
         fields = line.split()
         if len(fields) < 1:
-            raise ValueError('Too few fields in [ moleculetypes ] line: '+line);
+            raise ValueError('Too few fields in [ moleculetypes ] line: '+line)
         type = GromacsTopFile._MoleculeType()
         self._moleculeTypes[fields[0]] = type
         self._currentMoleculeType = type
@@ -281,7 +281,7 @@ class GromacsTopFile(object):
         """Process a line in the [ molecules ] category."""
         fields = line.split()
         if len(fields) < 2:
-            raise ValueError('Too few fields in [ molecules ] line: '+line);
+            raise ValueError('Too few fields in [ molecules ] line: '+line)
         self._molecules.append((fields[0], int(fields[1])))
 
     def _processAtom(self, line):
@@ -290,7 +290,7 @@ class GromacsTopFile(object):
             raise ValueError('Found [ atoms ] section before [ moleculetype ]')
         fields = line.split()
         if len(fields) < 5:
-            raise ValueError('Too few fields in [ atoms ] line: '+line);
+            raise ValueError('Too few fields in [ atoms ] line: '+line)
         self._currentMoleculeType.atoms.append(fields)
 
     def _processBond(self, line):
@@ -299,9 +299,9 @@ class GromacsTopFile(object):
             raise ValueError('Found [ bonds ] section before [ moleculetype ]')
         fields = line.split()
         if len(fields) < 3:
-            raise ValueError('Too few fields in [ bonds ] line: '+line);
+            raise ValueError('Too few fields in [ bonds ] line: '+line)
         if fields[2] != '1':
-            raise ValueError('Unsupported function type in [ bonds ] line: '+line);
+            raise ValueError('Unsupported function type in [ bonds ] line: '+line)
         self._currentMoleculeType.bonds.append(fields)
 
     def _processAngle(self, line):
@@ -310,9 +310,9 @@ class GromacsTopFile(object):
             raise ValueError('Found [ angles ] section before [ moleculetype ]')
         fields = line.split()
         if len(fields) < 4:
-            raise ValueError('Too few fields in [ angles ] line: '+line);
+            raise ValueError('Too few fields in [ angles ] line: '+line)
         if fields[3] not in ('1', '5'):
-            raise ValueError('Unsupported function type in [ angles ] line: '+line);
+            raise ValueError('Unsupported function type in [ angles ] line: '+line)
         self._currentMoleculeType.angles.append(fields)
 
     def _processDihedral(self, line):
@@ -321,9 +321,9 @@ class GromacsTopFile(object):
             raise ValueError('Found [ dihedrals ] section before [ moleculetype ]')
         fields = line.split()
         if len(fields) < 5:
-            raise ValueError('Too few fields in [ dihedrals ] line: '+line);
+            raise ValueError('Too few fields in [ dihedrals ] line: '+line)
         if fields[4] not in ('1', '2', '3', '4', '9'):
-            raise ValueError('Unsupported function type in [ dihedrals ] line: '+line);
+            raise ValueError('Unsupported function type in [ dihedrals ] line: '+line)
         self._currentMoleculeType.dihedrals.append(fields)
 
     def _processExclusion(self, line):
@@ -332,7 +332,7 @@ class GromacsTopFile(object):
             raise ValueError('Found [ exclusions ] section before [ moleculetype ]')
         fields = line.split()
         if len(fields) < 2:
-            raise ValueError('Too few fields in [ exclusions ] line: '+line);
+            raise ValueError('Too few fields in [ exclusions ] line: '+line)
         self._currentMoleculeType.exclusions.append(fields)
 
     def _processPair(self, line):
@@ -341,9 +341,9 @@ class GromacsTopFile(object):
             raise ValueError('Found [ pairs ] section before [ moleculetype ]')
         fields = line.split()
         if len(fields) < 3:
-            raise ValueError('Too few fields in [ pairs ] line: '+line);
+            raise ValueError('Too few fields in [ pairs ] line: '+line)
         if fields[2] != '1':
-            raise ValueError('Unsupported function type in [ pairs ] line: '+line);
+            raise ValueError('Unsupported function type in [ pairs ] line: '+line)
         self._currentMoleculeType.pairs.append(fields)
 
     def _processCmap(self, line):
@@ -352,14 +352,14 @@ class GromacsTopFile(object):
             raise ValueError('Found [ cmap ] section before [ moleculetype ]')
         fields = line.split()
         if len(fields) < 6:
-            raise ValueError('Too few fields in [ cmap ] line: '+line);
+            raise ValueError('Too few fields in [ cmap ] line: '+line)
         self._currentMoleculeType.cmaps.append(fields)
 
     def _processAtomType(self, line):
         """Process a line in the [ atomtypes ] category."""
         fields = line.split()
         if len(fields) < 6:
-            raise ValueError('Too few fields in [ atomtypes ] line: '+line);
+            raise ValueError('Too few fields in [ atomtypes ] line: '+line)
         if len(fields[3]) == 1:
             # Bonded type and atomic number are both missing.
             fields.insert(1, None)
@@ -377,27 +377,27 @@ class GromacsTopFile(object):
         """Process a line in the [ bondtypes ] category."""
         fields = line.split()
         if len(fields) < 5:
-            raise ValueError('Too few fields in [ bondtypes ] line: '+line);
+            raise ValueError('Too few fields in [ bondtypes ] line: '+line)
         if fields[2] != '1':
-            raise ValueError('Unsupported function type in [ bondtypes ] line: '+line);
+            raise ValueError('Unsupported function type in [ bondtypes ] line: '+line)
         self._bondTypes[tuple(fields[:2])] = fields
 
     def _processAngleType(self, line):
         """Process a line in the [ angletypes ] category."""
         fields = line.split()
         if len(fields) < 6:
-            raise ValueError('Too few fields in [ angletypes ] line: '+line);
+            raise ValueError('Too few fields in [ angletypes ] line: '+line)
         if fields[3] not in ('1', '5'):
-            raise ValueError('Unsupported function type in [ angletypes ] line: '+line);
+            raise ValueError('Unsupported function type in [ angletypes ] line: '+line)
         self._angleTypes[tuple(fields[:3])] = fields
 
     def _processDihedralType(self, line):
         """Process a line in the [ dihedraltypes ] category."""
         fields = line.split()
         if len(fields) < 7:
-            raise ValueError('Too few fields in [ dihedraltypes ] line: '+line);
+            raise ValueError('Too few fields in [ dihedraltypes ] line: '+line)
         if fields[4] not in ('1', '2', '3', '4', '9'):
-            raise ValueError('Unsupported function type in [ dihedraltypes ] line: '+line);
+            raise ValueError('Unsupported function type in [ dihedraltypes ] line: '+line)
         key = tuple(fields[:5])
         if fields[4] == '9' and key in self._dihedralTypes:
             # There are multiple dihedrals defined for these atom types.
@@ -409,25 +409,25 @@ class GromacsTopFile(object):
         """Process a line in the [ implicit_genborn_params ] category."""
         fields = line.split()
         if len(fields) < 6:
-            raise ValueError('Too few fields in [ implicit_genborn_params ] line: '+line);
+            raise ValueError('Too few fields in [ implicit_genborn_params ] line: '+line)
         self._implicitTypes[fields[0]] = fields
 
     def _processPairType(self, line):
         """Process a line in the [ pairtypes ] category."""
         fields = line.split()
         if len(fields) < 5:
-            raise ValueError('Too few fields in [ pairtypes] line: '+line);
+            raise ValueError('Too few fields in [ pairtypes] line: '+line)
         if fields[2] != '1':
-            raise ValueError('Unsupported function type in [ pairtypes ] line: '+line);
+            raise ValueError('Unsupported function type in [ pairtypes ] line: '+line)
         self._pairTypes[tuple(fields[:2])] = fields
 
     def _processCmapType(self, line):
         """Process a line in the [ cmaptypes ] category."""
         fields = line.split()
         if len(fields) < 8 or len(fields) < 8+int(fields[6])*int(fields[7]):
-            raise ValueError('Too few fields in [ cmaptypes ] line: '+line);
+            raise ValueError('Too few fields in [ cmaptypes ] line: '+line)
         if fields[5] != '1':
-            raise ValueError('Unsupported function type in [ cmaptypes ] line: '+line);
+            raise ValueError('Unsupported function type in [ cmaptypes ] line: '+line)
         self._cmapTypes[tuple(fields[:5])] = fields
 
     def __init__(self, file, periodicBoxVectors=None, unitCellDimensions=None, includeDir=None, defines=None):

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

@@ -34,8 +34,7 @@ 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.
 """
-from __future__ import absolute_import
-from __future__ import print_function
+from __future__ import absolute_import, print_function
 
 #=============================================================================================
 # GLOBAL IMPORTS
@@ -553,87 +552,6 @@ class PrmtopLoader(object):
             self._excludedAtoms.append(atomList)
         return self._excludedAtoms
 
-    def getGBParms(self, gbmodel, elements):
-        """Return list giving GB params, Radius and screening factor"""
-        gb_List=[]
-        radii=[float(r) for r in self._raw_data["RADII"]]
-        screen=[float(s) for s in self._raw_data["SCREEN"]]
-        if gbmodel == 'GBn2':
-            alpha = [0 for i in self._raw_data['RADII']]
-            beta = [0 for i in self._raw_data['RADII']]
-            gamma = [0 for i in self._raw_data['RADII']]
-        # Update screening parameters for GBn if specified
-        if gbmodel == 'GBn':
-            if elements is None:
-                raise Exception('GBn model requires element information')
-            for i, element in enumerate(elements):
-                if element is elem.carbon:
-                    screen[i] = 0.48435382330
-                elif element is elem.hydrogen:
-                    screen[i] = 1.09085413633
-                elif element is elem.nitrogen:
-                    screen[i] = 0.700147318409
-                elif element is elem.oxygen:
-                    screen[i] = 1.06557401132
-                elif element is elem.sulfur:
-                    screen[i] = 0.602256336067
-                else:
-                    screen[i] = 0.5
-                # radii is currently in Angstroms right now. GBn lookup tables
-                # only support radii between 1.0 and 2.0
-                if radii[i] < 1.0 or radii[i] > 2.0:
-                    raise ValueError('GBn requires intrinsic radii between 1 and '
-                                     '2 Angstroms (%.3f found for atom %d)' %
-                                     (radii[i], i))
-        if gbmodel == 'GBn2':
-            if elements is None:
-                raise Exception('GBn2 model requires element information')
-            for i, element in enumerate(elements):
-                if element is elem.carbon:
-                    screen[i] = 1.058554
-                    alpha[i] = 0.733756
-                    beta[i] = 0.506378
-                    gamma[i] = 0.205844
-                elif element is elem.hydrogen:
-                    screen[i] = 1.425952
-                    alpha[i] = 0.788440
-                    beta[i] = 0.798699
-                    gamma[i] = 0.437334
-                elif element is elem.nitrogen:
-                    screen[i] = 0.733599
-                    alpha[i] = 0.503364
-                    beta[i] = 0.316828
-                    gamma[i] = 0.192915
-                elif element is elem.oxygen:
-                    screen[i] = 1.061039
-                    alpha[i] = 0.867814
-                    beta[i] = 0.876635
-                    gamma[i] = 0.387882
-                elif element is elem.sulfur:
-                    screen[i] = -0.703469
-                    alpha[i] = 0.867814
-                    beta[i] = 0.876635
-                    gamma[i] = 0.387882
-                else: # not optimized
-                    screen[i] = 0.5
-                    alpha[i] = 1.0
-                    beta[i] = 0.8
-                    gamma[i] = 4.85
-                # radii is currently in Angstroms right now. GBn lookup tables
-                # only support radii between 1.0 and 2.0
-                if radii[i] < 1.0 or radii[i] > 2.0:
-                    raise ValueError('GBn2 requires intrinsic radii between 1 and '
-                                     '2 Angstroms (%.3f found for atom %d)' %
-                                     (radii[i], i))
-        lengthConversionFactor = units.angstrom.conversion_factor_to(units.nanometer)
-        if gbmodel == 'GBn2':
-            for rad, scr, alp, bet, gam in zip(radii, screen, alpha, beta, gamma):
-                gb_List.append((rad*lengthConversionFactor, scr, alp, bet, gam))
-        else:
-            for rad, scr in zip(radii, screen):
-                gb_List.append((rad*lengthConversionFactor, scr))
-        return gb_List
-
     def getBoxBetaAndDimensions(self):
         """Return periodic boundary box beta angle and dimensions"""
         beta=float(self._raw_data["BOX_DIMENSIONS"][0])
@@ -1053,9 +971,6 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
             cutoff = nonbondedCutoff
             if units.is_quantity(cutoff):
                 cutoff = cutoff.value_in_unit(units.nanometers)
-        gb_parms = prmtop.getGBParms(gbmodel, elements)
-        if gbmodel != 'OBC2' or implicitSolventKappa != 0:
-            gb_parms = customgb.convertParameters(gb_parms, gbmodel)
         if gbmodel == 'HCT':
             gb = customgb.GBSAHCTForce(solventDielectric, soluteDielectric, 'ACE', cutoff, implicitSolventKappa)
         elif gbmodel == 'OBC1':
@@ -1072,7 +987,29 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
         elif gbmodel == 'GBn2':
             gb = customgb.GBSAGBn2Force(solventDielectric, soluteDielectric, 'ACE', cutoff, implicitSolventKappa)
         else:
-            raise Exception("Illegal value specified for implicit solvent model")
+            raise ValueError("Illegal value specified for implicit solvent model")
+        if isinstance(gb, mm.GBSAOBCForce):
+            # Built-in GBSAOBCForce does not have getStandardParameters, so use
+            # the one from the equivalent CustomGBForce
+            gb_parms = customgb.GBSAOBC2Force.getStandardParameters(topology)
+        else:
+            gb_parms = type(gb).getStandardParameters(topology)
+        # Replace radii and screen, but screen *only* gets replaced by the
+        # prmtop contents for HCT, OBC1, and OBC2. GBn and GBn2 both override
+        # the prmtop screen factors from LEaP in sander and pmemd
+        if gbmodel in ('HCT', 'OBC1', 'OBC2'):
+            screen = [float(s) for s in prmtop._raw_data['SCREEN']]
+        else:
+            screen = [gb_parm[1] for gb_parm in gb_parms]
+        radii = [float(r)/10 for r in prmtop._raw_data['RADII']]
+        warned = False
+        for i, (r, s) in enumerate(zip(radii, screen)):
+            if abs(r - gb_parms[i][0]) > 1e-4 or abs(s - gb_parms[i][1]) > 1e-4:
+                if not warned:
+                    warnings.warn('Non-optimal GB parameters detected for GB '
+                                  'model %s' % gbmodel)
+                    warned = True
+            gb_parms[i][0], gb_parms[i][1] = r, s
 
         for charge, gb_parm in zip(charges, gb_parms):
             if gbmodel == 'OBC2' and implicitSolventKappa == 0:

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

@@ -29,10 +29,18 @@ OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
 USE OR OTHER DEALINGS IN THE SOFTWARE.
 """
 
-from __future__ import division
-from __future__ import absolute_import
+from __future__ import division, absolute_import
 
+from collections import defaultdict
+import copy
+from simtk.openmm.app import element as E
 from simtk.openmm import CustomGBForce, Continuous2DFunction
+import simtk.unit as u
+
+def strip_unit(value, unit):
+    if not u.is_quantity(value):
+        return value
+    return value.value_in_unit(unit)
 
 d0=[2.26685,2.32548,2.38397,2.44235,2.50057,2.55867,2.61663,2.67444,
     2.73212,2.78965,2.84705,2.9043,2.96141,3.0184,3.07524,3.13196,
@@ -189,6 +197,141 @@ for i in range (len(d0)):
     d0[i]=d0[i]/10
     m0[i]=m0[i]*10
 
+def _get_bonded_atom_list(topology):
+    """ Returns a list of atoms bonded to each other atom in a dict """
+    bondeds = defaultdict(list)
+    for a1, a2 in topology.bonds():
+        bondeds[a1].append(a2)
+        bondeds[a2].append(a1)
+    return bondeds
+
+def _is_carboxylateO(atom, all_bonds):
+    if atom is not E.oxygen: return False
+    bondeds = all_bonds[atom]
+    if len(bondeds) != 1:
+        return False
+    if bondeds[0].element is not E.carbon:
+        return False
+    bondedsC = all_bonds[bondeds[0]]
+    if len(bondedsC) != 3:
+        return False
+    for a3 in bondedsC:
+        if a3 is atom: continue
+        if a3.element is E.oxygen:
+            break
+    else:
+        return False
+    # If we got here, must be a carboxylate
+    return True
+
+def _bondi_radii(topology):
+    """ Sets the bondi radii """
+    radii = [0.0 for atom in topology.atoms()]
+    for i, atom in enumerate(topology.atoms()):
+        if atom.element is E.carbon:
+            radii[i] = 1.7
+        elif atom.element in (E.hydrogen, E.deuterium):
+            radii[i] = 1.2
+        elif atom.element is E.nitrogen:
+            radii[i] = 1.55
+        elif atom.element is E.oxygen:
+            radii[i] = 1.5
+        elif atom.element is E.fluorine:
+            radii[i] = 1.5
+        elif atom.element is E.silicon:
+            radii[i] = 2.1
+        elif atom.element is E.phosphorus:
+            radii[i] = 1.85
+        elif atom.element is E.sulfur:
+            radii[i] = 1.8
+        elif atom.element is E.chlorine:
+            radii[i] = 1.5
+        else:
+            radii[i] = 1.5
+    return radii  # converted to nanometers above
+
+def _mbondi_radii(topology):
+    """ Sets the mbondi radii """
+    radii = [0.0 for atom in topology.atoms()]
+    all_bonds = _get_bonded_atom_list(topology)
+    for i, atom in enumerate(topology.atoms()):
+        # Radius of H atom depends on element it is bonded to
+        if atom.element in (E.hydrogen, E.deuterium):
+            bondeds = all_bonds[atom]
+            if bondeds[0].element in (E.carbon, E.nitrogen):
+                radii[i] = 1.3
+            elif bondeds[0].type.atomic_number in (E.oxygen, E.sulfur):
+                radii[i] = 0.8
+            else:
+                radii[i] = 1.2
+        # Radius of C atom depends on what type it is
+        elif atom.element is E.carbon:
+            radii[i] = 1.7
+        # All other elements have fixed radii for all types/partners
+        elif atom.element is E.nitrogen:
+            radii[i] = 1.55
+        elif atom.element is E.oxygen:
+            radii[i] = 1.5
+        elif atom.element is E.fluorine:
+            radii[i] = 1.5
+        elif atom.element is E.silicon:
+            radii[i] = 2.1
+        elif atom.element is E.phosphorus:
+            radii[i] = 1.85
+        elif atom.element is E.sulfur:
+            radii[i] = 1.8
+        elif atom.element is E.chlorine:
+            radii[i] = 1.5
+        else:
+            radii[i] = 1.5
+    return radii  # converted to nanometers above
+
+def _mbondi2_radii(topology):
+    """ Sets the mbondi2 radii """
+    radii = [0.0 for atom in topology.atoms()]
+    all_bonds = _get_bonded_atom_list(topology)
+    for i, atom in enumerate(topology.atoms()):
+        # Radius of H atom depends on element it is bonded to
+        if atom.element in (E.hydrogen, E.deuterium):
+            bondeds = all_bonds[atom]
+            if bondeds[0].element is E.nitrogen:
+                radii[i] = 1.3
+            else:
+                radii[i] = 1.2
+        # Radius of C atom depends on what type it is
+        elif atom.element is E.carbon:
+            radii[i] = 1.7
+        # All other elements have fixed radii for all types/partners
+        elif atom.element is E.nitrogen:
+            radii[i] = 1.55
+        elif atom.element is E.oxygen:
+            radii[i] = 1.5
+        elif atom.element is E.fluorine:
+            radii[i] = 1.5
+        elif atom.element is E.silicon:
+            radii[i] = 2.1
+        elif atom.element == E.phosphorus:
+            radii[i] = 1.85
+        elif atom.element == E.sulfur:
+            radii[i] = 1.8
+        elif atom.element == E.chlorine:
+            radii[i] = 1.5
+        else:
+            radii[i] = 1.5
+    return radii  # Converted to nanometers above
+
+def _mbondi3_radii(topology):
+    """ Sets the mbondi3 radii """
+    radii = _mbondi2_radii(topology)
+    all_bonds = _get_bonded_atom_list(topology)
+    for i, atom in enumerate(topology.atoms()):
+        # carboxylate and HH/HE (ARG)
+        if _is_carboxylateO(atom, all_bonds):
+            radii[i] = 1.4
+        elif atom.residue.name == 'ARG':
+            if atom.name.startswith('HH') or atom.name.startswith('HE'):
+                radii[i] = 1.17
+    return radii  # Converted to nanometers above
 
 def _createEnergyTerms(force, solventDielectric, soluteDielectric, SA, cutoff, kappa, offset):
     # Add the energy terms to the CustomGBForce.  These are identical for all the GB models.
@@ -224,159 +367,239 @@ def _createEnergyTerms(force, solventDielectric, soluteDielectric, SA, cutoff, k
             force.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2*(1/f-"+str(1/cutoff)+");"
                                 "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))"+params, CustomGBForce.ParticlePairNoExclusions)
 
+_SCREEN_PARAMETERS = { # normal, GBn, GBn2
+        E.hydrogen : (0.85, 1.09085413633, 1.425952),
+        E.carbon : (0.72, 0.48435382330, 1.058554),
+        E.nitrogen : (0.79, 0.700147318409, 0.733599),
+        E.oxygen : (0.85, 1.06557401132, 1.061039),
+        E.fluorine : (0.88, 0.5, 0.5),
+        E.phosphorus : (0.86, 0.5, 0.5),
+        E.sulfur : (0.96, 0.602256336067, -0.703469),
+        None : (0.8, 0.5, 0.5) # default
+}
+_SCREEN_PARAMETERS[E.deuterium] = _SCREEN_PARAMETERS[E.hydrogen]
+def _screen_parameter(atom):
+    if atom.element in _SCREEN_PARAMETERS:
+        return _SCREEN_PARAMETERS[atom.element]
+    return _SCREEN_PARAMETERS[None]
+
+class CustomAmberGBForce(CustomGBForce):
+
+    OFFSET = 0.009
+    RADIUS_ARG_POSITION = 1
+    SCREEN_POSITION = 2
+
+    def __init__(self, *args, **kwargs):
+        raise NotImplementedError('Cannot instantiate ABC')
+
+    def addParticle(self, params):
+        params = copy.deepcopy(params)
+        params[self.RADIUS_ARG_POSITION] = strip_unit(params[self.RADIUS_ARG_POSITION], u.nanometer) - self.OFFSET
+        params[self.SCREEN_POSITION] *= params[self.RADIUS_ARG_POSITION]
+        CustomGBForce.addParticle(self, params)
+        return params
+
+    def setParticleParameters(self, idx, params):
+        params = copy.deepcopy(params)
+        params[self.RADIUS_ARG_POSITION] = strip_unit(params[self.RADIUS_ARG_POSITION], u.nanometer) - self.OFFSET
+        params[self.SCREEN_POSITION] *= params[self.RADIUS_ARG_POSITION]
+        CustomGBForce.addParticle(self, params)
+        return params
+
+    @staticmethod
+    def getStandardParameters(topology):
+        """ Gets list of standard parameters for this GB model based on an input Topology
+
+        Parameters
+        ----------
+        topology : simtk.openmm.app.Topology
+            Topology of the system to get parameters for
+
+        Returns
+        -------
+        list of float
+            List of all parameters needed for this GB model. These can be passed
+            to addParticle or setParticleParameters after the charge is inserted
+            at the beginning of the list
+        """
+        raise NotImplementedError('Must override getStandardParameters in subclasses')
+
 """
 Amber Equivalent: igb = 1
 """
+class GBSAHCTForce(CustomAmberGBForce):
 
-
-def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None,
+    def __init__(self, solventDielectric=78.5, soluteDielectric=1, SA=None,
                  cutoff=None, kappa=0.0):
-
-    custom = CustomGBForce()
-
-    custom.addPerParticleParameter("q")
-    custom.addPerParticleParameter("or") # Offset radius
-    custom.addPerParticleParameter("sr") # Scaled offset radius
-    custom.addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
-                                  "U=r+sr2;"
-                                  "L=max(or1, D);"
-                                  "D=abs(r-sr2)", CustomGBForce.ParticlePairNoExclusions)
-
-    custom.addComputedValue("B", "1/(1/or-I)", CustomGBForce.SingleParticle)
-    _createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
-    return custom
+        CustomGBForce.__init__(self)
+
+        self.addPerParticleParameter("q")
+        self.addPerParticleParameter("or") # Offset radius
+        self.addPerParticleParameter("sr") # Scaled offset radius
+        self.addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
+                                   "U=r+sr2;"
+                                   "L=max(or1, D);"
+                                   "D=abs(r-sr2)",
+                              CustomGBForce.ParticlePairNoExclusions)
+
+        self.addComputedValue("B", "1/(1/or-I)", CustomGBForce.SingleParticle)
+        _createEnergyTerms(self, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
+
+    @staticmethod
+    def getStandardParameters(topology):
+        radii = [[x/10] for x in _mbondi_radii(topology)]
+        for i, atom in enumerate(topology.atoms()):
+            radii[i].append(_screen_parameter(atom)[0])
+        return radii
 
 """
 Amber Equivalents: igb = 2
 """
-def GBSAOBC1Force(solventDielectric=78.5, soluteDielectric=1, SA=None,
-                  cutoff=None, kappa=0.0):
+class GBSAOBC1Force(CustomAmberGBForce):
+
+    def __init__(self, solventDielectric=78.5, soluteDielectric=1, SA=None,
+                 cutoff=None, kappa=0.0):
 
-    custom = CustomGBForce()
+        CustomGBForce.__init__(self)
 
-    custom.addPerParticleParameter("q")
-    custom.addPerParticleParameter("or") # Offset radius
-    custom.addPerParticleParameter("sr") # Scaled offset radius
-    custom.addComputedValue("I",  "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
-                                  "U=r+sr2;"
-                                  "L=max(or1, D);"
-                                  "D=abs(r-sr2)", CustomGBForce.ParticlePairNoExclusions)
+        self.addPerParticleParameter("q")
+        self.addPerParticleParameter("or") # Offset radius
+        self.addPerParticleParameter("sr") # Scaled offset radius
+        self.addComputedValue("I",  "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
+                                    "U=r+sr2;"
+                                    "L=max(or1, D);"
+                                    "D=abs(r-sr2)", CustomGBForce.ParticlePairNoExclusions)
 
-    custom.addComputedValue("B", "1/(1/or-tanh(0.8*psi+2.909125*psi^3)/radius);"
-                                 "psi=I*or; radius=or+offset; offset=0.009", CustomGBForce.SingleParticle)
-    _createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
-    return custom
+        self.addComputedValue("B", "1/(1/or-tanh(0.8*psi+2.909125*psi^3)/radius);"
+                                   "psi=I*or; radius=or+offset; offset=0.009", CustomGBForce.SingleParticle)
+        _createEnergyTerms(self, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
+
+    @staticmethod
+    def getStandardParameters(topology):
+        radii = [[x/10] for x in _mbondi2_radii(topology)]
+        for i, atom in enumerate(topology.atoms()):
+            radii[i].append(_screen_parameter(atom)[0])
+        return radii
 
 """
 Amber Equivalents: igb = 5
 """
-def GBSAOBC2Force(solventDielectric=78.5, soluteDielectric=1, SA=None,
-                  cutoff=None, kappa=0.0):
+class GBSAOBC2Force(GBSAOBC1Force):
+
+    def __init__(self, solventDielectric=78.5, soluteDielectric=1, SA=None,
+                 cutoff=None, kappa=0.0):
 
-    custom = CustomGBForce()
+        CustomGBForce.__init__(self)
 
-    custom.addPerParticleParameter("q")
-    custom.addPerParticleParameter("or") # Offset radius
-    custom.addPerParticleParameter("sr") # Scaled offset radius
-    custom.addComputedValue("I",  "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
-                                  "U=r+sr2;"
-                                  "L=max(or1, D);"
-                                  "D=abs(r-sr2)", CustomGBForce.ParticlePairNoExclusions)
+        self.addPerParticleParameter("q")
+        self.addPerParticleParameter("or") # Offset radius
+        self.addPerParticleParameter("sr") # Scaled offset radius
+        self.addComputedValue("I",  "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
+                                    "U=r+sr2;"
+                                    "L=max(or1, D);"
+                                    "D=abs(r-sr2)", CustomGBForce.ParticlePairNoExclusions)
 
-    custom.addComputedValue("B", "1/(1/or-tanh(psi-0.8*psi^2+4.85*psi^3)/radius);"
-                                 "psi=I*or; radius=or+offset; offset=0.009", CustomGBForce.SingleParticle)
-    _createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
-    return custom
+        self.addComputedValue("B", "1/(1/or-tanh(psi-0.8*psi^2+4.85*psi^3)/radius);"
+                                     "psi=I*or; radius=or+offset; offset=0.009", CustomGBForce.SingleParticle)
+        _createEnergyTerms(self, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
 
 """
 Amber Equivalents: igb = 7
 """
-def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None,
-                  cutoff=None, kappa=0.0):
-
+class GBSAGBnForce(CustomAmberGBForce):
 
-    """
-    Indexing for tables:
-        input: radius1, radius2
-        index = (radius2*200-20)*21 + (radius1*200-20)
-        output: index of desired value in row-by-row, 1D version of Tables 3 & 4
-    """
-
-
-    custom = CustomGBForce()
-
-    custom.addPerParticleParameter("q")
-    custom.addPerParticleParameter("or") # Offset radius
-    custom.addPerParticleParameter("sr") # Scaled offset radius
-
-    custom.addTabulatedFunction("getd0", Continuous2DFunction(21, 21, d0, 0.1, 0.2, 0.1, 0.2))
-    custom.addTabulatedFunction("getm0", Continuous2DFunction(21, 21, m0, 0.1, 0.2, 0.1, 0.2))
-
-    custom.addComputedValue("I",  "Ivdw+neckScale*Ineck;"
-                                  "Ineck=step(radius1+radius2+neckCut-r)*getm0(radius1,radius2)/(1+100*(r-getd0(radius1,radius2))^2+0.3*1000000*(r-getd0(radius1,radius2))^6);"
-                                  "Ivdw=step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
-                                  "U=r+sr2;"
-                                  "L=max(or1, D);"
-                                  "D=abs(r-sr2);"
-                                  "radius1=or1+offset; radius2=or2+offset;"
-                                  "neckScale=0.361825; neckCut=0.68; offset=0.009", CustomGBForce.ParticlePairNoExclusions)
+    def __init__(self, solventDielectric=78.5, soluteDielectric=1, SA=None,
+                 cutoff=None, kappa=0.0):
 
-    custom.addComputedValue("B", "1/(1/or-tanh(1.09511284*psi-1.907992938*psi^2+2.50798245*psi^3)/radius);"
-                              "psi=I*or; radius=or+offset; offset=0.009", CustomGBForce.SingleParticle)
-    _createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
-    return custom
+        CustomGBForce.__init__(self)
+
+        self.addPerParticleParameter("q")
+        self.addPerParticleParameter("or") # Offset radius
+        self.addPerParticleParameter("sr") # Scaled offset radius
+
+        self.addTabulatedFunction("getd0", Continuous2DFunction(21, 21, d0, 0.1, 0.2, 0.1, 0.2))
+        self.addTabulatedFunction("getm0", Continuous2DFunction(21, 21, m0, 0.1, 0.2, 0.1, 0.2))
+
+        self.addComputedValue("I",  "Ivdw+neckScale*Ineck;"
+                                    "Ineck=step(radius1+radius2+neckCut-r)*getm0(radius1,radius2)/(1+100*(r-getd0(radius1,radius2))^2+0.3*1000000*(r-getd0(radius1,radius2))^6);"
+                                    "Ivdw=step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
+                                    "U=r+sr2;"
+                                    "L=max(or1, D);"
+                                    "D=abs(r-sr2);"
+                                    "radius1=or1+offset; radius2=or2+offset;"
+                                    "neckScale=0.361825; neckCut=0.68; offset=0.009", CustomGBForce.ParticlePairNoExclusions)
+
+        self.addComputedValue("B", "1/(1/or-tanh(1.09511284*psi-1.907992938*psi^2+2.50798245*psi^3)/radius);"
+                                  "psi=I*or; radius=or+offset; offset=0.009", CustomGBForce.SingleParticle)
+        _createEnergyTerms(self, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
+
+    def addParticle(self, parameters):
+        parameters = CustomAmberGBForce.addParticle(self, parameters)
+        if parameters[1] < 0.1 or parameters[1] > 0.2:
+            raise ValueError('Radii must be between 1 and 2 Angstroms for neck lookup')
+
+    def setParticleParameters(self, idx, parameters):
+        parameters = CustomAmberGBForce.setParticleParameters(self, idx, parameters)
+        if parameters[1] < 0.1 or parameters[1] > 0.2:
+            raise ValueError('Radii must be between 1 and 2 Angstroms for neck lookup')
+
+    @staticmethod
+    def getStandardParameters(topology):
+        radii = [[x/10] for x in _bondi_radii(topology)]
+        for i, atom in enumerate(topology.atoms()):
+            radii[i].append(_screen_parameter(atom)[1])
+        return radii
 
 """
 Amber Equivalents: igb = 8
 """
-def GBSAGBn2Force(solventDielectric=78.5, soluteDielectric=1, SA=None,
-                  cutoff=None, kappa=0.0):
-
-
-    """
-    Indexing for tables:
-        input: radius1, radius2
-        index = (radius2*200-20)*21 + (radius1*200-20)
-        output: index of desired value in row-by-row, 1D version of Tables 3 & 4
-    """
+class GBSAGBn2Force(GBSAGBnForce):
 
+    OFFSET = 0.0195141
 
-    custom = CustomGBForce()
-
-    custom.addPerParticleParameter("q")
-    custom.addPerParticleParameter("or") # Offset radius
-    custom.addPerParticleParameter("sr") # Scaled offset radius
-    custom.addPerParticleParameter("alpha")
-    custom.addPerParticleParameter("beta")
-    custom.addPerParticleParameter("gamma")
-
-    custom.addTabulatedFunction("getd0", Continuous2DFunction(21, 21, d0, 0.1, 0.2, 0.1, 0.2))
-    custom.addTabulatedFunction("getm0", Continuous2DFunction(21, 21, m0, 0.1, 0.2, 0.1, 0.2))
-
-    custom.addComputedValue("I",  "Ivdw+neckScale*Ineck;"
-                                  "Ineck=step(radius1+radius2+neckCut-r)*getm0(radius1,radius2)/(1+100*(r-getd0(radius1,radius2))^2+0.3*1000000*(r-getd0(radius1,radius2))^6);"
-                                  "Ivdw=step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
-                                  "U=r+sr2;"
-                                  "L=max(or1, D);"
-                                  "D=abs(r-sr2);"
-                                  "radius1=or1+offset; radius2=or2+offset;"
-                                  "neckScale=0.826836; neckCut=0.68; offset=0.0195141", CustomGBForce.ParticlePairNoExclusions)
-
-    custom.addComputedValue("B", "1/(1/or-tanh(alpha*psi-beta*psi^2+gamma*psi^3)/radius);"
-                                 "psi=I*or; radius=or+offset; offset=0.0195141", CustomGBForce.SingleParticle)
-    _createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.0195141)
-    return custom
-
-
-def convertParameters(params, gbmodel):
-    """Convert the GB parameters from the file into the values expected by the appropriate CustomGBForce."""
-    if gbmodel == 'GBn2':
-        offset = 0.0195141
-    else:
-        offset = 0.009
+    def __init__(self, solventDielectric=78.5, soluteDielectric=1, SA=None,
+                 cutoff=None, kappa=0.0):
 
-    for p in params:
-        newParam = list(p)
-        newParam[0] -= offset
-        newParam[1] *= newParam[0]
-        yield newParam
+        CustomGBForce.__init__(self)
+
+        self.addPerParticleParameter("q")
+        self.addPerParticleParameter("or") # Offset radius
+        self.addPerParticleParameter("sr") # Scaled offset radius
+        self.addPerParticleParameter("alpha")
+        self.addPerParticleParameter("beta")
+        self.addPerParticleParameter("gamma")
+
+        self.addTabulatedFunction("getd0", Continuous2DFunction(21, 21, d0, 0.1, 0.2, 0.1, 0.2))
+        self.addTabulatedFunction("getm0", Continuous2DFunction(21, 21, m0, 0.1, 0.2, 0.1, 0.2))
+
+        self.addComputedValue("I",  "Ivdw+neckScale*Ineck;"
+                                    "Ineck=step(radius1+radius2+neckCut-r)*getm0(radius1,radius2)/(1+100*(r-getd0(radius1,radius2))^2+0.3*1000000*(r-getd0(radius1,radius2))^6);"
+                                    "Ivdw=step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
+                                    "U=r+sr2;"
+                                    "L=max(or1, D);"
+                                    "D=abs(r-sr2);"
+                                    "radius1=or1+offset; radius2=or2+offset;"
+                                    "neckScale=0.826836; neckCut=0.68; offset=0.0195141", CustomGBForce.ParticlePairNoExclusions)
+
+        self.addComputedValue("B", "1/(1/or-tanh(alpha*psi-beta*psi^2+gamma*psi^3)/radius);"
+                                     "psi=I*or; radius=or+offset; offset=0.0195141", CustomGBForce.SingleParticle)
+        _createEnergyTerms(self, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.0195141)
+
+    @staticmethod
+    def getStandardParameters(topology):
+        radii = [[x/10] for x in _mbondi3_radii(topology)]
+        for i, atom in enumerate(topology.atoms()):
+            radii[i].append(_screen_parameter(atom)[2])
+            if atom.element in (E.hydrogen, E.deuterium):
+                radii[i].extend([0.788440, 0.798699, 0.437334])
+            elif atom.element is E.carbon:
+                radii[i].extend([0.733756, 0.506378, 0.205844])
+            elif atom.element is E.nitrogen:
+                radii[i].extend([0.503364, 0.316828, 0.192915])
+            elif atom.element is E.oxygen:
+                radii[i].extend([0.867814, 0.876635, 0.387882])
+            elif atom.element is E.sulfur:
+                radii[i].extend([0.867814, 0.876635, 0.387882])
+            else:
+                radii[i].extend([0.8, 4.85, 0.5])
+        return radii

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

@@ -123,7 +123,7 @@ class PdbStructure(object):
     """
 
 
-    def __init__(self, input_stream, load_all_models=False):
+    def __init__(self, input_stream, load_all_models=False, extraParticleIdentifier='EP'):
         """Create a PDB model from a PDB file stream.
 
         Parameters
@@ -135,9 +135,12 @@ class PdbStructure(object):
         load_all_models : bool
             Whether to load every model of an NMR structure or trajectory, or
             just load the first model, to save memory.
+        extraParticleIdentifier : string='EP'
+            if this value appears in the element column for an ATOM record, the Atom's element will be set to 'EP' to mark it as an extra particle
         """
         # initialize models
         self.load_all_models = load_all_models
+        self.extraParticleIdentifier = extraParticleIdentifier
         self.models = []
         self._current_model = None
         self.default_model = None
@@ -157,7 +160,7 @@ class PdbStructure(object):
                 pdb_line = pdb_line.decode('utf-8')
             # Look for atoms
             if (pdb_line.find("ATOM  ") == 0) or (pdb_line.find("HETATM") == 0):
-                self._add_atom(Atom(pdb_line, self))
+                self._add_atom(Atom(pdb_line, self, self.extraParticleIdentifier))
             # Notice MODEL punctuation, for the next level of detail
             # in the structure->model->chain->residue->atom->position hierarchy
             elif (pdb_line.find("MODEL") == 0):
@@ -682,7 +685,7 @@ class Residue(object):
 class Atom(object):
     """Atom represents one atom in a PDB structure.
     """
-    def __init__(self, pdb_line, pdbstructure=None):
+    def __init__(self, pdb_line, pdbstructure=None, extraParticleIdentifier='EP'):
         """Create a new pdb.Atom from an ATOM or HETATM line.
 
         Example line:
@@ -795,11 +798,14 @@ class Atom(object):
         try: self.formal_charge = int(pdb_line[78:80])
         except ValueError: self.formal_charge = None
         # figure out atom element
-        try:
-            # Try to find a sensible element symbol from columns 76-77
-            self.element = element.get_by_symbol(self.element_symbol)
-        except KeyError:
-            self.element = None
+        if self.element_symbol == extraParticleIdentifier:
+            self.element = 'EP'
+        else:
+            try:
+                # Try to find a sensible element symbol from columns 76-77
+                self.element = element.get_by_symbol(self.element_symbol)
+            except KeyError:    
+                self.element = None
         if pdbstructure is not None:
             pdbstructure._next_atom_number = self.serial_number+1
             pdbstructure._next_residue_number = self.residue_number+1

wrappers/python/simtk/openmm/app/modeller.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:
 
@@ -264,7 +264,7 @@ class Modeller(object):
         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 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.
+        4. Ion pairs are added to give the requested total ionic strength.  Note that only monovalent ions are currently supported.
 
         The box size can be specified in any of several ways:
 
@@ -298,6 +298,7 @@ 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.
+            Note that only monovalent ions are currently supported.
         neutralize : bool=True
             whether to add ions to neutralize the system
         """
@@ -522,7 +523,7 @@ class Modeller(object):
         # Add ions based on the desired ionic strength.
 
         numIons = len(addedWaters)*ionicStrength/(55.4*molar) # Pure water is about 55.4 molar (depending on temperature)
-        numPairs = int(floor(numIons/2+0.5))
+        numPairs = int(floor(numIons+0.5))
         for i in range(numPairs):
             addIon(positiveElement)
         for i in range(numPairs):
@@ -614,6 +615,7 @@ class Modeller(object):
             HID: Neutral form with a hydrogen on the ND1 atom
             HIE: Neutral form with a hydrogen on the NE2 atom
             HIP: Positively charged form with hydrogens on both ND1 and NE2
+            HIN: Negatively charged form without a hydrogen on either ND1 or NE2
 
         Lysine:
             LYN: Neutral form with two hydrogens on the zeta nitrogen
@@ -625,9 +627,14 @@ class Modeller(object):
         2. Any Cysteine that participates in a disulfide bond uses the CYX variant regardless of pH.
         3. For a neutral Histidine residue, the HID or HIE variant is selected based on which one forms a better hydrogen bond.
 
-        You can override these rules by explicitly specifying a variant for any residue.  Also keep in mind that this
-        function will only add hydrogens.  It will never remove ones that are already present in the model, regardless
-        of the specified pH.
+        You can override these rules by explicitly specifying a variant for any residue.  To do that, provide a list for the
+        'variants' parameter, and set the corresponding element to the name of the variant to use.
+
+        A special case is when the model already contains a hydrogen that should not be present in the desired variant.
+        If you explicitly specify a variant using the 'variants' parameter, the residue will be modified to match the
+        desired variant, removing hydrogens if necessary.  On the other hand, for residues whose variant is selected
+        automatically, this function will only add hydrogens.  It will never remove ones that are already present in the
+        model, regardless of the specified pH.
 
         Definitions for standard amino acids and nucleotides are built in.  You can call loadHydrogenDefinitions() to load
         additional definitions for other residue types.
@@ -771,6 +778,7 @@ class Modeller(object):
                     if variant is not None and variant not in spec.variants:
                         raise ValueError('Illegal variant for %s residue: %s' % (residue.name, variant))
                     actualVariants[residue.index] = variant
+                    removeExtraHydrogens = (variants[residue.index] is not None)
 
                     # Make a list of hydrogens that should be present in the residue.
 
@@ -783,6 +791,11 @@ class Modeller(object):
                     # Loop over atoms in the residue, adding them to the new topology along with required hydrogens.
 
                     for parent in residue.atoms():
+                        # Check whether this is a hydrogen that should be removed.
+
+                        if removeExtraHydrogens and parent.element == elem.hydrogen and not any(parent.name == h.name for h in hydrogens):
+                            continue
+
                         # Add the atom.
 
                         newAtom = newTopology.addAtom(parent.name, parent.element, newResidue)
@@ -1091,7 +1104,7 @@ class Modeller(object):
                             else:
                                 a2 = newAtoms[matchingAtoms[atom2]]
                             newTopology.addBond(a1, a2)
-                            
+
         for bond in self.topology.bonds():
             if bond[0] in newAtoms and bond[1] in newAtoms:
                 newTopology.addBond(newAtoms[bond[0]], newAtoms[bond[1]])
@@ -1109,7 +1122,7 @@ class Modeller(object):
 
             # Now run a few iterations of SHAKE to try to select reasonable positions.
 
-            for iteration in range(10):
+            for iteration in range(15):
                 for atom1, atom2, distance in bonds:
                     if atom1 in missingPositions:
                         if atom2 in missingPositions:
@@ -1123,6 +1136,6 @@ class Modeller(object):
                     delta *= (distance-length)/length
                     newPositions[atom1] -= weights[0]*delta
                     newPositions[atom2] += weights[1]*delta
-                
+
         self.topology = newTopology
         self.positions = newPositions