Index: src/Actions/MoleculeAction/VerletIntegrationAction.cpp
===================================================================
--- src/Actions/MoleculeAction/VerletIntegrationAction.cpp	(revision eb33c4a913ab83436dbf30ab934797ec2d623036)
+++ src/Actions/MoleculeAction/VerletIntegrationAction.cpp	(revision 435065d3d9c23062fca0bb67ff882fade26acced)
@@ -20,9 +20,13 @@
 #include "CodePatterns/MemDebug.hpp"
 
+#include "atom.hpp"
+#include "AtomSet.hpp"
 #include "CodePatterns/Log.hpp"
+#include "CodePatterns/Verbose.hpp"
+#include "Dynamics/VerletForceIntegration.hpp"
 #include "molecule.hpp"
-#include "CodePatterns/Verbose.hpp"
 #include "World.hpp"
 
+#include <vector>
 #include <iostream>
 #include <fstream>
@@ -38,20 +42,21 @@
 /** =========== define the function ====================== */
 Action::state_ptr MoleculeVerletIntegrationAction::performCall() {
-  molecule *mol = NULL;
-
   // obtain information
   getParametersfromValueStorage();
 
-  for (World::MoleculeSelectionIterator iter = World::getInstance().beginMoleculeSelection(); iter != World::getInstance().endMoleculeSelection(); ++iter) {
-    mol = iter->second;
-    DoLog(1) && (Log() << Verbose(1) << "Parsing forces file and Verlet integrating." << endl);
-    // TODO: sollte besser stream nutzen, nicht filename direkt (es sei denn, ist prefix), besser fuer unit test
-    char outputname[MAXSTRINGSIZE];
-    strcpy(outputname, params.forcesfile.string().c_str());
-    if (!mol->VerletForceIntegration(outputname, *(World::getInstance().getConfig()), 0))
+
+  DoLog(1) && (Log() << Verbose(1) << "Parsing forces file and Verlet integrating." << endl);
+  // TODO: sollte besser stream nutzen, nicht filename direkt (es sei denn, ist prefix), besser fuer unit test
+  char outputname[MAXSTRINGSIZE];
+  strcpy(outputname, params.forcesfile.string().c_str());
+  AtomSetMixin<std::vector<atom *> > set(World::getInstance().getSelectedAtoms());
+  for (int step = 0; step < params.MDSteps; ++step) {
+    VerletForceIntegration<std::vector<atom *> > Verlet(set, params.Deltat, step, false);
+    if (!Verlet(outputname, 1, 0, params.FixedCenterOfMass))
       DoLog(2) && (Log() << Verbose(2) << "File " << params.forcesfile << " not found." << endl);
     else
       DoLog(2) && (Log() << Verbose(2) << "File " << params.forcesfile << " found and parsed." << endl);
   }
+
   return Action::success;
 }
Index: src/Actions/MoleculeAction/VerletIntegrationAction.def
===================================================================
--- src/Actions/MoleculeAction/VerletIntegrationAction.def	(revision eb33c4a913ab83436dbf30ab934797ec2d623036)
+++ src/Actions/MoleculeAction/VerletIntegrationAction.def	(revision 435065d3d9c23062fca0bb67ff882fade26acced)
@@ -12,9 +12,9 @@
 // ValueStorage by the token "Z" -> first column: int, Z, "Z"
 // "undefine" if no parameters are required, use (NODEFAULT) for each (undefined) default value
-#define paramtypes (boost::filesystem::path)
-#define paramtokens ("verlet-integration")
-#define paramdescriptions ("perform verlet integration of a given force file")
+#define paramtypes (boost::filesystem::path)(double)(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 paramreferences (forcesfile)
+#define paramreferences (forcesfile)(Deltat)(MDSteps)(FixedCenterOfMass)
 
 #undef statetypes
Index: src/Dynamics/VerletForceIntegration.hpp
===================================================================
--- src/Dynamics/VerletForceIntegration.hpp	(revision 435065d3d9c23062fca0bb67ff882fade26acced)
+++ src/Dynamics/VerletForceIntegration.hpp	(revision 435065d3d9c23062fca0bb67ff882fade26acced)
@@ -0,0 +1,152 @@
+/*
+ * VerletForceIntegration.hpp
+ *
+ *  Created on: Feb 23, 2011
+ *      Author: heber
+ */
+
+#ifndef VERLETFORCEINTEGRATION_HPP_
+#define VERLETFORCEINTEGRATION_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "atom.hpp"
+#include "AtomSet.hpp"
+#include "CodePatterns/Assert.hpp"
+#include "CodePatterns/Info.hpp"
+#include "CodePatterns/Log.hpp"
+#include "CodePatterns/Verbose.hpp"
+#include "Dynamics/MinimiseConstrainedPotential.hpp"
+#include "LinearAlgebra/Vector.hpp"
+#include "parser.hpp"
+#include "ThermoStatContainer.hpp"
+#include "Thermostats/Berendsen.hpp"
+#include "World.hpp"
+
+template <class T>
+class VerletForceIntegration
+{
+public:
+  VerletForceIntegration(AtomSetMixin<T> &_atoms, double _Deltat, int _startstep, bool _IsAngstroem) :
+    Deltat(_Deltat),
+    IsAngstroem(_IsAngstroem),
+    atoms(_atoms),
+    MDSteps(_startstep)
+  {}
+  ~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.
+   * \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)
+  {
+    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))) {
+      DoeLog(0) && (eLog()<< Verbose(0) << "Could not parse Force Matrix file " << file << "." << endl);
+      performCriticalExit();
+      return false;
+    }
+    input.close();
+    if (Force.RowCounter[0] != AtomCount) {
+      DoeLog(0) && (eLog()<< Verbose(0) << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << atoms.size() << "." << endl);
+      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.");
+    }
+
+    // solve a constrained potential if we are meant to
+    if (DoConstrainedMD) {
+      // calculate forces and potential
+      std::map<atom *, atom*> PermutationMap;
+      molecule::atomSet atoms_list;
+      copy(atoms.begin(), atoms.end(), atoms_list.begin());
+      MinimiseConstrainedPotential Minimiser(atoms_list, PermutationMap);
+      //double ConstrainedPotentialEnergy =
+      Minimiser(DoConstrainedMD, 0, IsAngstroem);
+      Minimiser.EvaluateConstrainedForces(&Force);
+    }
+
+    //std::cout << "Force before velocity verlet, " << Force << std::endl;
+    // and perform Verlet integration for each atom with position, velocity and force vector
+    // check size of vectors
+    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.");
+    }
+
+    // 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);
+    DoLog(1) && (Log() << Verbose(1) << "Kinetic energy is " << ekin << "." << endl);
+
+    // next step
+    MDSteps++;
+
+    // exit
+    return true;
+  };
+
+private:
+  double Deltat;
+  bool IsAngstroem;
+  AtomSetMixin<T> atoms;
+  int MDSteps;
+};
+
+#endif /* VERLETFORCEINTEGRATION_HPP_ */
Index: src/Makefile.am
===================================================================
--- src/Makefile.am	(revision eb33c4a913ab83436dbf30ab934797ec2d623036)
+++ src/Makefile.am	(revision 435065d3d9c23062fca0bb67ff882fade26acced)
@@ -125,5 +125,6 @@
 	Dynamics/LinearInterpolationBetweenSteps.hpp \
 	Dynamics/MinimiseConstrainedPotential.hpp \
