Fixed a bug in identifying molecules

platforms/cuda/src/kernels/gpu.cpp

@@ -103,6 +103,7 @@ struct Molecule {
     vector<int> periodicTorsions;
     vector<int> rbTorsions;
     vector<int> constraints;
+    vector<int> lj14s;
 };
 
 static const float dielectricOffset         =    0.009f;
@@ -2072,7 +2073,7 @@ static void findMoleculeGroups(gpuContext gpu)
         constraints.push_back(Constraint(atom1, atom2, distance2));
     }
 
-    // First make a list of every other atom to which each atom is connect by a bond or constraint.
+    // First make a list of every other atom to which each atom is connect by a bond, constraint, or exclusion.
 
     int numAtoms = gpu->natoms;
     vector<vector<int> > atomBonds(numAtoms);
@@ -2090,6 +2091,9 @@ static void findMoleculeGroups(gpuContext gpu)
         atomBonds[atom1].push_back(atom2);
         atomBonds[atom2].push_back(atom1);
     }
+    for (int i = 0; i < (int)gpu->exclusions.size(); i++)
+        for (int j = 0; j < (int)gpu->exclusions[i].size(); j++)
+            atomBonds[i].push_back(gpu->exclusions[i][j]);
 
     // Now tag atoms by which molecule they belong to.
 
@@ -2131,6 +2135,11 @@ static void findMoleculeGroups(gpuContext gpu)
     {
         molecules[atomMolecule[constraints[i].atom1]].constraints.push_back(i);
     }
+    for (int i = 0; i < (int)gpu->sim.LJ14s; i++)
+    {
+        int atom1 = (*gpu->psLJ14ID)[i].x;
+        molecules[atomMolecule[atom1]].lj14s.push_back(i);
+    }
 
     // Sort them into groups of identical molecules.
 
@@ -2149,7 +2158,8 @@ static void findMoleculeGroups(gpuContext gpu)
             bool identical = true;
             if (mol.atoms.size() != mol2.atoms.size() || mol.bonds.size() != mol2.bonds.size()
                     || mol.angles.size() != mol2.angles.size() || mol.periodicTorsions.size() != mol2.periodicTorsions.size()
-                    || mol.rbTorsions.size() != mol2.rbTorsions.size() || mol.constraints.size() != mol2.constraints.size())
+                    || mol.rbTorsions.size() != mol2.rbTorsions.size() || mol.constraints.size() != mol2.constraints.size()
+                    || mol.lj14s.size() != mol2.lj14s.size())
                 identical = false;
             int atomOffset = mol2.atoms[0]-mol.atoms[0];
             float4* posq = gpu->psPosq4->_pSysData;
@@ -2206,6 +2216,13 @@ static void findMoleculeGroups(gpuContext gpu)
                         constraints[mol.constraints[i]].atom2 != constraints[mol2.constraints[i]].atom2-atomOffset ||
                         constraints[mol.constraints[i]].distance2 != constraints[mol2.constraints[i]].distance2)
                     identical = false;
