source: src/molecule_dynamics.cpp@ 20943b

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

/*
 * molecule_dynamics.cpp
 *
 *  Created on: Oct 5, 2009
 *      Author: heber
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "CodePatterns/MemDebug.hpp"

#include "World.hpp"
#include "atom.hpp"
#include "config.hpp"
#include "Dynamics/MinimiseConstrainedPotential.hpp"
//#include "element.hpp"
#include "CodePatterns/Info.hpp"
//#include "CodePatterns/Verbose.hpp"
//#include "CodePatterns/Log.hpp"
#include "molecule.hpp"
#include "parser.hpp"
//#include "LinearAlgebra/Plane.hpp"
#include "ThermoStatContainer.hpp"
#include "Thermostats/Berendsen.hpp"

#include "CodePatterns/enumeration.hpp"

/************************************* Functions for class molecule *********************************/


/** 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
    atom **PermutationMap = NULL;
    MinimiseConstrainedPotential Minimiser(atoms);
    //double ConstrainedPotentialEnergy =
    Minimiser(PermutationMap, configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
    Minimiser.EvaluateConstrainedForces(&Force);
    delete[](PermutationMap);
  }

  // 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;
};