Reimplement constant acceleration groups
[alexxy/gromacs.git] / src / gromacs / modularsimulator / modularsimulator.cpp
index 4e91aa95d633aaf613b979739f3afc8df147f1ec..03e64cc6e36c6fdb4c95fd04625eab68fed98a4a 100644 (file)
 #include "andersentemperaturecoupling.h"
 #include "computeglobalselement.h"
 #include "constraintelement.h"
+#include "expandedensembleelement.h"
 #include "firstorderpressurecoupling.h"
 #include "forceelement.h"
 #include "mttk.h"
 #include "nosehooverchains.h"
 #include "parrinellorahmanbarostat.h"
+#include "pullelement.h"
 #include "simulatoralgorithm.h"
 #include "statepropagatordata.h"
 #include "velocityscalingtemperaturecoupling.h"
@@ -132,6 +134,12 @@ void ModularSimulator::addIntegrationElements(ModularSimulatorAlgorithmBuilder*
         {
             builder->add<ConstraintsElement<ConstraintVariable::Positions>>();
         }
+
+        if (legacySimulatorData_->inputrec->bPull)
+        {
+            builder->add<PullElement>();
+        }
+
         builder->add<ComputeGlobalsElement<ComputeGlobalsAlgorithm::LeapFrog>>();
         if (legacySimulatorData_->inputrec->epc == PressureCoupling::ParrinelloRahman)
         {
@@ -154,7 +162,12 @@ void ModularSimulator::addIntegrationElements(ModularSimulatorAlgorithmBuilder*
             builder->add<ConstraintsElement<ConstraintVariable::Velocities>>();
         }
         builder->add<ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>>();
+        // Here, we have x / v / f at the full time step
         builder->add<StatePropagatorData::Element>();
+        if (legacySimulatorData_->inputrec->bExpanded)
+        {
+            builder->add<ExpandedEnsembleElement>();
+        }
         if (legacySimulatorData_->inputrec->etc == TemperatureCoupling::VRescale
             || legacySimulatorData_->inputrec->etc == TemperatureCoupling::Berendsen)
         {
@@ -175,6 +188,12 @@ void ModularSimulator::addIntegrationElements(ModularSimulatorAlgorithmBuilder*
         {
             builder->add<ConstraintsElement<ConstraintVariable::Positions>>();
         }
+
+        if (legacySimulatorData_->inputrec->bPull)
+        {
+            builder->add<PullElement>();
+        }
+
         builder->add<ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>>();
         if (legacySimulatorData_->inputrec->epc == PressureCoupling::ParrinelloRahman)
         {
@@ -246,6 +265,10 @@ void ModularSimulator::addIntegrationElements(ModularSimulatorAlgorithmBuilder*
         }
         // We have a full state at time t here
         builder->add<StatePropagatorData::Element>();
+        if (legacySimulatorData_->inputrec->bExpanded)
+        {
+            builder->add<ExpandedEnsembleElement>();
+        }
 
         // Propagate extended system variables from t to t+dt/2
         if (legacySimulatorData_->inputrec->epc == PressureCoupling::Mttk)
@@ -294,6 +317,12 @@ void ModularSimulator::addIntegrationElements(ModularSimulatorAlgorithmBuilder*
         {
             builder->add<ConstraintsElement<ConstraintVariable::Positions>>();
         }
+
+        if (legacySimulatorData_->inputrec->bPull)
+        {
+            builder->add<PullElement>();
+        }
+
         builder->add<ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>>();
 
         // Propagate box from t to t+dt
@@ -372,17 +401,7 @@ bool ModularSimulator::isInputCompatible(bool                             exitOn
             && conditionalAssert(!doRerun, "Rerun is not supported by the modular simulator.");
     isInputCompatible =
             isInputCompatible
-            && conditionalAssert(inputrec->efep == FreeEnergyPerturbationType::No
-                                         || inputrec->efep == FreeEnergyPerturbationType::Yes
-                                         || inputrec->efep == FreeEnergyPerturbationType::SlowGrowth,
-                                 "Expanded ensemble free energy calculation is not "
-                                 "supported by the modular simulator.");
-    isInputCompatible = isInputCompatible
-                        && conditionalAssert(!inputrec->bPull,
-                                             "Pulling is not supported by the modular simulator.");
-    isInputCompatible =
-            isInputCompatible
-            && conditionalAssert(inputrec->cos_accel == 0.0,
+            && conditionalAssert(!inputrec->useConstantAcceleration && inputrec->cos_accel == 0.0,
                                  "Acceleration is not supported by the modular simulator.");
     isInputCompatible =
             isInputCompatible
@@ -447,10 +466,6 @@ bool ModularSimulator::isInputCompatible(bool                             exitOn
             isInputCompatible
             && conditionalAssert(!inputrec->bSimTemp,
                                  "Simulated tempering is not supported by the modular simulator.");
-    isInputCompatible = isInputCompatible
-                        && conditionalAssert(!inputrec->bExpanded,
-                                             "Expanded ensemble simulations are not supported by "
-                                             "the modular simulator.");
     isInputCompatible =
             isInputCompatible
             && conditionalAssert(!doEssentialDynamics,