+            int4* lj14ID = gpu->psLJ14ID->_pSysData;
+            float4* lj14Param = gpu->psLJ14Parameter->_pSysData;
+            for (int i = 0; i < (int)mol.lj14s.size() && identical; i++)
+                if (lj14ID[mol.lj14s[i]].x != lj14ID[mol2.lj14s[i]].x-atomOffset || lj14ID[mol.lj14s[i]].y != lj14ID[mol2.lj14s[i]].y-atomOffset ||
+                        lj14Param[mol.lj14s[i]].x != lj14Param[mol2.lj14s[i]].x || lj14Param[mol.lj14s[i]].y != lj14Param[mol2.lj14s[i]].y ||
+                        lj14Param[mol.lj14s[i]].z != lj14Param[mol2.lj14s[i]].z)
+                    identical = false;
             if (identical)
             {
                 moleculeInstances[j].push_back(mol.atoms[0]);

platforms/cuda/tests/TestCudaNonbondedForce.cpp

@@ -392,7 +392,7 @@ void testLargeSystem() {
     const int numMolecules = 600;
     const int numParticles = numMolecules*2;
     const double cutoff = 2.0;
-    const double boxSize = 10.0;
+    const double boxSize = 20.0;
     const double tol = 1e-3;
     CudaPlatform cuda;
     ReferencePlatform reference;
@@ -407,11 +407,11 @@ void testLargeSystem() {
     init_gen_rand(0);
     for (int i = 0; i < numMolecules; i++) {
         if (i < numMolecules/2) {
-            nonbonded->addParticle(1.0, 0.2, 0.1);
+            nonbonded->addParticle(-1.0, 0.2, 0.1);
             nonbonded->addParticle(1.0, 0.1, 0.1);
         }
         else {
-            nonbonded->addParticle(1.0, 0.2, 0.2);
+            nonbonded->addParticle(-1.0, 0.2, 0.2);
             nonbonded->addParticle(1.0, 0.1, 0.2);
         }
         positions[2*i] = Vec3(boxSize*genrand_real2(), boxSize*genrand_real2(), boxSize*genrand_real2());
@@ -419,6 +419,7 @@ void testLargeSystem() {
         velocities[2*i] = Vec3(genrand_real2(), genrand_real2(), genrand_real2());
         velocities[2*i+1] = Vec3(genrand_real2(), genrand_real2(), genrand_real2());
         bonds->addBond(2*i, 2*i+1, 1.0, 0.1);
+        nonbonded->addException(2*i, 2*i+1, 0.0, 0.15, 0.0);
     }
 
     // Try with cutoffs but not periodic boundary conditions, and make sure the Cuda and Reference
@@ -434,13 +435,14 @@ void testLargeSystem() {
     cudaContext.setVelocities(velocities);
     referenceContext.setPositions(positions);
     referenceContext.setVelocities(velocities);
-    State cudaState = cudaContext.getState(State::Positions | State::Velocities | State::Forces);
-    State referenceState = referenceContext.getState(State::Positions | State::Velocities | State::Forces);
+    State cudaState = cudaContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
+    State referenceState = referenceContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
     for (int i = 0; i < numParticles; i++) {
         ASSERT_EQUAL_VEC(cudaState.getPositions()[i], referenceState.getPositions()[i], tol);
         ASSERT_EQUAL_VEC(cudaState.getVelocities()[i], referenceState.getVelocities()[i], tol);
         ASSERT_EQUAL_VEC(cudaState.getForces()[i], referenceState.getForces()[i], tol);
     }
+    ASSERT_EQUAL_TOL(cudaState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
 
     // Now do the same thing with periodic boundary conditions.
 
@@ -452,8 +454,8 @@ void testLargeSystem() {
     cudaContext.setVelocities(velocities);
     referenceContext.setPositions(positions);
     referenceContext.setVelocities(velocities);
-    cudaState = cudaContext.getState(State::Positions | State::Velocities | State::Forces);
-    referenceState = referenceContext.getState(State::Positions | State::Velocities | State::Forces);
+    cudaState = cudaContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
+    referenceState = referenceContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
     for (int i = 0; i < numParticles; i++) {
         ASSERT_EQUAL_TOL(fmod(cudaState.getPositions()[i][0]-referenceState.getPositions()[i][0], boxSize), 0, tol);
         ASSERT_EQUAL_TOL(fmod(cudaState.getPositions()[i][1]-referenceState.getPositions()[i][1], boxSize), 0, tol);
@@ -461,6 +463,7 @@ void testLargeSystem() {
         ASSERT_EQUAL_VEC(cudaState.getVelocities()[i], referenceState.getVelocities()[i], tol);
         ASSERT_EQUAL_VEC(cudaState.getForces()[i], referenceState.getForces()[i], tol);
     }
+    ASSERT_EQUAL_TOL(cudaState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
 }
 
 void testBlockInteractions(bool periodic) {