source: src/Analysis/analysis_correlation.cpp@ a58c16

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

/*
 * analysis.cpp
 *
 *  Created on: Oct 13, 2009
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include <algorithm>
#include <iostream>
#include <iomanip>
#include <limits>

#include "Atom/atom.hpp"
#include "Bond/bond.hpp"
#include "Tesselation/BoundaryTriangleSet.hpp"
#include "Box.hpp"
#include "Element/element.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "Descriptors/AtomOfMoleculeSelectionDescriptor.hpp"
#include "Descriptors/MoleculeFormulaDescriptor.hpp"
#include "Descriptors/MoleculeOfAtomSelectionDescriptor.hpp"
#include "Formula.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "LinkedCell/LinkedCell_View.hpp"
#include "molecule.hpp"
#include "Tesselation/tesselation.hpp"
#include "Tesselation/tesselationhelpers.hpp"
#include "Tesselation/triangleintersectionlist.hpp"
#include "World.hpp"
#include "WorldTime.hpp"

#include "analysis_correlation.hpp"

/** Calculates the dipole vector of a given atomSet.
 *
 *  Note that we use the following procedure as rule of thumb:
 *   -# go through every bond of the atom
 *   -# calculate the difference of electronegativities \f$\Delta\mathrm{EN}\f$
 *   -# if \f$\Delta\mathrm{EN} > 0.5\f$, we align the bond vector in direction of the more negative element
 *   -# sum up all vectors
 *   -# finally, divide by the number of summed vectors
 *
 * @param atomsbegin begin iterator of atomSet
 * @param atomsend end iterator of atomset
 * @return dipole vector
 */
Vector getDipole(molecule::const_iterator atomsbegin, molecule::const_iterator atomsend)
{
  Vector DipoleVector;
  size_t SumOfVectors = 0;
  Box &domain = World::getInstance().getDomain();

      // go through all atoms
  for (molecule::const_iterator atomiter = atomsbegin;
      atomiter != atomsend;
      ++atomiter) {
    // go through all bonds
    const BondList& ListOfBonds = (*atomiter)->getListOfBonds();
    ASSERT(ListOfBonds.begin() != ListOfBonds.end(),
        "getDipole() - no bonds in molecule!");
    for (BondList::const_iterator bonditer = ListOfBonds.begin();
        bonditer != ListOfBonds.end();
        ++bonditer) {
      const atom * Otheratom = (*bonditer)->GetOtherAtom(*atomiter);
      if (Otheratom->getId() > (*atomiter)->getId()) {
        const double DeltaEN = (*atomiter)->getType()->getElectronegativity()
            -Otheratom->getType()->getElectronegativity();
        // get distance and correct for boundary conditions
        Vector BondDipoleVector = domain.periodicDistanceVector(
            (*atomiter)->getPosition(),
            Otheratom->getPosition());
        // DeltaEN is always positive, gives correct orientation of vector
        BondDipoleVector.Normalize();
        BondDipoleVector *= DeltaEN;
        LOG(3,"INFO: Dipole vector from bond " << **bonditer << " is " << BondDipoleVector);
        DipoleVector += BondDipoleVector;
        SumOfVectors++;
      }
    }
  }
  LOG(3,"INFO: Sum over all bond dipole vectors is "
      << DipoleVector << " with " << SumOfVectors << " in total.");
  if (SumOfVectors != 0)
    DipoleVector *= 1./(double)SumOfVectors;
  LOG(2, "INFO: Resulting dipole vector is " << DipoleVector);

  return DipoleVector;
};

/** Calculate minimum and maximum amount of trajectory steps by going through given atomic trajectories.
 * \param vector of atoms whose trajectories to check for [min,max]
 * \return range with [min, max]
 */
