source: src/Atom/atom_atominfo.cpp@ 4464ef

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * Copyright (C)  2014 Frederik Heber. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * atom_atominfo.cpp
 *
 *  Created on: Oct 19, 2009
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include "CodePatterns/Verbose.hpp"

#include "atom_atominfo.hpp"
#include "CodePatterns/Log.hpp"
#include "config.hpp"
#include "Element/element.hpp"
#include "Element/periodentafel.hpp"
#include "Fragmentation/ForceMatrix.hpp"
#include "World.hpp"
#include "WorldTime.hpp"

#include <iomanip>

/** Constructor of class AtomInfo.
 */
AtomInfo::AtomInfo() :
  AtomicElement(1),
  FixedIon(false),
  charge(0.)
{
  AtomicPosition.insert( std::make_pair(0, zeroVec) );
  AtomicVelocity.insert( std::make_pair(0, zeroVec) );
  AtomicForce.insert( std::make_pair(0, zeroVec) );
}

/** Copy constructor of class AtomInfo.
 */
AtomInfo::AtomInfo(const AtomInfo &_atom) :
    AtomicPosition(_atom.AtomicPosition),
    AtomicVelocity(_atom.AtomicVelocity),
    AtomicForce(_atom.AtomicForce),
    AtomicElement(_atom.AtomicElement),
    FixedIon(_atom.FixedIon),
    charge(_atom.charge)
{
}

AtomInfo::AtomInfo(const VectorInterface &_v) :
    AtomicElement(1),
    FixedIon(false),
    charge(0.)
{
  AtomicPosition.insert( std::make_pair(0, _v.getPosition()) );
  AtomicVelocity.insert( std::make_pair(0, zeroVec) );
  AtomicForce.insert( std::make_pair(0, zeroVec) );
};

/** Destructor of class AtomInfo.
 */
AtomInfo::~AtomInfo()
{
};

void AtomInfo::AppendTrajectoryStep(const unsigned int _step)
{
  ASSERT (WorldTime::getTime() != _step,
      "AtomInfo::AppendTrajectoryStep() - cannot append current time step.");
  NOTIFY(TrajectoryChanged);
  AtomicPosition.insert( std::make_pair(_step, zeroVec) );
  AtomicVelocity.insert( std::make_pair(_step, zeroVec) );
  AtomicForce.insert( std::make_pair(_step, zeroVec) );
  LOG(5,"AtomInfo::AppendTrajectoryStep() called, size is ("
      << AtomicPosition.size() << ","
      << AtomicVelocity.size() << ","
      << AtomicForce.size() << ")");
}

void AtomInfo::removeTrajectoryStep(const unsigned int _step)
{
  ASSERT (WorldTime::getTime() != _step,
      "AtomInfo::removeTrajectoryStep() - cannot remove current time step.");
  NOTIFY(TrajectoryChanged);
  AtomicPosition.erase(_step);
  AtomicVelocity.erase(_step);
  AtomicForce.erase(_step);
  LOG(5,"AtomInfo::removeTrajectoryStep() called, size is ("
      << AtomicPosition.size() << ","
      << AtomicVelocity.size() << ","
      << AtomicForce.size() << ")");
}

const element *AtomInfo::getType() const
{
  const element *elem = World::getInstance().getPeriode()->FindElement(AtomicElement);
  return elem;
}

const element &AtomInfo::getElement() const
{
  const element &elem = *World::getInstance().getPeriode()->FindElement(AtomicElement);
  return elem;
}

atomicNumber_t AtomInfo::getElementNo() const
{
  return AtomicElement;
}

const std::string &AtomInfo::getParticleName() const
{
  return particlename;
}

void AtomInfo::setParticleName(const std::string & _name)
{
  particlename = _name;
}

const double& AtomInfo::operator[](size_t i) const
{
  return atStep(i, WorldTime::getTime());
}

const double& AtomInfo::at(size_t i) const
{
  return atStep(i, WorldTime::getTime());
}

const double& AtomInfo::atStep(size_t i, unsigned int _step) const
{
  ASSERT(!AtomicPosition.empty(),
      "AtomInfo::operator[]() - AtomicPosition is empty.");
  VectorTrajectory_t::const_iterator iter =
      AtomicPosition.lower_bound(_step);
  return iter->second[i];
}

