Changeset 6625c3 for src

Timestamp:

Feb 24, 2011, 5:46:46 PM (14 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

056e70

Parents:

54b42e

git-author:

Frederik Heber <heber@…> (02/24/11 11:18:59)

git-committer:

Frederik Heber <heber@…> (02/24/11 17:46:46)

Message:

Removed atom_trajectoryparticle*, replaced by AtomInfo class now having std::vector<> for trajectories.

AtomInfo:

• AtomInfo stores position, velocity, and force in vectors.
• getter and setter implemented for velocity and force.
• all functions from TrajectoryParticle put into AtomInfo.
• Direct access to velocity and force exchanged by getter/setter everywhere.
• FIX: molecule::LinearInterpolationBetweenConfiguration() extends trajectory to MaxSteps+1 (for endstep).

Other changes:

• gsl_rng_gaussian() exchanged by boost::random specific type.

Location:

src

Files:

• AtomSet.hpp (modified) (7 diffs)
• Makefile.am (modified) (2 diffs)
• Parser/PcpParser.cpp (modified) (1 diff)
• Thermostats/Berendsen.cpp (modified) (1 diff)
• Thermostats/GaussianThermostat.cpp (modified) (2 diffs)
• Thermostats/Langevin.cpp (modified) (6 diffs)
• Thermostats/Langevin.hpp (modified) (2 diffs)
• Thermostats/NoseHoover.cpp (modified) (2 diffs)
• Thermostats/Woodcock.cpp (modified) (1 diff)
• atom.cpp (modified) (8 diffs)
• atom.hpp (modified) (3 diffs)
• atom_atominfo.cpp (modified) (10 diffs)
• atom_atominfo.hpp (modified) (5 diffs)
• atom_trajectoryparticle.cpp (deleted)
• atom_trajectoryparticle.hpp (deleted)
• atom_trajectoryparticleinfo.cpp (deleted)
• atom_trajectoryparticleinfo.hpp (deleted)
• config.cpp (modified) (3 diffs)
• molecule.cpp (modified) (5 diffs)
• molecule_dynamics.cpp (modified) (12 diffs)
• molecule_geometry.cpp (modified) (6 diffs)

src/AtomSet.hpp

@@ -82 +82 @@
   template<class T>
   struct valueSum {
-    valueSum(T (atom::*_f)() const,T startValue) :
+    valueSum(T (AtomInfo::*_f)() const,T startValue) :
       f(_f),
       value(startValue)
     {}
-    T operator+(atom *atom){
+    T operator+(AtomInfo *atom){
       return value + (atom->*f)();
     }
@@ -93 +93 @@
       return value;
     }
-    T (atom::*f)() const;
+    T (AtomInfo::*f)() const;
     T value;
   };
@@ -99 +99 @@
   template<class T>
   struct stepValueSum {
-    stepValueSum(unsigned int _step, T (atom::*_f)(unsigned int) const,T startValue) :
+    stepValueSum(unsigned int _step, T (AtomInfo::*_f)(unsigned int) const,T startValue) :
       step(_step),
       f(_f),
       value(startValue)
     {}
-    T operator+(atom *atom){
+    T operator+(AtomInfo *atom){
       return value + (atom->*f)(step);
     }
@@ -112 +112 @@
     }
     unsigned int step;
-    T (atom::*f)(unsigned int) const;
+    T (AtomInfo::*f)(unsigned int) const;
     T value;
   };
@@ -119 +119 @@
 template<class Set>
 inline void AtomSetMixin<Set>::translate(const Vector &translater){
-  BOOST_FOREACH(atom *atom,*this){
+  BOOST_FOREACH(AtomInfo *atom,*this){
     *(atom) += translater;
   }
@@ -126 +126 @@
 template<class Set>
 inline void AtomSetMixin<Set>::addVelocityAtStep(const Vector velocity, unsigned int step){
-  BOOST_FOREACH(atom *atom,*this){
-    atom->Trajectory.U.at(step) += velocity;
+  BOOST_FOREACH(AtomInfo *atom,*this){
+    atom->getAtomicVelocity(step) += velocity;
   }
 }
@@ -141 +141 @@
 template<class Set>
 inline double AtomSetMixin<Set>::totalMass() const{
-  return accumulate(this->begin(),this->end(),valueSum<double>(&atom::getMass,0)).value;
+  return accumulate(this->begin(),this->end(),valueSum<double>(&AtomInfo::getMass,0)).value;
 }
 
 template<class Set>
 inline double AtomSetMixin<Set>::totalTemperatureAtStep(unsigned int step) const{
-  return accumulate(this->begin(),this->end(),stepValueSum<double>(step,&atom::getKineticEnergy,0)).value;
+  return accumulate(this->begin(),this->end(),stepValueSum<double>(step,&AtomInfo::getKineticEnergy,0)).value;
 }
 
 template<class Set>
 inline Vector AtomSetMixin<Set>::totalMomentumAtStep(unsigned int step) const{
-  return accumulate(this->begin(),this->end(),stepValueSum<Vector>(step,&atom::getMomentum,Vector())).value;
+  return accumulate(this->begin(),this->end(),stepValueSum<Vector>(step,&AtomInfo::getMomentum,Vector())).value;
 }
 

src/Makefile.am

@@ -21 +21 @@
   atom_graphnode.cpp \
   atom_graphnodeinfo.cpp \
-  atom_particleinfo.cpp \
-  atom_trajectoryparticle.cpp \
-  atom_trajectoryparticleinfo.cpp
+  atom_particleinfo.cpp
 ATOMHEADER = \
   atom.hpp \
@@ -32 +30 @@
   atom_graphnode.hpp \
   atom_graphnodeinfo.hpp \
-  atom_particleinfo.hpp \
-  atom_trajectoryparticle.hpp \
-  atom_trajectoryparticleinfo.hpp
+  atom_particleinfo.hpp
 
 ANALYSISSOURCE = \

src/Parser/PcpParser.cpp

@@ -557 +557 @@
     *file << "Ion_Type" << ZtoIndexMap[Z] << "_" << ZtoCountMap[Z] << "\t"  << fixed << setprecision(9) << showpoint;
     *file << (*AtomRunner)->at(0) << "\t" << (*AtomRunner)->at(1) << "\t" << (*AtomRunner)->at(2);
-    *file << "\t" << (*AtomRunner)->FixedIon;
+    *file << "\t" << (int)((*AtomRunner)->getFixedIon());
     if ((*AtomRunner)->getAtomicVelocity().Norm() > MYEPSILON)
       *file << "\t" << scientific << setprecision(6) << (*AtomRunner)->getAtomicVelocity()[0] << "\t" << (*AtomRunner)->getAtomicVelocity()[1] << "\t" << (*AtomRunner)->getAtomicVelocity()[2] << "\t";

src/Thermostats/Berendsen.cpp

@@ -63 +63 @@
   double ScaleTempFactor = getContainer().TargetTemp/ActualTemp;
   for(ForwardIterator iter=begin;iter!=end;++iter){
-    Vector &U = (*iter)->Trajectory.U.at(step);
-    if ((*iter)->FixedIon == 0) { // even FixedIon moves, only not by other's forces
+    Vector &U = (*iter)->getAtomicVelocity(step);
+    if ((*iter)->getFixedIon() == 0) { // even FixedIon moves, only not by other's forces
       U *= sqrt(1+(World::getInstance().getConfig()->Deltat/TempFrequency)*(ScaleTempFactor-1));
       ekin += 0.5*(*iter)->getType()->getMass() * U.NormSquared();

src/Thermostats/GaussianThermostat.cpp

@@ -70 +70 @@
   double ekin =0;
   for(ForwardIterator iter=begin;iter!=end;++iter){
-    Vector &U = (*iter)->Trajectory.U.at(step);
-    if ((*iter)->FixedIon == 0) {// even FixedIon moves, only not by other's forces
+    Vector &U = (*iter)->getAtomicVelocity(step);
+    if ((*iter)->getFixedIon() == 0) {// even FixedIon moves, only not by other's forces
       U += World::getInstance().getConfig()->Deltat * G_over_E * U;
       ekin += (*iter)->getType()->getMass() * U.NormSquared();
@@ -84 +84 @@
   G=0;
   for(ForwardIterator iter=begin;iter!=end;++iter){
-    Vector &U = (*iter)->Trajectory.U.at(step);
-    Vector &F = (*iter)->Trajectory.F.at(step);
-    if ((*iter)->FixedIon == 0){ // even FixedIon moves, only not by other's forces
+    const Vector &U = (*iter)->getAtomicVelocity(step);
+    const Vector &F = (*iter)->getAtomicForce(step);
+    if ((*iter)->getFixedIon() == 0){ // even FixedIon moves, only not by other's forces
       G += U.ScalarProduct(F);
       E += U.NormSquared()*(*iter)->getType()->getMass();

src/Thermostats/Langevin.cpp

@@ -20 +20 @@
 #include "ThermoStatContainer.hpp"
 
+#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/normal_distribution.hpp>
+#include <boost/random/variate_generator.hpp>
+
 #include "RandomNumbers/RandomNumberGeneratorFactory.hpp"
 #include "RandomNumbers/RandomNumberGenerator.hpp"
@@ -25 +29 @@
 Langevin::Langevin(double _TempFrequency,double _alpha) :
   TempFrequency(_TempFrequency),
-  alpha(_alpha)
-{
-  gsl_rng_env_setup();
-  T = gsl_rng_default;
-  r = gsl_rng_alloc (T);
-}
+  alpha(_alpha),
+  rng_engine(new boost::mt19937),
+  rng_distribution(NULL)
+{}
 
 Langevin::Langevin() :
   TempFrequency(2.5),
-  alpha(0.)
+  alpha(0.),
+  rng_engine(new boost::mt19937),
+  rng_distribution(NULL)
 {}
 
 Langevin::~Langevin()
 {
-  gsl_rng_free (r);
+  delete rng_engine;
+  delete rng_distribution;
 }
 
@@ -76 +81 @@
 double Langevin::doScaleAtoms(unsigned int step,double ActualTemp,ForwardIterator begin, ForwardIterator end){
   DoLog(2) && (Log() << Verbose(2) <<  "Applying Langevin thermostat..." << endl);
-  RandomNumberGenerator &random = RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator();
+  RandomNumberGenerator &random = RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator("mt19937", "normal_distribution");
   const double rng_min = random.min();
   const double rng_max = random.max();
@@ -82 +87 @@
   for(ForwardIterator iter=begin;iter!=end;++iter){
     double sigma  = sqrt(getContainer().TargetTemp/(*iter)->getType()->getMass()); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
-    Vector &U = (*iter)->Trajectory.U.at(step);
-    if ((*iter)->FixedIon == 0) { // even FixedIon moves, only not by other's forces
+    rng_distribution = new boost::normal_distribution<>(0,sigma);
+    boost::variate_generator<boost::mt19937&, boost::normal_distribution<> > rng(*rng_engine, *rng_distribution);
+    Vector &U = (*iter)->getAtomicVelocity(step);
+    if ((*iter)->getFixedIon() == 0) { // even FixedIon moves, only not by other's forces
       // throw a dice to determine whether it gets hit by a heat bath particle
       if (((((random()/(rng_max-rng_min)))*TempFrequency) < 1.)) {
@@ -89 +96 @@
         // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis
         for (int d=0; d<NDIM; d++) {
-          U[d] = gsl_ran_gaussian (r, sigma);
+          U[d] = rng();
         }
         DoLog(2) && (Log() << Verbose(2) << U.Norm() << endl);
@@ -95 +102 @@
       ekin += 0.5*(*iter)->getType()->getMass() * U.NormSquared();
     }
+    delete rng_distribution;
+    rng_distribution = NULL;
   }
   return ekin;

src/Thermostats/Langevin.hpp

@@ -14 +14 @@
 #endif
 
+#include "Thermostats/Thermostat.hpp"
 
-#include "Thermostats/Thermostat.hpp"
+#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/normal_distribution.hpp>
 
 class Langevin : public Thermostat
@@ -38 +40 @@
   double TempFrequency;
   double alpha;
-  gsl_rng * r;
-  const gsl_rng_type * T;
+
+  boost::mt19937 *rng_engine;
+  boost::normal_distribution<> *rng_distribution;
 };
 

src/Thermostats/NoseHoover.cpp

@@ -68 +68 @@
   double ekin =0;
   for(ForwardIterator iter=begin;iter!=end;++iter){
-    Vector &U = (*iter)->Trajectory.U.at(step);
-    if ((*iter)->FixedIon == 0) { // even FixedIon moves, only not by other's forces
+    Vector &U = (*iter)->getAtomicVelocity(step);
+    if ((*iter)->getFixedIon() == 0) { // even FixedIon moves, only not by other's forces
       U += World::getInstance().getConfig()->Deltat/(*iter)->getType()->getMass() * (alpha * (U * (*iter)->getType()->getMass()));
       ekin += (0.5*(*iter)->getType()->getMass()) * U.NormSquared();
@@ -82 +82 @@
   count=0;
   for(ForwardIterator iter = begin;iter!=end;++iter){
-    Vector &U = (*iter)->Trajectory.U.at(step);
-    if ((*iter)->FixedIon == 0) { // even FixedIon moves, only not by other's forces
+    const Vector &U = (*iter)->getAtomicVelocity(step);
+    if ((*iter)->getFixedIon() == 0) { // even FixedIon moves, only not by other's forces
       delta_alpha += U.NormSquared()*(*iter)->getType()->getMass();
     }

src/Thermostats/Woodcock.cpp

@@ -67 +67 @@
     double ekin;
     for (ForwardIterator iter = begin; iter!=end;++iter){
-      Vector &U = (*iter)->Trajectory.U.at(step);
-      if ((*iter)->FixedIon == 0){ // even FixedIon moves, only not by other's forces
+      Vector &U = (*iter)->getAtomicVelocity(step);
+      if ((*iter)->getFixedIon() == 0){ // even FixedIon moves, only not by other's forces
         U *= ScaleTempFactor;
         ekin += 0.5*(*iter)->getType()->getMass() * U.NormSquared();

src/atom.cpp

@@ -55 +55 @@
   setPosition(pointer->getPosition()); // copy coordination
   setAtomicVelocity(pointer->getAtomicVelocity()); // copy velocity
-  FixedIon = pointer->FixedIon;
+  setFixedIon(pointer->getFixedIon());
 };
 
@@ -152 +152 @@
     *out << "Ion_Type" << ElementNo << "_" << AtomNo << "\t"  << fixed << setprecision(9) << showpoint;
     *out << at(0) << "\t" << at(1) << "\t" << at(2);
-    *out << "\t" << FixedIon;
+    *out << "\t" << (int)(getFixedIon());
     if (getAtomicVelocity().Norm() > MYEPSILON)
       *out << "\t" << scientific << setprecision(6) << getAtomicVelocity()[0] << "\t" << getAtomicVelocity()[1] << "\t" << getAtomicVelocity()[2] << "\t";
@@ -179 +179 @@
     *out << "Ion_Type" << elementLookup.there.find(elemental)->second << "_" << AtomNo[elemental->getAtomicNumber()] << "\t"  << fixed << setprecision(9) << showpoint;
     *out << at(0) << "\t" << at(1) << "\t" << at(2);
-    *out << "\t" << FixedIon;
+    *out << "\t" << getFixedIon();
     if (getAtomicVelocity().Norm() > MYEPSILON)
       *out << "\t" << scientific << setprecision(6) << getAtomicVelocity()[0] << "\t" << getAtomicVelocity()[1] << "\t" << getAtomicVelocity()[2] << "\t";
@@ -191 +191 @@
 };
 
-/** Output of a single atom as one lin in xyz file.
+/** Output of a single atom as one line in xyz file.
  * \param *out stream to output to
   * \return true - \a *out present, false - \a *out is NULL
@@ -204 +204 @@
 };
 
-/** Output of a single atom as one lin in xyz file.
+/** Output of a single atom as one line in xyz file.
  * \param *out stream to output to
  * \param *ElementNo array with ion type number in the config file this atom's element shall have
@@ -218 +218 @@
     ASSERT(elementLookup.there.find(elemental)!=elementLookup.there.end(),"Type of this atom was not in the formula upon enumeration");
     *out << "Ion_Type" << elementLookup.there.find(elemental)->second << "_" << AtomNo[getType()->getAtomicNumber()] << "\t"  << fixed << setprecision(9) << showpoint;
-    *out << Trajectory.R.at(step)[0] << "\t" << Trajectory.R.at(step)[1] << "\t" << Trajectory.R.at(step)[2];
-    *out << "\t" << FixedIon;
-    if (Trajectory.U.at(step).Norm() > MYEPSILON)
-      *out << "\t" << scientific << setprecision(6) << Trajectory.U.at(step)[0] << "\t" << Trajectory.U.at(step)[1] << "\t" << Trajectory.U.at(step)[2] << "\t";
-    if (Trajectory.F.at(step).Norm() > MYEPSILON)
-      *out << "\t" << scientific << setprecision(6) << Trajectory.F.at(step)[0] << "\t" << Trajectory.F.at(step)[1] << "\t" << Trajectory.F.at(step)[2] << "\t";
+    *out << getPosition(step)[0] << "\t" << getPosition(step)[1] << "\t" << getPosition(step)[2];
+    *out << "\t" << (int)(getFixedIon());
+    if (getAtomicVelocity(step).Norm() > MYEPSILON)
+      *out << "\t" << scientific << setprecision(6) << getAtomicVelocity(step)[0] << "\t" << getAtomicVelocity(step)[1] << "\t" << getAtomicVelocity(step)[2] << "\t";
+    if (getAtomicForce(step).Norm() > MYEPSILON)
+      *out << "\t" << scientific << setprecision(6) << getAtomicForce(step)[0] << "\t" << getAtomicForce(step)[1] << "\t" << getAtomicForce(step)[2] << "\t";
     *out << "\t# Number in molecule " << nr << endl;
     return true;
@@ -239 +239 @@
   if (out != NULL) {
     *out << getType()->getSymbol() << "\t";
-    *out << Trajectory.R.at(step)[0] << "\t";
-    *out << Trajectory.R.at(step)[1] << "\t";
-    *out << Trajectory.R.at(step)[2] << endl;
+    *out << getPosition(step)[0] << "\t";
+    *out << getPosition(step)[1] << "\t";
+    *out << getPosition(step)[2] << endl;
     return true;
   } else
@@ -383 +383 @@
 }
 
-double atom::getMass() const{
-  return getType()->getMass();
-}
-
 int atom::getNr() const{
   return nr;

src/atom.hpp

@@ -27 +27 @@
 #include "atom_graphnode.hpp"
 #include "atom_particleinfo.hpp"
-#include "atom_trajectoryparticle.hpp"
 #include "TesselPoint.hpp"
 #include "types.hpp"
@@ -45 +44 @@
  * Class incorporates position, type
  */
-class atom : public TrajectoryParticle, public GraphNode, public BondedParticle, public TesselPoint {
+class atom : public GraphNode, public BondedParticle, public TesselPoint {
   friend atom* NewAtom(atomId_t);
   friend void  DeleteAtom(atom*);
@@ -98 +97 @@
    void removeFromMolecule();
 
-   double getMass() const;
    int getNr() const;
 

src/atom_atominfo.cpp

@@ -20 +20 @@
 #include "CodePatterns/MemDebug.hpp"
 
+#include "CodePatterns/Log.hpp"
+#include "CodePatterns/Verbose.hpp"
+#include "config.hpp"
+#include "element.hpp"
+#include "parser.hpp"
 #include "periodentafel.hpp"
 #include "World.hpp"
-#include "element.hpp"
 #include "atom_atominfo.hpp"
 
@@ -28 +32 @@
  */
 AtomInfo::AtomInfo() :
-  AtomicElement(NULL)
+  AtomicElement(NULL),
+  FixedIon(false)
 {
   AtomicPosition.reserve(1);
@@ -42 +47 @@
 AtomInfo::AtomInfo(const AtomInfo &_atom) :
     AtomicPosition(_atom.AtomicPosition),
-    AtomicElement(_atom.AtomicElement)
-{};
+    AtomicElement(_atom.AtomicElement),
+    FixedIon(false)
+{
+  AtomicVelocity.reserve(1);
+  AtomicVelocity.push_back(zeroVec);
+  AtomicForce.reserve(1);
+  AtomicForce.push_back(zeroVec);
+};
 
 AtomInfo::AtomInfo(const VectorInterface &_v) :
-    AtomicElement(NULL)
+    AtomicElement(NULL),
+    FixedIon(false)
 {
   AtomicPosition[0] = _v.getPosition();
+  AtomicVelocity.reserve(1);
+  AtomicVelocity.push_back(zeroVec);
+  AtomicForce.reserve(1);
+  AtomicForce.push_back(zeroVec);
 };
 
@@ -80 +96 @@
 {
   return AtomicPosition[0];
+}
+
+const Vector& AtomInfo::getPosition(const int _step) const
+{
+  ASSERT(_step < AtomicPosition.size(),
+      "AtomInfo::getPosition() - Access out of range.");
+  return AtomicPosition[_step];
 }
 
@@ -96 +119 @@
 }
 
+Vector& AtomInfo::getAtomicVelocity(const int _step)
+{
+  ASSERT(_step < AtomicVelocity.size(),
+      "AtomInfo::getAtomicVelocity() - Access out of range.");
+  return AtomicVelocity[_step];
+}
+
 const Vector& AtomInfo::getAtomicVelocity() const
 {
@@ -101 +131 @@
 }
 
+const Vector& AtomInfo::getAtomicVelocity(const int _step) const
+{
+  ASSERT(_step < AtomicVelocity.size(),
+      "AtomInfo::getAtomicVelocity() - Access out of range.");
+  return AtomicVelocity[_step];
+}
+
 void AtomInfo::setAtomicVelocity(const Vector &_newvelocity)
 {
@@ -106 +143 @@
 }
 
+void AtomInfo::setAtomicVelocity(const int _step, const Vector &_newvelocity)
+{
+  ASSERT(_step <= AtomicVelocity.size(),
+      "AtomInfo::setAtomicVelocity() - Access out of range.");
+  if(_step < (int)AtomicVelocity.size()) {
+    AtomicVelocity[_step] = _newvelocity;
+  } else if (_step == (int)AtomicVelocity.size()) {
+    AtomicVelocity.push_back(_newvelocity);
+  }
+}
+
 const Vector& AtomInfo::getAtomicForce() const
 {
@@ -111 +159 @@
 }
 
+const Vector& AtomInfo::getAtomicForce(const int _step) const
+{
+  ASSERT(_step < AtomicForce.size(),
+      "AtomInfo::getAtomicForce() - Access out of range.");
+  return AtomicForce[_step];
+}
+
 void AtomInfo::setAtomicForce(const Vector &_newforce)
 {
@@ -116 +171 @@
 }
 
+void AtomInfo::setAtomicForce(const int _step, const Vector &_newforce)
+{
+  ASSERT(_step <= AtomicForce.size(),
+      "AtomInfo::setAtomicForce() - Access out of range.");
+  if(_step < (int)AtomicForce.size()) {
+    AtomicForce[_step] = _newforce;
+  } else if (_step == (int)AtomicForce.size()) {
+    AtomicForce.push_back(_newforce);
+  }
+}
+
+bool AtomInfo::getFixedIon() const
+{
+  return FixedIon;
+}
+
+void AtomInfo::setFixedIon(const bool _fixedion)
+{
+  FixedIon = _fixedion;
+}
+
 void AtomInfo::setPosition(const Vector& _vector)
 {
   AtomicPosition[0] = _vector;
+  //cout << "AtomInfo::setPosition: " << getType()->symbol << " at " << getPosition() << endl;
+}
+
+void AtomInfo::setPosition(int _step, const Vector& _vector)
+{
+  ASSERT(_step <= AtomicPosition.size(),
+      "AtomInfo::setPosition() - Access out of range.");
+  if(_step < (int)AtomicPosition.size()) {
+    AtomicPosition[_step] = _vector;
+  } else if (_step == (int)AtomicPosition.size()) {
+    AtomicPosition.push_back(_vector);
+  }
   //cout << "AtomInfo::setPosition: " << getType()->symbol << " at " << getPosition() << endl;
 }
@@ -204 +292 @@
 }
 
+/**
+ *  returns the kinetic energy of this atom at a given time step
+ */
+double AtomInfo::getKineticEnergy(unsigned int step) const{
+  return getMass() * AtomicVelocity[step].NormSquared();
+}
+
+Vector AtomInfo::getMomentum(unsigned int step) const{
+  return getMass()*AtomicVelocity[step];
+}
+
+/** Extends the trajectory STL vector to the new size.
+ * Does nothing if \a MaxSteps is smaller than current size.
+ * \param MaxSteps
+ */
+void AtomInfo::ResizeTrajectory(size_t MaxSteps)
+{
+  if (AtomicPosition.size() <= (unsigned int)(MaxSteps)) {
+    DoLog(0) && (Log() << Verbose(0) << "Increasing size for trajectory array from " << AtomicPosition.size() << " to " << (MaxSteps+1) << "." << endl);
+    AtomicPosition.resize(MaxSteps+1, zeroVec);
+    AtomicVelocity.resize(MaxSteps+1, zeroVec);
+    AtomicForce.resize(MaxSteps+1, zeroVec);
+  }
+};
+
+size_t AtomInfo::getTrajectorySize() const
+{
+  return AtomicPosition.size();
+}
+
+double AtomInfo::getMass() const{
+  return AtomicElement->getMass();
+}
+
+/** 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(int dest, int src)
+{
+  if (dest == src)  // self assignment check
+    return;
+
+  ASSERT(dest < AtomicPosition.size(),
+      "AtomInfo::CopyStepOnStep() - destination outside of current trajectory array.");
+  ASSERT(src < AtomicPosition.size(),
+      "AtomInfo::CopyStepOnStep() - source outside of current trajectory array.");
+  for (int n=NDIM;n--;) {
+    AtomicPosition.at(dest)[n] = AtomicPosition.at(src)[n];
+    AtomicVelocity.at(dest)[n] = AtomicVelocity.at(src)[n];
+    AtomicForce.at(dest)[n] = AtomicForce.at(src)[n];
+  }
+};
+
+/** Performs a velocity verlet update of the trajectory.
+ * Parameters are according to those in configuration class.
+ * \param NextStep index of sequential step to set
+ * \param *configuration pointer to configuration with parameters
+ * \param *Force matrix with forces
+ */
+void AtomInfo::VelocityVerletUpdate(int nr, int NextStep, config *configuration, ForceMatrix *Force, const size_t offset)
+{
+  //a = configuration.Deltat*0.5/walker->type->mass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
+  if ((int)AtomicPosition.size() <= NextStep)
+    ResizeTrajectory(NextStep+10);
+  for (int d=0; d<NDIM; d++) {
+    AtomicForce.at(NextStep)[d] = -Force->Matrix[0][nr][d+offset]*(configuration->GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
+    AtomicPosition.at(NextStep)[d] = AtomicPosition.at(NextStep-1)[d];
+    AtomicPosition.at(NextStep)[d] += configuration->Deltat*(AtomicVelocity.at(NextStep-1)[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
+    AtomicPosition.at(NextStep)[d] += 0.5*configuration->Deltat*configuration->Deltat*(AtomicForce.at(NextStep)[d]/getMass());     // F = m * a and s =
+  }
+  // Update U
+  for (int d=0; d<NDIM; d++) {
+    AtomicVelocity.at(NextStep)[d] = AtomicVelocity.at(NextStep-1)[d];
+    AtomicVelocity.at(NextStep)[d] += configuration->Deltat * (AtomicForce.at(NextStep)[d]+AtomicForce.at(NextStep-1)[d]/getMass()); // v = F/m * t
+  }
+  // Update R (and F)
+//      out << "Integrated position&velocity of step " << (NextStep) << ": (";
+//      for (int d=0;d<NDIM;d++)
+//        out << AtomicPosition.at(NextStep).x[d] << " ";          // next step
+//      out << ")\t(";
+//      for (int d=0;d<NDIM;d++)
+//        Log() << Verbose(0) << AtomicVelocity.at(NextStep).x[d] << " ";          // next step
+//      out << ")" << endl;
+};
+
 const AtomInfo& operator*=(AtomInfo& a, const double m)
 {

src/atom_atominfo.hpp

@@ -27 +27 @@
 
 class AtomInfo;
+class config;
 class element;
+class ForceMatrix;
 class RealSpaceMatrix;
 
@@ -58 +60 @@
   /** Getter for AtomicVelocity.
    *
+   * Current time step is used.
+   *
    * @return constant reference to AtomicVelocity
    */
@@ -63 +67 @@
   /** Getter for AtomicVelocity.
    *
+   * @param _step time step to return
+   * @return constant reference to AtomicVelocity
+   */
+  Vector& getAtomicVelocity(const int _step);
+  /** Getter for AtomicVelocity.
+   *
+   * Current time step is used.
+   *
    * @return constant reference to AtomicVelocity
    */
   const Vector& getAtomicVelocity() const;
+  /** Getter for AtomicVelocity.
+   *
+   * @param _step time step to return
+   * @return constant reference to AtomicVelocity
+   */
+  const Vector& getAtomicVelocity(const int _step) const;
   /** Setter for AtomicVelocity.
    *
+   * Current time step is used.
+   *
    * @param _newvelocity new velocity to set
    */
   void setAtomicVelocity(const Vector &_newvelocity);
+  /** Setter for AtomicVelocity.
+   *
+   * @param _step time step to set
+   * @param _newvelocity new velocity to set
+   */
+  void setAtomicVelocity(const int _step, const Vector &_newvelocity);
 
   /** Getter for AtomicForce.
    *
+   * Current time step is used.
+   *
    * @return constant reference to AtomicForce
    */
   const Vector& getAtomicForce() const;
+  /** Getter for AtomicForce.
+   *
+   * @param _step time step to return
+   * @return constant reference to AtomicForce
+   */
+  const Vector& getAtomicForce(const int _step) const;
   /** Setter for AtomicForce.
    *
+   * Current time step is used.
+   *
    * @param _newvelocity new force vector to set
    */
   void setAtomicForce(const Vector &_newforce);
+  /** Setter for AtomicForce.
+   *
+   * @param _step time step to set
+   * @param _newvelocity new force vector to set
+   */
+  void setAtomicForce(const int _step, const Vector &_newforce);
+
+  /** Getter for FixedIon.
+   *
+   * @return constant reference to FixedIon
+   */
+  bool getFixedIon() const;
+  /** Setter for FixedIon.
+   *
+   * @param _fixedion new state of FixedIon
+   */
+  void setFixedIon(const bool _fixedion);
 
   ///// manipulation of the atomic position
 
   // Accessors ussually come in pairs... and sometimes even more than that
+  /** Getter for AtomicPosition.
+   *
+   * Current time step is used.
+   *
+   * @param i component of vector
+   * @return i-th component of atomic position
+   */
   const double& operator[](size_t i) const;
+  /** Getter for AtomicPosition.
+   *
+   * Current time step is used.
+   *
+   * \sa operator[], this is if instance is a reference.
+   *
+   * @param i component of vector
+   * @return i-th component of atomic position
+   */
   const double& at(size_t i) const;
+  /** Getter for AtomicPosition.
+   *
+   * \sa operator[], this is if instance is a reference.
+   *
+   * @param i index of component of AtomicPosition
+   * @param _step time step to return
+   * @return atomic position at time step _step
+   */
+  const Vector& atStep(size_t i, int _step) const;
+  /** Setter for AtomicPosition.
+   *
+   * Current time step is used.
+   *
+   * @param i component to set
+   * @param value value to set to
+   */
   void set(size_t i, const double value);
+  /** Setter for AtomicPosition.
+   *
+   * @param i component to set
+   * @param _step time step to set
+   * @param value value to set to
+   */
+  void setAtStep(size_t i, int _step, const double value);
+  /** Getter for AtomicPosition.
+   *
+   * Current time step is used.
+   *
+   * @return atomic position
+   */
   const Vector& getPosition() const;
+  /** Getter for AtomicPosition.
+   *
+   * @param _step time step to return
+   * @return atomic position at time step _step
+   */
+  const Vector& getPosition(int _step) const;
 
   // Assignment operator
+  /** Setter for AtomicPosition.
+   *
+   * Current time step is used.
+   *
+   * @param _vector new position to set
+   */
   void setPosition(const Vector& _vector);
+  /** Setter for AtomicPosition.
+   *
+   * @param _step time step to set
+   * @param _vector new position to set for time step _step
+   */
+  void setPosition(int _step, const Vector& _vector);
   class VectorInterface &operator=(const Vector& _vector);
 
@@ -114 +230 @@
   void Scale(const double factor);
 
+  // operations for trajectories
+  void ResizeTrajectory(size_t MaxSteps);
+  size_t getTrajectorySize() const;
+  void CopyStepOnStep(int dest, int src);
+  void VelocityVerletUpdate(int nr, int MDSteps, config *configuration, ForceMatrix *Force, const size_t offset);
+  double getKineticEnergy(unsigned int step) const;
+  Vector getMomentum(unsigned int step) const;
+  double getMass() const;
+
   std::ostream & operator << (std::ostream &ost) const;
 
@@ -122 +247 @@
 
   const element *AtomicElement;  //!< pointing to element
+  bool FixedIon;
 };
 

src/config.cpp

@@ -462 +462 @@
   int No[MAX_ELEMENTS];
   int verbose = 0;
+  Vector tempVector;
+  int _fixedion;
   double value[3];
 
@@ -525 +527 @@
                     ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &position[1], 1, (repetition == 0) ? critical : optional) &&
                     ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &position[2], 1, (repetition == 0) ? critical : optional) &&
-                    ParseForParameter(verbose,FileBuffer, keyword, 0, 4, 1, int_type, &neues->FixedIon, 1, (repetition == 0) ? critical : optional));
+                    ParseForParameter(verbose,FileBuffer, keyword, 0, 4, 1, int_type, &_fixedion, 1, (repetition == 0) ? critical : optional));
             if (!status)
               break;
-            neues ->setPosition(position);
+            neues->setFixedIon(_fixedion);
+            neues->setPosition(position);
 
             // check size of vectors
-            if (neues->Trajectory.R.size() <= (unsigned int)(repetition)) {
+            if (neues->getTrajectorySize() <= (unsigned int)(repetition)) {
               //Log() << Verbose(0) << "Increasing size for trajectory array of " << keyword << " to " << (repetition+10) << "." << endl;
-              neues->Trajectory.R.resize(repetition+10);
-              neues->Trajectory.U.resize(repetition+10);
-              neues->Trajectory.F.resize(repetition+10);
+              neues->ResizeTrajectory(repetition+1);
+              neues->ResizeTrajectory(repetition+1);
+              neues->ResizeTrajectory(repetition+1);
             }
 
-            // put into trajectories list
-            for (int d=0;d<NDIM;d++)
-              neues->Trajectory.R.at(repetition)[d] = neues->at(d);
-
             // parse velocities if present
-            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &neues->getAtomicVelocity()[0], 1,optional))
-              neues->getAtomicVelocity()[0] = 0.;
-            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &neues->getAtomicVelocity()[1], 1,optional))
-              neues->getAtomicVelocity()[1] = 0.;
-            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &neues->getAtomicVelocity()[2], 1,optional))
-              neues->getAtomicVelocity()[2] = 0.;
-            for (int d=0;d<NDIM;d++)
-              neues->Trajectory.U.at(repetition)[d] = neues->getAtomicVelocity()[d];
+            tempVector.Zero();
+            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &tempVector[0], 1,optional))
+              tempVector[0] = 0.;
+            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &tempVector[1], 1,optional))
+              tempVector[1] = 0.;
+            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &tempVector[2], 1,optional))
+              tempVector[2] = 0.;
+            neues->setAtomicVelocity(repetition, tempVector);
 
             // parse forces if present
-            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 8, 1, double_type, &value[0], 1,optional))
-              value[0] = 0.;
-            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 9, 1, double_type, &value[1], 1,optional))
-              value[1] = 0.;
-            if(!ParseForParameter(verbose,FileBuffer, keyword, 1, 10, 1, double_type, &value[2], 1,optional))
-              value[2] = 0.;
-            for (int d=0;d<NDIM;d++)
-              neues->Trajectory.F.at(repetition)[d] = value[d];
+            tempVector.Zero();
+            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 8, 1, double_type, &tempVector[0], 1,optional))
+              tempVector[0] = 0.;
+            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 9, 1, double_type, &tempVector[1], 1,optional))
+              tempVector[1] = 0.;
+            if(!ParseForParameter(verbose,FileBuffer, keyword, 1, 10, 1, double_type, &tempVector[2], 1,optional))
+              tempVector[2] = 0.;
+            neues->setAtomicForce(repetition, tempVector);
 
   //            Log() << Verbose(0) << "Parsed position of step " << (repetition) << ": (";
-  //            for (int d=0;d<NDIM;d++)
-  //              Log() << Verbose(0) << neues->Trajectory.R.at(repetition).x[d] << " ";          // next step
-  //            Log() << Verbose(0) << ")\t(";
-  //            for (int d=0;d<NDIM;d++)
-  //              Log() << Verbose(0) << neues->Trajectory.U.at(repetition).x[d] << " ";          // next step
-  //            Log() << Verbose(0) << ")\t(";
-  //            for (int d=0;d<NDIM;d++)
-  //              Log() << Verbose(0) << neues->Trajectory.F.at(repetition).x[d] << " ";          // next step
+  //            Log() << Verbose(0) << neues->getPosition(repetition) << " ";          // next step
+  //            Log() << Verbose(0) << neues->getAtomicVelocity(repetition) << " ";          // next step
+  //            Log() << Verbose(0) << neues->getAtomicForce(repetition) << " ";          // next step
   //            Log() << Verbose(0) << ")" << endl;
           }
@@ -603 +597 @@
             neues = AtomList[i][j];
           // then parse for each atom the coordinates as often as present
-          ParseForParameter(verbose,FileBuffer, keyword, 0, 1, 1, double_type, &position[0], repetition,critical);
-          ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &position[1], repetition,critical);
-          ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &position[2], repetition,critical);
-          neues->setPosition(position);
-          ParseForParameter(verbose,FileBuffer, keyword, 0, 4, 1, int_type, &neues->FixedIon, repetition,critical);
-          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &neues->getAtomicVelocity()[0], repetition,optional))
-            neues->getAtomicVelocity()[0] = 0.;
-          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &neues->getAtomicVelocity()[1], repetition,optional))
-            neues->getAtomicVelocity()[1] = 0.;
-          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &neues->getAtomicVelocity()[2], repetition,optional))
-            neues->getAtomicVelocity()[2] = 0.;
+          tempVector.Zero();
+          ParseForParameter(verbose,FileBuffer, keyword, 0, 1, 1, double_type, &tempVector[0], repetition,critical);
+          ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &tempVector[1], repetition,critical);
+          ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &tempVector[2], repetition,critical);
+          neues->setPosition(tempVector);
+          // fixed ion
+          ParseForParameter(verbose,FileBuffer, keyword, 0, 4, 1, int_type, &_fixedion, repetition,critical);
+          neues->setFixedIon(_fixedion);
+          // and velocity
+          tempVector.Zero();
+          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &tempVector[0], repetition,optional))
+            tempVector[0] = 0.;
+          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &tempVector[1], repetition,optional))
+            tempVector[1] = 0.;
+          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &tempVector[2], repetition,optional))
+            tempVector[2] = 0.;
+          neues->setAtomicVelocity(tempVector);
           // here we don't care if forces are present (last in trajectories is always equal to current position)
           neues->setType(elementhash[i]); // find element type