range<size_t> getMaximumTrajectoryBounds(const std::vector<atom *> &atoms)
{
  // get highest trajectory size
  LOG(0,"STATUS: Retrieving maximum amount of time steps ...");
  if (atoms.size() == 0)
    return range<size_t>(0,0);
  size_t max_timesteps = std::numeric_limits<size_t>::min();
  size_t min_timesteps = std::numeric_limits<size_t>::max();
  BOOST_FOREACH(atom *_atom, atoms) {
    if (_atom->getTrajectorySize() > max_timesteps)
      max_timesteps  = _atom->getTrajectorySize();
    if (_atom->getTrajectorySize() < min_timesteps)
      min_timesteps = _atom->getTrajectorySize();
  }
  LOG(1,"INFO: Minimum number of time steps found is " << min_timesteps);
  LOG(1,"INFO: Maximum number of time steps found is " << max_timesteps);

  return range<size_t>(min_timesteps, max_timesteps);
}

/** Calculates the angular dipole zero orientation from current time step.
 * \param molecules vector of molecules to calculate dipoles of
 * \return map with orientation vector for each atomic id given in \a atoms.
 */
std::map<atomId_t, Vector> CalculateZeroAngularDipole(const std::vector<molecule *> &molecules)
{
  // get zero orientation for each molecule.
  LOG(0,"STATUS: Calculating dipoles for current time step ...");
  std::map<atomId_t, Vector> ZeroVector;
  BOOST_FOREACH(molecule *_mol, molecules) {
    const Vector Dipole =
        getDipole(
            const_cast<const molecule *>(_mol)->begin(),
            const_cast<const molecule *>(_mol)->end());
    for(molecule::const_iterator iter = const_cast<const molecule *>(_mol)->begin();
        iter != const_cast<const molecule *>(_mol)->end();
        ++iter)
      ZeroVector[(*iter)->getId()] = Dipole;
    LOG(2,"INFO: Zero alignment for molecule " << _mol->getId() << " is " << Dipole);
  }
  LOG(1,"INFO: We calculated zero orientation for a total of " << molecules.size() << " molecule(s).");

  return ZeroVector;
}

/** Calculates the dipole angular correlation for given molecule type.
 * Calculate the change of the dipole orientation angle over time.
 * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
 * Angles are given in degrees.
 * \param &atoms list of atoms of the molecules taking part (Note: molecules may
 * change over time as bond structure is recalculated, hence we need the atoms)
 * \param timestep time step to calculate angular correlation for (relative to
 *  \a ZeroVector)
 * \param ZeroVector map with Zero orientation vector for each atom in \a atoms.
 * \param DontResetTime don't reset time to old value (triggers re-creation of bond system)
 * \return Map of doubles with values the pair of the two atoms.
 */
DipoleAngularCorrelationMap *DipoleAngularCorrelation(
    const Formula &DipoleFormula,
    const size_t timestep,
    const std::map<atomId_t, Vector> &ZeroVector,
    const enum ResetWorldTime DoTimeReset
    )
{
  Info FunctionInfo(__func__);
  DipoleAngularCorrelationMap *outmap = new DipoleAngularCorrelationMap;

  unsigned int oldtime = 0;
  if (DoTimeReset == DoResetTime) {
    // store original time step
    oldtime = WorldTime::getTime();
  }

  // set time step
  LOG(0,"STATUS: Stepping onto to time step " << timestep << ".");
  World::getInstance().setTime(timestep);

  // get all molecules for this time step
  World::getInstance().clearMoleculeSelection();
  World::getInstance().selectAllMolecules(MoleculeByFormula(DipoleFormula));
  std::vector<molecule *> molecules = World::getInstance().getSelectedMolecules();
  LOG(1,"INFO: There are " << molecules.size() << " molecules for time step " << timestep << ".");

  // calculate dipoles for each
  LOG(0,"STATUS: Calculating dipoles for time step " << timestep << " ...");
  size_t i=0;
  size_t Counter_rejections = 0;
  BOOST_FOREACH(molecule *_mol, molecules) {
    const Vector Dipole =
        getDipole(
            const_cast<const molecule *>(_mol)->begin(),
            const_cast<const molecule *>(_mol)->end());
    LOG(3,"INFO: Dipole vector at time step " << timestep << " for for molecule "
        << _mol->getId() << " is " << Dipole);
    // check that all atoms are valid (zeroVector known)
    molecule::const_iterator iter = const_cast<const molecule *>(_mol)->begin();
    for(; iter != const_cast<const molecule *>(_mol)->end(); ++iter) {
      if (!ZeroVector.count((*iter)->getId()))
        break;
    }
    if (iter != const_cast<const molecule *>(_mol)->end()) {
      ELOG(2, "Skipping molecule " << _mol->getName() << " as not all atoms have a valid zeroVector.");
      ++Counter_rejections;
      continue;
    } else
      iter = const_cast<const molecule *>(_mol)->begin();
    std::map<atomId_t, Vector>::const_iterator zeroValue = ZeroVector.find((*iter)->getId()); //due to iter is const
    double angle = 0.;
    LOG(2, "INFO: ZeroVector of first atom " << **iter << " is "
        << zeroValue->second << ".");
    LOG(4, "INFO: Squared norm of difference vector is "
        << (zeroValue->second - Dipole).NormSquared() << ".");
    if ((zeroValue->second - Dipole).NormSquared() > MYEPSILON)
      angle = Dipole.Angle(zeroValue->second) * (180./M_PI);
    else
      LOG(2, "INFO: Both vectors (almost) coincide, numerically unstable, angle set to zero.");
    // we print six digits, hence round up to six digit precision
    const double precision = 1e-6;
    angle = precision*floor(angle/precision);
    LOG(1,"INFO: Resulting relative angle for molecule " << _mol->getName()
        << " is " << angle << ".");
    outmap->insert ( std::make_pair (angle, *iter ) );
    ++i;
  }
  ASSERT(Counter_rejections <= molecules.size(),
      "DipoleAngularCorrelation() - more rejections ("+toString(Counter_rejections)
      +") than there are molecules ("+toString(molecules.size())+").");
  LOG(1,"INFO: " << Counter_rejections << " molecules have been rejected in time step " << timestep << ".");

  LOG(0,"STATUS: Done with calculating dipoles.");

  if (DoTimeReset == DoResetTime) {
    // re-set to original time step again
    World::getInstance().setTime(oldtime);
  }

  // and return results
  return outmap;
};

/** Calculates the dipole correlation for given molecule type.
 * I.e. we calculate how the angle between any two given dipoles in the
 * systems behaves. Sort of pair correlation but distance is replaced by
 * the orientation distance, i.e. an angle.
 * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
 * Angles are given in degrees.
 * \param *molecules vector of molecules
 * \return Map of doubles with values the pair of the two atoms.
 */
DipoleCorrelationMap *DipoleCorrelation(std::vector<molecule *> &molecules)
{
  Info FunctionInfo(__func__);
  DipoleCorrelationMap *outmap = new DipoleCorrelationMap;
//  double distance = 0.;
//  Box &domain = World::getInstance().getDomain();
//
  if (molecules.empty()) {
    ELOG(1, "No molecule given.");
    return outmap;
  }

  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin();
      MolWalker != molecules.end(); ++MolWalker) {
    LOG(2, "INFO: Current molecule is " << (*MolWalker)->getId() << ".");
    const Vector Dipole =
        getDipole(
            const_cast<const molecule *>(*MolWalker)->begin(),
            const_cast<const molecule *>(*MolWalker)->end());
    std::vector<molecule *>::const_iterator MolOtherWalker = MolWalker;
    for (++MolOtherWalker;
        MolOtherWalker != molecules.end();
        ++MolOtherWalker) {
      LOG(2, "INFO: Current other molecule is " << (*MolOtherWalker)->getId() << ".");
      const Vector OtherDipole = getDipole(
          const_cast<const molecule *>(*MolOtherWalker)->begin(),
          const_cast<const molecule *>(*MolOtherWalker)->end());
      const double angle = Dipole.Angle(OtherDipole) * (180./M_PI);
      LOG(1, "Angle is " << angle << ".");
      outmap->insert ( make_pair (angle, make_pair ((*MolWalker), (*MolOtherWalker)) ) );
    }
  }
  return outmap;
};

