Changeset c78d44 for src


Ignore:
Timestamp:
Jun 9, 2010, 11:12:08 AM (15 years ago)
Author:
Frederik Heber <heber@…>
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, 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:
104524
Parents:
8540f0
Message:

Change to PairCorrelation...() functions, accept now vector<element *> instead of ptrs to element *.

  • modified ParseCommndLineOptions for all three cases 'C': Pair (E), Point (P) and Surface (S).
  • modified all Analysis(Pair)Correlation...() unit tests to have a vector of elements be delivered.
Location:
src
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • src/analysis_correlation.cpp

    r8540f0 rc78d44  
    2323/** Calculates the pair correlation between given elements.
    2424 * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
    25  * \param *out output stream for debugging
    2625 * \param *molecules list of molecules structure
    27  * \param *type1 first element or NULL (if any element)
    28  * \param *type2 second element or NULL (if any element)
     26 * \param &elements vector of elements to correlate
    2927 * \return Map of doubles with values the pair of the two atoms.
    3028 */
    31 PairCorrelationMap *PairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const type2 )
     29PairCorrelationMap *PairCorrelation(MoleculeListClass * const &molecules, const std::vector<element *> &elements)
    3230{
    3331  Info FunctionInfo(__func__);
    3432  PairCorrelationMap *outmap = NULL;
    3533  double distance = 0.;
     34  double *domain = World::getInstance().getDomain();
    3635
    3736  if (molecules->ListOfMolecules.empty()) {
    … …  
    4140  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    4241    (*MolWalker)->doCountAtoms();
     42
     43  // create all possible pairs of elements
     44  set <pair<element *, element *> > PairsOfElements;
     45  if (elements.size() >= 2) {
     46    for (vector<element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
     47      for (vector<element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
     48        if (type1 != type2) {
     49          PairsOfElements.insert( pair<element *, element*>(*type1,*type2) );
     50          DoLog(1) && (Log() << Verbose(1) << "Creating element pair " << (*type1)->symbol << " and " << (*type2)->symbol << "." << endl);
     51        }
     52  } else if (elements.size() == 1) { // one to all are valid
     53    element *elemental = *elements.begin();
     54    PairsOfElements.insert( pair<element *, element*>(elemental,(element *)NULL) );
     55    PairsOfElements.insert( pair<element *, element*>((element *)NULL,elemental) );
     56  } else { // all elements valid
     57    PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
     58  }
     59
    4360  outmap = new PairCorrelationMap;
    4461  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
    … …  
    4865      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    4966        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
    50         if ((type1 == NULL) || ((*iter)->type == type1)) {
    51           for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++){
    52             if ((*MolOtherWalker)->ActiveFlag) {
    53               DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
    54               for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
    55                 DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
    56                 if ((*iter)->getId() < (*runner)->getId()){
    57                   if ((type2 == NULL) || ((*runner)->type == type2)) {
    58                     distance = (*iter)->node->PeriodicDistance(*(*runner)->node,  World::getInstance().getDomain());
     67        for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++){
     68          if ((*MolOtherWalker)->ActiveFlag) {
     69            DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
     70            for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
     71              DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
     72              if ((*iter)->getId() < (*runner)->getId()){
     73                for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
     74                  if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
     75                    distance = (*iter)->node->PeriodicDistance(*(*runner)->node,  domain);
    5976                    //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
    6077                    outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
    6178                  }
    62                 }
    6379              }
    6480            }
    … …  
    7389/** Calculates the pair correlation between given elements.
    7490 * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
    75  * \param *out output stream for debugging
    7691 * \param *molecules list of molecules structure
    77  * \param *type1 first element or NULL (if any element)
    78  * \param *type2 second element or NULL (if any element)
     92 * \param &elements vector of elements to correlate
    7993 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
    8094 * \return Map of doubles with values the pair of the two atoms.
    8195 */
    82 PairCorrelationMap *PeriodicPairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const type2, const int ranges[NDIM] )
     96PairCorrelationMap *PeriodicPairCorrelation(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const int ranges[NDIM] )
    8397{
    8498  Info FunctionInfo(__func__);
    … …  
    98112  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    99113    (*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
    100132  outmap = new PairCorrelationMap;
    101   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     133  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
    102134    if ((*MolWalker)->ActiveFlag) {
    103135      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
    104136      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    105137      DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
     138      eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
    106139      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    107140        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
    108         if ((type1 == NULL) || ((*iter)->type == type1)) {
    109           periodicX = *(*iter)->node;
    110           periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    111           // go through every range in xyz and get distance
    112           for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
    113             for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
    114               for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
    115                 checkX = Vector(n[0], n[1], n[2]) + periodicX;
    116                 checkX.MatrixMultiplication(FullMatrix);
    117                 for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
    118                   if ((*MolOtherWalker)->ActiveFlag) {
    119                     DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
    120                     for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
    121                       DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
    122                       if ((*iter)->nr < (*runner)->nr)
    123                         if ((type2 == NULL) || ((*runner)->type == type2)) {
     141        periodicX = *(*iter)->node;
     142        periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
     143        // go through every range in xyz and get distance
     144        for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
     145          for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
     146            for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
     147              checkX = Vector(n[0], n[1], n[2]) + periodicX;
     148              checkX.MatrixMultiplication(FullMatrix);
     149              for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++){
     150                if ((*MolOtherWalker)->ActiveFlag) {
     151                  DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
     152                  for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
     153                    DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
     154                    if ((*iter)->getId() < (*runner)->getId()){
     155                      for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
     156                        if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
    124157                          periodicOtherX = *(*runner)->node;
    125158                          periodicOtherX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    … …  
    135168                              }
    136169                        }
     170                    }
    137171                  }
     172                }
    138173              }
    139174            }
    140         }
    141175      }
    142176      delete[](FullMatrix);
    143177      delete[](FullInverseMatrix);
    144178    }
     179  }
     180
     181//  outmap = new PairCorrelationMap;
     182//  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     183//    if ((*MolWalker)->ActiveFlag) {
     184//      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
     185//      double * FullInverseMatrix = InverseMatrix(FullMatrix);
     186//      DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
     187//      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
     188//        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
     189//        if ((type1 == NULL) || ((*iter)->type == type1)) {
     190//          periodicX = *(*iter)->node;
     191//          periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
     192//          // go through every range in xyz and get distance
     193//          for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
     194//            for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
     195//              for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
     196//                checkX = Vector(n[0], n[1], n[2]) + periodicX;
     197//                checkX.MatrixMultiplication(FullMatrix);
     198//                for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
     199//                  if ((*MolOtherWalker)->ActiveFlag) {
     200//                    DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
     201//                    for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
     202//                      DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
     203//                      if ((*iter)->nr < (*runner)->nr)
     204//                        if ((type2 == NULL) || ((*runner)->type == type2)) {
     205//                          periodicOtherX = *(*runner)->node;
     206//                          periodicOtherX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
     207//                          // go through every range in xyz and get distance
     208//                          for (Othern[0]=-ranges[0]; Othern[0] <= ranges[0]; Othern[0]++)
     209//                            for (Othern[1]=-ranges[1]; Othern[1] <= ranges[1]; Othern[1]++)
     210//                              for (Othern[2]=-ranges[2]; Othern[2] <= ranges[2]; Othern[2]++) {
     211//                                checkOtherX = Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX;
     212//                                checkOtherX.MatrixMultiplication(FullMatrix);
     213//                                distance = checkX.distance(checkOtherX);
     214//                                //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
     215//                                outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
     216//                              }
     217//                        }
     218//                  }
     219//              }
     220//            }
     221//        }
     222//      }
     223//      delete[](FullMatrix);
     224//      delete[](FullInverseMatrix);
     225//    }
    145226
    146227  return outmap;
    … …  
    148229
    149230/** Calculates the distance (pair) correlation between a given element and a point.
    150  * \param *out output stream for debugging
    151231 * \param *molecules list of molecules structure
    152  * \param *type element or NULL (if any element)
     232 * \param &elements vector of elements to correlate with point
    153233 * \param *point vector to the correlation point
    154234 * \return Map of dobules with values as pairs of atom and the vector
    155235 */
    156 CorrelationToPointMap *CorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point )
     236CorrelationToPointMap *CorrelationToPoint(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Vector *point )
    157237{
    158238  Info FunctionInfo(__func__);
    … …  
    172252      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    173253        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
    174         if ((type == NULL) || ((*iter)->type == type)) {
    175           distance = (*iter)->node->PeriodicDistance(*point, World::getInstance().getDomain());
    176           DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
    177           outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
    178         }
     254        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
     255          if ((*type == NULL) || ((*iter)->type == *type)) {
     256            distance = (*iter)->node->PeriodicDistance(*point, World::getInstance().getDomain());
     257            DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
     258            outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
     259          }
    179260      }
    180261    }
    … …  
    184265
    185266/** Calculates the distance (pair) correlation between a given element, all its periodic images and a point.
    186  * \param *out output stream for debugging
    187267 * \param *molecules list of molecules structure
    188  * \param *type element or NULL (if any element)
     268 * \param &elements vector of elements to correlate to point
    189269 * \param *point vector to the correlation point
    190270 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
    191271 * \return Map of dobules with values as pairs of atom and the vector
    192272 */
    193 CorrelationToPointMap *PeriodicCorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point, const int ranges[NDIM] )
     273CorrelationToPointMap *PeriodicCorrelationToPoint(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Vector *point, const int ranges[NDIM] )
    194274{
    195275  Info FunctionInfo(__func__);
    … …  
    214294      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    215295        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
    216         if ((type == NULL) || ((*iter)->type == type)) {
    217           periodicX = *(*iter)->node;
    218           periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    219           // go through every range in xyz and get distance
    220           for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
    221             for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
    222               for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
    223                 checkX = Vector(n[0], n[1], n[2]) + periodicX;
    224                 checkX.MatrixMultiplication(FullMatrix);
    225                 distance = checkX.distance(*point);
    226                 DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
    227                 outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (*iter, point) ) );
    228               }
    229         }
     296        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
     297          if ((*type == NULL) || ((*iter)->type == *type)) {
     298            periodicX = *(*iter)->node;
     299            periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
     300            // go through every range in xyz and get distance
     301            for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
     302              for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
     303                for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
     304                  checkX = Vector(n[0], n[1], n[2]) + periodicX;
     305                  checkX.MatrixMultiplication(FullMatrix);
     306                  distance = checkX.distance(*point);
     307                  DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
     308                  outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (*iter, point) ) );
     309                }
     310          }
    230311      }
    231312      delete[](FullMatrix);
    … …  
    237318
    238319/** Calculates the distance (pair) correlation between a given element and a surface.
    239  * \param *out output stream for debugging
    240320 * \param *molecules list of molecules structure
    241  * \param *type element or NULL (if any element)
     321 * \param &elements vector of elements to correlate to surface
    242322 * \param *Surface pointer to Tesselation class surface
    243323 * \param *LC LinkedCell structure to quickly find neighbouring atoms
    244324 * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
    245325 */
    246 CorrelationToSurfaceMap *CorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC )
     326CorrelationToSurfaceMap *CorrelationToSurface(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Tesselation * const Surface, const LinkedCell *LC )
    247327{
    248328  Info FunctionInfo(__func__);
    … …  
    266346      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    267347        DoLog(1) && (Log() << Verbose(1) << "\tCurrent atom is " << *(*iter) << "." << endl);
    268         if ((type == NULL) || ((*iter)->type == type)) {
    269           TriangleIntersectionList Intersections((*iter)->node,Surface,LC);
    270           distance = Intersections.GetSmallestDistance();
    271           triangle = Intersections.GetClosestTriangle();
    272           outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> ((*iter), triangle) ) );
    273         }
     348        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
     349          if ((*type == NULL) || ((*iter)->type == *type)) {
     350            TriangleIntersectionList Intersections((*iter)->node,Surface,LC);
     351            distance = Intersections.GetSmallestDistance();
     352            triangle = Intersections.GetClosestTriangle();
     353            outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> ((*iter), triangle) ) );
     354          }
    274355      }
    275356    } else {
    … …  
    285366 * axis an integer from [ -ranges[i], ranges[i] ] onto it and multiply with the domain matrix to bring it back into
    286367 * the real space. Then, we Tesselation::FindClosestTriangleToPoint() and DistanceToTrianglePlane().
    287  * \param *out output stream for debugging
    288368 * \param *molecules list of molecules structure
    289  * \param *type element or NULL (if any element)
     369 * \param &elements vector of elements to correlate to surface
    290370 * \param *Surface pointer to Tesselation class surface
    291371 * \param *LC LinkedCell structure to quickly find neighbouring atoms
    … …  
    293373 * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
    294374 */
    295 CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
     375CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
    296376{
    297377  Info FunctionInfo(__func__);
    … …  
    320400      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    321401        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
    322         if ((type == NULL) || ((*iter)->type == type)) {
    323           periodicX = *(*iter)->node;
    324           periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    325           // go through every range in xyz and get distance
    326           ShortestDistance = -1.;
    327           for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
    328             for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
    329               for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
    330                 checkX = Vector(n[0], n[1], n[2]) + periodicX;
    331                 checkX.MatrixMultiplication(FullMatrix);
    332                 TriangleIntersectionList Intersections(&checkX,Surface,LC);
    333                 distance = Intersections.GetSmallestDistance();
    334                 triangle = Intersections.GetClosestTriangle();
    335                 if ((ShortestDistance == -1.) || (distance < ShortestDistance)) {
    336                   ShortestDistance = distance;
    337                   ShortestTriangle = triangle;
     402        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
     403          if ((*type == NULL) || ((*iter)->type == *type)) {
     404            periodicX = *(*iter)->node;
     405            periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
     406            // go through every range in xyz and get distance
     407            ShortestDistance = -1.;
     408            for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
     409              for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
     410                for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
     411                  checkX = Vector(n[0], n[1], n[2]) + periodicX;
     412                  checkX.MatrixMultiplication(FullMatrix);
     413                  TriangleIntersectionList Intersections(&checkX,Surface,LC);
     414                  distance = Intersections.GetSmallestDistance();
     415                  triangle = Intersections.GetClosestTriangle();
     416                  if ((ShortestDistance == -1.) || (distance < ShortestDistance)) {
     417                    ShortestDistance = distance;
     418                    ShortestTriangle = triangle;
     419                  }
    338420                }
    339               }
    340           // insert
    341           outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (*iter, ShortestTriangle) ) );
    342           //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
    343         }
     421            // insert
     422            outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (*iter, ShortestTriangle) ) );
     423            //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
     424          }
    344425      }
    345426      delete[](FullMatrix);

  • src/analysis_correlation.hpp

    r8540f0 rc78d44  
    4545/********************************************** declarations *******************************/
    4646
    47 PairCorrelationMap *PairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const  type2 );
    48 CorrelationToPointMap *CorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point );
    49 CorrelationToSurfaceMap *CorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC );
    50 PairCorrelationMap *PeriodicPairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const  type2, const int ranges[NDIM] );
    51 CorrelationToPointMap *PeriodicCorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point, const int ranges[NDIM] );
    52 CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] );
     47PairCorrelationMap *PairCorrelation(MoleculeListClass * const &molecules, const std::vector<element *> &elements);
     48CorrelationToPointMap *CorrelationToPoint(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Vector *point );
     49CorrelationToSurfaceMap *CorrelationToSurface(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Tesselation * const Surface, const LinkedCell *LC );
     50PairCorrelationMap *PeriodicPairCorrelation(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const int ranges[NDIM] );
     51CorrelationToPointMap *PeriodicCorrelationToPoint(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Vector *point, const int ranges[NDIM] );
     52CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] );
    5353int GetBin ( const double value, const double BinWidth, const double BinStart );
    5454void OutputCorrelation( ofstream * const file, const BinPairMap * const map );

  • src/builder.cpp

    r8540f0 rc78d44  
    17741774                          const double BinEnd = atof(argv[argptr+6]);
    17751775
    1776                           const element *elemental = periode->FindElement((const int) atoi(argv[argptr+1]));
    1777                           const element *elemental2 = periode->FindElement((const int) atoi(argv[argptr+2]));
     1776                          std::vector<element *> elements;
     1777                          elements.push_back(periode->FindElement((const int) atoi(argv[argptr+1])));
     1778                          elements.push_back(periode->FindElement((const int) atoi(argv[argptr+2])));
    17781779                          PairCorrelationMap *correlationmap = NULL;
    17791780                          if (periodic)
    1780                             correlationmap = PeriodicPairCorrelation(molecules, elemental, elemental2, ranges);
     1781                            correlationmap = PeriodicPairCorrelation(molecules, elements, ranges);
    17811782                          else
    1782                             correlationmap = PairCorrelation(molecules, elemental, elemental2);
     1783                            correlationmap = PairCorrelation(molecules, elements);
    17831784                          OutputPairCorrelation(&output, correlationmap);
    17841785                          BinPairMap *binmap = BinData( correlationmap, 0.5, BinStart, BinEnd );
    … …  
    18051806                          const double BinEnd = atof(argv[argptr+8]);
    18061807
    1807                           const element *elemental = periode->FindElement((const int) atoi(argv[argptr+1]));
     1808                          std::vector<element *> elements;
     1809                          elements.push_back(periode->FindElement((const int) atoi(argv[argptr+1])));
    18081810                          Vector *Point = new Vector((const double) atof(argv[argptr+1]),(const double) atof(argv[argptr+2]),(const double) atof(argv[argptr+3]));
    18091811                          CorrelationToPointMap *correlationmap = NULL;
    18101812                          if (periodic)
    1811                             correlationmap  = PeriodicCorrelationToPoint(molecules, elemental, Point, ranges);
     1813                            correlationmap  = PeriodicCorrelationToPoint(molecules, elements, Point, ranges);
    18121814                          else
    1813                             correlationmap = CorrelationToPoint(molecules, elemental, Point);
     1815                            correlationmap = CorrelationToPoint(molecules, elements, Point);
    18141816                          OutputCorrelationToPoint(&output, correlationmap);
    18151817                          BinPairMap *binmap = BinData( correlationmap, 0.5, BinStart, BinEnd );
    … …  
    18651867                          }
    18661868                          LCList = new LinkedCell(Boundary, LCWidth);
    1867                           const element *elemental = periode->FindElement((const int) atoi(argv[argptr+1]));
     1869                          std::vector<element *> elements;
     1870                          elements.push_back(periode->FindElement((const int) atoi(argv[argptr+1])));
    18681871                          FindNonConvexBorder(Boundary, TesselStruct, LCList, radius, NULL);
    18691872                          CorrelationToSurfaceMap *surfacemap = NULL;
    18701873                          if (periodic)
    1871                             surfacemap = PeriodicCorrelationToSurface( molecules, elemental, TesselStruct, LCList, ranges);
     1874                            surfacemap = PeriodicCorrelationToSurface( molecules, elements, TesselStruct, LCList, ranges);
    18721875                          else
    1873                             surfacemap = CorrelationToSurface( molecules, elemental, TesselStruct, LCList);
     1876                            surfacemap = CorrelationToSurface( molecules, elements, TesselStruct, LCList);
    18741877                          OutputCorrelationToSurface(&output, surfacemap);
    18751878                          // check whether radius was appropriate

  • src/unittests/AnalysisCorrelationToPointUnitTest.cpp

    r8540f0 rc78d44  
    4444  point = NULL;
    4545
    46   // construct molecule (tetraeder of hydrogens)
     46  // construct element list
     47  std::vector<element *> elements;
    4748  hydrogen = World::getInstance().getPeriode()->FindElement(1);
    4849  CPPUNIT_ASSERT(hydrogen != NULL && "hydrogen element not found");
     50  elements.push_back(hydrogen);
     51  // construct molecule (tetraeder of hydrogens)
    4952  TestMolecule = World::getInstance().createMolecule();
    5053  Walker = World::getInstance().createAtom();
    … …  
    7679
    7780  // init maps
    78   pointmap = CorrelationToPoint( (MoleculeListClass * const)TestList, (const element * const)hydrogen, (const Vector *)point );
     81  pointmap = CorrelationToPoint( (MoleculeListClass * const)TestList, elements, (const Vector *)point );
    7982  binmap = NULL;
    8083

  • src/unittests/AnalysisCorrelationToPointUnitTest.hpp

    r8540f0 rc78d44  
    3737      MoleculeListClass *TestList;
    3838      molecule *TestMolecule;
    39       const element *hydrogen;
     39      element *hydrogen;
    4040
    4141      CorrelationToPointMap *pointmap;

  • src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp

    r8540f0 rc78d44  
    5252  LC = NULL;
    5353
     54  // prepare element list
     55  hydrogen = World::getInstance().getPeriode()->FindElement(1);
     56  CPPUNIT_ASSERT(hydrogen != NULL && "hydrogen element not found");
     57  elements.clear();
     58
    5459  // construct molecule (tetraeder of hydrogens) base
    55   hydrogen = World::getInstance().getPeriode()->FindElement(1);
    5660  TestSurfaceMolecule = World::getInstance().createMolecule();
    5761
    … …  
    151155{
    152156  // do the pair correlation
    153   surfacemap = CorrelationToSurface( TestList, hydrogen, Surface, LC );
     157  elements.push_back(hydrogen);
     158  surfacemap = CorrelationToSurface( TestList, elements, Surface, LC );
    154159//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
    155160  CPPUNIT_ASSERT( surfacemap != NULL );
    … …  
    160165{
    161166  BinPairMap::iterator tester;
    162   surfacemap = CorrelationToSurface( TestList, hydrogen, Surface, LC );
     167  elements.push_back(hydrogen);
     168  surfacemap = CorrelationToSurface( TestList, elements, Surface, LC );
    163169  // put pair correlation into bins and check with no range
    164170//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
    … …  
    175181{
    176182  BinPairMap::iterator tester;
    177   surfacemap = CorrelationToSurface( TestList, hydrogen, Surface, LC );
     183  elements.push_back(hydrogen);
     184  surfacemap = CorrelationToSurface( TestList, elements, Surface, LC );
    178185//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
    179186  // ... and check with [0., 2.] range
    … …  
    193200{
    194201  BinPairMap::iterator tester;
    195   surfacemap = CorrelationToSurface( TestList, carbon, Surface, LC );
     202  elements.push_back(carbon);
     203  surfacemap = CorrelationToSurface( TestList, elements, Surface, LC );
    196204//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
    197205  // put pair correlation into bins and check with no range
    … …  
    212220{
    213221  BinPairMap::iterator tester;
    214   surfacemap = CorrelationToSurface( TestList, carbon, Surface, LC );
     222  elements.push_back(carbon);
     223  surfacemap = CorrelationToSurface( TestList, elements, Surface, LC );
    215224//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
    216225  // ... and check with [0., 2.] range

  • src/unittests/AnalysisCorrelationToSurfaceUnitTest.hpp

    r8540f0 rc78d44  
    4545      MoleculeListClass *TestList;
    4646      molecule *TestSurfaceMolecule;
    47       const element *hydrogen;
    48       const element *carbon;
     47      element *hydrogen;
     48      element *carbon;
     49      std::vector<element *> elements;
    4950
    5051      CorrelationToSurfaceMap *surfacemap;

  • src/unittests/AnalysisPairCorrelationUnitTest.cpp

    r8540f0 rc78d44  
    4747  binmap = NULL;
    4848
     49  // construct element list
     50  std::vector<element *> elements;
     51  hydrogen = World::getInstance().getPeriode()->FindElement(1);
     52  CPPUNIT_ASSERT(hydrogen != NULL && "hydrogen element not found");
     53  elements.push_back(hydrogen);
     54  elements.push_back(hydrogen);
     55
    4956  // construct molecule (tetraeder of hydrogens)
    50   hydrogen = World::getInstance().getPeriode()->FindElement(1);
    5157  TestMolecule = World::getInstance().createMolecule();
    5258  Walker = World::getInstance().createAtom();
    … …  
    7581
    7682  // init maps
    77   correlationmap = PairCorrelation( TestList, hydrogen, hydrogen );
     83  correlationmap = PairCorrelation( TestList, elements);
    7884  binmap = NULL;
    7985

  • src/unittests/AnalysisPairCorrelationUnitTest.hpp

    r8540f0 rc78d44  
    3737      MoleculeListClass *TestList;
    3838      molecule *TestMolecule;
    39       const element *hydrogen;
     39      element *hydrogen;
    4040
    4141      PairCorrelationMap *correlationmap;
Note: See TracChangeset for help on using the changeset viewer.