Add bonded interaction TorsionBond.

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)

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 {

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}

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 <http://www.gnu.org/licenses/>.
+ */
+/** \file
+ *
+ *  Implementation of \ref torsion_bond.hpp
+ */
+
+#include "torsion_bond.hpp"
+
+TorsionBond::TorsionBond(double k) { this->k = k; }

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 <http://www.gnu.org/licenses/>.
+ */
+#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 <utils/Vector.hpp>
+
+#include <boost/optional.hpp>
+
+/** 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<Utils::Vector3d> torque(Particle const &p1,
+                                          Particle const &p2) const;
+  boost::optional<double> energy(Particle const &p1, Particle const &p2) const;
+
+private:
+  friend boost::serialization::access;
+  template <typename Archive>
+  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<Utils::Vector3d>
+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<double> TorsionBond::energy(Particle const &p1,
+                                                   Particle const &p2) const {
+#ifdef ROTATION
+  auto const cos_alpha = p1.calc_director() * p2.calc_director();
+  return k * (1. - cos_alpha);
+#else
+  return {};
+#endif
+}
+
+#endif // TORSION_BOND_HPP

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<TabulatedDistanceBond>(&iaparams)) {
       return iap->energy(dx);
     }
+#endif
+#ifdef ROTATION
+    if (auto const *iap = boost::get<TorsionBond>(&iaparams)) {
+      return iap->energy(p1, *p2);
+    }
 #endif
     if (boost::get<VirtualBond>(&iaparams)) {
       return {0.};

src/core/forces_inline.hpp

@@ -311,7 +311,21 @@ inline bool add_bonded_two_body_force(
 
       return false;
     }
-  } else {
+  }
+#ifdef ROTATION
+  else if (auto const *iap = boost::get<TorsionBond>(&iaparams)) {
+    auto result = iap->torque(p1, p2);
+    if (result) {
+      auto const &torque = result.get();
+
+      p1.torque() -= torque;
+      p2.torque() += torque;
+
+      return false;
+    }
+  }
+#endif
+  else {
     auto result = calc_bond_pair_force(p1, p2, iaparams, dx, kernel);
     if (result) {
       p1.force() += result.get();

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,30 @@ 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"
+    _so_feature = "ROTATION"
+    _type_number = BONDED_IA.TORSION_BOND
+
+    def get_default_params(self):
+        """Gets default values of optional parameters.
+
+        """
+        return {}
+
+
 @script_interface_register
 class BondedCoulomb(BondedInteraction):
 

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 &params) override {
+    m_bonded_ia = std::make_shared<::Bonded_IA_Parameters>(
+        CoreBondedInteraction(get_value<double>(params, "k")));
+  }
+};
+
 class DihedralBond : public BondedInteractionImpl<::DihedralBond> {
 public:
   DihedralBond() {

src/script_interface/interactions/initialize.cpp

@@ -44,6 +44,7 @@ void initialize(Utils::Factory<ObjectHandle> *om) {
   om->register_new<TabulatedDihedralBond>(
       "Interactions::TabulatedDihedralBond");
   om->register_new<ThermalizedBond>("Interactions::ThermalizedBond");
+  om->register_new<TorsionBond>("Interactions::TorsionBond");
   om->register_new<RigidBond>("Interactions::RigidBond");
   om->register_new<IBMTriel>("Interactions::IBMTriel");
   om->register_new<IBMVolCons>("Interactions::IBMVolCons");