Changeset e65878 for src/Analysis

Timestamp:

Oct 31, 2011, 11:49:19 AM (13 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:

870b4b

Parents:

505d05

git-author:

Frederik Heber <heber@…> (05/18/11 19:55:05)

git-committer:

Frederik Heber <heber@…> (10/31/11 11:49:19)

Message:

Selection is now done inside DipoleAngularCorrelation.

• as molecules may change over time it makes more sense to select molecules internally right now.
• DipoleAngularCorrelation() now takes Formula instead of std::vector<atom *>.
• also adjusted const'ness of some parameters.

Location:

src/Analysis

Files:

• analysis_correlation.cpp (modified) (7 diffs)
• analysis_correlation.hpp (modified) (1 diff)

src/Analysis/analysis_correlation.cpp

@@ -32 +32 @@
 #include "CodePatterns/Log.hpp"
 #include "CodePatterns/Verbose.hpp"
+#include "Descriptors/AtomOfMoleculeSelectionDescriptor.hpp"
+#include "Descriptors/MoleculeFormulaDescriptor.hpp"
 #include "Descriptors/MoleculeOfAtomSelectionDescriptor.hpp"
 #include "Formula.hpp"
@@ -100 +102 @@
  * \return range with [min, max]
  */
-range<size_t> getMaximumTrajectoryBounds(std::vector<atom *> &atoms)
+range<size_t> getMaximumTrajectoryBounds(const std::vector<atom *> &atoms)
 {
   // get highest trajectory size
@@ -121 +123 @@
 
 /** Calculates the angular dipole zero orientation from current time step.
- * \param atoms vector of atoms to calculate it for
+ * \param molecules vector of molecules to calculate dipoles of
  * \return map with orientation vector for each atomic id given in \a atoms.
  */
-std::map<atomId_t, Vector> CalculateZeroAngularDipole(std::vector<atom *> &atoms)
-{
-  // calculate molecules for this time step
-  std::set<molecule *> molecules;
-  BOOST_FOREACH(atom *_atom, atoms)
-    molecules.insert(_atom->getMolecule());
-
+std::map<atomId_t, Vector> CalculateZeroAngularDipole(const std::vector<molecule *> &molecules)
+{
   // get zero orientation for each molecule.
-  LOG(0,"STATUS: Calculating dipoles for first time step ...");
+  LOG(0,"STATUS: Calculating dipoles for current time step ...");
   std::map<atomId_t, Vector> ZeroVector;
   BOOST_FOREACH(molecule *_mol, molecules) {
@@ -154 +151 @@
  *  \a ZeroVector)
  * \param ZeroVector map with Zero orientation vector for each atom in \a atoms.
- *   Is filled from initial time step if size of map does not match size of \a atoms.
  * \param DontResetTime don't reset time to old value (triggers re-creation of bond system)
  * \return Map of doubles with values the pair of the two atoms.
  */
 DipoleAngularCorrelationMap *DipoleAngularCorrelation(
-    std::vector<atom *> &atoms,
+    const Formula &DipoleFormula,
     const size_t timestep,
-    std::map<atomId_t, Vector> &ZeroVector,
+    const std::map<atomId_t, Vector> &ZeroVector,
     const enum ResetWorldTime DoTimeReset
     )
@@ -167 +163 @@
   Info FunctionInfo(__func__);
   DipoleAngularCorrelationMap *outmap = new DipoleAngularCorrelationMap;
-
-  // get zero orientation for each molecule if not given
-  if (ZeroVector.size() != atoms.size()) {
-    LOG(1, "INFO: Mismatch size of zero orientation map ("
-        << ZeroVector.size() << ") and atom vector ("<< atoms.size() << +").");
-    ZeroVector.clear();
-    ZeroVector = CalculateZeroAngularDipole(atoms);
-  }
 
   unsigned int oldtime = 0;
@@ -187 +175 @@
 
   // get all molecules for this time step
-  LOG(0,"STATUS: Gathering molecules from " << atoms.size() << " atoms for time step " << timestep << " ...");
-  std::set<molecule *> molecules;
-  BOOST_FOREACH(atom *_atom, atoms)
-    molecules.insert(_atom->getMolecule());
+  World::getInstance().clearMoleculeSelection();
+  World::getInstance().selectAllMolecules(MoleculeByFormula(DipoleFormula));
+  World::getInstance().clearAtomSelection();
+  World::getInstance().selectAllAtoms(AtomsByMoleculeSelection());
+  std::vector<molecule *> molecules = World::getInstance().getSelectedMolecules();
+  LOG(1,"STATUS: Gathering " << molecules.size() << " molecules for time step " << timestep << ".");
 
   // calculate dipoles for each
-  LOG(0,"STATUS: Calculating dipoles for " << molecules.size() << " molecules for time step " << timestep << " ...");
+  LOG(1,"STATUS: Calculating dipoles for time step " << timestep << " ...");
   size_t i=0;
   BOOST_FOREACH(molecule *_mol, molecules) {
     const Vector Dipole = getDipole(_mol->begin(), _mol->end());
-    LOG(2,"INFO: Dipole vector at time step " << timestep << " for molecule "
+    LOG(3,"INFO: Dipole vector at time step " << timestep << " for for molecule "
         << _mol->getId() << " is " << Dipole);
+    // check that all atoms are valid (zeroVector known)
     molecule::const_iterator iter = _mol->begin();
-    ASSERT(ZeroVector.count((*iter)->getId()),
-        "DipoleAngularCorrelation() - ZeroVector for atom "+toString(**iter)+" not present.");
+    for(; iter != _mol->end(); ++iter) {
+      if (!ZeroVector.count((*iter)->getId()))
+        break;
+    }
+    if (iter != _mol->end()) {
+      ELOG(2, "Skipping molecule " << _mol->getName() << " as not all atoms have a valid zeroVector.");
+      continue;
+    } else
+      iter = _mol->begin();
+    std::map<atomId_t, Vector>::const_iterator zeroValue = ZeroVector.find((*iter)->getId()); //due to iter is const
     double angle = 0.;
     LOG(2, "INFO: ZeroVector of first atom " << **iter << " is "
-        << ZeroVector[(*iter)->getId()] << ".");
+        << zeroValue->second << ".");
     LOG(4, "INFO: Squared norm of difference vector is "
-        << (ZeroVector[(*iter)->getId()] - Dipole).NormSquared() << ".");
-    if ((ZeroVector[(*iter)->getId()] - Dipole).NormSquared() > MYEPSILON)
-      angle = Dipole.Angle(ZeroVector[(*iter)->getId()]) * (180./M_PI);
+        << (zeroValue->second - Dipole).NormSquared() << ".");
+    if ((zeroValue->second - Dipole).NormSquared() > MYEPSILON)
+      angle = Dipole.Angle(zeroValue->second) * (180./M_PI);
     else
       LOG(2, "INFO: Both vectors (almost) coincide, numerically unstable, angle set to zero.");
@@ -216 +215 @@
     ++i;
   }
-  LOG(0,"STATUS: Done with calculating dipoles.");
+  LOG(1,"STATUS: Done with calculating dipoles.");
 
   if (DoTimeReset == DoResetTime) {

src/Analysis/analysis_correlation.hpp

@@ -60 +60 @@
 /********************************************** declarations *******************************/
 
-range<size_t> getMaximumTrajectoryBounds(std::vector<atom *> &atoms);
-std::map<atomId_t, Vector> CalculateZeroAngularDipole(std::vector<atom *> &atoms);
-
-DipoleAngularCorrelationMap *DipoleAngularCorrelation(std::vector<atom *> &atoms, const size_t timestep, std::map<atomId_t, Vector> &ZeroVector, const enum ResetWorldTime DoTimeReset = DontResetTime);
+range<size_t> getMaximumTrajectoryBounds(const std::vector<atom *> &atoms);
+std::map<atomId_t, Vector> CalculateZeroAngularDipole(const std::vector<molecule *> &molecules);
+
+DipoleAngularCorrelationMap *DipoleAngularCorrelation(const Formula &DipoleFormula, const size_t timestep, const std::map<atomId_t, Vector> &ZeroVector, const enum ResetWorldTime DoTimeReset = DontResetTime);
 DipoleCorrelationMap *DipoleCorrelation(std::vector<molecule *> &molecules);
 PairCorrelationMap *PairCorrelation(std::vector<molecule *> &molecules, const std::vector<const element *> &elements);