-	Dynamics/OutputTemperature.hpp
+	Dynamics/OutputTemperature.hpp \
+	Dynamics/VerletForceIntegration.hpp
 
 GRAPHSOURCE = \
Index: src/atom_atominfo.cpp
===================================================================
--- src/atom_atominfo.cpp	(revision eb33c4a913ab83436dbf30ab934797ec2d623036)
+++ src/atom_atominfo.cpp	(revision 435065d3d9c23062fca0bb67ff882fade26acced)
@@ -30,4 +30,6 @@
 #include "atom_atominfo.hpp"
 
+#include <iomanip>
+
 /** Constructor of class AtomInfo.
  */
@@ -492,32 +494,48 @@
  * Parameters are according to those in configuration class.
  * \param NextStep index of sequential step to set
- * \param *configuration pointer to configuration with parameters
+ * \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, config *configuration, ForceMatrix *Force, const size_t offset)
-{
+void AtomInfo::VelocityVerletUpdate(int nr, const unsigned int NextStep, double Deltat, bool IsAngstroem, ForceMatrix *Force, const size_t offset)
+{
+  LOG(2, "INFO: Particle that currently " << *this);
+  LOG(2, "INFO: Integrating with mass=" << getMass() << " and Deltat="
+      << Deltat << " at NextStep=" << NextStep);
   // 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++)
-    tempVector[d] = -Force->Matrix[0][nr][d+offset]*(configuration->GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
+  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+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);
 
   // update position
   tempVector = getPositionAtStep(NextStep-1);
-  tempVector += configuration->Deltat*(getAtomicVelocityAtStep(NextStep-1));     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
-  tempVector += 0.5*configuration->Deltat*configuration->Deltat*(getAtomicForceAtStep(NextStep))*(1./getMass());     // F = m * a and s =
+  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);
 
   // Update U
   tempVector = getAtomicVelocityAtStep(NextStep-1);
-  tempVector += configuration->Deltat * (getAtomicForceAtStep(NextStep)+getAtomicForceAtStep(NextStep-1))*(1./getMass()); // v = F/m * t
+  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);
-
-  // some info for debugging
-  DoLog(2) && (Log() << Verbose(2)
-      << "Integrated position&velocity of step " << (NextStep) << ": ("
-      << getPositionAtStep(NextStep) << ")\t("
-      << getAtomicVelocityAtStep(NextStep) << ")" << std::endl);
+  LOG(3, "INFO: Velocity at step " << NextStep << " set to " << tempVector);
 };
 
