Merge branch 'master' into membrane

docs-source/usersguide/application.rst

@@ -1794,7 +1794,7 @@ Here is the definition of the :class:`ForceReporter` class:
 
             def describeNextReport(self, simulation):
                 steps = self._reportInterval - simulation.currentStep%self._reportInterval
-                return (steps, False, False, True, False)
+                return (steps, False, False, True, False, None)
 
             def report(self, simulation, state):
                 forces = state.getForces().value_in_unit(kilojoules/mole/nanometer)
@@ -1814,7 +1814,7 @@ We then have two methods that every reporter must implement:
 :meth:`describeNextReport()` and :meth:`report()`.  A Simulation object
 periodically calls :meth:`describeNextReport()` on each of its reporters to
 find out when that reporter will next generate a report, and what information
-will be needed to generate it.  The return value should be a five element tuple,
+will be needed to generate it.  The return value should be a six element tuple,
 whose elements are as follows:
 
 * The number of time steps until the next report.  We calculate this as
@@ -1825,6 +1825,9 @@ whose elements are as follows:
 * Whether the next report will need particle velocities.
 * Whether the next report will need forces.
 * Whether the next report will need energies.
+* Whether the positions should be wrapped to the periodic box.  If None, it will
+  automatically decide whether to wrap positions based on whether the System uses
+  periodic boundary conditions.
 
 
 When the time comes for the next scheduled report, the :class:`Simulation` calls

openmmapi/include/openmm/Context.h

@@ -217,15 +217,25 @@ public:
      */
     void computeVirtualSites();
     /**
-     * When a Context is created, it may cache information about the System being simulated
+     * When a Context is created, it caches information about the System being simulated
      * and the Force objects contained in it.  This means that, if the System or Forces are then
-     * modified, the Context might not see all of the changes.  Call reinitialize() to force
+     * modified, the Context does not see the changes.  Call reinitialize() to force
      * the Context to rebuild its internal representation of the System and pick up any changes
      * that have been made.
      * 
      * This is an expensive operation, so you should try to avoid calling it too frequently.
+     * Most Force classes have an updateParametersInContext() method that provides a less expensive
+     * way of updating certain types of information.  However, this method is the only way to
+     * make some types of changes, so it is sometimes necessary to call it.
+     * 
+     * By default, reinitializing a Context causes all state information (positions, velocities,
+     * etc.) to be discarded.  You can optionally tell it to try to preserve state information.
+     * It does this by internally creating a checkpoint, then reinitializing the Context, then
+     * loading the checkpoint.  Be aware that if the System has changed in a way that prevents
+     * the checkpoint from being loaded (such as changing the number of particles), this will
+     * throw an exception and the state information will be lost.
      */