void AtomInfo::set(size_t i, const double value)
{
  OBSERVE;
  NOTIFY(AtomObservable::PositionChanged);
  VectorTrajectory_t::iterator iter = AtomicPosition.find(WorldTime::getTime());
  if (iter !=  AtomicPosition.end()) {
    iter->second[i] = value;
  } else {
    Vector newPos;
    newPos[i] = value;
#ifndef NDEBUG
    std::pair<VectorTrajectory_t::iterator, bool> inserter =
#endif
        AtomicPosition.insert( std::make_pair(WorldTime::getTime(), newPos) );
    ASSERT( inserter.second,
        "AtomInfo::set() - time step "+toString(WorldTime::getTime())
        +" present after all?");
  }
}

/** Helps to determine whether the current step really exists or getPosition() has just
 * delivered the one closest to it in the past.
 *
 * \param _step step to check
 * \param true - step exists, false - step does not exist, getPosition() delivers closest
 */
bool AtomInfo::isStepPresent(const unsigned int _step) const
{
  VectorTrajectory_t::const_iterator iter =
      AtomicPosition.find(_step);
  return iter != AtomicPosition.end();
}

const Vector& AtomInfo::getPosition() const
{
  return getPositionAtStep(WorldTime::getTime());
}

const Vector& AtomInfo::getPositionAtStep(const unsigned int _step) const
{
  ASSERT(!AtomicPosition.empty(),
      "AtomInfo::operator[]() - AtomicPosition is empty.");
  VectorTrajectory_t::const_iterator iter =
      AtomicPosition.lower_bound(_step);
  return iter->second;
}

void AtomInfo::setType(const element* _type)
{
  OBSERVE;
  NOTIFY(AtomObservable::ElementChanged);
  AtomicElement = _type->getAtomicNumber();
}

void AtomInfo::setType(const int Z)
{
  const element *elem = World::getInstance().getPeriode()->FindElement(Z);
  setType(elem);
}

const Vector& AtomInfo::getAtomicVelocity() const
{
  return getAtomicVelocityAtStep(WorldTime::getTime());
}

const Vector& AtomInfo::getAtomicVelocityAtStep(const unsigned int _step) const
{
  ASSERT(!AtomicVelocity.empty(),
      "AtomInfo::operator[]() - AtomicVelocity is empty.");
  VectorTrajectory_t::const_iterator iter =
      AtomicVelocity.lower_bound(_step);
  // special, we only interpolate between present time steps not into the future
  if (_step > AtomicVelocity.begin()->first)
    return zeroVec;
  else
    return iter->second;
}

void AtomInfo::setAtomicVelocity(const Vector &_newvelocity)
{
  setAtomicVelocityAtStep(WorldTime::getTime(), _newvelocity);
}

void AtomInfo::setAtomicVelocityAtStep(const unsigned int _step, const Vector &_newvelocity)
{
  OBSERVE;
  VectorTrajectory_t::iterator iter = AtomicVelocity.find(_step);
  if (iter !=  AtomicVelocity.end()) {
    iter->second = _newvelocity;
  } else {
#ifndef NDEBUG
    std::pair<VectorTrajectory_t::iterator, bool> inserter =
#endif
        AtomicVelocity.insert( std::make_pair(_step, _newvelocity) );
    ASSERT( inserter.second,
        "AtomInfo::set() - time step "+toString(_step)
        +" present after all?");
  }
  if (WorldTime::getTime() == _step)
    NOTIFY(AtomObservable::VelocityChanged);
}

const Vector& AtomInfo::getAtomicForce() const
{
  return getAtomicForceAtStep(WorldTime::getTime());
}

const Vector& AtomInfo::getAtomicForceAtStep(const unsigned int _step) const
{
  ASSERT(!AtomicForce.empty(),
      "AtomInfo::operator[]() - AtomicForce is empty.");
  VectorTrajectory_t::const_iterator iter =
      AtomicForce.lower_bound(_step);
  // special, we only interpolate between present time steps not into the future
  if (_step > AtomicForce.begin()->first)
    return zeroVec;
  else
    return iter->second;
}

void AtomInfo::setAtomicForce(const Vector &_newforce)
{
  setAtomicForceAtStep(WorldTime::getTime(), _newforce);
}

void AtomInfo::setAtomicForceAtStep(const unsigned int _step, const Vector &_newforce)
{
  OBSERVE;
  VectorTrajectory_t::iterator iter = AtomicForce.find(_step);
  if (iter !=  AtomicForce.end()) {
    iter->second = _newforce;
  } else {
#ifndef NDEBUG
    std::pair<VectorTrajectory_t::iterator, bool> inserter =
#endif
        AtomicForce.insert( std::make_pair(_step, _newforce) );
    ASSERT( inserter.second,
        "AtomInfo::set() - time step "+toString(_step)
        +" present after all?");
  }
  if (WorldTime::getTime() == _step)
    NOTIFY(AtomObservable::ForceChanged);
}

