source: src/molecule.cpp@ ebcade

/** \file molecules.cpp
 *
 * Functions for the class molecule.
 *
 */

#include "config.hpp"
#include "memoryallocator.hpp"
#include "molecule.hpp"

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

/** Constructor of class molecule.
 * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
 */
molecule::molecule(periodentafel *teil)
{
  // init atom chain list
  start = new atom;
  end = new atom;
  start->father = NULL;
  end->father = NULL;
  link(start,end);
  InternalPointer = start;
  // init bond chain list
  first = new bond(start, end, 1, -1);
  last = new bond(start, end, 1, -1);
  link(first,last);
  // other stuff
  MDSteps = 0;
  last_atom = 0;
  elemente = teil;
  AtomCount = 0;
  BondCount = 0;
  NoNonBonds = 0;
  NoNonHydrogen = 0;
  NoCyclicBonds = 0;
  ListOfBondsPerAtom = NULL;
  NumberOfBondsPerAtom = NULL;
  ElementCount = 0;
  for(int i=MAX_ELEMENTS;i--;)
    ElementsInMolecule[i] = 0;
  cell_size[0] = cell_size[2] = cell_size[5]= 20.;
  cell_size[1] = cell_size[3] = cell_size[4]= 0.;
  strcpy(name,"none");
  IndexNr  = -1;
  ActiveFlag = false;
  TesselStruct = NULL;
};

/** Destructor of class molecule.
 * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
 */
molecule::~molecule()
{
  if (ListOfBondsPerAtom != NULL)
    for(int i=AtomCount;i--;)
      Free(&ListOfBondsPerAtom[i]);
  Free(&ListOfBondsPerAtom);
  Free(&NumberOfBondsPerAtom);
  if (TesselStruct != NULL)
    delete(TesselStruct);
  CleanupMolecule();
  delete(first);
  delete(last);
  delete(end);
  delete(start);
};


/** Adds given atom \a *pointer from molecule list.
 * Increases molecule::last_atom and gives last number to added atom and names it according to its element::abbrev and molecule::AtomCount
 * \param *pointer allocated and set atom
 * \return true - succeeded, false - atom not found in list
 */
bool molecule::AddAtom(atom *pointer)
{
  if (pointer != NULL) {
    pointer->sort = &pointer->nr;
    pointer->nr = last_atom++;  // increase number within molecule
    AtomCount++;
    if (pointer->type != NULL) {
      if (ElementsInMolecule[pointer->type->Z] == 0)
        ElementCount++;
      ElementsInMolecule[pointer->type->Z]++; // increase number of elements
      if (pointer->type->Z != 1)
        NoNonHydrogen++;
      if (pointer->Name == NULL) {
        Free(&pointer->Name);
        pointer->Name = Malloc<char>(6, "molecule::AddAtom: *pointer->Name");
        sprintf(pointer->Name, "%2s%02d", pointer->type->symbol, pointer->nr+1);
      }
    }
    return add(pointer, end);
  } else
    return false;
};

/** Adds a copy of the given atom \a *pointer from molecule list.
 * Increases molecule::last_atom and gives last number to added atom.
 * \param *pointer allocated and set atom
 * \return pointer to the newly added atom
 */
atom * molecule::AddCopyAtom(atom *pointer)
{
  if (pointer != NULL) {
    atom *walker = new atom(pointer);
    walker->Name = Malloc<char>(strlen(pointer->Name) + 1, "atom::atom: *Name");
    strcpy (walker->Name, pointer->Name);
    walker->nr = last_atom++;  // increase number within molecule
    add(walker, end);
    if ((pointer->type != NULL) && (pointer->type->Z != 1))
      NoNonHydrogen++;
    AtomCount++;
    return walker;
  } else
    return NULL;
};

/** Adds a Hydrogen atom in replacement for the given atom \a *partner in bond with a *origin.
 * Here, we have to distinguish between single, double or triple bonds as stated by \a BondDegree, that each demand
 * a different scheme when adding \a *replacement atom for the given one.
 * -# Single Bond: Simply add new atom with bond distance rescaled to typical hydrogen one
 * -# Double Bond: Here, we need the **BondList of the \a *origin atom, by scanning for the other bonds instead of
 *    *Bond, we use the through these connected atoms to determine the plane they lie in, vector::MakeNormalvector().
 *    The orthonormal vector to this plane along with the vector in *Bond direction determines the plane the two
 *    replacing hydrogens shall lie in. Now, all remains to do is take the usual hydrogen double bond angle for the
 *    element of *origin and form the sin/cos admixture of both plane vectors for the new coordinates of the two
 *    hydrogens forming this angle with *origin.
 * -# Triple Bond: The idea is to set up a tetraoid (C1-H1-H2-H3) (however the lengths \f$b\f$ of the sides of the base
 *    triangle formed by the to be added hydrogens are not equal to the typical bond distance \f$l\f$ but have to be
 *    determined from the typical angle \f$\alpha\f$ for a hydrogen triple connected to the element of *origin):
 *    We have the height \f$d\f$ as the vector in *Bond direction (from triangle C1-H1-H2).
 *    \f[ h = l \cdot \cos{\left (\frac{\alpha}{2} \right )} \qquad b = 2l \cdot \sin{\left (\frac{\alpha}{2} \right)} \quad \rightarrow \quad d = l \cdot \sqrt{\cos^2{\left (\frac{\alpha}{2} \right)}-\frac{1}{3}\cdot\sin^2{\left (\frac{\alpha}{2}\right )}}
 *    \f]
 *    vector::GetNormalvector() creates one orthonormal vector from this *Bond vector and vector::MakeNormalvector creates
 *    the third one from the former two vectors. The latter ones form the plane of the base triangle mentioned above.
 *    The lengths for these are \f$f\f$ and \f$g\f$ (from triangle H1-H2-(center of H1-H2-H3)) with knowledge that
 *    the median lines in an isosceles triangle meet in the center point with a ratio 2:1.
 *    \f[ f = \frac{b}{\sqrt{3}} \qquad g = \frac{b}{2}
 *    \f]
 *    as the coordination of all three atoms in the coordinate system of these three vectors:
 *    \f$\pmatrix{d & f & 0}\f$, \f$\pmatrix{d & -0.5 \cdot f & g}\f$ and \f$\pmatrix{d & -0.5 \cdot f & -g}\f$.
 *
 * \param *out output stream for debugging
 * \param *Bond pointer to bond between \a *origin and \a *replacement
 * \param *TopOrigin son of \a *origin of upper level molecule (the atom added to this molecule as a copy of \a *origin)
 * \param *origin pointer to atom which acts as the origin for scaling the added hydrogen to correct bond length
 * \param *replacement pointer to the atom which shall be copied as a hydrogen atom in this molecule
 * \param **BondList list of bonds \a *replacement has (necessary to determine plane for double and triple bonds)
 * \param NumBond  number of bonds in \a **BondList
 * \param isAngstroem whether the coordination of the given atoms is in AtomicLength (false) or Angstrom(true)
 * \return number of atoms added, if < bond::BondDegree then something went wrong
 * \todo double and triple bonds splitting (always use the tetraeder angle!)
 */
bool molecule::AddHydrogenReplacementAtom(ofstream *out, bond *TopBond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bond **BondList, int NumBond, bool IsAngstroem)
{
  double bondlength;  // bond length of the bond to be replaced/cut
  double bondangle;  // bond angle of the bond to be replaced/cut
  double BondRescale;   // rescale value for the hydrogen bond length
  bool AllWentWell = true;    // flag gathering the boolean return value of molecule::AddAtom and other functions, as return value on exit
  bond *FirstBond = NULL, *SecondBond = NULL; // Other bonds in double bond case to determine "other" plane
  atom *FirstOtherAtom = NULL, *SecondOtherAtom = NULL, *ThirdOtherAtom = NULL; // pointer to hydrogen atoms to be added
  double b,l,d,f,g, alpha, factors[NDIM];    // hold temporary values in triple bond case for coordination determination
  Vector Orthovector1, Orthovector2;  // temporary vectors in coordination construction
  Vector InBondvector;    // vector in direction of *Bond
  bond *Binder = NULL;
  double *matrix;

//  *out << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
  // create vector in direction of bond
  InBondvector.CopyVector(&TopReplacement->x);
  InBondvector.SubtractVector(&TopOrigin->x);
  bondlength = InBondvector.Norm();

   // is greater than typical bond distance? Then we have to correct periodically
   // the problem is not the H being out of the box, but InBondvector have the wrong direction
   // due to TopReplacement or Origin being on the wrong side!
  if (bondlength > BondDistance) {
//    *out << Verbose(4) << "InBondvector is: ";
//    InBondvector.Output(out);
//    *out << endl;
    Orthovector1.Zero();
    for (int i=NDIM;i--;) {
      l = TopReplacement->x.x[i] - TopOrigin->x.x[i];
      if (fabs(l) > BondDistance) { // is component greater than bond distance
        Orthovector1.x[i] = (l < 0) ? -1. : +1.;
      } // (signs are correct, was tested!)
    }
    matrix = ReturnFullMatrixforSymmetric(cell_size);
    Orthovector1.MatrixMultiplication(matrix);
    InBondvector.SubtractVector(&Orthovector1); // subtract just the additional translation
    Free(&matrix);
    bondlength = InBondvector.Norm();
//    *out << Verbose(4) << "Corrected InBondvector is now: ";
//    InBondvector.Output(out);
//    *out << endl;
  } // periodic correction finished

  InBondvector.Normalize();
  // get typical bond length and store as scale factor for later
  BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1];
  if (BondRescale == -1) {
    cerr << Verbose(3) << "ERROR: There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
    return false;
    BondRescale = bondlength;
  } else {
    if (!IsAngstroem)
      BondRescale /= (1.*AtomicLengthToAngstroem);
  }

  // discern single, double and triple bonds
  switch(TopBond->BondDegree) {
    case 1:
      FirstOtherAtom = new atom();    // new atom
      FirstOtherAtom->type = elemente->FindElement(1);  // element is Hydrogen
      FirstOtherAtom->v.CopyVector(&TopReplacement->v); // copy velocity
      FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
      if (TopReplacement->type->Z == 1) { // neither rescale nor replace if it's already hydrogen
        FirstOtherAtom->father = TopReplacement;
        BondRescale = bondlength;
      } else {
        FirstOtherAtom->father = NULL;  // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
      }
      InBondvector.Scale(&BondRescale);   // rescale the distance vector to Hydrogen bond length
      FirstOtherAtom->x.CopyVector(&TopOrigin->x); // set coordination to origin ...
      FirstOtherAtom->x.AddVector(&InBondvector);  // ... and add distance vector to replacement atom
      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
//      *out << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
//      FirstOtherAtom->x.Output(out);
//      *out << endl;
      Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
      Binder->Cyclic = false;
      Binder->Type = TreeEdge;
      break;
    case 2:
      // determine two other bonds (warning if there are more than two other) plus valence sanity check
      for (int i=0;i<NumBond;i++) {
        if (BondList[i] != TopBond) {
          if (FirstBond == NULL) {
            FirstBond = BondList[i];
            FirstOtherAtom = BondList[i]->GetOtherAtom(TopOrigin);
          } else if (SecondBond == NULL) {
            SecondBond = BondList[i];
            SecondOtherAtom = BondList[i]->GetOtherAtom(TopOrigin);
          } else {
            *out << Verbose(3) << "WARNING: Detected more than four bonds for atom " << TopOrigin->Name;
          }
        }
      }
      if (SecondOtherAtom == NULL) {  // then we have an atom with valence four, but only 3 bonds: one to replace and one which is TopBond (third is FirstBond)
        SecondBond = TopBond;
        SecondOtherAtom = TopReplacement;
      }
      if (FirstOtherAtom != NULL) { // then we just have this double bond and the plane does not matter at all
//        *out << Verbose(3) << "Regarding the double bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") to be constructed: Taking " << FirstOtherAtom->Name << " and " << SecondOtherAtom->Name << " along with " << TopOrigin->Name << " to determine orthogonal plane." << endl;

        // determine the plane of these two with the *origin
        AllWentWell = AllWentWell && Orthovector1.MakeNormalVector(&TopOrigin->x, &FirstOtherAtom->x, &SecondOtherAtom->x);
      } else {
        Orthovector1.GetOneNormalVector(&InBondvector);
      }
      //*out << Verbose(3)<< "Orthovector1: ";
      //Orthovector1.Output(out);
      //*out << endl;
      // orthogonal vector and bond vector between origin and replacement form the new plane
      Orthovector1.MakeNormalVector(&InBondvector);
      Orthovector1.Normalize();
      //*out << Verbose(3) << "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << "." << endl;

      // create the two Hydrogens ...
      FirstOtherAtom = new atom();
      SecondOtherAtom = new atom();
      FirstOtherAtom->type = elemente->FindElement(1);
      SecondOtherAtom->type = elemente->FindElement(1);
      FirstOtherAtom->v.CopyVector(&TopReplacement->v); // copy velocity
      FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
      SecondOtherAtom->v.CopyVector(&TopReplacement->v); // copy velocity
      SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
      FirstOtherAtom->father = NULL;  // we are just an added hydrogen with no father
      SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
      bondangle = TopOrigin->type->HBondAngle[1];
      if (bondangle == -1) {
        *out << Verbose(3) << "ERROR: There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
        return false;
        bondangle = 0;
      }
      bondangle *= M_PI/180./2.;
//      *out << Verbose(3) << "ReScaleCheck: InBondvector ";
//      InBondvector.Output(out);
//      *out << endl;
//      *out << Verbose(3) << "ReScaleCheck: Orthovector ";
//      Orthovector1.Output(out);
//      *out << endl;
//      *out << Verbose(3) << "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle) << endl;
      FirstOtherAtom->x.Zero();
      SecondOtherAtom->x.Zero();
      for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction)
        FirstOtherAtom->x.x[i] = InBondvector.x[i] * cos(bondangle) + Orthovector1.x[i] * (sin(bondangle));
        SecondOtherAtom->x.x[i] = InBondvector.x[i] * cos(bondangle) + Orthovector1.x[i] * (-sin(bondangle));
      }
      FirstOtherAtom->x.Scale(&BondRescale);  // rescale by correct BondDistance
      SecondOtherAtom->x.Scale(&BondRescale);
      //*out << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl;
      for(int i=NDIM;i--;) { // and make relative to origin atom
        FirstOtherAtom->x.x[i] += TopOrigin->x.x[i];
        SecondOtherAtom->x.x[i] += TopOrigin->x.x[i];
      }
      // ... and add to molecule
      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
      AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