-    void reinitialize();
+    void reinitialize(bool preserveState=false);
     /**
      * Create a checkpoint recording the current state of the Context.  This should be treated
      * as an opaque block of binary data.  See loadCheckpoint() for more details.

openmmapi/src/Context.cpp

@@ -36,6 +36,8 @@
 #include "SimTKOpenMMRealType.h"
 #include "sfmt/SFMT.h"
 #include <cmath>
+#include <iostream>
+#include <sstream>
 
 using namespace OpenMM;
 using namespace std;
@@ -242,14 +244,19 @@ void Context::computeVirtualSites() {
     impl->computeVirtualSites();
 }
 
-void Context::reinitialize() {
+void Context::reinitialize(bool preserveState) {
     const System& system = impl->getSystem();
     Integrator& integrator = impl->getIntegrator();
     Platform& platform = impl->getPlatform();
+    stringstream checkpoint(ios_base::out | ios_base::in | ios_base::binary);
+    if (preserveState)
+        createCheckpoint(checkpoint);
     integrator.cleanup();
     delete impl;
     impl = new ContextImpl(*this, system, integrator, &platform, properties);
     impl->initialize();
+    if (preserveState)
+        loadCheckpoint(checkpoint);
 }
 
 void Context::createCheckpoint(ostream& stream) {

platforms/cpu/src/CpuLangevinDynamics.cpp

@@ -112,7 +112,6 @@ void CpuLangevinDynamics::threadUpdate2(int threadIndex) {
 
     for (int i = start; i < end; i++) {
         if (inverseMasses[i] != 0.0) {
-            double sqrtInvMass = sqrt(inverseMasses[i]);
             xPrime[i] = atomCoordinates[i]+velocities[i]*dt;
         }
    }

platforms/cuda/include/CudaKernels.h

@@ -778,7 +778,7 @@ private:
     std::vector<CudaArray*> tabulatedFunctions;
     double longRangeCoefficient;
     std::vector<double> longRangeCoefficientDerivs;
-    bool hasInitializedLongRangeCorrection, hasInitializedKernel;
+    bool hasInitializedLongRangeCorrection, hasInitializedKernel, hasParamDerivs;
     int numGroupThreadBlocks;
     CustomNonbondedForce* forceCopy;
     const System& system;

platforms/cuda/src/CudaContext.cpp

@@ -236,12 +236,6 @@ CudaContext::CudaContext(const System& system, int deviceIndex, bool useBlocking
             minor = 3;
         }
     }
-    if (major == 7) {
-        // Don't generate Volta-specific code until we've made the changes needed
-        // to support it properly.
-        major = 6;
-        minor = 0;
-    }
     gpuArchitecture = intToString(major)+intToString(minor);
     computeCapability = major+0.1*minor;
 
@@ -263,6 +257,16 @@ CudaContext::CudaContext(const System& system, int deviceIndex, bool useBlocking
     int multiprocessors;
     CHECK_RESULT(cuDeviceGetAttribute(&multiprocessors, CU_DEVICE_ATTRIBUTE_MULTIPROCESSOR_COUNT, device));
     numThreadBlocks = numThreadBlocksPerComputeUnit*multiprocessors;
+    if (cudaDriverVersion >= 9000) {
+        compilationDefines["SYNC_WARPS"] = "__syncwarp();";
+        compilationDefines["SHFL(var, srcLane)"] = "__shfl_sync(0xffffffff, var, srcLane);";
+        compilationDefines["BALLOT(var)"] = "__ballot_sync(0xffffffff, var);";
+    }
+    else {
+        compilationDefines["SYNC_WARPS"] = "";
+        compilationDefines["SHFL(var, srcLane)"] = "__shfl(var, srcLane);";
+        compilationDefines["BALLOT(var)"] = "__ballot(var);";
+    }
     if (useDoublePrecision) {
         posq = CudaArray::create<double4>(*this, paddedNumAtoms, "posq");
         velm = CudaArray::create<double4>(*this, paddedNumAtoms, "velm");

platforms/cuda/src/CudaKernels.cpp

@@ -2664,6 +2664,7 @@ void CudaCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNonbo
     
     // Create the kernel.
     
+    hasParamDerivs = (force.getNumEnergyParameterDerivatives() > 0);
     map<string, string> replacements;
     replacements["COMPUTE_INTERACTION"] = interactionSource;
     const string suffixes[] = {"x", "y", "z", "w"};
@@ -2687,6 +2688,8 @@ void CudaCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNonbo
         args << ", const " << tableTypes[i]<< "* __restrict__ table" << i;
     if (globals != NULL)
         args<<", const float* __restrict__ globals";
+    if (hasParamDerivs)
+        args << ", mixed* __restrict__ energyParamDerivs";
     replacements["PARAMETER_ARGUMENTS"] = args.str();
     stringstream load1;
     for (int i = 0; i < (int) buffers.size(); i++)
@@ -2718,6 +2721,19 @@ void CudaCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNonbo
         }
     }
     replacements["LOAD_ATOM2_PARAMETERS"] = load2.str();
+    stringstream initDerivs, saveDerivs;
+    const vector<string>& allParamDerivNames = cu.getEnergyParamDerivNames();
+    int numDerivs = allParamDerivNames.size();
+    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
+        string paramName = force.getEnergyParameterDerivativeName(i);
+        string derivVariable = cu.getNonbondedUtilities().addEnergyParameterDerivative(paramName);
+        initDerivs<<"mixed "<<derivVariable<<" = 0;\n";
+        for (int index = 0; index < numDerivs; index++)
+            if (allParamDerivNames[index] == paramName)
+                saveDerivs<<"energyParamDerivs[(blockIdx.x*blockDim.x+threadIdx.x)*"<<numDerivs<<"+"<<index<<"] += "<<derivVariable<<";\n";
+    }
+    replacements["INIT_DERIVATIVES"] = initDerivs.str();
+    replacements["SAVE_DERIVATIVES"] = saveDerivs.str();
     map<string, string> defines;
     if (force.getNonbondedMethod() != CustomNonbondedForce::NoCutoff)
         defines["USE_CUTOFF"] = "1";
@@ -2779,6 +2795,8 @@ double CudaCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool in
                 interactionGroupArgs.push_back(&function->getDevicePointer());
             if (globals != NULL)
                 interactionGroupArgs.push_back(&globals->getDevicePointer());
+            if (hasParamDerivs)
+                interactionGroupArgs.push_back(&cu.getEnergyParamDerivBuffer().getDevicePointer());
         }
         int forceThreadBlockSize = cu.getNonbondedUtilities().getForceThreadBlockSize();
         cu.executeKernel(interactionGroupKernel, &interactionGroupArgs[0], numGroupThreadBlocks*forceThreadBlockSize, forceThreadBlockSize);

platforms/cuda/src/kernels/customManyParticle.cu

@@ -225,7 +225,7 @@ extern "C" __global__ void findNeighbors(real4 periodicBoxSize, real4 invPeriodi
             
             // Loop over any blocks we identified as potentially containing neighbors.
             
-            int includeBlockFlags = __ballot(includeBlock2);
+            int includeBlockFlags = BALLOT(includeBlock2);
             while (includeBlockFlags != 0) {
                 int i = __ffs(includeBlockFlags)-1;
                 includeBlockFlags &= includeBlockFlags-1;

platforms/cuda/src/kernels/customNonbondedGroups.cu

@@ -17,6 +17,7 @@ extern "C" __global__ void computeInteractionGroups(
     const unsigned int tgx = threadIdx.x & (TILE_SIZE-1); // index within the warp
     const unsigned int tbx = threadIdx.x - tgx;           // block warpIndex
     mixed energy = 0;
+    INIT_DERIVATIVES
     __shared__ AtomData localData[LOCAL_MEMORY_SIZE];
 
     const unsigned int startTile = FIRST_TILE+warp*(LAST_TILE-FIRST_TILE)/totalWarps;
@@ -58,6 +59,7 @@ extern "C" __global__ void computeInteractionGroups(
                 LOAD_ATOM2_PARAMETERS
                 real dEdR = 0.0f;
                 real tempEnergy = 0.0f;
+                const real interactionScale = 1.0f;
                 COMPUTE_INTERACTION
                 energy += tempEnergy;
                 delta *= dEdR;
@@ -82,4 +84,5 @@ extern "C" __global__ void computeInteractionGroups(
         atomicAdd(&forceBuffers[atom2+2*PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (localData[threadIdx.x].fz*0x100000000)));
     }
     energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] += energy;
+    SAVE_DERIVATIVES
 }

platforms/cuda/src/kernels/findInteractingBlocks.cu

@@ -115,7 +115,7 @@ __device__ int saveSinglePairs(int x, int* atoms, int* flags, int length, unsign
         int atom = atoms[i];
         int flag = flags[i];
         bool include = (i < length && __popc(flags[i]) > MAX_BITS_FOR_PAIRS);
-        int includeFlags = __ballot(include);
+        int includeFlags = BALLOT(include);
         if (include) {
             int index = numCompacted+__popc(includeFlags&warpMask);
             atoms[index] = atom;
@@ -271,7 +271,7 @@ extern "C" __global__ void findBlocksWithInteractions(real4 periodicBoxSize, rea
             
             // Loop over any blocks we identified as potentially containing neighbors.
             
-            int includeBlockFlags = __ballot(includeBlock2);
+            int includeBlockFlags = BALLOT(includeBlock2);
             while (includeBlockFlags != 0) {
                 int i = __ffs(includeBlockFlags)-1;
                 includeBlockFlags &= includeBlockFlags-1;
@@ -291,7 +291,7 @@ extern "C" __global__ void findBlocksWithInteractions(real4 periodicBoxSize, rea
 #ifdef USE_PERIODIC
                 APPLY_PERIODIC_TO_DELTA(atomDelta)
 #endif
-                int atomFlags = ballot(atomDelta.x*atomDelta.x+atomDelta.y*atomDelta.y+atomDelta.z*atomDelta.z < (PADDED_CUTOFF+blockCenterY.w)*(PADDED_CUTOFF+blockCenterY.w));
+                int atomFlags = BALLOT(atomDelta.x*atomDelta.x+atomDelta.y*atomDelta.y+atomDelta.z*atomDelta.z < (PADDED_CUTOFF+blockCenterY.w)*(PADDED_CUTOFF+blockCenterY.w));
                 int interacts = 0;
                 if (atom2 < NUM_ATOMS && atomFlags != 0) {
                     int first = __ffs(atomFlags)-1;
@@ -317,7 +317,7 @@ extern "C" __global__ void findBlocksWithInteractions(real4 periodicBoxSize, rea
                 
                 // Add any interacting atoms to the buffer.
                 
-                int includeAtomFlags = __ballot(interacts);
+                int includeAtomFlags = BALLOT(interacts);
                 if (interacts) {
                     int index = neighborsInBuffer+__popc(includeAtomFlags&warpMask);
                     buffer[index] = atom2;

platforms/cuda/src/kernels/nonbonded.cu

@@ -15,22 +15,22 @@ typedef struct {
 #ifdef ENABLE_SHUFFLE
 //support for 64 bit shuffles
 static __inline__ __device__ float real_shfl(float var, int srcLane) {
-    return __shfl(var, srcLane);
+    return SHFL(var, srcLane);
 }
 
 static __inline__ __device__ double real_shfl(double var, int srcLane) {
     int hi, lo;
     asm volatile("mov.b64 { %0, %1 }, %2;" : "=r"(lo), "=r"(hi) : "d"(var));
-    hi = __shfl(hi, srcLane);
-    lo = __shfl(lo, srcLane);
+    hi = SHFL(hi, srcLane);
+    lo = SHFL(lo, srcLane);
     return __hiloint2double( hi, lo );
 }
 
 static __inline__ __device__ long long real_shfl(long long var, int srcLane) {
     int hi, lo;
     asm volatile("mov.b64 { %0, %1 }, %2;" : "=r"(lo), "=r"(hi) : "l"(var));
-    hi = __shfl(hi, srcLane);
-    lo = __shfl(lo, srcLane);
+    hi = SHFL(hi, srcLane);
+    lo = SHFL(lo, srcLane);
     // unforunately there isn't an __nv_hiloint2long(hi,lo) intrinsic cast
     int2 fuse; fuse.x = lo; fuse.y = hi;
     return *reinterpret_cast<long long*>(&fuse);

platforms/opencl/include/OpenCLKernels.h

@@ -758,7 +758,7 @@ private:
     std::vector<OpenCLArray*> tabulatedFunctions;
     double longRangeCoefficient;
     std::vector<double> longRangeCoefficientDerivs;
-    bool hasInitializedLongRangeCorrection, hasInitializedKernel;
+    bool hasInitializedLongRangeCorrection, hasInitializedKernel, hasParamDerivs;
     int numGroupThreadBlocks;
     CustomNonbondedForce* forceCopy;
     const System& system;

platforms/opencl/src/OpenCLKernels.cpp

@@ -2789,6 +2789,7 @@ void OpenCLCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNon
     
     // Create the kernel.
     
+    hasParamDerivs = (force.getNumEnergyParameterDerivatives() > 0);
     map<string, string> replacements;
     replacements["COMPUTE_INTERACTION"] = interactionSource;
     const string suffixes[] = {"x", "y", "z", "w"};
@@ -2812,6 +2813,8 @@ void OpenCLCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNon
         args << ", __global const " << tableTypes[i]<< "* restrict table" << i;
     if (globals != NULL)
         args<<", __global const float* restrict globals";
+    if (hasParamDerivs)
+        args << ", __global mixed* restrict energyParamDerivs";
     replacements["PARAMETER_ARGUMENTS"] = args.str();
     stringstream load1;
     for (int i = 0; i < (int) buffers.size(); i++)
@@ -2843,6 +2846,19 @@ void OpenCLCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNon
         }
     }
     replacements["LOAD_ATOM2_PARAMETERS"] = load2.str();
+    stringstream initDerivs, saveDerivs;
+    const vector<string>& allParamDerivNames = cl.getEnergyParamDerivNames();
+    int numDerivs = allParamDerivNames.size();
+    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
+        string paramName = force.getEnergyParameterDerivativeName(i);
+        string derivVariable = cl.getNonbondedUtilities().addEnergyParameterDerivative(paramName);
+        initDerivs<<"mixed "<<derivVariable<<" = 0;\n";
+        for (int index = 0; index < numDerivs; index++)
+            if (allParamDerivNames[index] == paramName)
+                saveDerivs<<"energyParamDerivs[get_global_id(0)*"<<numDerivs<<"+"<<index<<"] += "<<derivVariable<<";\n";
+    }
+    replacements["INIT_DERIVATIVES"] = initDerivs.str();
+    replacements["SAVE_DERIVATIVES"] = saveDerivs.str();
     map<string, string> defines;
     if (force.getNonbondedMethod() != CustomNonbondedForce::NoCutoff)
         defines["USE_CUTOFF"] = "1";
@@ -2902,6 +2918,8 @@ double OpenCLCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool
                 interactionGroupKernel.setArg<cl::Memory>(index++, function->getDeviceBuffer());
             if (globals != NULL)
                 interactionGroupKernel.setArg<cl::Buffer>(index++, globals->getDeviceBuffer());
+            if (hasParamDerivs)
+                interactionGroupKernel.setArg<cl::Memory>(index++, cl.getEnergyParamDerivBuffer().getDeviceBuffer());
         }
         setPeriodicBoxArgs(cl, interactionGroupKernel, 4);
         int forceThreadBlockSize = max(32, cl.getNonbondedUtilities().getForceThreadBlockSize());

platforms/opencl/src/kernels/customNonbondedGroups.cl

@@ -50,6 +50,7 @@ __kernel void computeInteractionGroups(
     const unsigned int tgx = get_local_id(0) & (TILE_SIZE-1); // index within the warp
     const unsigned int tbx = get_local_id(0) - tgx;           // block warpIndex
     mixed energy = 0;
+    INIT_DERIVATIVES
     __local AtomData localData[LOCAL_MEMORY_SIZE];
 
     const unsigned int startTile = FIRST_TILE+warp*(LAST_TILE-FIRST_TILE)/totalWarps;
@@ -93,6 +94,7 @@ __kernel void computeInteractionGroups(
                     LOAD_ATOM2_PARAMETERS
                     real dEdR = 0.0f;
                     real tempEnergy = 0.0f;
+                    const real interactionScale = 1.0f;
                     COMPUTE_INTERACTION
                     energy += tempEnergy;
                     delta *= dEdR;
@@ -125,4 +127,5 @@ __kernel void computeInteractionGroups(
 #endif
     }
     energyBuffer[get_global_id(0)] += energy;
+    SAVE_DERIVATIVES
 }

tests/TestCustomNonbondedForce.h

@@ -1188,6 +1188,47 @@ void testEnergyParameterDerivatives2() {
     ASSERT_EQUAL_TOL((energy1-energy2)/(2*delta), derivs["a"], 1e-4);
 }
 
+void testEnergyParameterDerivativesWithGroups() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomNonbondedForce* nonbonded = new CustomNonbondedForce("k*(r-r0)^2");
+    nonbonded->addGlobalParameter("r0", 0.0);
+    nonbonded->addGlobalParameter("k", 0.0);
+    nonbonded->addEnergyParameterDerivative("k");
+    nonbonded->addEnergyParameterDerivative("r0");
+    vector<double> parameters;
+    nonbonded->addParticle(parameters);
+    nonbonded->addParticle(parameters);
+    nonbonded->addParticle(parameters);
+    set<int> set1, set2;
+    set1.insert(1);
+    set2.insert(0);
+    set2.insert(2);
+    nonbonded->addInteractionGroup(set1, set2);
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 2, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    context.setPositions(positions);
+    for (int i = 0; i < 10; i++) {
+        double r0 = 0.1*i;
+        double k = 10-i;
+        context.setParameter("r0", r0);
+        context.setParameter("k", k);
+        State state = context.getState(State::ParameterDerivatives);
+        map<string, double> derivs = state.getEnergyParameterDerivatives();
+        double dEdr0 = -2*k*((2-r0)+(1-r0));
+        double dEdk = (2-r0)*(2-r0) + (1-r0)*(1-r0);
+        ASSERT_EQUAL_TOL(dEdr0, derivs["r0"], 1e-5);
+        ASSERT_EQUAL_TOL(dEdk, derivs["k"], 1e-5);
+    }
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -1217,6 +1258,7 @@ int main(int argc, char* argv[]) {
         testIllegalVariable();
         testEnergyParameterDerivatives();
         testEnergyParameterDerivatives2();
+        testEnergyParameterDerivativesWithGroups();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestNonbondedForce.h

@@ -261,6 +261,7 @@ void testCutoff14() {
     positions[2] = Vec3(2, 0, 0);
     positions[3] = Vec3(3, 0, 0);
     positions[4] = Vec3(4, 0, 0);
+    context.setPositions(positions);
     for (int i = 1; i < 5; ++i) {
 
         // Test LJ forces
@@ -271,8 +272,7 @@ void testCutoff14() {
         nonbonded->setParticleParameters(i, 0, 1.5, 1);
         nonbonded->setExceptionParameters(first14, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
         nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0.0);
-        context.reinitialize();
-        context.setPositions(positions);
+        context.reinitialize(true);
         State state = context.getState(State::Forces | State::Energy);
         const vector<Vec3>& forces = state.getForces();
         double r = positions[i][0];
@@ -299,8 +299,7 @@ void testCutoff14() {
         nonbonded->setParticleParameters(i, q, 1.5, 0);
         nonbonded->setExceptionParameters(first14, 0, 3, i == 3 ? q*q/1.2 : 0, 1.5, 0);
         nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0);
-        context.reinitialize();
-        context.setPositions(positions);
+        context.reinitialize(true);
         state = context.getState(State::Forces | State::Energy);
         const vector<Vec3>& forces2 = state.getForces();
         force = ONE_4PI_EPS0*q*q/(r*r);

wrappers/python/CMakeLists.txt

@@ -55,6 +55,7 @@ foreach(SUBDIR ${SUBDIRS})
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.xml"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.pdb"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.pdbx"
+        "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.cif"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.prmtop"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.prm"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.inpcrd"

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

@@ -42,7 +42,7 @@ class DCDReporter(object):
     To use it, create a DCDReporter, then add it to the Simulation's list of reporters.
     """
 
-    def __init__(self, file, reportInterval, append=False):
+    def __init__(self, file, reportInterval, append=False, enforcePeriodicBox=None):
         """Create a DCDReporter.
 
         Parameters
@@ -53,9 +53,15 @@ class DCDReporter(object):
             The interval (in time steps) at which to write frames
         append : bool=False
             If True, open an existing DCD file to append to.  If False, create a new file.
+        enforcePeriodicBox: bool
+            Specifies whether particle positions should be translated so the center of every molecule
+            lies in the same periodic box.  If None (the default), it will automatically decide whether
+            to translate molecules based on whether the system being simulated uses periodic boundary
+            conditions.
         """
         self._reportInterval = reportInterval
         self._append = append
+        self._enforcePeriodicBox = enforcePeriodicBox
         if append:
             mode = 'r+b'
         else:
@@ -74,13 +80,14 @@ class DCDReporter(object):
         Returns
         -------
         tuple
-            A five element tuple. The first element is the number of steps
-            until the next report. The remaining elements specify whether
+            A six element tuple. The first element is the number of steps
+            until the next report. The next four elements specify whether
             that report will require positions, velocities, forces, and
-            energies respectively.
+            energies respectively.  The final element specifies whether
+            positions should be wrapped to lie in a single periodic box.
         """
         steps = self._reportInterval - simulation.currentStep%self._reportInterval
-        return (steps, True, False, False, False)
+        return (steps, True, False, False, False, self._enforcePeriodicBox)
 
     def report(self, simulation, state):
         """Generate a report.

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

@@ -426,8 +426,8 @@ class ForceField(object):
     def registerTemplatePatch(self, residue, patch, patchResidueIndex):
         """Register that a particular patch can be used with a particular residue."""
         if residue not in self._templatePatches:
-            self._templatePatches[residue] = []
-        self._templatePatches[residue].append((patch, patchResidueIndex))
+            self._templatePatches[residue] = set()
+        self._templatePatches[residue].add((patch, patchResidueIndex))
 
     def registerScript(self, script):
         """Register a new script to be executed after building the System."""

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

@@ -41,7 +41,7 @@ class PDBReporter(object):
     To use it, create a PDBReporter, then add it to the Simulation's list of reporters.
     """
 
-    def __init__(self, file, reportInterval):
+    def __init__(self, file, reportInterval, enforcePeriodicBox=None):
         """Create a PDBReporter.
 
         Parameters
@@ -50,8 +50,14 @@ class PDBReporter(object):
             The file to write to
         reportInterval : int
             The interval (in time steps) at which to write frames
+        enforcePeriodicBox: bool
+            Specifies whether particle positions should be translated so the center of every molecule
+            lies in the same periodic box.  If None (the default), it will automatically decide whether
+            to translate molecules based on whether the system being simulated uses periodic boundary
+            conditions.
         """
         self._reportInterval = reportInterval
+        self._enforcePeriodicBox = enforcePeriodicBox
         self._out = open(file, 'w')
         self._topology = None
         self._nextModel = 0
@@ -67,13 +73,14 @@ class PDBReporter(object):
         Returns
         -------
         tuple
-            A five element tuple. The first element is the number of steps
-            until the next report. The remaining elements specify whether
+            A six element tuple. The first element is the number of steps
+            until the next report. The next four elements specify whether
             that report will require positions, velocities, forces, and
-            energies respectively.
+            energies respectively.  The final element specifies whether
+            positions should be wrapped to lie in a single periodic box.
         """
         steps = self._reportInterval - simulation.currentStep%self._reportInterval
-        return (steps, True, False, False, False)
+        return (steps, True, False, False, False, self._enforcePeriodicBox)
 
     def report(self, simulation, state):
         """Generate a report.