Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/analysis_correlation.cpp

    r112b09 re65de8  
    99
    1010#include <iostream>
     11#include <iomanip>
    1112
    1213#include "analysis_correlation.hpp"
    … …  
    1920#include "triangleintersectionlist.hpp"
    2021#include "vector.hpp"
     22#include "Matrix.hpp"
    2123#include "verbose.hpp"
    2224#include "World.hpp"
     25#include "Box.hpp"
    2326
    2427
    2528/** Calculates the pair correlation between given elements.
    2629 * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
    27  * \param *out output stream for debugging
    28  * \param *molecules list of molecules structure
    29  * \param *type1 first element or NULL (if any element)
    30  * \param *type2 second element or NULL (if any element)
     30 * \param *molecules vector of molecules
     31 * \param &elements vector of elements to correlate
    3132 * \return Map of doubles with values the pair of the two atoms.
    3233 */
    33 PairCorrelationMap *PairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const type2 )
     34PairCorrelationMap *PairCorrelation(std::vector<molecule *> &molecules, const std::vector<element *> &elements)
    3435{
    3536  Info FunctionInfo(__func__);
    3637  PairCorrelationMap *outmap = NULL;
    3738  double distance = 0.;
    38 
    39   if (molecules->ListOfMolecules.empty()) {
     39  Box &domain = World::getInstance().getDomain();
     40
     41  if (molecules.empty()) {
    4042    DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
    4143    return outmap;
    4244  }
    43   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    44     (*MolWalker)->doCountAtoms();
     45  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
     46    (*MolWalker)->doCountAtoms();
     47
     48  // create all possible pairs of elements
     49  set <pair<element *, element *> > PairsOfElements;
     50  if (elements.size() >= 2) {
     51    for (vector<element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
     52      for (vector<element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
     53        if (type1 != type2) {
     54          PairsOfElements.insert( pair<element *, element*>(*type1,*type2) );
     55          DoLog(1) && (Log() << Verbose(1) << "Creating element pair " << (*type1)->symbol << " and " << (*type2)->symbol << "." << endl);
     56        }
     57  } else if (elements.size() == 1) { // one to all are valid
     58    element *elemental = *elements.begin();
     59    PairsOfElements.insert( pair<element *, element*>(elemental,(element *)NULL) );
     60    PairsOfElements.insert( pair<element *, element*>((element *)NULL,elemental) );
     61  } else { // all elements valid
     62    PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
     63  }
     64
    4565  outmap = new PairCorrelationMap;
    46   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
    47     if ((*MolWalker)->ActiveFlag) {
    48       DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    49       eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
    50       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    51         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
    52         if ((type1 == NULL) || ((*iter)->type == type1)) {
    53           for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++){
    54             if ((*MolOtherWalker)->ActiveFlag) {
    55               DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
    56               for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
    57                 DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
    58                 if ((*iter)->getId() < (*runner)->getId()){
    59                   if ((type2 == NULL) || ((*runner)->type == type2)) {
    60                     distance = (*iter)->node->PeriodicDistance(*(*runner)->node,  World::getInstance().getDomain());
    61                     //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
    62                     outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
    63                   }
    64                 }
     66  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++){
     67    DoLog(2) && (Log()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
     68    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
     69      DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
     70      for (std::vector<molecule *>::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules.end(); MolOtherWalker++){
     71        DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
     72        for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
     73          DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
     74          if ((*iter)->getId() < (*runner)->getId()){
     75            for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
     76              if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
     77                distance = domain.periodicDistance(*(*iter)->node,*(*runner)->node);
     78                //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
     79                outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
    6580              }
    66             }
    6781          }
    6882        }
    … …  
    7589/** Calculates the pair correlation between given elements.
    7690 * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
    77  * \param *out output stream for debugging
    7891 * \param *molecules list of molecules structure
    79  * \param *type1 first element or NULL (if any element)
    80  * \param *type2 second element or NULL (if any element)
     92 * \param &elements vector of elements to correlate
    8193 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
    8294 * \return Map of doubles with values the pair of the two atoms.
    8395 */
    84 PairCorrelationMap *PeriodicPairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const type2, const int ranges[NDIM] )
     96PairCorrelationMap *PeriodicPairCorrelation(std::vector<molecule *> &molecules, const std::vector<element *> &elements, const int ranges[NDIM] )
    8597{
    8698  Info FunctionInfo(__func__);
    … …  
    94106  Vector periodicOtherX;
    95107
    96   if (molecules->ListOfMolecules.empty()) {
     108  if (molecules.empty()) {
    97109    DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
    98110    return outmap;
    99111  }
    100   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    101     (*MolWalker)->doCountAtoms();
     112  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
     113    (*MolWalker)->doCountAtoms();
     114
     115  // create all possible pairs of elements
     116  set <pair<element *, element *> > PairsOfElements;
     117  if (elements.size() >= 2) {
     118    for (vector<element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
     119      for (vector<element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
     120        if (type1 != type2) {
     121          PairsOfElements.insert( pair<element *, element*>(*type1,*type2) );
     122          DoLog(1) && (Log() << Verbose(1) << "Creating element pair " << (*type1)->symbol << " and " << (*type2)->symbol << "." << endl);
     123        }
     124  } else if (elements.size() == 1) { // one to all are valid
     125    element *elemental = *elements.begin();
     126    PairsOfElements.insert( pair<element *, element*>(elemental,(element *)NULL) );
     127    PairsOfElements.insert( pair<element *, element*>((element *)NULL,elemental) );
     128  } else { // all elements valid
     129    PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
     130  }
     131
    102132  outmap = new PairCorrelationMap;
    103   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    104     if ((*MolWalker)->ActiveFlag) {
    105       double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
    106       double * FullInverseMatrix = InverseMatrix(FullMatrix);
    107       DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    108       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    109         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
    110         if ((type1 == NULL) || ((*iter)->type == type1)) {
    111           periodicX = *(*iter)->node;
    112           periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    113           // go through every range in xyz and get distance
    114           for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
    115             for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
    116               for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
    117                 checkX = Vector(n[0], n[1], n[2]) + periodicX;
    118                 checkX.MatrixMultiplication(FullMatrix);
    119                 for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
    120                   if ((*MolOtherWalker)->ActiveFlag) {
    121                     DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
    122                     for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
    123                       DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
    124                       if ((*iter)->nr < (*runner)->nr)
    125                         if ((type2 == NULL) || ((*runner)->type == type2)) {
    126                           periodicOtherX = *(*runner)->node;
    127                           periodicOtherX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    128                           // go through every range in xyz and get distance
    129                           for (Othern[0]=-ranges[0]; Othern[0] <= ranges[0]; Othern[0]++)
    130                             for (Othern[1]=-ranges[1]; Othern[1] <= ranges[1]; Othern[1]++)
    131                               for (Othern[2]=-ranges[2]; Othern[2] <= ranges[2]; Othern[2]++) {
    132                                 checkOtherX = Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX;
    133                                 checkOtherX.MatrixMultiplication(FullMatrix);
    134                                 distance = checkX.distance(checkOtherX);
    135                                 //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
    136                                 outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
    137                               }
    138                         }
     133  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++){
     134    Matrix FullMatrix = World::getInstance().getDomain().getM();
     135    Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
     136    DoLog(2) && (Log()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
     137    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
     138      DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
     139      periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
     140      // go through every range in xyz and get distance
     141      for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
     142        for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
     143          for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
     144            checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
     145            for (std::vector<molecule *>::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules.end(); MolOtherWalker++){
     146                DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
     147                for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
     148                  DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
     149                  if ((*iter)->getId() < (*runner)->getId()){
     150                    for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
     151                      if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
     152                        periodicOtherX = FullInverseMatrix * (*(*runner)->node); // x now in [0,1)^3
     153                        // go through every range in xyz and get distance
     154                        for (Othern[0]=-ranges[0]; Othern[0] <= ranges[0]; Othern[0]++)
     155                          for (Othern[1]=-ranges[1]; Othern[1] <= ranges[1]; Othern[1]++)
     156                            for (Othern[2]=-ranges[2]; Othern[2] <= ranges[2]; Othern[2]++) {
     157                              checkOtherX = FullMatrix * (Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX);
     158                              distance = checkX.distance(checkOtherX);
     159                              //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
     160                              outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
     161                            }
     162                      }
     163                    }
    139164                  }
    140               }
    141             }
    142         }
     165                }
    143166      }
    144       delete[](FullMatrix);
    145       delete[](FullInverseMatrix);
    146     }
     167    }
     168  }
    147169
    148170  return outmap;
    … …  
    150172
    151173/** Calculates the distance (pair) correlation between a given element and a point.
    152  * \param *out output stream for debugging
    153174 * \param *molecules list of molecules structure
    154  * \param *type element or NULL (if any element)
     175 * \param &elements vector of elements to correlate with point
    155176 * \param *point vector to the correlation point
    156177 * \return Map of dobules with values as pairs of atom and the vector
    157178 */
    158 CorrelationToPointMap *CorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point )
     179CorrelationToPointMap *CorrelationToPoint(std::vector<molecule *> &molecules, const std::vector<element *> &elements, const Vector *point )
    159180{
    160181  Info FunctionInfo(__func__);
    161182  CorrelationToPointMap *outmap = NULL;
    162183  double distance = 0.;
    163 
    164   if (molecules->ListOfMolecules.empty()) {
     184  Box &domain = World::getInstance().getDomain();
     185
     186  if (molecules.empty()) {
    165187    DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
    166188    return outmap;
    167189  }
    168   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     190  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
    169191    (*MolWalker)->doCountAtoms();
    170192  outmap = new CorrelationToPointMap;
    171   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    172     if ((*MolWalker)->ActiveFlag) {
    173       DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    174       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    175         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
    176         if ((type == NULL) || ((*iter)->type == type)) {
    177           distance = (*iter)->node->PeriodicDistance(*point, World::getInstance().getDomain());
     193  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
     194    DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
     195    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
     196      DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
     197      for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
     198        if ((*type == NULL) || ((*iter)->type == *type)) {
     199          distance = domain.periodicDistance(*(*iter)->node,*point);
    178200          DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
    179201          outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
    180202        }
    181       }
    182     }
     203    }
     204  }
    183205
    184206  return outmap;
    … …  
    186208
    187209/** Calculates the distance (pair) correlation between a given element, all its periodic images and a point.
    188  * \param *out output stream for debugging
    189210 * \param *molecules list of molecules structure
    190  * \param *type element or NULL (if any element)
     211 * \param &elements vector of elements to correlate to point
    191212 * \param *point vector to the correlation point
    192213 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
    193214 * \return Map of dobules with values as pairs of atom and the vector
    194215 */
    195 CorrelationToPointMap *PeriodicCorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point, const int ranges[NDIM] )
     216CorrelationToPointMap *PeriodicCorrelationToPoint(std::vector<molecule *> &molecules, const std::vector<element *> &elements, const Vector *point, const int ranges[NDIM] )
    196217{
    197218  Info FunctionInfo(__func__);
    … …  
    202223  Vector checkX;
    203224
    204   if (molecules->ListOfMolecules.empty()) {
     225  if (molecules.empty()) {
    205226    DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
    206227    return outmap;
    207228  }
    208   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     229  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
    209230    (*MolWalker)->doCountAtoms();
    210231  outmap = new CorrelationToPointMap;
    211   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    212     if ((*MolWalker)->ActiveFlag) {
    213       double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
    214       double * FullInverseMatrix = InverseMatrix(FullMatrix);
    215       DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    216       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    217         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
    218         if ((type == NULL) || ((*iter)->type == type)) {
    219           periodicX = *(*iter)->node;
    220           periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
     232  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
     233    Matrix FullMatrix = World::getInstance().getDomain().getM();
     234    Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
     235    DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
     236    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
     237      DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
     238      for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
     239        if ((*type == NULL) || ((*iter)->type == *type)) {
     240          periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
    221241          // go through every range in xyz and get distance
    222242          for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
    223243            for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
    224244              for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
    225                 checkX = Vector(n[0], n[1], n[2]) + periodicX;
    226                 checkX.MatrixMultiplication(FullMatrix);
     245                checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
    227246                distance = checkX.distance(*point);
    228247                DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
    … …  
    230249              }
    231250        }
    232       }
    233       delete[](FullMatrix);
    234       delete[](FullInverseMatrix);
    235     }
     251    }
     252  }
    236253
    237254  return outmap;
    … …  
    239256
    240257/** Calculates the distance (pair) correlation between a given element and a surface.
    241  * \param *out output stream for debugging
    242258 * \param *molecules list of molecules structure
    243  * \param *type element or NULL (if any element)
     259 * \param &elements vector of elements to correlate to surface
    244260 * \param *Surface pointer to Tesselation class surface
    245261 * \param *LC LinkedCell structure to quickly find neighbouring atoms
    246262 * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
    247263 */
    248 CorrelationToSurfaceMap *CorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC )
     264CorrelationToSurfaceMap *CorrelationToSurface(std::vector<molecule *> &molecules, const std::vector<element *> &elements, const Tesselation * const Surface, const LinkedCell *LC )
    249265{
    250266  Info FunctionInfo(__func__);
    … …  
    254270  Vector centroid;
    255271
    256   if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
     272  if ((Surface == NULL) || (LC == NULL) || (molecules.empty())) {
    257273    DoeLog(1) && (eLog()<< Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
    258274    return outmap;
    259275  }
    260   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     276  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
    261277    (*MolWalker)->doCountAtoms();
    262278  outmap = new CorrelationToSurfaceMap;
    263   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    264     if ((*MolWalker)->ActiveFlag) {
    265       DoLog(1) && (Log() << Verbose(1) << "Current molecule is " << (*MolWalker)->name << "." << endl);
    266       if ((*MolWalker)->empty())
    267         DoLog(1) && (1) && (Log() << Verbose(1) << "\t is empty." << endl);
    268       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    269         DoLog(1) && (Log() << Verbose(1) << "\tCurrent atom is " << *(*iter) << "." << endl);
    270         if ((type == NULL) || ((*iter)->type == type)) {
     279  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
     280    DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << (*MolWalker)->name << "." << endl);
     281    if ((*MolWalker)->empty())
     282      DoLog(2) && (2) && (Log() << Verbose(2) << "\t is empty." << endl);
     283    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
     284      DoLog(3) && (Log() << Verbose(3) << "\tCurrent atom is " << *(*iter) << "." << endl);
     285      for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
     286        if ((*type == NULL) || ((*iter)->type == *type)) {
    271287          TriangleIntersectionList Intersections((*iter)->node,Surface,LC);
    272288          distance = Intersections.GetSmallestDistance();
    … …  
    274290          outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> ((*iter), triangle) ) );
    275291        }
    276       }
    277     } else {
    278       DoLog(1) && (Log() << Verbose(1) << "molecule " << (*MolWalker)->name << " is not active." << endl);
    279     }
     292    }
     293  }
    280294
    281295  return outmap;
    … …  
    287301 * axis an integer from [ -ranges[i], ranges[i] ] onto it and multiply with the domain matrix to bring it back into
    288302 * the real space. Then, we Tesselation::FindClosestTriangleToPoint() and DistanceToTrianglePlane().
    289  * \param *out output stream for debugging
    290303 * \param *molecules list of molecules structure
    291  * \param *type element or NULL (if any element)
     304 * \param &elements vector of elements to correlate to surface
    292305 * \param *Surface pointer to Tesselation class surface
    293306 * \param *LC LinkedCell structure to quickly find neighbouring atoms
    … …  
    295308 * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
    296309 */
    297 CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
     310CorrelationToSurfaceMap *PeriodicCorrelationToSurface(std::vector<molecule *> &molecules, const std::vector<element *> &elements, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
    298311{
    299312  Info FunctionInfo(__func__);
    … …  
    306319  Vector checkX;
    307320
    308   if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
     321  if ((Surface == NULL) || (LC == NULL) || (molecules.empty())) {
    309322    DoLog(1) && (Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
    310323    return outmap;
    311324  }
    312   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     325  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
    313326    (*MolWalker)->doCountAtoms();
    314327  outmap = new CorrelationToSurfaceMap;
    315328  double ShortestDistance = 0.;
    316329  BoundaryTriangleSet *ShortestTriangle = NULL;
    317   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    318     if ((*MolWalker)->ActiveFlag) {
    319       double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
    320       double * FullInverseMatrix = InverseMatrix(FullMatrix);
    321       DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    322       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    323         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
    324         if ((type == NULL) || ((*iter)->type == type)) {
    325           periodicX = *(*iter)->node;
    326           periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
     330  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++) {
     331    Matrix FullMatrix = World::getInstance().getDomain().getM();
     332    Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
     333    DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
     334    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
     335      DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
     336      for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
     337        if ((*type == NULL) || ((*iter)->type == *type)) {
     338          periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
    327339          // go through every range in xyz and get distance
    328340          ShortestDistance = -1.;
    … …  
    330342            for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
    331343              for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
    332                 checkX = Vector(n[0], n[1], n[2]) + periodicX;
    333                 checkX.MatrixMultiplication(FullMatrix);
     344                checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
    334345                TriangleIntersectionList Intersections(&checkX,Surface,LC);
    335346                distance = Intersections.GetSmallestDistance();
    … …  
    344355          //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
    345356        }
    346       }
    347       delete[](FullMatrix);
    348       delete[](FullInverseMatrix);
    349     }
     357    }
     358  }
    350359
    351360  return outmap;
Note: See TracChangeset for help on using the changeset viewer.