/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2011 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

/*
 * 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 <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 "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;
  // 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();
        Vector BondDipoleVector = (*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(1, "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(_mol->begin(), _mol->end());
    for(molecule::const_iterator iter = _mol->begin(); iter != _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(_mol->begin(), _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 = _mol->begin();
    for(; iter != _mol->end(); ++iter) {
      if (!ZeroVector.count((*iter)->getId()))
        break;
    }
    if (iter != _mol->end()) {
      ELOG(2, "Skipping molecule " << _mol->getName() << " as not all atoms have a valid zeroVector.");
      ++Counter_rejections;
      continue;
    } else
      iter = _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.");
    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((*MolWalker)->begin(), (*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((*MolOtherWalker)->begin(), (*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 elements.
 * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
 * \param *molecules vector of molecules
 * \param &elements vector of elements to correlate
 * \return Map of doubles with values the pair of the two atoms.
 */
PairCorrelationMap *PairCorrelation(std::vector<molecule *> &molecules, const std::vector<const element *> &elements)
{
  Info FunctionInfo(__func__);
  PairCorrelationMap *outmap = new PairCorrelationMap;
  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++)
    (*MolWalker)->doCountAtoms();

  // create all possible pairs of elements
  set <pair<const element *,const element *> > PairsOfElements;
  if (elements.size() >= 2) {
    for (vector<const element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
      for (vector<const element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
        if (type1 != type2) {
          PairsOfElements.insert( make_pair(*type1,*type2) );
          LOG(1, "Creating element pair " << *(*type1) << " and " << *(*type2) << ".");
        }
  } else if (elements.size() == 1) { // one to all are valid
    const element *elemental = *elements.begin();
    PairsOfElements.insert( pair<const element *,const element*>(elemental,0) );
    PairsOfElements.insert( pair<const element *,const element*>(0,elemental) );
  } else { // all elements valid
    PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
  }

  outmap = new PairCorrelationMap;
  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++){
    LOG(2, "Current molecule is " << (*MolWalker)->getName() << ".");
    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
      LOG(3, "Current atom is " << **iter << ".");
      for (std::vector<molecule *>::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules.end(); MolOtherWalker++){
        LOG(2, "Current other molecule is " << (*MolOtherWalker)->getName() << ".");
        for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
          LOG(3, "Current otheratom is " << **runner << ".");
          if ((*iter)->getId() < (*runner)->getId()){
            for (set <pair<const element *, const element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
              if ((PairRunner->first == (**iter).getType()) && (PairRunner->second == (**runner).getType())) {
                distance = domain.periodicDistance((*iter)->getPosition(),(*runner)->getPosition());
                //LOG(1, "Inserting " << *(*iter) << " and " << *(*runner));
                outmap->insert (
                    std::pair<double, std::pair <const atom *, const atom*> > (
                        distance,
                        std::pair<const atom *, const atom*> ((*iter), (*runner))
                        )
                    );
              }
          }
        }
      }
    }
  }
  return outmap;
};

/** Calculates the pair correlation between given elements.
 * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
 * \param *molecules list of molecules structure
 * \param &elements vector of elements to correlate
 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
 * \return Map of doubles with values the pair of the two atoms.
 */
PairCorrelationMap *PeriodicPairCorrelation(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const int ranges[NDIM] )
{
  Info FunctionInfo(__func__);
  PairCorrelationMap *outmap = new PairCorrelationMap;
  double distance = 0.;
  int n[NDIM];
  Vector checkX;
  Vector periodicX;
  int Othern[NDIM];
  Vector checkOtherX;
  Vector periodicOtherX;

  if (molecules.empty()) {
    ELOG(1, "No molecule given.");
    return outmap;
  }
  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
    (*MolWalker)->doCountAtoms();

  // create all possible pairs of elements
  set <pair<const element *,const element *> > PairsOfElements;
  if (elements.size() >= 2) {
    for (vector<const element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
      for (vector<const element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
        if (type1 != type2) {
          PairsOfElements.insert( make_pair(*type1,*type2) );
          LOG(1, "Creating element pair " << *(*type1) << " and " << *(*type2) << ".");
        }
  } else if (elements.size() == 1) { // one to all are valid
    const element *elemental = *elements.begin();
    PairsOfElements.insert( pair<const element *,const element*>(elemental,0) );
    PairsOfElements.insert( pair<const element *,const element*>(0,elemental) );
  } else { // all elements valid
    PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
  }

  outmap = new PairCorrelationMap;
  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 = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
      LOG(3, "Current atom is " << **iter << ".");
      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);
            for (std::vector<molecule *>::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules.end(); MolOtherWalker++){
                LOG(2, "Current other molecule is " << *MolOtherWalker << ".");
                for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
                  LOG(3, "Current otheratom is " << **runner << ".");
                  if ((*iter)->getId() < (*runner)->getId()){
                    for (set <pair<const element *,const element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
                      if ((PairRunner->first == (**iter).getType()) && (PairRunner->second == (**runner).getType())) {
                        periodicOtherX = FullInverseMatrix * ((*runner)->getPosition()); // x now in [0,1)^3
                        // go through every range in xyz and get distance
                        for (Othern[0]=-ranges[0]; Othern[0] <= ranges[0]; Othern[0]++)
                          for (Othern[1]=-ranges[1]; Othern[1] <= ranges[1]; Othern[1]++)
                            for (Othern[2]=-ranges[2]; Othern[2] <= ranges[2]; Othern[2]++) {
                              checkOtherX = FullMatrix * (Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX);
                              distance = checkX.distance(checkOtherX);
                              //LOG(1, "Inserting " << *(*iter) << " and " << *(*runner));
                              outmap->insert (
                                  std::pair<double, std::pair <const atom *, const atom*> > (
                                      distance,
                                      std::pair<const atom *, const atom*> (
                                          (*iter),
                                          (*runner))
                                      )
                                  );
                            }
                      }
                    }
                  }
                }
      }
    }
  }

  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++)
    (*MolWalker)->doCountAtoms();
  outmap = new CorrelationToPointMap;
  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
    LOG(2, "Current molecule is " << *MolWalker << ".");
    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*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++)
    (*MolWalker)->doCountAtoms();
  outmap = new CorrelationToPointMap;
  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 = (*MolWalker)->begin(); iter != (*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++)
    (*MolWalker)->doCountAtoms();
  outmap = new CorrelationToSurfaceMap;
  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 = (*MolWalker)->begin(); iter != (*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;
  }
  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
    (*MolWalker)->doCountAtoms();
  outmap = new CorrelationToSurfaceMap;
  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 = (*MolWalker)->begin(); iter != (*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);
};