Index: src/Actions/AnalysisAction/PrincipalAxisSystemAction.cpp
===================================================================
--- src/Actions/AnalysisAction/PrincipalAxisSystemAction.cpp	(revision d223d5e07c484fa4ef7f54e88e97fa78b5c0b504)
+++ src/Actions/AnalysisAction/PrincipalAxisSystemAction.cpp	(revision 6e5084ba8942d5b07318fb32c642c0078c9dafe5)
@@ -10,7 +10,10 @@
 #include "Actions/AnalysisAction/PrincipalAxisSystemAction.hpp"
 #include "Actions/ActionRegistry.hpp"
-#include "molecule.hpp"
 #include "Helpers/Log.hpp"
 #include "Helpers/Verbose.hpp"
+#include "LinearAlgebra/Matrix.hpp"
+#include "LinearAlgebra/Vector.hpp"
+#include "element.hpp"
+#include "molecule.hpp"
 
 #include <iostream>
@@ -46,4 +49,5 @@
 Action::state_ptr AnalysisPrincipalAxisSystemAction::performCall() {
   molecule *mol = NULL;
+  Matrix InertiaTensor;
 
   ValueStorage::getInstance().queryCurrentValue(NAME, mol);
@@ -51,5 +55,25 @@
   for (World::MoleculeSelectionIterator iter = World::getInstance().beginMoleculeSelection(); iter != World::getInstance().endMoleculeSelection(); ++iter) {
     molecule *mol = iter->second;
-    mol->PrincipalAxisSystem(false);
+    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);
   }
   return Action::success;
Index: src/Actions/MoleculeAction/RotateToPrincipalAxisSystemAction.cpp
===================================================================
--- src/Actions/MoleculeAction/RotateToPrincipalAxisSystemAction.cpp	(revision d223d5e07c484fa4ef7f54e88e97fa78b5c0b504)
+++ src/Actions/MoleculeAction/RotateToPrincipalAxisSystemAction.cpp	(revision 6e5084ba8942d5b07318fb32c642c0078c9dafe5)
@@ -11,6 +11,10 @@
 #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,5 +52,6 @@
 {}
 
-void MoleculeRotateToPrincipalAxisSystem() {
+void MoleculeRotateToPrincipalAxisSystem(Vector &Axis) {
+  ValueStorage::getInstance().setCurrentValue(MoleculeRotateToPrincipalAxisSystemAction::NAME, Axis);
   ActionRegistry::getInstance().getActionByName(MoleculeRotateToPrincipalAxisSystemAction::NAME)->call(Action::NonInteractive);
 };
@@ -55,5 +60,5 @@
   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,9 +67,87 @@
 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;
Index: src/Actions/MoleculeAction/RotateToPrincipalAxisSystemAction.hpp
===================================================================
--- src/Actions/MoleculeAction/RotateToPrincipalAxisSystemAction.hpp	(revision d223d5e07c484fa4ef7f54e88e97fa78b5c0b504)
+++ src/Actions/MoleculeAction/RotateToPrincipalAxisSystemAction.hpp	(revision 6e5084ba8942d5b07318fb32c642c0078c9dafe5)
@@ -13,9 +13,10 @@
 
 class MoleculeListClass;
+class Vector;
 
-void MoleculeRotateToPrincipalAxisSystem();
+void MoleculeRotateToPrincipalAxisSystem(Vector &Axis);
 
 class MoleculeRotateToPrincipalAxisSystemAction : public Action {
-  friend void MoleculeRotateToPrincipalAxisSystem();
+  friend void MoleculeRotateToPrincipalAxisSystem(Vector &Axis);
 
 public:
