Merge branch 'master' of github.com:SimTk/openmm

examples/CMakeLists.txt

@@ -93,7 +93,7 @@ FOREACH(EX_ROOT ${F_EXAMPLES})
     INSTALL(FILES ${EX_ROOT}.f90 DESTINATION examples)
 ENDFOREACH(EX_ROOT ${F_EXAMPLES})
 
-INSTALL(FILES simulateAmber.py simulatePdb.py testInstallation.py argon-chemical-potential.py input.inpcrd input.prmtop input.pdb
+INSTALL(FILES simulateAmber.py simulatePdb.py simulateGromacs.py testInstallation.py argon-chemical-potential.py input.inpcrd input.prmtop input.pdb input.gro input.top
         DESTINATION examples)
 
 INSTALL(FILES VisualStudio/HelloArgon.vcproj 

examples/argon-chemical-potential.py

@@ -15,6 +15,7 @@
 # [1] Michael R. Shirts and John D. Chodera. Statistically optimal analysis of samples from multiple equilibrium states.
 # J. Chem. Phys. 129:124105 (2008)  http://dx.doi.org/10.1063/1.2978177
 
+from __future__ import print_function
 from simtk.openmm import *
 from simtk.unit import *
 import numpy
@@ -51,15 +52,15 @@ nlambda = len(lambda_values)
 volume = nparticles*(sigma**3)/reduced_density
 box_edge = volume**(1.0/3.0)
 cutoff = min(box_edge*0.49, 2.5*sigma) # Compute cutoff
-print "sigma = %s" % sigma
-print "box_edge = %s" % box_edge
-print "cutoff = %s" % cutoff
+print("sigma = %s" % sigma)
+print("box_edge = %s" % box_edge)
+print("cutoff = %s" % cutoff)
 
 # =============================================================================
 # Build systems at each alchemical lambda value.
 # =============================================================================
 
-print "Building alchemically-modified systems..."
+print("Building alchemically-modified systems...")
 alchemical_systems = list() # alchemical_systems[i] is the alchemically-modified System object corresponding to lambda_values[i]
 for lambda_value in lambda_values:
    # Create argon system where first particle is alchemically modified by lambda_value.
@@ -112,17 +113,17 @@ for (lambda_index, lambda_value) in enumerate(lambda_values):
    context.setPositions(positions)
 
    # Minimize energy from coordinates.
-   print "Lambda %d/%d : minimizing..." % (lambda_index+1, nlambda)
+   print("Lambda %d/%d : minimizing..." % (lambda_index+1, nlambda))
    LocalEnergyMinimizer.minimize(context)
 
    # Equilibrate.
-   print "Lambda %d/%d : equilibrating..." % (lambda_index+1, nlambda)
+   print("Lambda %d/%d : equilibrating..." % (lambda_index+1, nlambda))
    integrator.step(nequil_steps)
 
    # Sample.
    position_history = list() # position_history[i] is the set of positions after iteration i
    for iteration in range(nprod_iterations):
-      print "Lambda %d/%d : production iteration %d/%d" % (lambda_index+1, nlambda, iteration+1, nprod_iterations)
+      print("Lambda %d/%d : production iteration %d/%d" % (lambda_index+1, nlambda, iteration+1, nprod_iterations))
 
       # Run dynamics.
       integrator.step(nprod_steps)
@@ -138,7 +139,7 @@ for (lambda_index, lambda_value) in enumerate(lambda_values):
    del context, integrator
 
    # Compute reduced potentials of all snapshots at all alchemical states for MBAR.
-   print "Lambda %d/%d : computing energies at all states..." % (lambda_index+1, nlambda)
+   print("Lambda %d/%d : computing energies at all states..." % (lambda_index+1, nlambda))
    beta = 1.0 / (BOLTZMANN_CONSTANT_kB * AVOGADRO_CONSTANT_NA * temperature) # inverse temperature
    for l in range(nlambda):
       # Set up Context just to evaluate energies.
@@ -159,13 +160,13 @@ try:
    from timeseries import subsampleCorrelatedData
    from pymbar import MBAR
 except:
-   raise "pymbar [https://simtk.org/home/pymbar] must be installed to complete analysis of free energies."
+   raise ImportError("pymbar [https://simtk.org/home/pymbar] must be installed to complete analysis of free energies.")
    
 # =============================================================================
 # Subsample correlated samples to generate uncorrelated subsample.
 # =============================================================================
 
-print "Subsampling data to remove correlation..."
+print("Subsampling data to remove correlation...")
 K = nlambda # number of states
 N_k = nprod_iterations*numpy.ones([K], numpy.int32) # N_k[k] is the number of uncorrelated samples at state k
 u_kln_subsampled = numpy.zeros([K,K,nprod_iterations], numpy.float64) # subsampled data
@@ -176,24 +177,24 @@ for k in range(K):
    N_k[k] = len(indices)
    for l in range(K):
       u_kln_subsampled[k,l,0:len(indices)] = u_kln[k,l,indices]
-print "Number of uncorrelated samples per state:"
-print N_k
+print("Number of uncorrelated samples per state:")
+print(N_k)
 
 # =============================================================================
 # Analyze with MBAR to compute free energy differences and statistical errors.
 # =============================================================================
 
-print "Analyzing with MBAR..."
+print("Analyzing with MBAR...")
 mbar = MBAR(u_kln_subsampled, N_k)
-[Deltaf_ij, dDeltaf_ij] = mbar.getFreeEnergyDifferences()
-print "Free energy differences (in kT)"
-print Deltaf_ij
-print "Statistical errors (in kT)"
-print dDeltaf_ij
+Deltaf_ij, dDeltaf_ij = mbar.getFreeEnergyDifferences()
+print("Free energy differences (in kT)")
+print(Deltaf_ij)
+print("Statistical errors (in kT)")
+print(dDeltaf_ij)
 
 # =============================================================================
 # Report result.
 # =============================================================================
 
-print "Free energy of inserting argon particle: %.3f +- %.3f kT" % (Deltaf_ij[0,K-1], dDeltaf_ij[0,K-1])
+print("Free energy of inserting argon particle: %.3f +- %.3f kT" % (Deltaf_ij[0,K-1], dDeltaf_ij[0,K-1]))
 

examples/testInstallation.py

+from __future__ import print_function
 # First make sure OpenMM is installed.
 
+
 import sys
 try:
     from simtk.openmm.app import *
     from simtk.openmm import *
     from simtk.unit import *
 except ImportError as err:
-    print "Failed to import OpenMM packages:", err.message
-    print "Make sure OpenMM is installed and the library path is set correctly."
+    print("Failed to import OpenMM packages:", err.message)
+    print("Make sure OpenMM is installed and the library path is set correctly.")
     sys.exit()
 
 # Create a System for the tests.
@@ -19,28 +21,28 @@ system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, nonbondedCut
 # List all installed platforms and compute forces with each one.
 
 numPlatforms = Platform.getNumPlatforms()
-print "There are", numPlatforms, "Platforms available:"
-print
+print("There are", numPlatforms, "Platforms available:")
+print()
 forces = [None]*numPlatforms
 for i in range(numPlatforms):
     platform = Platform.getPlatform(i)
-    print i+1, platform.getName(),
+    print(i+1, platform.getName(), end=" ")
     integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
     try:
         simulation = Simulation(pdb.topology, system, integrator, platform)
         simulation.context.setPositions(pdb.positions)
         forces[i] = simulation.context.getState(getForces=True).getForces()
         del simulation
-        print "- Successfully computed forces"
+        print("- Successfully computed forces")
     except:
-        print "- Error computing forces with", platform.getName(), "platform"
+        print("- Error computing forces with", platform.getName(), "platform")
 
 # See how well the platforms agree.
 
 if numPlatforms > 1:
-    print
-    print "Median difference in forces between platforms:"
-    print
+    print()
+    print("Median difference in forces between platforms:")
+    print()
     for i in range(numPlatforms):
         for j in range(i):
             if forces[i] is not None and forces[j] is not None:
@@ -49,4 +51,6 @@ if numPlatforms > 1:
                     d = f1-f2
                     error = sqrt((d[0]*d[0]+d[1]*d[1]+d[2]*d[2])/(f1[0]*f1[0]+f1[1]*f1[1]+f1[2]*f1[2]))
                     errors.append(error)
-                print "%s vs. %s: %g" % (Platform.getPlatform(j).getName(), Platform.getPlatform(i).getName(), sorted(errors)[len(errors)/2])
+                print("{} vs. {}: {:g}".format(Platform.getPlatform(j).getName(),
+                                              Platform.getPlatform(i).getName(),
+                                              sorted(errors)[len(errors)//2]))

platforms/opencl/include/OpenCLContext.h

@@ -191,6 +191,12 @@ public:
     int getDeviceIndex() {
         return deviceIndex;
     }
+    /**
+     * Get the index of the cl::Platform associated with this object.
+     */
+    int getPlatformIndex() {
+        return platformIndex;
+    }
     /**
      * Get the PlatformData object this context is part of.
      */
@@ -604,6 +610,7 @@ private:
     double time;
     OpenCLPlatform::PlatformData& platformData;
     int deviceIndex;
+    int platformIndex;
     int contextIndex;
     int stepCount;
     int computeForceCount;

platforms/opencl/src/OpenCLContext.cpp

@@ -88,17 +88,25 @@ OpenCLContext::OpenCLContext(const System& system, int platformIndex, int device
         contextIndex = platformData.contexts.size();
         std::vector<cl::Platform> platforms;
         cl::Platform::get(&platforms);
-        if (platformIndex < 0 || platformIndex >= (int) platforms.size())
-            throw OpenMMException("Illegal value for OpenCL platform index");
-        string platformVendor = platforms[platformIndex].getInfo<CL_PLATFORM_VENDOR>();
-        vector<cl::Device> devices;
-        platforms[platformIndex].getDevices(CL_DEVICE_TYPE_ALL, &devices);
         const int minThreadBlockSize = 32;
-        if (deviceIndex < 0 || deviceIndex >= (int) devices.size()) {
-            // Try to figure out which device is the fastest.
 
-            int bestSpeed = -1;
+        int bestSpeed = -1;
+        int bestDevice = -1;
+        int bestPlatform = -1;
+        for (int j = 0; j < platforms.size(); j++) {
+            // if they supplied a valid platformIndex, we only look through that platform
+            if (j != platformIndex && platformIndex >= 0 && platformIndex < (int) platforms.size())
+                continue;
+
+            string platformVendor = platforms[j].getInfo<CL_PLATFORM_VENDOR>();
+            vector<cl::Device> devices;
+            platforms[j].getDevices(CL_DEVICE_TYPE_ALL, &devices);
+
             for (int i = 0; i < (int) devices.size(); i++) {
+                // if they supplied a valid deviceIndex, we only look through that one
+                if (i != deviceIndex && deviceIndex >= 0 && deviceIndex < (int) devices.size())
+                    continue;
+
                 if (platformVendor == "Apple" && devices[i].getInfo<CL_DEVICE_VENDOR>() == "AMD")
                     continue; // Don't use AMD GPUs on OS X due to serious bugs.
                 int maxSize = devices[i].getInfo<CL_DEVICE_MAX_WORK_ITEM_SIZES>()[0];
@@ -137,15 +145,26 @@ OpenCLContext::OpenCLContext(const System& system, int platformIndex, int device
                 }
                 int speed = devices[i].getInfo<CL_DEVICE_MAX_COMPUTE_UNITS>()*processingElementsPerComputeUnit*devices[i].getInfo<CL_DEVICE_MAX_CLOCK_FREQUENCY>();
                 if (maxSize >= minThreadBlockSize && speed > bestSpeed) {
-                    deviceIndex = i;
+                    bestDevice = i;
                     bestSpeed = speed;
+                    bestPlatform = j;
                 }
             }
         }
-        if (deviceIndex == -1)
+
+        if (bestPlatform == -1)
+            throw OpenMMException("No compatible OpenCL platform is available");
+
+        if (bestDevice == -1)
             throw OpenMMException("No compatible OpenCL device is available");
-        device = devices[deviceIndex];
-        this->deviceIndex = deviceIndex;
+
+        vector<cl::Device> devices;
+        platforms[bestPlatform].getDevices(CL_DEVICE_TYPE_ALL, &devices);
+        string platformVendor = platforms[bestPlatform].getInfo<CL_PLATFORM_VENDOR>();
+        device = devices[bestDevice];
+
+        this->deviceIndex = bestDevice;
+        this->platformIndex = bestPlatform;
         if (device.getInfo<CL_DEVICE_MAX_WORK_GROUP_SIZE>() < minThreadBlockSize)
             throw OpenMMException("The specified OpenCL device is not compatible with OpenMM");
         compilationDefines["WORK_GROUP_SIZE"] = intToString(ThreadBlockSize);
@@ -227,7 +246,7 @@ OpenCLContext::OpenCLContext(const System& system, int platformIndex, int device
             compilationDefines["SYNC_WARPS"] = "barrier(CLK_LOCAL_MEM_FENCE)";
         vector<cl::Device> contextDevices;
         contextDevices.push_back(device);
-        cl_context_properties cprops[] = {CL_CONTEXT_PLATFORM, (cl_context_properties) platforms[platformIndex](), 0};
+        cl_context_properties cprops[] = {CL_CONTEXT_PLATFORM, (cl_context_properties) platforms[bestPlatform](), 0};
         context = cl::Context(contextDevices, cprops, errorCallback);
         queue = cl::CommandQueue(context, device);
         numAtoms = system.getNumParticles();

platforms/opencl/src/OpenCLPlatform.cpp

@@ -133,7 +133,7 @@ void OpenCLPlatform::contextDestroyed(ContextImpl& context) const {
 
 OpenCLPlatform::PlatformData::PlatformData(const System& system, const string& platformPropValue, const string& deviceIndexProperty,
         const string& precisionProperty, const string& cpuPmeProperty) : removeCM(false), stepCount(0), computeForceCount(0), time(0.0)  {
-    int platformIndex = 0;
+    int platformIndex = -1;
     if (platformPropValue.length() > 0)
         stringstream(platformPropValue) >> platformIndex;
     vector<string> devices;
@@ -161,6 +161,8 @@ OpenCLPlatform::PlatformData::PlatformData(const System& system, const string& p
         deviceIndex << contexts[i]->getDeviceIndex();
         deviceName << contexts[i]->getDevice().getInfo<CL_DEVICE_NAME>();
     }
+    platformIndex = contexts[0]->getPlatformIndex();
+
     useCpuPme = (cpuPmeProperty == "true" && !contexts[0]->getUseDoublePrecision());
     propertyValues[OpenCLPlatform::OpenCLDeviceIndex()] = deviceIndex.str();
     propertyValues[OpenCLPlatform::OpenCLDeviceName()] = deviceName.str();

plugins/cpupme/src/CpuPmeKernels.cpp

@@ -95,6 +95,7 @@ static int getNumProcessors() {
 #define cpuid __cpuid
 #else
 static void cpuid(int cpuInfo[4], int infoType){
+#ifdef __LP64__
     __asm__ __volatile__ (
         "cpuid":
         "=a" (cpuInfo[0]),
@@ -103,6 +104,19 @@ static void cpuid(int cpuInfo[4], int infoType){
         "=d" (cpuInfo[3]) :
         "a" (infoType)
     );
+#else
+    __asm__ __volatile__ (
+        "pushl %%ebx\n"
+        "cpuid\n"
+        "movl %%ebx, %1\n"
+        "popl %%ebx\n" :
+        "=a" (cpuInfo[0]),
+        "=r" (cpuInfo[1]),
+        "=c" (cpuInfo[2]),
+        "=d" (cpuInfo[3]) :
+        "a" (infoType)
+    );
+#endif
 }
 #endif
 

plugins/rpmd/platforms/cuda/src/kernels/rpmdContraction.cu

@@ -114,8 +114,8 @@ extern "C" __global__ void contractForces(long long* force, long long* contracte
         
         // Store results.
         
-        force[forceIndex] = (long long) (FORCE_SCALE*freal[indexInBlock].x);
-        force[forceIndex+PADDED_NUM_ATOMS] = (long long) (FORCE_SCALE*freal[indexInBlock].y);
-        force[forceIndex+PADDED_NUM_ATOMS*2] = (long long) (FORCE_SCALE*freal[indexInBlock].z);
+        force[forceIndex] += (long long) (FORCE_SCALE*freal[indexInBlock].x);
+        force[forceIndex+PADDED_NUM_ATOMS] += (long long) (FORCE_SCALE*freal[indexInBlock].y);
+        force[forceIndex+PADDED_NUM_ATOMS*2] += (long long) (FORCE_SCALE*freal[indexInBlock].z);
     }
 }

plugins/rpmd/platforms/opencl/src/kernels/rpmdContraction.cl

@@ -113,6 +113,6 @@ __kernel void contractForces(__global real4* force, __global real4* contracted)
         
         // Store results.
 
-        force[index] = convert_real4(FORCE_SCALE*freal[indexInBlock]);
+        force[index] += convert_real4(FORCE_SCALE*freal[indexInBlock]);
     }
 }

plugins/rpmd/platforms/reference/src/ReferenceRpmdKernels.cpp

@@ -336,7 +336,7 @@ void ReferenceIntegrateRPMDStepKernel::computeForces(ContextImpl& context, const
                     q[k] = t_complex(0, 0);
                 fftpack_exec_1d(fft, FFTPACK_BACKWARD, &q[0], &q[0]);
                 for (int k = 0; k < totalCopies; k++)
-                    forces[k][particle][component] = scale2*q[k].re;
+                    forces[k][particle][component] += scale2*q[k].re;
             }
         }
     }

wrappers/python/CMakeLists.txt

@@ -13,6 +13,7 @@ mark_as_advanced(OPENMM_PYTHON_STAGING_DIR)
 file(MAKE_DIRECTORY ${OPENMM_PYTHON_STAGING_DIR}/simtk/openmm)
 file(MAKE_DIRECTORY ${OPENMM_PYTHON_STAGING_DIR}/simtk/unit)
 file(MAKE_DIRECTORY ${OPENMM_PYTHON_STAGING_DIR}/src/swig_doxygen/swig_lib/python)
+file(MAKE_DIRECTORY ${OPENMM_PYTHON_STAGING_DIR}/tests)
 
 ##############################################################################
 ### Identify files that need to be copied from source area to staging area ###
@@ -46,7 +47,7 @@ file(WRITE "${OPENMM_PYTHON_STAGING_DIR}/simtk/openmm/version.py" "git_revision
 #     set(temp2 "${temp2}\n${f}")
 # endforeach()
 
-set(SUBDIRS src simtk)
+set(SUBDIRS src simtk tests)
 foreach(SUBDIR ${SUBDIRS})
     file(GLOB_RECURSE STAGING_INPUT_FILES1 RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*README.txt"
@@ -55,6 +56,8 @@ foreach(SUBDIR ${SUBDIRS})
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.sh"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.xml"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.pdb"
+        "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.prmtop"
+        "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.top"
     )
     foreach(file ${STAGING_INPUT_FILES1})
         set(STAGING_INPUT_FILES ${STAGING_INPUT_FILES} "${file}")

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

@@ -335,8 +335,8 @@ class PrmtopLoader(object):
                                  % ((bondPointers[ii],
                                      bondPointers[ii+1]),))
              iType=int(bondPointers[ii+2])-1
-             returnList.append((int(bondPointers[ii])/3,
-                                int(bondPointers[ii+1])/3,
+             returnList.append((int(bondPointers[ii])//3,
+                                int(bondPointers[ii+1])//3,
                                 float(forceConstant[iType])*forceConstConversionFactor,
                                 float(bondEquil[iType])*lengthConversionFactor))
         return returnList
@@ -383,9 +383,9 @@ class PrmtopLoader(object):
                                      anglePointers[ii+1],
                                      anglePointers[ii+2]),))
              iType=int(anglePointers[ii+3])-1
-             self._angleList.append((int(anglePointers[ii])/3,
-                                int(anglePointers[ii+1])/3,
-                                int(anglePointers[ii+2])/3,
+             self._angleList.append((int(anglePointers[ii])//3,
+                                int(anglePointers[ii+1])//3,
+                                int(anglePointers[ii+2])//3,
                                 float(forceConstant[iType])*forceConstConversionFactor,
                                 float(angleEquil[iType])))
         return self._angleList
@@ -411,10 +411,10 @@ class PrmtopLoader(object):
                                     dihedralPointers[ii+2],
                                     dihedralPointers[ii+3]),))
              iType=int(dihedralPointers[ii+4])-1
-             self._dihedralList.append((int(dihedralPointers[ii])/3,
-                                int(dihedralPointers[ii+1])/3,
-                                abs(int(dihedralPointers[ii+2]))/3,
-                                abs(int(dihedralPointers[ii+3]))/3,
+             self._dihedralList.append((int(dihedralPointers[ii])//3,
+                                int(dihedralPointers[ii+1])//3,
+                                abs(int(dihedralPointers[ii+2]))//3,
+                                abs(int(dihedralPointers[ii+3]))//3,
                                 float(forceConstant[iType])*forceConstConversionFactor,
                                 float(phase[iType]),
                                 int(0.5+float(periodicity[iType]))))
@@ -429,8 +429,8 @@ class PrmtopLoader(object):
         nonbondTerms = self.getNonbondTerms()
         for ii in range(0,len(dihedralPointers),5):
              if int(dihedralPointers[ii+2])>0 and int(dihedralPointers[ii+3])>0:
-                 iAtom = int(dihedralPointers[ii])/3
-                 lAtom = int(dihedralPointers[ii+3])/3
+                 iAtom = int(dihedralPointers[ii])//3
+                 lAtom = int(dihedralPointers[ii+3])//3
                  chargeProd = charges[iAtom]*charges[lAtom]
                  (rVdwI, epsilonI) = nonbondTerms[iAtom]
                  (rVdwL, epsilonL) = nonbondTerms[lAtom]

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

@@ -769,6 +769,7 @@ class Modeller(object):
         LocalEnergyMinimizer.minimize(context)
         self.topology = newTopology
         self.positions = context.getState(getPositions=True).getPositions()
+        del context
         return actualVariants
 
     def addExtraParticles(self, forcefield):

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

@@ -275,6 +275,7 @@ Parameters:
   }
 
   %feature(docstring, "This method exists only for backward compatibility. @deprecated Use deserialize() instead.") deserializeSystem;
+  %newobject deserializeSystem;
   static OpenMM::System* deserializeSystem(const char* inputString) {
       std::stringstream ss;
       ss << inputString;
@@ -287,6 +288,7 @@ Parameters:
       return ss.str();
   }
 
+  %newobject _deserializeForce;
   static OpenMM::Force* _deserializeForce(const char* inputString) {
       std::stringstream ss;
       ss << inputString;
@@ -299,6 +301,7 @@ Parameters:
       return ss.str();
   }
 
+  %newobject _deserializeIntegrator;
   static OpenMM::Integrator* _deserializeIntegrator(const char* inputString) {
       std::stringstream ss;
       ss << inputString;