//      *out << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
//      FirstOtherAtom->x.Output(out);
//      *out << endl;
//      *out << Verbose(4) << "Added " << *SecondOtherAtom << " at: ";
//      SecondOtherAtom->x.Output(out);
//      *out << endl;
      Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
      Binder->Cyclic = false;
      Binder->Type = TreeEdge;
      Binder = AddBond(BottomOrigin, SecondOtherAtom, 1);
      Binder->Cyclic = false;
      Binder->Type = TreeEdge;
      break;
    case 3:
      // take the "usual" tetraoidal angle and add the three Hydrogen in direction of the bond (height of the tetraoid)
      FirstOtherAtom = new atom();
      SecondOtherAtom = new atom();
      ThirdOtherAtom = new atom();
      FirstOtherAtom->type = elemente->FindElement(1);
      SecondOtherAtom->type = elemente->FindElement(1);
      ThirdOtherAtom->type = elemente->FindElement(1);
      FirstOtherAtom->v.CopyVector(&TopReplacement->v); // copy velocity
      FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
      SecondOtherAtom->v.CopyVector(&TopReplacement->v); // copy velocity
      SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
      ThirdOtherAtom->v.CopyVector(&TopReplacement->v); // copy velocity
      ThirdOtherAtom->FixedIon = TopReplacement->FixedIon;
      FirstOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
      SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
      ThirdOtherAtom->father = NULL;  //  we are just an added hydrogen with no father

      // we need to vectors orthonormal the InBondvector
      AllWentWell = AllWentWell && Orthovector1.GetOneNormalVector(&InBondvector);
//      *out << Verbose(3) << "Orthovector1: ";
//      Orthovector1.Output(out);
//      *out << endl;
      AllWentWell = AllWentWell && Orthovector2.MakeNormalVector(&InBondvector, &Orthovector1);
//      *out << Verbose(3) << "Orthovector2: ";
//      Orthovector2.Output(out);
//      *out << endl;

      // create correct coordination for the three atoms
      alpha = (TopOrigin->type->HBondAngle[2])/180.*M_PI/2.;  // retrieve triple bond angle from database
      l = BondRescale;        // desired bond length
      b = 2.*l*sin(alpha);    // base length of isosceles triangle
      d = l*sqrt(cos(alpha)*cos(alpha) - sin(alpha)*sin(alpha)/3.);   // length for InBondvector
      f = b/sqrt(3.);   // length for Orthvector1
      g = b/2.;         // length for Orthvector2
//      *out << Verbose(3) << "Bond length and half-angle: " << l << ", " << alpha << "\t (b,d,f,g) = " << b << ", " << d << ", " << f << ", " << g << ", " << endl;
//      *out << Verbose(3) << "The three Bond lengths: " << sqrt(d*d+f*f) << ", " << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g) << ", "  << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g) << endl;
      factors[0] = d;
      factors[1] = f;
      factors[2] = 0.;
      FirstOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
      factors[1] = -0.5*f;
      factors[2] = g;
      SecondOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
      factors[2] = -g;
      ThirdOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);

      // rescale each to correct BondDistance
//      FirstOtherAtom->x.Scale(&BondRescale);
//      SecondOtherAtom->x.Scale(&BondRescale);
//      ThirdOtherAtom->x.Scale(&BondRescale);

      // and relative to *origin atom
      FirstOtherAtom->x.AddVector(&TopOrigin->x);
      SecondOtherAtom->x.AddVector(&TopOrigin->x);
      ThirdOtherAtom->x.AddVector(&TopOrigin->x);

      // ... and add to molecule
      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
      AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
      AllWentWell = AllWentWell && AddAtom(ThirdOtherAtom);
//      *out << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
//      FirstOtherAtom->x.Output(out);
//      *out << endl;
//      *out << Verbose(4) << "Added " << *SecondOtherAtom << " at: ";
//      SecondOtherAtom->x.Output(out);
//      *out << endl;
//      *out << Verbose(4) << "Added " << *ThirdOtherAtom << " at: ";
//      ThirdOtherAtom->x.Output(out);
//      *out << endl;
      Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
      Binder->Cyclic = false;
      Binder->Type = TreeEdge;
      Binder = AddBond(BottomOrigin, SecondOtherAtom, 1);
      Binder->Cyclic = false;
      Binder->Type = TreeEdge;
      Binder = AddBond(BottomOrigin, ThirdOtherAtom, 1);
      Binder->Cyclic = false;
      Binder->Type = TreeEdge;
      break;
    default:
      cerr << "ERROR: BondDegree does not state single, double or triple bond!" << endl;
      AllWentWell = false;
      break;
  }

//  *out << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;
  return AllWentWell;
};

/** Adds given atom \a *pointer from molecule list.
 * Increases molecule::last_atom and gives last number to added atom.
 * \param filename name and path of xyz file
 * \return true - succeeded, false - file not found
 */
