source: src/Dynamics/VerletForceIntegration.hpp@ 41a467

/*
 * 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 "Fragmentation/ForceMatrix.hpp"
#include "Helpers/helpers.hpp"
#include "LinearAlgebra/Vector.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))) {
      ELOG(0, "Could not parse Force Matrix file " << file << ".");
      performCriticalExit();
      return false;
    }
    input.close();
    if (Force.RowCounter[0] != AtomCount) {
      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.");
    }

    // 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);
    LOG(1, "Kinetic energy is " << ekin << ".");

    // next step
    MDSteps++;

    // exit
    return true;
  };

private:
  double Deltat;
  bool IsAngstroem;
  AtomSetMixin<T> atoms;
  int MDSteps;
};

#endif /* VERLETFORCEINTEGRATION_HPP_ */