src/molecule.cpp

@@ -365 +365 @@
       FirstOtherAtom->setType(1);  // element is Hydrogen
       FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
-      FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
+      FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
       if (TopReplacement->getType()->getAtomicNumber() == 1) { // neither rescale nor replace if it's already hydrogen
         FirstOtherAtom->father = TopReplacement;
@@ -430 +430 @@
       SecondOtherAtom->setType(1);
       FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
-      FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
+      FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
       SecondOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
-      SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
+      SecondOtherAtom->setFixedIon(TopReplacement->getFixedIon());
       FirstOtherAtom->father = NULL;  // we are just an added hydrogen with no father
       SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
@@ -485 +485 @@
       ThirdOtherAtom->setType(1);
       FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
-      FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
+      FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
       SecondOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
-      SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
+      SecondOtherAtom->setFixedIon(TopReplacement->getFixedIon());
       ThirdOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
-      ThirdOtherAtom->FixedIon = TopReplacement->FixedIon;
+      ThirdOtherAtom->setFixedIon(TopReplacement->getFixedIon());
       FirstOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
       SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
@@ -582 +582 @@
   istringstream *input = NULL;
   int NumberOfAtoms = 0; // atom number in xyz read
-  int i, j; // loop variables
+  int i; // loop variables
   atom *Walker = NULL;  // pointer to added atom
   char shorthand[3];  // shorthand for atom name
@@ -618 +618 @@
       Walker->setType(1);
     }
