Fixed several bugs in OpenMM

openmmapi/src/NonbondedForceImpl.cpp

@@ -51,21 +51,35 @@ void NonbondedForceImpl::initialize(OpenMMContextImpl& context) {
     
     System& system = context.getSystem();
     vector<set<int> > exclusions(owner.getNumParticles());
+    vector<vector<int> > bondIndices;
+    set<pair<int, int> > bonded14set;
     for (int i = 0; i < system.getNumForces(); i++) {
         if (dynamic_cast<HarmonicBondForce*>(&system.getForce(i)) != NULL) {
             const HarmonicBondForce& force = dynamic_cast<const HarmonicBondForce&>(system.getForce(i));
-            vector<vector<int> > bondIndices(force.getNumBonds());
-            set<pair<int, int> > bonded14set;
-            for (int i = 0; i < force.getNumBonds(); ++i) {
+            bondIndices.resize(force.getNumBonds());
+            for (int j = 0; j < force.getNumBonds(); ++j) {
                 int particle1, particle2;
                 double length, k;
-                force.getBondParameters(i, particle1, particle2, length, k);
-                bondIndices[i].push_back(particle1);
-                bondIndices[i].push_back(particle2);
+                force.getBondParameters(j, particle1, particle2, length, k);
+                bondIndices[j].push_back(particle1);
+                bondIndices[j].push_back(particle2);
             }
-            findExclusions(bondIndices, exclusions, bonded14set);
+            break;
         }
     }
+
+    // Also treat constrained distances as bonds.
+
+    int numBonds = bondIndices.size();
+    bondIndices.resize(numBonds+system.getNumConstraints());
+    for (int j = 0; j < system.getNumConstraints(); j++) {
+        int particle1, particle2;
+        double distance;
+        system.getConstraintParameters(j, particle1, particle2, distance);
+        bondIndices[j+numBonds].push_back(particle1);
+        bondIndices[j+numBonds].push_back(particle2);
+    }
+    findExclusions(bondIndices, exclusions, bonded14set);
     dynamic_cast<CalcNonbondedForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner, exclusions);
 }
 

platforms/cuda/src/CudaKernels.cpp

@@ -269,9 +269,8 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
             force.getNonbonded14Parameters(i, particle1[i], particle2[i], charge, sig, eps);
             c6[i] = (float) (4*eps*pow(sig, 6.0));
             c12[i] = (float) (4*eps*pow(sig, 12.0));
-            float q = (float) std::sqrt(charge);
-            q1[i] = q;
-            q2[i] = q;
+            q1[i] = (float) charge;
+            q2[i] = 1.0f;
         }
         gpuSetLJ14Parameters(gpu, 138.935485f, 1.0f, particle1, particle2, c6, c12, q1, q2);
     }

platforms/cuda/src/CudaPlatform.cpp

@@ -40,7 +40,8 @@
 using namespace OpenMM;
 
 extern "C" void initOpenMMPlugin() {
-    Platform::registerPlatform(new CudaPlatform());
+    if (gpuIsAvailable())
+        Platform::registerPlatform(new CudaPlatform());
 }
 
 CudaPlatform::CudaPlatform() {

platforms/reference/src/ReferenceKernels.cpp

@@ -378,8 +378,8 @@ double ReferenceCalcNonbondedForceKernel::executeEnergy(OpenMMContextImpl& conte
     ReferenceLJCoulomb14 nonbonded14;
     if (nonbondedMethod != NoCutoff)
         nonbonded14.setUseCutoff(nonbondedCutoff, 78.3f);
-    RealOpenMM* energyArray = new RealOpenMM[numParticles];
-    for (int i = 0; i < numParticles; ++i)
+    RealOpenMM* energyArray = new RealOpenMM[num14];
+    for (int i = 0; i < num14; ++i)
         energyArray[i] = 0;
     refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, energyArray, 0, &energy, nonbonded14);
     disposeRealArray(forceData, numParticles);

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp

@@ -243,7 +243,7 @@ int ReferenceLJCoulombIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM** ato
 
      --------------------------------------------------------------------------------------- */
 
-  int ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, RealOpenMM** atomCoordinates,
+int ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, RealOpenMM** atomCoordinates,
                         RealOpenMM** atomParameters, RealOpenMM** forces,
                         RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
 
@@ -294,7 +294,7 @@ int ReferenceLJCoulombIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM** ato
                else
                    dEdR += atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
                dEdR     *= inverseR*inverseR;
-         
+               
     // accumulate forces
 
     for( int kk = 0; kk < 3; kk++ ){