Index: src/atom_atominfo.hpp
===================================================================
--- src/atom_atominfo.hpp	(revision eb33c4a913ab83436dbf30ab934797ec2d623036)
+++ src/atom_atominfo.hpp	(revision 435065d3d9c23062fca0bb67ff882fade26acced)
@@ -27,5 +27,4 @@
 
 class AtomInfo;
-class config;
 class element;
 class ForceMatrix;
@@ -242,5 +241,5 @@
   size_t getTrajectorySize() const;
   void CopyStepOnStep(const unsigned int dest, const unsigned int src);
-  void VelocityVerletUpdate(int nr, const unsigned int MDSteps, config *configuration, ForceMatrix *Force, const size_t offset);
+  void VelocityVerletUpdate(int nr, const unsigned int NextStep, double Deltat, bool IsAngstroem, ForceMatrix *Force, const size_t offset);
   double getKineticEnergy(const unsigned int step) const;
   Vector getMomentum(const unsigned int step) const;
Index: src/molecule.hpp
===================================================================
--- src/molecule.hpp	(revision eb33c4a913ab83436dbf30ab934797ec2d623036)
+++ src/molecule.hpp	(revision 435065d3d9c23062fca0bb67ff882fade26acced)
@@ -176,5 +176,4 @@
   void SetBoxDimension(Vector *dim);
   bool ScanForPeriodicCorrection();
-  bool VerletForceIntegration(char *file, config &configuration, const size_t offset);
   double VolumeOfConvexEnvelope(bool IsAngstroem);
   RealSpaceMatrix getInertiaTensor() const;
Index: src/molecule_dynamics.cpp
===================================================================
--- src/molecule_dynamics.cpp	(revision eb33c4a913ab83436dbf30ab934797ec2d623036)
+++ src/molecule_dynamics.cpp	(revision 435065d3d9c23062fca0bb67ff882fade26acced)
@@ -39,79 +39,2 @@
 
 
-/** 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.
- * \param *file filename
- * \param config structure with config::Deltat, config::IsAngstroem, config::DoConstrained
- * \param offset offset in matrix file to the first force component
- * \return true - file found and parsed, false - file not found or imparsable
- * \todo This is not yet checked if it is correctly working with DoConstrained set to true.
- */
-bool molecule::VerletForceIntegration(char *file, config &configuration, const size_t offset)
-{
-  Info FunctionInfo(__func__);
-  string token;
-  stringstream item;
-  double IonMass, ConstrainedPotentialEnergy, ActualTemp;
-  Vector Velocity;
-  ForceMatrix Force;
-
-  const int AtomCount = getAtomCount();
-  // parse file into ForceMatrix
-  std::ifstream input(file);
-  if ((input.good()) && (!Force.ParseMatrix(input, 0,0,0))) {
-    DoeLog(0) && (eLog()<< Verbose(0) << "Could not parse Force Matrix file " << file << "." << endl);
-    performCriticalExit();
-    return false;
-  }
-  input.close();
-  if (Force.RowCounter[0] != AtomCount) {
-    DoeLog(0) && (eLog()<< Verbose(0) << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << getAtomCount() << "." << endl);
-    performCriticalExit();
-    return false;
-  }
-  // correct Forces
-  Velocity.Zero();
-  for(int i=0;i<AtomCount;i++)
-    for(int d=0;d<NDIM;d++) {
-      Velocity[d] += Force.Matrix[0][i][d+offset];
-    }
-  for(int i=0;i<AtomCount;i++)
-    for(int d=0;d<NDIM;d++) {
-      Force.Matrix[0][i][d+offset] -= Velocity[d]/static_cast<double>(AtomCount);
-    }
-  // solve a constrained potential if we are meant to
-  if (configuration.DoConstrainedMD) {
-    // calculate forces and potential
-    std::map<atom*, atom *> PermutationMap;
-    MinimiseConstrainedPotential Minimiser(atoms, PermutationMap);
-    //double ConstrainedPotentialEnergy =
-    Minimiser(configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
-    Minimiser.EvaluateConstrainedForces(&Force);
-  }
-
-  // and perform Verlet integration for each atom with position, velocity and force vector
-  // check size of vectors
-  BOOST_FOREACH(atom *_atom, atoms) {
-    _atom->VelocityVerletUpdate(_atom->getNr(), MDSteps+1, &configuration, &Force, (const size_t) 0);
-  }
-
-  // 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);
-  ActualTemp = atoms.totalTemperatureAtStep(MDSteps+1);
-  Berendsen berendsen = Berendsen();
-  berendsen.addToContainer(World::getInstance().getThermostats());
-  double ekin = berendsen.scaleAtoms(MDSteps,ActualTemp,atoms);
-  DoLog(1) && (Log() << Verbose(1) << "Kinetic energy is " << ekin << "." << endl);
-  MDSteps++;
-
-  // exit
-  return true;
-};
