Changeset 6e5084


Ignore:
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.

Files:
8 edited

Legend:

Unmodified
Added
Removed

  • src/Actions/AnalysisAction/PrincipalAxisSystemAction.cpp

    ra439e5 r6e5084  
    1010#include "Actions/AnalysisAction/PrincipalAxisSystemAction.hpp"
    1111#include "Actions/ActionRegistry.hpp"
    12 #include "molecule.hpp"
    1312#include "Helpers/Log.hpp"
    1413#include "Helpers/Verbose.hpp"
     14#include "LinearAlgebra/Matrix.hpp"
     15#include "LinearAlgebra/Vector.hpp"
     16#include "element.hpp"
     17#include "molecule.hpp"
    1518
    1619#include <iostream>
    … …  
    4649Action::state_ptr AnalysisPrincipalAxisSystemAction::performCall() {
    4750  molecule *mol = NULL;
     51  Matrix InertiaTensor;
    4852
    4953  ValueStorage::getInstance().queryCurrentValue(NAME, mol);
    … …  
    5155  for (World::MoleculeSelectionIterator iter = World::getInstance().beginMoleculeSelection(); iter != World::getInstance().endMoleculeSelection(); ++iter) {
    5256    molecule *mol = iter->second;
    53     mol->PrincipalAxisSystem(false);
     57    Vector *CenterOfGravity = mol->DetermineCenterOfGravity();
     58
     59    // reset inertia tensor
     60    InertiaTensor.zero();
     61
     62    // sum up inertia tensor
     63    for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     64      Vector x = (*iter)->x;
     65      x -= *CenterOfGravity;
     66      InertiaTensor.at(0,0) += (*iter)->type->mass*(x[1]*x[1] + x[2]*x[2]);
     67      InertiaTensor.at(0,1) += (*iter)->type->mass*(-x[0]*x[1]);
     68      InertiaTensor.at(0,2) += (*iter)->type->mass*(-x[0]*x[2]);
     69      InertiaTensor.at(1,0) += (*iter)->type->mass*(-x[1]*x[0]);
     70      InertiaTensor.at(1,1) += (*iter)->type->mass*(x[0]*x[0] + x[2]*x[2]);
     71      InertiaTensor.at(1,2) += (*iter)->type->mass*(-x[1]*x[2]);
     72      InertiaTensor.at(2,0) += (*iter)->type->mass*(-x[2]*x[0]);
     73      InertiaTensor.at(2,1) += (*iter)->type->mass*(-x[2]*x[1]);
     74      InertiaTensor.at(2,2) += (*iter)->type->mass*(x[0]*x[0] + x[1]*x[1]);
     75    }
     76    // print InertiaTensor for debugging
     77    DoLog(0) && (Log() << Verbose(0) << "The inertia tensor is:" << InertiaTensor << endl);
    5478  }
    5579  return Action::success;

  • src/Actions/MoleculeAction/RotateToPrincipalAxisSystemAction.cpp

    ra439e5 r6e5084  
    1111#include "Actions/ActionRegistry.hpp"
    1212#include "Helpers/Log.hpp"
     13#include "Helpers/Verbose.hpp"
     14#include "LinearAlgebra/Line.hpp"
     15#include "LinearAlgebra/Matrix.hpp"
     16#include "LinearAlgebra/Vector.hpp"
     17#include "element.hpp"
    1318#include "molecule.hpp"
    14 #include "Helpers/Verbose.hpp"
    1519
    1620
    … …  
    4852{}
    4953
    50 void MoleculeRotateToPrincipalAxisSystem() {
     54void MoleculeRotateToPrincipalAxisSystem(Vector &Axis) {
     55  ValueStorage::getInstance().setCurrentValue(MoleculeRotateToPrincipalAxisSystemAction::NAME, Axis);
    5156  ActionRegistry::getInstance().getActionByName(MoleculeRotateToPrincipalAxisSystemAction::NAME)->call(Action::NonInteractive);
    5257};
    … …  
    5560  ASSERT(dialog,"No Dialog given when filling action dialog");
    5661
    57   dialog->queryEmpty(NAME, MapOfActions::getInstance().getDescription(NAME));
     62  dialog->queryVector(NAME, false, MapOfActions::getInstance().getDescription(NAME));
    5863
    5964  return dialog;
    … …  
    6267Action::state_ptr MoleculeRotateToPrincipalAxisSystemAction::performCall() {
    6368  molecule *mol = NULL;
     69  Vector Axis;
     70
     71  // obtain axis to rotate to
     72  ValueStorage::getInstance().queryCurrentValue(NAME, Axis);
    6473
    6574  for (World::MoleculeSelectionIterator iter = World::getInstance().beginMoleculeSelection(); iter != World::getInstance().endMoleculeSelection(); ++iter) {
    6675    mol = iter->second;
    6776    DoLog(0) && (Log() << Verbose(0) << "Converting to prinicipal axis system." << endl);
    68     mol->PrincipalAxisSystem(true);
     77    Matrix InertiaTensor;
     78    Vector *CenterOfGravity = mol->DetermineCenterOfGravity();
     79
     80    // reset inertia tensor
     81    InertiaTensor.zero();
     82
     83    // sum up inertia tensor
     84    for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     85      Vector x = (*iter)->x;
     86      x -= *CenterOfGravity;
     87      InertiaTensor.at(0,0) += (*iter)->type->mass*(x[1]*x[1] + x[2]*x[2]);
     88      InertiaTensor.at(0,1) += (*iter)->type->mass*(-x[0]*x[1]);
     89      InertiaTensor.at(0,2) += (*iter)->type->mass*(-x[0]*x[2]);
     90      InertiaTensor.at(1,0) += (*iter)->type->mass*(-x[1]*x[0]);
     91      InertiaTensor.at(1,1) += (*iter)->type->mass*(x[0]*x[0] + x[2]*x[2]);
     92      InertiaTensor.at(1,2) += (*iter)->type->mass*(-x[1]*x[2]);
     93      InertiaTensor.at(2,0) += (*iter)->type->mass*(-x[2]*x[0]);
     94      InertiaTensor.at(2,1) += (*iter)->type->mass*(-x[2]*x[1]);
     95      InertiaTensor.at(2,2) += (*iter)->type->mass*(x[0]*x[0] + x[1]*x[1]);
     96    }
     97    // print InertiaTensor for debugging
     98    DoLog(0) && (Log() << Verbose(0) << "The inertia tensor is:" << InertiaTensor << endl);
     99
     100    // diagonalize to determine principal axis system
     101    Vector Eigenvalues = InertiaTensor.transformToEigenbasis();
     102
     103    for(int i=0;i<NDIM;i++)
     104      DoLog(0) && (Log() << Verbose(0) << "eigenvalue = " << Eigenvalues[i] << ", eigenvector = " << InertiaTensor.column(i) << endl);
     105
     106    // check whether we rotate or not
     107    DoLog(0) && (Log() << Verbose(0) << "Transforming molecule into PAS ... ");
     108
     109    // obtain first column, eigenvector to biggest eigenvalue
     110    Vector BiggestEigenvector(InertiaTensor.column(Eigenvalues.SmallestComponent()));
     111    Vector DesiredAxis(Axis);
     112
     113    // Creation Line that is the rotation axis
     114    DesiredAxis.VectorProduct(BiggestEigenvector);
     115    Line RotationAxis(Vector(0.,0.,0.), DesiredAxis);
     116
     117    // determine angle
     118    const double alpha = BiggestEigenvector.Angle(Axis);
     119
     120    DoLog(0) && (Log() << Verbose(0) << alpha << endl);
     121
     122    for (molecule::iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     123      *((*iter)->node) -= *CenterOfGravity;
     124      *((*iter)->node) = RotationAxis.rotateVector(*((*iter)->node), alpha);
     125      *((*iter)->node) += *CenterOfGravity;
     126    }
     127    DoLog(0) && (Log() << Verbose(0) << "done." << endl);
     128
     129    // summing anew for debugging (resulting matrix has to be diagonal!)
     130    // reset inertia tensor
     131    InertiaTensor.zero();
     132
     133    // sum up inertia tensor
     134    for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     135      Vector x = (*iter)->x;
     136      x -= *CenterOfGravity;
     137      InertiaTensor.at(0,0) += (*iter)->type->mass*(x[1]*x[1] + x[2]*x[2]);
     138      InertiaTensor.at(0,1) += (*iter)->type->mass*(-x[0]*x[1]);
     139      InertiaTensor.at(0,2) += (*iter)->type->mass*(-x[0]*x[2]);
     140      InertiaTensor.at(1,0) += (*iter)->type->mass*(-x[1]*x[0]);
     141      InertiaTensor.at(1,1) += (*iter)->type->mass*(x[0]*x[0] + x[2]*x[2]);
     142      InertiaTensor.at(1,2) += (*iter)->type->mass*(-x[1]*x[2]);
     143      InertiaTensor.at(2,0) += (*iter)->type->mass*(-x[2]*x[0]);
     144      InertiaTensor.at(2,1) += (*iter)->type->mass*(-x[2]*x[1]);
     145      InertiaTensor.at(2,2) += (*iter)->type->mass*(x[0]*x[0] + x[1]*x[1]);
     146      // print InertiaTensor for debugging
     147      DoLog(0) && (Log() << Verbose(0) << "The inertia tensor is:" << InertiaTensor << endl);
     148    }
     149
     150    // free everything
     151    delete(CenterOfGravity);
    69152  }
    70153  return Action::success;

  • src/Actions/MoleculeAction/RotateToPrincipalAxisSystemAction.hpp

    ra439e5 r6e5084  
    1313
    1414class MoleculeListClass;
     15class Vector;
    1516
    16 void MoleculeRotateToPrincipalAxisSystem();
     17void MoleculeRotateToPrincipalAxisSystem(Vector &Axis);
    1718
    1819class MoleculeRotateToPrincipalAxisSystemAction : public Action {
    19   friend void MoleculeRotateToPrincipalAxisSystem();
     20  friend void MoleculeRotateToPrincipalAxisSystem(Vector &Axis);
    2021
    2122public:

  • src/Makefile.am

    ra439e5 r6e5084  
    303303joiner_SOURCES = joiner.cpp datacreator.cpp datacreator.hpp periodentafel.hpp
    304304joiner_LDADD = \
     305        Actions/libMolecuilderActions-@MOLECUILDER_API_VERSION@.la \
    305306        libMolecuilder-@MOLECUILDER_API_VERSION@.la \
     307        Parser/libMolecuilderParser-@MOLECUILDER_API_VERSION@.la \
    306308        LinearAlgebra/libMolecuilderLinearAlgebra-@MOLECUILDER_API_VERSION@.la \
    307309        Exceptions/libMolecuilderExceptions-@MOLECUILDER_API_VERSION@.la \
    … …  
    313315analyzer_SOURCES = analyzer.cpp datacreator.cpp periodentafel.hpp datacreator.hpp
    314316analyzer_LDADD = \
     317        Actions/libMolecuilderActions-@MOLECUILDER_API_VERSION@.la \
    315318        libMolecuilder-@MOLECUILDER_API_VERSION@.la \
     319        Parser/libMolecuilderParser-@MOLECUILDER_API_VERSION@.la \
    316320        LinearAlgebra/libMolecuilderLinearAlgebra-@MOLECUILDER_API_VERSION@.la \
    317321        Exceptions/libMolecuilderExceptions-@MOLECUILDER_API_VERSION@.la \

  • src/boundary.cpp

    ra439e5 r6e5084  
    66#include "Helpers/MemDebug.hpp"
    77
     8#include "Actions/MoleculeAction/RotateToPrincipalAxisSystemAction.hpp"
    89#include "World.hpp"
    910#include "atom.hpp"
    … …  
    674675  double cellvolume = 0.;
    675676
    676   // transform to PAS
    677   mol->PrincipalAxisSystem(true);
     677  // transform to PAS by Action
     678  Vector MainAxis(0.,0.,1.);
     679  MoleculeRotateToPrincipalAxisSystem(MainAxis);
    678680
    679681  IsAngstroem = configuration->GetIsAngstroem();

  • src/molecule.hpp

    ra439e5 r6e5084  
    278278  bool VerletForceIntegration(char *file, config &configuration, const size_t offset);
    279279  void Thermostats(config &configuration, double ActualTemp, int Thermostat);
    280   void PrincipalAxisSystem(bool DoRotate);
    281280  double VolumeOfConvexEnvelope(bool IsAngstroem);
    282281

  • src/molecule_geometry.cpp

    ra439e5 r6e5084  
    1010#endif
    1111
     12#include "Helpers/helpers.hpp"
     13#include "Helpers/Log.hpp"
    1214#include "Helpers/MemDebug.hpp"
     15#include "Helpers/Verbose.hpp"
     16#include "LinearAlgebra/Line.hpp"
     17#include "LinearAlgebra/Matrix.hpp"
     18#include "LinearAlgebra/Plane.hpp"
    1319
    1420#include "atom.hpp"
    … …  
    1622#include "config.hpp"
    1723#include "element.hpp"
    18 #include "Helpers/helpers.hpp"
    1924#include "leastsquaremin.hpp"
    20 #include "Helpers/Verbose.hpp"
    21 #include "Helpers/Log.hpp"
    2225#include "molecule.hpp"
    2326#include "World.hpp"
    24 #include "LinearAlgebra/Plane.hpp"
    25 #include "LinearAlgebra/Matrix.hpp"
     27
    2628#include "Box.hpp"
     29
    2730#include <boost/foreach.hpp>
    2831
    … …  
    316319};
    317320
    318 /** Transforms/Rotates the given molecule into its principal axis system.
    319  * \param *out output stream for debugging
    320  * \param DoRotate whether to rotate (true) or only to determine the PAS.
    321  * TODO treatment of trajetories missing
    322  */
    323 void molecule::PrincipalAxisSystem(bool DoRotate)
    324 {
    325   double InertiaTensor[NDIM*NDIM];
    326   Vector *CenterOfGravity = DetermineCenterOfGravity();
    327 
    328   CenterPeriodic();
    329 
    330   // reset inertia tensor
    331   for(int i=0;i<NDIM*NDIM;i++)
    332     InertiaTensor[i] = 0.;
    333 
    334   // sum up inertia tensor
    335   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    336     Vector x = (*iter)->x;
    337     //x.SubtractVector(CenterOfGravity);
    338     InertiaTensor[0] += (*iter)->type->mass*(x[1]*x[1] + x[2]*x[2]);
    339     InertiaTensor[1] += (*iter)->type->mass*(-x[0]*x[1]);
    340     InertiaTensor[2] += (*iter)->type->mass*(-x[0]*x[2]);
    341     InertiaTensor[3] += (*iter)->type->mass*(-x[1]*x[0]);
    342     InertiaTensor[4] += (*iter)->type->mass*(x[0]*x[0] + x[2]*x[2]);
    343     InertiaTensor[5] += (*iter)->type->mass*(-x[1]*x[2]);
    344     InertiaTensor[6] += (*iter)->type->mass*(-x[2]*x[0]);
    345     InertiaTensor[7] += (*iter)->type->mass*(-x[2]*x[1]);
    346     InertiaTensor[8] += (*iter)->type->mass*(x[0]*x[0] + x[1]*x[1]);
    347   }
    348   // print InertiaTensor for debugging
    349   DoLog(0) && (Log() << Verbose(0) << "The inertia tensor is:" << endl);
    350   for(int i=0;i<NDIM;i++) {
    351     for(int j=0;j<NDIM;j++)
    352       DoLog(0) && (Log() << Verbose(0) << InertiaTensor[i*NDIM+j] << " ");
    353     DoLog(0) && (Log() << Verbose(0) << endl);
    354   }
    355   DoLog(0) && (Log() << Verbose(0) << endl);
    356 
    357   // diagonalize to determine principal axis system
    358   gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(NDIM);
    359   gsl_matrix_view m = gsl_matrix_view_array(InertiaTensor, NDIM, NDIM);
    360   gsl_vector *eval = gsl_vector_alloc(NDIM);
    361   gsl_matrix *evec = gsl_matrix_alloc(NDIM, NDIM);
    362   gsl_eigen_symmv(&m.matrix, eval, evec, T);
    363   gsl_eigen_symmv_free(T);
    364   gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC);
    365 
    366   for(int i=0;i<NDIM;i++) {
    367     DoLog(1) && (Log() << Verbose(1) << "eigenvalue = " << gsl_vector_get(eval, i));
    368     DoLog(0) && (Log() << Verbose(0) << ", eigenvector = (" << evec->data[i * evec->tda + 0] << "," << evec->data[i * evec->tda + 1] << "," << evec->data[i * evec->tda + 2] << ")" << endl);
    369   }
    370 
    371   // check whether we rotate or not
    372   if (DoRotate) {
    373     DoLog(1) && (Log() << Verbose(1) << "Transforming molecule into PAS ... ");
    374     // the eigenvectors specify the transformation matrix
    375     Matrix M = Matrix(evec->data);
    376 
    377     BOOST_FOREACH(atom* iter, atoms){
    378       (*iter->node) *= M;
    379     }
    380     DoLog(0) && (Log() << Verbose(0) << "done." << endl);
    381 
    382     // summing anew for debugging (resulting matrix has to be diagonal!)
    383     // reset inertia tensor
    384     for(int i=0;i<NDIM*NDIM;i++)
    385       InertiaTensor[i] = 0.;
    386 
    387     // sum up inertia tensor
    388     for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    389       Vector x = (*iter)->x;
    390       InertiaTensor[0] += (*iter)->type->mass*(x[1]*x[1] + x[2]*x[2]);
    391       InertiaTensor[1] += (*iter)->type->mass*(-x[0]*x[1]);
    392       InertiaTensor[2] += (*iter)->type->mass*(-x[0]*x[2]);
    393       InertiaTensor[3] += (*iter)->type->mass*(-x[1]*x[0]);
    394       InertiaTensor[4] += (*iter)->type->mass*(x[0]*x[0] + x[2]*x[2]);
    395       InertiaTensor[5] += (*iter)->type->mass*(-x[1]*x[2]);
    396       InertiaTensor[6] += (*iter)->type->mass*(-x[2]*x[0]);
    397       InertiaTensor[7] += (*iter)->type->mass*(-x[2]*x[1]);
    398       InertiaTensor[8] += (*iter)->type->mass*(x[0]*x[0] + x[1]*x[1]);
    399     }
    400     // print InertiaTensor for debugging
    401     DoLog(0) && (Log() << Verbose(0) << "The inertia tensor is:" << endl);
    402     for(int i=0;i<NDIM;i++) {
    403       for(int j=0;j<NDIM;j++)
    404         DoLog(0) && (Log() << Verbose(0) << InertiaTensor[i*NDIM+j] << " ");
    405       DoLog(0) && (Log() << Verbose(0) << endl);
    406     }
    407     DoLog(0) && (Log() << Verbose(0) << endl);
    408   }
    409 
    410   // free everything
    411   delete(CenterOfGravity);
    412   gsl_vector_free(eval);
    413   gsl_matrix_free(evec);
    414 };
    415 
    416 
    417321/** Align all atoms in such a manner that given vector \a *n is along z axis.
    418322 * \param n[] alignment vector.

  • tests/regression/testsuite-molecules.at

    ra439e5 r6e5084  
    6464AT_KEYWORDS([Molecules])
    6565AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/Molecules/8/pre/test.* .], 0)
    66 AT_CHECK([../../molecuilder -i test.conf -e ${abs_top_srcdir}/src/ --select-molecule-by-id 0 -m 0], 0, [stdout], [stderr])
     66AT_CHECK([../../molecuilder -i test.conf -e ${abs_top_srcdir}/src/ --select-molecule-by-id 0 -m "0,0,1"], 0, [stdout], [stderr])
    6767#AT_CHECK([file=test.conf; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/Molecules/8/post/$file], 0, [ignore], [ignore])
    6868AT_CLEANUP
Note: See TracChangeset for help on using the changeset viewer.