api/nblib/listed_forces/bondtypes.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2020, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \inpublicapi \file
  36  * \brief
  37  * Implements nblib supported bondtypes
  38  *
  39  * We choose to forward comparison operations to the
  40  * corresponding std::tuple comparison operations.
  41  * In order to do that without temporary copies,
  42  * we employ std::tie, which requires lvalues as input.
  43  * For this reason, bond type parameter getters are implemented
  44  * with a const lvalue reference return.
  45  *
  46  * \author Victor Holanda <victor.holanda@cscs.ch>
  47  * \author Joe Jordan <ejjordan@kth.se>
  48  * \author Prashanth Kanduri <kanduri@cscs.ch>
  49  * \author Sebastian Keller <keller@cscs.ch>
  50  * \author Artem Zhmurov <zhmurov@gmail.com>
  51  */
  52 #ifndef NBLIB_LISTEDFORCES_BONDTYPES_H
  53 #define NBLIB_LISTEDFORCES_BONDTYPES_H
  54
  55 #include <array>
  56
  57 #include "nblib/particletype.h"
  58 #include "nblib/ppmap.h"
  59 #include "nblib/util/user.h"
  60
  61 namespace nblib
  62 {
  63 using Name          = std::string;
  64 using ForceConstant = real;
  65 using EquilDistance = real;
  66 using Exponent      = real;
  67
  68 using Degrees = StrongType<real, struct DegreeParameter>;
  69 using Radians = StrongType<real, struct RadianParameter>;
  70
  71 /*! \brief Basic template for interactions with 2 parameters named forceConstant and equilDistance
  72  *
  73  * \tparam Phantom unused template parameter for type distinction
  74  *
  75  * Distinct bond types can be generated from this template with using declarations
  76  * and declared, but undefined structs. For example:
  77  * using HarmonicBondType = TwoParameterInteraction<struct HarmonicBondTypeParameter>;
  78  * Note that HarmonicBondTypeParameter does not have to be defined.
  79  */
  80 template<class Phantom>
  81 class TwoParameterInteraction
  82 {
  83 public:
  84     TwoParameterInteraction() = default;
  85     TwoParameterInteraction(ForceConstant f, EquilDistance d) : forceConstant_(f), equilDistance_(d)
  86     {
  87     }
  88
  89     [[nodiscard]] const ForceConstant& forceConstant() const { return forceConstant_; }
  90     [[nodiscard]] const EquilDistance& equilDistance() const { return equilDistance_; }
  91
  92 private:
  93     ForceConstant forceConstant_;
  94     EquilDistance equilDistance_;
  95 };
  96
  97 template<class Phantom>
  98 inline bool operator<(const TwoParameterInteraction<Phantom>& a, const TwoParameterInteraction<Phantom>& b)
  99 {
 100     return std::tie(a.forceConstant(), a.equilDistance())
 101            < std::tie(b.forceConstant(), b.equilDistance());
 102 }
 103
 104 template<class Phantom>
 105 inline bool operator==(const TwoParameterInteraction<Phantom>& a, const TwoParameterInteraction<Phantom>& b)
 106 {
 107     return std::tie(a.forceConstant(), a.equilDistance())
 108            == std::tie(b.forceConstant(), b.equilDistance());
 109 }
 110
 111 /*! \brief harmonic bond type
 112  *
 113  *  It represents the interaction of the form
 114  *  V(r; forceConstant, equilDistance) = 0.5 * forceConstant * (r - equilDistance)^2
 115  */
 116 using HarmonicBondType = TwoParameterInteraction<struct HarmonicBondTypeParameter>;
 117
 118
 119 /*! \brief GROMOS bond type
 120  *
 121  * It represents the interaction of the form
 122  * V(r; forceConstant, equilDistance) = 0.25 * forceConstant * (r^2 - equilDistance^2)^2
 123  */
 124 using G96BondType = TwoParameterInteraction<struct G96BondTypeParameter>;
 125
 126
 127 /*! \brief FENE bond type
 128  *
 129  * It represents the interaction of the form
 130  * V(r; forceConstant, equilDistance) = - 0.5 * forceConstant * equilDistance^2 * log( 1 - (r / equilDistance)^2)
 131  */
 132 using FENEBondType = TwoParameterInteraction<struct FENEBondTypeParameter>;
 133
 134
 135 /*! \brief Half-attractive quartic bond type
 136  *
 137  * It represents the interaction of the form
 138  * V(r; forceConstant, equilDistance) = 0.5 * forceConstant * (r - equilDistance)^4
 139  */
 140 using HalfAttractiveQuarticBondType =
 141         TwoParameterInteraction<struct HalfAttractiveQuarticBondTypeParameter>;
 142
 143
 144 /*! \brief Cubic bond type
 145  *
 146  * It represents the interaction of the form
 147  * V(r; quadraticForceConstant, cubicForceConstant, equilDistance) = quadraticForceConstant * (r -
 148  * equilDistance)^2 + quadraticForceConstant * cubicForceConstant * (r - equilDistance)
 149  */
 150 struct CubicBondType
 151 {
 152     CubicBondType() = default;
 153     CubicBondType(ForceConstant fq, ForceConstant fc, EquilDistance d) :
 154         quadraticForceConstant_(fq),
 155         cubicForceConstant_(fc),
 156         equilDistance_(d)
 157     {
 158     }
 159
 160     [[nodiscard]] const ForceConstant& quadraticForceConstant() const
 161     {
 162         return quadraticForceConstant_;
 163     }
 164     [[nodiscard]] const ForceConstant& cubicForceConstant() const { return cubicForceConstant_; }
 165     [[nodiscard]] const EquilDistance& equilDistance() const { return equilDistance_; }
 166
 167 private:
 168     ForceConstant quadraticForceConstant_;
 169     ForceConstant cubicForceConstant_;
 170     EquilDistance equilDistance_;
 171 };
 172
 173 inline bool operator<(const CubicBondType& a, const CubicBondType& b)
 174 {
 175     return std::tie(a.quadraticForceConstant(), a.cubicForceConstant(), a.equilDistance())
 176            < std::tie(b.quadraticForceConstant(), b.cubicForceConstant(), b.equilDistance());
 177 }
 178
 179 inline bool operator==(const CubicBondType& a, const CubicBondType& b)
 180 {
 181     return std::tie(a.quadraticForceConstant(), a.cubicForceConstant(), a.equilDistance())
 182            == std::tie(b.quadraticForceConstant(), b.cubicForceConstant(), b.equilDistance());
 183 }
 184
 185 /*! \brief Morse bond type
 186  *
 187  * It represents the interaction of the form
 188  * V(r; forceConstant, exponent, equilDistance) = forceConstant * ( 1 - exp( -exponent * (r - equilDistance))
 189  */
 190 class MorseBondType
 191 {
 192 public:
 193     MorseBondType() = default;
 194     MorseBondType(ForceConstant f, Exponent e, EquilDistance d) :
 195         forceConstant_(f),
 196         exponent_(e),
 197         equilDistance_(d)
 198     {
 199     }
 200
 201     [[nodiscard]] const ForceConstant& forceConstant() const { return forceConstant_; }
 202     [[nodiscard]] const Exponent&      exponent() const { return exponent_; }
 203     [[nodiscard]] const EquilDistance& equilDistance() const { return equilDistance_; }
 204
 205 private:
 206     ForceConstant forceConstant_;
 207     Exponent      exponent_;
 208     EquilDistance equilDistance_;
 209 };
 210
 211 inline bool operator<(const MorseBondType& a, const MorseBondType& b)
 212 {
 213     return std::tie(a.forceConstant(), a.exponent(), a.equilDistance())
 214            < std::tie(b.forceConstant(), b.exponent(), b.equilDistance());
 215 }
 216
 217 inline bool operator==(const MorseBondType& a, const MorseBondType& b)
 218 {
 219     return std::tie(a.forceConstant(), a.exponent(), a.equilDistance())
 220            == std::tie(b.forceConstant(), b.exponent(), b.equilDistance());
 221 }
 222
 223
 224 /*! \brief default angle type
 225  *
 226  * Note: the angle is always stored as radians internally
 227  */
 228 struct DefaultAngle : public TwoParameterInteraction<struct DefaultAngleParameter>
 229 {
 230     DefaultAngle() = default;
 231     //! \brief construct from angle given in radians
 232     DefaultAngle(Radians angle, ForceConstant f) :
 233         TwoParameterInteraction<struct DefaultAngleParameter>{ f, angle }
 234     {
 235     }
 236
 237     //! \brief construct from angle given in degrees
 238     DefaultAngle(Degrees angle, ForceConstant f) :
 239         TwoParameterInteraction<struct DefaultAngleParameter>{ f, angle * DEG2RAD }
 240     {
 241     }
 242 };
 243
 244 /*! \brief Proper Dihedral Implementation
 245  */
 246 class ProperDihedral
 247 {
 248 public:
 249     using Multiplicity = int;
 250
 251     ProperDihedral() = default;
 252     ProperDihedral(Radians phi, ForceConstant f, Multiplicity m) :
 253         phi_(phi),
 254         forceConstant_(f),
 255         multiplicity_(m)
 256     {
 257     }
 258     ProperDihedral(Degrees phi, ForceConstant f, Multiplicity m) :
 259         phi_(phi * DEG2RAD),
 260         forceConstant_(f),
 261         multiplicity_(m)
 262     {
 263     }
 264
 265     [[nodiscard]] const EquilDistance& equilDistance() const { return phi_; }
 266     [[nodiscard]] const ForceConstant& forceConstant() const { return forceConstant_; }
 267     [[nodiscard]] const Multiplicity&  multiplicity() const { return multiplicity_; }
 268
 269 private:
 270     EquilDistance phi_;
 271     ForceConstant forceConstant_;
 272     Multiplicity  multiplicity_;
 273 };
 274
 275 inline bool operator<(const ProperDihedral& a, const ProperDihedral& b)
 276 {
 277     return std::tie(a.equilDistance(), a.forceConstant(), a.multiplicity())
 278            < std::tie(b.equilDistance(), b.forceConstant(), b.multiplicity());
 279 }
 280
 281 inline bool operator==(const ProperDihedral& a, const ProperDihedral& b)
 282 {
 283     return std::tie(a.equilDistance(), a.forceConstant(), a.multiplicity())
 284            == std::tie(b.equilDistance(), b.forceConstant(), b.multiplicity());
 285 }
 286
 287
 288 /*! \brief Improper Dihedral Implementation
 289  */
 290 struct ImproperDihedral : public TwoParameterInteraction<struct ImproperDihdedralParameter>
 291 {
 292     ImproperDihedral() = default;
 293     ImproperDihedral(Radians phi, ForceConstant f) :
 294         TwoParameterInteraction<struct ImproperDihdedralParameter>{ f, phi }
 295     {
 296     }
 297     ImproperDihedral(Degrees phi, ForceConstant f) :
 298         TwoParameterInteraction<struct ImproperDihdedralParameter>{ f, phi * DEG2RAD }
 299     {
 300     }
 301 };
 302
 303 /*! \brief Ryckaert-Belleman Dihedral Implementation
 304  */
 305 class RyckaertBellemanDihedral
 306 {
 307 public:
 308     RyckaertBellemanDihedral() = default;
 309     RyckaertBellemanDihedral(real p1, real p2, real p3, real p4, real p5, real p6) :
 310         parameters_{ p1, p2, p3, p4, p5, p6 }
 311     {
 312     }
 313
 314     const real& operator[](std::size_t i) const { return parameters_[i]; }
 315
 316     [[nodiscard]] const std::array<real, 6>& parameters() const { return parameters_; }
 317
 318     [[nodiscard]] std::size_t size() const { return parameters_.size(); }
 319
 320 private:
 321     std::array<real, 6> parameters_;
 322 };
 323
 324 inline bool operator<(const RyckaertBellemanDihedral& a, const RyckaertBellemanDihedral& b)
 325 {
 326     return a.parameters() < b.parameters();
 327 }
 328
 329 inline bool operator==(const RyckaertBellemanDihedral& a, const RyckaertBellemanDihedral& b)
 330 {
 331     return a.parameters() == b.parameters();
 332 }
 333
 334
 335 /*! \brief Type for 5-center interaction (C-MAP)
 336  *
 337  *  Note: no kernels currently implemented
 338  */
 339 class Default5Center
 340 {
 341 public:
 342     Default5Center() = default;
 343     Default5Center(Radians phi, Radians psi, ForceConstant fphi, ForceConstant fpsi) :
 344         phi_(phi),
 345         psi_(psi),
 346         fphi_(fphi),
 347         fpsi_(fpsi)
 348     {
 349     }
 350
 351     [[nodiscard]] const Radians&       phi() const { return phi_; }
 352     [[nodiscard]] const Radians&       psi() const { return psi_; }
 353     [[nodiscard]] const ForceConstant& fphi() const { return fphi_; }
 354     [[nodiscard]] const ForceConstant& fpsi() const { return fpsi_; }
 355
 356 private:
 357     Radians       phi_, psi_;
 358     ForceConstant fphi_, fpsi_;
 359 };
 360
 361 inline bool operator<(const Default5Center& a, const Default5Center& b)
 362 {
 363     return std::tie(a.phi(), a.psi(), a.fphi(), a.fpsi())
 364            < std::tie(b.phi(), b.psi(), b.fphi(), b.fpsi());
 365 }
 366
 367 inline bool operator==(const Default5Center& a, const Default5Center& b)
 368 {
 369     return std::tie(a.phi(), a.psi(), a.fphi(), a.fpsi())
 370            == std::tie(b.phi(), b.psi(), b.fphi(), b.fpsi());
 371 }
 372
 373
 374 } // namespace nblib
 375 #endif // NBLIB_LISTEDFORCES_BONDTYPES_H