-    if (Walker->Trajectory.R.size() <= (unsigned int)MDSteps) {
-      Walker->Trajectory.R.resize(MDSteps+10);
-      Walker->Trajectory.U.resize(MDSteps+10);
-      Walker->Trajectory.F.resize(MDSteps+10);
+    if (Walker->getTrajectorySize() <= (unsigned int)MDSteps) {
+      Walker->ResizeTrajectory(MDSteps+10);
     }
     Walker->setPosition(Vector(x));
-    for(j=NDIM;j--;) {
-      Walker->Trajectory.R.at(MDSteps-1)[j] = x[j];
-      Walker->Trajectory.U.at(MDSteps-1)[j] = 0;
-      Walker->Trajectory.F.at(MDSteps-1)[j] = 0;
-    }
+    Walker->setPosition(MDSteps-1, Vector(x));
+    Walker->setAtomicVelocity(MDSteps-1, zeroVec);
+    Walker->setAtomicForce(MDSteps-1, zeroVec);
     AddAtom(Walker);  // add to molecule
     delete(item);

src/molecule_dynamics.cpp

@@ -60 +60 @@
     // determine normalized trajectories direction vector (n1, n2)
     Sprinter = Params.PermutationMap[Walker->nr];   // find first target point
-    trajectory1 = Sprinter->Trajectory.R.at(Params.endstep) - Walker->Trajectory.R.at(Params.startstep);
+    trajectory1 = Sprinter->getPosition(Params.endstep) - Walker->getPosition(Params.startstep);
     trajectory1.Normalize();
     Norm1 = trajectory1.Norm();
     Sprinter = Params.PermutationMap[(*iter)->nr];   // find second target point
-    trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - (*iter)->Trajectory.R.at(Params.startstep);
+    trajectory2 = Sprinter->getPosition(Params.endstep) - (*iter)->getPosition(Params.startstep);
     trajectory2.Normalize();
     Norm2 = trajectory1.Norm();
     // check whether either is zero()
     if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
-      tmp = Walker->Trajectory.R.at(Params.startstep).distance((*iter)->Trajectory.R.at(Params.startstep));
+      tmp = Walker->getPosition(Params.startstep).distance((*iter)->getPosition(Params.startstep));
     } else if (Norm1 < MYEPSILON) {
       Sprinter = Params.PermutationMap[Walker->nr];   // find first target point
-      trajectory1 = Sprinter->Trajectory.R.at(Params.endstep) - (*iter)->Trajectory.R.at(Params.startstep);
+      trajectory1 = Sprinter->getPosition(Params.endstep) - (*iter)->getPosition(Params.startstep);
       trajectory2 *= trajectory1.ScalarProduct(trajectory2); // trajectory2 is scaled to unity, hence we don't need to divide by anything
       trajectory1 -= trajectory2;   // project the part in norm direction away
@@ -78 +78 @@
     } else if (Norm2 < MYEPSILON) {
       Sprinter = Params.PermutationMap[(*iter)->nr];   // find second target point
-      trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - Walker->Trajectory.R.at(Params.startstep);  // copy second offset
+      trajectory2 = Sprinter->getPosition(Params.endstep) - Walker->getPosition(Params.startstep);  // copy second offset
       trajectory1 *= trajectory2.ScalarProduct(trajectory1); // trajectory1 is scaled to unity, hence we don't need to divide by anything
       trajectory2 -= trajectory1;   // project the part in norm direction away
@@ -87 +87 @@
   //        Log() << Verbose(0) << " and ";
   //        Log() << Verbose(0) << trajectory2;
-      tmp = Walker->Trajectory.R.at(Params.startstep).distance((*iter)->Trajectory.R.at(Params.startstep));
+      tmp = Walker->getPosition(Params.startstep).distance((*iter)->getPosition(Params.startstep));
   //        Log() << Verbose(0) << " with distance " << tmp << "." << endl;
     } else { // determine distance by finding minimum distance
@@ -106 +106 @@
         gsl_matrix_set(A, 1, i, trajectory2[i]);
         gsl_matrix_set(A, 2, i, normal[i]);
-        gsl_vector_set(x,i, (Walker->Trajectory.R.at(Params.startstep)[i] - (*iter)->Trajectory.R.at(Params.startstep)[i]));
+        gsl_vector_set(x,i, (Walker->getPosition(Params.startstep)[i] - (*iter)->getPosition(Params.startstep)[i]));
       }
       // solve the linear system by Householder transformations