bool molecule::AddXYZFile(string filename)
{
  istringstream *input = NULL;
  int NumberOfAtoms = 0; // atom number in xyz read
  int i, j; // loop variables
  atom *Walker = NULL;  // pointer to added atom
  char shorthand[3];  // shorthand for atom name
  ifstream xyzfile;   // xyz file
  string line;    // currently parsed line
  double x[3];    // atom coordinates

  xyzfile.open(filename.c_str());
  if (!xyzfile)
    return false;

  getline(xyzfile,line,'\n'); // Read numer of atoms in file
  input = new istringstream(line);
  *input >> NumberOfAtoms;
  cout << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl;
  getline(xyzfile,line,'\n'); // Read comment
  cout << Verbose(1) << "Comment: " << line << endl;

  if (MDSteps == 0) // no atoms yet present
    MDSteps++;
  for(i=0;i<NumberOfAtoms;i++){
    Walker = new atom;
    getline(xyzfile,line,'\n');
    istringstream *item = new istringstream(line);
    //istringstream input(line);
    //cout << Verbose(1) << "Reading: " << line << endl;
    *item >> shorthand;
    *item >> x[0];
    *item >> x[1];
    *item >> x[2];
    Walker->type = elemente->FindElement(shorthand);
    if (Walker->type == NULL) {
      cerr << "Could not parse the element at line: '" << line << "', setting to H.";
      Walker->type = elemente->FindElement(1);
    }
    if (Trajectories[Walker].R.size() <= (unsigned int)MDSteps) {
      Trajectories[Walker].R.resize(MDSteps+10);
      Trajectories[Walker].U.resize(MDSteps+10);
      Trajectories[Walker].F.resize(MDSteps+10);
    }
    for(j=NDIM;j--;) {
      Walker->x.x[j] = x[j];
      Trajectories[Walker].R.at(MDSteps-1).x[j] = x[j];
      Trajectories[Walker].U.at(MDSteps-1).x[j] = 0;
      Trajectories[Walker].F.at(MDSteps-1).x[j] = 0;
    }
    AddAtom(Walker);  // add to molecule
    delete(item);
  }
  xyzfile.close();
  delete(input);
  return true;
};

/** Creates a copy of this molecule.
 * \return copy of molecule
 */
molecule *molecule::CopyMolecule()
{
  molecule *copy = new molecule(elemente);
  atom *CurrentAtom = NULL;
  atom *LeftAtom = NULL, *RightAtom = NULL;
  atom *Walker = NULL;

  // copy all atoms
  Walker = start;
  while(Walker->next != end) {
    Walker = Walker->next;
    CurrentAtom = copy->AddCopyAtom(Walker);
  }

  // copy all bonds
  bond *Binder = first;
  bond *NewBond = NULL;
  while(Binder->next != last) {
    Binder = Binder->next;

    // get the pendant atoms of current bond in the copy molecule
    copy->ActOnAllAtoms( &atom::EqualsFather, Binder->leftatom, &LeftAtom );
    copy->ActOnAllAtoms( &atom::EqualsFather, Binder->rightatom, &RightAtom );

    NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->BondDegree);
    NewBond->Cyclic = Binder->Cyclic;
    if (Binder->Cyclic)
      copy->NoCyclicBonds++;
    NewBond->Type = Binder->Type;
  }
  // correct fathers
  ActOnAllAtoms( &atom::CorrectFather );

  // copy values
  copy->CountAtoms((ofstream *)&cout);
  copy->CountElements();
  if (first->next != last) {  // if adjaceny list is present
    copy->BondDistance = BondDistance;
    copy->CreateListOfBondsPerAtom((ofstream *)&cout);
  }

  return copy;
};


/**
 * Copies all atoms of a molecule which are within the defined parallelepiped.
 *
 * @param offest for the origin of the parallelepiped
 * @param three vectors forming the matrix that defines the shape of the parallelpiped
 */
molecule* molecule::CopyMoleculeFromSubRegion(Vector offset, double *parallelepiped) {
  molecule *copy = new molecule(elemente);
  atom *Walker = start;

  while(Walker->next != end) {
    Walker = Walker->next;
    if (Walker->x.IsInParallelepiped(offset, parallelepiped)) {
      cout << "Adding atom " << *Walker << endl;
      copy->AddCopyAtom(Walker);
    }
  }

  //TODO: copy->BuildInducedSubgraph((ofstream *)&cout, this);

  return copy;
}

/** Adds a bond to a the molecule specified by two atoms, \a *first and \a *second.
 * Also updates molecule::BondCount and molecule::NoNonBonds.
 * \param *first first atom in bond
 * \param *second atom in bond
 * \return pointer to bond or NULL on failure
 */
bond * molecule::AddBond(atom *atom1, atom *atom2, int degree)
{
  bond *Binder = NULL;
  if ((atom1 != NULL) && (FindAtom(atom1->nr) != NULL) && (atom2 != NULL) && (FindAtom(atom2->nr) != NULL)) {
    Binder = new bond(atom1, atom2, degree, BondCount++);
    if ((atom1->type != NULL) && (atom1->type->Z != 1) && (atom2->type != NULL) && (atom2->type->Z != 1))
      NoNonBonds++;
    add(Binder, last);
  } else {
    cerr << Verbose(1) << "ERROR: Could not add bond between " << atom1->Name << " and " << atom2->Name << " as one or both are not present in the molecule." << endl;
  }
  return Binder;
};

/** Remove bond from bond chain list.
 * \todo Function not implemented yet
 * \param *pointer bond pointer
 * \return true - bound found and removed, false - bond not found/removed
 */
bool molecule::RemoveBond(bond *pointer)
{
  //cerr << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
  removewithoutcheck(pointer);
  return true;
};

/** Remove every bond from bond chain list that atom \a *BondPartner is a constituent of.
 * \todo Function not implemented yet
 * \param *BondPartner atom to be removed
 * \return true - bounds found and removed, false - bonds not found/removed
 */
bool molecule::RemoveBonds(atom *BondPartner)
{
  cerr << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
  return false;
};

/** Set molecule::name from the basename without suffix in the given \a *filename.
 * \param *filename filename
 */
void molecule::SetNameFromFilename(const char *filename)
{
  int length = 0;
  const char *molname = strrchr(filename, '/');
  if (molname != NULL)
    molname += sizeof(char);  // search for filename without dirs
  else
    molname = filename; // contains no slashes
  char *endname = strchr(molname, '.');
  if ((endname == NULL) || (endname < molname))
    length = strlen(molname);
  else
    length = strlen(molname) - strlen(endname);
  strncpy(name, molname, length);
  name[length]='\0';
};

