Index: molecuilder/src/Makefile.am
===================================================================
--- molecuilder/src/Makefile.am	(revision 3b0ba203af31dbee5c5a38092340126009136ed1)
+++ molecuilder/src/Makefile.am	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -1,4 +1,7 @@
-SOURCE = analysis_correlation.cpp atom.cpp bond.cpp boundary.cpp config.cpp element.cpp ellipsoid.cpp graph.cpp helpers.cpp leastsquaremin.cpp linkedcell.cpp memoryusageobserver.cpp moleculelist.cpp molecule.cpp molecule_dynamics.cpp molecule_fragmentation.cpp molecule_geometry.cpp molecule_graph.cpp molecule_pointcloud.cpp parser.cpp periodentafel.cpp tesselation.cpp tesselationhelpers.cpp vector.cpp verbose.cpp
-HEADER = analysis_correlation.hpp atom.hpp bond.hpp boundary.hpp config.hpp defs.hpp element.hpp ellipsoid.hpp graph.hpp helpers.hpp leastsquaremin.hpp linkedcell.hpp lists.hpp memoryallocator.hpp memoryusageobserver.hpp molecule.hpp molecule_template.hpp parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp vector.hpp verbose.hpp
+ATOMSOURCE = atom.cpp atom_atominfo.cpp atom_bondedparticle.cpp atom_bondedparticleinfo.cpp atom_graphnode.cpp atom_graphnodeinfo.cpp atom_particleinfo.cpp atom_trajectoryparticle.cpp atom_trajectoryparticleinfo.cpp
+ATOMHEADER = atom.hpp atom_atominfo.hpp atom_bondedparticle.hpp atom_bondedparticleinfo.hpp atom_graphnode.hpp atom_graphnodeinfo.hpp atom_particleinfo.hpp atom_trajectoryparticle.hpp atom_trajectoryparticleinfo.hpp
+
+SOURCE = analysis_correlation.cpp ${ATOMSOURCE} bond.cpp boundary.cpp config.cpp element.cpp ellipsoid.cpp graph.cpp helpers.cpp leastsquaremin.cpp linkedcell.cpp memoryusageobserver.cpp moleculelist.cpp molecule.cpp molecule_dynamics.cpp molecule_fragmentation.cpp molecule_geometry.cpp molecule_graph.cpp molecule_pointcloud.cpp parser.cpp periodentafel.cpp tesselation.cpp tesselationhelpers.cpp vector.cpp verbose.cpp
+HEADER = analysis_correlation.hpp ${ATOMHEADER} bond.hpp boundary.hpp config.hpp defs.hpp element.hpp ellipsoid.hpp graph.hpp helpers.hpp leastsquaremin.hpp linkedcell.hpp lists.hpp memoryallocator.hpp memoryusageobserver.hpp molecule.hpp molecule_template.hpp parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp vector.hpp verbose.hpp
 
 BOOST_LIB = $(BOOST_LDFLAGS) $(BOOST_MPL_LIB)
Index: molecuilder/src/atom.cpp
===================================================================
--- molecuilder/src/atom.cpp	(revision 3b0ba203af31dbee5c5a38092340126009136ed1)
+++ molecuilder/src/atom.cpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -15,87 +15,4 @@
 
 /************************************* Functions for class atom *************************************/