@@ -117 +117 @@
       trajectory2.Scale(gsl_vector_get(x,1));
       normal.Scale(gsl_vector_get(x,2));
-      TestVector = (*iter)->Trajectory.R.at(Params.startstep) + trajectory2 + normal
-                   - (Walker->Trajectory.R.at(Params.startstep) + trajectory1);
+      TestVector = (*iter)->getPosition(Params.startstep) + trajectory2 + normal
+                   - (Walker->getPosition(Params.startstep) + trajectory1);
       if (TestVector.Norm() < MYEPSILON) {
   //          Log() << Verbose(2) << "Test: ok.\tDistance of " << tmp << " is correct." << endl;
@@ -150 +150 @@
   //    atom *Sprinter = PermutationMap[Walker->nr];
   //        Log() << Verbose(0) << *Walker << " and " << *(*iter) << " are heading to the same target at ";
-  //        Log() << Verbose(0) << Sprinter->Trajectory.R.at(endstep);
+  //        Log() << Verbose(0) << Sprinter->getPosition(endstep);
   //        Log() << Verbose(0) << ", penalting." << endl;
       result += Params.PenaltyConstants[2];
@@ -183 +183 @@
     // first term: distance to target
     Runner = Params.PermutationMap[(*iter)->nr];   // find target point
-    tmp = ((*iter)->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.endstep)));
+    tmp = ((*iter)->getPosition(Params.startstep).distance(Runner->getPosition(Params.endstep)));
     tmp *= Params.IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
     result += Params.PenaltyConstants[0] * tmp;
