Changeset be945c

Timestamp:

Feb 3, 2011, 9:51:18 AM (14 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:

caa30b

Parents:

67c92b

git-author:

Frederik Heber <heber@…> (12/17/10 18:52:44)

git-committer:

Frederik Heber <heber@…> (02/03/11 09:51:18)

Message:

New DipoleAngularCorrelationAction that calculates angular correlation between dipoles.

Location:

src

Files:

• Actions/AnalysisAction/DipoleAngularCorrelationAction.cpp (added)
• Actions/AnalysisAction/DipoleAngularCorrelationAction.def (added)
• Actions/AnalysisAction/DipoleAngularCorrelationAction.hpp (added)
• Actions/Makefile.am (modified) (2 diffs)
• analysis_correlation.cpp (modified) (5 diffs)
• analysis_correlation.hpp (modified) (2 diffs)

src/Actions/Makefile.am ¶

@@ -37 +37 @@
   
 ANALYSISACTIONSOURCE = \
+  AnalysisAction/DipoleAngularCorrelationAction.cpp \
   AnalysisAction/MolecularVolumeAction.cpp \
   AnalysisAction/PairCorrelationAction.cpp \
@@ -43 +44 @@
   AnalysisAction/SurfaceCorrelationAction.cpp
 ANALYSISACTIONHEADER = \
+  AnalysisAction/DipoleAngularCorrelationAction.hpp \
   AnalysisAction/MolecularVolumeAction.hpp \
   AnalysisAction/PairCorrelationAction.hpp \

src/analysis_correlation.cpp ¶

@@ -23 +23 @@
 #include <iomanip>
 
+#include "atom.hpp"
+#include "bond.hpp"
 #include "BoundaryTriangleSet.hpp"
-#include "analysis_correlation.hpp"
+#include "Box.hpp"
 #include "element.hpp"
 #include "CodePatterns/Info.hpp"
@@ -33 +35 @@
 #include "tesselationhelpers.hpp"
 #include "triangleintersectionlist.hpp"
+#include "World.hpp"
 #include "LinearAlgebra/Vector.hpp"
 #include "LinearAlgebra/RealSpaceMatrix.hpp"
@@ -39 +42 @@
 #include "Box.hpp"
 
+#include "analysis_correlation.hpp"
+
+/** Calculates the dipole vector of a given atomSet.
+ *
+ *  Note that we use the following procedure as rule of thumb:
+ *   -# go through every bond of the atom
+ *   -# calculate the difference of electronegativities \f$\Delta\text{EN}\f$
+ *   -# if \f$\Delta\text{EN} > 0.5\f$, we align the bond vector in direction of the more negative element
+ *   -# sum up all vectors
+ *   -# finally, divide by the number of summed vectors
+ *
+ * @param atomsbegin begin iterator of atomSet
+ * @param atomsend end iterator of atomset
+ * @return dipole vector
+ */
+Vector getDipole(molecule::const_iterator atomsbegin, molecule::const_iterator atomsend)
+{
+  Vector DipoleVector;
+  size_t SumOfVectors = 0;
+  // go through all atoms
+  for (molecule::const_iterator atomiter = atomsbegin;
+      atomiter != atomsend;
+      ++atomiter) {
+    // go through all bonds
+    for (BondList::const_iterator bonditer = (*atomiter)->ListOfBonds.begin();
+        bonditer != (*atomiter)->ListOfBonds.end();
+        ++bonditer) {
+      const atom * Otheratom = (*bonditer)->GetOtherAtom(*atomiter);
+      if (Otheratom->getId() > (*atomiter)->getId()) {
+        const double DeltaEN = (*atomiter)->getType()->getElectronegativity()
+            -Otheratom->getType()->getElectronegativity();
+        Vector BondDipoleVector = (*atomiter)->getPosition() - Otheratom->getPosition();
+        // DeltaEN is always positive, gives correct orientation of vector
+        BondDipoleVector.Normalize();
+        BondDipoleVector *= DeltaEN;
+        DipoleVector += BondDipoleVector;
+        SumOfVectors++;
+      }
+    }
+  }
+  DipoleVector *= 1./(double)SumOfVectors;
+  DoLog(1) && (Log() << Verbose(1) << "Resulting dipole vector is " << DipoleVector << std::endl);
+
+  return DipoleVector;
+};
+
 /** Calculates the dipole angular correlation for given molecule type.
  * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
+ * Angles are given in degrees.
  * \param *molecules vector of molecules
- * \param &elements vector of elements to correlate
  * \return Map of doubles with values the pair of the two atoms.
  */
-DipoleAngularCorrelationMap *DipoleAngularCorrelation(std::vector<molecule *> &molecules, const Formula &formula)
+DipoleAngularCorrelationMap *DipoleAngularCorrelation(std::vector<molecule *> &molecules)
 {
   Info FunctionInfo(__func__);
@@ -52 +101 @@
 //  Box &domain = World::getInstance().getDomain();
 //
-//  if (molecules.empty()) {
-//    DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
-//    return outmap;
-//  }
-//  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
-//    (*MolWalker)->doCountAtoms();
-//
-//  // create all possible pairs of elements
-//  set <pair<const element *,const element *> > PairsOfElements;
-//  if (elements.size() >= 2) {
-//    for (vector<const element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
-//      for (vector<const element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
-//        if (type1 != type2) {
-//          PairsOfElements.insert( make_pair(*type1,*type2) );
-//          DoLog(1) && (Log() << Verbose(1) << "Creating element pair " << *(*type1) << " and " << *(*type2) << "." << endl);
-//        }
-//  } else if (elements.size() == 1) { // one to all are valid
-//    const element *elemental = *elements.begin();
-//    PairsOfElements.insert( pair<const element *,const element*>(elemental,0) );
-//    PairsOfElements.insert( pair<const element *,const element*>(0,elemental) );
-//  } else { // all elements valid
-//    PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
-//  }
-//
-//  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);
-//    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);
-//        for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
-//          DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
-//          if ((*iter)->getId() < (*runner)->getId()){
-//            for (set <pair<const element *, const element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
-//              if ((PairRunner->first == (**iter).getType()) && (PairRunner->second == (**runner).getType())) {
-//                distance = domain.periodicDistance((*iter)->getPosition(),(*runner)->getPosition());
-//                //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
-//                outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
-//              }
-//          }
-//        }
-//      }
-//    }
-//  }
+  if (molecules.empty()) {
+    DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
+    return outmap;
+  }
+
+  outmap = new DipoleAngularCorrelationMap;
+  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin();
+      MolWalker != molecules.end();) {
+    DoLog(2) && (Log()<< Verbose(2) << "Current molecule is "
+        << (*MolWalker)->getId() << "." << endl);
+    const Vector Dipole = getDipole((*MolWalker)->begin(), (*MolWalker)->end());
+    for (std::vector<molecule *>::const_iterator MolOtherWalker = ++MolWalker;
+        MolOtherWalker != molecules.end();
+        MolOtherWalker++) {
+      DoLog(2) && (Log() << Verbose(2) << "Current other molecule is "
+          << (*MolOtherWalker)->getId() << "." << endl);
+      const Vector OtherDipole = getDipole((*MolOtherWalker)->begin(), (*MolOtherWalker)->end());
+      const double angle = Dipole.Angle(OtherDipole) * (180./M_PI);
+      DoLog(1) && (Log() << Verbose(1) << "Angle is " << angle << "." << endl);
+      outmap->insert ( make_pair (angle, make_pair ((*MolWalker), (*MolOtherWalker)) ) );
+    }
+  }
   return outmap;
 };
