Rotation calculated for spherical particles (#274)

* Optimization of Quaternion for vector {0,0,0}

* Minor optimization: Avoid normalizing identity quaternions

* Profile with gprof; Reverted minor changes

* Feedback from PR

* CMake option and documentation for profiling

---------

Co-authored-by: Ruben Horn <>

cmake/modules/options.cmake

@@ -17,7 +17,14 @@ else()
     message(FATAL_ERROR "wrong precision option ")
 endif()
 
-
+# ---- PROFILING ----
+option(ENABLE_GPROF "Use the GNU profiler gprof (Only supported by GNU/Clang compiler)" OFF)
+if(ENABLE_GPROF)
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg")
+    if(NOT (CMAKE_CXX_COMPILER_ID MATCHES "Clang" OR CMAKE_CXX_COMPILER_ID MATCHES "GNU"))
+        message(WARNING "${CMAKE_CXX_COMPILER_ID} may not support the -pg option.\n(Only use GNU/Clang with the GNU profiler!)\n")
+    endif()
+endif()
 
 #set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DQUICKSCHED")
 #set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DTASKTIMINGPROFILE")

doc/profiling.dox

@@ -0,0 +1,10 @@
+/** \page Profiling
+
+In order to assess the runtime for parts of the code or even individual functions, one may use the <a href="https://ftp.gnu.org/old-gnu/Manuals/gprof-2.9.1/html_chapter/gprof_toc.html">GNU profiler</a> by following the steps below:
+
+-# Set the option `ENABLE_GPROF` to `ON` and compile ls1. (Make sure to use either the GNU or Clang compiler.)
+-# Run an experiment and make sure that the program exits normally.
+-# Rename the output file `gmon.out` generated by the profiler as it is overwritten by subsequent runs.
+-# Extract the human-readable report using `gprof build/src/MarDyn gmon.out`.
+
+*/

src/molecules/Quaternion.cpp

@@ -41,44 +41,48 @@ void Quaternion::multiply_left(const Quaternion& q) {
 }
 
 std::array<double, 3> Quaternion::rotate(const std::array<double, 3>& d) const {
-	std::array<double, 3> drot;
-	double ww = m_qw*m_qw;
-	double xx = m_qx*m_qx;
-	double yy = m_qy*m_qy;
-	double zz = m_qz*m_qz;
-	double wx = m_qw*m_qx;
-	double wy = m_qw*m_qy;
-	double wz = m_qw*m_qz;
-	double xy = m_qx*m_qy;
-	double xz = m_qx*m_qz;
-	double yz = m_qy*m_qz;
-	//          1-2*(yy+zz)
-	drot[0] = (ww+xx-yy-zz)*d[0] + 2.*(xy-wz)*d[1] + 2.*(wy+xz)*d[2];
-	//                            1-2*(xx+zz)
-	drot[1] = 2.*(wz+xy)*d[0] + (ww-xx+yy-zz)*d[1] + 2.*(yz-wx)*d[2];
-	//                                              1-2*(xx+yy)
-	drot[2] = 2.*(xz-wy)*d[0] + 2.*(wx+yz)*d[1] + (ww-xx-yy+zz)*d[2];
+	std::array<double, 3> drot{};
+	if(std::array<double, 3>{} != d) {
+		double ww = m_qw*m_qw;
+		double xx = m_qx*m_qx;
+		double yy = m_qy*m_qy;
+		double zz = m_qz*m_qz;
+		double wx = m_qw*m_qx;
+		double wy = m_qw*m_qy;
+		double wz = m_qw*m_qz;
+		double xy = m_qx*m_qy;
+		double xz = m_qx*m_qz;
+		double yz = m_qy*m_qz;
+		//          1-2*(yy+zz)
+		drot[0] = (ww+xx-yy-zz)*d[0] + 2.*(xy-wz)*d[1] + 2.*(wy+xz)*d[2];
+		//                            1-2*(xx+zz)
+		drot[1] = 2.*(wz+xy)*d[0] + (ww-xx+yy-zz)*d[1] + 2.*(yz-wx)*d[2];
+		//                                              1-2*(xx+yy)
+		drot[2] = 2.*(xz-wy)*d[0] + 2.*(wx+yz)*d[1] + (ww-xx-yy+zz)*d[2];
+	}
 	return drot;
 }
 
 std::array<double, 3> Quaternion::rotateinv(const std::array<double, 3>& d) const {
-	std::array<double, 3> drot;
-	double ww = m_qw*m_qw;
-	double xx = m_qx*m_qx;
-	double yy = m_qy*m_qy;
-	double zz = m_qz*m_qz;
-	double wx = m_qw*m_qx;
-	double wy = m_qw*m_qy;
-	double wz = m_qw*m_qz;
-	double xy = m_qx*m_qy;
-	double xz = m_qx*m_qz;
-	double yz = m_qy*m_qz;
-	//          1-2*(yy+zz)
-	drot[0] = (ww+xx-yy-zz)*d[0] + 2.*(xy+wz)*d[1] + 2.*(xz-wy)*d[2];
-	//                            1-2*(xx+zz)
-	drot[1] = 2.*(xy-wz)*d[0] + (ww-xx+yy-zz)*d[1] + 2.*(yz+wx)*d[2];
-	//                                              1-2*(xx+yy)
-	drot[2] = 2.*(xz+wy)*d[0] + 2.*(yz-wx)*d[1] + (ww-xx-yy+zz)*d[2];
+	std::array<double, 3> drot{};
+	if (std::array<double, 3>{} != d) {
+		double ww = m_qw*m_qw;
+		double xx = m_qx*m_qx;
+		double yy = m_qy*m_qy;
+		double zz = m_qz*m_qz;
+		double wx = m_qw*m_qx;
+		double wy = m_qw*m_qy;
+		double wz = m_qw*m_qz;
+		double xy = m_qx*m_qy;
+		double xz = m_qx*m_qz;
+		double yz = m_qy*m_qz;
+		//          1-2*(yy+zz)
+		drot[0] = (ww+xx-yy-zz)*d[0] + 2.*(xy+wz)*d[1] + 2.*(xz-wy)*d[2];
+		//                            1-2*(xx+zz)
+		drot[1] = 2.*(xy-wz)*d[0] + (ww-xx+yy-zz)*d[1] + 2.*(yz+wx)*d[2];
+		//                                              1-2*(xx+yy)
+		drot[2] = 2.*(xz+wy)*d[0] + 2.*(yz-wx)*d[1] + (ww-xx-yy+zz)*d[2];
+	}
 	return drot;
 }
 

src/molecules/Quaternion.h

@@ -100,7 +100,6 @@ class Quaternion {
 
 private:
 	double m_qw, m_qx, m_qy, m_qz; // components
-
 };
 
 #endif /*QUATERNION_H_*/