/** Sets the molecule::cell_size to the components of \a *dim (rectangular box)
 * \param *dim vector class
 */
void molecule::SetBoxDimension(Vector *dim)
{
  cell_size[0] = dim->x[0];
  cell_size[1] = 0.;
  cell_size[2] = dim->x[1];
  cell_size[3] = 0.;
  cell_size[4] = 0.;
  cell_size[5] = dim->x[2];
};

/** Removes atom from molecule list and deletes it.
 * \param *pointer atom to be removed
 * \return true - succeeded, false - atom not found in list
 */
bool molecule::RemoveAtom(atom *pointer)
{
  if (ElementsInMolecule[pointer->type->Z] != 0)  { // this would indicate an error
    ElementsInMolecule[pointer->type->Z]--;  // decrease number of atom of this element
    AtomCount--;
  } else
    cerr << "ERROR: Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
  if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
    ElementCount--;
  Trajectories.erase(pointer);
  return remove(pointer, start, end);
};

/** Removes atom from molecule list, but does not delete it.
 * \param *pointer atom to be removed
 * \return true - succeeded, false - atom not found in list
 */
bool molecule::UnlinkAtom(atom *pointer)
{
  if (pointer == NULL)
    return false;
  if (ElementsInMolecule[pointer->type->Z] != 0)  // this would indicate an error
    ElementsInMolecule[pointer->type->Z]--; // decrease number of atom of this element
  else
    cerr << "ERROR: Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
  if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
    ElementCount--;
  Trajectories.erase(pointer);
  unlink(pointer);
  return true;
};

/** Removes every atom from molecule list.
 * \return true - succeeded, false - atom not found in list
 */
bool molecule::CleanupMolecule()
{
  return (cleanup(start,end) && cleanup(first,last));
};

/** Finds an atom specified by its continuous number.
 * \param Nr number of atom withim molecule
 * \return pointer to atom or NULL
 */
atom * molecule::FindAtom(int Nr)  const{
  atom * walker = find(&Nr, start,end);
  if (walker != NULL) {
    //cout << Verbose(0) << "Found Atom Nr. " << walker->nr << endl;
    return walker;
  } else {
    cout << Verbose(0) << "Atom not found in list." << endl;
    return NULL;
  }
};

/** Asks for atom number, and checks whether in list.
 * \param *text question before entering
 */
atom * molecule::AskAtom(string text)
{
  int No;
  atom *ion = NULL;
  do {
    //cout << Verbose(0) << "============Atom list==========================" << endl;
    //mol->Output((ofstream *)&cout);
    //cout << Verbose(0) << "===============================================" << endl;
    cout << Verbose(0) << text;
    cin >> No;
    ion = this->FindAtom(No);
  } while (ion == NULL);
  return ion;
};

/** Checks if given coordinates are within cell volume.
 * \param *x array of coordinates
 * \return true - is within, false - out of cell
 */
bool molecule::CheckBounds(const Vector *x) const
{
  bool result = true;
  int j =-1;
  for (int i=0;i<NDIM;i++) {
    j += i+1;
    result = result && ((x->x[i] >= 0) && (x->x[i] < cell_size[j]));
  }
  //return result;
  return true; /// probably not gonna use the check no more
};

/** Prints molecule to *out.
 * \param *out output stream
 */
bool molecule::Output(ofstream *out)
{
  atom *walker = NULL;
  int ElementNo[MAX_ELEMENTS], AtomNo[MAX_ELEMENTS];
  CountElements();

  for (int i=0;i<MAX_ELEMENTS;++i) {
    AtomNo[i] = 0;
    ElementNo[i] = 0;
  }
  if (out == NULL) {
    return false;
  } else {
    *out << "#Ion_TypeNr._Nr.R[0]    R[1]    R[2]    MoveType (0 MoveIon, 1 FixedIon)" << endl;
    walker = start;
    while (walker->next != end) { // go through every atom of this element
      walker = walker->next;
      ElementNo[walker->type->Z] = 1;
    }
    int current=1;
    for (int i=0;i<MAX_ELEMENTS;++i) {
      if (ElementNo[i] == 1)
        ElementNo[i] = current++;
    }
    walker = start;
    while (walker->next != end) { // go through every atom of this element
      walker = walker->next;
      AtomNo[walker->type->Z]++;
      walker->Output(ElementNo[walker->type->Z], AtomNo[walker->type->Z], out); // removed due to trajectories
    }
    return true;
  }
};

/** Prints molecule with all atomic trajectory positions to *out.
 * \param *out output stream
 */
bool molecule::OutputTrajectories(ofstream *out)
{
  atom *walker = NULL;
  int ElementNo[MAX_ELEMENTS], AtomNo[MAX_ELEMENTS];
  CountElements();

  if (out == NULL) {
    return false;
  } else {
    for (int step = 0; step < MDSteps; step++) {
      if (step == 0) {
        *out << "#Ion_TypeNr._Nr.R[0]    R[1]    R[2]    MoveType (0 MoveIon, 1 FixedIon)" << endl;
      } else {
        *out << "# ====== MD step " << step << " =========" << endl;
      }
      for (int i=0;i<MAX_ELEMENTS;++i) {
        AtomNo[i] = 0;
        ElementNo[i] = 0;
      }
    walker = start;
    while (walker->next != end) { // go through every atom of this element
      walker = walker->next;
      ElementNo[walker->type->Z] = 1;
    }
    int current=1;
    for (int i=0;i<MAX_ELEMENTS;++i) {
      if (ElementNo[i] == 1)
        ElementNo[i] = current++;
    }
    walker = start;
      while (walker->next != end) { // go through every atom of this element
        walker = walker->next;
        AtomNo[walker->type->Z]++;
        *out << "Ion_Type" << ElementNo[walker->type->Z] << "_" << AtomNo[walker->type->Z] << "\t"  << fixed << setprecision(9) << showpoint;
        *out << Trajectories[walker].R.at(step).x[0] << "\t" << Trajectories[walker].R.at(step).x[1] << "\t" << Trajectories[walker].R.at(step).x[2];
        *out << "\t" << walker->FixedIon;
        if (Trajectories[walker].U.at(step).Norm() > MYEPSILON)
          *out << "\t" << scientific << setprecision(6) << Trajectories[walker].U.at(step).x[0] << "\t" << Trajectories[walker].U.at(step).x[1] << "\t" << Trajectories[walker].U.at(step).x[2] << "\t";
        if (Trajectories[walker].F.at(step).Norm() > MYEPSILON)
          *out << "\t" << scientific << setprecision(6) << Trajectories[walker].F.at(step).x[0] << "\t" << Trajectories[walker].F.at(step).x[1] << "\t" << Trajectories[walker].F.at(step).x[2] << "\t";
        *out << "\t# Number in molecule " << walker->nr << endl;
      }
    }
    return true;
  }
};

