Changeset e65878

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

Files:

• Actions/AnalysisAction/DipoleAngularCorrelationAction.cpp (modified) (3 diffs)
• Actions/AnalysisAction/DipoleAngularCorrelationAction.def (modified) (1 diff)
• Analysis/analysis_correlation.cpp (modified) (7 diffs)
• Analysis/analysis_correlation.hpp (modified) (1 diff)

src/Actions/AnalysisAction/DipoleAngularCorrelationAction.cpp

@@ -24 +24 @@
 #include "linkedcell.hpp"
 #include "CodePatterns/Log.hpp"
+#include "Descriptors/AtomOfMoleculeSelectionDescriptor.hpp"
+#include "Descriptors/MoleculeFormulaDescriptor.hpp"
 #include "Element/element.hpp"
-#include "molecule.hpp"
 #include "Element/periodentafel.hpp"
 #include "LinearAlgebra/Vector.hpp"
+#include "molecule.hpp"
 #include "World.hpp"
 #include "WorldTime.hpp"
@@ -53 +55 @@
 
   // get selected atoms
-  std::vector<atom*> atoms = World::getInstance().getSelectedAtoms();
-  ASSERT(atoms.size() != 0, "AnalysisDipoleAngularCorrelationAction() - not atoms selected.");
+  std::vector<atom*> old_atom_selection = World::getInstance().getSelectedAtoms();
+  std::vector<molecule*> old_molecule_selection = World::getInstance().getSelectedMolecules();
 
   // get current time step
   const unsigned int oldtime = WorldTime::getTime();
 
-  // obtain zero dipole orientation
+  // select atoms and obtain zero dipole orientation
+  Formula DipoleFormula(params.DipoleFormula);
   World::getInstance().setTime(0);
-  std::map<atomId_t, Vector> ZeroVector = CalculateZeroAngularDipole(atoms);
+  World::getInstance().clearMoleculeSelection(); // TODO: This should be done in setTime or where molecules are re-done
+  World::getInstance().selectAllMolecules(MoleculeByFormula(DipoleFormula));
+  std::vector<molecule *> molecules = World::getInstance().getSelectedMolecules();
+  std::map<atomId_t, Vector> ZeroVector = CalculateZeroAngularDipole(molecules);
 
   // go through each step of common trajectory of all atoms in set
+  World::getInstance().clearAtomSelection();
+  World::getInstance().selectAllAtoms(AtomsByMoleculeSelection());
+  std::vector<atom *> atoms = World::getInstance().getSelectedAtoms();
   range<size_t> timesteps = getMaximumTrajectoryBounds(atoms);
   for (size_t step = 0; step < timesteps.first; ++step) {
     // calculate dipoles relative to zero orientation
     DipoleAngularCorrelationMap *correlationmap = NULL;
-    correlationmap = DipoleAngularCorrelation(atoms, step, ZeroVector, DontResetTime);
+    correlationmap = DipoleAngularCorrelation(DipoleFormula, step, ZeroVector, DontResetTime);
 
     // prepare step string in filename
@@ -102 +111 @@
   World::getInstance().setTime(oldtime);
 
+  // reset to old selections
+  World::getInstance().clearAtomSelection();
+  BOOST_FOREACH(atom *_atom, old_atom_selection) {
+    World::getInstance().selectAtom(_atom);
+  }
+  World::getInstance().clearMoleculeSelection();
+  BOOST_FOREACH(molecule *_mol, old_molecule_selection) {
+    World::getInstance().selectMolecule(_mol);
+  }
+
   // exit
   return Action::success;

src/Actions/AnalysisAction/DipoleAngularCorrelationAction.def

@@ -13 +13 @@
 // ValueStorage by the token "Z" -> first column: int, Z, "Z"
 // "undefine" if no parameters are required, use (NODEFAULT) for each (undefined) default value
-#define paramtypes (double)(double)(double)(boost::filesystem::path)(boost::filesystem::path)(bool)
-#define paramreferences (BinStart)(BinWidth)(BinEnd)(outputname)(binoutputname)(periodic)
-#define paramtokens ("bin-start")("bin-width")("bin-end")("output-file")("bin-output-file")("periodic")
-#define paramdescriptions ("start of the first bin")("width of the bins")("start of the last bin")("name of the output file")("name of the bin output file")("system is constraint to periodic boundary conditions")
-#define paramdefaults (NODEFAULT)("0.5")(NODEFAULT)(NODEFAULT)(NODEFAULT)("0")
+#define paramtypes (std::string)(double)(double)(double)(boost::filesystem::path)(boost::filesystem::path)(bool)
+#define paramreferences (DipoleFormula)(BinStart)(BinWidth)(BinEnd)(outputname)(binoutputname)(periodic)
+#define paramtokens ("dipole-angular-correlation")("bin-start")("bin-width")("bin-end")("output-file")("bin-output-file")("periodic")
+#define paramdescriptions ("formula of molecules to calculate dipole of")("start of the first bin")("width of the bins")("start of the last bin")("name of the output file")("name of the bin output file")("system is constraint to periodic boundary conditions")
+#define paramdefaults (NODEFAULT)(NODEFAULT)("0.5")(NODEFAULT)(NODEFAULT)(NODEFAULT)("0")
 
 // some defines for all the names, you may use ACTION, STATE and PARAMS

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);