Changeset 4882d5 for src

Timestamp:

Dec 14, 2012, 5:39:30 PM (12 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

5251af

Parents:

8b886f

git-author:

Frederik Heber <heber@…> (10/05/12 07:59:22)

git-committer:

Frederik Heber <heber@…> (12/14/12 17:39:30)

Message:

Rewrote VerletIntegrationAction to diminish dependence on ForceMatrix.

• atom::VelocityVerletUpdate() now works on the forces stored in AtomicForces.
• ThermoStatContainer::getActive() getter for the currently active thermostat.
• VerletForceIntegration::operator() split up into several smaller functions.

Location:

src

Files:

• Actions/MoleculeAction/VerletIntegrationAction.cpp (modified) (1 diff)
• Actions/MoleculeAction/VerletIntegrationAction.def (modified) (2 diffs)
• Atom/atom_atominfo.cpp (modified) (2 diffs)
• Atom/atom_atominfo.hpp (modified) (1 diff)
• Dynamics/VerletForceIntegration.hpp (modified) (4 diffs)
• Thermostats/ThermoStatContainer.hpp (modified) (1 diff)

src/Actions/MoleculeAction/VerletIntegrationAction.cpp

@@ -57 +57 @@
 /** =========== define the function ====================== */
 Action::state_ptr MoleculeVerletIntegrationAction::performCall() {
-  LOG(1, "Parsing forces file and Verlet integrating.");
   // TODO: sollte besser stream nutzen, nicht filename direkt (es sei denn, ist prefix), besser fuer unit test
-  char outputname[MAXSTRINGSIZE];
-  strcpy(outputname, params.forcesfile.get().string().c_str());
   AtomSetMixin<std::vector<atom *> > set(World::getInstance().getSelectedAtoms());
   for (unsigned int step = 0; step < params.MDSteps.get(); ++step) {
-    VerletForceIntegration<std::vector<atom *> > Verlet(set, params.Deltat.get(), step, false);
-    if (!Verlet(outputname, 1, 0, params.FixedCenterOfMass.get()))
+    VerletForceIntegration<std::vector<atom *> > Verlet(set, params.Deltat.get(), false);
+    // parse forces into next step
+    if (!params.forcesfile.get().string().empty()) {
+      LOG(1, "Parsing forces file.");
+      if (!Verlet.parseForcesFile(params.forcesfile.get().string().c_str(), step))
+        LOG(2, "File " << params.forcesfile.get() << " not found.");
+      else
+        LOG(2, "File " << params.forcesfile.get() << " found and parsed.");
+    }
+    // perform velocity verlet update
+    LOG(1, "Verlet integrating.");
+    if (!Verlet(step+1, 1, 0, params.FixedCenterOfMass.get()))
       LOG(2, "File " << params.forcesfile.get() << " not found.");
     else

src/Actions/MoleculeAction/VerletIntegrationAction.def

@@ -11 +11 @@
 #include "Parameters/Validators/DummyValidator.hpp"
 #include "Parameters/Validators/GenericValidators.hpp"
+#include "Parameters/Validators/Ops_Validator.hpp"
 #include "Parameters/Validators/Specific/FilePresentValidator.hpp"
 
@@ -17 +18 @@
 // "undefine" if no parameters are required, use (NOPARAM_DEFAULT) for each (undefined) default value
 #define paramtypes (boost::filesystem::path)(double)(unsigned int)(bool)
-#define paramtokens ("verlet-integration")("deltat")("MDSteps")("keep-fixed-CenterOfMass")
-#define paramdescriptions ("perform verlet integration of a given force file")("time step width")("number of MDSteps to integrate")("whether forces and velocities shall be corrected such that center of mass remains at rest")
-#undef paramdefaults
+#define paramtokens ("forces-file")("deltat")("MDSteps")("keep-fixed-CenterOfMass")
+#define paramdescriptions ("file containing")("time step width")("number of MDSteps to integrate")("whether forces and velocities shall be corrected such that center of mass remains at rest")
+#define paramdefaults (PARAM_DEFAULT(""))(NOPARAM_DEFAULT)(NOPARAM_DEFAULT)(NOPARAM_DEFAULT)
 #define paramreferences (forcesfile)(Deltat)(MDSteps)(FixedCenterOfMass)
 #define paramvalids \
-(FilePresentValidator()) \
+(DummyValidator< boost::filesystem::path >()) \
 (PositiveValidator< double >()) \
 (NotZeroValidator< unsigned int >()) \

src/Atom/atom_atominfo.cpp

@@ -584 +584 @@
  * \param Deltat time step width
  * \param IsAngstroem whether the force's underlying unit of length is angstroem or bohr radii
- * \param *Force matrix with forces
- */
-void AtomInfo::VelocityVerletUpdate(int nr, const unsigned int NextStep, double Deltat, bool IsAngstroem, ForceMatrix *Force, const size_t offset)
-{
+ */
+void AtomInfo::VelocityVerletUpdate(int nr, const unsigned int NextStep, double Deltat, bool IsAngstroem)
+{
+  const unsigned int LastStep = NextStep == 0 ? 0 : NextStep-1;
+
   LOG(2, "INFO: Particle that currently " << *this);
   LOG(2, "INFO: Integrating with mass=" << getMass() << " and Deltat="
@@ -593 +594 @@
   // update force
   // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
-  Vector tempVector;
-  for (int d=0; d<NDIM; d++) {
-    ASSERT(Force->RowCounter[0] > nr,
-        "AtomInfo::VelocityVerletUpdate() - "+toString(nr)
-        +" out of bounds [0,"+toString(Force->RowCounter[0])
-        +"] access on force matrix.");
-    ASSERT(Force->ColumnCounter[0] > d+(int)offset,
-        "AtomInfo::VelocityVerletUpdate() - "+toString(d+offset)
-        +" out of bounds [0,"+toString(Force->ColumnCounter[0])
-        +"] access on force matrix.");
-    tempVector[d] = -Force->Matrix[0][nr][d+offset]*(IsAngstroem ? AtomicLengthToAngstroem : 1.);
-  }
-  LOG(3, "INFO: Force at step " << NextStep << " set to " << tempVector);
-  setAtomicForceAtStep(NextStep, tempVector);
+  LOG(3, "INFO: Force at step " << NextStep << " is " << getAtomicForceAtStep(NextStep));
 
   // update position
-  tempVector = getPositionAtStep(NextStep-1);
-  LOG(4, "INFO: initial position from last step " << setprecision(4) << tempVector);
-  tempVector += Deltat*(getAtomicVelocityAtStep(NextStep-1));     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
-  LOG(4, "INFO: position with velocity " << getAtomicVelocityAtStep(NextStep-1) << " from last step " << tempVector);
-  tempVector += 0.5*Deltat*Deltat*(getAtomicForceAtStep(NextStep))*(1./getMass());     // F = m * a and s =
-  LOG(4, "INFO: position with force " << getAtomicForceAtStep(NextStep) << " from last step " << tempVector);
-  setPositionAtStep(NextStep, tempVector);
-  LOG(3, "INFO: Position at step " << NextStep << " set to " << tempVector);
+  {
+    Vector tempVector = getPositionAtStep(LastStep);
+    LOG(4, "INFO: initial position from last step " << setprecision(4) << tempVector);
+    tempVector += Deltat*(getAtomicVelocityAtStep(LastStep));     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
+    LOG(4, "INFO: position with velocity " << getAtomicVelocityAtStep(LastStep) << " from last step " << tempVector);
+    tempVector += .5*Deltat*Deltat*(getAtomicForceAtStep(LastStep))*(1./getMass());     // F = m * a and s =
+    LOG(4, "INFO: position with force " << getAtomicForceAtStep(NextStep) << " from last step " << tempVector);
+    setPositionAtStep(NextStep, tempVector);
+    LOG(3, "INFO: Position at step " << NextStep << " set to " << tempVector);
+  }
 
   // Update U
-  tempVector = getAtomicVelocityAtStep(NextStep-1);
-  LOG(4, "INFO: initial velocity from last step " << tempVector);
-  tempVector += Deltat * (getAtomicForceAtStep(NextStep)+getAtomicForceAtStep(NextStep-1))*(1./getMass()); // v = F/m * t
-  LOG(4, "INFO: Velocity with force from last " << getAtomicForceAtStep(NextStep-1)
-      << " and present " << getAtomicForceAtStep(NextStep) << " step " << tempVector);
-  setAtomicVelocityAtStep(NextStep, tempVector);
-  LOG(3, "INFO: Velocity at step " << NextStep << " set to " << tempVector);
+  {
+    Vector tempVector = getAtomicVelocityAtStep(LastStep);
+    LOG(4, "INFO: initial velocity from last step " << tempVector);
+    tempVector += Deltat * .5*(getAtomicForceAtStep(LastStep)+getAtomicForceAtStep(NextStep))*(1./getMass()); // v = F/m * t
+    LOG(4, "INFO: Velocity with force from last " << getAtomicForceAtStep(LastStep)
+        << " and present " << getAtomicForceAtStep(NextStep) << " step " << tempVector);
+    setAtomicVelocityAtStep(NextStep, tempVector);
+    LOG(3, "INFO: Velocity at step " << NextStep << " set to " << tempVector);
+  }
 };
 

src/Atom/atom_atominfo.hpp

@@ -267 +267 @@
   size_t getTrajectorySize() const;
   void CopyStepOnStep(const unsigned int dest, const unsigned int src);
-  void VelocityVerletUpdate(int nr, const unsigned int NextStep, double Deltat, bool IsAngstroem, ForceMatrix *Force, const size_t offset);
+  void VelocityVerletUpdate(int nr, const unsigned int NextStep, double Deltat, bool IsAngstroem);
   double getKineticEnergy(const unsigned int step) const;
   Vector getMomentum(const unsigned int step) const;

src/Dynamics/VerletForceIntegration.hpp

@@ -23 +23 @@
 #include "Fragmentation/ForceMatrix.hpp"
 #include "Helpers/helpers.hpp"
+#include "Helpers/defs.hpp"
 #include "LinearAlgebra/Vector.hpp"
 #include "Thermostats/ThermoStatContainer.hpp"
-#include "Thermostats/Berendsen.hpp"
+#include "Thermostats/Thermostat.hpp"
 #include "World.hpp"
 
@@ -32 +33 @@
 {
 public:
-  VerletForceIntegration(AtomSetMixin<T> &_atoms, double _Deltat, int _startstep, bool _IsAngstroem) :
+  /** Constructor of class VerletForceIntegration.
+   *
+   * \param _atoms set of atoms to integrate
+   * \param _Deltat time step width in atomic units
+   * \param _IsAngstroem whether length units are in angstroem or bohr radii
+   */
+  VerletForceIntegration(AtomSetMixin<T> &_atoms, double _Deltat, bool _IsAngstroem) :
     Deltat(_Deltat),
     IsAngstroem(_IsAngstroem),
-    atoms(_atoms),
-    MDSteps(_startstep)
+    atoms(_atoms)
   {}
+  /** Destructor of class VerletForceIntegration.
+   *
+   */
   ~VerletForceIntegration()
   {}
 
-  /** Parses nuclear forces from file and performs Verlet integration.
-   * Note that we assume the parsed forces to be in atomic units (hence, if coordinates are in angstroem, we
-   * have to transform them).
-   * This adds a new MD step to the config file.
+  /** Parses nuclear forces from file.
+   * Forces are stored in the time step \a TimeStep in the atomicForces in \a atoms.
    * \param *file filename
-   * \param offset offset in matrix file to the first force component
-   * \param DoConstrainedMD whether a constrained MD shall be done
-   * \param FixedCenterOfMass whether forces and velocities are correct to have fixed center of mass
-   * \return true - file found and parsed, false - file not found or imparsable
-   * \todo This is not yet checked if it is correctly working with DoConstrainedMD set >0.
-   */
-  bool operator()(char *file, const size_t offset, int DoConstrainedMD, bool FixedCenterOfMass)
+   * \param TimeStep time step to parse forces file into
+   * \return true - file parsed, false - file not found or imparsable
+   */
+  bool parseForcesFile(const char *file, const int TimeStep)
   {
     Info FunctionInfo(__func__);
-    string token;
-    stringstream item;
-    double IonMass, ActualTemp;
     ForceMatrix Force;
 
-    const int AtomCount = atoms.size();
-    ASSERT(AtomCount != 0, "VerletForceIntegration::operator() - no atoms to integrate.");
     // parse file into ForceMatrix
     std::ifstream input(file);
     if ((input.good()) && (!Force.ParseMatrix(input, 0,0,0))) {
       ELOG(0, "Could not parse Force Matrix file " << file << ".");
-      performCriticalExit();
       return false;
     }
     input.close();
-    if (Force.RowCounter[0] != AtomCount) {
+    if (Force.RowCounter[0] != (int)atoms.size()) {
       ELOG(0, "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << atoms.size() << ".");
-      performCriticalExit();
       return false;
     }
 
-    if (FixedCenterOfMass) {
-      Vector ForceVector;
-      // correct Forces
-      //std::cout << "Force before correction, " << Force << std::endl;
-      ForceVector.Zero();
-      for(int i=0;i<AtomCount;i++)
-        for(int d=0;d<NDIM;d++) {
-          ForceVector[d] += Force.Matrix[0][i][d+offset];
-        }
-      ForceVector.Scale(1./static_cast<double>(AtomCount));
-      //std::cout << "Force before second correction, " << Force << std::endl;
-      for(int i=0;i<AtomCount;i++)
-        for(int d=0;d<NDIM;d++) {
-          Force.Matrix[0][i][d+offset] -= ForceVector[d];
-        }
-      LOG(3, "INFO: forces correct by " << ForceVector << "each.");
-    }
+    addForceMatrixToAtomicForce(Force, TimeStep, 1);
+    return true;
+  }
+
+  /** Performs Verlet integration.
+   * Note that we assume the parsed forces to be in atomic units (hence, if coordinates are in angstroem, we
+   * have to transform them).
+   * This adds a new MD step \f$ t + \Delta t \f$ to the config file.
+   * \param NextStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
+   * \param offset offset in matrix file to the first force component
+   * \param DoConstrainedMD whether a constrained MD shall be done
+   * \param FixedCenterOfMass whether forces and velocities are correct to have fixed center of mass
+   * \return true - file found and parsed, false - no atoms, file not found or imparsable
+   * \todo This is not yet checked if it is correctly working with DoConstrainedMD set >0.
+   */
+  bool operator()(const int NextStep, const size_t offset, const int DoConstrainedMD, const bool FixedCenterOfMass)
+  {
+    Info FunctionInfo(__func__);
+
+    // check that atoms are present at all
+    if (atoms.size() == 0) {
+      ELOG(2, "VerletForceIntegration::operator() - no atoms to integrate.");
+      return false;
+    }
+
+    // make sum of forces equal zero
+    if (FixedCenterOfMass)
+      correctForceMatrixForFixedCenterOfMass(offset,NextStep);
 
     // solve a constrained potential if we are meant to
     if (DoConstrainedMD) {
-      // calculate forces and potential
-      std::map<atom *, atom*> PermutationMap;
-      MinimiseConstrainedPotential Minimiser(atoms, PermutationMap);
-      //double ConstrainedPotentialEnergy =
-      Minimiser(DoConstrainedMD, 0, IsAngstroem);
-      Minimiser.EvaluateConstrainedForces(&Force);
+      performConstraintMinimization(DoConstrainedMD,NextStep);
     }
 
@@ -109 +112 @@
     for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
       //std::cout << "Id of atom is " << (*iter)->getId() << std::endl;
-      (*iter)->VelocityVerletUpdate((*iter)->getId(), MDSteps+1, Deltat, IsAngstroem, &Force, (const size_t) 1);
-    }
-
-    if (FixedCenterOfMass) {
-      Vector Velocity;
-      // correct velocities (rather momenta) so that center of mass remains motionless
-      Velocity = atoms.totalMomentumAtStep(MDSteps+1);
-      IonMass = atoms.totalMass();
-
-      // correct velocities (rather momenta) so that center of mass remains motionless
-      Velocity *= 1./IonMass;
-      atoms.addVelocityAtStep(-1.*Velocity,MDSteps+1);
-
-      LOG(3, "INFO: Velocities corrected by " << Velocity << " each.");
-    }
+      (*iter)->VelocityVerletUpdate((*iter)->getId(), NextStep, Deltat, IsAngstroem);
+    }
+
+    // make sum of velocities equal zero
+    if (FixedCenterOfMass)
+      correctVelocitiesForFixedCenterOfMass(NextStep);
 
     // thermostat
-    ActualTemp = atoms.totalTemperatureAtStep(MDSteps+1);
-    LOG(3, "INFO: Current temperature is " << ActualTemp);
-    Berendsen berendsen = Berendsen();
-    berendsen.addToContainer(World::getInstance().getThermostats());
-    double ekin = berendsen.scaleAtoms(MDSteps,ActualTemp,atoms);
-    ActualTemp = atoms.totalTemperatureAtStep(MDSteps+1);
-    LOG(3, "INFO: New temperature after thermostat is " << ActualTemp);
-    LOG(1, "Kinetic energy is " << ekin << ".");
-
-    // next step
-    MDSteps++;
+    performThermostatControl(NextStep);
 
     // exit
     return true;
   };
+
+private:
+  void addForceMatrixToAtomicForce(const ForceMatrix &Force, const int &TimeStep, const int offset)
+  {
+    // place forces from matrix into atoms
+    Vector tempVector;
+    size_t i=0;
+    for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter,++i) {
+      for(size_t d=0;d<NDIM;d++) {
+        tempVector[d] = Force.Matrix[0][i][d+offset]*(IsAngstroem ? AtomicLengthToAngstroem : 1.);
+      }
+      tempVector += (*iter)->getAtomicForceAtStep(TimeStep);
+      (*iter)->setAtomicForceAtStep(TimeStep, tempVector);
+    }
+  }
+
+  void correctForceMatrixForFixedCenterOfMass(const size_t offset, const int &TimeStep) {
+    Vector ForceVector;
+    // correct Forces
+    //std::cout << "Force before correction, " << Force << std::endl;
+    ForceVector.Zero();
+    for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
+      ForceVector += (*iter)->getAtomicForceAtStep(TimeStep);
+    }
+    ForceVector.Scale(1./(double)atoms.size());
+    //std::cout << "Force before second correction, " << Force << std::endl;
+    for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
+      const Vector tempVector = (*iter)->getAtomicForceAtStep(TimeStep) - ForceVector;
+      (*iter)->setAtomicForceAtStep(TimeStep, tempVector);
+    }
+    LOG(3, "INFO: forces correct by " << ForceVector << "each.");
+  }
+
+  void correctVelocitiesForFixedCenterOfMass(const int &TimeStep) {
+    Vector Velocity;
+    double IonMass;
+    // correct velocities (rather momenta) so that center of mass remains motionless
+    Velocity = atoms.totalMomentumAtStep(TimeStep);
+    IonMass = atoms.totalMass();
+
+    // correct velocities (rather momenta) so that center of mass remains motionless
+    Velocity *= 1./IonMass;
+    atoms.addVelocityAtStep(-1.*Velocity,TimeStep);
+
+    LOG(3, "INFO: Velocities corrected by " << Velocity << " each.");
+  }
+
+  void performConstraintMinimization(const int &DoConstrainedMD, const int &TimeStep) {
+    // calculate forces and potential
+    ForceMatrix Force;
+    std::map<atom *, atom*> PermutationMap;
+    MinimiseConstrainedPotential Minimiser(atoms, PermutationMap);
+    //double ConstrainedPotentialEnergy =
+    Minimiser(DoConstrainedMD, 0, IsAngstroem);
+    Minimiser.EvaluateConstrainedForces(&Force);
+    addForceMatrixToAtomicForce(Force, TimeStep, 1);
+  }
+
+  void performThermostatControl(const int &TimeStep) {
+    double ActualTemp;
+
+    // calculate current temperature
+    ActualTemp = atoms.totalTemperatureAtStep(TimeStep);
+    LOG(3, "INFO: Current temperature is " << ActualTemp);
+
+    // rescale to desired value
+    double ekin = World::getInstance().getThermostats()->getActive()->scaleAtoms(TimeStep,ActualTemp,atoms);
+    ActualTemp = atoms.totalTemperatureAtStep(TimeStep);
+    LOG(3, "INFO: New temperature after thermostat is " << ActualTemp);
+    LOG(1, "Kinetic energy is " << ekin << ".");
+  }
 
 private:
@@ -146 +201 @@
   bool IsAngstroem;
   AtomSetMixin<T> atoms;
-  int MDSteps;
 };
 

src/Thermostats/ThermoStatContainer.hpp

@@ -25 +25 @@
   Thermostat *makeByName(const std::string name,class ConfigFileBuffer * const fb);
   void makeActive(const std::string name,class ConfigFileBuffer * const fb);
+  Thermostat *getActive() {
+    return activeThermostat;
+  }
   void chooseNone();