/** Calculates the pair correlation between given atom sets.
 *
 * Note we correlate each of the \a &atomsfirst with each of the second set
 * \a &atoms_second. However, we are aware of double counting. If an atom is
 * in either set, the pair is counted only once.
 *
 * \param &atoms_first vector of atoms
 * \param &atoms_second vector of atoms
 * \param max_distance maximum distance for the correlation
 * \return Map of doubles with values the pair of the two atoms.
 */
PairCorrelationMap *PairCorrelation(
    const World::ConstAtomComposite &atoms_first,
    const World::ConstAtomComposite &atoms_second,
    const double max_distance)
{
  Info FunctionInfo(__func__);
  PairCorrelationMap *outmap = new PairCorrelationMap;
  //double distance = 0.;
  Box &domain = World::getInstance().getDomain();

  if (atoms_first.empty() || atoms_second.empty()) {
    ELOG(1, "No atoms given.");
    return outmap;
  }

  //!> typedef for an unsorted container, (output) compatible with STL algorithms
  typedef std::vector<const TesselPoint *> LinkedVector;

  // create intersection (to know when to check for double-counting)
  LinkedVector intersected_atoms(atoms_second.size(), NULL);
  LinkedVector::iterator intersected_atoms_end =
      std::set_intersection(
          atoms_first.begin(),atoms_first.end(),
          atoms_second.begin(), atoms_second.end(),
          intersected_atoms.begin());
  const LinkedCell::LinkedList intersected_atoms_set(intersected_atoms.begin(), intersected_atoms_end);

  // get linked cell view
  LinkedCell::LinkedCell_View LC = World::getInstance().getLinkedCell(max_distance);

  // convert second to _sorted_ set
  LinkedCell::LinkedList atoms_second_set(atoms_second.begin(), atoms_second.end());
  LOG(2, "INFO: first set has " << atoms_first.size()
      << " and second set has " << atoms_second_set.size() << " atoms.");

  // fill map
  for (World::ConstAtomComposite::const_iterator iter = atoms_first.begin();
      iter != atoms_first.end();
      ++iter) {
    const TesselPoint * const Walker = *iter;
    LOG(3, "INFO: Current point is " << Walker->getName() << ".");
    // obtain all possible neighbors (that is a sorted set)
    LinkedCell::LinkedList ListOfNeighbors = LC.getPointsInsideSphere(
        max_distance,
        Walker->getPosition());
    LOG(2, "INFO: There are " << ListOfNeighbors.size() << " neighbors.");

    // create intersection with second set
    // NOTE: STL algorithms do mostly not work on sorted container because reassignment
    // of a value may also require changing its position.
    LinkedVector intersected_set(atoms_second.size(), NULL);
    LinkedVector::iterator intersected_end =
        std::set_intersection(
            ListOfNeighbors.begin(),ListOfNeighbors.end(),
            atoms_second_set.begin(), atoms_second_set.end(),
            intersected_set.begin());
    // count remaining elements
    LOG(2, "INFO: Intersection with second set has " << int(intersected_end - intersected_set.begin()) << " elements.");
    // we have some possible candidates, go through each
    for (LinkedVector::const_iterator neighboriter = intersected_set.begin();
        neighboriter != intersected_end;
        ++neighboriter) {
      const TesselPoint * const OtherWalker = (*neighboriter);
      LinkedCell::LinkedList::const_iterator equaliter = intersected_atoms_set.find(OtherWalker);
      if ((equaliter !=  intersected_atoms_set.end()) && (OtherWalker <= Walker)) {
        // present in both sets, assure that we are larger
        continue;
      }
      LOG(3, "INFO: Current other point is " << *OtherWalker << ".");
      const double distance = domain.periodicDistance(OtherWalker->getPosition(),Walker->getPosition());
      LOG(3, "INFO: Resulting distance is " << distance << ".");
      outmap->insert (
          std::pair<double, std::pair <const TesselPoint *, const TesselPoint*> > (
              distance,
              std::make_pair (Walker, OtherWalker)
              )
          );
    }
  }
  // and return
  return outmap;
};

