Skip to content

Commit

Permalink
core: Consider rotational kinetic energy of virtual particles.
Browse files Browse the repository at this point in the history
  • Loading branch information
fweik committed Apr 3, 2021
1 parent 0c342e6 commit a4d5b76
Show file tree
Hide file tree
Showing 3 changed files with 100 additions and 16 deletions.
37 changes: 21 additions & 16 deletions src/core/energy_inline.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
1 change: 1 addition & 0 deletions src/core/unit_tests/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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)
78 changes: 78 additions & 0 deletions src/core/unit_tests/energy_test.cpp
Original file line number Diff line number Diff line change
@@ -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 <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;
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);
}

0 comments on commit a4d5b76

Please sign in to comment.