source: src/atom.cpp@ 4455f4

/** \file atom.cpp
 *
 * Function implementations for the class atom.
 *
 */

#include "atom.hpp"
#include "bond.hpp"
#include "config.hpp"
#include "element.hpp"
#include "lists.hpp"
#include "memoryallocator.hpp"
#include "parser.hpp"
#include "vector.hpp"

/************************************* 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.
 */
atom::atom()
{
  father = this;  // generally, father is itself
  previous = NULL;
  next = NULL;
  sort = &nr;

  // set LCNode::Vector to our Vector
  node = &x;
};

/** Constructor of class atom.
 */
atom::atom(atom *pointer)
{
  previous = NULL;
  next = NULL;
  father = pointer;  // generally, father is itself

  type = pointer->type;  // copy element of atom
  x.CopyVector(&pointer->x); // copy coordination
  v.CopyVector(&pointer->v); // copy velocity
  FixedIon = pointer->FixedIon;
  sort = &nr;
  node = &x;
}


/** Destructor of class atom.
 */
atom::~atom()
{
  unlink(this);
};


/** Climbs up the father list until NULL, last is returned.
 * \return true father, i.e. whose father points to itself, NULL if it could not be found or has none (added hydrogen)
 */
atom *atom::GetTrueFather()
{
  atom *walker = this;
  do {
    if (walker == walker->father) // top most father is the one that points on itself
      break;
    walker = walker->father;
  } while (walker != NULL);
  return walker;
};

/** Sets father to itself or its father in case of copying a molecule.
 */
void atom::CorrectFather()
{
  if (father->father == father)   // same atom in copy's father points to itself
    father = this;  // set father to itself (copy of a whole molecule)
  else
   father = father->father;  // set father to original's father

};

/** Check whether father is equal to given atom.
 * \param *ptr atom to compare father to
 * \param **res return value (only set if atom::father is equal to \a *ptr)
 */
void atom::EqualsFather ( atom *ptr, atom **res )
{
  if ( ptr == father )
    *res = this;
};

/** Checks whether atom is within the given box.
 * \param offset offset to box origin
 * \param *parallelepiped box matrix
 * \return true - is inside, false - is not
 */
bool atom::IsInParallelepiped(Vector offset, double *parallelepiped)
{
  return (node->IsInParallelepiped(offset, parallelepiped));
};

/** Counts the number of bonds weighted by bond::BondDegree.
 * \param bonds times bond::BondDegree
 */
int BondedParticle::CountBonds() const
{
  int NoBonds = 0;
  for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner))
    NoBonds += (*Runner)->BondDegree;
  return NoBonds;
};

/** Output of a single atom.
 * \param ElementNo cardinal number of the element
 * \param AtomNo cardinal number among these atoms of the same element
 * \param *out stream to output to
 * \param *comment commentary after '#' sign
  * \return true - \a *out present, false - \a *out is NULL
 */
bool atom::Output(ofstream *out, int ElementNo, int AtomNo, const char *comment) const
{
  if (out != NULL) {
    *out << "Ion_Type" << ElementNo << "_" << AtomNo << "\t"  << fixed << setprecision(9) << showpoint;
    *out << x.x[0] << "\t" << x.x[1] << "\t" << x.x[2];
    *out << "\t" << FixedIon;
    if (v.Norm() > MYEPSILON)
      *out << "\t" << scientific << setprecision(6) << v.x[0] << "\t" << v.x[1] << "\t" << v.x[2] << "\t";
    if (comment != NULL)
      *out << " # " << comment << endl;
    else
      *out << " # molecule nr " << nr << endl;
    return true;
  } else
    return false;
};
bool atom::Output(ofstream *out, int *ElementNo, int *AtomNo, const char *comment)
{
  AtomNo[type->Z]++;  // increment number
  if (out != NULL) {
    *out << "Ion_Type" << ElementNo[type->Z] << "_" << AtomNo[type->Z] << "\t"  << fixed << setprecision(9) << showpoint;
    *out << x.x[0] << "\t" << x.x[1] << "\t" << x.x[2];
    *out << "\t" << FixedIon;
    if (v.Norm() > MYEPSILON)
      *out << "\t" << scientific << setprecision(6) << v.x[0] << "\t" << v.x[1] << "\t" << v.x[2] << "\t";
    if (comment != NULL)
      *out << " # " << comment << endl;
    else
      *out << " # molecule nr " << nr << endl;
    return true;
  } else
    return false;
};