@@ -466 +492 @@
  * \param *map map to write
  */
+void OutputDipoleAngularCorrelation( ofstream * const file, const DipoleAngularCorrelationMap * const map )
+{
+  Info FunctionInfo(__func__);
+  *file << "BinStart\tMolecule1\tMolecule2" << endl;
+  for (DipoleAngularCorrelationMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
+    *file << setprecision(8) << runner->first << "\t" << runner->second.first->getId() << "\t" << runner->second.second->getId() << endl;
+  }
+};
+
+/** Prints correlation (double, (atom*,atom*) ) pairs to file.
+ * \param *file file to write to
+ * \param *map map to write
+ */
 void OutputPairCorrelation( ofstream * const file, const PairCorrelationMap * const map )
 {

src/analysis_correlation.hpp ¶

@@ -49 +49 @@
 /********************************************** declarations *******************************/
 
-DipoleAngularCorrelationMap *DipoleAngularCorrelation(std::vector<molecule *> &molecules, const Formula &formula);
+DipoleAngularCorrelationMap *DipoleAngularCorrelation(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 );
@@ -58 +58 @@
 int GetBin ( const double value, const double BinWidth, const double BinStart );
 void OutputCorrelation( ofstream * const file, const BinPairMap * const map );
+void OutputDipoleAngularCorrelation( ofstream * const file, const DipoleAngularCorrelationMap * const map );
 void OutputPairCorrelation( ofstream * const file, const PairCorrelationMap * const map );
 void OutputCorrelationToPoint( ofstream * const file, const CorrelationToPointMap * const map );