@@ -239 +239 @@
   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     for (molecule::const_iterator runner = mol->begin(); runner != mol->end(); ++runner) {
-      Params.DistanceList[(*iter)->nr]->insert( DistancePair((*iter)->Trajectory.R.at(Params.startstep).distance((*runner)->Trajectory.R.at(Params.endstep)), (*runner)) );
+      Params.DistanceList[(*iter)->nr]->insert( DistancePair((*iter)->getPosition(Params.startstep).distance((*runner)->getPosition(Params.endstep)), (*runner)) );
     }
   }
@@ -489 +489 @@
     // set forces
     for (int i=NDIM;i++;)
-      Force->Matrix[0][_atom->nr][5+i] += 2.*constant*sqrt(_atom->Trajectory.R.at(startstep).distance(Sprinter->Trajectory.R.at(endstep)));
+      Force->Matrix[0][_atom->nr][5+i] += 2.*constant*sqrt(_atom->getPosition(startstep).distance(Sprinter->getPosition(endstep)));
   }
   DoLog(1) && (Log() << Verbose(1) << "done." << endl);
@@ -525 +525 @@
 
   // check whether we have sufficient space in Trajectories for each atom
-  for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::ResizeTrajectory),MaxSteps));
+  LOG(1, "STATUS: Extending each trajectory size to " << MaxSteps+1 << ".");
+  for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::ResizeTrajectory),MaxSteps+1));
   // push endstep to last one
   for_each(atoms.begin(),atoms.end(),boost::bind(&atom::CopyStepOnStep,_1,MaxSteps,endstep));
