/*
 * 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/atom.hpp"
#include "Atom/AtomSet.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "Dynamics/AtomicForceManipulator.hpp"
#include "Fragmentation/ForceMatrix.hpp"
#include "Helpers/helpers.hpp"
#include "Helpers/defs.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "Thermostats/ThermoStatContainer.hpp"
#include "Thermostats/Thermostat.hpp"
#include "World.hpp"

template <class T>
class VerletForceIntegration : public AtomicForceManipulator<T>
{
public:
  /** 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) :
    AtomicForceManipulator<T>(_atoms, _Deltat, _IsAngstroem)
  {}
  /** Destructor of class VerletForceIntegration.
   *
   */
  ~VerletForceIntegration()
  {}

  /** Performs Verlet integration.
   *
   * We assume that forces have just been calculated. Then, we perform the velocity
   * and the position calculation for \f$ t + \Delta t \f$, such that forces may be
   * again calculated with respect to the new position.
   *
   * \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
   * \todo This is not yet checked if it is correctly working with DoConstrainedMD set >0.
   */
  void operator()(const int NextStep, const size_t offset, const int DoConstrainedMD, const bool FixedCenterOfMass)
  {
    Info FunctionInfo(__func__);

    // make sum of forces equal zero
    if (FixedCenterOfMass)
      AtomicForceManipulator<T>::correctForceMatrixForFixedCenterOfMass(offset,NextStep-1);

    // solve a constrained potential if we are meant to
//    if (DoConstrainedMD)
//      performConstraintMinimization(DoConstrainedMD,NextStep-1);

    if (NextStep > 0) {
      for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
          iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
        //std::cout << "Id of atom is " << (*iter)->getId() << std::endl;
        (*iter)->VelocityVerletUpdateU(
            (*iter)->getId(),
            NextStep-1,
            AtomicForceManipulator<T>::Deltat,
            AtomicForceManipulator<T>::IsAngstroem);
      }

      // make sum of velocities equal zero
      if (FixedCenterOfMass)
        AtomicForceManipulator<T>::correctVelocitiesForFixedCenterOfMass(NextStep-1);

      // thermostat
      AtomicForceManipulator<T>::performThermostatControl(NextStep-1);
    }

    //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 = AtomicForceManipulator<T>::atoms.begin();
        iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
      //std::cout << "Id of atom is " << (*iter)->getId() << std::endl;
      (*iter)->VelocityVerletUpdateX(
          (*iter)->getId(),
          NextStep,
          AtomicForceManipulator<T>::Deltat,
          AtomicForceManipulator<T>::IsAngstroem);
    }
  }
};

#endif /* VERLETFORCEINTEGRATION_HPP_ */