Changeset 6e5084 for src/Actions/MoleculeAction/RotateToPrincipalAxisSystemAction.cpp

Timestamp:

Aug 6, 2010, 1:32:00 PM (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:

eaf4ae

Parents:

a439e5

git-author:

Frederik Heber <heber@…> (08/05/10 14:34:18)

git-committer:

Frederik Heber <heber@…> (08/06/10 13:32:00)

Message:

Fixed PrincipalAxisSystemAction.

• PrincipalAxisSystemAction:performCall() now has initial parts of former molecule::PrincipalAxisSystem()
• RotateToPrincipalAxisSystemAction::performCall () now has former molecule::PrincipalAxisSystem()
• call in PrepareClustersinWater() replaced by call to this Action
• changes in PrincipalAxisSystemAction/RotateToPrincipalAxisSystemAction:
  • use Matrix class for InertiaTensor
  • use Vector class for axis
  • use Line class for rotation axis

File:

• src/Actions/MoleculeAction/RotateToPrincipalAxisSystemAction.cpp (modified) (4 diffs)

src/Actions/MoleculeAction/RotateToPrincipalAxisSystemAction.cpp

@@ -11 +11 @@
 #include "Actions/ActionRegistry.hpp"
 #include "Helpers/Log.hpp"
+#include "Helpers/Verbose.hpp"
+#include "LinearAlgebra/Line.hpp"
+#include "LinearAlgebra/Matrix.hpp"
+#include "LinearAlgebra/Vector.hpp"
+#include "element.hpp"
 #include "molecule.hpp"
-#include "Helpers/Verbose.hpp"
 
 
@@ -48 +52 @@
 {}
 
-void MoleculeRotateToPrincipalAxisSystem() {
+void MoleculeRotateToPrincipalAxisSystem(Vector &Axis) {
+  ValueStorage::getInstance().setCurrentValue(MoleculeRotateToPrincipalAxisSystemAction::NAME, Axis);
   ActionRegistry::getInstance().getActionByName(MoleculeRotateToPrincipalAxisSystemAction::NAME)->call(Action::NonInteractive);
 };
@@ -55 +60 @@
   ASSERT(dialog,"No Dialog given when filling action dialog");
 
-  dialog->queryEmpty(NAME, MapOfActions::getInstance().getDescription(NAME));
+  dialog->queryVector(NAME, false, MapOfActions::getInstance().getDescription(NAME));
 
   return dialog;
@@ -62 +67 @@
 Action::state_ptr MoleculeRotateToPrincipalAxisSystemAction::performCall() {
   molecule *mol = NULL;
+  Vector Axis;
+
+  // obtain axis to rotate to
+  ValueStorage::getInstance().queryCurrentValue(NAME, Axis);
 
   for (World::MoleculeSelectionIterator iter = World::getInstance().beginMoleculeSelection(); iter != World::getInstance().endMoleculeSelection(); ++iter) {
     mol = iter->second;
     DoLog(0) && (Log() << Verbose(0) << "Converting to prinicipal axis system." << endl);
-    mol->PrincipalAxisSystem(true);
+    Matrix InertiaTensor;
+    Vector *CenterOfGravity = mol->DetermineCenterOfGravity();
+
+    // reset inertia tensor
+    InertiaTensor.zero();
+
+    // sum up inertia tensor
+    for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+      Vector x = (*iter)->x;
+      x -= *CenterOfGravity;
+      InertiaTensor.at(0,0) += (*iter)->type->mass*(x[1]*x[1] + x[2]*x[2]);
+      InertiaTensor.at(0,1) += (*iter)->type->mass*(-x[0]*x[1]);
+      InertiaTensor.at(0,2) += (*iter)->type->mass*(-x[0]*x[2]);
+      InertiaTensor.at(1,0) += (*iter)->type->mass*(-x[1]*x[0]);
+      InertiaTensor.at(1,1) += (*iter)->type->mass*(x[0]*x[0] + x[2]*x[2]);
+      InertiaTensor.at(1,2) += (*iter)->type->mass*(-x[1]*x[2]);
+      InertiaTensor.at(2,0) += (*iter)->type->mass*(-x[2]*x[0]);
+      InertiaTensor.at(2,1) += (*iter)->type->mass*(-x[2]*x[1]);
+      InertiaTensor.at(2,2) += (*iter)->type->mass*(x[0]*x[0] + x[1]*x[1]);
+    }
+    // print InertiaTensor for debugging
+    DoLog(0) && (Log() << Verbose(0) << "The inertia tensor is:" << InertiaTensor << endl);
+
+    // diagonalize to determine principal axis system
+    Vector Eigenvalues = InertiaTensor.transformToEigenbasis();
+
+    for(int i=0;i<NDIM;i++)
+      DoLog(0) && (Log() << Verbose(0) << "eigenvalue = " << Eigenvalues[i] << ", eigenvector = " << InertiaTensor.column(i) << endl);
+
+    // check whether we rotate or not
+    DoLog(0) && (Log() << Verbose(0) << "Transforming molecule into PAS ... ");
+
+    // obtain first column, eigenvector to biggest eigenvalue
+    Vector BiggestEigenvector(InertiaTensor.column(Eigenvalues.SmallestComponent()));
+    Vector DesiredAxis(Axis);
+
+    // Creation Line that is the rotation axis
+    DesiredAxis.VectorProduct(BiggestEigenvector);
+    Line RotationAxis(Vector(0.,0.,0.), DesiredAxis);
+
+    // determine angle
+    const double alpha = BiggestEigenvector.Angle(Axis);
+
+    DoLog(0) && (Log() << Verbose(0) << alpha << endl);
+
+    for (molecule::iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+      *((*iter)->node) -= *CenterOfGravity;
+      *((*iter)->node) = RotationAxis.rotateVector(*((*iter)->node), alpha);
+      *((*iter)->node) += *CenterOfGravity;
+    }
+    DoLog(0) && (Log() << Verbose(0) << "done." << endl);
+
+    // summing anew for debugging (resulting matrix has to be diagonal!)
+    // reset inertia tensor
+    InertiaTensor.zero();
+
+    // sum up inertia tensor
+    for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
+      Vector x = (*iter)->x;
+      x -= *CenterOfGravity;
+      InertiaTensor.at(0,0) += (*iter)->type->mass*(x[1]*x[1] + x[2]*x[2]);
+      InertiaTensor.at(0,1) += (*iter)->type->mass*(-x[0]*x[1]);
+      InertiaTensor.at(0,2) += (*iter)->type->mass*(-x[0]*x[2]);
+      InertiaTensor.at(1,0) += (*iter)->type->mass*(-x[1]*x[0]);
+      InertiaTensor.at(1,1) += (*iter)->type->mass*(x[0]*x[0] + x[2]*x[2]);
+      InertiaTensor.at(1,2) += (*iter)->type->mass*(-x[1]*x[2]);
+      InertiaTensor.at(2,0) += (*iter)->type->mass*(-x[2]*x[0]);
+      InertiaTensor.at(2,1) += (*iter)->type->mass*(-x[2]*x[1]);
+      InertiaTensor.at(2,2) += (*iter)->type->mass*(x[0]*x[0] + x[1]*x[1]);
+      // print InertiaTensor for debugging
+      DoLog(0) && (Log() << Verbose(0) << "The inertia tensor is:" << InertiaTensor << endl);
+    }
+
+    // free everything
+    delete(CenterOfGravity);
   }
   return Action::success;