@@ -538 +539 @@
       // add to molecule list
       Sprinter = mol->AddCopyAtom((*iter));
-      for (int n=NDIM;n--;) {
-        Sprinter->set(n, (*iter)->Trajectory.R.at(startstep)[n] + (PermutationMap[(*iter)->nr]->Trajectory.R.at(endstep)[n] - (*iter)->Trajectory.R.at(startstep)[n])*((double)step/(double)MaxSteps));
-        // add to Trajectories
-        //Log() << Verbose(3) << step << ">=" << MDSteps-1 << endl;
-        if (step < MaxSteps) {
-          (*iter)->Trajectory.R.at(step)[n] = (*iter)->Trajectory.R.at(startstep)[n] + (PermutationMap[(*iter)->nr]->Trajectory.R.at(endstep)[n] - (*iter)->Trajectory.R.at(startstep)[n])*((double)step/(double)MaxSteps);
-          (*iter)->Trajectory.U.at(step)[n] = 0.;
-          (*iter)->Trajectory.F.at(step)[n] = 0.;
-        }
-      }
+      // add to Trajectories
+      Vector temp = (*iter)->getPosition(startstep) + (PermutationMap[(*iter)->nr]->getPosition(endstep) - (*iter)->getPosition(startstep))*((double)step/(double)MaxSteps);
+      Sprinter->setPosition(temp);
+      (*iter)->setAtomicVelocity(step, zeroVec);
+      (*iter)->setAtomicForce(step, zeroVec);
+      //Log() << Verbose(3) << step << ">=" << MDSteps-1 << endl;
     }
   }