/** Calculates the distance (pair) correlation between a given element and a point.
 * \param *molecules list of molecules structure
 * \param &elements vector of elements to correlate with point
 * \param *point vector to the correlation point
 * \return Map of dobules with values as pairs of atom and the vector
 */
CorrelationToPointMap *CorrelationToPoint(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Vector *point )
{
  Info FunctionInfo(__func__);
  CorrelationToPointMap *outmap = new CorrelationToPointMap;
  double distance = 0.;
  Box &domain = World::getInstance().getDomain();

  if (molecules.empty()) {
    LOG(1, "No molecule given.");
    return outmap;
  }

  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
    LOG(2, "Current molecule is " << *MolWalker << ".");
    for (molecule::const_iterator iter = const_cast<const molecule *>(*MolWalker)->begin();
        iter != const_cast<const molecule *>(*MolWalker)->end();
        ++iter) {
      LOG(3, "Current atom is " << **iter << ".");
      for (vector<const element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
        if ((*type == NULL) || ((*iter)->getType() == *type)) {
          distance = domain.periodicDistance((*iter)->getPosition(),*point);
          LOG(4, "Current distance is " << distance << ".");
          outmap->insert (
              std::pair<double, std::pair<const atom *, const Vector*> >(
                  distance,
                  std::pair<const atom *, const Vector*> (
                      (*iter),
                      point)
                  )
              );
        }
    }
  }

  return outmap;
};

/** Calculates the distance (pair) correlation between a given element, all its periodic images and a point.
 * \param *molecules list of molecules structure
 * \param &elements vector of elements to correlate to point
 * \param *point vector to the correlation point
 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
 * \return Map of dobules with values as pairs of atom and the vector
 */
CorrelationToPointMap *PeriodicCorrelationToPoint(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Vector *point, const int ranges[NDIM] )
{
  Info FunctionInfo(__func__);
  CorrelationToPointMap *outmap = new CorrelationToPointMap;
  double distance = 0.;
  int n[NDIM];
  Vector periodicX;
  Vector checkX;

  if (molecules.empty()) {
    LOG(1, "No molecule given.");
    return outmap;
  }

  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
    RealSpaceMatrix FullMatrix = World::getInstance().getDomain().getM();
    RealSpaceMatrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
    LOG(2, "Current molecule is " << *MolWalker << ".");
    for (molecule::const_iterator iter = const_cast<const molecule *>(*MolWalker)->begin();
        iter != const_cast<const molecule *>(*MolWalker)->end();
        ++iter) {
      LOG(3, "Current atom is " << **iter << ".");
      for (vector<const element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
        if ((*type == NULL) || ((*iter)->getType() == *type)) {
          periodicX = FullInverseMatrix * ((*iter)->getPosition()); // x now in [0,1)^3
          // go through every range in xyz and get distance
          for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
            for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
              for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
                checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
                distance = checkX.distance(*point);
                LOG(4, "Current distance is " << distance << ".");
                outmap->insert (
                    std::pair<double,
                    std::pair<const atom *, const Vector*> >(
                        distance,
                        std::pair<const atom *, const Vector*> (
                            *iter,
                            point)
                        )
                    );
              }
        }
    }
  }

  return outmap;
};

/** Calculates the distance (pair) correlation between a given element and a surface.
 * \param *molecules list of molecules structure
 * \param &elements vector of elements to correlate to surface
 * \param *Surface pointer to Tesselation class surface
 * \param *LC LinkedCell_deprecated structure to quickly find neighbouring atoms
 * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
 */