/** Outputs contents of molecule::ListOfBondsPerAtom.
 * \param *out output stream
 */
void molecule::OutputListOfBonds(ofstream *out) const
{
  *out << Verbose(2) << endl << "From Contents of ListOfBondsPerAtom, all non-hydrogen atoms:" << endl;
  atom *Walker = start;
  while (Walker->next != end) {
    Walker = Walker->next;
#ifdef ADDHYDROGEN
    if (Walker->type->Z != 1) {   // regard only non-hydrogen
#endif
      *out << Verbose(2) << "Atom " << Walker->Name << " has Bonds: "<<endl;
      for(int j=0;j<NumberOfBondsPerAtom[Walker->nr];j++) {
        *out << Verbose(3) << *(ListOfBondsPerAtom)[Walker->nr][j] << endl;
      }
#ifdef ADDHYDROGEN
    }
#endif
  }
  *out << endl;
};

/** Output of element before the actual coordination list.
 * \param *out stream pointer
 */
bool molecule::Checkout(ofstream *out)  const
{
  return elemente->Checkout(out, ElementsInMolecule);
};

/** Prints molecule with all its trajectories to *out as xyz file.
 * \param *out output stream
 */
bool molecule::OutputTrajectoriesXYZ(ofstream *out)
{
  atom *walker = NULL;
  int No = 0;
  time_t now;

  now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
  walker = start;
  while (walker->next != end) { // go through every atom and count
    walker = walker->next;
    No++;
  }
  if (out != NULL) {
    for (int step=0;step<MDSteps;step++) {
      *out << No << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
      walker = start;
      while (walker->next != end) { // go through every atom of this element
        walker = walker->next;
        *out << walker->type->symbol << "\t" << Trajectories[walker].R.at(step).x[0] << "\t" << Trajectories[walker].R.at(step).x[1] << "\t" << Trajectories[walker].R.at(step).x[2] << endl;
      }
    }
    return true;
  } else
    return false;
};

/** Prints molecule to *out as xyz file.
* \param *out output stream
 */
bool molecule::OutputXYZ(ofstream *out) const
{
  atom *walker = NULL;
  int AtomNo = 0;
  time_t now;

  now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
  walker = start;
  while (walker->next != end) { // go through every atom and count
    walker = walker->next;
    AtomNo++;
  }
  if (out != NULL) {
    *out << AtomNo << "\n\tCreated by molecuilder on " << ctime(&now);
    walker = start;
    while (walker->next != end) { // go through every atom of this element
      walker = walker->next;
      walker->OutputXYZLine(out);
    }
    return true;
  } else
    return false;
};

/** Brings molecule::AtomCount and atom::*Name up-to-date.
 * \param *out output stream for debugging
 */
void molecule::CountAtoms(ofstream *out)
{
  int i = 0;
  atom *Walker = start;
  while (Walker->next != end) {
    Walker = Walker->next;
    i++;
  }
  if ((AtomCount == 0) || (i != AtomCount)) {
    *out << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl;
    AtomCount = i;

    // count NonHydrogen atoms and give each atom a unique name
    if (AtomCount != 0) {
      i=0;
      NoNonHydrogen = 0;
      Walker = start;
      while (Walker->next != end) {
        Walker = Walker->next;
        Walker->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
        if (Walker->type->Z != 1) // count non-hydrogen atoms whilst at it
          NoNonHydrogen++;
        Free(&Walker->Name);
        Walker->Name = Malloc<char>(6, "molecule::CountAtoms: *walker->Name");
        sprintf(Walker->Name, "%2s%02d", Walker->type->symbol, Walker->nr+1);
        *out << "Naming atom nr. " << Walker->nr << " " << Walker->Name << "." << endl;
        i++;
      }
    } else
      *out << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl;
  }
};

/** Brings molecule::ElementCount and molecule::ElementsInMolecule up-to-date.
 */
void molecule::CountElements()
{
  int i = 0;
  for(i=MAX_ELEMENTS;i--;)
    ElementsInMolecule[i] = 0;
  ElementCount = 0;

  atom *walker = start;
  while (walker->next != end) {
    walker = walker->next;
    ElementsInMolecule[walker->type->Z]++;
    i++;
  }
  for(i=MAX_ELEMENTS;i--;)
    ElementCount += (ElementsInMolecule[i] != 0 ? 1 : 0);
};



/** Counts necessary number of valence electrons and returns number and SpinType.
 * \param configuration containing everything
 */
void molecule::CalculateOrbitals(class config &configuration)
{
  configuration.MaxPsiDouble = configuration.PsiMaxNoDown = configuration.PsiMaxNoUp = configuration.PsiType = 0;
  for(int i=MAX_ELEMENTS;i--;) {
    if (ElementsInMolecule[i] != 0) {
      //cout << "CalculateOrbitals: " << elemente->FindElement(i)->name << " has a valence of " << (int)elemente->FindElement(i)->Valence << " and there are " << ElementsInMolecule[i] << " of it." << endl;
      configuration.MaxPsiDouble += ElementsInMolecule[i]*((int)elemente->FindElement(i)->Valence);
    }
  }
  configuration.PsiMaxNoDown = configuration.MaxPsiDouble/2 + (configuration.MaxPsiDouble % 2);
  configuration.PsiMaxNoUp = configuration.MaxPsiDouble/2;
  configuration.MaxPsiDouble /= 2;
  configuration.PsiType = (configuration.PsiMaxNoDown == configuration.PsiMaxNoUp) ? 0 : 1;
  if ((configuration.PsiType == 1) && (configuration.ProcPEPsi < 2)) {
    configuration.ProcPEGamma /= 2;
    configuration.ProcPEPsi *= 2;
  } else {
    configuration.ProcPEGamma *= configuration.ProcPEPsi;
    configuration.ProcPEPsi = 1;
  }
  configuration.InitMaxMinStopStep = configuration.MaxMinStopStep = configuration.MaxPsiDouble;
};


