Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/analysis_correlation.cpp

    ra67d19 rbd61b41  
    1515#include "tesselation.hpp"
    1616#include "tesselationhelpers.hpp"
    17 #include "triangleintersectionlist.hpp"
    1817#include "vector.hpp"
    1918#include "verbose.hpp"
    20 #include "World.hpp"
    2119
    2220
    … …  
    3634
    3735  if (molecules->ListOfMolecules.empty()) {
    38     DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
     36    eLog() << Verbose(1) <<"No molecule given." << endl;
    3937    return outmap;
    4038  }
    … …  
    4240  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    4341    if ((*MolWalker)->ActiveFlag) {
    44       DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    45       atom *Walker = (*MolWalker)->start;
    46       while (Walker->next != (*MolWalker)->end) {
    47         Walker = Walker->next;
    48         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
     42      eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
     43      atom *Walker = (*MolWalker)->start;
     44      while (Walker->next != (*MolWalker)->end) {
     45        Walker = Walker->next;
     46        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
    4947        if ((type1 == NULL) || (Walker->type == type1)) {
    5048          for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
    5149            if ((*MolOtherWalker)->ActiveFlag) {
    52               DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
     50              Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl;
    5351              atom *OtherWalker = (*MolOtherWalker)->start;
    5452              while (OtherWalker->next != (*MolOtherWalker)->end) { // only go up to Walker
    5553                OtherWalker = OtherWalker->next;
    56                 DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl);
     54                Log() << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl;
    5755                if (Walker->nr < OtherWalker->nr)
    5856                  if ((type2 == NULL) || (OtherWalker->type == type2)) {
    59                     distance = Walker->node->PeriodicDistance(OtherWalker->node, World::get()->cell_size);
     57                    distance = Walker->node->PeriodicDistance(OtherWalker->node, (*MolWalker)->cell_size);
    6058                    //Log() << Verbose(1) <<"Inserting " << *Walker << " and " << *OtherWalker << endl;
    6159                    outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> (Walker, OtherWalker) ) );
    … …  
    9290
    9391  if (molecules->ListOfMolecules.empty()) {
    94     DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
     92    eLog() << Verbose(1) <<"No molecule given." << endl;
    9593    return outmap;
    9694  }
    … …  
    9896  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    9997    if ((*MolWalker)->ActiveFlag) {
    100       double * FullMatrix = ReturnFullMatrixforSymmetric(World::get()->cell_size);
     98      double * FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
    10199      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    102       DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    103       atom *Walker = (*MolWalker)->start;
    104       while (Walker->next != (*MolWalker)->end) {
    105         Walker = Walker->next;
    106         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
     100      eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
     101      atom *Walker = (*MolWalker)->start;
     102      while (Walker->next != (*MolWalker)->end) {
     103        Walker = Walker->next;
     104        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
    107105        if ((type1 == NULL) || (Walker->type == type1)) {
    108106          periodicX.CopyVector(Walker->node);
    … …  
    117115                for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
    118116                  if ((*MolOtherWalker)->ActiveFlag) {
    119                     DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
     117                    Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl;
    120118                    atom *OtherWalker = (*MolOtherWalker)->start;
    121119                    while (OtherWalker->next != (*MolOtherWalker)->end) { // only go up to Walker
    122120                      OtherWalker = OtherWalker->next;
    123                       DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl);
     121                      Log() << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl;
    124122                      if (Walker->nr < OtherWalker->nr)
    125123                        if ((type2 == NULL) || (OtherWalker->type == type2)) {
    … …  
    164162
    165163  if (molecules->ListOfMolecules.empty()) {
    166     DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
     164    Log() << Verbose(1) <<"No molecule given." << endl;
    167165    return outmap;
    168166  }
    … …  
    170168  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    171169    if ((*MolWalker)->ActiveFlag) {
    172       DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    173       atom *Walker = (*MolWalker)->start;
    174       while (Walker->next != (*MolWalker)->end) {
    175         Walker = Walker->next;
    176         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
     170      Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
     171      atom *Walker = (*MolWalker)->start;
     172      while (Walker->next != (*MolWalker)->end) {
     173        Walker = Walker->next;
     174        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
    177175        if ((type == NULL) || (Walker->type == type)) {
    178           distance = Walker->node->PeriodicDistance(point, World::get()->cell_size);
    179           DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
     176          distance = Walker->node->PeriodicDistance(point, (*MolWalker)->cell_size);
     177          Log() << Verbose(4) << "Current distance is " << distance << "." << endl;
    180178          outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
    181179        }
    … …  
    204202
    205203  if (molecules->ListOfMolecules.empty()) {
    206     DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
     204    Log() << Verbose(1) <<"No molecule given." << endl;
    207205    return outmap;
    208206  }
    … …  
    210208  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    211209    if ((*MolWalker)->ActiveFlag) {
    212       double * FullMatrix = ReturnFullMatrixforSymmetric(World::get()->cell_size);
     210      double * FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
    213211      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    214       DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    215       atom *Walker = (*MolWalker)->start;
    216       while (Walker->next != (*MolWalker)->end) {
    217         Walker = Walker->next;
    218         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
     212      Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
     213      atom *Walker = (*MolWalker)->start;
     214      while (Walker->next != (*MolWalker)->end) {
     215        Walker = Walker->next;
     216        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
    219217        if ((type == NULL) || (Walker->type == type)) {
    220218          periodicX.CopyVector(Walker->node);
    … …  
    228226                checkX.MatrixMultiplication(FullMatrix);
    229227                distance = checkX.Distance(point);
    230                 DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
     228                Log() << Verbose(4) << "Current distance is " << distance << "." << endl;
    231229                outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
    232230              }
    … …  
    257255
    258256  if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
    259     DoeLog(1) && (eLog()<< Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
     257    Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl;
    260258    return outmap;
    261259  }
    … …  
    263261  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    264262    if ((*MolWalker)->ActiveFlag) {
    265       DoLog(1) && (Log() << Verbose(1) << "Current molecule is " << (*MolWalker)->name << "." << endl);
    266       atom *Walker = (*MolWalker)->start;
    267       while (Walker->next != (*MolWalker)->end) {
    268         Walker = Walker->next;
    269         //Log() << Verbose(1) << "Current atom is " << *Walker << "." << endl;
     263      Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
     264      atom *Walker = (*MolWalker)->start;
     265      while (Walker->next != (*MolWalker)->end) {
     266        Walker = Walker->next;
     267        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
    270268        if ((type == NULL) || (Walker->type == type)) {
    271           TriangleIntersectionList Intersections(Walker->node,Surface,LC);
    272           distance = Intersections.GetSmallestDistance();
    273           triangle = Intersections.GetClosestTriangle();
    274           outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
    275         }
    276       }
    277     } else
    278       DoLog(1) && (Log() << Verbose(1) << "molecule " << (*MolWalker)->name << " is not active." << endl);
    279 
     269          triangle = Surface->FindClosestTriangleToPoint(Walker->node, LC );
     270          if (triangle != NULL) {
     271            distance = DistanceToTrianglePlane(Walker->node, triangle);
     272            outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
     273          }
     274        }
     275      }
     276    }
    280277
    281278  return outmap;
    … …  
    307304
    308305  if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
    309     DoLog(1) && (Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
     306    Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl;
    310307    return outmap;
    311308  }
    312309  outmap = new CorrelationToSurfaceMap;
    313   double ShortestDistance = 0.;
    314   BoundaryTriangleSet *ShortestTriangle = NULL;
    315   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    316     if ((*MolWalker)->ActiveFlag) {
    317       double * FullMatrix = ReturnFullMatrixforSymmetric(World::get()->cell_size);
     310  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     311    if ((*MolWalker)->ActiveFlag) {
     312      double * FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
    318313      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    319       DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    320       atom *Walker = (*MolWalker)->start;
    321       while (Walker->next != (*MolWalker)->end) {
    322         Walker = Walker->next;
    323         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
     314      Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
     315      atom *Walker = (*MolWalker)->start;
     316      while (Walker->next != (*MolWalker)->end) {
     317        Walker = Walker->next;
     318        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
    324319        if ((type == NULL) || (Walker->type == type)) {
    325320          periodicX.CopyVector(Walker->node);
    326321          periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    327322          // go through every range in xyz and get distance
    328           ShortestDistance = -1.;
    329323          for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
    330324            for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
    … …  
    333327                checkX.AddVector(&periodicX);
    334328                checkX.MatrixMultiplication(FullMatrix);
    335                 TriangleIntersectionList Intersections(&checkX,Surface,LC);
    336                 distance = Intersections.GetSmallestDistance();
    337                 triangle = Intersections.GetClosestTriangle();
    338                 if ((ShortestDistance == -1.) || (distance < ShortestDistance)) {
    339                   ShortestDistance = distance;
    340                   ShortestTriangle = triangle;
     329                triangle = Surface->FindClosestTriangleToPoint(&checkX, LC );
     330                if (triangle != NULL) {
     331                  distance = DistanceToTrianglePlane(&checkX, triangle);
     332                  outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
    341333                }
    342               }
    343           // insert
    344           outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (Walker, ShortestTriangle) ) );
    345           //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
     334          }
    346335        }
    347336      }
    … …  
    373362{
    374363  Info FunctionInfo(__func__);
    375   *file << "BinStart\tCount" << endl;
     364  *file << "# BinStart\tCount" << endl;
    376365  for (BinPairMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    377     *file << setprecision(8) << runner->first << "\t" << runner->second << endl;
     366    *file << runner->first << "\t" << runner->second << endl;
    378367  }
    379368};
    … …  
    386375{
    387376  Info FunctionInfo(__func__);
    388   *file << "BinStart\tAtom1\tAtom2" << endl;
     377  *file << "# BinStart\tAtom1\tAtom2" << endl;
    389378  for (PairCorrelationMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    390     *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
     379    *file << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
    391380  }
    392381};
    … …  
    399388{
    400389  Info FunctionInfo(__func__);
    401   *file << "BinStart\tAtom::x[i]-point.x[i]" << endl;
     390  *file << "# BinStart\tAtom::x[i]-point.x[i]" << endl;
    402391  for (CorrelationToPointMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    403392    *file << runner->first;
    404393    for (int i=0;i<NDIM;i++)
    405       *file << "\t" << setprecision(8) << (runner->second.first->node->x[i] - runner->second.second->x[i]);
     394      *file << "\t" << (runner->second.first->node->x[i] - runner->second.second->x[i]);
    406395    *file << endl;
    407396  }
    … …  
    415404{
    416405  Info FunctionInfo(__func__);
    417   *file << "BinStart\tTriangle" << endl;
    418   if (!map->empty())
    419     for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    420       *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
    421     }
    422 };
    423 
     406  *file << "# BinStart\tTriangle" << endl;
     407  for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
     408    *file << runner->first << "\t" << *(runner->second.second) << endl;
     409  }
     410};
     411
Note: See TracChangeset for help on using the changeset viewer.