CorrelationToSurfaceMap *CorrelationToSurface(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Tesselation * const Surface, const LinkedCell_deprecated *LC )
{
  Info FunctionInfo(__func__);
  CorrelationToSurfaceMap *outmap = new CorrelationToSurfaceMap;
  double distance = 0;
  class BoundaryTriangleSet *triangle = NULL;
  Vector centroid;

  if ((Surface == NULL) || (LC == NULL) || (molecules.empty())) {
    ELOG(1, "No Tesselation, no LinkedCell or no molecule given.");
    return outmap;
  }

  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
    LOG(2, "Current molecule is " << (*MolWalker)->name << ".");
    if ((*MolWalker)->empty())
      LOG(2, "\t is empty.");
    for (molecule::const_iterator iter = const_cast<const molecule *>(*MolWalker)->begin();
        iter != const_cast<const molecule *>(*MolWalker)->end();
        ++iter) {
      LOG(3, "\tCurrent atom is " << *(*iter) << ".");
      for (vector<const element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
        if ((*type == NULL) || ((*iter)->getType() == *type)) {
          TriangleIntersectionList Intersections((*iter)->getPosition(),Surface,LC);
          distance = Intersections.GetSmallestDistance();
          triangle = Intersections.GetClosestTriangle();
          outmap->insert (
              std::pair<double,
              std::pair<const atom *, BoundaryTriangleSet*> >(
                  distance,
                  std::pair<const atom *, BoundaryTriangleSet*> (
                      (*iter),
                      triangle)
                  )
              );
        }
    }
  }

  return outmap;
};

/** Calculates the distance (pair) correlation between a given element, all its periodic images and and a surface.
 * Note that we also put all periodic images found in the cells given by [ -ranges[i], ranges[i] ] and i=0,...,NDIM-1.
 * I.e. We multiply the atom::node with the inverse of the domain matrix, i.e. transform it to \f$[0,0^3\f$, then add per
 * axis an integer from [ -ranges[i], ranges[i] ] onto it and multiply with the domain matrix to bring it back into
 * the real space. Then, we Tesselation::FindClosestTriangleToPoint() and DistanceToTrianglePlane().
 * \param *molecules list of molecules structure
 * \param &elements vector of elements to correlate to surface
 * \param *Surface pointer to Tesselation class surface
 * \param *LC LinkedCell_deprecated structure to quickly find neighbouring atoms
 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
 * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
 */
CorrelationToSurfaceMap *PeriodicCorrelationToSurface(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Tesselation * const Surface, const LinkedCell_deprecated *LC, const int ranges[NDIM] )
{
  Info FunctionInfo(__func__);
  CorrelationToSurfaceMap *outmap = new CorrelationToSurfaceMap;
  double distance = 0;
  class BoundaryTriangleSet *triangle = NULL;
  Vector centroid;
  int n[NDIM];
  Vector periodicX;
  Vector checkX;

  if ((Surface == NULL) || (LC == NULL) || (molecules.empty())) {
    LOG(1, "No Tesselation, no LinkedCell or no molecule given.");
    return outmap;
  }

  double ShortestDistance = 0.;
  BoundaryTriangleSet *ShortestTriangle = NULL;
  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
    RealSpaceMatrix FullMatrix = World::getInstance().getDomain().getM();
    RealSpaceMatrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
    LOG(2, "Current molecule is " << *MolWalker << ".");
    for (molecule::const_iterator iter = const_cast<const molecule *>(*MolWalker)->begin();
        iter != const_cast<const molecule *>(*MolWalker)->end();
        ++iter) {
      LOG(3, "Current atom is " << **iter << ".");
      for (vector<const element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
        if ((*type == NULL) || ((*iter)->getType() == *type)) {
          periodicX = FullInverseMatrix * ((*iter)->getPosition()); // x now in [0,1)^3
          // go through every range in xyz and get distance
          ShortestDistance = -1.;
          for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
            for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
              for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
                checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
                TriangleIntersectionList Intersections(checkX,Surface,LC);
                distance = Intersections.GetSmallestDistance();
                triangle = Intersections.GetClosestTriangle();
                if ((ShortestDistance == -1.) || (distance < ShortestDistance)) {
                  ShortestDistance = distance;
                  ShortestTriangle = triangle;
                }
              }
          // insert
          outmap->insert (
              std::pair<double,
              std::pair<const atom *, BoundaryTriangleSet*> >(
                  ShortestDistance,
                  std::pair<const atom *, BoundaryTriangleSet*> (
                      *iter,
                      ShortestTriangle)
                  )
              );
          //LOG(1, "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << ".");
        }
    }
  }

  return outmap;
};