/** Creates an 2d array of pointer with an entry for each atom and each bond it has.
 * Updates molecule::ListOfBondsPerAtom, molecule::NumberOfBondsPerAtom by parsing through
 * bond chain list, using molecule::AtomCount and molecule::BondCount.
 * Allocates memory, fills the array and exits
 * \param *out output stream for debugging
 */
void molecule::CreateListOfBondsPerAtom(ofstream *out)
{
  bond *Binder = NULL;
  atom *Walker = NULL;
  int TotalDegree;
  *out << Verbose(1) << "Begin of Creating ListOfBondsPerAtom: AtomCount = " << AtomCount << "\tBondCount = " << BondCount << "\tNoNonBonds = " << NoNonBonds << "." << endl;

  // re-allocate memory
  *out << Verbose(2) << "(Re-)Allocating memory." << endl;
  if (ListOfBondsPerAtom != NULL) {
    for(int i=AtomCount;i--;)
      Free(&ListOfBondsPerAtom[i]);
    Free(&ListOfBondsPerAtom);
  }
  if (NumberOfBondsPerAtom != NULL)
    Free(&NumberOfBondsPerAtom);
  ListOfBondsPerAtom = Malloc<bond**>(AtomCount, "molecule::CreateListOfBondsPerAtom: ***ListOfBondsPerAtom");
  NumberOfBondsPerAtom = Malloc<int>(AtomCount, "molecule::CreateListOfBondsPerAtom: *NumberOfBondsPerAtom");

  // reset bond counts per atom
  for(int i=AtomCount;i--;)
    NumberOfBondsPerAtom[i] = 0;
  // count bonds per atom
  Binder = first;
  while (Binder->next != last) {
    Binder = Binder->next;
    NumberOfBondsPerAtom[Binder->leftatom->nr]++;
    NumberOfBondsPerAtom[Binder->rightatom->nr]++;
  }
  for(int i=AtomCount;i--;) {
    // allocate list of bonds per atom
    ListOfBondsPerAtom[i] = Malloc<bond*>(NumberOfBondsPerAtom[i], "molecule::CreateListOfBondsPerAtom: **ListOfBondsPerAtom[]");
    // clear the list again, now each NumberOfBondsPerAtom marks current free field
    NumberOfBondsPerAtom[i] = 0;
  }
  // fill the list
  Binder = first;
  while (Binder->next != last) {
    Binder = Binder->next;
    ListOfBondsPerAtom[Binder->leftatom->nr][NumberOfBondsPerAtom[Binder->leftatom->nr]++] = Binder;
    ListOfBondsPerAtom[Binder->rightatom->nr][NumberOfBondsPerAtom[Binder->rightatom->nr]++] = Binder;
  }

  // output list for debugging
  *out << Verbose(3) << "ListOfBondsPerAtom for each atom:" << endl;
  Walker = start;
  while (Walker->next != end) {
    Walker = Walker->next;
    *out << Verbose(4) << "Atom " << Walker->Name << "/" << Walker->nr << " with " << NumberOfBondsPerAtom[Walker->nr] << " bonds: ";
    TotalDegree = 0;
    for (int j=0;j<NumberOfBondsPerAtom[Walker->nr];j++) {
      *out << *ListOfBondsPerAtom[Walker->nr][j] << "\t";
      TotalDegree += ListOfBondsPerAtom[Walker->nr][j]->BondDegree;
    }
    *out << " -- TotalDegree: " << TotalDegree << endl;
  }
  *out << Verbose(1) << "End of Creating ListOfBondsPerAtom." << endl << endl;
};


/** Blows the 6-dimensional \a cell_size array up to a full NDIM by NDIM matrix.
 * \param *symm 6-dim array of unique symmetric matrix components
 * \return allocated NDIM*NDIM array with the symmetric matrix
 */
double * molecule::ReturnFullMatrixforSymmetric(double *symm)
{
  double *matrix = Malloc<double>(NDIM * NDIM, "molecule::ReturnFullMatrixforSymmetric: *matrix");
  matrix[0] = symm[0];
  matrix[1] = symm[1];
  matrix[2] = symm[3];
  matrix[3] = symm[1];
  matrix[4] = symm[2];
  matrix[5] = symm[4];
  matrix[6] = symm[3];
  matrix[7] = symm[4];
  matrix[8] = symm[5];
  return matrix;
};


/** Comparison function for GSL heapsort on distances in two molecules.
 * \param *a
 * \param *b
 * \return <0, \a *a less than \a *b, ==0 if equal, >0 \a *a greater than \a *b
 */
inline int CompareDoubles (const void * a, const void * b)
{
  if (*(double *)a > *(double *)b)
    return -1;
  else if (*(double *)a < *(double *)b)
    return 1;
  else
    return 0;
};

/** Determines whether two molecules actually contain the same atoms and coordination.
 * \param *out output stream for debugging
 * \param *OtherMolecule the molecule to compare this one to
 * \param threshold upper limit of difference when comparing the coordination.
 * \return NULL - not equal, otherwise an allocated (molecule::AtomCount) permutation map of the atom numbers (which corresponds to which)
 */