-
-/** Constructor of class AtomInfo.
- */
-AtomInfo::AtomInfo() : type(NULL) {};
-
-/** Destructor of class AtomInfo.
- */
-AtomInfo::~AtomInfo()
-{
-};
-
-/** Constructor of class TrajectoryParticleInfo.
- */
-TrajectoryParticleInfo::TrajectoryParticleInfo() : FixedIon(0) {};
-
-/** Destructor of class TrajectoryParticleInfo.
- */
-TrajectoryParticleInfo::~TrajectoryParticleInfo()
-{
-  Trajectory.R.clear();
-  Trajectory.U.clear();
-  Trajectory.F.clear();
-};
-
-/** Constructor of class TrajectoryParticle.
- */
-TrajectoryParticle::TrajectoryParticle()
-{
-};
-
-/** Destructor of class TrajectoryParticle.
- */
-TrajectoryParticle::~TrajectoryParticle()
-{
-};
-
-/** Constructor of class GraphNodeInfo.
- */
-GraphNodeInfo::GraphNodeInfo() : GraphNr(-1), ComponentNr(NULL), LowpointNr(-1), SeparationVertex(false), IsCyclic(false), Ancestor(NULL) {};
-
-/** Destructor of class GraphNodeInfo.
- */
-GraphNodeInfo::~GraphNodeInfo()
-{
-  Free<int>(&ComponentNr, "atom::~atom: *ComponentNr");
-};
-
-/** Constructor of class GraphNode.
- */
-GraphNode::GraphNode()
-{
-};
-
-/** Destructor of class GraphNode.
- */
-GraphNode::~GraphNode()
-{
-};
-
-/** Constructor of class BondedParticleInfo.
- */
-BondedParticleInfo::BondedParticleInfo() : AdaptiveOrder(0), MaxOrder(false) {};
-
-/** Destructor of class BondedParticleInfo.
- */
-BondedParticleInfo::~BondedParticleInfo()
-{
-};
-
-/** Constructor of class BondedParticle.
- */
-BondedParticle::BondedParticle(){};
-
-/** Destructor of class BondedParticle.
- */
-BondedParticle::~BondedParticle()
-{
-  BondList::const_iterator Runner;
-  while (!ListOfBonds.empty()) {
-    Runner = ListOfBonds.begin();
-    removewithoutcheck(*Runner);
-  }
-};
 
 /** Constructor of class atom.
@@ -360,374 +277,2 @@
 };
 
-/** Output graph info of this atom.
- * \param *out output stream
- */
-void GraphNode::OutputGraphInfo(ofstream *out) const
-{
-  *out << Verbose(2) << "Atom " << Name << " is " << ((SeparationVertex) ? "a" : "not a") << " separation vertex, components are ";
-  OutputComponentNumber(out);
-  *out << " with Lowpoint " << LowpointNr << " and Graph Nr. " << GraphNr << "." << endl;
-};
-
-/** Output a list of flags, stating whether the bond was visited or not.
- * Note, we make use of the last entry of the ComponentNr always being -1 if allocated.
- * \param *out output stream for debugging
- */
-void GraphNode::OutputComponentNumber(ofstream *out) const
-{
-  if (ComponentNr != NULL) {
-    for (int i=0; ComponentNr[i] != -1; i++)
-      *out << ComponentNr[i] << " ";
-  }
-};
-
-/** Outputs the current atom::AdaptiveOrder and atom::MaxOrder to \a *file.
- * \param *file output stream
- */
-void BondedParticle::OutputOrder(ofstream *file)
-{
-  *file << nr << "\t" << (int)AdaptiveOrder << "\t" << (int)MaxOrder << endl;
-  //cout << Verbose(2) << "Storing: " << nr << "\t" << (int)AdaptiveOrder << "\t" << (int)MaxOrder << "." << endl;
-};
-
-/** Prints all bonds of this atom with total degree.
- * \param *out stream to output to
- * \return true - \a *out present, false - \a *out is NULL
- */
-bool BondedParticle::OutputBondOfAtom(ofstream *out) const
-{
-  if (out != NULL) {
-    *out << Verbose(4) << "Atom " << Name << "/" << nr << " with " << ListOfBonds.size() << " bonds: ";
-    int TotalDegree = 0;
-    for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); ++Runner) {
-      *out << **Runner << "\t";
-      TotalDegree += (*Runner)->BondDegree;
-    }
-    *out << " -- TotalDegree: " << TotalDegree << endl;
-    return true;
-  } else
-    return false;
-};
-
-/** Output of atom::nr along with all bond partners.
- * \param *AdjacencyFile output stream
- */
-void BondedParticle::OutputAdjacency(ofstream *AdjacencyFile) const
-{
-  *AdjacencyFile << nr << "\t";
-  for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner))
-    *AdjacencyFile << (*Runner)->GetOtherAtom(this)->nr << "\t";
-  *AdjacencyFile << endl;
-};
-
-/** Puts a given bond into atom::ListOfBonds.
- * \param *Binder bond to insert
- */
-bool BondedParticle::RegisterBond(bond *Binder)
-{
-  bool status = false;
-  if (Binder != NULL) {
-    if (Binder->Contains(this)) {
-      ListOfBonds.push_back(Binder);
-      status = true;
-    } else {
-      cout << Verbose(1) << "ERROR: " << *Binder << " does not contain " << *this << "." << endl;
-    }
-  } else {
-    cout << Verbose(1) << "ERROR: Binder is " << Binder << "." << endl;
-  }
-  return status;
-};
-
-/** Removes a given bond from atom::ListOfBonds.
- * \param *Binder bond to remove
- */
-bool BondedParticle::UnregisterBond(bond *Binder)
-{
-  bool status = false;
-  if (Binder != NULL) {
-    if (Binder->Contains(this)) {
-      ListOfBonds.remove(Binder);
-      status = true;
-    } else {
-      cout << Verbose(1) << "ERROR: " << *Binder << " does not contain " << *this << "." << endl;
-    }
-  } else {
-    cout << Verbose(1) << "ERROR: Binder is " << Binder << "." << endl;
-  }
-  return status;
-};
-
-/** Removes all bonds from atom::ListOfBonds.
- * \note Does not do any memory de-allocation.
- */
-void BondedParticle::UnregisterAllBond()
-{
-  ListOfBonds.clear();
-};
-
-/** Corrects the bond degree by one at most if necessary.
- * \param *out output stream for debugging
- */
-int BondedParticle::CorrectBondDegree(ofstream *out)
-{
-  int NoBonds = 0;
-  int OtherNoBonds = 0;
-  int FalseBondDegree = 0;
-  atom *OtherWalker = NULL;
-  bond *CandidateBond = NULL;
-
-  *out << Verbose(3) << "Walker " << *this << ": " << (int)this->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
-  NoBonds = CountBonds();
-  if ((int)(type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check all bonding partners for mismatch
-    for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner)) {
-      OtherWalker = (*Runner)->GetOtherAtom(this);
-      OtherNoBonds = OtherWalker->CountBonds();
-      *out << Verbose(3) << "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
-      if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) { // check if possible candidate
-        if ((CandidateBond == NULL) || (ListOfBonds.size() > OtherWalker->ListOfBonds.size())) { // pick the one with fewer number of bonds first
-          CandidateBond = (*Runner);
-          *out << Verbose(3) << "New candidate is " << *CandidateBond << "." << endl;
-        }
-      }
-    }
-    if ((CandidateBond != NULL)) {
-      CandidateBond->BondDegree++;
-      *out << Verbose(2) << "Increased bond degree for bond " << *CandidateBond << "." << endl;
-    } else {
-      *out << Verbose(2) << "Could not find correct degree for atom " << *this << "." << endl;
-      FalseBondDegree++;
-    }
-  }
-  return FalseBondDegree;
-};
-
-/** Adds kinetic energy of this atom to given temperature value.
- * \param *temperature add on this value
- * \param step given step of trajectory to add
- */
-void TrajectoryParticle::AddKineticToTemperature(double *temperature, int step) const
-{
-  for (int i=NDIM;i--;)
-    *temperature += type->mass * Trajectory.U.at(step).x[i]* Trajectory.U.at(step).x[i];
-};
-
-/** Evaluates some constraint potential if atom moves from \a startstep at once to \endstep in trajectory.
- * \param startstep trajectory begins at
- * \param endstep trajectory ends at
- * \param **PermutationMap if atom switches places with some other atom, there is no translation but a permutaton noted here (not in the trajectories of each).
- * \param *Force Force matrix to store result in
- */
-void TrajectoryParticle::EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force)
-{
-  double constant = 10.;
-  TrajectoryParticle *Sprinter = PermutationMap[nr];
-  // set forces
-  for (int i=NDIM;i++;)
-    Force->Matrix[0][nr][5+i] += 2.*constant*sqrt(Trajectory.R.at(startstep).Distance(&Sprinter->Trajectory.R.at(endstep)));
-};
-
-/** Correct velocity against the summed \a CoGVelocity for \a step.
- * \param *ActualTemp sum up actual temperature meanwhile
- * \param Step MD step in atom::Tracjetory
- * \param *CoGVelocity remnant velocity (i.e. vector sum of all atom velocities)
- */
-void TrajectoryParticle::CorrectVelocity(double *ActualTemp, int Step, Vector *CoGVelocity)
-{
-  for(int d=0;d<NDIM;d++) {
-    Trajectory.U.at(Step).x[d] -= CoGVelocity->x[d];
-    *ActualTemp += 0.5 * type->mass * Trajectory.U.at(Step).x[d] * Trajectory.U.at(Step).x[d];
-  }
-};
-
-/** Extends the trajectory STL vector to the new size.
- * Does nothing if \a MaxSteps is smaller than current size.
- * \param MaxSteps
- */
-void TrajectoryParticle::ResizeTrajectory(int MaxSteps)
-{
-  if (Trajectory.R.size() <= (unsigned int)(MaxSteps)) {
-    //cout << "Increasing size for trajectory array of " << keyword << " to " << (MaxSteps+1) << "." << endl;
-    Trajectory.R.resize(MaxSteps+1);
-    Trajectory.U.resize(MaxSteps+1);
-    Trajectory.F.resize(MaxSteps+1);
-  }
-};
-
-/** Copies a given trajectory step \a src onto another \a dest
- * \param dest index of destination step
- * \param src index of source step
- */
-void TrajectoryParticle::CopyStepOnStep(int dest, int src)
-{
-  if (dest == src)  // self assignment check
-    return;
-
-  for (int n=NDIM;n--;) {
-    Trajectory.R.at(dest).x[n] = Trajectory.R.at(src).x[n];
-    Trajectory.U.at(dest).x[n] = Trajectory.U.at(src).x[n];
-    Trajectory.F.at(dest).x[n] = Trajectory.F.at(src).x[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 TrajectoryParticle::VelocityVerletUpdate(int NextStep, config *configuration, ForceMatrix *Force)
-{
-  //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
-  for (int d=0; d<NDIM; d++) {
-    Trajectory.F.at(NextStep).x[d] = -Force->Matrix[0][nr][d+5]*(configuration->GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
-    Trajectory.R.at(NextStep).x[d] = Trajectory.R.at(NextStep-1).x[d];
-    Trajectory.R.at(NextStep).x[d] += configuration->Deltat*(Trajectory.U.at(NextStep-1).x[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
-    Trajectory.R.at(NextStep).x[d] += 0.5*configuration->Deltat*configuration->Deltat*(Trajectory.F.at(NextStep).x[d]/type->mass);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
-  }
-  // Update U
-  for (int d=0; d<NDIM; d++) {
-    Trajectory.U.at(NextStep).x[d] = Trajectory.U.at(NextStep-1).x[d];
-    Trajectory.U.at(NextStep).x[d] += configuration->Deltat * (Trajectory.F.at(NextStep).x[d]+Trajectory.F.at(NextStep-1).x[d]/type->mass); // v = F/m * t
-  }
-  // Update R (and F)
-//      out << "Integrated position&velocity of step " << (NextStep) << ": (";
-//      for (int d=0;d<NDIM;d++)
-//        out << Trajectory.R.at(NextStep).x[d] << " ";          // next step
-//      out << ")\t(";
-//      for (int d=0;d<NDIM;d++)
-//        cout << Trajectory.U.at(NextStep).x[d] << " ";          // next step
-//      out << ")" << endl;
-};
-
-/** Sums up mass and kinetics.
- * \param Step step to sum for
- * \param *TotalMass pointer to total mass sum
- * \param *TotalVelocity pointer to tota velocity sum
- */
-void TrajectoryParticle::SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity )
-{
-  *TotalMass += type->mass;  // sum up total mass
-  for(int d=0;d<NDIM;d++) {
-    TotalVelocity->x[d] += Trajectory.U.at(Step).x[d]*type->mass;
-  }
-};
-
-/** Scales velocity of atom according to Woodcock thermostat.
- * \param ScaleTempFactor factor to scale the velocities with (i.e. sqrt of energy scale factor)
- * \param Step MD step to scale
- * \param *ekin sum of kinetic energy
- */
-void TrajectoryParticle::Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin)
-{
-  double *U = Trajectory.U.at(Step).x;
-  if (FixedIon == 0) // even FixedIon moves, only not by other's forces
-    for (int d=0; d<NDIM; d++) {
-      U[d] *= ScaleTempFactor;
-      *ekin += 0.5*type->mass * U[d]*U[d];
-    }
-};
-
-/** Scales velocity of atom according to Gaussian thermostat.
- * \param Step MD step to scale
- * \param *G
- * \param *E
- */
-void TrajectoryParticle::Thermostat_Gaussian_init(int Step, double *G, double *E)
-{
-  double *U = Trajectory.U.at(Step).x;
-  double *F = Trajectory.F.at(Step).x;
-  if (FixedIon == 0) // even FixedIon moves, only not by other's forces
-    for (int d=0; d<NDIM; d++) {
-      *G += U[d] * F[d];
-      *E += U[d]*U[d]*type->mass;
-    }
-};
-
-/** Determines scale factors according to Gaussian thermostat.
- * \param Step MD step to scale
- * \param GE G over E ratio
- * \param *ekin sum of kinetic energy
- * \param *configuration configuration class with TempFrequency and TargetTemp
- */
-void TrajectoryParticle::Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration)
-{
-  double *U = Trajectory.U.at(Step).x;
-  if (FixedIon == 0) // even FixedIon moves, only not by other's forces
-    for (int d=0; d<NDIM; d++) {
-      U[d] += configuration->Deltat/type->mass * ( (G_over_E) * (U[d]*type->mass) );
-      *ekin += type->mass * U[d]*U[d];
-    }
-};
-
-/** Scales velocity of atom according to Langevin thermostat.
- * \param Step MD step to scale
- * \param *r random number generator
- * \param *ekin sum of kinetic energy
- * \param *configuration configuration class with TempFrequency and TargetTemp
- */
-void TrajectoryParticle::Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration)
-{
-  double sigma  = sqrt(configuration->TargetTemp/type->mass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
-  double *U = Trajectory.U.at(Step).x;
-  if (FixedIon == 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 (((((rand()/(double)RAND_MAX))*configuration->TempFrequency) < 1.)) {
-      cout << Verbose(3) << "Particle " << *this << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> ";
-      // 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);
-      }
-      cout << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << endl;
-    }
-    for (int d=0; d<NDIM; d++)
-      *ekin += 0.5*type->mass * U[d]*U[d];
-  }
-};
-
-/** Scales velocity of atom according to Berendsen thermostat.
- * \param Step MD step to scale
- * \param ScaleTempFactor factor to scale energy (not velocity!) with
- * \param *ekin sum of kinetic energy
- * \param *configuration configuration class with TempFrequency and Deltat
- */
-void TrajectoryParticle::Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration)
-{
-  double *U = Trajectory.U.at(Step).x;
-  if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
-    for (int d=0; d<NDIM; d++) {
-      U[d] *= sqrt(1+(configuration->Deltat/configuration->TempFrequency)*(ScaleTempFactor-1));
-      *ekin += 0.5*type->mass * U[d]*U[d];
-    }
-  }
-};
-
-/** Initializes current run of NoseHoover thermostat.
- * \param Step MD step to scale
- * \param *delta_alpha additional sum of kinetic energy on return
- */
-void TrajectoryParticle::Thermostat_NoseHoover_init(int Step, double *delta_alpha)
-{
-  double *U = Trajectory.U.at(Step).x;
-  if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
-    for (int d=0; d<NDIM; d++) {
-      *delta_alpha += U[d]*U[d]*type->mass;
-    }
-  }
-};
-
-/** Initializes current run of NoseHoover thermostat.
- * \param Step MD step to scale
- * \param *ekin sum of kinetic energy
- * \param *configuration configuration class with TempFrequency and Deltat
- */
-void TrajectoryParticle::Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration)
-{
-  double *U = Trajectory.U.at(Step).x;
-  if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
-    for (int d=0; d<NDIM; d++) {
-        U[d] += configuration->Deltat/type->mass * (configuration->alpha * (U[d] * type->mass));
-        *ekin += (0.5*type->mass) * U[d]*U[d];
-      }
-  }
-};
Index: molecuilder/src/atom.hpp
===================================================================
--- molecuilder/src/atom.hpp	(revision 3b0ba203af31dbee5c5a38092340126009136ed1)
+++ molecuilder/src/atom.hpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -22,142 +22,16 @@
 #include <vector>
 
