Additional test case (#4408)

* Improve numerical stability for edge cases

* Undo checking of derivative factors

* Undo checking for NaNs

tests/TestATMForce.h

@@ -47,6 +47,7 @@
 #include "openmm/serialization/XmlSerializer.h"
 #include "sfmt/SFMT.h"
 #include <algorithm>
+#include <random>
 #include <iostream>
 #include <vector>
 
@@ -240,47 +241,64 @@ void test2ParticlesSoftCore() {
     ASSERT_EQUAL_VEC(Vec3(-lmbd*df*displ[0], 0.0, 0.0), state.getForces()[1], 1e-6);
 }
 
-void test2ParticlesNonbonded() {
-
+void testNonbonded() {
     System system;
-    system.addParticle(1.0);
-    system.addParticle(1.0);
+    double u0, u1, energy;
+    double lambda = 0.5;
+    int numParticles = 216;
+    double width = 4.0;
 
+    system.setDefaultPeriodicBoxVectors(Vec3(width, 0, 0), Vec3(0, width, 0), Vec3(0, 0, width));
     NonbondedForce* nbforce = new NonbondedForce();
-    nbforce->addParticle( 1.0, 1.0, 1.0);
-    nbforce->addParticle(-1.0, 1.0, 1.0);
+    nbforce->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    nbforce->setCutoffDistance(0.7);
+    ATMForce* atm = new ATMForce(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0);
+    double spacing = width/6.0;
+    double offset = spacing/5.0;
+    vector<Vec3> positions;
+    for (int i = 0; i < 6; i++)
+        for (int j = 0; j < 6; j++)
+            for (int k = 0; k < 6; k++) {
+	        positions.push_back(Vec3(spacing*i+offset, spacing*j+offset, spacing*k+offset));
+		system.addParticle(10.0);
+		nbforce->addParticle(0, 0.3, 1.0);
+		atm->addParticle(Vec3());
+            }
+    auto rng = std::default_random_engine {};
+    std::shuffle(std::begin(positions), std::end(positions), rng);
+    atm->setParticleParameters(0, Vec3(0.5, 0, 0), Vec3(0.0, 0, 0));
 
+    //in this scenario the non-bonded force is added to the System, a copy is added to ATMForce and
+    //the System's copy is disabled by giving it a force group that is not evaluated.
+    //This should cause atom reordering in the main context
     system.addForce(nbforce);
-
-    ATMForce* atm = new ATMForce(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0);
-    Vec3 nodispl = Vec3(0., 0., 0.);
-    Vec3   displ = Vec3(1., 0., 0.);
-    atm->addParticle( nodispl );
-    atm->addParticle(   displ );
     atm->addForce(XmlSerializer::clone<Force>(*nbforce));
-    system.removeForce(0);
+    nbforce->setForceGroup(1);
     system.addForce(atm);
-
-    vector<Vec3> positions;
-    positions.push_back(Vec3(0., 0., 0.));
-    positions.push_back(Vec3(1., 0., 0.));
-
-    VerletIntegrator integrator(1.0);
-    Context context(system, integrator, platform);
-    context.setPositions(positions);
-
-    double lambda = 0.5;
-    context.setParameter(ATMForce::Lambda1(), lambda);
-    context.setParameter(ATMForce::Lambda2(), lambda);
-    State state = context.getState( State::Energy );
-    double epot = state.getPotentialEnergy();
-    double u0, u1, energy;
-    atm->getPerturbationEnergy(context, u1, u0, energy);
-    double epert = u1 - u0;
-    ASSERT_EQUAL_TOL(-104.2320, epot,  1e-3);
-    ASSERT_EQUAL_TOL(  69.4062, epert, 1e-3);
-    // std::cout << "Nonbonded: epot = " << epot << std::endl;
-    // std::cout << "Nonbonded: epert = " << epert << std::endl;
+    LangevinMiddleIntegrator integrator1(300, 1.0, 0.004);
+    Context context1(system, integrator1, platform);
+    context1.setPositions(positions);
+    context1.setParameter(ATMForce::Lambda1(), lambda);
+    context1.setParameter(ATMForce::Lambda2(), lambda);
+    State state1 = context1.getState( State::Energy, false, 0 );
+    double epot1 = state1.getPotentialEnergy();
+    atm->getPerturbationEnergy(context1, u1, u0, energy);
+    double epert1 = u1 - u0;
+
+    //in this second scenario the non-bonded force is remove from the System
+    //and atom reordering is disabled.
+    system.removeForce(0);
+    LangevinMiddleIntegrator integrator2(300, 1.0, 0.004);
+    Context context2(system, integrator2, platform);
+    context2.setPositions(positions);
+    context2.setParameter(ATMForce::Lambda1(), lambda);
+    context2.setParameter(ATMForce::Lambda2(), lambda);
+    State state2 = context2.getState( State::Energy );
+    double epot2 = state2.getPotentialEnergy();
+    atm->getPerturbationEnergy(context2, u1, u0, energy);
+    double epert2 = u1 - u0;
+
+    ASSERT_EQUAL_TOL(epert1, epert2,  1e-3);
 }
 
 
@@ -520,7 +538,7 @@ int main(int argc, char* argv[]) {
         test2Particles();
         test2Particles2Displacement0();
         test2ParticlesSoftCore();
-        test2ParticlesNonbonded();
+        testNonbonded();
         testParticlesCustomExpressionLinear();
         testParticlesCustomExpressionSoftplus();
         testLargeSystem();