Ignore:
Timestamp:
Jun 16, 2010, 12:24:21 PM (16 years ago)
Author:
Tillmann Crueger <crueger@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, Candidate_v1.7.0, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
Children:
492279
Parents:
f8e486 (diff), 980dd6 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Broken: Merge commit 'Gitosis/stable' into stable

Conflicts:

molecuilder/src/Actions/AnalysisAction/PairCorrelationToPointAction.cpp
molecuilder/src/Actions/AnalysisAction/PairCorrelationToSurfaceAction.cpp
molecuilder/src/Makefile.am

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/analysis_correlation.cpp

    rf8e486 re6317b  
    2525/** Calculates the pair correlation between given elements.
    2626 * 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)
     27 * \param *molecules list of molecules structure
     28 * \param &elements vector of elements to correlate
    3129 * \return Map of doubles with values the pair of the two atoms.
    3230 */
    33 PairCorrelationMap *PairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const type2 )
     31PairCorrelationMap *PairCorrelation(MoleculeListClass * const &molecules, const std::vector<element *> &elements)
    3432{
    3533  Info FunctionInfo(__func__);
    3634  PairCorrelationMap *outmap = NULL;
    3735  double distance = 0.;
     36  double *domain = World::getInstance().getDomain();
    3837
    3938  if (molecules->ListOfMolecules.empty()) {
    … …  
    4342  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    4443    (*MolWalker)->doCountAtoms();
     44
     45  // create all possible pairs of elements
     46  set <pair<element *, element *> > PairsOfElements;
     47  if (elements.size() >= 2) {
     48    for (vector<element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
     49      for (vector<element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
     50        if (type1 != type2) {
     51          PairsOfElements.insert( pair<element *, element*>(*type1,*type2) );
     52          DoLog(1) && (Log() << Verbose(1) << "Creating element pair " << (*type1)->symbol << " and " << (*type2)->symbol << "." << endl);
     53        }
     54  } else if (elements.size() == 1) { // one to all are valid
     55    element *elemental = *elements.begin();
     56    PairsOfElements.insert( pair<element *, element*>(elemental,(element *)NULL) );
     57    PairsOfElements.insert( pair<element *, element*>((element *)NULL,elemental) );
     58  } else { // all elements valid
     59    PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
     60  }
     61
    4562  outmap = new PairCorrelationMap;
    4663  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
    … …  
    5067      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    5168        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());
     69        for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++){
     70          if ((*MolOtherWalker)->ActiveFlag) {
     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 = (*iter)->node->PeriodicDistance(*(*runner)->node,  domain);
    6178                    //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
    6279                    outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
    6380                  }
    64                 }
    6581              }
    6682            }
    … …  
    7591/** Calculates the pair correlation between given elements.
    7692 * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
    77  * \param *out output stream for debugging
    78  * \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)
     93 * \param *molecules list of molecules structure
     94 * \param &elements vector of elements to correlate
    8195 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
    8296 * \return Map of doubles with values the pair of the two atoms.
    8397 */
    84 PairCorrelationMap *PeriodicPairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const type2, const int ranges[NDIM] )
     98PairCorrelationMap *PeriodicPairCorrelation(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const int ranges[NDIM] )
    8599{
    86100  Info FunctionInfo(__func__);
    … …  
    100114  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    101115    (*MolWalker)->doCountAtoms();
     116
     117  // create all possible pairs of elements
     118  set <pair<element *, element *> > PairsOfElements;
     119  if (elements.size() >= 2) {
     120    for (vector<element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
     121      for (vector<element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
     122        if (type1 != type2) {
     123          PairsOfElements.insert( pair<element *, element*>(*type1,*type2) );
     124          DoLog(1) && (Log() << Verbose(1) << "Creating element pair " << (*type1)->symbol << " and " << (*type2)->symbol << "." << endl);
     125        }
     126  } else if (elements.size() == 1) { // one to all are valid
     127    element *elemental = *elements.begin();
     128    PairsOfElements.insert( pair<element *, element*>(elemental,(element *)NULL) );
     129    PairsOfElements.insert( pair<element *, element*>((element *)NULL,elemental) );
     130  } else { // all elements valid
     131    PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
     132  }
     133
    102134  outmap = new PairCorrelationMap;
    103   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     135  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
    104136    if ((*MolWalker)->ActiveFlag) {
    105137      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
    106138      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    107139      DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
     140      eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
    108141      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    109142        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)) {
     143        periodicX = *(*iter)->node;
     144        periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
     145        // go through every range in xyz and get distance
     146        for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
     147          for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
     148            for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
     149              checkX = Vector(n[0], n[1], n[2]) + periodicX;
     150              checkX.MatrixMultiplication(FullMatrix);
     151              for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++){
     152                if ((*MolOtherWalker)->ActiveFlag) {
     153                  DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
     154                  for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
     155                    DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
     156                    if ((*iter)->getId() < (*runner)->getId()){
     157                      for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
     158                        if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
    126159                          periodicOtherX = *(*runner)->node;
    127160                          periodicOtherX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    … …  
    137170                              }
    138171                        }
     172                    }
    139173                  }
     174                }
    140175              }
    141176            }
    142         }
    143177      }
    144178      delete[](FullMatrix);
    145179      delete[](FullInverseMatrix);
    146180    }
     181  }
    147182
    148183  return outmap;
    … …  
    150185
    151186/** Calculates the distance (pair) correlation between a given element and a point.
    152  * \param *out output stream for debugging
    153  * \param *molecules list of molecules structure
    154  * \param *type element or NULL (if any element)
     187 * \param *molecules list of molecules structure
     188 * \param &elements vector of elements to correlate with point
    155189 * \param *point vector to the correlation point
    156190 * \return Map of dobules with values as pairs of atom and the vector
    157191 */
    158 CorrelationToPointMap *CorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point )
     192CorrelationToPointMap *CorrelationToPoint(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Vector *point )
    159193{
    160194  Info FunctionInfo(__func__);
    161195  CorrelationToPointMap *outmap = NULL;
    162196  double distance = 0.;
     197  double *cell_size = World::getInstance().getDomain();
    163198
    164199  if (molecules->ListOfMolecules.empty()) {
    … …  
    174209      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    175210        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());
    178           DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
    179           outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
    180         }
     211        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
     212          if ((*type == NULL) || ((*iter)->type == *type)) {
     213            distance = (*iter)->node->PeriodicDistance(*point, cell_size);
     214            DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
     215            outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
     216          }
    181217      }
    182218    }
    … …  
    186222
    187223/** Calculates the distance (pair) correlation between a given element, all its periodic images and a point.
    188  * \param *out output stream for debugging
    189  * \param *molecules list of molecules structure
    190  * \param *type element or NULL (if any element)
     224 * \param *molecules list of molecules structure
     225 * \param &elements vector of elements to correlate to point
    191226 * \param *point vector to the correlation point
    192227 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
    193228 * \return Map of dobules with values as pairs of atom and the vector
    194229 */
    195 CorrelationToPointMap *PeriodicCorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point, const int ranges[NDIM] )
     230CorrelationToPointMap *PeriodicCorrelationToPoint(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Vector *point, const int ranges[NDIM] )
    196231{
    197232  Info FunctionInfo(__func__);
    … …  
    216251      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    217252        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
    221           // go through every range in xyz and get distance
    222           for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
    223             for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
    224               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);
    227                 distance = checkX.distance(*point);
    228                 DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
    229                 outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (*iter, point) ) );
    230               }
    231         }
     253        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
     254          if ((*type == NULL) || ((*iter)->type == *type)) {
     255            periodicX = *(*iter)->node;
     256            periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
     257            // go through every range in xyz and get distance
     258            for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
     259              for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
     260                for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
     261                  checkX = Vector(n[0], n[1], n[2]) + periodicX;
     262                  checkX.MatrixMultiplication(FullMatrix);
     263                  distance = checkX.distance(*point);
     264                  DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
     265                  outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (*iter, point) ) );
     266                }
     267          }
    232268      }
    233269      delete[](FullMatrix);
    … …  
    239275
    240276/** Calculates the distance (pair) correlation between a given element and a surface.
    241  * \param *out output stream for debugging
    242  * \param *molecules list of molecules structure
    243  * \param *type element or NULL (if any element)
     277 * \param *molecules list of molecules structure
     278 * \param &elements vector of elements to correlate to surface
    244279 * \param *Surface pointer to Tesselation class surface
    245280 * \param *LC LinkedCell structure to quickly find neighbouring atoms
    246281 * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
    247282 */
    248 CorrelationToSurfaceMap *CorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC )
     283CorrelationToSurfaceMap *CorrelationToSurface(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Tesselation * const Surface, const LinkedCell *LC )
    249284{
    250285  Info FunctionInfo(__func__);
    … …  
    268303      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    269304        DoLog(1) && (Log() << Verbose(1) << "\tCurrent atom is " << *(*iter) << "." << endl);
    270         if ((type == NULL) || ((*iter)->type == type)) {
    271           TriangleIntersectionList Intersections((*iter)->node,Surface,LC);
    272           distance = Intersections.GetSmallestDistance();
    273           triangle = Intersections.GetClosestTriangle();
    274           outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> ((*iter), triangle) ) );
    275         }
     305        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
     306          if ((*type == NULL) || ((*iter)->type == *type)) {
     307            TriangleIntersectionList Intersections((*iter)->node,Surface,LC);
     308            distance = Intersections.GetSmallestDistance();
     309            triangle = Intersections.GetClosestTriangle();
     310            outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> ((*iter), triangle) ) );
     311          }
    276312      }
    277313    } else {
    … …  
    287323 * axis an integer from [ -ranges[i], ranges[i] ] onto it and multiply with the domain matrix to bring it back into
    288324 * the real space. Then, we Tesselation::FindClosestTriangleToPoint() and DistanceToTrianglePlane().
    289  * \param *out output stream for debugging
    290  * \param *molecules list of molecules structure
    291  * \param *type element or NULL (if any element)
     325 * \param *molecules list of molecules structure
     326 * \param &elements vector of elements to correlate to surface
    292327 * \param *Surface pointer to Tesselation class surface
    293328 * \param *LC LinkedCell structure to quickly find neighbouring atoms
    … …  
    295330 * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
    296331 */
    297 CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
     332CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
    298333{
    299334  Info FunctionInfo(__func__);
    … …  
    322357      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    323358        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
    327           // go through every range in xyz and get distance
    328           ShortestDistance = -1.;
    329           for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
    330             for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
    331               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);
    334                 TriangleIntersectionList Intersections(&checkX,Surface,LC);
    335                 distance = Intersections.GetSmallestDistance();
    336                 triangle = Intersections.GetClosestTriangle();
    337                 if ((ShortestDistance == -1.) || (distance < ShortestDistance)) {
    338                   ShortestDistance = distance;
    339                   ShortestTriangle = triangle;
     359        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
     360          if ((*type == NULL) || ((*iter)->type == *type)) {
     361            periodicX = *(*iter)->node;
     362            periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
     363            // go through every range in xyz and get distance
     364            ShortestDistance = -1.;
     365            for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
     366              for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
     367                for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
     368                  checkX = Vector(n[0], n[1], n[2]) + periodicX;
     369                  checkX.MatrixMultiplication(FullMatrix);
     370                  TriangleIntersectionList Intersections(&checkX,Surface,LC);
     371                  distance = Intersections.GetSmallestDistance();
     372                  triangle = Intersections.GetClosestTriangle();
     373                  if ((ShortestDistance == -1.) || (distance < ShortestDistance)) {
     374                    ShortestDistance = distance;
     375                    ShortestTriangle = triangle;
     376                  }
    340377                }
    341               }
    342           // insert
    343           outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (*iter, ShortestTriangle) ) );
    344           //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
    345         }
     378            // insert
     379            outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (*iter, ShortestTriangle) ) );
     380            //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
     381          }
    346382      }
    347383      delete[](FullMatrix);
Note: See TracChangeset for help on using the changeset viewer.