source: src/analysis_correlation.cpp@ 98a8b4

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

#include <iostream>

#include "analysis_correlation.hpp"
#include "element.hpp"
#include "info.hpp"
#include "log.hpp"
#include "molecule.hpp"
#include "tesselation.hpp"
#include "tesselationhelpers.hpp"
#include "triangleintersectionlist.hpp"
#include "vector.hpp"
#include "verbose.hpp"
#include "World.hpp"


/** Calculates the pair correlation between given elements.
 * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
 * \param *out output stream for debugging
 * \param *molecules list of molecules structure
 * \param *type1 first element or NULL (if any element)
 * \param *type2 second element or NULL (if any element)
 * \return Map of doubles with values the pair of the two atoms.
 */
PairCorrelationMap *PairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const type2 )
{
  Info FunctionInfo(__func__);
  PairCorrelationMap *outmap = NULL;
  double distance = 0.;

  if (molecules->ListOfMolecules.empty()) {
    DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
    return outmap;
  }
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    (*MolWalker)->doCountAtoms();
  outmap = new PairCorrelationMap;
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
    if ((*MolWalker)->ActiveFlag) {
      DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
      eLog() << 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);
        if ((type1 == NULL) || ((*iter)->type == type1)) {
          for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++){
            if ((*MolOtherWalker)->ActiveFlag) {
              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()){
                  if ((type2 == NULL) || ((*runner)->type == type2)) {
                    distance = (*iter)->node->PeriodicDistance(*(*runner)->node,  World::getInstance().getDomain());
                    //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 *out output stream for debugging
 * \param *molecules list of molecules structure
 * \param *type1 first element or NULL (if any element)
 * \param *type2 second element or NULL (if any element)
 * \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(MoleculeListClass * const &molecules, const element * const type1, const element * const type2, 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->ListOfMolecules.empty()) {
    DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
    return outmap;
  }
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    (*MolWalker)->doCountAtoms();
  outmap = new PairCorrelationMap;
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    if ((*MolWalker)->ActiveFlag) {
      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
      double * FullInverseMatrix = InverseMatrix(FullMatrix);
      DoeLog(2) && (eLog()<< 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);
        if ((type1 == NULL) || ((*iter)->type == type1)) {
          periodicX = *(*iter)->node;
          periodicX.MatrixMultiplication(FullInverseMatrix);  // 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 = Vector(n[0], n[1], n[2]) + periodicX;
                checkX.MatrixMultiplication(FullMatrix);
                for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
                  if ((*MolOtherWalker)->ActiveFlag) {
                    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)->nr < (*runner)->nr)
                        if ((type2 == NULL) || ((*runner)->type == type2)) {
                          periodicOtherX = *(*runner)->node;
                          periodicOtherX.MatrixMultiplication(FullInverseMatrix);  // 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 = Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX;
                                checkOtherX.MatrixMultiplication(FullMatrix);
                                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)) ) );
                              }
                        }
                  }
              }
            }
        }
      }
      delete[](FullMatrix);
      delete[](FullInverseMatrix);
    }

  return outmap;
};

/** Calculates the distance (pair) correlation between a given element and a point.
 * \param *out output stream for debugging
 * \param *molecules list of molecules structure
 * \param *type element or NULL (if any element)
 * \param *point vector to the correlation point
 * \return Map of dobules with values as pairs of atom and the vector
 */
CorrelationToPointMap *CorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point )
{
  Info FunctionInfo(__func__);
  CorrelationToPointMap *outmap = NULL;
  double distance = 0.;

  if (molecules->ListOfMolecules.empty()) {
    DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
    return outmap;
  }
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    (*MolWalker)->doCountAtoms();
  outmap = new CorrelationToPointMap;
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    if ((*MolWalker)->ActiveFlag) {
      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);
        if ((type == NULL) || ((*iter)->type == type)) {
          distance = (*iter)->node->PeriodicDistance(*point, World::getInstance().getDomain());
          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 *out output stream for debugging
 * \param *molecules list of molecules structure
 * \param *type element or NULL (if any element)
 * \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(MoleculeListClass * const &molecules, const element * const type, 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->ListOfMolecules.empty()) {
    DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
    return outmap;
  }
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    (*MolWalker)->doCountAtoms();
  outmap = new CorrelationToPointMap;
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    if ((*MolWalker)->ActiveFlag) {
      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
      double * FullInverseMatrix = InverseMatrix(FullMatrix);
      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);
        if ((type == NULL) || ((*iter)->type == type)) {
          periodicX = *(*iter)->node;
          periodicX.MatrixMultiplication(FullInverseMatrix);  // 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 = Vector(n[0], n[1], n[2]) + periodicX;
                checkX.MatrixMultiplication(FullMatrix);
                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) ) );
              }
        }
      }
      delete[](FullMatrix);
      delete[](FullInverseMatrix);
    }

  return outmap;
};

/** Calculates the distance (pair) correlation between a given element and a surface.
 * \param *out output stream for debugging
 * \param *molecules list of molecules structure
 * \param *type element or NULL (if any element)
 * \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(MoleculeListClass * const &molecules, const element * const type, 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->ListOfMolecules.empty())) {
    DoeLog(1) && (eLog()<< Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
    return outmap;
  }
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    (*MolWalker)->doCountAtoms();
  outmap = new CorrelationToSurfaceMap;
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    if ((*MolWalker)->ActiveFlag) {
      DoLog(1) && (Log() << Verbose(1) << "Current molecule is " << (*MolWalker)->name << "." << endl);
      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
        //Log() << Verbose(3) << "Current atom is " << *(*iter) << "." << endl;
        if ((type == NULL) || ((*iter)->type == type)) {
          TriangleIntersectionList Intersections((*iter)->node,Surface,LC);
          distance = Intersections.GetSmallestDistance();
          triangle = Intersections.GetClosestTriangle();
          outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> ((*iter), triangle) ) );
        }
      }
    } else
      DoLog(1) && (Log() << Verbose(1) << "molecule " << (*MolWalker)->name << " is not active." << endl);

  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 *out output stream for debugging
 * \param *molecules list of molecules structure
 * \param *type element or NULL (if any element)
 * \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(MoleculeListClass * const &molecules, const element * const type, 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->ListOfMolecules.empty())) {
    DoLog(1) && (Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
    return outmap;
  }
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    (*MolWalker)->doCountAtoms();
  outmap = new CorrelationToSurfaceMap;
  double ShortestDistance = 0.;
  BoundaryTriangleSet *ShortestTriangle = NULL;
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    if ((*MolWalker)->ActiveFlag) {
      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
      double * FullInverseMatrix = InverseMatrix(FullMatrix);
      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);
        if ((type == NULL) || ((*iter)->type == type)) {
          periodicX = *(*iter)->node;
          periodicX.MatrixMultiplication(FullInverseMatrix);  // 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 = Vector(n[0], n[1], n[2]) + periodicX;
                checkX.MatrixMultiplication(FullMatrix);
                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;
        }
      }
      delete[](FullMatrix);
      delete[](FullInverseMatrix);
    }

  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->node->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" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
    }
};