Changeset 7ba268 for src/Analysis

Timestamp:

Nov 3, 2011, 9:52:29 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:

bbff92

Parents:

13510b (diff), 5837dd (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:

Merge branch 'DipoleAngularCorrelation' into mergeMaster

Conflicts:

src/World.cpp

• tiny conflict in World.cpp

Location:

src/Analysis

Files:

• analysis_correlation.cpp (modified) (12 diffs)
• analysis_correlation.hpp (modified) (3 diffs)

src/Analysis/analysis_correlation.cpp

@@ -22 +22 @@
 #include <iostream>
 #include <iomanip>
+#include <limits>
 
 #include "atom.hpp"
@@ -30 +31 @@
 #include "CodePatterns/Info.hpp"
 #include "CodePatterns/Log.hpp"
+#include "CodePatterns/Verbose.hpp"
+#include "Descriptors/AtomOfMoleculeSelectionDescriptor.hpp"
+#include "Descriptors/MoleculeFormulaDescriptor.hpp"
+#include "Descriptors/MoleculeOfAtomSelectionDescriptor.hpp"
 #include "Formula.hpp"
+#include "LinearAlgebra/Vector.hpp"
+#include "LinearAlgebra/RealSpaceMatrix.hpp"
 #include "molecule.hpp"
 #include "Tesselation/tesselation.hpp"
@@ -36 +43 @@
 #include "Tesselation/triangleintersectionlist.hpp"
 #include "World.hpp"
-#include "LinearAlgebra/Vector.hpp"
-#include "LinearAlgebra/RealSpaceMatrix.hpp"
-#include "CodePatterns/Verbose.hpp"
-#include "World.hpp"
-#include "Box.hpp"
+#include "WorldTime.hpp"
 
 #include "analysis_correlation.hpp"
@@ -67 +70 @@
     // go through all bonds
     const BondList& ListOfBonds = (*atomiter)->getListOfBonds();
+    ASSERT(ListOfBonds.begin() != ListOfBonds.end(),
+        "getDipole() - no bonds in molecule!");
     for (BondList::const_iterator bonditer = ListOfBonds.begin();
         bonditer != ListOfBonds.end();
@@ -78 +83 @@
         BondDipoleVector.Normalize();
         BondDipoleVector *= DeltaEN;
+        LOG(3,"INFO: Dipole vector from bond " << **bonditer << " is " << BondDipoleVector);
         DipoleVector += BondDipoleVector;
         SumOfVectors++;
@@ -83 +89 @@
     }
   }
-  DipoleVector *= 1./(double)SumOfVectors;
+  LOG(3,"INFO: Sum over all bond dipole vectors is "
+      << DipoleVector << " with " << SumOfVectors << " in total.");
+  if (SumOfVectors != 0)
+    DipoleVector *= 1./(double)SumOfVectors;
   DoLog(1) && (Log() << Verbose(1) << "Resulting dipole vector is " << DipoleVector << std::endl);
 
@@ -89 +98 @@
 };
 