@@ -620 +617 @@
   // and perform Verlet integration for each atom with position, velocity and force vector
   // check size of vectors
-  for_each(atoms.begin(),
-           atoms.end(),
-           boost::bind(&atom::VelocityVerletUpdate,_1,MDSteps+1, &configuration, &Force, (const size_t) 0));
+  BOOST_FOREACH(atom *_atom, atoms) {
+    _atom->VelocityVerletUpdate(_atom->nr, MDSteps+1, &configuration, &Force, (const size_t) 0);
+  }
 
   // correct velocities (rather momenta) so that center of mass remains motionless

src/molecule_geometry.cpp

@@ -58 +58 @@
   // go through all atoms
   BOOST_FOREACH(atom* iter, atoms){
+    std::cout << "atom before is at " << *iter << std::endl;
     *iter -= *Center;
-    *iter -= *CenterBox;
+    *iter += *CenterBox;
+    std::cout << "atom after is at " << *iter << std::endl;
   }
   atoms.transformNodes(boost::bind(&Box::WrapPeriodically,domain,_1));
@@ -314 +316 @@
 {
   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    for (int j=0;j<MDSteps;j++)
-      (*iter)->Trajectory.R.at(j).ScaleAll(*factor);
-    (*iter)->ScaleAll(*factor);
+    for (size_t j=0;j<(*iter)->getTrajectorySize();j++) {
+      Vector temp = (*iter)->getPosition(j);
+      temp.ScaleAll(*factor);
+      (*iter)->setPosition(j,temp);
+    }
   }
 };
