Merge remote-tracking branch 'upstream/master'

platforms/cuda/src/CudaExpressionUtilities.cpp

@@ -463,6 +463,9 @@ void CudaExpressionUtilities::processExpression(stringstream& out, const Express
         case Operation::CEIL:
             out << "ceil(" << getTempName(node.getChildren()[0], temps) << ")";
             break;
+        case Operation::SELECT:
+            out << "(" << getTempName(node.getChildren()[0], temps) << " != 0 ? " << getTempName(node.getChildren()[1], temps) << " : " << getTempName(node.getChildren()[2], temps) << ")";
+            break;
         default:
             throw OpenMMException("Internal error: Unknown operation in user-defined expression: "+node.getOperation().getName());
     }

platforms/cuda/tests/TestCudaCustomIntegrator.cpp

@@ -334,7 +334,7 @@ void testMonteCarlo() {
     integrator.addComputePerDof("oldx", "x");
     integrator.addComputePerDof("x", "x+dt*gaussian");
     integrator.addComputeGlobal("accept", "step(exp((oldE-energy)/kT)-uniform)");
-    integrator.addComputePerDof("x", "accept*x + (1-accept)*oldx");
+    integrator.addComputePerDof("x", "select(accept, x, oldx)");
     HarmonicBondForce* forceField = new HarmonicBondForce();
     forceField->addBond(0, 1, 2.0, 10.0);
     system.addForce(forceField);

platforms/cuda/tests/TestCudaLocalEnergyMinimizer.cpp

@@ -7,7 +7,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2010-2014 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -81,7 +81,7 @@ void testLargeSystem() {
     const int numParticles = numMolecules*2;
     const double cutoff = 2.0;
     const double boxSize = 4.0;
-    const double tolerance = 5;
+    const double tolerance = 10;
     System system;
     system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     NonbondedForce* nonbonded = new NonbondedForce();
@@ -114,7 +114,7 @@ void testLargeSystem() {
     State finalState = context.getState(State::Forces | State::Energy | State::Positions);
     ASSERT(finalState.getPotentialEnergy() < initialState.getPotentialEnergy());
 
-    // Compute the force magnitude, substracting off any component parallel to a constraint, and
+    // Compute the force magnitude, subtracting off any component parallel to a constraint, and
     // check that it satisfies the requested tolerance.
 
     double forceNorm = 0.0;
@@ -129,8 +129,8 @@ void testLargeSystem() {
         f -= dir*dir.dot(f);
         forceNorm += f.dot(f);
     }
-    forceNorm = sqrt(forceNorm/(4*numMolecules));
-    ASSERT(forceNorm < 3*tolerance);
+    forceNorm = sqrt(forceNorm/(5*numMolecules));
+    ASSERT(forceNorm < 2*tolerance);
 }
 
 void testVirtualSites() {
@@ -138,7 +138,7 @@ void testVirtualSites() {
     const int numParticles = numMolecules*3;
     const double cutoff = 2.0;
     const double boxSize = 4.0;
-    const double tolerance = 5;
+    const double tolerance = 10;
     System system;
     system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     NonbondedForce* nonbonded = new NonbondedForce();
@@ -195,8 +195,8 @@ void testVirtualSites() {
         
         ASSERT_EQUAL_VEC((finalState.getPositions()[i+1]+finalState.getPositions()[i])*0.5, finalState.getPositions()[i+2], 1e-5);
     }
-    forceNorm = sqrt(forceNorm/(4*numMolecules));
-    ASSERT(forceNorm < 3*tolerance);
+    forceNorm = sqrt(forceNorm/(5*numMolecules));
+    ASSERT(forceNorm < 2*tolerance);
 }
 
 int main(int argc, char* argv[]) {

platforms/opencl/src/OpenCLExpressionUtilities.cpp

@@ -455,6 +455,9 @@ void OpenCLExpressionUtilities::processExpression(stringstream& out, const Expre
         case Operation::CEIL:
             out << "ceil(" << getTempName(node.getChildren()[0], temps) << ")";
             break;
+        case Operation::SELECT:
+            out << "(" << getTempName(node.getChildren()[0], temps) << " != 0 ? " << getTempName(node.getChildren()[1], temps) << " : " << getTempName(node.getChildren()[2], temps) << ")";
+            break;
         default:
             throw OpenMMException("Internal error: Unknown operation in user-defined expression: "+node.getOperation().getName());
     }

platforms/opencl/src/OpenCLNonbondedUtilities.cpp

@@ -303,6 +303,7 @@ void OpenCLNonbondedUtilities::initialize(const System& system) {
         defines["PADDED_CUTOFF_SQUARED"] = context.doubleToString(paddedCutoff*paddedCutoff);
         defines["NUM_TILES_WITH_EXCLUSIONS"] = context.intToString(exclusionTiles->getSize());
         defines["NUM_BLOCKS"] = context.intToString(context.getNumAtomBlocks());
+        defines["SIMD_WIDTH"] = context.intToString(context.getSIMDWidth());
         if (usePeriodic)
             defines["USE_PERIODIC"] = "1";
         int maxExclusions = 0;

platforms/opencl/src/kernels/findInteractingBlocks.cl

 #pragma OPENCL EXTENSION cl_khr_global_int32_base_atomics : enable
 #pragma OPENCL EXTENSION cl_khr_byte_addressable_store : enable
-#define BUFFER_SIZE 256
 
 /**
  * Find a bounding box for the atoms in each block.
@@ -65,6 +64,10 @@ __kernel void sortBoxData(__global const real2* restrict sortedBlock, __global c
     }
 }
 
+#if SIMD_WIDTH <= 32
+
+#define BUFFER_SIZE 256
+
 __kernel void findBlocksWithInteractions(real4 periodicBoxSize, real4 invPeriodicBoxSize, real4 periodicBoxVecX, real4 periodicBoxVecY, real4 periodicBoxVecZ,
         __global unsigned int* restrict interactionCount, __global int* restrict interactingTiles, __global unsigned int* restrict interactingAtoms,
         __global const real4* restrict posq, unsigned int maxTiles, unsigned int startBlockIndex, unsigned int numBlocks, __global real2* restrict sortedBlocks,
@@ -133,8 +136,8 @@ __kernel void findBlocksWithInteractions(real4 periodicBoxSize, real4 invPeriodi
             int block2 = block2Base+indexInWarp;
             bool includeBlock2 = (block2 < NUM_BLOCKS);
             if (includeBlock2) {
-                real4 blockCenterY = (block2 < NUM_BLOCKS ? sortedBlockCenter[block2] : (real4) (0));
-                real4 blockSizeY = (block2 < NUM_BLOCKS ? sortedBlockBoundingBox[block2] : (real4) (0));
+                real4 blockCenterY = sortedBlockCenter[block2];
+                real4 blockSizeY = sortedBlockBoundingBox[block2];
                 real4 blockDelta = blockCenterX-blockCenterY;
 #ifdef USE_PERIODIC
                 APPLY_PERIODIC_TO_DELTA(blockDelta)
@@ -252,3 +255,279 @@ __kernel void findBlocksWithInteractions(real4 periodicBoxSize, real4 invPeriodi
     for (int i = get_global_id(0); i < NUM_ATOMS; i += get_global_size(0))
         oldPositions[i] = posq[i];
 }
+
+#else
+// This is the old implementation of finding interacting blocks.  It is quite a bit more complicated,
+// and slower on most GPUs.  On AMD, however, it is faster, so we keep it around to use there.
+
+#define BUFFER_SIZE BUFFER_GROUPS*GROUP_SIZE
+#define WARP_SIZE 32
+#define INVALID 0xFFFF
+
+/**
+ * Perform a parallel prefix sum over an array.  The input values are all assumed to be 0 or 1.
+ */
+void prefixSum(__local short* sum, __local ushort2* temp) {
+    for (int i = get_local_id(0); i < BUFFER_SIZE; i += get_local_size(0))
+        temp[i].x = sum[i];
+    barrier(CLK_LOCAL_MEM_FENCE);
+    int whichBuffer = 0;
+    for (int offset = 1; offset < BUFFER_SIZE; offset *= 2) {
+        if (whichBuffer == 0)
+            for (int i = get_local_id(0); i < BUFFER_SIZE; i += get_local_size(0))
+                temp[i].y = (i < offset ? temp[i].x : temp[i].x+temp[i-offset].x);
+        else
+            for (int i = get_local_id(0); i < BUFFER_SIZE; i += get_local_size(0))
+                temp[i].x = (i < offset ? temp[i].y : temp[i].y+temp[i-offset].y);
+        whichBuffer = 1-whichBuffer;
+        barrier(CLK_LOCAL_MEM_FENCE);
+    }
+    if (whichBuffer == 0)
+        for (int i = get_local_id(0); i < BUFFER_SIZE; i += get_local_size(0))
+            sum[i] = temp[i].x;
+    else
+        for (int i = get_local_id(0); i < BUFFER_SIZE; i += get_local_size(0))
+            sum[i] = temp[i].y;
+    barrier(CLK_LOCAL_MEM_FENCE);
+}
+
+/**
+ * This is called by findBlocksWithInteractions().  It compacts the list of blocks, identifies interactions
+ * in them, and writes the result to global memory.
+ */
+void storeInteractionData(int x, __local unsigned short* buffer, __local short* sum, __local ushort2* temp, __local int* atoms, __local int* numAtoms,
+            __local int* baseIndex, __global unsigned int* interactionCount, __global int* interactingTiles, __global unsigned int* interactingAtoms, real4 periodicBoxSize,
+            real4 invPeriodicBoxSize, real4 periodicBoxVecX, real4 periodicBoxVecY, real4 periodicBoxVecZ, __global const real4* posq, __local real4* posBuffer,
+            real4 blockCenterX, real4 blockSizeX, unsigned int maxTiles, bool finish) {
+    const bool singlePeriodicCopy = (0.5f*periodicBoxSize.x-blockSizeX.x >= PADDED_CUTOFF &&
+                                     0.5f*periodicBoxSize.y-blockSizeX.y >= PADDED_CUTOFF &&
+                                     0.5f*periodicBoxSize.z-blockSizeX.z >= PADDED_CUTOFF);
+    if (get_local_id(0) < TILE_SIZE) {
+        real4 pos = posq[x*TILE_SIZE+get_local_id(0)];
+#ifdef USE_PERIODIC
+        if (singlePeriodicCopy) {
+            // The box is small enough that we can just translate all the atoms into a single periodic
+            // box, then skip having to apply periodic boundary conditions later.
+            
+            APPLY_PERIODIC_TO_POS_WITH_CENTER(pos, blockCenterX)
+        }
+#endif
+        posBuffer[get_local_id(0)] = pos;
+    }
+    
+    // The buffer is full, so we need to compact it and write out results.  Start by doing a parallel prefix sum.
+
+    barrier(CLK_LOCAL_MEM_FENCE);
+    for (int i = get_local_id(0); i < BUFFER_SIZE; i += get_local_size(0))
+        sum[i] = (buffer[i] == INVALID ? 0 : 1);
+    barrier(CLK_LOCAL_MEM_FENCE);
+    prefixSum(sum, temp);
+    int numValid = sum[BUFFER_SIZE-1];
+
+    // Compact the buffer.
+
+    for (int i = get_local_id(0); i < BUFFER_SIZE; i += get_local_size(0))
+        if (buffer[i] != INVALID)
+            temp[sum[i]-1].x = buffer[i];
+    barrier(CLK_LOCAL_MEM_FENCE);
+    for (int i = get_local_id(0); i < BUFFER_SIZE; i += get_local_size(0))
+        buffer[i] = temp[i].x;
+    barrier(CLK_LOCAL_MEM_FENCE);
+
+    // Loop over the tiles and find specific interactions in them.
+
+    const int indexInWarp = get_local_id(0)%WARP_SIZE;
+    for (int base = 0; base < numValid; base += BUFFER_SIZE/WARP_SIZE) {
+        for (int i = get_local_id(0)/WARP_SIZE; i < BUFFER_SIZE/WARP_SIZE && base+i < numValid; i += GROUP_SIZE/WARP_SIZE) {
+            // Check each atom in block Y for interactions.
+            
+            real4 pos = posq[buffer[base+i]*TILE_SIZE+indexInWarp];
+#ifdef USE_PERIODIC
+            if (singlePeriodicCopy)
+                APPLY_PERIODIC_TO_POS_WITH_CENTER(pos, blockCenterX)
+#endif
+            bool interacts = false;
+#ifdef USE_PERIODIC
+            if (!singlePeriodicCopy) {
+                for (int j = 0; j < TILE_SIZE; j++) {
+                    real4 delta = pos-posBuffer[j];
+                    APPLY_PERIODIC_TO_DELTA(delta)
+                    interacts |= (delta.x*delta.x+delta.y*delta.y+delta.z*delta.z < PADDED_CUTOFF_SQUARED);
+                }
+            }
+            else {
+#endif
+                for (int j = 0; j < TILE_SIZE; j++) {
+                    real4 delta = pos-posBuffer[j];
+                    interacts |= (delta.x*delta.x+delta.y*delta.y+delta.z*delta.z < PADDED_CUTOFF_SQUARED);
+                }
+#ifdef USE_PERIODIC
+            }
+#endif
+            sum[i*WARP_SIZE+indexInWarp] = (interacts ? 1 : 0);
+        }
+        for (int i = numValid-base+get_local_id(0)/WARP_SIZE; i < BUFFER_SIZE/WARP_SIZE; i += GROUP_SIZE/WARP_SIZE)
+            sum[i*WARP_SIZE+indexInWarp] = 0;
+
+        // Compact the list of atoms.
+
+        barrier(CLK_LOCAL_MEM_FENCE);
+        prefixSum(sum, temp);
+        for (int i = get_local_id(0); i < BUFFER_SIZE; i += get_local_size(0))
+            if (sum[i] != (i == 0 ? 0 : sum[i-1]))
+                atoms[*numAtoms+sum[i]-1] = buffer[base+i/WARP_SIZE]*TILE_SIZE+indexInWarp;
+
+        // Store them to global memory.
+
+        int atomsToStore = *numAtoms+sum[BUFFER_SIZE-1];
+        bool storePartialTile = (finish && base >= numValid-BUFFER_SIZE/WARP_SIZE);
+        int tilesToStore = (storePartialTile ? (atomsToStore+TILE_SIZE-1)/TILE_SIZE : atomsToStore/TILE_SIZE);
+        if (tilesToStore > 0) {
+            if (get_local_id(0) == 0)
+                *baseIndex = atom_add(interactionCount, tilesToStore);
+            barrier(CLK_LOCAL_MEM_FENCE);
+            if (get_local_id(0) == 0)
+                *numAtoms = atomsToStore-tilesToStore*TILE_SIZE;
+            if (*baseIndex+tilesToStore <= maxTiles) {
+                if (get_local_id(0) < tilesToStore)
+                    interactingTiles[*baseIndex+get_local_id(0)] = x;
+                for (int i = get_local_id(0); i < tilesToStore*TILE_SIZE; i += get_local_size(0))
+                    interactingAtoms[*baseIndex*TILE_SIZE+i] = (i < atomsToStore ? atoms[i] : NUM_ATOMS);
+            }
+        }
+        else {
+            barrier(CLK_LOCAL_MEM_FENCE);
+            if (get_local_id(0) == 0)
+                *numAtoms += sum[BUFFER_SIZE-1];
+        }
+        barrier(CLK_LOCAL_MEM_FENCE);
+        if (get_local_id(0) < *numAtoms && !storePartialTile)
+            atoms[get_local_id(0)] = atoms[tilesToStore*TILE_SIZE+get_local_id(0)];
+    }
+
+    if (numValid == 0 && *numAtoms > 0 && finish) {
+        // We didn't have any more tiles to process, but there were some atoms left over from a
+        // previous call to this function.  Save them now.
+
+        if (get_local_id(0) == 0)
+            *baseIndex = atom_add(interactionCount, 1);
+        barrier(CLK_LOCAL_MEM_FENCE);
+        if (*baseIndex < maxTiles) {
+            if (get_local_id(0) == 0)
+                interactingTiles[*baseIndex] = x;
+            if (get_local_id(0) < TILE_SIZE)
+                interactingAtoms[*baseIndex*TILE_SIZE+get_local_id(0)] = (get_local_id(0) < *numAtoms ? atoms[get_local_id(0)] : NUM_ATOMS);
+        }
+    }
+
+    // Reset the buffer for processing more tiles.
+
+    for (int i = get_local_id(0); i < BUFFER_SIZE; i += get_local_size(0))
+        buffer[i] = INVALID;
+    barrier(CLK_LOCAL_MEM_FENCE);
+}
+
+/**
+ * Compare the bounding boxes for each pair of blocks.  If they are sufficiently far apart,
+ * mark them as non-interacting.
+ */
+__kernel void findBlocksWithInteractions(real4 periodicBoxSize, real4 invPeriodicBoxSize, real4 periodicBoxVecX, real4 periodicBoxVecY, real4 periodicBoxVecZ,
+        __global unsigned int* restrict interactionCount, __global int* restrict interactingTiles, __global unsigned int* restrict interactingAtoms,
+        __global const real4* restrict posq, unsigned int maxTiles, unsigned int startBlockIndex, unsigned int numBlocks, __global real2* restrict sortedBlocks,
+        __global const real4* restrict sortedBlockCenter, __global const real4* restrict sortedBlockBoundingBox,
+        __global const unsigned int* restrict exclusionIndices, __global const unsigned int* restrict exclusionRowIndices, __global real4* restrict oldPositions,
+        __global const int* restrict rebuildNeighborList) {
+    __local unsigned short buffer[BUFFER_SIZE];
+    __local short sum[BUFFER_SIZE];
+    __local ushort2 temp[BUFFER_SIZE];
+    __local int atoms[BUFFER_SIZE+TILE_SIZE];
+    __local real4 posBuffer[TILE_SIZE];
+    __local int exclusionsForX[MAX_EXCLUSIONS];
+    __local int bufferFull;
+    __local int globalIndex;
+    __local int numAtoms;
+#ifdef AMD_ATOMIC_WORK_AROUND
+    // Do a byte write to force all memory accesses to interactionCount to use the complete path.
+    // This avoids the atomic access from causing all word accesses to other buffers from using the slow complete path.
+    // The IF actually causes the write to never be executed, its presence is all that is needed.
+    // AMD APP SDK 2.4 has this problem.
+    if (get_global_id(0) == get_local_id(0)+1)
+        ((__global char*)interactionCount)[sizeof(unsigned int)+1] = 0;
+#endif
+
+    if (rebuildNeighborList[0] == 0)
+        return; // The neighbor list doesn't need to be rebuilt.
+
+    int valuesInBuffer = 0;
+    if (get_local_id(0) == 0)
+        bufferFull = false;
+    for (int i = 0; i < BUFFER_GROUPS; ++i)
+        buffer[i*GROUP_SIZE+get_local_id(0)] = INVALID;
+    barrier(CLK_LOCAL_MEM_FENCE);
+    
+    // Loop over blocks sorted by size.
+    
+    for (int i = startBlockIndex+get_group_id(0); i < startBlockIndex+numBlocks; i += get_num_groups(0)) {
+        if (get_local_id(0) == get_local_size(0)-1)
+            numAtoms = 0;
+        real2 sortedKey = sortedBlocks[i];
+        int x = (int) sortedKey.y;
+        real4 blockCenterX = sortedBlockCenter[i];
+        real4 blockSizeX = sortedBlockBoundingBox[i];
+
+        // Load exclusion data for block x.
+        
+        const int exclusionStart = exclusionRowIndices[x];
+        const int exclusionEnd = exclusionRowIndices[x+1];
+        const int numExclusions = exclusionEnd-exclusionStart;
+        for (int j = get_local_id(0); j < numExclusions; j += get_local_size(0))
+            exclusionsForX[j] = exclusionIndices[exclusionStart+j];
+        barrier(CLK_LOCAL_MEM_FENCE);
+        
+        // Compare it to other blocks after this one in sorted order.
+        
+        for (int base = i+1; base < NUM_BLOCKS; base += get_local_size(0)) {
+            int j = base+get_local_id(0);
+            real2 sortedKey2 = (j < NUM_BLOCKS ? sortedBlocks[j] : (real2) 0);
+            real4 blockCenterY = (j < NUM_BLOCKS ? sortedBlockCenter[j] : (real4) 0);
+            real4 blockSizeY = (j < NUM_BLOCKS ? sortedBlockBoundingBox[j] : (real4) 0);
+            unsigned short y = (unsigned short) sortedKey2.y;
+            real4 delta = blockCenterX-blockCenterY;
+#ifdef USE_PERIODIC
+            APPLY_PERIODIC_TO_DELTA(delta)
+#endif
+            delta.x = max((real) 0, fabs(delta.x)-blockSizeX.x-blockSizeY.x);
+            delta.y = max((real) 0, fabs(delta.y)-blockSizeX.y-blockSizeY.y);
+            delta.z = max((real) 0, fabs(delta.z)-blockSizeX.z-blockSizeY.z);
+            bool hasExclusions = false;
+            for (int k = 0; k < numExclusions; k++)
+                hasExclusions |= (exclusionsForX[k] == y);
+            if (j < NUM_BLOCKS && delta.x*delta.x+delta.y*delta.y+delta.z*delta.z < PADDED_CUTOFF_SQUARED && !hasExclusions) {
+                // Add this tile to the buffer.
+
+                int bufferIndex = valuesInBuffer*GROUP_SIZE+get_local_id(0);
+                buffer[bufferIndex] = y;
+                valuesInBuffer++;
+                if (!bufferFull && valuesInBuffer == BUFFER_GROUPS)
+                    bufferFull = true;
+            }
+            barrier(CLK_LOCAL_MEM_FENCE);
+            if (bufferFull) {
+                storeInteractionData(x, buffer, sum, temp, atoms, &numAtoms, &globalIndex, interactionCount, interactingTiles, interactingAtoms, periodicBoxSize, invPeriodicBoxSize, periodicBoxVecX, periodicBoxVecY, periodicBoxVecZ, posq, posBuffer, blockCenterX, blockSizeX, maxTiles, false);
+                valuesInBuffer = 0;
+                if (get_local_id(0) == 0)
+                    bufferFull = false;
+            }
+            barrier(CLK_LOCAL_MEM_FENCE);
+        }
+        storeInteractionData(x, buffer, sum, temp, atoms, &numAtoms, &globalIndex, interactionCount, interactingTiles, interactingAtoms, periodicBoxSize, invPeriodicBoxSize, periodicBoxVecX, periodicBoxVecY, periodicBoxVecZ, posq, posBuffer, blockCenterX, blockSizeX, maxTiles, true);
+    }
+    
+    // Record the positions the neighbor list is based on.
+    
+    for (int i = get_global_id(0); i < NUM_ATOMS; i += get_global_size(0))
+        oldPositions[i] = posq[i];
+}
+
+#endif
+

platforms/opencl/tests/TestOpenCLCustomIntegrator.cpp

@@ -334,7 +334,7 @@ void testMonteCarlo() {
     integrator.addComputePerDof("oldx", "x");
     integrator.addComputePerDof("x", "x+dt*gaussian");
     integrator.addComputeGlobal("accept", "step(exp((oldE-energy)/kT)-uniform)");
-    integrator.addComputePerDof("x", "accept*x + (1-accept)*oldx");
+    integrator.addComputePerDof("x", "select(accept, x, oldx)");
     HarmonicBondForce* forceField = new HarmonicBondForce();
     forceField->addBond(0, 1, 2.0, 10.0);
     system.addForce(forceField);

platforms/opencl/tests/TestOpenCLLocalEnergyMinimizer.cpp

@@ -7,7 +7,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2010-2014 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -81,7 +81,7 @@ void testLargeSystem() {
     const int numParticles = numMolecules*2;
     const double cutoff = 2.0;
     const double boxSize = 4.0;
-    const double tolerance = 5;
+    const double tolerance = 10;
     System system;
     system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     NonbondedForce* nonbonded = new NonbondedForce();
@@ -114,7 +114,7 @@ void testLargeSystem() {
     State finalState = context.getState(State::Forces | State::Energy | State::Positions);
     ASSERT(finalState.getPotentialEnergy() < initialState.getPotentialEnergy());
 
-    // Compute the force magnitude, substracting off any component parallel to a constraint, and
+    // Compute the force magnitude, subtracting off any component parallel to a constraint, and
     // check that it satisfies the requested tolerance.
 
     double forceNorm = 0.0;
@@ -129,8 +129,8 @@ void testLargeSystem() {
         f -= dir*dir.dot(f);
         forceNorm += f.dot(f);
     }
-    forceNorm = sqrt(forceNorm/(4*numMolecules));
-    ASSERT(forceNorm < 3*tolerance);
+    forceNorm = sqrt(forceNorm/(5*numMolecules));
+    ASSERT(forceNorm < 2*tolerance);
 }
 
 void testVirtualSites() {
@@ -138,7 +138,7 @@ void testVirtualSites() {
     const int numParticles = numMolecules*3;
     const double cutoff = 2.0;
     const double boxSize = 4.0;
-    const double tolerance = 5;
+    const double tolerance = 10;
     System system;
     system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     NonbondedForce* nonbonded = new NonbondedForce();
@@ -195,8 +195,8 @@ void testVirtualSites() {
         
         ASSERT_EQUAL_VEC((finalState.getPositions()[i+1]+finalState.getPositions()[i])*0.5, finalState.getPositions()[i+2], 1e-5);
     }
-    forceNorm = sqrt(forceNorm/(4*numMolecules));
-    ASSERT(forceNorm < 3*tolerance);
+    forceNorm = sqrt(forceNorm/(5*numMolecules));
+    ASSERT(forceNorm < 2*tolerance);
 }
 
 int main(int argc, char* argv[]) {

platforms/reference/tests/TestReferenceCustomIntegrator.cpp

@@ -334,7 +334,7 @@ void testMonteCarlo() {
     integrator.addComputePerDof("oldx", "x");
     integrator.addComputePerDof("x", "x+dt*gaussian");
     integrator.addComputeGlobal("accept", "step(exp((oldE-energy)/kT)-uniform)");
-    integrator.addComputePerDof("x", "accept*x + (1-accept)*oldx");
+    integrator.addComputePerDof("x", "select(accept, x, oldx)");
     HarmonicBondForce* forceField = new HarmonicBondForce();
     forceField->addBond(0, 1, 2.0, 10.0);
     system.addForce(forceField);

platforms/reference/tests/TestReferenceLocalEnergyMinimizer.cpp

@@ -7,7 +7,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2010-2014 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -129,8 +129,8 @@ void testLargeSystem() {
         f -= dir*dir.dot(f);
         forceNorm += f.dot(f);
     }
-    forceNorm = sqrt(forceNorm/(4*numMolecules));
-    ASSERT(forceNorm < 3*tolerance);
+    forceNorm = sqrt(forceNorm/(5*numMolecules));
+    ASSERT(forceNorm < 2*tolerance);
 }
 
 void testVirtualSites() {
@@ -195,8 +195,8 @@ void testVirtualSites() {
         
         ASSERT_EQUAL_VEC((finalState.getPositions()[i+1]+finalState.getPositions()[i])*0.5, finalState.getPositions()[i+2], 1e-5);
     }
-    forceNorm = sqrt(forceNorm/(4*numMolecules));
-    ASSERT(forceNorm < 3*tolerance);
+    forceNorm = sqrt(forceNorm/(5*numMolecules));
+    ASSERT(forceNorm < 2*tolerance);
 }
 
 int main() {

plugins/cpupme/src/CpuPmeKernels.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2013-2014 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -38,6 +38,7 @@
 #include "openmm/internal/vectorize.h"
 #include <cmath>
 #include <cstring>
+#include <sstream>
 
 using namespace OpenMM;
 using namespace std;
@@ -353,6 +354,9 @@ static void* threadBody(void* args) {
 void CpuCalcPmeReciprocalForceKernel::initialize(int xsize, int ysize, int zsize, int numParticles, double alpha) {
     if (!hasInitializedThreads) {
         numThreads = getNumProcessors();
+        char* threadsEnv = getenv("OPENMM_CPU_THREADS");
+        if (threadsEnv != NULL)
+            stringstream(threadsEnv) >> numThreads;
         fftwf_init_threads();
         hasInitializedThreads = true;
     }

tests/TestParser.cpp

@@ -178,6 +178,8 @@ void verifyDerivative(const string& expression, const string& expectedDeriv) {
     verifySameValue(computed, expected, 2.0, 3.0);
     verifySameValue(computed, expected, -2.0, 3.0);
     verifySameValue(computed, expected, 2.0, -3.0);
+    verifySameValue(computed, expected, 0.0, -3.0);
+    verifySameValue(computed, expected, 2.0, 0.0);
 }
 
 /**
@@ -249,6 +251,8 @@ int main() {
         verifyEvaluation("step(x-3)+y*step(x)", 2.0, 3.0, 3.0);
         verifyEvaluation("floor(x)", -2.1, 3.0, -3.0);
         verifyEvaluation("ceil(x)", -2.1, 3.0, -2.0);
+        verifyEvaluation("select(x, 1.0, y)", 0.3, 2.0, 1.0);
+        verifyEvaluation("select(x, 1.0, y)", 0.0, 2.0, 2.0);
         verifyInvalidExpression("1..2");
         verifyInvalidExpression("1*(2+3");
         verifyInvalidExpression("5++4");
@@ -282,6 +286,7 @@ int main() {
         verifyDerivative("max(5, x^2)", "(1-step(5-x^2))*2*x");
         verifyDerivative("abs(3*x)", "step(3*x)*3+(1-step(3*x))*-3");
         verifyDerivative("floor(x)+0.5*x*ceil(x)", "0.5*ceil(x)");
+        verifyDerivative("select(x, x^2, 3*x)", "select(x, 2*x, 3)");
         testCustomFunction("custom(x, y)/2", "x*y");
         testCustomFunction("custom(x^2, 1)+custom(2, y-1)", "2*x^2+4*(y-1)");
         cout << Parser::parse("x*x").optimize() << endl;

wrappers/python/setup.py

@@ -85,7 +85,7 @@ version = '%(version)s'
 full_version = '%(full_version)s'
 git_revision = '%(git_revision)s'
 release = %(isrelease)s
-openmm_library_path = '%(path)s'
+openmm_library_path = r'%(path)s'
 
 if not release:
     version = full_version

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

@@ -127,6 +127,10 @@ class PrmtopLoader(object):
                         tag, self._prmtopVersion = line.rstrip().split(None, 1)
                     elif line.startswith('%FLAG'):
                         tag, flag = line.rstrip().split(None, 1)
+                        if flag == 'CTITLE':
+                            raise TypeError('CHAMBER-style topology files are not supported here. '
+                                            'Consider using the CHARMM files directly with CharmmPsfFile '
+                                            'or ParmEd (where CHAMBER topologies are supported)')
                         self._flags.append(flag)
                         self._raw_data[flag] = []
                     elif line.startswith('%FORMAT'):

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

@@ -240,7 +240,7 @@ class Modeller(object):
         self.topology = newTopology
         self.positions = newPositions
 
-    def addSolvent(self, forcefield, model='tip3p', boxSize=None, boxVectors=None, padding=None, positiveIon='Na+', negativeIon='Cl-', ionicStrength=0*molar):
+    def addSolvent(self, forcefield, model='tip3p', boxSize=None, boxVectors=None, padding=None, numAdded=None, positiveIon='Na+', negativeIon='Cl-', ionicStrength=0*molar):
         """Add solvent (both water and ions) to the model to fill a rectangular box.
 
         The algorithm works as follows:
@@ -250,11 +250,15 @@ class Modeller(object):
            randomly selecting a water molecule and replacing it with the ion.
         4. Ion pairs are added to give the requested total ionic strength.
 
-        The box size can be specified in four ways.  First, you can explicitly give the vectors defining the periodic box to
-        use.  Alternatively, for a rectangular box you can simply give the dimensions of the unit cell.  Third, you can
-        give a padding distance.  The largest dimension of the solute (along the x, y, or z axis) is determined, and a cubic
-        box of size (largest dimension)+2*padding is used.  Finally, if neither box vectors, box size, nor padding distance is specified,
-        the existing Topology's box vectors are used.
+        The box size can be specified in any of several ways:
+        
+        1. You can explicitly give the vectors defining the periodic box to use.
+        2. Alternatively, for a rectangular box you can simply give the dimensions of the unit cell.
+        3. You can give a padding distance.  The largest dimension of the solute (along the x, y, or z axis) is determined, and a cubic
+        box of size (largest dimension)+2*padding is used.
+        4. You can specify the total number of molecules (both waters and ions) to add.  A cubic box is then created whose size is
+        just large enough hold the specified amount of solvent.
+        5. Finally, if none of the above options is specified, the existing Topology's box vectors are used.
 
         Parameters:
          - forcefield (ForceField) the ForceField to use for determining van der Waals radii and atomic charges
@@ -262,14 +266,43 @@ class Modeller(object):
          - boxSize (Vec3=None) the size of the box to fill with water
          - boxVectors (tuple of Vec3=None) the vectors defining the periodic box to fill with water
          - padding (distance=None) the padding distance to use
+         - numAdded (int=None) the total number of molecules (waters and ions) to add
          - positiveIon (string='Na+') the type of positive ion to add.  Allowed values are 'Cs+', 'K+', 'Li+', 'Na+', and 'Rb+'
          - negativeIon (string='Cl-') the type of negative ion to add.  Allowed values are 'Cl-', 'Br-', 'F-', and 'I-'. Be aware
            that not all force fields support all ion types.
          - ionicStrength (concentration=0*molar) the total concentration of ions (both positive and negative) to add.  This
            does not include ions that are added to neutralize the system.
         """
+        if len([x for x in (boxSize, boxVectors, padding, numAdded) if x is not None]) > 1:
+            raise ValueError('At most one of the following arguments may be specified: boxSize, boxVectors, padding, numAdded')
+
+        # Load the pre-equilibrated water box.
+
+        vdwRadiusPerSigma = 0.5612310241546864907
+        if model == 'tip3p':
+            waterRadius = 0.31507524065751241*vdwRadiusPerSigma
+        elif model == 'spce':
+            waterRadius = 0.31657195050398818*vdwRadiusPerSigma
+        elif model == 'tip4pew':
+            waterRadius = 0.315365*vdwRadiusPerSigma
+        elif model == 'tip5p':
+            waterRadius = 0.312*vdwRadiusPerSigma
+        else:
+            raise ValueError('Unknown water model: %s' % model)
+        pdb = PDBFile(os.path.join(os.path.dirname(__file__), 'data', model+'.pdb'))
+        pdbTopology = pdb.getTopology()
+        pdbPositions = pdb.getPositions().value_in_unit(nanometer)
+        pdbResidues = list(pdbTopology.residues())
+        pdbBoxSize = pdbTopology.getUnitCellDimensions().value_in_unit(nanometer)            
+        
         # Pick a unit cell size.
 
+        if numAdded is not None:
+            # Select a padding distance which is guaranteed to give more than the specified number of molecules.
+            
+            padding = 1.1*(numAdded/((len(pdbResidues)/pdbBoxSize[0]**3)*8))**(1.0/3.0)
+            if padding < 0.5:
+                padding = 0.5 # Ensure we have enough when adding very small numbers of molecules
         if boxVectors is not None:
             if is_quantity(boxVectors[0]):
                 boxVectors = (boxVectors[0].value_in_unit(nanometer), boxVectors[1].value_in_unit(nanometer), boxVectors[2].value_in_unit(nanometer))
@@ -305,25 +338,6 @@ class Modeller(object):
         positiveElement = posIonElements[positiveIon]
         negativeElement = negIonElements[negativeIon]
 
-        # Load the pre-equilibrated water box.
-
-        vdwRadiusPerSigma = 0.5612310241546864907
-        if model == 'tip3p':
-            waterRadius = 0.31507524065751241*vdwRadiusPerSigma
-        elif model == 'spce':
-            waterRadius = 0.31657195050398818*vdwRadiusPerSigma
-        elif model == 'tip4pew':
-            waterRadius = 0.315365*vdwRadiusPerSigma
-        elif model == 'tip5p':
-            waterRadius = 0.312*vdwRadiusPerSigma
-        else:
-            raise ValueError('Unknown water model: %s' % model)
-        pdb = PDBFile(os.path.join(os.path.dirname(__file__), 'data', model+'.pdb'))
-        pdbTopology = pdb.getTopology()
-        pdbPositions = pdb.getPositions().value_in_unit(nanometer)
-        pdbResidues = list(pdbTopology.residues())
-        pdbBoxSize = pdbTopology.getUnitCellDimensions().value_in_unit(nanometer)
-
         # Have the ForceField build a System for the solute from which we can determine van der Waals radii.
 
         system = forcefield.createSystem(self.topology)
@@ -424,6 +438,21 @@ class Modeller(object):
 
                                 addedWaters.append((residue.index, atomPos))
 
+        if numAdded is not None:
+            # We added many more waters than we actually want.  Sort them based on distance to the nearest box edge and
+            # only keep the ones in the middle.
+            
+            lowerBound = center-box/2
+            upperBound = center+box/2
+            distToEdge = (min(min(pos-lowerBound), min(upperBound-pos)) for index, pos in addedWaters)
+            sortedIndex = [i[0] for i in sorted(enumerate(distToEdge), key=lambda x: -x[1])]
+            addedWaters = [addedWaters[i] for i in sortedIndex[:numAdded]]
+            
+            # Compute a new periodic box size.
+            
+            maxSize = max(max((pos[i] for index, pos in addedWaters))-min((pos[i] for index, pos in addedWaters)) for i in range(3))
+            newTopology.setUnitCellDimensions(Vec3(maxSize, maxSize, maxSize))
+        else:
             # There could be clashes between water molecules at the box edges.  Find ones to remove.
     
             upperCutoff = center+box/2-Vec3(waterCutoff, waterCutoff, waterCutoff)

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

@@ -78,11 +78,11 @@ class Simulation(object):
         else:
             self.context = mm.Context(system, integrator, platform, platformProperties)
 
-    def minimizeEnergy(self, tolerance=1*unit.kilojoule/unit.mole, maxIterations=0):
+    def minimizeEnergy(self, tolerance=10*unit.kilojoule/unit.mole, maxIterations=0):
         """Perform a local energy minimization on the system.
 
         Parameters:
-         - tolerance (energy=1*kilojoule/mole) The energy tolerance to which the system should be minimized
+         - tolerance (energy=10*kilojoules/mole) The energy tolerance to which the system should be minimized
          - maxIterations (int=0) The maximum number of iterations to perform.  If this is 0, minimization is continued
            until the results converge without regard to how many iterations it takes.
         """

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

@@ -26,7 +26,7 @@
   }
 
 
-  %pythoncode {
+  %pythoncode %{
     def getState(self,
                  getPositions=False,
                  getVelocities=False,
@@ -107,7 +107,7 @@
         if state._paramMap is not None:
              for param in state._paramMap:
                  self.setParameter(param, state._paramMap[param])
-  }
+  %}
   
   %feature("docstring") createCheckpoint "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.
@@ -175,7 +175,7 @@ Parameters:
   }
 
 
-  %pythoncode {
+  %pythoncode %{
     def getState(self,
                  copy,
                  getPositions=False,
@@ -235,11 +235,11 @@ Parameters:
                       periodicBoxVectorsList=periodicBoxVectorsList,
                       paramMap=paramMap)
         return state
-  }
+  %}
 }
 
 %extend OpenMM::NonbondedForce {
-  %pythoncode {
+  %pythoncode %{
     def addParticle_usingRVdw(self, charge, rVDW, epsilon):
         """Add particle using elemetrary charge.  Rvdw and epsilon,
            which is consistent with AMBER parameter file usage.
@@ -260,11 +260,11 @@ Parameters:
         """
         return self.addException(particle1, particle2,
                                  chargeProd, rMin/RMIN_PER_SIGMA, epsilon)
-  }
+  %}
 }
 
 %extend OpenMM::System {
-  %pythoncode {
+  %pythoncode %{
     def __getstate__(self):
         serializationString = XmlSerializer.serializeSystem(self)
         return serializationString
@@ -275,7 +275,7 @@ Parameters:
     def getForces(self):
         """Get the list of Forces in this System"""
         return [self.getForce(i) for i in range(self.getNumForces())]
-  }
+  %}
 }
 
 %extend OpenMM::XmlSerializer {
@@ -345,7 +345,7 @@ Parameters:
     return obj;
   }
 
-  %pythoncode {
+  %pythoncode %{
     @staticmethod
     def _serializeState(pythonState):
       positions = []
@@ -431,7 +431,7 @@ Parameters:
       if type == "State":
         return XmlSerializer._deserializeState(inputString)
       raise ValueError("Unsupported object type")
-  }
+  %}
 }
 
 %extend OpenMM::CustomIntegrator {
@@ -443,7 +443,7 @@ Parameters:
 }
 
 %extend OpenMM::Force {
-  %pythoncode {
+  %pythoncode %{
     def __copy__(self):
         copy = self.__class__.__new__(self.__class__)
         copy.__init__(self)
@@ -451,11 +451,11 @@ Parameters:
 
     def __deepcopy__(self, memo):
         return self.__copy__()
-  }
+  %}
 }
 
 %extend OpenMM::Integrator {
-  %pythoncode {
+  %pythoncode %{
     def __getstate__(self):
         serializationString = XmlSerializer.serialize(self)
         return serializationString
@@ -463,5 +463,5 @@ Parameters:
     def __setstate__(self, serializationString):
         system = XmlSerializer.deserialize(serializationString)
         self.this = system.this
-  }
+  %}
 }

wrappers/python/tests/TestAmberPrmtopFile.py

@@ -317,6 +317,16 @@ class TestAmberPrmtopFile(unittest.TestCase):
         simulation.step(5)
         os.remove(fname)
 
+    def testChamber(self):
+        """ Tests that Chamber prmtops fail with proper error message """
+        self.assertRaises(TypeError, lambda: AmberPrmtopFile('systems/ala3_solv.parm7'))
+        try:
+            parm = AmberPrmtopFile('systems/ala3_solv.parm7')
+            # Should not make it past here
+            self.assertTrue(False)
+        except TypeError as e:
+            # Make sure it says something about chamber
+            self.assertTrue('chamber' in str(e).lower())
 
 if __name__ == '__main__':
     unittest.main()

wrappers/python/tests/TestModeller.py

@@ -308,7 +308,7 @@ class TestModeller(unittest.TestCase):
                         self.assertTrue(len(matoms)==0 and len(m1atoms)==1 and len(m2atoms)==1)
     
     def test_addSolventPeriodicBox(self):
-        """ Test the addSolvent() method; test that the four ways of passing in the periodic box all work. """
+        """ Test the addSolvent() method; test that the five ways of passing in the periodic box all work. """
     
         # First way of passing in periodic box vectors:  set it in the original topology.
         topology_start = self.pdb.topology
@@ -359,6 +359,14 @@ class TestModeller(unittest.TestCase):
         self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0, 2.8802, 0))
         self.assertVecAlmostEqual(dim3[2]/nanometers, Vec3(0, 0, 2.8802))
         
+        # Fifth way: specify a number of molecules to add instead of a box size
+        topology_start = self.pdb.topology
+        modeller = Modeller(topology_start, self.positions)
+        modeller.deleteWater()
+        numInitial = len(list(modeller.topology.residues()))
+        modeller.addSolvent(self.forcefield, numAdded=1000)
+        self.assertEqual(numInitial+1000, len(list(modeller.topology.residues())))
+
     def test_addSolventNeutralSolvent(self):
         """ Test the addSolvent() method; test adding ions to neutral solvent. """