diff --git a/src/core/energy_inline.hpp b/src/core/energy_inline.hpp
index 66ccfb19d00..4797ee91994 100644
--- a/src/core/energy_inline.hpp
+++ b/src/core/energy_inline.hpp
@@ -275,26 +275,31 @@ calc_bonded_energy(Bonded_IA_Parameters const &iaparams, Particle const &p1,
throw BondInvalidSizeError(n_partners);
}
-/** Calculate kinetic energies for one particle.
- * @param[in] p1 particle for which to calculate energies
+/** Calculate kinetic energies from translation for one particle.
+ * @param p particle for which to calculate energies
*/
-inline double calc_kinetic_energy(Particle const &p1) {
- if (p1.p.is_virtual)
- return 0.0;
-
- /* kinetic energy */
- auto res = 0.5 * p1.p.mass * p1.m.v.norm2();
+inline double translational_kinetic_energy(Particle const &p) {
+ return p.p.is_virtual ? 0. : 0.5 * p.p.mass * p.m.v.norm2();
+}
- // Note that rotational degrees of virtual sites are integrated
- // and therefore can contribute to kinetic energy
+/** Calculate kinetic energies from rotation for one particle.
+ * @param p particle for which to calculate energies
+ */
+inline double rotational_kinetic_energy(Particle const &p) {
#ifdef ROTATION
- if (p1.p.rotation) {
- /* the rotational part is added to the total kinetic energy;
- * Here we use the rotational inertia */
- res += 0.5 * (hadamard_product(p1.m.omega, p1.m.omega) * p1.p.rinertia);
- }
+ return p.p.rotation
+ ? 0.5 * (hadamard_product(p.m.omega, p.m.omega) * p.p.rinertia)
+ : 0.0;
+#else
+ return 0.0;
#endif
- return res;
+}
+
+/** Calculate kinetic energies for one particle.
+ * @param p particle for which to calculate energies
+ */
+inline double calc_kinetic_energy(Particle const &p) {
+ return translational_kinetic_energy(p) + rotational_kinetic_energy(p);
}
#endif // ENERGY_INLINE_HPP
diff --git a/src/core/unit_tests/CMakeLists.txt b/src/core/unit_tests/CMakeLists.txt
index 754601ac1da..f42a4b64aa4 100644
--- a/src/core/unit_tests/CMakeLists.txt
+++ b/src/core/unit_tests/CMakeLists.txt
@@ -51,3 +51,4 @@ unit_test(NAME reaction_ensemble_utils_test SRC
unit_test(NAME reaction_ensemble_classes_test SRC
reaction_ensemble_classes_test.cpp DEPENDS EspressoCore Boost::mpi
MPI::MPI_CXX)
+unit_test(NAME energy_test SRC energy_test.cpp DEPENDS EspressoCore)
diff --git a/src/core/unit_tests/energy_test.cpp b/src/core/unit_tests/energy_test.cpp
new file mode 100644
index 00000000000..84311f624be
--- /dev/null
+++ b/src/core/unit_tests/energy_test.cpp
@@ -0,0 +1,78 @@
+/*
+ * Copyright (C) 2021 The ESPResSo project
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see .
+ */
+
+#define BOOST_TEST_MODULE tests
+#define BOOST_TEST_DYN_LINK
+#include
+
+#include "energy_inline.hpp"
+
+BOOST_AUTO_TEST_CASE(translational_kinetic_energy_) {
+
+ Particle p;
+ p.p.mass = 2.;
+ p.m.v = {3., 4., 5.};
+
+ auto const expected = 0.5 * p.p.mass * p.m.v.norm2();
+ BOOST_TEST(translational_kinetic_energy(p) == expected);
+
+/* virtual */
+#ifdef VIRTUAL_SITES
+ {
+ Particle p;
+ p.p.mass = 2.;
+ p.p.is_virtual = true;
+ p.m.v = {3., 4., 5.};
+
+ auto const expected = 0.;
+ BOOST_TEST(translational_kinetic_energy(p) == expected);
+ }
+#endif
+}
+
+BOOST_AUTO_TEST_CASE(rotational_kinetic_energy_) {
+ BOOST_TEST(rotational_kinetic_energy(Particle{}) == 0.);
+
+#ifdef ROTATION
+ {
+ Particle p;
+ p.m.omega = {1., 2., 3.};
+ p.p.rotation = 1;
+
+ auto const expected =
+ 0.5 * (hadamard_product(p.m.omega, p.m.omega) * p.p.rinertia);
+ BOOST_TEST(rotational_kinetic_energy(p) == expected);
+ }
+#endif
+}
+
+BOOST_AUTO_TEST_CASE(kinetic_energy_) {
+ Particle p;
+ p.p.mass = 2.;
+ p.m.v = {3., 4., 5.};
+
+#ifdef ROTATION
+ p.m.omega = {1., 2., 3.};
+ p.p.rotation = 1;
+#endif
+
+ auto const expected =
+ translational_kinetic_energy(p) + rotational_kinetic_energy(p);
+ BOOST_TEST(calc_kinetic_energy(p) == expected);
+}
\ No newline at end of file