From 0bfef142d52169fd8a1d53a317be0a81dfbd007e Mon Sep 17 00:00:00 2001 From: Patrick Kreissl Date: Thu, 30 Mar 2023 12:29:32 +0200 Subject: [PATCH 1/2] Add bonded interaction TorsionBond. --- src/core/bonded_interactions/CMakeLists.txt | 1 + .../bonded_interaction_data.hpp | 3 +- .../bonded_interactions.dox | 2 +- src/core/bonded_interactions/torsion_bond.cpp | 28 ++++++ src/core/bonded_interactions/torsion_bond.hpp | 86 +++++++++++++++++++ src/core/energy_inline.hpp | 5 ++ src/core/forces_inline.hpp | 10 +++ src/python/espressomd/interactions.py | 24 ++++++ .../interactions/BondedInteraction.hpp | 15 ++++ .../interactions/initialize.cpp | 1 + 10 files changed, 173 insertions(+), 2 deletions(-) create mode 100644 src/core/bonded_interactions/torsion_bond.cpp create mode 100644 src/core/bonded_interactions/torsion_bond.hpp diff --git a/src/core/bonded_interactions/CMakeLists.txt b/src/core/bonded_interactions/CMakeLists.txt index 3c4ad3afa5d..6749a14617c 100644 --- a/src/core/bonded_interactions/CMakeLists.txt +++ b/src/core/bonded_interactions/CMakeLists.txt @@ -26,4 +26,5 @@ target_sources( ${CMAKE_CURRENT_SOURCE_DIR}/fene.cpp ${CMAKE_CURRENT_SOURCE_DIR}/rigid_bond.cpp ${CMAKE_CURRENT_SOURCE_DIR}/thermalized_bond.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/torsion_bond.cpp ${CMAKE_CURRENT_SOURCE_DIR}/thermalized_bond_utils.cpp) diff --git a/src/core/bonded_interactions/bonded_interaction_data.hpp b/src/core/bonded_interactions/bonded_interaction_data.hpp index 2f006186cec..53d787f04f1 100644 --- a/src/core/bonded_interactions/bonded_interaction_data.hpp +++ b/src/core/bonded_interactions/bonded_interaction_data.hpp @@ -41,6 +41,7 @@ #include "quartic.hpp" #include "rigid_bond.hpp" #include "thermalized_bond.hpp" +#include "torsion_bond.hpp" #include "TabulatedPotential.hpp" @@ -96,7 +97,7 @@ using Bonded_IA_Parameters = BondedCoulombSR, AngleHarmonicBond, AngleCosineBond, AngleCossquareBond, DihedralBond, TabulatedDistanceBond, TabulatedAngleBond, TabulatedDihedralBond, ThermalizedBond, - RigidBond, IBMTriel, IBMVolCons, IBMTribend, + TorsionBond, RigidBond, IBMTriel, IBMVolCons, IBMTribend, OifGlobalForcesBond, OifLocalForcesBond, VirtualBond>; class BondedInteractionsMap { diff --git a/src/core/bonded_interactions/bonded_interactions.dox b/src/core/bonded_interactions/bonded_interactions.dox index ab7eb29a285..c6a66855f8a 100644 --- a/src/core/bonded_interactions/bonded_interactions.dox +++ b/src/core/bonded_interactions/bonded_interactions.dox @@ -61,7 +61,7 @@ * @endcode * The return value of \c cutoff() should be as large as the interaction * range of the new interaction. This is only relevant to pairwise bonds - * and dihedral bonds. In all other cases, the return value should be -1, + * and dihedral bonds. In all other cases, the return value should be 0, * namely angle bonds as well as other bonds that don't have an interaction * range. * * @code{.cpp} diff --git a/src/core/bonded_interactions/torsion_bond.cpp b/src/core/bonded_interactions/torsion_bond.cpp new file mode 100644 index 00000000000..3cb3b98cd04 --- /dev/null +++ b/src/core/bonded_interactions/torsion_bond.cpp @@ -0,0 +1,28 @@ +/* + * Copyright (C) 2010-2022 The ESPResSo project + * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 + * Max-Planck-Institute for Polymer Research, Theory Group + * + * 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 . + */ +/** \file + * + * Implementation of \ref torsion_bond.hpp + */ + +#include "torsion_bond.hpp" + +TorsionBond::TorsionBond(double k) { this->k = k; } diff --git a/src/core/bonded_interactions/torsion_bond.hpp b/src/core/bonded_interactions/torsion_bond.hpp new file mode 100644 index 00000000000..1e296c1c7d3 --- /dev/null +++ b/src/core/bonded_interactions/torsion_bond.hpp @@ -0,0 +1,86 @@ +/* + * Copyright (C) 2010-2022 The ESPResSo project + * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 + * Max-Planck-Institute for Polymer Research, Theory Group + * + * 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 . + */ +#ifndef TORSION_BOND_HPP +#define TORSION_BOND_HPP +/** \file + * Routines to calculate the torsion_bond potential between particle pairs. + * + * Implementation in \ref torsion_bond.cpp. + */ + +#include "config/config.hpp" + +#include "Particle.hpp" +#include + +#include + +/** Parameters for torsion_bond Potential. */ +struct TorsionBond { + /** spring constant */ + double k; + + static constexpr int num = 1; + + double cutoff() const { return 0; } + + TorsionBond(double k); + + boost::optional torque(Particle const &p1, + Particle const &p2) const; + boost::optional energy(Particle const &p1, Particle const &p2) const; + +private: + friend boost::serialization::access; + template + void serialize(Archive &ar, long int /* version */) { + ar &k; + } +}; + +/** Compute the torque resulting from the torsion_bond between p1 and p2. + * @param[in] p1 %First particle. + * @param[in] p2 %Second particle. + */ +inline boost::optional +TorsionBond::torque(Particle const &p1, Particle const &p2) const { +#ifdef ROTATION + return -k * vector_product(p1.calc_director(), p2.calc_director()); +#else + return {}; +#endif +} + +/** Compute the torsion_bond energy. + * @param[in] p1 %First particle. + * @param[in] p2 %Second particle. + */ +inline boost::optional TorsionBond::energy(Particle const &p1, + Particle const &p2) const { +#ifdef ROTATION + double const angle = std::acos(p1.calc_director() * p2.calc_director()); + return 0.5 * k * angle * angle; +#else + return {}; +#endif +} + +#endif // TORSION_BOND_HPP diff --git a/src/core/energy_inline.hpp b/src/core/energy_inline.hpp index eaa24054fd3..08faae0c9e0 100644 --- a/src/core/energy_inline.hpp +++ b/src/core/energy_inline.hpp @@ -238,6 +238,11 @@ calc_bonded_energy(Bonded_IA_Parameters const &iaparams, Particle const &p1, if (auto const *iap = boost::get(&iaparams)) { return iap->energy(dx); } +#endif +#ifdef ROTATION + if (auto const *iap = boost::get(&iaparams)) { + return iap->energy(p1, *p2); + } #endif if (boost::get(&iaparams)) { return {0.}; diff --git a/src/core/forces_inline.hpp b/src/core/forces_inline.hpp index 935cf1632c4..adaaaac8750 100644 --- a/src/core/forces_inline.hpp +++ b/src/core/forces_inline.hpp @@ -309,6 +309,16 @@ inline bool add_bonded_two_body_force( p1.force() += std::get<0>(forces); p2.force() += std::get<1>(forces); + return false; + } + } else if (auto const *iap = boost::get(&iaparams)) { + auto result = iap->torque(p1, p2); + if (result) { + auto const &torque = result.get(); + + p1.torque() -= torque; + p2.torque() += torque; + return false; } } else { diff --git a/src/python/espressomd/interactions.py b/src/python/espressomd/interactions.py index f6e9578cf80..f960f2a3de3 100644 --- a/src/python/espressomd/interactions.py +++ b/src/python/espressomd/interactions.py @@ -783,6 +783,7 @@ class BONDED_IA(enum.IntEnum): TABULATED_ANGLE = enum.auto() TABULATED_DIHEDRAL = enum.auto() THERMALIZED_DIST = enum.auto() + TORSION_BOND = enum.auto() RIGID_BOND = enum.auto() IBM_TRIEL = enum.auto() IBM_VOLUME_CONSERVATION = enum.auto() @@ -944,6 +945,29 @@ def get_default_params(self): return {"r_cut": 0.} +@script_interface_register +class TorsionBond(BondedInteraction): + + """ + Torsion bond. + + Parameters + ---------- + k : :obj:`float` + Spring constant for the bond interaction. + + """ + + _so_name = "Interactions::TorsionBond" + _type_number = BONDED_IA.TORSION_BOND + + def get_default_params(self): + """Gets default values of optional parameters. + + """ + return {} + + @script_interface_register class BondedCoulomb(BondedInteraction): diff --git a/src/script_interface/interactions/BondedInteraction.hpp b/src/script_interface/interactions/BondedInteraction.hpp index fa22352491e..6e737a1b133 100644 --- a/src/script_interface/interactions/BondedInteraction.hpp +++ b/src/script_interface/interactions/BondedInteraction.hpp @@ -298,6 +298,21 @@ class AngleCossquareBond : public BondedInteractionImpl<::AngleCossquareBond> { } }; +class TorsionBond : public BondedInteractionImpl<::TorsionBond> { +public: + TorsionBond() { + add_parameters({ + {"k", AutoParameter::read_only, [this]() { return get_struct().k; }}, + }); + } + +private: + void construct_bond(VariantMap const ¶ms) override { + m_bonded_ia = std::make_shared<::Bonded_IA_Parameters>( + CoreBondedInteraction(get_value(params, "k"))); + } +}; + class DihedralBond : public BondedInteractionImpl<::DihedralBond> { public: DihedralBond() { diff --git a/src/script_interface/interactions/initialize.cpp b/src/script_interface/interactions/initialize.cpp index 57129779939..fa2f4a726e0 100644 --- a/src/script_interface/interactions/initialize.cpp +++ b/src/script_interface/interactions/initialize.cpp @@ -44,6 +44,7 @@ void initialize(Utils::Factory *om) { om->register_new( "Interactions::TabulatedDihedralBond"); om->register_new("Interactions::ThermalizedBond"); + om->register_new("Interactions::TorsionBond"); om->register_new("Interactions::RigidBond"); om->register_new("Interactions::IBMTriel"); om->register_new("Interactions::IBMVolCons"); From 35e09cbcadf024f94e0248796394976a25b4cb2f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 30 Mar 2023 19:02:37 +0200 Subject: [PATCH 2/2] core: Implement TorsionBond::energy --- src/core/bonded_interactions/torsion_bond.hpp | 4 ++-- src/core/forces_inline.hpp | 8 ++++++-- src/python/espressomd/interactions.py | 1 + 3 files changed, 9 insertions(+), 4 deletions(-) diff --git a/src/core/bonded_interactions/torsion_bond.hpp b/src/core/bonded_interactions/torsion_bond.hpp index 1e296c1c7d3..1e68b975fcf 100644 --- a/src/core/bonded_interactions/torsion_bond.hpp +++ b/src/core/bonded_interactions/torsion_bond.hpp @@ -76,8 +76,8 @@ TorsionBond::torque(Particle const &p1, Particle const &p2) const { inline boost::optional TorsionBond::energy(Particle const &p1, Particle const &p2) const { #ifdef ROTATION - double const angle = std::acos(p1.calc_director() * p2.calc_director()); - return 0.5 * k * angle * angle; + auto const cos_alpha = p1.calc_director() * p2.calc_director(); + return k * (1. - cos_alpha); #else return {}; #endif diff --git a/src/core/forces_inline.hpp b/src/core/forces_inline.hpp index adaaaac8750..771ea7e1a9e 100644 --- a/src/core/forces_inline.hpp +++ b/src/core/forces_inline.hpp @@ -311,7 +311,9 @@ inline bool add_bonded_two_body_force( return false; } - } else if (auto const *iap = boost::get(&iaparams)) { + } +#ifdef ROTATION + else if (auto const *iap = boost::get(&iaparams)) { auto result = iap->torque(p1, p2); if (result) { auto const &torque = result.get(); @@ -321,7 +323,9 @@ inline bool add_bonded_two_body_force( return false; } - } else { + } +#endif + else { auto result = calc_bond_pair_force(p1, p2, iaparams, dx, kernel); if (result) { p1.force() += result.get(); diff --git a/src/python/espressomd/interactions.py b/src/python/espressomd/interactions.py index f960f2a3de3..680d36eed6b 100644 --- a/src/python/espressomd/interactions.py +++ b/src/python/espressomd/interactions.py @@ -959,6 +959,7 @@ class TorsionBond(BondedInteraction): """ _so_name = "Interactions::TorsionBond" + _so_feature = "ROTATION" _type_number = BONDED_IA.TORSION_BOND def get_default_params(self):