source: src/analysis_correlation.cpp@ 7cb450

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010 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 "Helpers/MemDebug.hpp"

#include <iostream>
#include <iomanip>

#include "BoundaryTriangleSet.hpp"
#include "analysis_correlation.hpp"
#include "element.hpp"
#include "Helpers/Info.hpp"
#include "Helpers/Log.hpp"
#include "molecule.hpp"
#include "tesselation.hpp"
#include "tesselationhelpers.hpp"
#include "triangleintersectionlist.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "LinearAlgebra/Matrix.hpp"
#include "Helpers/Verbose.hpp"
#include "World.hpp"
#include "Box.hpp"


/** 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 = NULL;
  double distance = 0.;
  Box &domain = World::getInstance().getDomain();

  if (molecules.empty()) {
    DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
    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) );
          DoLog(1) && (Log() << Verbose(1) << "Creating element pair " << *(*type1) << " and " << *(*type2) << "." << endl);
        }
  } 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++){
    DoLog(2) && (Log()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
      DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
      for (std::vector<molecule *>::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules.end(); MolOtherWalker++){
        DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
        for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
          DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
          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() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
                outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, 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 = NULL;
  double distance = 0.;
  int n[NDIM];
  Vector checkX;
  Vector periodicX;
  int Othern[NDIM];
  Vector checkOtherX;
  Vector periodicOtherX;

  if (molecules.empty()) {
    DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
    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) );
          DoLog(1) && (Log() << Verbose(1) << "Creating element pair " << *(*type1) << " and " << *(*type2) << "." << endl);
        }
  } 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++){
    Matrix FullMatrix = World::getInstance().getDomain().getM();
    Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
    DoLog(2) && (Log()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
      DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
      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++){
                DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
                for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
                  DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
                  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() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
                              outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, 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 = NULL;
  double distance = 0.;
  Box &domain = World::getInstance().getDomain();

  if (molecules.empty()) {
    DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
    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++) {
    DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
      DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
      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);
          DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
          outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<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 = NULL;
  double distance = 0.;
  int n[NDIM];
  Vector periodicX;
  Vector checkX;

  if (molecules.empty()) {
    DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
    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++) {
    Matrix FullMatrix = World::getInstance().getDomain().getM();
    Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
    DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
      DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
      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);
                DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
                outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<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 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 *LC )
{
  Info FunctionInfo(__func__);
  CorrelationToSurfaceMap *outmap = NULL;
  double distance = 0;
  class BoundaryTriangleSet *triangle = NULL;
  Vector centroid;

  if ((Surface == NULL) || (LC == NULL) || (molecules.empty())) {
    DoeLog(1) && (eLog()<< Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
    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++) {
    DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << (*MolWalker)->name << "." << endl);
    if ((*MolWalker)->empty())
      DoLog(2) && (2) && (Log() << Verbose(2) << "\t is empty." << endl);
    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
      DoLog(3) && (Log() << Verbose(3) << "\tCurrent atom is " << *(*iter) << "." << endl);
      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 ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<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 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 *LC, const int ranges[NDIM] )
{
  Info FunctionInfo(__func__);
  CorrelationToSurfaceMap *outmap = NULL;
  double distance = 0;
  class BoundaryTriangleSet *triangle = NULL;
  Vector centroid;
  int n[NDIM];
  Vector periodicX;
  Vector checkX;

  if ((Surface == NULL) || (LC == NULL) || (molecules.empty())) {
    DoLog(1) && (Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
    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++) {
    Matrix FullMatrix = World::getInstance().getDomain().getM();
    Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
    DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
      DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
      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 ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (*iter, ShortestTriangle) ) );
          //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
        }
    }
  }

  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);
};


/** Prints correlation (double, int) pairs to file.
 * \param *file file to write to
 * \param *map map to write
 */
void OutputCorrelation( ofstream * const file, const BinPairMap * const map )
{
  Info FunctionInfo(__func__);
  *file << "BinStart\tCount" << endl;
  for (BinPairMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    *file << setprecision(8) << runner->first << "\t" << runner->second << endl;
  }
};

/** Prints correlation (double, (atom*,atom*) ) pairs to file.
 * \param *file file to write to
 * \param *map map to write
 */
void OutputPairCorrelation( ofstream * const file, const PairCorrelationMap * const map )
{
  Info FunctionInfo(__func__);
  *file << "BinStart\tAtom1\tAtom2" << endl;
  for (PairCorrelationMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
  }
};

/** Prints correlation (double, int) pairs to file.
 * \param *file file to write to
 * \param *map map to write
 */
void OutputCorrelationToPoint( ofstream * const file, const CorrelationToPointMap * const map )
{
  Info FunctionInfo(__func__);
  *file << "BinStart\tAtom::x[i]-point.x[i]" << endl;
  for (CorrelationToPointMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    *file << runner->first;
    for (int i=0;i<NDIM;i++)
      *file << "\t" << setprecision(8) << (runner->second.first->at(i) - runner->second.second->at(i));
    *file << endl;
  }
};

/** Prints correlation (double, int) pairs to file.
 * \param *file file to write to
 * \param *map map to write
 */
void OutputCorrelationToSurface( ofstream * const file, const CorrelationToSurfaceMap * const map )
{
  Info FunctionInfo(__func__);
  *file << "BinStart\tTriangle" << endl;
  if (!map->empty())
    for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
      *file << setprecision(8) << runner->first << "\t";
      *file << *(runner->second.first) << "\t";
      *file << *(runner->second.second) << endl;
    }
};