+/** Calculate minimum and maximum amount of trajectory steps by going through given atomic trajectories.
+ * \param vector of atoms whose trajectories to check for [min,max]
+ * \return range with [min, max]
+ */
+range<size_t> getMaximumTrajectoryBounds(const std::vector<atom *> &atoms)
+{
+  // get highest trajectory size
+  LOG(0,"STATUS: Retrieving maximum amount of time steps ...");
+  if (atoms.size() == 0)
+    return range<size_t>(0,0);
+  size_t max_timesteps = std::numeric_limits<size_t>::min();
+  size_t min_timesteps = std::numeric_limits<size_t>::max();
+  BOOST_FOREACH(atom *_atom, atoms) {
+    if (_atom->getTrajectorySize() > max_timesteps)
+      max_timesteps  = _atom->getTrajectorySize();
+    if (_atom->getTrajectorySize() < min_timesteps)
+      min_timesteps = _atom->getTrajectorySize();
+  }
+  LOG(1,"INFO: Minimum number of time steps found is " << min_timesteps);
+  LOG(1,"INFO: Maximum number of time steps found is " << max_timesteps);
+
+  return range<size_t>(min_timesteps, max_timesteps);
+}
+
+/** Calculates the angular dipole zero orientation from current time step.
+ * \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(const std::vector<molecule *> &molecules)
+{
+  // get zero orientation for each molecule.
+  LOG(0,"STATUS: Calculating dipoles for current time step ...");
+  std::map<atomId_t, Vector> ZeroVector;
+  BOOST_FOREACH(molecule *_mol, molecules) {
+    const Vector Dipole = getDipole(_mol->begin(), _mol->end());
+    for(molecule::const_iterator iter = _mol->begin(); iter != _mol->end(); ++iter)
+      ZeroVector[(*iter)->getId()] = Dipole;
+    LOG(2,"INFO: Zero alignment for molecule " << _mol->getId() << " is " << Dipole);
+  }
+  LOG(1,"INFO: We calculated zero orientation for a total of " << molecules.size() << " molecule(s).");
+
+  return ZeroVector;
+}
+
 /** Calculates the dipole angular correlation for given molecule type.
+ * Calculate the change of the dipole orientation angle over time.
+ * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
+ * Angles are given in degrees.
+ * \param &atoms list of atoms of the molecules taking part (Note: molecules may
+ * change over time as bond structure is recalculated, hence we need the atoms)
+ * \param timestep time step to calculate angular correlation for (relative to
+ *  \a ZeroVector)
+ * \param ZeroVector map with Zero orientation vector for each atom in \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(
+    const Formula &DipoleFormula,
+    const size_t timestep,
+    const std::map<atomId_t, Vector> &ZeroVector,
+    const enum ResetWorldTime DoTimeReset
+    )
+{
+  Info FunctionInfo(__func__);
+  DipoleAngularCorrelationMap *outmap = new DipoleAngularCorrelationMap;
+
+  unsigned int oldtime = 0;
+  if (DoTimeReset == DoResetTime) {
+    // store original time step
+    oldtime = WorldTime::getTime();
+  }
+
+  // set time step
+  LOG(0,"STATUS: Stepping onto to time step " << timestep << ".");
+  World::getInstance().setTime(timestep);
+
+  // get all molecules for this time step
+  World::getInstance().clearMoleculeSelection();
+  World::getInstance().selectAllMolecules(MoleculeByFormula(DipoleFormula));
+  std::vector<molecule *> molecules = World::getInstance().getSelectedMolecules();
+  LOG(1,"INFO: There are " << molecules.size() << " molecules for time step " << timestep << ".");
+
+  // calculate dipoles for each
+  LOG(0,"STATUS: Calculating dipoles for time step " << timestep << " ...");
+  size_t i=0;
+  size_t Counter_rejections = 0;
+  BOOST_FOREACH(molecule *_mol, molecules) {
+    const Vector Dipole = getDipole(_mol->begin(), _mol->end());
+    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();
+    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.");
+      ++Counter_rejections;
+      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 "
+        << zeroValue->second << ".");
+    LOG(4, "INFO: Squared norm of difference vector is "
+        << (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.");
+    LOG(1,"INFO: Resulting relative angle for molecule " << _mol->getName()
+        << " is " << angle << ".");
+    outmap->insert ( make_pair (angle, *iter ) );
+    ++i;
+  }
+  ASSERT(Counter_rejections <= molecules.size(),
+      "DipoleAngularCorrelation() - more rejections ("+toString(Counter_rejections)
+      +") than there are molecules ("+toString(molecules.size())+").");
+  LOG(1,"INFO: " << Counter_rejections << " molecules have been rejected in time step " << timestep << ".");
+
+  LOG(0,"STATUS: Done with calculating dipoles.");
+
+  if (DoTimeReset == DoResetTime) {
+    // re-set to original time step again
+    World::getInstance().setTime(oldtime);
+  }
+
+  // and return results
+  return outmap;
+};
+
+/** Calculates the dipole correlation for given molecule type.
+ * I.e. we calculate how the angle between any two given dipoles in the
+ * systems behaves. Sort of pair correlation but distance is replaced by
+ * the orientation distance, i.e. an angle.
  * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
  * Angles are given in degrees.
@@ -95 +240 @@
  * \return Map of doubles with values the pair of the two atoms.
  */