@@ -326 +330 @@
 {
   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    for (int j=0;j<MDSteps;j++)
-      (*iter)->Trajectory.R.at(j) += (*trans);
-    *(*iter) += (*trans);
+    for (size_t j=0;j<(*iter)->getTrajectorySize();j++) {
+      Vector temp = (*iter)->getPosition(j);
+      temp += (*trans);
+      (*iter)->setPosition(j, temp);
+    }
   }
 };
@@ -334 +340 @@
 /** Translate the molecule periodically in the box.
  * \param trans[] translation vector.
- * TODO treatment of trajetories missing
+ * TODO treatment of trajectories missing
  */
 void molecule::TranslatePeriodically(const Vector *trans)
@@ -438 +444 @@
     (*iter)->set(2, -sin(alpha) * tmp + cos(alpha) * (*iter)->at(2));
     for (int j=0;j<MDSteps;j++) {
-      tmp = (*iter)->Trajectory.R.at(j)[0];
-      (*iter)->Trajectory.R.at(j)[0] =  cos(alpha) * tmp + sin(alpha) * (*iter)->Trajectory.R.at(j)[2];
-      (*iter)->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->Trajectory.R.at(j)[2];
+      Vector temp;
+      temp[0] =  cos(alpha) * (*iter)->getPosition(j)[0] + sin(alpha) * (*iter)->getPosition(j)[2];
+      temp[2] = -sin(alpha) * (*iter)->getPosition(j)[0] + cos(alpha) * (*iter)->getPosition(j)[2];
+      (*iter)->setPosition(j,temp);
     }
   }
@@ -457 +464 @@
     (*iter)->set(2, -sin(alpha) * tmp + cos(alpha) * (*iter)->at(2));
     for (int j=0;j<MDSteps;j++) {
-      tmp = (*iter)->Trajectory.R.at(j)[1];
-      (*iter)->Trajectory.R.at(j)[1] =  cos(alpha) * tmp + sin(alpha) * (*iter)->Trajectory.R.at(j)[2];
-      (*iter)->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->Trajectory.R.at(j)[2];
+      Vector temp;
+      temp[1] =  cos(alpha) * (*iter)->getPosition(j)[1] + sin(alpha) * (*iter)->getPosition(j)[2];
+      temp[2] = -sin(alpha) * (*iter)->getPosition(j)[1] + cos(alpha) * (*iter)->getPosition(j)[2];
+      (*iter)->setPosition(j,temp);
     }
   }