int * molecule::IsEqualToWithinThreshold(ofstream *out, molecule *OtherMolecule, double threshold)
{
  int flag;
  double *Distances = NULL, *OtherDistances = NULL;
  Vector CenterOfGravity, OtherCenterOfGravity;
  size_t *PermMap = NULL, *OtherPermMap = NULL;
  int *PermutationMap = NULL;
  atom *Walker = NULL;
  bool result = true; // status of comparison

  *out << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl;
  /// first count both their atoms and elements and update lists thereby ...
  //*out << Verbose(0) << "Counting atoms, updating list" << endl;
  CountAtoms(out);
  OtherMolecule->CountAtoms(out);
  CountElements();
  OtherMolecule->CountElements();

  /// ... and compare:
  /// -# AtomCount
  if (result) {
    if (AtomCount != OtherMolecule->AtomCount) {
      *out << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
      result = false;
    } else *out << Verbose(4) << "AtomCounts match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
  }
  /// -# ElementCount
  if (result) {
    if (ElementCount != OtherMolecule->ElementCount) {
      *out << Verbose(4) << "ElementCount don't match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;
      result = false;
    } else *out << Verbose(4) << "ElementCount match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;
  }
  /// -# ElementsInMolecule
  if (result) {
    for (flag=MAX_ELEMENTS;flag--;) {
      //*out << Verbose(5) << "Element " <<  flag << ": " << ElementsInMolecule[flag] << " <-> " << OtherMolecule->ElementsInMolecule[flag] << "." << endl;
      if (ElementsInMolecule[flag] != OtherMolecule->ElementsInMolecule[flag])
        break;
    }
    if (flag < MAX_ELEMENTS) {
      *out << Verbose(4) << "ElementsInMolecule don't match." << endl;
      result = false;
    } else *out << Verbose(4) << "ElementsInMolecule match." << endl;
  }
  /// then determine and compare center of gravity for each molecule ...
  if (result) {
    *out << Verbose(5) << "Calculating Centers of Gravity" << endl;
    DeterminePeriodicCenter(CenterOfGravity);
    OtherMolecule->DeterminePeriodicCenter(OtherCenterOfGravity);
    *out << Verbose(5) << "Center of Gravity: ";
    CenterOfGravity.Output(out);
    *out << endl << Verbose(5) << "Other Center of Gravity: ";
    OtherCenterOfGravity.Output(out);
    *out << endl;
    if (CenterOfGravity.DistanceSquared(&OtherCenterOfGravity) > threshold*threshold) {
      *out << Verbose(4) << "Centers of gravity don't match." << endl;
      result = false;
    }
  }

  /// ... then make a list with the euclidian distance to this center for each atom of both molecules
  if (result) {
    *out << Verbose(5) << "Calculating distances" << endl;
    Distances = Malloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: Distances");
    OtherDistances = Malloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: OtherDistances");
    Walker = start;
    while (Walker->next != end) {
      Walker = Walker->next;
      Distances[Walker->nr] = CenterOfGravity.DistanceSquared(&Walker->x);
    }
    Walker = OtherMolecule->start;
    while (Walker->next != OtherMolecule->end) {
      Walker = Walker->next;
      OtherDistances[Walker->nr] = OtherCenterOfGravity.DistanceSquared(&Walker->x);
    }

    /// ... sort each list (using heapsort (o(N log N)) from GSL)
    *out << Verbose(5) << "Sorting distances" << endl;
    PermMap = Malloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermMap");
    OtherPermMap = Malloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *OtherPermMap");
    gsl_heapsort_index (PermMap, Distances, AtomCount, sizeof(double), CompareDoubles);
    gsl_heapsort_index (OtherPermMap, OtherDistances, AtomCount, sizeof(double), CompareDoubles);
    PermutationMap = Malloc<int>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermutationMap");
    *out << Verbose(5) << "Combining Permutation Maps" << endl;
    for(int i=AtomCount;i--;)
      PermutationMap[PermMap[i]] = (int) OtherPermMap[i];

    /// ... and compare them step by step, whether the difference is individually(!) below \a threshold for all
    *out << Verbose(4) << "Comparing distances" << endl;
    flag = 0;
    for (int i=0;i<AtomCount;i++) {
      *out << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " <<  threshold << endl;
      if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold*threshold)
        flag = 1;
    }

    // free memory
    Free(&PermMap);
    Free(&OtherPermMap);
    Free(&Distances);
    Free(&OtherDistances);
    if (flag) { // if not equal
      Free(&PermutationMap);
      result = false;
    }
  }
  /// return pointer to map if all distances were below \a threshold
  *out << Verbose(3) << "End of IsEqualToWithinThreshold." << endl;
  if (result) {
    *out << Verbose(3) << "Result: Equal." << endl;
    return PermutationMap;
  } else {
    *out << Verbose(3) << "Result: Not equal." << endl;
    return NULL;
  }
};

/** Returns an index map for two father-son-molecules.
 * The map tells which atom in this molecule corresponds to which one in the other molecul with their fathers.
 * \param *out output stream for debugging
 * \param *OtherMolecule corresponding molecule with fathers
 * \return allocated map of size molecule::AtomCount with map
 * \todo make this with a good sort O(n), not O(n^2)
 */
int * molecule::GetFatherSonAtomicMap(ofstream *out, molecule *OtherMolecule)
{
  atom *Walker = NULL, *OtherWalker = NULL;
  *out << Verbose(3) << "Begin of GetFatherAtomicMap." << endl;
  int *AtomicMap = Malloc<int>(AtomCount, "molecule::GetAtomicMap: *AtomicMap");
  for (int i=AtomCount;i--;)
    AtomicMap[i] = -1;
  if (OtherMolecule == this) {  // same molecule
    for (int i=AtomCount;i--;) // no need as -1 means already that there is trivial correspondence
      AtomicMap[i] = i;
    *out << Verbose(4) << "Map is trivial." << endl;
  } else {
    *out << Verbose(4) << "Map is ";
    Walker = start;
    while (Walker->next != end) {
      Walker = Walker->next;
      if (Walker->father == NULL) {
        AtomicMap[Walker->nr] = -2;
      } else {
        OtherWalker = OtherMolecule->start;
        while (OtherWalker->next != OtherMolecule->end) {
          OtherWalker = OtherWalker->next;
      //for (int i=0;i<AtomCount;i++) { // search atom
        //for (int j=0;j<OtherMolecule->AtomCount;j++) {
          //*out << Verbose(4) << "Comparing father " << Walker->father << " with the other one " << OtherWalker->father << "." << endl;
          if (Walker->father == OtherWalker)
            AtomicMap[Walker->nr] = OtherWalker->nr;
        }
      }
      *out << AtomicMap[Walker->nr] << "\t";
    }
    *out << endl;
  }
  *out << Verbose(3) << "End of GetFatherAtomicMap." << endl;
  return AtomicMap;
};

/** Stores the temperature evaluated from velocities in molecule::Trajectories.
 * We simply use the formula equivaleting temperature and kinetic energy:
 * \f$k_B T = \sum_i m_i v_i^2\f$
 * \param *out output stream for debugging
 * \param startstep first MD step in molecule::Trajectories
 * \param endstep last plus one MD step in molecule::Trajectories
 * \param *output output stream of temperature file
 * \return file written (true), failure on writing file (false)
 */
bool molecule::OutputTemperatureFromTrajectories(ofstream *out, int startstep, int endstep, ofstream *output)
{
  double temperature;
  atom *Walker = NULL;
  // test stream
  if (output == NULL)
    return false;
  else
    *output << "# Step Temperature [K] Temperature [a.u.]" << endl;
  for (int step=startstep;step < endstep; step++) { // loop over all time steps
    temperature = 0.;
    Walker = start;
    while (Walker->next != end) {
      Walker = Walker->next;
      for (int i=NDIM;i--;)
        temperature += Walker->type->mass * Trajectories[Walker].U.at(step).x[i]* Trajectories[Walker].U.at(step).x[i];
    }
    *output << step << "\t" << temperature*AtomicEnergyToKelvin << "\t" << temperature << endl;
  }
  return true;
};