/** Output of a single atom as one lin in xyz file.
 * \param *out stream to output to
  * \return true - \a *out present, false - \a *out is NULL
 */
bool atom::OutputXYZLine(ofstream *out) const
{
  if (out != NULL) {
    *out << type->symbol << "\t" << x.x[0] << "\t" << x.x[1] << "\t" << x.x[2] << "\t" << endl;
    return true;
  } else
    return false;
};

/** Output of a single atom as one lin 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
 * \param *AtomNo array with atom number in the config file this atom shall have, is increase by one automatically
 * \param step Trajectory time step to output
  * \return true - \a *out present, false - \a *out is NULL
 */
bool atom::OutputTrajectory(ofstream *out, int *ElementNo, int *AtomNo, int step) const
{
  AtomNo[type->Z]++;
  if (out != NULL) {
    *out << "Ion_Type" << ElementNo[type->Z] << "_" << AtomNo[type->Z] << "\t"  << fixed << setprecision(9) << showpoint;
    *out << Trajectory.R.at(step).x[0] << "\t" << Trajectory.R.at(step).x[1] << "\t" << Trajectory.R.at(step).x[2];
    *out << "\t" << FixedIon;
    if (Trajectory.U.at(step).Norm() > MYEPSILON)
      *out << "\t" << scientific << setprecision(6) << Trajectory.U.at(step).x[0] << "\t" << Trajectory.U.at(step).x[1] << "\t" << Trajectory.U.at(step).x[2] << "\t";
    if (Trajectory.F.at(step).Norm() > MYEPSILON)
      *out << "\t" << scientific << setprecision(6) << Trajectory.F.at(step).x[0] << "\t" << Trajectory.F.at(step).x[1] << "\t" << Trajectory.F.at(step).x[2] << "\t";
    *out << "\t# Number in molecule " << nr << endl;
    return true;
  } else
    return false;
};

/** Output of a single atom as one lin in xyz file.
 * \param *out stream to output to
 * \param step Trajectory time step to output
 * \return true - \a *out present, false - \a *out is NULL
 */
bool atom::OutputTrajectoryXYZ(ofstream *out, int step) const
{
  if (out != NULL) {
    *out << type->symbol << "\t";
    *out << Trajectory.R.at(step).x[0] << "\t";
    *out << Trajectory.R.at(step).x[1] << "\t";
    *out << Trajectory.R.at(step).x[2] << endl;
    return true;
  } else
    return false;
};

/** Outputs the MPQC configuration line for this atom.
 * \param *out output stream
 * \param *center center of molecule subtracted from position
 * \param *AtomNo pointer to atom counter that is increased by one
 */
void atom::OutputMPQCLine(ofstream *out, Vector *center, int *AtomNo = NULL) const
{
  *out << "\t\t" << type->symbol << " [ " << x.x[0]-center->x[0] << "\t" << x.x[1]-center->x[1] << "\t" << x.x[2]-center->x[2] << " ]" << endl;
  if (AtomNo != NULL)
    *AtomNo++;
};

ostream & operator << (ostream &ost, const ParticleInfo &a)
{
  ost << "[" << a.Name << "|" << &a << "]";
  return ost;
};

ostream & ParticleInfo::operator << (ostream &ost)
{
  ost << "[" << Name << "|" << this << "]";
  return ost;
};

/** Compares the indices of \a this atom with a given \a ptr.
 * \param ptr atom to compare index against
 * \return true - this one's is smaller, false - not
 */
bool atom::Compare(const atom &ptr)
{
  if (nr < ptr.nr)
    return true;
  else
    return false;
};

/** Returns squared distance to a given vector.
 * \param origin vector to calculate distance to
 * \return distance squared
 */
double atom::DistanceSquaredToVector(Vector &origin)
{
  return origin.DistanceSquared(&x);
};

/** Returns distance to a given vector.
 * \param origin vector to calculate distance to
 * \return distance
 */
double atom::DistanceToVector(Vector &origin)
{
  return origin.Distance(&x);
};

/** Initialises the component number array.
 * Size is set to atom::ListOfBonds.size()+1 (last is th encode end by -1)
 */
void atom::InitComponentNr()
{
  if (ComponentNr != NULL)
    Free(&ComponentNr);
  ComponentNr = Malloc<int>(ListOfBonds.size()+1, "atom::InitComponentNumbers: *ComponentNr");
  for (int i=ListOfBonds.size()+1;i--;)
    ComponentNr[i] = -1;
};


bool operator < (atom &a, atom &b)
{
  return a.Compare(b);
};

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