/** Returns the index of the bin for a given value.
 * \param value value whose bin to look for
 * \param BinWidth width of bin
 * \param BinStart first bin
 */
int GetBin ( const double value, const double BinWidth, const double BinStart )
{
  //Info FunctionInfo(__func__);
  int bin =(int) (floor((value - BinStart)/BinWidth));
  return (bin);
};


/** Adds header part that is unique to BinPairMap.
 *
 * @param file stream to print to
 */
void OutputCorrelation_Header( ofstream * const file )
{
  *file << "\tCount";
};

/** Prints values stored in BinPairMap iterator.
 *
 * @param file stream to print to
 * @param runner iterator pointing at values to print
 */
void OutputCorrelation_Value( ofstream * const file, BinPairMap::const_iterator &runner )
{
  *file << runner->second;
};


/** Adds header part that is unique to DipoleAngularCorrelationMap.
 *
 * @param file stream to print to
 */
void OutputDipoleAngularCorrelation_Header( ofstream * const file )
{
  *file << "\tFirstAtomOfMolecule";
};

/** Prints values stored in DipoleCorrelationMap iterator.
 *
 * @param file stream to print to
 * @param runner iterator pointing at values to print
 */
void OutputDipoleAngularCorrelation_Value( ofstream * const file, DipoleAngularCorrelationMap::const_iterator &runner )
{
  *file << *(runner->second);
};


/** Adds header part that is unique to DipoleAngularCorrelationMap.
 *
 * @param file stream to print to
 */
void OutputDipoleCorrelation_Header( ofstream * const file )
{
  *file << "\tMolecule";
};

/** Prints values stored in DipoleCorrelationMap iterator.
 *
 * @param file stream to print to
 * @param runner iterator pointing at values to print
 */
void OutputDipoleCorrelation_Value( ofstream * const file, DipoleCorrelationMap::const_iterator &runner )
{
  *file << runner->second.first->getId() << "\t" << runner->second.second->getId();
};


/** Adds header part that is unique to PairCorrelationMap.
 *
 * @param file stream to print to
 */
void OutputPairCorrelation_Header( ofstream * const file )
{
  *file << "\tAtom1\tAtom2";
};

/** Prints values stored in PairCorrelationMap iterator.
 *
 * @param file stream to print to
 * @param runner iterator pointing at values to print
 */
void OutputPairCorrelation_Value( ofstream * const file, PairCorrelationMap::const_iterator &runner )
{
  *file << *(runner->second.first) << "\t" << *(runner->second.second);
};


/** Adds header part that is unique to CorrelationToPointMap.
 *
 * @param file stream to print to
 */
void OutputCorrelationToPoint_Header( ofstream * const file )
{
  *file << "\tAtom::x[i]-point.x[i]";
};

/** Prints values stored in CorrelationToPointMap iterator.
 *
 * @param file stream to print to
 * @param runner iterator pointing at values to print
 */
void OutputCorrelationToPoint_Value( ofstream * const file, CorrelationToPointMap::const_iterator &runner )
{
  for (int i=0;i<NDIM;i++)
    *file << "\t" << setprecision(8) << (runner->second.first->at(i) - runner->second.second->at(i));
};


/** Adds header part that is unique to CorrelationToSurfaceMap.
 *
 * @param file stream to print to
 */
void OutputCorrelationToSurface_Header( ofstream * const file )
{
  *file << "\tTriangle";
};

/** Prints values stored in CorrelationToSurfaceMap iterator.
 *
 * @param file stream to print to
 * @param runner iterator pointing at values to print
 */
void OutputCorrelationToSurface_Value( ofstream * const file, CorrelationToSurfaceMap::const_iterator &runner )
{
  *file << *(runner->second.first) << "\t" << *(runner->second.second);
};