bool AtomInfo::getFixedIon() const
{
  return FixedIon;
}

void AtomInfo::setFixedIon(const bool _fixedion)
{
  OBSERVE;
  NOTIFY(AtomObservable::PropertyChanged);
  FixedIon = _fixedion;
}

void AtomInfo::setPosition(const Vector& _vector)
{
  setPositionAtStep(WorldTime::getTime(), _vector);
}

void AtomInfo::setPositionAtStep(unsigned int _step, const Vector& _vector)
{
  OBSERVE;
  VectorTrajectory_t::iterator iter = AtomicPosition.find(_step);
  if (iter !=  AtomicPosition.end()) {
    iter->second = _vector;
  } else {
#ifndef NDEBUG
    std::pair<VectorTrajectory_t::iterator, bool> inserter =
#endif
        AtomicPosition.insert( std::make_pair(_step, _vector) );
    ASSERT( inserter.second,
        "AtomInfo::set() - time step "+toString(_step)
        +" present after all?");
  }
  if (WorldTime::getTime() == _step)
    NOTIFY(AtomObservable::PositionChanged);
}

const VectorInterface& AtomInfo::operator+=(const Vector& b)
{
  setPosition(getPosition()+b);
  return *this;
}

const VectorInterface& AtomInfo::operator-=(const Vector& b)
{
  setPosition(getPosition()-b);
  return *this;
}

Vector const AtomInfo::operator+(const Vector& b) const
{
  Vector a(getPosition());
  a += b;
  return a;
}

Vector const AtomInfo::operator-(const Vector& b) const
{
  Vector a(getPosition());
  a -= b;
  return a;
}

double AtomInfo::distance(const Vector &point) const
{
  return getPosition().distance(point);
}

double AtomInfo::DistanceSquared(const Vector &y) const
{
  return getPosition().DistanceSquared(y);
}

double AtomInfo::distance(const VectorInterface &_atom) const
{
  return _atom.distance(getPosition());
}

double AtomInfo::DistanceSquared(const VectorInterface &_atom) const
{
  return _atom.DistanceSquared(getPosition());
}

VectorInterface &AtomInfo::operator=(const Vector& _vector)
{
  setPosition(_vector);
  return *this;
}

void AtomInfo::ScaleAll(const double *factor)
{
  Vector temp(getPosition());
  temp.ScaleAll(factor);
  setPosition(temp);
}

void AtomInfo::ScaleAll(const Vector &factor)
{
  Vector temp(getPosition());
  temp.ScaleAll(factor);
  setPosition(temp);
}

void AtomInfo::Scale(const double factor)
{
  Vector temp(getPosition());
  temp.Scale(factor);
  setPosition(temp);
}

void AtomInfo::Zero()
{
  setPosition(zeroVec);
}

void AtomInfo::One(const double one)
{
  setPosition(Vector(one,one,one));
}

void AtomInfo::LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors)
{
  Vector newPos;
  newPos.LinearCombinationOfVectors(x1,x2,x3,factors);
  setPosition(newPos);
}

/**
 *  returns the kinetic energy of this atom at a given time step
 */
double AtomInfo::getKineticEnergy(const unsigned int _step) const
{
  return getMass() * getAtomicVelocityAtStep(_step).NormSquared();
}

Vector AtomInfo::getMomentum(const unsigned int _step) const
{
  return getMass() * getAtomicVelocityAtStep(_step);
}

/** Decrease the trajectory if given \a MaxSteps is smaller.
 * Does nothing if \a MaxSteps is larger than current size.
 *
 * \param MaxSteps
 */
void AtomInfo::ResizeTrajectory(size_t MaxSteps)
{
  // mind the reverse ordering due to std::greater, latest time steps are at beginning
  VectorTrajectory_t::iterator positer = AtomicPosition.lower_bound(MaxSteps);
  if (positer != AtomicPosition.begin()) {
    if (positer->first == MaxSteps)
      --positer;
    AtomicPosition.erase(AtomicPosition.begin(), positer);
  }
  VectorTrajectory_t::iterator veliter = AtomicVelocity.lower_bound(MaxSteps);
  if (veliter != AtomicVelocity.begin()) {
    if (veliter->first == MaxSteps)
      --veliter;
    AtomicVelocity.erase(AtomicVelocity.begin(), veliter);
  }
  VectorTrajectory_t::iterator forceiter = AtomicForce.lower_bound(MaxSteps);
  if (forceiter != AtomicForce.begin()) {
    if (forceiter->first == MaxSteps)
      --forceiter;
    AtomicForce.erase(AtomicForce.begin(), forceiter);
  }
}

