Merge pull request #4 from peastman/master

Most of the features to support Drude particles are now complete and functional

plugins/drude/CMakeLists.txt

@@ -205,9 +205,9 @@ ELSE (EXECUTABLE_OUTPUT_PATH)
   SET (TEST_PATH .)
 ENDIF (EXECUTABLE_OUTPUT_PATH)
 
-#IF (OPENMM_BUILD_SERIALIZATION_SUPPORT)
-#    ADD_SUBDIRECTORY(serialization)
-#ENDIF (OPENMM_BUILD_SERIALIZATION_SUPPORT)
+IF (OPENMM_BUILD_SERIALIZATION_SUPPORT)
+    ADD_SUBDIRECTORY(serialization)
+ENDIF (OPENMM_BUILD_SERIALIZATION_SUPPORT)
 
 #INCLUDE(ApiDoxygen.cmake)
 

plugins/drude/platforms/cuda/src/CudaDrudeKernelFactory.cpp

 /* -------------------------------------------------------------------------- *
- *                              OpenMMAmoeba                                  *
+ *                              OpenMMDrude                                   *
  * -------------------------------------------------------------------------- *
  * This is part of the OpenMM molecular simulation toolkit originating from   *
  * Simbios, the NIH National Center for Physics-Based Simulation of           *

plugins/drude/platforms/opencl/src/OpenCLDrudeKernelFactory.cpp

 /* -------------------------------------------------------------------------- *
- *                              OpenMMAmoeba                                  *
+ *                              OpenMMDrude                                   *
  * -------------------------------------------------------------------------- *
  * This is part of the OpenMM molecular simulation toolkit originating from   *
  * Simbios, the NIH National Center for Physics-Based Simulation of           *

plugins/drude/platforms/reference/src/ReferenceDrudeKernelFactory.cpp

 /* -------------------------------------------------------------------------- *
- *                              OpenMMAmoeba                                  *
+ *                              OpenMMDrude                                   *
  * -------------------------------------------------------------------------- *
  * 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) 2011 Stanford University and the Authors.           *
+ * Portions copyright (c) 2011-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *

plugins/drude/serialization/CMakeLists.txt

+#---------------------------------------------------
+# OpenMMDrude Serialization Library
+#
+# Creates OpenMMDrude serializatin library, base name=OpenMMDrudeSerialization.
+# Default libraries are shared & optimized. Variants
+# are created for static (_static) and debug (_d).
+#
+# Windows:
+#   OpenMMDrudeSerialization[_d].dll
+#   OpenMMDrudeSerialization[_d].lib
+# Unix:
+#   libOpenMMDrudeSerialization[_d].so
+#----------------------------------------------------
+
+# The source is organized into subdirectories, but we handle them all from
+# this CMakeLists file rather than letting CMake visit them as SUBDIRS.
+SET(OPENMM_SOURCE_SUBDIRS . ../../../serialization)
+SET(OPENMM_DRUDE_SOURCE_SUBDIRS . )
+
+# Collect up information about the version of the OpenMM library we're building
+# and make it available to the code so it can be built into the binaries.
+
+SET(OPENMM_DRUDE_LIBRARY_NAME OpenMMDrude)
+SET(OPENMM_SERIALIZATION_LIBRARY_NAME OpenMMSerialization)
+SET(OPENMM_DRUDE_SERIALIZATION_LIBRARY_NAME OpenMMDrudeSerialization)
+
+# Ensure that debug libraries have "_d" appended to their names.
+# CMake gets this right on Windows automatically with this definition.
+IF (${CMAKE_GENERATOR} MATCHES "Visual Studio")
+    SET(CMAKE_DEBUG_POSTFIX "_d" CACHE INTERNAL "" FORCE)
+ENDIF (${CMAKE_GENERATOR} MATCHES "Visual Studio")
+
+# But on Unix or Cygwin we have to add the suffix manually
+IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
+    SET(OPENMM_DRUDE_LIBRARY_NAME ${OPENMM_DRUDE_LIBRARY_NAME}_d)
+    SET(OPENMM_SERIALIZATION_LIBRARY_NAME ${OPENMM_SERIALIZATION_LIBRARY_NAME}_d)
+    SET(OPENMM_DRUDE_SERIALIZATION_LIBRARY_NAME ${OPENMM_DRUDE_SERIALIZATION_LIBRARY_NAME}_d)
+ENDIF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
+
+# These are all the places to search for header files which are
+# to be part of the API.
+SET(API_INCLUDE_DIRS) # start empty
+FOREACH(subdir ${OPENMM_SOURCE_SUBDIRS})
+    # append
+    SET(API_INCLUDE_DIRS ${API_INCLUDE_DIRS}
+                         ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include
+                         ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include/internal)
+ENDFOREACH(subdir)
+
+# We'll need both *relative* path names, starting with their API_INCLUDE_DIRS,
+# and absolute pathnames.
+SET(API_REL_INCLUDE_FILES)   # start these out empty
+SET(API_ABS_INCLUDE_FILES)
+
+FOREACH(dir ${API_INCLUDE_DIRS})
+    FILE(GLOB fullpaths ${dir}/*.h)	# returns full pathnames
+    SET(API_ABS_INCLUDE_FILES ${API_ABS_INCLUDE_FILES} ${fullpaths})
+
+    FOREACH(pathname ${fullpaths})
+        GET_FILENAME_COMPONENT(filename ${pathname} NAME)
+        SET(API_REL_INCLUDE_FILES ${API_REL_INCLUDE_FILES} ${dir}/${filename})
+    ENDFOREACH(pathname)
+ENDFOREACH(dir)
+
+# collect up source files
+SET(SOURCE_FILES) # empty
+SET(SOURCE_INCLUDE_FILES)
+
+FOREACH(subdir ${OPENMM_DRUDE_SOURCE_SUBDIRS})
+    FILE(GLOB_RECURSE src_files   ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*.cpp  ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*.c)
+    FILE(GLOB incl_files  ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*.h  ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*.hpp)
+    SET(SOURCE_FILES         ${SOURCE_FILES}         ${src_files})   #append
+    SET(SOURCE_INCLUDE_FILES ${SOURCE_INCLUDE_FILES} ${incl_files})
+    INCLUDE_DIRECTORIES(BEFORE  ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include)
+    INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/../../../serialization/${subdir}/include)
+ENDFOREACH(subdir)
+#Message( "API_REL_INCLUDE_FILES=${API_REL_INCLUDE_FILES}\n" )
+
+INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/src )
+
+# Create the library
+ADD_LIBRARY(${OPENMM_DRUDE_SERIALIZATION_LIBRARY_NAME} SHARED ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES} ${API_ABS_INCLUDE_FILES})
+TARGET_LINK_LIBRARIES(${OPENMM_DRUDE_SERIALIZATION_LIBRARY_NAME} ${OPENMM_DRUDE_LIBRARY_NAME} ${OPENMM_SERIALIZATION_LIBRARY_NAME} ${SHARED_TARGET})
+SET_TARGET_PROPERTIES(${OPENMM_DRUDE_SERIALIZATION_LIBRARY_NAME} PROPERTIES COMPILE_FLAGS "-DOPENMM_DRUDE_SERIALIZATION_BUILDING_SHARED_LIBRARY -DTIXML_USE_STL -DIEEE_8087")
+INSTALL_TARGETS(/lib/plugins RUNTIME_DIRECTORY /lib/plugins ${OPENMM_DRUDE_SERIALIZATION_LIBRARY_NAME})
+
+ADD_SUBDIRECTORY(tests)

plugins/drude/serialization/include/openmm/serialization/DrudeForceProxy.h

+#ifndef OPENMM_DRUDE_FORCE_PROXY_H_
+#define OPENMM_DRUDE_FORCE_PROXY_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                OpenMMDrude                                 *
+ * -------------------------------------------------------------------------- *
+ * 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) 2013 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/internal/windowsExportDrudeSerialization.h"
+#include "openmm/serialization/SerializationProxy.h"
+
+namespace OpenMM {
+
+/**
+ * This is a proxy for serializing DrudeForce objects.
+ */
+
+class OPENMM_EXPORT_DRUDE_SERIALIZATION DrudeForceProxy : public SerializationProxy {
+public:
+    DrudeForceProxy();
+    void serialize(const void* object, SerializationNode& node) const;
+    void* deserialize(const SerializationNode& node) const;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_DRUDE_FORCE_PROXY_H_*/

plugins/drude/serialization/include/openmm/serialization/DrudeLangevinIntegratorProxy.h

+#ifndef OPENMM_DRUDE_LANGEVIN_INTEGRATOR_PROXY_H_
+#define OPENMM_DRUDE_LANGEVIN_INTEGRATOR_PROXY_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                OpenMMDrude                                 *
+ * -------------------------------------------------------------------------- *
+ * 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) 2013 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/internal/windowsExportDrudeSerialization.h"
+#include "openmm/serialization/SerializationProxy.h"
+
+namespace OpenMM {
+
+/**
+ * This is a proxy for serializing DrudeLangevinIntegrator objects.
+ */
+
+class OPENMM_EXPORT_DRUDE_SERIALIZATION DrudeLangevinIntegratorProxy : public SerializationProxy {
+public:
+    DrudeLangevinIntegratorProxy();
+    void serialize(const void* object, SerializationNode& node) const;
+    void* deserialize(const SerializationNode& node) const;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_DRUDE_LANGEVIN_INTEGRATOR_PROXY_H_*/

plugins/drude/serialization/include/openmm/serialization/internal/windowsExportDrudeSerialization.h

+#ifndef OPENMM_WINDOWSEXPORTDRUDESERIALIZATION_H_
+#define OPENMM_WINDOWSEXPORTDRUDESERIALIZATION_H_
+
+/*
+ * Shared libraries are messy in Visual Studio. We have to distinguish three
+ * cases:
+ *   (1) this header is being used to build the OpenMM shared library
+ *       (dllexport)
+ *   (2) this header is being used by a *client* of the OpenMM shared
+ *       library (dllimport)
+ *   (3) we are building the OpenMM static library, or the client is
+ *       being compiled with the expectation of linking with the
+ *       OpenMM static library (nothing special needed)
+ * In the CMake script for building this library, we define one of the symbols
+ *     OPENMM_DRUDE_SERIALIZATION_BUILDING_{SHARED|STATIC}_LIBRARY
+ * Client code normally has no special symbol defined, in which case we'll
+ * assume it wants to use the shared library. However, if the client defines
+ * the symbol OPENMM_USE_STATIC_LIBRARIES we'll suppress the dllimport so
+ * that the client code can be linked with static libraries. Note that
+ * the client symbol is not library dependent, while the library symbols
+ * affect only the OpenMM library, meaning that other libraries can
+ * be clients of this one. However, we are assuming all-static or all-shared.
+ */
+
+#ifdef _MSC_VER
+    // We don't want to hear about how sprintf is "unsafe".
+    #pragma warning(disable:4996)
+    // Keep MS VC++ quiet about lack of dll export of private members.
+    #pragma warning(disable:4251)
+    #if defined(OPENMM_DRUDE_SERIALIZATION_BUILDING_SHARED_LIBRARY)
+        #define OPENMM_EXPORT_DRUDE_SERIALIZATION __declspec(dllexport)
+    #elif defined(OPENMM_DRUDE_SERIALIZATION_BUILDING_STATIC_LIBRARY) || defined(OPENMM_DRUDE_SERIALIZATION_USE_STATIC_LIBRARIES)
+        #define OPENMM_EXPORT_DRUDE_SERIALIZATION
+    #else
+        #define OPENMM_EXPORT_DRUDE_SERIALIZATION __declspec(dllimport)   // i.e., a client of a shared library
+    #endif
+#else
+    #define OPENMM_EXPORT_DRUDE_SERIALIZATION // Linux, Mac
+#endif
+
+#endif // OPENMM_WINDOWSEXPORTDRUDESERIALIZATION_H_

plugins/drude/serialization/src/DrudeForceProxy.cpp

+/* -------------------------------------------------------------------------- *
+ *                                OpenMMDrude                                 *
+ * -------------------------------------------------------------------------- *
+ * 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) 2013 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/DrudeForceProxy.h"
+#include "openmm/serialization/SerializationNode.h"
+#include "openmm/Force.h"
+#include "openmm/DrudeForce.h"
+#include <sstream>
+
+using namespace OpenMM;
+using namespace std;
+
+DrudeForceProxy::DrudeForceProxy() : SerializationProxy("DrudeForce") {
+}
+
+void DrudeForceProxy::serialize(const void* object, SerializationNode& node) const {
+    node.setIntProperty("version", 1);
+    const DrudeForce& force = *reinterpret_cast<const DrudeForce*>(object);
+    SerializationNode& particles = node.createChildNode("Particles");
+    for (int i = 0; i < force.getNumParticles(); i++) {
+        int p, p1, p2, p3, p4;
+        double charge, polarizability, aniso12, aniso34;
+        force.getParticleParameters(i, p, p1, p2, p3, p4, charge, polarizability, aniso12, aniso34);
+        particles.createChildNode("Particle").setIntProperty("p", p).setIntProperty("p1", p1).setIntProperty("p2", p2).setIntProperty("p3", p3).setIntProperty("p4", p4)
+            .setDoubleProperty("charge", charge).setDoubleProperty("polarizability", polarizability).setDoubleProperty("a12", aniso12).setDoubleProperty("a34", aniso34);
+    }
+    SerializationNode& pairs = node.createChildNode("ScreenedPairs");
+    for (int i = 0; i < force.getNumScreenedPairs(); i++) {
+        int p1, p2;
+        double thole;
+        force.getScreenedPairParameters(i, p1, p2, thole);
+        pairs.createChildNode("Pair").setIntProperty("p1", p1).setIntProperty("p2", p2).setDoubleProperty("thole", thole);
+    }
+}
+
+void* DrudeForceProxy::deserialize(const SerializationNode& node) const {
+    if (node.getIntProperty("version") != 1)
+        throw OpenMMException("Unsupported version number");
+    DrudeForce* force = new DrudeForce();
+    try {
+        const SerializationNode& particles = node.getChildNode("Particles");
+        for (int i = 0; i < (int) particles.getChildren().size(); i++) {
+            const SerializationNode& particle = particles.getChildren()[i];
+            force->addParticle(particle.getIntProperty("p"), particle.getIntProperty("p1"), particle.getIntProperty("p2"), particle.getIntProperty("p3"), particle.getIntProperty("p4"),
+                    particle.getDoubleProperty("charge"), particle.getDoubleProperty("polarizability"), particle.getDoubleProperty("a12"), particle.getDoubleProperty("a34"));
+        }
+        const SerializationNode& pairs = node.getChildNode("ScreenedPairs");
+        for (int i = 0; i < (int) pairs.getChildren().size(); i++) {
+            const SerializationNode& pair = pairs.getChildren()[i];
+            force->addScreenedPair(pair.getIntProperty("p1"), pair.getIntProperty("p2"), pair.getDoubleProperty("thole"));
+        }
+    }
+    catch (...) {
+        delete force;
+        throw;
+    }
+    return force;
+}

plugins/drude/serialization/src/DrudeLangevinIntegratorProxy.cpp

+/* -------------------------------------------------------------------------- *
+ *                                OpenMMDrude                                 *
+ * -------------------------------------------------------------------------- *
+ * 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) 2013 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/DrudeLangevinIntegratorProxy.h"
+#include "openmm/serialization/SerializationNode.h"
+#include "openmm/DrudeLangevinIntegrator.h"
+#include <sstream>
+
+using namespace OpenMM;
+using namespace std;
+
+DrudeLangevinIntegratorProxy::DrudeLangevinIntegratorProxy() : SerializationProxy("DrudeLangevinIntegrator") {
+}
+
+void DrudeLangevinIntegratorProxy::serialize(const void* object, SerializationNode& node) const {
+    node.setIntProperty("version", 1);
+    const DrudeLangevinIntegrator& integrator = *reinterpret_cast<const DrudeLangevinIntegrator*>(object);
+    node.setDoubleProperty("stepSize", integrator.getStepSize());
+    node.setDoubleProperty("constraintTolerance", integrator.getConstraintTolerance());
+    node.setDoubleProperty("temperature", integrator.getTemperature());
+    node.setDoubleProperty("friction", integrator.getFriction());
+    node.setDoubleProperty("drudeTemperature", integrator.getDrudeTemperature());
+    node.setDoubleProperty("drudeFriction", integrator.getDrudeFriction());
+    node.setIntProperty("randomSeed", integrator.getRandomNumberSeed());
+}
+
+void* DrudeLangevinIntegratorProxy::deserialize(const SerializationNode& node) const {
+    if (node.getIntProperty("version") != 1)
+        throw OpenMMException("Unsupported version number");
+    DrudeLangevinIntegrator *integrator = new DrudeLangevinIntegrator(node.getDoubleProperty("temperature"),
+            node.getDoubleProperty("friction"), node.getDoubleProperty("drudeTemperature"),
+            node.getDoubleProperty("drudeFriction"), node.getDoubleProperty("stepSize"));
+    integrator->setConstraintTolerance(node.getDoubleProperty("constraintTolerance"));
+    integrator->setRandomNumberSeed(node.getIntProperty("randomSeed"));
+    return integrator;
+}
\ No newline at end of file

plugins/drude/serialization/src/DrudeSerializationProxyRegistration.cpp

+/* -------------------------------------------------------------------------- *
+ *                                OpenMMDrude                                 *
+ * -------------------------------------------------------------------------- *
+ * 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) 2013 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+#ifdef WIN32
+#include <windows.h>
+#include <sstream>
+#else
+#include <dlfcn.h>
+#include <dirent.h>
+#include <cstdlib>
+#endif
+
+#include "openmm/OpenMMException.h"
+
+#include "openmm/DrudeForce.h"
+#include "openmm/DrudeLangevinIntegrator.h"
+
+#include "openmm/serialization/SerializationProxy.h"
+
+#include "openmm/serialization/DrudeForceProxy.h"
+#include "openmm/serialization/DrudeLangevinIntegratorProxy.h"
+#include "openmm/serialization/internal/windowsExportDrudeSerialization.h"
+
+#if defined(WIN32)
+    #include <windows.h>
+    extern "C" OPENMM_EXPORT_DRUDE_SERIALIZATION void registerDrudeSerializationProxies();
+    BOOL WINAPI DllMain(HANDLE hModule, DWORD  ul_reason_for_call, LPVOID lpReserved) {
+        if (ul_reason_for_call == DLL_PROCESS_ATTACH)
+            registerDrudeSerializationProxies();
+        return TRUE;
+    }
+#else
+    extern "C" void __attribute__((constructor)) registerDrudeSerializationProxies();
+#endif
+
+using namespace OpenMM;
+
+extern "C" OPENMM_EXPORT_DRUDE_SERIALIZATION void registerDrudeSerializationProxies() {
+    SerializationProxy::registerProxy(typeid(DrudeForce), new DrudeForceProxy());
+    SerializationProxy::registerProxy(typeid(DrudeLangevinIntegrator), new DrudeLangevinIntegratorProxy());
+}

plugins/drude/serialization/tests/CMakeLists.txt

+#
+# Testing
+#
+
+ENABLE_TESTING()
+
+# Automatically create tests using files named "Test*.cpp"
+FILE(GLOB TEST_PROGS "*Test*.cpp")
+FOREACH(TEST_PROG ${TEST_PROGS})
+    GET_FILENAME_COMPONENT(TEST_ROOT ${TEST_PROG} NAME_WE)
+
+    # All tests use shared libraries
+    ADD_EXECUTABLE(${TEST_ROOT} ${TEST_PROG})
+    TARGET_LINK_LIBRARIES(${TEST_ROOT} ${OPENMM_SERIALIZATION_LIBRARY_NAME} ${OPENMM_DRUDE_SERIALIZATION_LIBRARY_NAME} ${OPENMM_DRUDE_LIBRARY_NAME})
+    ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT})
+
+ENDFOREACH(TEST_PROG ${TEST_PROGS})

plugins/drude/serialization/tests/TestSerializeDrudeForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                OpenMMDrude                                 *
+ * -------------------------------------------------------------------------- *
+ * 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) 2013 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/Platform.h"
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/DrudeForce.h"
+#include "openmm/serialization/XmlSerializer.h"
+#include <iostream>
+#include <sstream>
+
+using namespace OpenMM;
+using namespace std;
+
+extern "C" void registerDrudeSerializationProxies();
+
+void testSerialization() {
+    // Create a Force.
+
+    DrudeForce force1;
+    force1.addParticle(0, 1, 2, 3, 4, 0.5, 1.0, 1.5, 2.0);
+    force1.addParticle(2, 3, 7, 8, 9, 0.1, 1e-4, 1.0, 0.9);
+    force1.addParticle(5, 6, -1, -1, -1, 0.2, 0.1, 1.0, 1.0);
+    force1.addScreenedPair(0, 1, 2.6);
+    force1.addScreenedPair(1, 2, 2.2);
+
+    // Serialize and then deserialize it.
+
+    stringstream buffer;
+    XmlSerializer::serialize<DrudeForce>(&force1, "Force", buffer);
+    DrudeForce* copy = XmlSerializer::deserialize<DrudeForce>(buffer);
+
+    // Compare the two forces to see if they are identical.  
+    DrudeForce& force2 = *copy;
+    ASSERT_EQUAL(force1.getNumParticles(), force2.getNumParticles());
+    for (int i = 0; i < (int) force1.getNumParticles(); i++) {
+        int a1, a2, a3, a4, a5, b1, b2, b3, b4, b5;
+        double charge1, charge2;
+        double polar1, polar2;
+        double aa12, ba12, aa34, ba34;
+        force1.getParticleParameters(i, a1, a2, a3, a4, a5, charge1, polar1, aa12, aa34);
+        force2.getParticleParameters(i, b1, b2, b3, b4, b5, charge2, polar2, ba12, ba34);
+        ASSERT_EQUAL(a1, b1);
+        ASSERT_EQUAL(a2, b2);
+        ASSERT_EQUAL(a3, b3);
+        ASSERT_EQUAL(a4, b4);
+        ASSERT_EQUAL(a5, b5);
+        ASSERT_EQUAL(charge1, charge2);
+        ASSERT_EQUAL(polar1, polar2);
+        ASSERT_EQUAL(aa12, ba12);
+        ASSERT_EQUAL(aa34, ba34);
+    }
+    ASSERT_EQUAL(force1.getNumScreenedPairs(), force2.getNumScreenedPairs());
+    for (int i = 0; i < (int) force1.getNumScreenedPairs(); i++) {
+        int a1, a2, b1, b2;
+        double thole1, thole2;
+        force1.getScreenedPairParameters(i, a1, a2, thole1);
+        force1.getScreenedPairParameters(i, b1, b2, thole2);
+        ASSERT_EQUAL(a1, b1);
+        ASSERT_EQUAL(a2, b2);
+        ASSERT_EQUAL(thole1, thole2);
+    }
+}
+
+int main() {
+    try {
+        registerDrudeSerializationProxies();
+        testSerialization();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+

plugins/drude/serialization/tests/TestSerializeDrudeLangevinIntegrator.cpp

+/* -------------------------------------------------------------------------- *
+ *                                OpenMMDrude                                 *
+ * -------------------------------------------------------------------------- *
+ * 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) 2013 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/Platform.h"
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/DrudeLangevinIntegrator.h"
+#include "openmm/serialization/XmlSerializer.h"
+#include <iostream>
+#include <sstream>
+
+using namespace OpenMM;
+using namespace std;
+
+extern "C" void registerDrudeSerializationProxies();
+
+void testSerialization() {
+    // Create an Integrator.
+
+    DrudeLangevinIntegrator integ1(301.1, 0.95, 10.5, 1.2, 0.001);
+    integ1.setRandomNumberSeed(18);
+
+    // Serialize and then deserialize it.
+
+    stringstream buffer;
+    XmlSerializer::serialize<DrudeLangevinIntegrator>(&integ1, "Integrator", buffer);
+    DrudeLangevinIntegrator* copy = XmlSerializer::deserialize<DrudeLangevinIntegrator>(buffer);
+
+    // Compare the two integrators to see if they are identical.
+    
+    DrudeLangevinIntegrator& integ2 = *copy;
+    ASSERT_EQUAL(integ1.getTemperature(), integ2.getTemperature());
+    ASSERT_EQUAL(integ1.getFriction(), integ2.getFriction());
+    ASSERT_EQUAL(integ1.getDrudeTemperature(), integ2.getDrudeTemperature());
+    ASSERT_EQUAL(integ1.getDrudeFriction(), integ2.getDrudeFriction());
+    ASSERT_EQUAL(integ1.getRandomNumberSeed(), integ2.getRandomNumberSeed());
+}
+
+int main() {
+    try {
+        registerDrudeSerializationProxies();
+        testSerialization();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+

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

@@ -114,15 +114,11 @@ class ForceField(object):
                         template.externalBonds.append(b)
                         template.atoms[b].externalBonds += 1
             for template in self._templates.values():
-                template.signature = _createResidueSignature([atom.element for atom in template.atoms])
-                if template.signature is None:
-                    sigString = None
+                signature = _createResidueSignature([atom.element for atom in template.atoms])
+                if signature in self._templateSignatures:
+                    self._templateSignatures[signature].append(template)
                 else:
-                    sigString = _signatureToString(template.signature)
-                if sigString in self._templateSignatures:
-                    self._templateSignatures[sigString].append(template)
-                else:
-                    self._templateSignatures[sigString] = [template]
+                    self._templateSignatures[signature] = [template]
             
             # Build sets of every atom type belonging to each class
             
@@ -162,18 +158,20 @@ class ForceField(object):
                 if typeAttrib in attrib:
                     raise ValueError('Tag specifies both a type and a class for the same atom: '+etree.tostring(node))
                 if attrib[classAttrib] not in self._atomClasses:
-                    return None # Unknown atom class
-                types.append(self._atomClasses[attrib[classAttrib]])
+                    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:
-                    return None # Unknown atom type
-                types.append([attrib[typeAttrib]])
+                    types.append(None) # Unknown atom type
+                else:
+                    types.append([attrib[typeAttrib]])
         return types
 
     def _parseTorsion(self, node):
         """Parse the node defining a torsion."""
         types = self._findAtomTypes(node, 4)
-        if types is None:
+        if None in types:
             return None
         torsion = PeriodicTorsion(types)
         attrib = node.attrib
@@ -260,20 +258,13 @@ class ForceField(object):
            particular force fields.
         Returns: the newly created System
         """
-        
-        # Record atom indices
-        
         data = ForceField._SystemData()
-        atomIndices = {}
-        for index, atom in enumerate(topology.atoms()):
-            data.atoms.append(atom)
-            atomIndices[atom] = index
+        data.atoms = list(topology.atoms())
 
         # Make a list of all bonds
         
         for bond in topology.bonds():
-            if bond[0] in atomIndices and bond[1] in atomIndices:
-                data.bonds.append(ForceField._BondData(atomIndices[bond[0]], atomIndices[bond[1]]))
+            data.bonds.append(ForceField._BondData(bond[0].index, bond[1].index))
 
         # Record which atoms are bonded to each other atom
         
@@ -294,19 +285,10 @@ class ForceField(object):
             for res in chain.residues():
                 template = None
                 matches = None
-                sig = _createResidueSignature([atom.element for atom in res.atoms()])
-                if sig is not None:
-                    signature = _signatureToString(sig)
-                    if signature in self._templateSignatures:
-                        for t in self._templateSignatures[signature]:
-                            matches = _matchResidue(res, t, bondedToAtom, atomIndices)
-                            if matches is not None:
-                                template = t
-                                break
-                if matches is None:
-                    # Check templates involving virtual sites
-                    for t in self._templateSignatures[None]:
-                        matches = _matchResidue(res, t, bondedToAtom, atomIndices)
+                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
@@ -421,7 +403,7 @@ class ForceField(object):
         
         for atom in data.virtualSites:
             site = data.virtualSites[atom]
-            index = atomIndices[atom]
+            index = atom.index
             if site.type == 'average2':
                 sys.setVirtualSite(index, mm.TwoParticleAverageSite(index+site.atoms[0], index+site.atoms[1], site.weights[0], site.weights[1]))
             elif site.type == 'average3':
@@ -436,6 +418,12 @@ class ForceField(object):
         if removeCMMotion:
             sys.addForce(mm.CMMotionRemover())
         
+        # Let generators do postprocessing
+        
+        for force in self._forces:
+            if 'postprocessSystem' in dir(force):
+                force.postprocessSystem(sys, data, args)
+        
         # Execute scripts found in the XML files.
         
         for script in self._scripts:
@@ -448,8 +436,8 @@ def _createResidueSignature(elements):
     counts = {}
     for element in elements:
         if element is None:
-            return None # This residue contains "atoms" (probably virtual sites) that should match any element
-        if element in counts:
+            pass # This residue contains "atoms" (probably virtual sites) that should match any element
+        elif element in counts:
             counts[element] += 1
         else:
             counts[element] = 1
@@ -457,33 +445,28 @@ def _createResidueSignature(elements):
     for c in counts:
         sig.append((c, counts[c]))
     sig.sort(key=lambda x: -x[0].mass)
-    return sig
-
+    
+    # Convert it to a string.
 
-def _signatureToString(signature):
-    """Convert the signature returned by _createResidueSignature() to a string."""
     s = ''
-    for element, count in signature:
+    for element, count in sig:
         s += element.symbol+str(count)
     return s
 
 
-def _matchResidue(res, template, bondedToAtom, atomIndices):
+def _matchResidue(res, template, bondedToAtom):
     """Determine whether a residue matches a template and return a list of corresponding atoms.
     
     Parameters:
      - res (Residue) The residue to check
      - template (_TemplateData) The template to compare it to
      - bondedToAtom (list) Enumerates which other atoms each atom is bonded to
-     - atomIndices (map) Maps from atoms to their indices in the System
     Returns: a list specifying which atom of the template each atom of the residue corresponds to,
     or None if it does not match the template
     """
     atoms = list(res.atoms())
     if len(atoms) != len(template.atoms):
         return None
-    residueAtomBonds = []
-    templateAtomBonds = []
     matches = len(atoms)*[0]
     hasMatch = len(atoms)*[False]
     
@@ -491,13 +474,13 @@ def _matchResidue(res, template, bondedToAtom, atomIndices):
     
     renumberAtoms = {}
     for i in range(len(atoms)):
-        renumberAtoms[atomIndices[atoms[i]]] = i
+        renumberAtoms[atoms[i].index] = i
     bondedTo = []
     externalBonds = []
     for atom in atoms:
-        bonds = [renumberAtoms[x] for x in bondedToAtom[atomIndices[atom]] if x in renumberAtoms]
+        bonds = [renumberAtoms[x] for x in bondedToAtom[atom.index] if x in renumberAtoms]
         bondedTo.append(bonds)
-        externalBonds.append(len([x for x in bondedToAtom[atomIndices[atom]] if x not in renumberAtoms]))
+        externalBonds.append(len([x for x in bondedToAtom[atom.index] if x not in renumberAtoms]))
     if _findAtomMatches(atoms, template, bondedTo, externalBonds, matches, hasMatch, 0):
         return matches
     return None
@@ -508,9 +491,10 @@ def _findAtomMatches(atoms, template, bondedTo, externalBonds, matches, hasMatch
     if position == len(atoms):
         return True
     elem = atoms[position].element
+    name = atoms[position].name
     for i in range(len(atoms)):
         atom = template.atoms[i]
-        if (atom.element == elem or atom.element is None) and not hasMatch[i] and len(atom.bondedTo) == len(bondedTo[position]) and atom.externalBonds == externalBonds[position]:
+        if (atom.element == elem or (atom.element is None and atom.name == name)) and not hasMatch[i] and len(atom.bondedTo) == len(bondedTo[position]) and atom.externalBonds == externalBonds[position]:
             # See if the bonds for this identification are consistent
             
             allBondsMatch = all((bonded > position or matches[bonded] in atom.bondedTo for bonded in bondedTo[position]))
@@ -546,7 +530,7 @@ class HarmonicBondGenerator:
         ff._forces.append(generator)
         for bond in element.findall('Bond'):
             types = ff._findAtomTypes(bond, 2)
-            if types is not None:
+            if None not in types:
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
                 generator.length.append(float(bond.attrib['length']))
@@ -594,7 +578,7 @@ class HarmonicAngleGenerator:
         ff._forces.append(generator)
         for angle in element.findall('Angle'):
             types = ff._findAtomTypes(angle, 3)
-            if types is not None:
+            if None not in types:
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
                 generator.types3.append(types[2])
@@ -773,11 +757,11 @@ class RBTorsionGenerator:
         ff._forces.append(generator)
         for torsion in element.findall('Proper'):
             types = ff._findAtomTypes(torsion, 4)
-            if types is not None:
+            if None not in types:
                 generator.proper.append(RBTorsion(types, [float(torsion.attrib['c'+str(i)]) for i in range(6)]))
         for torsion in element.findall('Improper'):
             types = ff._findAtomTypes(torsion, 4)
-            if types is not None:
+            if None not in types:
                 generator.improper.append(RBTorsion(types, [float(torsion.attrib['c'+str(i)]) for i in range(6)]))
     
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
@@ -875,7 +859,7 @@ class CMAPTorsionGenerator:
             generator.maps.append(values)
         for torsion in element.findall('Torsion'):
             types = ff._findAtomTypes(torsion, 5)
-            if types is not None:
+            if None not in types:
                 generator.torsions.append(CMAPTorsion(types, int(torsion.attrib['map'])))
     
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
@@ -956,7 +940,7 @@ class NonbondedGenerator:
                 raise ValueError('Found multiple NonbondedForce tags with different 1-4 scales') 
         for atom in element.findall('Atom'):
             types = ff._findAtomTypes(atom, 1)
-            if types is not None:
+            if None not in types:
                 values = (float(atom.attrib['charge']), float(atom.attrib['sigma']), float(atom.attrib['epsilon']))
                 for t in types[0]:
                     generator.typeMap[t] = values
@@ -977,7 +961,14 @@ class NonbondedGenerator:
                 force.addParticle(values[0], values[1], values[2])
             else:
                 raise ValueError('No nonbonded parameters defined for atom type '+t)
-        # Create exceptions based on bonds and virtual sites.
+        force.setNonbondedMethod(methodMap[nonbondedMethod])
+        force.setCutoffDistance(nonbondedCutoff)
+        if 'ewaldErrorTolerance' in args:
+            force.setEwaldErrorTolerance(args['ewaldErrorTolerance'])
+        sys.addForce(force)
+    
+    def postprocessSystem(self, sys, data, args):
+        # Create exceptions based on bonds, virtual sites, and Drude particles.
         bondIndices = []
         for bond in data.bonds:
             bondIndices.append((bond.atom1, bond.atom2))
@@ -986,12 +977,26 @@ class NonbondedGenerator:
                 site = sys.getVirtualSite(i)
                 for j in range(site.getNumParticles()):
                     bondIndices.append((i, site.getParticle(j)))
-        force.createExceptionsFromBonds(bondIndices, self.coulomb14scale, self.lj14scale)
-        force.setNonbondedMethod(methodMap[nonbondedMethod])
-        force.setCutoffDistance(nonbondedCutoff)
-        if 'ewaldErrorTolerance' in args:
-            force.setEwaldErrorTolerance(args['ewaldErrorTolerance'])
-        sys.addForce(force)
+        drude = [f for f in sys.getForces() if isinstance(f, mm.DrudeForce)]
+        if len(drude) > 0:
+            drude = drude[0]
+            # For purposes of creating exceptions, a Drude particle is "bonded" to anything
+            # its parent atom is bonded to.
+            drudeMap = {}
+            for i in range(drude.getNumParticles()):
+                params = drude.getParticleParameters(i)
+                drudeMap[params[1]] = params[0]
+            for atom1, atom2 in bondIndices:
+                drude1 = drudeMap[atom1] if atom1 in drudeMap else None
+                drude2 = drudeMap[atom2] if atom2 in drudeMap else None
+                if drude1 is not None:
+                    bondIndices.append((drude1, atom2))
+                    if drude2 is not None:
+                        bondIndices.append((drude1, drude2))
+                if drude2 is not None:
+                    bondIndices.append((atom1, drude2))
+        nonbonded = [f for f in sys.getForces() if isinstance(f, mm.NonbondedForce)][0]
+        nonbonded.createExceptionsFromBonds(bondIndices, self.coulomb14scale, self.lj14scale)        
 
 parsers["NonbondedForce"] = NonbondedGenerator.parseElement
 
@@ -1014,7 +1019,7 @@ class GBSAOBCGenerator:
             generator = existing[0]
         for atom in element.findall('Atom'):
             types = ff._findAtomTypes(atom, 1)
-            if types is not None:
+            if None not in types:
                 values = (float(atom.attrib['charge']), float(atom.attrib['radius']), float(atom.attrib['scale']))
                 for t in types[0]:
                     generator.typeMap[t] = values
@@ -1071,7 +1076,7 @@ class GBVIGenerator:
         ff._forces.append(generator)
         for atom in element.findall('Atom'):
             types = ff._findAtomTypes(atom, 1)
-            if types is not None:
+            if None not in types:
                 values = (float(atom.attrib['charge']), float(atom.attrib['radius']), float(atom.attrib['gamma']))
                 for t in types[0]:
                     generator.typeMap[t] = values
@@ -1154,7 +1159,7 @@ class CustomBondGenerator:
             generator.perBondParams.append(param.attrib['name'])
         for bond in element.findall('Bond'):
             types = ff._findAtomTypes(bond, 2)
-            if types is not None:
+            if None not in types:
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
                 generator.paramValues.append([float(bond.attrib[param]) for param in generator.perBondParams])
@@ -1202,7 +1207,7 @@ class CustomAngleGenerator:
             generator.perAngleParams.append(param.attrib['name'])
         for angle in element.findall('Angle'):
             types = ff._findAtomTypes(angle, 3)
-            if types is not None:
+            if None not in types:
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
                 generator.types3.append(types[2])
@@ -1263,11 +1268,11 @@ class CustomTorsionGenerator:
             generator.perTorsionParams.append(param.attrib['name'])
         for torsion in element.findall('Proper'):
             types = ff._findAtomTypes(torsion, 4)
-            if types is not None:
+            if None not in types:
                 generator.proper.append(CustomTorsion(types, [float(torsion.attrib[param]) for param in generator.perTorsionParams]))
         for torsion in element.findall('Improper'):
             types = ff._findAtomTypes(torsion, 4)
-            if types is not None:
+            if None not in types:
                 generator.improper.append(CustomTorsion(types, [float(torsion.attrib[param]) for param in generator.perTorsionParams]))
     
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
@@ -1354,7 +1359,7 @@ class CustomGBGenerator:
             generator.perParticleParams.append(param.attrib['name'])
         for atom in element.findall('Atom'):
             types = ff._findAtomTypes(atom, 1)
-            if types is not None:
+            if None not in types:
                 values = [float(atom.attrib[param]) for param in generator.perParticleParams]
                 for t in types[0]:
                     generator.typeMap[t] = values
@@ -1456,7 +1461,7 @@ class AmoebaBondGenerator:
         forceField._forces.append(generator)
         for bond in element.findall('Bond'):
             types = forceField._findAtomTypes(bond, 2)
-            if types is not None:
+            if None not in types:
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
                 generator.length.append(float(bond.attrib['length']))
@@ -1569,7 +1574,7 @@ class AmoebaAngleGenerator:
         forceField._forces.append(generator)
         for angle in element.findall('Angle'):
             types = forceField._findAtomTypes(angle, 3)
-            if types is not None:
+            if None not in types:
 
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
@@ -2054,7 +2059,7 @@ class AmoebaTorsionGenerator:
 
         for torsion in element.findall('Torsion'):
             types = forceField._findAtomTypes(torsion, 4)
-            if types is not None:
+            if None not in types:
 
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
@@ -2165,7 +2170,7 @@ class AmoebaPiTorsionGenerator:
 
         for piTorsion in element.findall('PiTorsion'):
             types = forceField._findAtomTypes(piTorsion, 2)
-            if types is not None:
+            if None not in types:
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
                 generator.k.append(float(piTorsion.attrib['k']))
@@ -2284,7 +2289,7 @@ class AmoebaTorsionTorsionGenerator:
 
         for torsionTorsion in element.findall('TorsionTorsion'):
             types = forceField._findAtomTypes(torsionTorsion, 5)
-            if types is not None:
+            if None not in types:
 
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
@@ -2521,7 +2526,7 @@ class AmoebaStretchBendGenerator:
 
         for stretchBend in element.findall('StretchBend'):
             types = forceField._findAtomTypes(stretchBend, 3)
-            if types is not None:
+            if None not in types:
 
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
@@ -2678,7 +2683,7 @@ class AmoebaVdwGenerator:
 
         for atom in element.findall('Vdw'):
             types = forceField._findAtomTypes(atom, 1)
-            if types is not None:
+            if None not in types:
 
                 values = [float(atom.attrib['sigma']), float(atom.attrib['epsilon']), float(atom.attrib['reduction'])]
 
@@ -2951,7 +2956,7 @@ class AmoebaMultipoleGenerator:
 
         for atom in element.findall('Multipole'):
             types = forceField._findAtomTypes(atom, 1)
-            if types is not None:
+            if None not in types:
 
                 # k-indices not provided default to 0
 
@@ -3008,7 +3013,7 @@ class AmoebaMultipoleGenerator:
  
         for atom in element.findall('Polarize'):
             types = forceField._findAtomTypes(atom, 1)
-            if types is not None:
+            if None not in types:
 
                 classIndex = atom.attrib['type']
                 polarizability = float(atom.attrib['polarizability'])
@@ -3536,7 +3541,7 @@ class AmoebaWcaDispersionGenerator:
 
         for atom in element.findall('WcaDispersion'):
             types = forceField._findAtomTypes(atom, 1)
-            if types is not None:
+            if None not in types:
 
                 values = [float(atom.attrib['radius']), float(atom.attrib['epsilon'])]
                 for t in types[0]:
@@ -3880,7 +3885,7 @@ class AmoebaUreyBradleyGenerator:
         forceField._forces.append(generator)
         for bond in element.findall('UreyBradley'):
             types = forceField._findAtomTypes(bond, 3)
-            if types is not None:
+            if None not in types:
 
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
@@ -3924,3 +3929,87 @@ class AmoebaUreyBradleyGenerator:
 parsers["AmoebaUreyBradleyForce"] = AmoebaUreyBradleyGenerator.parseElement
 
 #=============================================================================================
+
+
+## @private
+class DrudeGenerator:
+    """A DrudeGenerator constructs a DrudeForce."""
+    
+    def __init__(self):
+        self.typeMap = {}
+
+    @staticmethod
+    def parseElement(element, ff):
+        existing = [f for f in ff._forces if isinstance(f, DrudeGenerator)]
+        if len(existing) == 0:
+            generator = DrudeGenerator()
+            ff._forces.append(generator)
+        else:
+            # Multiple <DrudeForce> tags were found, probably in different files.  Simply add more types to the existing one.
+            generator = existing[0]
+        for particle in element.findall('Particle'):
+            types = ff._findAtomTypes(particle, 5)
+            if None not in types[:2]:
+                aniso12 = 0.0
+                aniso34 = 0.0
+                if 'aniso12' in particle.attrib:
+                    aniso12 = float(particle.attrib['aniso12'])
+                if 'aniso34' in particle.attrib:
+                    aniso34 = float(particle.attrib['aniso34'])
+                values = (types[1], types[2], types[3], types[4], float(particle.attrib['charge']), float(particle.attrib['polarizability']), aniso12, aniso34, float(particle.attrib['thole']))
+                for t in types[0]:
+                    generator.typeMap[t] = values
+    
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
+        force = mm.DrudeForce()
+        if not any(isinstance(f, mm.NonbondedForce) for f in sys.getForces()):
+            raise ValueError('<DrudeForce> must come after <NonbondedForce> in XML file')
+        
+        # Add Drude particles.
+        
+        drudeMap = {}
+        parentMap = {}
+        for atom in data.atoms:
+            t = data.atomType[atom]
+            if t in self.typeMap:
+                # Find other atoms in the residue that affect the Drude particle.
+                p = [-1, -1, -1, -1]
+                values = self.typeMap[t]
+                for atom2 in atom.residue.atoms():
+                    type2 = data.atomType[atom2]
+                    if type2 in values[0]:
+                        p[0] = atom2.index
+                    elif values[1] is not None and type2 in values[1]:
+                        p[1] = atom2.index
+                    elif values[2] is not None and type2 in values[2]:
+                        p[2] = atom2.index
+                    elif values[3] is not None and type2 in values[3]:
+                        p[3] = atom2.index
+                drudeIndex = force.addParticle(atom.index, p[0], p[1], p[2], p[3], values[4], values[5], values[6], values[7])
+                drudeMap[atom.index] = p[0]
+                parentMap[p[0]] = (atom.index, drudeIndex)
+        sys.addForce(force)
+    
+    def postprocessSystem(self, sys, data, args):
+        # For every nonbonded exclusion between Drude particles, add a screened pair.
+        
+        drude = [f for f in sys.getForces() if isinstance(f, mm.DrudeForce)][0]
+        nonbonded = [f for f in sys.getForces() if isinstance(f, mm.NonbondedForce)][0]
+        particleMap = {}
+        for i in range(drude.getNumParticles()):
+            particleMap[drude.getParticleParameters(i)[0]] = i
+        for i in range(nonbonded.getNumExceptions()):
+            (particle1, particle2, charge, sigma, epsilon) = nonbonded.getExceptionParameters(i)
+            if charge == 0 and epsilon == 0:
+                # This is an exclusion.
+                if particle1 in particleMap and particle2 in particleMap:
+                    # It connects two Drude particles, so add a screened pair.
+                    drude1 = particleMap[particle1]
+                    drude2 = particleMap[particle2]
+                    type1 = data.atomType[data.atoms[drude1]]
+                    type2 = data.atomType[data.atoms[drude2]]
+                    thole1 = self.typeMap[type1][8]
+                    thole2 = self.typeMap[type2][8]
+                    drude.addScreenedPair(drude1, drude2, thole1+thole2)
+
+parsers["DrudeForce"] = DrudeGenerator.parseElement
\ No newline at end of file

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

@@ -33,11 +33,12 @@ from __future__ import division
 __author__ = "Peter Eastman"
 __version__ = "1.0"
 
-from simtk.openmm.app import Topology, PDBFile
-from simtk.openmm.app.forcefield import HAngles
+from simtk.openmm.app import Topology, PDBFile, ForceField
+from simtk.openmm.app.forcefield import HAngles, _createResidueSignature, _matchResidue, DrudeGenerator
+from simtk.openmm.app.topology import Residue
 from simtk.openmm.vec3 import Vec3
 from simtk.openmm import System, Context, NonbondedForce, VerletIntegrator, LocalEnergyMinimizer
-from simtk.unit import nanometer, molar, elementary_charge, amu, gram, liter, degree, sqrt, acos, is_quantity, dot, norm
+from simtk.unit import nanometer, molar, elementary_charge, amu, gram, liter, degree, sqrt, acos, is_quantity, dot, norm, sum
 import element as elem
 import os
 import random
@@ -516,7 +517,7 @@ class Modeller(object):
                     terminal = hydrogen.attrib['terminal']
                 data.hydrogens.append(Modeller._Hydrogen(hydrogen.attrib['name'], hydrogen.attrib['parent'], maxph, atomVariants, terminal))
     
-    def addHydrogens(self, forcefield, pH=7.0, variants=None, hydrogenDefinitions=None):
+    def addHydrogens(self, forcefield, pH=7.0, variants=None):
         """Add missing hydrogens to the model.
         
         Some residues can exist in multiple forms depending on the pH and properties of the local environment.  These
@@ -763,3 +764,160 @@ class Modeller(object):
         self.topology = newTopology
         self.positions = context.getState(getPositions=True).getPositions()
         return actualVariants
+    
+    def addExtraParticles(self, forcefield):
+        """Add missing extra particles to the model that are required by a force field.
+        
+        Some force fields use "extra particles" that do not represent actual atoms, but still need to be included in
+        the System.  Examples include lone pairs, Drude particles, and the virtual sites used in some water models
+        to adjust the charge distribution.  Extra particles can be recognized by the fact that their element is None.
+        
+        This method is primarily used to add extra particles, but it can also remove them.  It tries to match every
+        residue in the Topology to a template in the force field.  If there is no match, it will both add and remove
+        extra particles as necessary to make it match.
+        
+        Parameters:
+         - forcefield (ForceField) the ForceField defining what extra particles should be present
+        """
+        # Create copies of all residue templates that have had all extra points removed.
+        
+        templatesNoEP = {}
+        for resName, template in forcefield._templates.iteritems():
+            if any(atom.element is None for atom in template.atoms):
+                index = 0
+                newIndex = {}
+                newTemplate = ForceField._TemplateData(resName)
+                for i, atom in enumerate(template.atoms):
+                    if atom.element is not None:
+                        newIndex[i] = index
+                        index += 1
+                        newTemplate.atoms.append(ForceField._TemplateAtomData(atom.name, atom.type, atom.element))
+                for b1, b2 in template.bonds:
+                    if b1 in newIndex and b2 in newIndex:
+                        newTemplate.bonds.append((newIndex[b1], newIndex[b2]))
+                        newTemplate.atoms[newIndex[b1]].bondedTo.append(newIndex[b2])
+                        newTemplate.atoms[newIndex[b2]].bondedTo.append(newIndex[b1])
+                for b in template.externalBonds:
+                    if b in newIndex:
+                        newTemplate.externalBonds.append(newIndex[b])
+                templatesNoEP[template] = newTemplate
+        
+        # Record which atoms are bonded to each other atom, with and without extra particles.
+        
+        bondedToAtom = []
+        bondedToAtomNoEP = []
+        for atom in self.topology.atoms():
+            bondedToAtom.append(set())
+            if atom.element is not None:
+                bondedToAtomNoEP.append(set())
+        for atom1, atom2 in self.topology.bonds():
+            bondedToAtom[atom1.index].add(atom2.index)
+            bondedToAtom[atom2.index].add(atom1.index)
+            if atom1.element is not None and atom2.element is not None:
+                bondedToAtomNoEP[atom1.index].add(atom2.index)
+                bondedToAtomNoEP[atom2.index].add(atom1.index)
+        
+        # If the force field has a DrudeForce, record the types of Drude particles and their parents since we'll
+        # need them for picking particle positions.
+        
+        drudeTypeMap = {}
+        for force in forcefield._forces:
+            if isinstance(force, DrudeGenerator):
+                for type in force.typeMap:
+                    drudeTypeMap[type] = force.typeMap[type][0]
+        
+        # Create the new Topology.
+        
+        newTopology = Topology()
+        newTopology.setUnitCellDimensions(deepcopy(self.topology.getUnitCellDimensions()))
+        newAtoms = {}
+        newPositions = []*nanometer
+        for chain in self.topology.chains():
+            newChain = newTopology.addChain()
+            for residue in chain.residues():
+                newResidue = newTopology.addResidue(residue.name, newChain)
+                
+                # Look for a matching template.
+                
+                matchFound = False
+                signature = _createResidueSignature([atom.element for atom in residue.atoms()])
+                if signature in forcefield._templateSignatures:
+                    for t in forcefield._templateSignatures[signature]:
+                        if _matchResidue(residue, t, bondedToAtom) is not None:
+                            matchFound = True
+                if matchFound:
+                    # Just copy the residue over.
+                    
+                    for atom in residue.atoms():
+                        newAtom = newTopology.addAtom(atom.name, atom.element, newResidue)
+                        newAtoms[atom] = newAtom
+                        newPositions.append(deepcopy(self.positions[atom.index]))
+                else:
+                    # There's no matching template.  Try to find one that matches based on everything except
+                    # extra points.
+                    
+                    template = None
+                    residueNoEP = Residue(residue.name, residue.index, residue.chain)
+                    residueNoEP._atoms = [atom for atom in residue.atoms() if atom.element is not None]
+                    if signature in forcefield._templateSignatures:
+                        for t in forcefield._templateSignatures[signature]:
+                            if t in templatesNoEP:
+                                matches = _matchResidue(residueNoEP, templatesNoEP[t], bondedToAtomNoEP)
+                                if matches is not None:
+                                    template = t;
+                                    # Record the corresponding atoms.
+                                    matchingAtoms = {}
+                                    for atom, match in zip(residueNoEP.atoms(), matches):
+                                        templateAtomName = t.atoms[match].name
+                                        for templateAtom in template.atoms:
+                                            if templateAtom.name == templateAtomName:
+                                                matchingAtoms[templateAtom] = atom
+                                    break
+                    if template is None:
+                        raise ValueError('Residue %d (%s) does not match any template defined by the ForceField.' % (residue.index+1, residue.name))
+                    
+                    # Add the regular atoms.
+                    
+                    for atom in residue.atoms():
+                        if atom.element is not None:
+                            newAtoms[atom] = newTopology.addAtom(atom.name, atom.element, newResidue)
+                            newPositions.append(deepcopy(self.positions[atom.index]))
+                    
+                    # Add the extra points.
+                    
+                    templateAtomPositions = len(template.atoms)*[None]
+                    for index, atom in enumerate(template.atoms):
+                        if atom in matchingAtoms:
+                            templateAtomPositions[index] = self.positions[matchingAtoms[atom].index]
+                    for index, atom in enumerate(template.atoms):
+                        if atom.element is None:
+                            newTopology.addAtom(atom.name, None, newResidue)
+                            position = None
+                            for site in template.virtualSites:
+                                if site.index == index:
+                                    # This is a virtual site.  Compute its position by the correct rule.
+                                    
+                                    if site.type == 'average2':
+                                        position = site.weights[0]*templateAtomPositions[index+site.atoms[0]] + site.weights[1]*templateAtomPositions[index+site.atoms[1]]
+                                    elif site.type == 'average3':
+                                        position = site.weights[0]*templateAtomPositions[index+site.atoms[0]] + site.weights[1]*templateAtomPositions[index+site.atoms[1]] + site.weights[2]*templateAtomPositions[index+site.atoms[2]]
+                                    elif site.type == 'outOfPlane':
+                                        position = site.weights[0]*templateAtomPositions[index+site.atoms[0]] + site.weights[1]*templateAtomPositions[index+site.atoms[1]] + site.weights[2]*templateAtomPositions[index+site.atoms[2]]
+                            if position is None and atom.type in drudeTypeMap:
+                                # This is a Drude particle.  Put it on top of its parent atom.
+                                
+                                for atom2, pos in zip(template.atoms, templateAtomPositions):
+                                    if atom2.type in drudeTypeMap[atom.type]:
+                                        position = deepcopy(pos)
+                            if position is None:
+                                # We couldn't figure out the correct position.  As a wild guess, just put it at the center of the residue
+                                # and hope that energy minimization will fix it.
+                                
+                                knownPositions = [x for x in templateAtomPositions if x is not None]
+                                position = sum(knownPositions)/len(knownPositions)
+                            newPositions.append(position)
+        for bond in self.topology.bonds():
+            if bond[0] in newAtoms and bond[1] in newAtoms:
+                newTopology.addBond(newAtoms[bond[0]], newAtoms[bond[1]])
+        self.topology = newTopology
+        self.positions = newPositions
\ No newline at end of file