-DipoleAngularCorrelationMap *DipoleAngularCorrelation(std::vector<molecule *> &molecules)
+DipoleCorrelationMap *DipoleCorrelation(std::vector<molecule *> &molecules)
 {
   Info FunctionInfo(__func__);
-  DipoleAngularCorrelationMap *outmap = new DipoleAngularCorrelationMap;
+  DipoleCorrelationMap *outmap = new DipoleCorrelationMap;
 //  double distance = 0.;
 //  Box &domain = World::getInstance().getDomain();
@@ -107 +252 @@
   }
 
-  outmap = new DipoleAngularCorrelationMap;
   for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin();
       MolWalker != molecules.end(); ++MolWalker) {
@@ -168 +312 @@
   outmap = new PairCorrelationMap;
   for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++){
-    DoLog(2) && (Log()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
+    DoLog(2) && (Log()<< Verbose(2) << "Current molecule is " << (*MolWalker)->getName() << "." << endl);
     for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
       DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
       for (std::vector<molecule *>::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules.end(); MolOtherWalker++){
-        DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
+        DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << (*MolOtherWalker)->getName() << "." << endl);
         for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
           DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
@@ -503 +647 @@
 void OutputDipoleAngularCorrelation_Header( ofstream * const file )
 {
-  *file << "\tAtom1\tAtom2";
-};
-
-/** Prints values stored in DipoleAngularCorrelationMap iterator.
+  *file << "\tFirstAtomOfMolecule";
+};
+
+/** Prints values stored in DipoleCorrelationMap iterator.
  *
  * @param file stream to print to
@@ -512 +656 @@
  */
 void OutputDipoleAngularCorrelation_Value( ofstream * const file, DipoleAngularCorrelationMap::const_iterator &runner )
+{
+  *file << *(runner->second);
+};
+
+
+/** Adds header part that is unique to DipoleAngularCorrelationMap.
+ *
+ * @param file stream to print to
+ */
+void OutputDipoleCorrelation_Header( ofstream * const file )
+{
+  *file << "\tMolecule";
+};
+
+/** Prints values stored in DipoleCorrelationMap iterator.
+ *
+ * @param file stream to print to
+ * @param runner iterator pointing at values to print
+ */
+void OutputDipoleCorrelation_Value( ofstream * const file, DipoleCorrelationMap::const_iterator &runner )
 {
   *file << runner->second.first->getId() << "\t" << runner->second.second->getId();

src/Analysis/analysis_correlation.hpp

@@ -31 +31 @@
 #include "CodePatterns/Info.hpp"
 #include "CodePatterns/Log.hpp"
+#include "CodePatterns/Range.hpp"
 #include "CodePatterns/Verbose.hpp"
 #include "Helpers/helpers.hpp"
@@ -46 +47 @@
 
 typedef multimap<double, pair<atom *, atom *> > PairCorrelationMap;
-typedef multimap<double, pair<molecule *, molecule *> > DipoleAngularCorrelationMap;
+typedef multimap<double, atom * > DipoleAngularCorrelationMap;
+typedef multimap<double, pair<molecule *, molecule *> > DipoleCorrelationMap;
 typedef multimap<double, pair<atom *, const Vector *> > CorrelationToPointMap;
 typedef multimap<double, pair<atom *, BoundaryTriangleSet *> > CorrelationToSurfaceMap;
 typedef map<double, int> BinPairMap;
 
+enum ResetWorldTime {
+  DontResetTime,
+  DoResetTime
+};
+
 /********************************************** declarations *******************************/
 
-DipoleAngularCorrelationMap *DipoleAngularCorrelation(std::vector<molecule *> &molecules);
+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);
 CorrelationToPointMap *CorrelationToPoint(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Vector *point );
@@ -65 +76 @@
 void OutputDipoleAngularCorrelation_Header( ofstream * const file );
 void OutputDipoleAngularCorrelation_Value( ofstream * const file, DipoleAngularCorrelationMap::const_iterator &runner );
+void OutputDipoleCorrelation_Header( ofstream * const file );
+void OutputDipoleCorrelation_Value( ofstream * const file, DipoleCorrelationMap::const_iterator &runner );
 void OutputPairCorrelation_Header( ofstream * const file );
 void OutputPairCorrelation_Value( ofstream * const file, PairCorrelationMap::const_iterator &runner );