size_t AtomInfo::getTrajectorySize() const
{
  // mind greater comp for map here: first element is latest in time steps!
  return AtomicPosition.begin()->first+1;
}

double AtomInfo::getMass() const
{
  return getType()->getMass();
}

/** Helper function to either insert or assign, depending on the element being
 * present already.
 *
 * \param _trajectory vector of Vectors to assign
 * \param dest step to insert/assign to
 * \param _newvalue new Vector value
 */
void assignTrajectoryElement(
    std::map<unsigned int, Vector, std::greater<unsigned int> > &_trajectory,
    const unsigned int dest,
    const Vector &_newvalue)
{
  std::pair<std::map<unsigned int, Vector, std::greater<unsigned int> >::iterator, bool> inserter =
      _trajectory.insert( std::make_pair(dest, _newvalue) );
  if (!inserter.second)
    inserter.first->second = _newvalue;
}

/** Copies a given trajectory step \a src onto another \a dest
 * \param dest index of destination step
 * \param src index of source step
 */
void AtomInfo::CopyStepOnStep(const unsigned int dest, const unsigned int src)
{
  if (dest == src)  // self assignment check
    return;

  if (WorldTime::getTime() == dest){
    NOTIFY(AtomObservable::PositionChanged);
    NOTIFY(AtomObservable::VelocityChanged);
    NOTIFY(AtomObservable::ForceChanged);
  }

  VectorTrajectory_t::iterator positer = AtomicPosition.find(src);
  ASSERT( positer != AtomicPosition.end(),
      "AtomInfo::CopyStepOnStep() - step "
      +toString(src)+" to copy from not present in AtomicPosition.");
  assignTrajectoryElement(AtomicPosition, dest, positer->second);
  VectorTrajectory_t::iterator veliter = AtomicVelocity.find(src);
  if (veliter != AtomicVelocity.end())
    assignTrajectoryElement(AtomicVelocity, dest, veliter->second);
  VectorTrajectory_t::iterator forceiter = AtomicForce.find(src);
  if (forceiter != AtomicForce.end())
    assignTrajectoryElement(AtomicForce, dest, forceiter->second);
};

/** Performs a velocity verlet update of the position at \a NextStep from \a LastStep information only.
 *
 * We calculate \f$x(t + \delta t) = x(t) + v(t)* \delta t + .5 * \delta t * \delta t * F(t)/m \f$.
 *
 *
 * \param NextStep index of sequential step to set
 * \param Deltat time step width
 * \param IsAngstroem whether the force's underlying unit of length is angstroem or bohr radii
 */
void AtomInfo::VelocityVerletUpdateX(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 position with mass=" << getMass() << " and Deltat="
      << Deltat << " at NextStep=" << NextStep);

  // update position
  {
    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(LastStep) << " from last step " << tempVector);
    setPositionAtStep(NextStep, tempVector);
    LOG(3, "INFO: Position at step " << NextStep << " set to " << tempVector);
  }
};

/** Performs a velocity verlet update of the velocity at \a NextStep.
 *
 * \note forces at NextStep should have been calculated based on position at NextStep prior
 * to calling this function.
 *
 * We calculate \f$v(t) = v(t - \delta t) + \delta _t * .5 * (F(t - \delta t) + F(t))/m \f$.
 *
 * Parameters are according to those in configuration class.
 * \param NextStep index of sequential step to set
 * \param Deltat time step width
 * \param IsAngstroem whether the force's underlying unit of length is angstroem or bohr radii
 */
void AtomInfo::VelocityVerletUpdateU(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 velocity with mass=" << getMass() << " and Deltat="
      << Deltat << " at NextStep=" << NextStep);

  // Update U
  {
    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);
  }
};

std::ostream & AtomInfo::operator << (std::ostream &ost) const
{
  return (ost << getPosition());
}

std::ostream & operator << (std::ostream &ost, const AtomInfo &a)
{
  const size_t terminalstep = a.getTrajectorySize()-1;
  if (terminalstep) {
    ost << "starts at "
        << a.getPositionAtStep(0) << " and ends at "
        << a.getPositionAtStep(terminalstep)
        << " at time step " << terminalstep;
  } else {
    ost << "is at "
        << a.getPositionAtStep(0) << " with a single time step only";
  }
  return ost;
}