core: Consider rotational kinetic energy of virtual particles.

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

src/core/unit_tests/CMakeLists.txt

@@ -46,3 +46,4 @@ unit_test(NAME LocalBox_test SRC LocalBox_test.cpp DEPENDS EspressoCore)
 unit_test(NAME thermostats_test SRC thermostats_test.cpp DEPENDS EspressoCore)
 unit_test(NAME random_test SRC random_test.cpp DEPENDS EspressoUtils Random123)
 unit_test(NAME BondList_test SRC BondList_test.cpp DEPENDS EspressoCore)
+unit_test(NAME energy_test SRC energy_test.cpp DEPENDS EspressoCore)

src/core/unit_tests/energy_test.cpp

@@ -0,0 +1,84 @@
+/*
+ * 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 <http://www.gnu.org/licenses/>.
+ */
+
+#define BOOST_TEST_MODULE tests
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+
+#include "energy_inline.hpp"
+
+BOOST_AUTO_TEST_CASE(translational_kinetic_energy_) {
+
+  Particle p;
+#ifdef MASS
+  p.p.mass = 2.;
+#endif
+  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;
+#ifdef MASS
+    p.p.mass = 2.;
+#endif
+    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;
+#ifdef MASS
+  p.p.mass = 2.;
+#endif
+  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);
+}