source: src/analysis_correlation.cpp@ 7326b2

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

#include <iostream>

#include "analysis_correlation.hpp"
#include "element.hpp"
#include "molecule.hpp"
#include "tesselation.hpp"
#include "tesselationhelpers.hpp"
#include "vector.hpp"
#include "verbose.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( ofstream * const out, MoleculeListClass * const &molecules, const element * const type1, const element * const type2 )
{
  PairCorrelationMap *outmap = NULL;
  double distance = 0.;

  if (molecules->ListOfMolecules.empty()) {
    cerr << Verbose(1) <<"No molecule given." << endl;
    return outmap;
  }
  outmap = new PairCorrelationMap;
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    if ((*MolWalker)->ActiveFlag) {
      cerr << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
      atom *Walker = (*MolWalker)->start;
      while (Walker->next != (*MolWalker)->end) {
        Walker = Walker->next;
        *out << Verbose(3) << "Current atom is " << *Walker << "." << endl;
        if ((type1 == NULL) || (Walker->type == type1)) {
          for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
            if ((*MolOtherWalker)->ActiveFlag) {
              *out << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl;
              atom *OtherWalker = (*MolOtherWalker)->start;
              while (OtherWalker->next != (*MolOtherWalker)->end) { // only go up to Walker
                OtherWalker = OtherWalker->next;
                *out << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl;
                if (Walker->nr < OtherWalker->nr)
                  if ((type2 == NULL) || (OtherWalker->type == type2)) {
                    distance = Walker->node->PeriodicDistance(OtherWalker->node, (*MolWalker)->cell_size);
                    //*out << Verbose(1) <<"Inserting " << *Walker << " and " << *OtherWalker << endl;
                    outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> (Walker, OtherWalker) ) );
                  }
              }
          }
        }
      }
    }

  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( ofstream * const out, MoleculeListClass * const &molecules, const element * const type1, const element * const type2, const int ranges[NDIM] )
{
  PairCorrelationMap *outmap = NULL;
  double distance = 0.;
  int n[NDIM];
  Vector checkX;
  Vector periodicX;
  int Othern[NDIM];
  Vector checkOtherX;
  Vector periodicOtherX;

  if (molecules->ListOfMolecules.empty()) {
    cerr << Verbose(1) <<"No molecule given." << endl;
    return outmap;
  }
  outmap = new PairCorrelationMap;
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    if ((*MolWalker)->ActiveFlag) {
      const double * const FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
      const double * const FullInverseMatrix = InverseMatrix(FullMatrix);
      cerr << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
      atom *Walker = (*MolWalker)->start;
      while (Walker->next != (*MolWalker)->end) {
        Walker = Walker->next;
        *out << Verbose(3) << "Current atom is " << *Walker << "." << endl;
        if ((type1 == NULL) || (Walker->type == type1)) {
          periodicX.CopyVector(Walker->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.Init(n[0], n[1], n[2]);
                checkX.AddVector(&periodicX);
                checkX.MatrixMultiplication(FullMatrix);
                for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
                  if ((*MolOtherWalker)->ActiveFlag) {
                    *out << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl;
                    atom *OtherWalker = (*MolOtherWalker)->start;
                    while (OtherWalker->next != (*MolOtherWalker)->end) { // only go up to Walker
                      OtherWalker = OtherWalker->next;
                      *out << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl;
                      if (Walker->nr < OtherWalker->nr)
                        if ((type2 == NULL) || (OtherWalker->type == type2)) {
                          periodicOtherX.CopyVector(OtherWalker->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.Init(Othern[0], Othern[1], Othern[2]);
                                checkOtherX.AddVector(&periodicOtherX);
                                checkOtherX.MatrixMultiplication(FullMatrix);
                                distance = checkX.Distance(&checkOtherX);
                                //*out << Verbose(1) <<"Inserting " << *Walker << " and " << *OtherWalker << endl;
                                outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> (Walker, OtherWalker) ) );
                              }
                        }
                  }
              }
            }
        }
      }
    }

  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(  ofstream * const out, MoleculeListClass * const &molecules, const element * const type, const Vector *point )
{
  CorrelationToPointMap *outmap = NULL;
  double distance = 0.;

  if (molecules->ListOfMolecules.empty()) {
    *out << Verbose(1) <<"No molecule given." << endl;
    return outmap;
  }
  outmap = new CorrelationToPointMap;
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    if ((*MolWalker)->ActiveFlag) {
      *out << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
      atom *Walker = (*MolWalker)->start;
      while (Walker->next != (*MolWalker)->end) {
        Walker = Walker->next;
        *out << Verbose(3) << "Current atom is " << *Walker << "." << endl;
        if ((type == NULL) || (Walker->type == type)) {
          distance = Walker->node->PeriodicDistance(point, (*MolWalker)->cell_size);
          *out << Verbose(4) << "Current distance is " << distance << "." << endl;
          outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, 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(  ofstream * const out, MoleculeListClass * const &molecules, const element * const type, const Vector *point, const int ranges[NDIM] )
{
  CorrelationToPointMap *outmap = NULL;
  double distance = 0.;
  int n[NDIM];
  Vector periodicX;
  Vector checkX;

  if (molecules->ListOfMolecules.empty()) {
    *out << Verbose(1) <<"No molecule given." << endl;
    return outmap;
  }
  outmap = new CorrelationToPointMap;
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    if ((*MolWalker)->ActiveFlag) {
      const double * const FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
      const double * const FullInverseMatrix = InverseMatrix(FullMatrix);
      *out << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
      atom *Walker = (*MolWalker)->start;
      while (Walker->next != (*MolWalker)->end) {
        Walker = Walker->next;
        *out << Verbose(3) << "Current atom is " << *Walker << "." << endl;
        if ((type == NULL) || (Walker->type == type)) {
          periodicX.CopyVector(Walker->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.Init(n[0], n[1], n[2]);
                checkX.AddVector(&periodicX);
                checkX.MatrixMultiplication(FullMatrix);
                distance = checkX.Distance(point);
                *out << Verbose(4) << "Current distance is " << distance << "." << endl;
                outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
              }
        }
      }
    }

  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( ofstream * const out, MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC )
{
  CorrelationToSurfaceMap *outmap = NULL;
  double distance = 0;
  class BoundaryTriangleSet *triangle = NULL;
  Vector centroid;

  if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
    *out << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl;
    return outmap;
  }
  outmap = new CorrelationToSurfaceMap;
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    if ((*MolWalker)->ActiveFlag) {
      *out << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
      atom *Walker = (*MolWalker)->start;
      while (Walker->next != (*MolWalker)->end) {
        Walker = Walker->next;
        *out << Verbose(3) << "Current atom is " << *Walker << "." << endl;
        if ((type == NULL) || (Walker->type == type)) {
          triangle = Surface->FindClosestTriangleToPoint(out, Walker->node, LC );
          if (triangle != NULL) {
            distance = DistanceToTrianglePlane(out, Walker->node, triangle);
            outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, 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 *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( ofstream * const out, MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
{
  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())) {
    *out << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl;
    return outmap;
  }
  outmap = new CorrelationToSurfaceMap;
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    if ((*MolWalker)->ActiveFlag) {
      const double * const FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
      const double * const FullInverseMatrix = InverseMatrix(FullMatrix);
      *out << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
      atom *Walker = (*MolWalker)->start;
      while (Walker->next != (*MolWalker)->end) {
        Walker = Walker->next;
        *out << Verbose(3) << "Current atom is " << *Walker << "." << endl;
        if ((type == NULL) || (Walker->type == type)) {
          periodicX.CopyVector(Walker->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.Init(n[0], n[1], n[2]);
                checkX.AddVector(&periodicX);
                checkX.MatrixMultiplication(FullMatrix);
                triangle = Surface->FindClosestTriangleToPoint(out, &checkX, LC );
                if (triangle != NULL) {
                  distance = DistanceToTrianglePlane(out, &checkX, triangle);
                  outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
                }
          }
        }
      }
    }

  return outmap;
};

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


/** 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 )
{
  *file << "# BinStart\tCount" << endl;
  for (BinPairMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    *file << 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 )
{
  *file << "# BinStart\tAtom1\tAtom2" << endl;
  for (PairCorrelationMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    *file << 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 )
{
  *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" << (runner->second.first->node->x[i] - runner->second.second->x[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 )
{
  *file << "# BinStart\tTriangle" << endl;
  for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    *file << runner->first << "\t" << *(runner->second.second) << endl;
  }
};