-#include <gsl/gsl_randist.h>
-
+#include "atom_atominfo.hpp"
+#include "atom_bondedparticle.hpp"
+#include "atom_graphnode.hpp"
+#include "atom_particleinfo.hpp"
+#include "atom_trajectoryparticle.hpp"
 #include "tesselation.hpp"
 
 /****************************************** forward declarations *****************************/
 
-class atom;
-class bond;
-class config;
-class element;
-class ForceMatrix;
 class Vector;
 
-#define BondList list<bond *>
-
 /********************************************** declarations *******************************/
-
-class AtomInfo {
-public:
-  Vector x;       //!< coordinate vector of atom, giving last position within cell
-  Vector v;       //!< velocity vector of atom, giving last velocity within cell
-  Vector F;       //!< Force vector of atom, giving last force within cell
-  element *type;  //!< pointing to element
-
-  AtomInfo();
-  ~AtomInfo();
-
-private:
-};
-
-class TrajectoryParticleInfo {
-public:
-  struct
-  {
-    vector<Vector> R;  //!< position vector
-    vector<Vector> U;  //!< velocity vector
-    vector<Vector> F;  //!< last force vector
-  } Trajectory;
-  int FixedIon;   //!< config variable that states whether forces act on the ion or not
-
-  TrajectoryParticleInfo();
-  ~TrajectoryParticleInfo();
-
-private:
-};
-
-class TrajectoryParticle :  public TrajectoryParticleInfo, public virtual AtomInfo, public virtual ParticleInfo {
-public:
-
-  TrajectoryParticle();
-  virtual ~TrajectoryParticle();
-
-  // constraint potential and dynamics stuff
-  void AddKineticToTemperature(double *temperature, int step) const;
-  void EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
-  void CorrectVelocity(double *ActualTemp, int Step, Vector *CoGVelocity);
-
-  // trajectory stuff
-  void ResizeTrajectory(int MaxSteps);
-  void CopyStepOnStep(int dest, int src);
-  void VelocityVerletUpdate(int MDSteps, config *configuration, ForceMatrix *Force);
-  void SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity );
-
-  // thermostats
-  void Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin);
-  void Thermostat_Gaussian_init(int Step, double *G, double *E);
-  void Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration);
-  void Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration);
-  void Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration);
-  void Thermostat_NoseHoover_init(int Step, double *delta_alpha);
-  void Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration);
-
-private:
-
-};
-
-class GraphNode;
-
-class GraphNodeInfo {
-public:
-
-  int GraphNr;      //!< unique number, given in DepthFirstSearchAnalysis()
-  int *ComponentNr;//!< belongs to this non-separable components, given in DepthFirstSearchAnalysis() (if more than one, then is SeparationVertex)
-  int LowpointNr; //!< needed in DepthFirstSearchAnalysis() to detect non-separable components, is the lowest possible number of an atom to reach via tree edges only followed by at most one back edge.
-  bool SeparationVertex; //!< whether this atom separates off subsets of atoms or not, determined in DepthFirstSearchAnalysis()
-  bool IsCyclic;        //!< whether atom belong to as cycle or not, determined in DepthFirstSearchAnalysis()
-  atom *Ancestor; //!< "Father" in Depth-First-Search
-
-  GraphNodeInfo();
-  ~GraphNodeInfo();
-private:
-
-};
-
-
-class GraphNode : public GraphNodeInfo, public virtual ParticleInfo {
-public:
-
-  GraphNode();
-  virtual ~GraphNode();
-
-  void OutputGraphInfo(ofstream *out) const;
-  void OutputComponentNumber(ofstream *out) const;
-
-private:
-
-};
-
-class BondedParticleInfo {
-public:
-  unsigned char AdaptiveOrder;  //!< current present bond order at site (0 means "not set")
-  bool MaxOrder;  //!< whether this atom as a root in fragmentation still creates more fragments on higher orders or not
-  BondList ListOfBonds; //!< list of all bonds
-
-  BondedParticleInfo();
-  ~BondedParticleInfo();
-
-private:
-
-};
-
-class BondedParticle : public BondedParticleInfo, public virtual ParticleInfo, public virtual AtomInfo {
-public:
-  BondedParticle();
-  virtual ~BondedParticle();
-
-  bool RegisterBond(bond *Binder);
-  bool UnregisterBond(bond *Binder);
-  void UnregisterAllBond();
-  int CountBonds() const;
-  int CorrectBondDegree(ofstream *out);
-  bool OutputBondOfAtom(ofstream *out) const;
-  void OutputAdjacency(ofstream *AdjacencyFile) const;
-  void OutputOrder(ofstream *file);
-
-private:
-
-};
 
 /** Single atom.
Index: molecuilder/src/atom_atominfo.cpp
===================================================================
--- molecuilder/src/atom_atominfo.cpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
+++ molecuilder/src/atom_atominfo.cpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -0,0 +1,19 @@
+/*
+ * atom_atominfo.cpp
+ *
+ *  Created on: Oct 19, 2009
+ *      Author: heber
+ */
+
+#include "atom_atominfo.hpp"
+
+/** Constructor of class AtomInfo.
+ */
+AtomInfo::AtomInfo() : type(NULL) {};
+
+/** Destructor of class AtomInfo.
+ */
+AtomInfo::~AtomInfo()
+{
+};
+
Index: molecuilder/src/atom_atominfo.hpp
===================================================================
--- molecuilder/src/atom_atominfo.hpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
+++ molecuilder/src/atom_atominfo.hpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -0,0 +1,42 @@
+/*
+ * atom_atominfo.hpp
+ *
+ *  Created on: Oct 19, 2009
+ *      Author: heber
+ */
+
+#ifndef ATOM_ATOMINFO_HPP_
+#define ATOM_ATOMINFO_HPP_
+
+
+using namespace std;
+
+/*********************************************** includes ***********************************/
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "vector.hpp"
+
+/****************************************** forward declarations *****************************/
+
+class element;
+
+/********************************************** declarations *******************************/
+
+class AtomInfo {
+public:
+  Vector x;       //!< coordinate vector of atom, giving last position within cell
+  Vector v;       //!< velocity vector of atom, giving last velocity within cell
+  Vector F;       //!< Force vector of atom, giving last force within cell
+  element *type;  //!< pointing to element
+
+  AtomInfo();
+  ~AtomInfo();
+
+private:
+};
+
+#endif /* ATOM_ATOMINFO_HPP_ */
Index: molecuilder/src/atom_bondedparticle.cpp
===================================================================
--- molecuilder/src/atom_bondedparticle.cpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
+++ molecuilder/src/atom_bondedparticle.cpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -0,0 +1,150 @@
+/*
+ * atom_bondedparticle.cpp
+ *
+ *  Created on: Oct 19, 2009
+ *      Author: heber
+ */
+
+#include "atom.hpp"
+#include "atom_bondedparticle.hpp"
+#include "bond.hpp"
+#include "element.hpp"
+#include "lists.hpp"
+#include "verbose.hpp"
+
+/** Constructor of class BondedParticle.
+ */
+BondedParticle::BondedParticle(){};
+
+/** Destructor of class BondedParticle.
+ */
+BondedParticle::~BondedParticle()
+{
+  BondList::const_iterator Runner;
+  while (!ListOfBonds.empty()) {
+    Runner = ListOfBonds.begin();
+    removewithoutcheck(*Runner);
+  }
+};
+
+/** Outputs the current atom::AdaptiveOrder and atom::MaxOrder to \a *file.
+ * \param *file output stream
+ */
+void BondedParticle::OutputOrder(ofstream *file)
+{
+  *file << nr << "\t" << (int)AdaptiveOrder << "\t" << (int)MaxOrder << endl;
+  //cout << Verbose(2) << "Storing: " << nr << "\t" << (int)AdaptiveOrder << "\t" << (int)MaxOrder << "." << endl;
+};
+
+/** Prints all bonds of this atom with total degree.
+ * \param *out stream to output to
+ * \return true - \a *out present, false - \a *out is NULL
+ */
+bool BondedParticle::OutputBondOfAtom(ofstream *out) const
+{
+  if (out != NULL) {
+    *out << Verbose(4) << "Atom " << Name << "/" << nr << " with " << ListOfBonds.size() << " bonds: ";
+    int TotalDegree = 0;
+    for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); ++Runner) {
+      *out << **Runner << "\t";
+      TotalDegree += (*Runner)->BondDegree;
+    }
+    *out << " -- TotalDegree: " << TotalDegree << endl;
+    return true;
+  } else
+    return false;
+};
+
+/** Output of atom::nr along with all bond partners.
+ * \param *AdjacencyFile output stream
+ */
+void BondedParticle::OutputAdjacency(ofstream *AdjacencyFile) const
+{
+  *AdjacencyFile << nr << "\t";
+  for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner))
+    *AdjacencyFile << (*Runner)->GetOtherAtom(this)->nr << "\t";
+  *AdjacencyFile << endl;
+};
+
+/** Puts a given bond into atom::ListOfBonds.
+ * \param *Binder bond to insert
+ */
+bool BondedParticle::RegisterBond(bond *Binder)
+{
+  bool status = false;
+  if (Binder != NULL) {
+    if (Binder->Contains(this)) {
+      ListOfBonds.push_back(Binder);
+      status = true;
+    } else {
+      cout << Verbose(1) << "ERROR: " << *Binder << " does not contain " << *this << "." << endl;
+    }
+  } else {
+    cout << Verbose(1) << "ERROR: Binder is " << Binder << "." << endl;
+  }
+  return status;
+};
+
+/** Removes a given bond from atom::ListOfBonds.
+ * \param *Binder bond to remove
+ */
+bool BondedParticle::UnregisterBond(bond *Binder)
+{
+  bool status = false;
+  if (Binder != NULL) {
+    if (Binder->Contains(this)) {
+      ListOfBonds.remove(Binder);
+      status = true;
+    } else {
+      cout << Verbose(1) << "ERROR: " << *Binder << " does not contain " << *this << "." << endl;
+    }
+  } else {
+    cout << Verbose(1) << "ERROR: Binder is " << Binder << "." << endl;
+  }
+  return status;
+};
+
+/** Removes all bonds from atom::ListOfBonds.
+ * \note Does not do any memory de-allocation.
+ */
+void BondedParticle::UnregisterAllBond()
+{
+  ListOfBonds.clear();
+};
+
+/** Corrects the bond degree by one at most if necessary.
+ * \param *out output stream for debugging
+ */
+int BondedParticle::CorrectBondDegree(ofstream *out)
+{
+  int NoBonds = 0;
+  int OtherNoBonds = 0;
+  int FalseBondDegree = 0;
+  atom *OtherWalker = NULL;
+  bond *CandidateBond = NULL;
+
+  *out << Verbose(3) << "Walker " << *this << ": " << (int)this->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
+  NoBonds = CountBonds();
+  if ((int)(type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check all bonding partners for mismatch
+    for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner)) {
+      OtherWalker = (*Runner)->GetOtherAtom(this);
+      OtherNoBonds = OtherWalker->CountBonds();
+      *out << Verbose(3) << "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
+      if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) { // check if possible candidate
+        if ((CandidateBond == NULL) || (ListOfBonds.size() > OtherWalker->ListOfBonds.size())) { // pick the one with fewer number of bonds first
+          CandidateBond = (*Runner);
+          *out << Verbose(3) << "New candidate is " << *CandidateBond << "." << endl;
+        }
+      }
+    }
+    if ((CandidateBond != NULL)) {
+      CandidateBond->BondDegree++;
+      *out << Verbose(2) << "Increased bond degree for bond " << *CandidateBond << "." << endl;
+    } else {
+      *out << Verbose(2) << "Could not find correct degree for atom " << *this << "." << endl;
+      FalseBondDegree++;
+    }
+  }
+  return FalseBondDegree;
+};
+
Index: molecuilder/src/atom_bondedparticle.hpp
===================================================================
--- molecuilder/src/atom_bondedparticle.hpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
+++ molecuilder/src/atom_bondedparticle.hpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -0,0 +1,53 @@
+/*
+ * atom_bondedparticle.hpp
+ *
+ *  Created on: Oct 19, 2009
+ *      Author: heber
+ */
+
+#ifndef ATOM_BONDEDPARTICLE_HPP_
+#define ATOM_BONDEDPARTICLE_HPP_
+
+
+using namespace std;
+
+/*********************************************** includes ***********************************/
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <fstream>
+#include <list>
+
+#include "atom_atominfo.hpp"
+#include "atom_bondedparticleinfo.hpp"
+#include "atom_particleinfo.hpp"
+
+/****************************************** forward declarations *****************************/
+
+class bond;
+
+/********************************************** declarations *******************************/
+
+class BondedParticle : public BondedParticleInfo, public virtual ParticleInfo, public virtual AtomInfo {
+public:
+  BondedParticle();
+  virtual ~BondedParticle();
+
+  bool RegisterBond(bond *Binder);
+  bool UnregisterBond(bond *Binder);
+  void UnregisterAllBond();
+  int CountBonds() const;
+  int CorrectBondDegree(ofstream *out);
+  bool OutputBondOfAtom(ofstream *out) const;
+  void OutputAdjacency(ofstream *AdjacencyFile) const;
+  void OutputOrder(ofstream *file);
+
+private:
+
+};
+
+
+#endif /* ATOM_BONDEDPARTICLE_HPP_ */
Index: molecuilder/src/atom_bondedparticleinfo.cpp
===================================================================
--- molecuilder/src/atom_bondedparticleinfo.cpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
+++ molecuilder/src/atom_bondedparticleinfo.cpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -0,0 +1,19 @@
+/*
+ * atom_bondedparticleinfo.cpp
+ *
+ *  Created on: Oct 19, 2009
+ *      Author: heber
+ */
+
+#include "atom_bondedparticleinfo.hpp"
+
+/** Constructor of class BondedParticleInfo.
+ */
+BondedParticleInfo::BondedParticleInfo() : AdaptiveOrder(0), MaxOrder(false) {};
+
+/** Destructor of class BondedParticleInfo.
+ */
+BondedParticleInfo::~BondedParticleInfo()
+{
+};
+
Index: molecuilder/src/atom_bondedparticleinfo.hpp
===================================================================
--- molecuilder/src/atom_bondedparticleinfo.hpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
+++ molecuilder/src/atom_bondedparticleinfo.hpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -0,0 +1,45 @@
+/*
+ * atom_bondedparticleinfo.hpp
+ *
+ *  Created on: Oct 19, 2009
+ *      Author: heber
+ */
+
+#ifndef ATOM_BONDEDPARTICLEINFO_HPP_
+#define ATOM_BONDEDPARTICLEINFO_HPP_
+
+
+using namespace std;
+
+/*********************************************** includes ***********************************/
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <list>
+
+/****************************************** forward declarations *****************************/
+
+class bond;
+
+#define BondList list<bond *>
+
+/********************************************** declarations *******************************/
+
+class BondedParticleInfo {
+public:
+  unsigned char AdaptiveOrder;  //!< current present bond order at site (0 means "not set")
+  bool MaxOrder;  //!< whether this atom as a root in fragmentation still creates more fragments on higher orders or not
+  BondList ListOfBonds; //!< list of all bonds
+
+  BondedParticleInfo();
+  ~BondedParticleInfo();
+
+private:
+
+};
+
+
+#endif /* ATOM_BONDEDPARTICLEINFO_HPP_ */
Index: molecuilder/src/atom_graphnode.cpp
===================================================================
--- molecuilder/src/atom_graphnode.cpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
+++ molecuilder/src/atom_graphnode.cpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -0,0 +1,44 @@
+/*
+ * atom_graphnode.cpp
+ *
+ *  Created on: Oct 19, 2009
+ *      Author: heber
+ */
+
+#include "atom_graphnode.hpp"
+#include "verbose.hpp"
+
+/** Constructor of class GraphNode.
+ */
+GraphNode::GraphNode()
+{
+};
+
+/** Destructor of class GraphNode.
+ */
+GraphNode::~GraphNode()
+{
+};
+
+/** Output graph info of this atom.
+ * \param *out output stream
+ */
+void GraphNode::OutputGraphInfo(ofstream *out) const
+{
+  *out << Verbose(2) << "Atom " << Name << " is " << ((SeparationVertex) ? "a" : "not a") << " separation vertex, components are ";
+  OutputComponentNumber(out);
+  *out << " with Lowpoint " << LowpointNr << " and Graph Nr. " << GraphNr << "." << endl;
+};
+
+/** Output a list of flags, stating whether the bond was visited or not.
+ * Note, we make use of the last entry of the ComponentNr always being -1 if allocated.
+ * \param *out output stream for debugging
+ */
+void GraphNode::OutputComponentNumber(ofstream *out) const
+{
+  if (ComponentNr != NULL) {
+    for (int i=0; ComponentNr[i] != -1; i++)
+      *out << ComponentNr[i] << " ";
+  }
+};
+
Index: molecuilder/src/atom_graphnode.hpp
===================================================================
--- molecuilder/src/atom_graphnode.hpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
+++ molecuilder/src/atom_graphnode.hpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -0,0 +1,43 @@
+/*
+ * atom_graphnode.hpp
+ *
+ *  Created on: Oct 19, 2009
+ *      Author: heber
+ */
+
+#ifndef ATOM_GRAPHNODE_HPP_
+#define ATOM_GRAPHNODE_HPP_
+
+
+using namespace std;
+
+/*********************************************** includes ***********************************/
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <fstream>
+
+#include "atom_graphnodeinfo.hpp"
+#include "atom_particleinfo.hpp"
+
+/****************************************** forward declarations *****************************/
+
+/********************************************** declarations *******************************/
+
+class GraphNode : public GraphNodeInfo, public virtual ParticleInfo {
+public:
+
+  GraphNode();
+  virtual ~GraphNode();
+
+  void OutputGraphInfo(ofstream *out) const;
+  void OutputComponentNumber(ofstream *out) const;
+
+private:
+
+};
+
+#endif /* ATOM_GRAPHNODE_HPP_ */
Index: molecuilder/src/atom_graphnodeinfo.cpp
===================================================================
--- molecuilder/src/atom_graphnodeinfo.cpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
+++ molecuilder/src/atom_graphnodeinfo.cpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -0,0 +1,20 @@
+/*
+ * atom_graphnodeinfo.cpp
+ *
+ *  Created on: Oct 19, 2009
+ *      Author: heber
+ */
+
+#include "atom_graphnodeinfo.hpp"
+#include "memoryallocator.hpp"
+
+/** Constructor of class GraphNodeInfo.
+ */
+GraphNodeInfo::GraphNodeInfo() : GraphNr(-1), ComponentNr(NULL), LowpointNr(-1), SeparationVertex(false), IsCyclic(false), Ancestor(NULL) {};
+
+/** Destructor of class GraphNodeInfo.
+ */
+GraphNodeInfo::~GraphNodeInfo()
+{
+  Free<int>(&ComponentNr, "atom::~atom: *ComponentNr");
+};
Index: molecuilder/src/atom_graphnodeinfo.hpp
===================================================================
--- molecuilder/src/atom_graphnodeinfo.hpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
+++ molecuilder/src/atom_graphnodeinfo.hpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -0,0 +1,44 @@
+/*
+ * atom_graphnodeinfo.hpp
+ *
+ *  Created on: Oct 19, 2009
+ *      Author: heber
+ */
+
+#ifndef ATOM_GRAPHNODEINFO_HPP_
+#define ATOM_GRAPHNODEINFO_HPP_
+
+using namespace std;
+
+/*********************************************** includes ***********************************/
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+/****************************************** forward declarations *****************************/
+
+class atom;
+
+/********************************************** declarations *******************************/
+
+class GraphNode;
+
+class GraphNodeInfo {
+public:
+
+  int GraphNr;      //!< unique number, given in DepthFirstSearchAnalysis()
+  int *ComponentNr;//!< belongs to this non-separable components, given in DepthFirstSearchAnalysis() (if more than one, then is SeparationVertex)
+  int LowpointNr; //!< needed in DepthFirstSearchAnalysis() to detect non-separable components, is the lowest possible number of an atom to reach via tree edges only followed by at most one back edge.
+  bool SeparationVertex; //!< whether this atom separates off subsets of atoms or not, determined in DepthFirstSearchAnalysis()
+  bool IsCyclic;        //!< whether atom belong to as cycle or not, determined in DepthFirstSearchAnalysis()
+  atom *Ancestor; //!< "Father" in Depth-First-Search
+
+  GraphNodeInfo();
+  ~GraphNodeInfo();
+private:
+
+};
+
+#endif /* ATOM_GRAPHNODEINFO_HPP_ */
Index: molecuilder/src/atom_particleinfo.cpp
===================================================================
--- molecuilder/src/atom_particleinfo.cpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
+++ molecuilder/src/atom_particleinfo.cpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -0,0 +1,21 @@
+/*
+ * atom_particleinfo.cpp
+ *
+ *  Created on: Oct 19, 2009
+ *      Author: heber
+ */
+
+#include "atom_particleinfo.hpp"
+#include "memoryallocator.hpp"
+
+/** Constructor of ParticleInfo.
+ */
+ParticleInfo::ParticleInfo() : nr(-1), Name(NULL) {};
+
+/** Destructor of ParticleInfo.
+ */
+ParticleInfo::~ParticleInfo()
+{
+  Free(&Name);
+};
+
Index: molecuilder/src/atom_particleinfo.hpp
===================================================================
--- molecuilder/src/atom_particleinfo.hpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
+++ molecuilder/src/atom_particleinfo.hpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -0,0 +1,42 @@
+/*
+ * atom_particleinfo.hpp
+ *
+ *  Created on: Oct 19, 2009
+ *      Author: heber
+ */
+
+#ifndef ATOM_PARTICLEINFO_HPP_
+#define ATOM_PARTICLEINFO_HPP_
+
+using namespace std;
+
+/*********************************************** includes ***********************************/
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <iostream>
+
+/****************************************** forward declarations *****************************/
+
+
+/********************************************** declarations *******************************/
+
+class ParticleInfo {
+public:
+  int nr;       // index to easierly identify
+  char *Name;   // some name to reference to on output
+
+  ParticleInfo();
+  ~ParticleInfo();
+
+  ostream & operator << (ostream &ost);
+
+private:
+};
+
+ostream & operator << (ostream &ost, const ParticleInfo &a);
+
+#endif /* ATOM_PARTICLEINFO_HPP_ */
Index: molecuilder/src/atom_trajectoryparticle.cpp
===================================================================
--- molecuilder/src/atom_trajectoryparticle.cpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
+++ molecuilder/src/atom_trajectoryparticle.cpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -0,0 +1,256 @@
+/*
+ * atom_trajectoryparticle.cpp
+ *
+ *  Created on: Oct 19, 2009
+ *      Author: heber
+ */
+
+#include "atom.hpp"
+#include "atom_trajectoryparticle.hpp"
+#include "config.hpp"
+#include "element.hpp"
+#include "parser.hpp"
+#include "verbose.hpp"
+
+/** Constructor of class TrajectoryParticle.
+ */
+TrajectoryParticle::TrajectoryParticle()
+{
+};
+
+/** Destructor of class TrajectoryParticle.
+ */
+TrajectoryParticle::~TrajectoryParticle()
+{
+};
+
+
+/** Adds kinetic energy of this atom to given temperature value.
+ * \param *temperature add on this value
+ * \param step given step of trajectory to add
+ */
+void TrajectoryParticle::AddKineticToTemperature(double *temperature, int step) const
+{
+  for (int i=NDIM;i--;)
+    *temperature += type->mass * Trajectory.U.at(step).x[i]* Trajectory.U.at(step).x[i];
+};
+
+/** Evaluates some constraint potential if atom moves from \a startstep at once to \endstep in trajectory.
+ * \param startstep trajectory begins at
+ * \param endstep trajectory ends at
+ * \param **PermutationMap if atom switches places with some other atom, there is no translation but a permutaton noted here (not in the trajectories of ea
+ * \param *Force Force matrix to store result in
+ */
+void TrajectoryParticle::EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force)
+{
+  double constant = 10.;
+  TrajectoryParticle *Sprinter = PermutationMap[nr];
+  // set forces
+  for (int i=NDIM;i++;)
+    Force->Matrix[0][nr][5+i] += 2.*constant*sqrt(Trajectory.R.at(startstep).Distance(&Sprinter->Trajectory.R.at(endstep)));
+};
+
+/** Correct velocity against the summed \a CoGVelocity for \a step.
+ * \param *ActualTemp sum up actual temperature meanwhile
+ * \param Step MD step in atom::Tracjetory
+ * \param *CoGVelocity remnant velocity (i.e. vector sum of all atom velocities)
+ */
+void TrajectoryParticle::CorrectVelocity(double *ActualTemp, int Step, Vector *CoGVelocity)
+{
+  for(int d=0;d<NDIM;d++) {
+    Trajectory.U.at(Step).x[d] -= CoGVelocity->x[d];
+    *ActualTemp += 0.5 * type->mass * Trajectory.U.at(Step).x[d] * Trajectory.U.at(Step).x[d];
+  }
+};
+
+/** Extends the trajectory STL vector to the new size.
+ * Does nothing if \a MaxSteps is smaller than current size.
+ * \param MaxSteps
+ */
+void TrajectoryParticle::ResizeTrajectory(int MaxSteps)
+{
+  if (Trajectory.R.size() <= (unsigned int)(MaxSteps)) {
+    //cout << "Increasing size for trajectory array of " << keyword << " to " << (MaxSteps+1) << "." << endl;
+    Trajectory.R.resize(MaxSteps+1);
+    Trajectory.U.resize(MaxSteps+1);
+    Trajectory.F.resize(MaxSteps+1);
+  }
+};
+
+/** Copies a given trajectory step \a src onto another \a dest
+ * \param dest index of destination step
+ * \param src index of source step
+ */
+void TrajectoryParticle::CopyStepOnStep(int dest, int src)
+{
+  if (dest == src)  // self assignment check
+    return;
+
+  for (int n=NDIM;n--;) {
+    Trajectory.R.at(dest).x[n] = Trajectory.R.at(src).x[n];
+    Trajectory.U.at(dest).x[n] = Trajectory.U.at(src).x[n];
+    Trajectory.F.at(dest).x[n] = Trajectory.F.at(src).x[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 TrajectoryParticle::VelocityVerletUpdate(int NextStep, config *configuration, ForceMatrix *Force)
+{
+  //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
+  for (int d=0; d<NDIM; d++) {
+    Trajectory.F.at(NextStep).x[d] = -Force->Matrix[0][nr][d+5]*(configuration->GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
+    Trajectory.R.at(NextStep).x[d] = Trajectory.R.at(NextStep-1).x[d];
+    Trajectory.R.at(NextStep).x[d] += configuration->Deltat*(Trajectory.U.at(NextStep-1).x[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
+    Trajectory.R.at(NextStep).x[d] += 0.5*configuration->Deltat*configuration->Deltat*(Trajectory.F.at(NextStep).x[d]/type->mass);     // F = m * a and s =
+  }
+  // Update U
+  for (int d=0; d<NDIM; d++) {
+    Trajectory.U.at(NextStep).x[d] = Trajectory.U.at(NextStep-1).x[d];
+    Trajectory.U.at(NextStep).x[d] += configuration->Deltat * (Trajectory.F.at(NextStep).x[d]+Trajectory.F.at(NextStep-1).x[d]/type->mass); // v = F/m * t
+  }
+  // Update R (and F)
+//      out << "Integrated position&velocity of step " << (NextStep) << ": (";
+//      for (int d=0;d<NDIM;d++)
+//        out << Trajectory.R.at(NextStep).x[d] << " ";          // next step
+//      out << ")\t(";
+//      for (int d=0;d<NDIM;d++)
+//        cout << Trajectory.U.at(NextStep).x[d] << " ";          // next step
+//      out << ")" << endl;
+};
+
+/** Sums up mass and kinetics.
+ * \param Step step to sum for
+ * \param *TotalMass pointer to total mass sum
+ * \param *TotalVelocity pointer to tota velocity sum
+ */
+void TrajectoryParticle::SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity )
+{
+  *TotalMass += type->mass;  // sum up total mass
+  for(int d=0;d<NDIM;d++) {
+    TotalVelocity->x[d] += Trajectory.U.at(Step).x[d]*type->mass;
+  }
+};
+
+/** Scales velocity of atom according to Woodcock thermostat.
+ * \param ScaleTempFactor factor to scale the velocities with (i.e. sqrt of energy scale factor)
+ * \param Step MD step to scale
+ * \param *ekin sum of kinetic energy
+ */
+void TrajectoryParticle::Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin)
+{
+  double *U = Trajectory.U.at(Step).x;
+  if (FixedIon == 0) // even FixedIon moves, only not by other's forces
+    for (int d=0; d<NDIM; d++) {
+      U[d] *= ScaleTempFactor;
+      *ekin += 0.5*type->mass * U[d]*U[d];
+    }
+};
+
+/** Scales velocity of atom according to Gaussian thermostat.
+ * \param Step MD step to scale
+ * \param *G
+ * \param *E
+ */
+void TrajectoryParticle::Thermostat_Gaussian_init(int Step, double *G, double *E)
+{
+  double *U = Trajectory.U.at(Step).x;
+  double *F = Trajectory.F.at(Step).x;
+  if (FixedIon == 0) // even FixedIon moves, only not by other's forces
+    for (int d=0; d<NDIM; d++) {
+      *G += U[d] * F[d];
+      *E += U[d]*U[d]*type->mass;
+    }
+};
+
+/** Determines scale factors according to Gaussian thermostat.
+ * \param Step MD step to scale
+ * \param GE G over E ratio
+ * \param *ekin sum of kinetic energy
+ * \param *configuration configuration class with TempFrequency and TargetTemp
+ */
+void TrajectoryParticle::Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration)
+{
+  double *U = Trajectory.U.at(Step).x;
+  if (FixedIon == 0) // even FixedIon moves, only not by other's forces
+    for (int d=0; d<NDIM; d++) {
+      U[d] += configuration->Deltat/type->mass * ( (G_over_E) * (U[d]*type->mass) );
+      *ekin += type->mass * U[d]*U[d];
+    }
+};
+
+/** Scales velocity of atom according to Langevin thermostat.
+ * \param Step MD step to scale
+ * \param *r random number generator
+ * \param *ekin sum of kinetic energy
+ * \param *configuration configuration class with TempFrequency and TargetTemp
+ */
+void TrajectoryParticle::Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration)
+{
+  double sigma  = sqrt(configuration->TargetTemp/type->mass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
+  double *U = Trajectory.U.at(Step).x;
+  if (FixedIon == 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 (((((rand()/(double)RAND_MAX))*configuration->TempFrequency) < 1.)) {
+      cout << Verbose(3) << "Particle " << *this << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> ";
+      // 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);
+      }
+      cout << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << endl;
+    }
+    for (int d=0; d<NDIM; d++)
+      *ekin += 0.5*type->mass * U[d]*U[d];
+  }
+};
+
+/** Scales velocity of atom according to Berendsen thermostat.
+ * \param Step MD step to scale
+ * \param ScaleTempFactor factor to scale energy (not velocity!) with
+ * \param *ekin sum of kinetic energy
+ * \param *configuration configuration class with TempFrequency and Deltat
+ */
+void TrajectoryParticle::Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration)
+{
+  double *U = Trajectory.U.at(Step).x;
+  if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
+    for (int d=0; d<NDIM; d++) {
+      U[d] *= sqrt(1+(configuration->Deltat/configuration->TempFrequency)*(ScaleTempFactor-1));
+      *ekin += 0.5*type->mass * U[d]*U[d];
+    }
+  }
+};
+
+/** Initializes current run of NoseHoover thermostat.
+ * \param Step MD step to scale
+ * \param *delta_alpha additional sum of kinetic energy on return
+ */
+void TrajectoryParticle::Thermostat_NoseHoover_init(int Step, double *delta_alpha)
+{
+  double *U = Trajectory.U.at(Step).x;
+  if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
+    for (int d=0; d<NDIM; d++) {
+      *delta_alpha += U[d]*U[d]*type->mass;
+    }
+  }
+};
+
+/** Initializes current run of NoseHoover thermostat.
+ * \param Step MD step to scale
+ * \param *ekin sum of kinetic energy
+ * \param *configuration configuration class with TempFrequency and Deltat
+ */
+void TrajectoryParticle::Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration)
+{
+  double *U = Trajectory.U.at(Step).x;
+  if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
+    for (int d=0; d<NDIM; d++) {
+        U[d] += configuration->Deltat/type->mass * (configuration->alpha * (U[d] * type->mass));
+        *ekin += (0.5*type->mass) * U[d]*U[d];
+      }
+  }
+};
Index: molecuilder/src/atom_trajectoryparticle.hpp
===================================================================
--- molecuilder/src/atom_trajectoryparticle.hpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
+++ molecuilder/src/atom_trajectoryparticle.hpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -0,0 +1,69 @@
+/*
+ * atom_trajectoryparticle.hpp
+ *
+ *  Created on: Oct 19, 2009
+ *      Author: heber
+ */
+
+#ifndef ATOM_TRAJECTORYPARTICLE_HPP_
+#define ATOM_TRAJECTORYPARTICLE_HPP_
+
+using namespace std;
+
+/*********************************************** includes ***********************************/
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <fstream>
+
+#include <gsl/gsl_randist.h>
+
+#include "atom_atominfo.hpp"
+#include "atom_particleinfo.hpp"
+#include "atom_trajectoryparticleinfo.hpp"
+
+/****************************************** forward declarations *****************************/
+
+class atom;
+class config;
+class ForceMatrix;
+class Vector;
+
+/********************************************** declarations *******************************/
+
+
+class TrajectoryParticle :  public TrajectoryParticleInfo, public virtual AtomInfo, public virtual ParticleInfo {
+public:
+
+  TrajectoryParticle();
+  virtual ~TrajectoryParticle();
+
+  // constraint potential and dynamics stuff
+  void AddKineticToTemperature(double *temperature, int step) const;
+  void EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
+  void CorrectVelocity(double *ActualTemp, int Step, Vector *CoGVelocity);
+
+  // trajectory stuff
+  void ResizeTrajectory(int MaxSteps);
+  void CopyStepOnStep(int dest, int src);
+  void VelocityVerletUpdate(int MDSteps, config *configuration, ForceMatrix *Force);
+  void SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity );
+
+  // thermostats
+  void Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin);
+  void Thermostat_Gaussian_init(int Step, double *G, double *E);
+  void Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration);
+  void Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration);
+  void Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration);
+  void Thermostat_NoseHoover_init(int Step, double *delta_alpha);
+  void Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration);
+
+private:
+
+};
+
+
+#endif /* ATOM_TRAJECTORYPARTICLE_HPP_ */
Index: molecuilder/src/atom_trajectoryparticleinfo.cpp
===================================================================
--- molecuilder/src/atom_trajectoryparticleinfo.cpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
+++ molecuilder/src/atom_trajectoryparticleinfo.cpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -0,0 +1,21 @@
+/*
+ * atom_trajectoryparticleinfo.cpp
+ *
+ *  Created on: Oct 19, 2009
+ *      Author: heber
+ */
+
+#include "atom_trajectoryparticleinfo.hpp"
+
+/** Constructor of class TrajectoryParticleInfo.
+ */
+TrajectoryParticleInfo::TrajectoryParticleInfo() : FixedIon(0) {};
+
+/** Destructor of class TrajectoryParticleInfo.
+ */
+TrajectoryParticleInfo::~TrajectoryParticleInfo()
+{
+  Trajectory.R.clear();
+  Trajectory.U.clear();
+  Trajectory.F.clear();
+};
Index: molecuilder/src/atom_trajectoryparticleinfo.hpp
===================================================================
--- molecuilder/src/atom_trajectoryparticleinfo.hpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
+++ molecuilder/src/atom_trajectoryparticleinfo.hpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -0,0 +1,42 @@
+/*
+ * atom_trajectoryparticleinfo.hpp
+ *
+ *  Created on: Oct 19, 2009
+ *      Author: heber
+ */
+
+#ifndef ATOM_TRAJECTORYPARTICLEINFO_HPP_
+#define ATOM_TRAJECTORYPARTICLEINFO_HPP_
+
+using namespace std;
+
+/*********************************************** includes ***********************************/
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <vector>
+
+#include "vector.hpp"
+
+/********************************************** declarations *******************************/
+
+class TrajectoryParticleInfo {
+public:
+  struct
+  {
+    vector<Vector> R;  //!< position vector
+    vector<Vector> U;  //!< velocity vector
+    vector<Vector> F;  //!< last force vector
+  } Trajectory;
+  int FixedIon;   //!< config variable that states whether forces act on the ion or not
+
+  TrajectoryParticleInfo();
+  ~TrajectoryParticleInfo();
+
+private:
+};
+
+#endif /* ATOM_TRAJECTORYPARTICLEINFO_HPP_ */
Index: molecuilder/src/tesselation.cpp
===================================================================
--- molecuilder/src/tesselation.cpp	(revision 3b0ba203af31dbee5c5a38092340126009136ed1)
+++ molecuilder/src/tesselation.cpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -16,17 +16,4 @@
 
 class molecule;
-
-// =========================================== ParticleInfo ====================================
-
-/** Constructor of ParticleInfo.
- */
-ParticleInfo::ParticleInfo() : nr(-1), Name(NULL) {};
-
-/** Destructor of ParticleInfo.
- */
-ParticleInfo::~ParticleInfo()
-{
-  Free(&Name);
-};
 
 // ======================================== Points on Boundary =================================
Index: molecuilder/src/tesselation.hpp
===================================================================
--- molecuilder/src/tesselation.hpp	(revision 3b0ba203af31dbee5c5a38092340126009136ed1)
+++ molecuilder/src/tesselation.hpp	(revision 2e2a70645a2820252e3e6c47dc1d3c15afb2cc94)
@@ -24,6 +24,7 @@
 #include <set>
 
+#include "atom_particleinfo.hpp"
+#include "helpers.hpp"
 #include "vector.hpp"
-#include "helpers.hpp"
 
 /****************************************** forward declarations *****************************/
@@ -144,19 +145,4 @@
 
 // =========================================================== class TESSELPOINT ===========================================
-
-class ParticleInfo {
-public:
-  int nr;       // index to easierly identify
-  char *Name;   // some name to reference to on output
-
-  ParticleInfo();
-  ~ParticleInfo();
-
-  ostream & operator << (ostream &ost);
-
-private:
-};
-
-ostream & operator << (ostream &ost, const ParticleInfo &a);
 
 /** Is a single point of the set of Vectors, also a super-class to be inherited and and its functions to be implemented.
