api/nblib/listed_forces/bondtypes.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2020,2021, 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/exception.h"
  58 #include "nblib/particletype.h"
  59 #include "nblib/util/util.hpp"
  60
  61 namespace nblib
  62 {
  63 using Name          = std::string;
  64 using ForceConstant = real;
  65 using EquilConstant = 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 equilConstant
  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, EquilConstant d) : forceConstant_(f), equilConstant_(d)
  86     {
  87     }
  88
  89     [[nodiscard]] const ForceConstant& forceConstant() const { return forceConstant_; }
  90     [[nodiscard]] const EquilConstant& equilConstant() const { return equilConstant_; }
  91
  92 private:
  93     ForceConstant forceConstant_;
  94     EquilConstant equilConstant_;
  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.equilConstant())
 101            < std::tie(b.forceConstant(), b.equilConstant());
 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.equilConstant())
 108            == std::tie(b.forceConstant(), b.equilConstant());
 109 }
 110
 111 /*! \brief harmonic bond type
 112  *
 113  *  It represents the interaction of the form
 114  *  V(r, forceConstant, equilConstant) = 0.5 * forceConstant * (r - equilConstant)^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, equilConstant) = 0.25 * forceConstant * (r^2 - equilConstant^2)^2
 123  */
 124 class G96BondType : public TwoParameterInteraction<struct G96BondTypeParameter>
 125 {
 126 public:
 127     G96BondType() = default;
 128     //! \brief Store square of equilibrium distance
 129     G96BondType(ForceConstant f, EquilConstant equilConstant) :
 130         TwoParameterInteraction<struct G96BondTypeParameter>{ f, equilConstant * equilConstant }
 131     {
 132     }
 133 };
 134
 135
 136 /*! \brief FENE bond type
 137  *
 138  * It represents the interaction of the form
 139  * V(r, forceConstant, equilConstant) = - 0.5 * forceConstant * equilConstant^2 * log( 1 - (r / equilConstant)^2)
 140  */
 141 using FENEBondType = TwoParameterInteraction<struct FENEBondTypeParameter>;
 142
 143
 144 /*! \brief Half-attractive quartic bond type
 145  *
 146  * It represents the interaction of the form
 147  * V(r, forceConstant, equilConstant) = 0.5 * forceConstant * (r - equilConstant)^4
 148  */
 149 using HalfAttractiveQuarticBondType =
 150         TwoParameterInteraction<struct HalfAttractiveQuarticBondTypeParameter>;
 151
 152
 153 /*! \brief Cubic bond type
 154  *
 155  * It represents the interaction of the form
 156  * V(r, quadraticForceConstant, cubicForceConstant, equilConstant) = quadraticForceConstant * (r -
 157  * equilConstant)^2 + quadraticForceConstant * cubicForceConstant * (r - equilConstant)
 158  */
 159 class CubicBondType
 160 {
 161 public:
 162     CubicBondType() = default;
 163     CubicBondType(ForceConstant fq, ForceConstant fc, EquilConstant d) :
 164         quadraticForceConstant_(fq), cubicForceConstant_(fc), equilDistance_(d)
 165     {
 166     }
 167
 168     [[nodiscard]] const ForceConstant& quadraticForceConstant() const
 169     {
 170         return quadraticForceConstant_;
 171     }
 172     [[nodiscard]] const ForceConstant& cubicForceConstant() const { return cubicForceConstant_; }
 173     [[nodiscard]] const EquilConstant& equilDistance() const { return equilDistance_; }
 174
 175 private:
 176     ForceConstant quadraticForceConstant_;
 177     ForceConstant cubicForceConstant_;
 178     EquilConstant equilDistance_;
 179 };
 180
 181 inline bool operator<(const CubicBondType& a, const CubicBondType& b)
 182 {
 183     return std::tie(a.quadraticForceConstant(), a.cubicForceConstant(), a.equilDistance())
 184            < std::tie(b.quadraticForceConstant(), b.cubicForceConstant(), b.equilDistance());
 185 }
 186
 187 inline bool operator==(const CubicBondType& a, const CubicBondType& b)
 188 {
 189     return std::tie(a.quadraticForceConstant(), a.cubicForceConstant(), a.equilDistance())
 190            == std::tie(b.quadraticForceConstant(), b.cubicForceConstant(), b.equilDistance());
 191 }
 192
 193 /*! \brief Morse bond type
 194  *
 195  * It represents the interaction of the form
 196  * V(r, forceConstant, exponent, equilConstant) = forceConstant * ( 1 - exp( -exponent * (r - equilConstant))
 197  */
 198 class MorseBondType
 199 {
 200 public:
 201     MorseBondType() = default;
 202     MorseBondType(ForceConstant f, Exponent e, EquilConstant d) :
 203         forceConstant_(f), exponent_(e), equilDistance_(d)
 204     {
 205     }
 206
 207     [[nodiscard]] const ForceConstant& forceConstant() const { return forceConstant_; }
 208     [[nodiscard]] const Exponent&      exponent() const { return exponent_; }
 209     [[nodiscard]] const EquilConstant& equilDistance() const { return equilDistance_; }
 210
 211 private:
 212     ForceConstant forceConstant_;
 213     Exponent      exponent_;
 214     EquilConstant 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 inline bool operator==(const MorseBondType& a, const MorseBondType& b)
 224 {
 225     return std::tie(a.forceConstant(), a.exponent(), a.equilDistance())
 226            == std::tie(b.forceConstant(), b.exponent(), b.equilDistance());
 227 }
 228
 229 /*! \brief Non-Bonded Pair Interaction Type
 230  *
 231  * It represents the interaction of the form
 232  * of LJ interactions, but occur between atoms
 233  * of the same bonded chain
 234  * V(r, c_6, c_12) = c_12*(r^-12) - c_6*(r^-6)
 235  */
 236 class PairLJType
 237 {
 238 public:
 239     PairLJType() = default;
 240     PairLJType(C6 c6, C12 c12) : c6_(c6), c12_(c12) {}
 241
 242     [[nodiscard]] const C6&  c6() const { return c6_; }
 243     [[nodiscard]] const C12& c12() const { return c12_; }
 244
 245 private:
 246     C6  c6_;
 247     C12 c12_;
 248 };
 249
 250 inline bool operator<(const PairLJType& a, const PairLJType& b)
 251 {
 252     return std::tie(a.c6(), a.c12()) < std::tie(b.c6(), b.c12());
 253 }
 254
 255 inline bool operator==(const PairLJType& a, const PairLJType& b)
 256 {
 257     return std::tie(a.c6(), a.c12()) == std::tie(b.c6(), b.c12());
 258 }
 259
 260 /*! \brief Basic template for interactions with 2 parameters named forceConstant and equilAngle
 261  *
 262  * \tparam Phantom unused template parameter for type distinction
 263  *
 264  * Distinct angle types can be generated from this template with using declarations
 265  * and declared, but undefined structs. For example:
 266  * using HarmonicAngleType = AngleInteractionType<struct HarmonicAngleParameter>;
 267  * HarmonicAngleParameter does not have to be defined.
 268  *
 269  * Note: the angle is always stored as radians internally
 270  */
 271 template<class Phantom>
 272 class AngleInteractionType : public TwoParameterInteraction<Phantom>
 273 {
 274 public:
 275     AngleInteractionType() = default;
 276     //! \brief construct from angle given in radians
 277     AngleInteractionType(ForceConstant f, Radians angle) :
 278         TwoParameterInteraction<Phantom>{ f, angle }
 279     {
 280     }
 281
 282     //! \brief construct from angle given in degrees
 283     AngleInteractionType(ForceConstant f, Degrees angle) :
 284         TwoParameterInteraction<Phantom>{ f, angle * DEG2RAD }
 285     {
 286     }
 287 };
 288
 289 /*! \brief Harmonic angle type
 290  *
 291  * It represents the interaction of the form
 292  * V(theta, forceConstant, equilConstant) = 0.5 * forceConstant * (theta - equilConstant)^2
 293  */
 294 using HarmonicAngle = AngleInteractionType<struct HarmonicAngleParameter>;
 295
 296 /*! \brief linear angle type
 297  *
 298  * It represents the interaction of the form
 299  * V(theta, forceConstant, a) = 0.5 * forceConstant * (dr)^2
 300  * where dr = - a * r_ij - (1 - a) * r_kj
 301  */
 302 using LinearAngle = TwoParameterInteraction<struct LinearAngleParameter>;
 303
 304 /*! \brief Basic template for angle types that use the cosines of their equilibrium angles
 305  *         in their potential expression
 306  */
 307 template<class Phantom>
 308 class CosineParamAngle : public TwoParameterInteraction<Phantom>
 309 {
 310 public:
 311     CosineParamAngle() = default;
 312     //! \brief construct from angle given in radians
 313     CosineParamAngle(ForceConstant f, Radians angle) :
 314         TwoParameterInteraction<Phantom>{ f, std::cos(angle) }
 315     {
 316     }
 317
 318     //! \brief construct from angle given in degrees
 319     CosineParamAngle(ForceConstant f, Degrees angle) :
 320         TwoParameterInteraction<Phantom>{ f, std::cos(angle * DEG2RAD) }
 321     {
 322     }
 323 };
 324
 325 /*! \brief G96 or Cosine-based angle type
 326  *
 327  * This represents the interaction of the form
 328  * V(cos(theta), forceConstant, cos(equilConstant)) = 0.5 * forceConstant * (cos(theta) - cos(equilConstant))^2
 329  */
 330 using G96Angle = CosineParamAngle<struct G96AngleParameter>;
 331
 332 /*! \brief Restricted angle type
 333  *
 334  * This represents the interaction of the form
 335  * V(cos(theta), forceConstant, cos(equilConstant)) =
 336  *     0.5 * forceConstant * (cos(theta) - cos(equilConstant))^2 / (sin(theta))^2
 337  */
 338 using RestrictedAngle = CosineParamAngle<struct RestrictedAngleParameter>;
 339
 340 /*! \brief Quartic angle type
 341  *
 342  * It represents the interaction of the form of a fourth order polynomial
 343  * V(theta, forceConstant, equilConstant) = sum[i = 0 -> 4](forceConstant_i * (theta - equilConstant)^i
 344  */
 345 class QuarticAngle
 346 {
 347 public:
 348     QuarticAngle() = default;
 349     //! \brief construct from given angle in radians
 350     QuarticAngle(ForceConstant f0, ForceConstant f1, ForceConstant f2, ForceConstant f3, ForceConstant f4, Radians angle) :
 351         forceConstants_{ f0, f1, f2, f3, f4 }, equilConstant_(angle)
 352     {
 353     }
 354
 355     //! \brief construct from given angle in degrees
 356     QuarticAngle(ForceConstant f0, ForceConstant f1, ForceConstant f2, ForceConstant f3, ForceConstant f4, Degrees angle) :
 357         forceConstants_{ f0, f1, f2, f3, f4 }, equilConstant_(angle * DEG2RAD)
 358     {
 359     }
 360
 361     [[nodiscard]] const std::array<ForceConstant, 5>& forceConstants() const
 362     {
 363         return forceConstants_;
 364     }
 365
 366     ForceConstant forceConstant(int order) const
 367     {
 368         switch (order)
 369         {
 370             case 0: return forceConstants_[0];
 371             case 1: return forceConstants_[1];
 372             case 2: return forceConstants_[2];
 373             case 3: return forceConstants_[3];
 374             case 4: return forceConstants_[4];
 375             default:
 376                 throw InputException(
 377                         "Please enter a value between 0-4 for the Quartic Angle force constants");
 378         }
 379     }
 380
 381     Radians equilConstant() const { return equilConstant_; }
 382
 383 private:
 384     std::array<ForceConstant, 5> forceConstants_;
 385     Radians                      equilConstant_;
 386 };
 387
 388 inline bool operator<(const QuarticAngle& a, const QuarticAngle& b)
 389 {
 390     return (a.forceConstants() < b.forceConstants()) && (a.equilConstant() < b.equilConstant());
 391 }
 392
 393 inline bool operator==(const QuarticAngle& a, const QuarticAngle& b)
 394 {
 395     return (a.forceConstants() == b.forceConstants()) && (a.equilConstant() == b.equilConstant());
 396 }
 397
 398
 399 /*! \brief Cross bond-bond interaction type
 400  */
 401 class CrossBondBond
 402 {
 403 public:
 404     CrossBondBond() = default;
 405     CrossBondBond(ForceConstant f, EquilConstant r0ij, EquilConstant r0kj) :
 406         forceConstant_(f), r0ij_(r0ij), r0kj_(r0kj)
 407     {
 408     }
 409
 410     [[nodiscard]] const ForceConstant& forceConstant() const { return forceConstant_; }
 411     [[nodiscard]] const EquilConstant& equilDistanceIJ() const { return r0ij_; }
 412     [[nodiscard]] const EquilConstant& equilDistanceKJ() const { return r0kj_; }
 413
 414 private:
 415     ForceConstant forceConstant_;
 416     EquilConstant r0ij_;
 417     EquilConstant r0kj_;
 418 };
 419
 420 inline bool operator<(const CrossBondBond& a, const CrossBondBond& b)
 421 {
 422     return std::tie(a.forceConstant(), a.equilDistanceIJ(), a.equilDistanceKJ())
 423            < std::tie(b.forceConstant(), b.equilDistanceIJ(), b.equilDistanceKJ());
 424 }
 425
 426 inline bool operator==(const CrossBondBond& a, const CrossBondBond& b)
 427 {
 428     return std::tie(a.forceConstant(), a.equilDistanceIJ(), a.equilDistanceKJ())
 429            == std::tie(b.forceConstant(), b.equilDistanceIJ(), b.equilDistanceKJ());
 430 }
 431
 432 /*! \brief Cross bond-angle interaction type
 433  */
 434 class CrossBondAngle
 435 {
 436 public:
 437     CrossBondAngle() = default;
 438     CrossBondAngle(ForceConstant f, EquilConstant r0ij, EquilConstant r0kj, EquilConstant r0ik) :
 439         forceConstant_(f), r0ij_(r0ij), r0kj_(r0kj), r0ik_(r0ik)
 440     {
 441     }
 442
 443     [[nodiscard]] const ForceConstant& forceConstant() const { return forceConstant_; }
 444     [[nodiscard]] const EquilConstant& equilDistanceIJ() const { return r0ij_; }
 445     [[nodiscard]] const EquilConstant& equilDistanceKJ() const { return r0kj_; }
 446     [[nodiscard]] const EquilConstant& equilDistanceIK() const { return r0ik_; }
 447
 448 private:
 449     ForceConstant forceConstant_;
 450     EquilConstant r0ij_;
 451     EquilConstant r0kj_;
 452     EquilConstant r0ik_;
 453 };
 454
 455 inline bool operator<(const CrossBondAngle& a, const CrossBondAngle& b)
 456 {
 457     return std::tie(a.forceConstant(), a.equilDistanceIJ(), a.equilDistanceKJ(), a.equilDistanceIK())
 458            < std::tie(b.forceConstant(), b.equilDistanceIJ(), b.equilDistanceKJ(), b.equilDistanceIK());
 459 }
 460
 461 inline bool operator==(const CrossBondAngle& a, const CrossBondAngle& b)
 462 {
 463     return std::tie(a.forceConstant(), a.equilDistanceIJ(), a.equilDistanceKJ(), a.equilDistanceIK())
 464            == std::tie(b.forceConstant(), b.equilDistanceIJ(), b.equilDistanceKJ(), b.equilDistanceIK());
 465 }
 466
 467 /*! \brief Proper Dihedral Implementation
 468  */
 469 class ProperDihedral
 470 {
 471 public:
 472     using Multiplicity = int;
 473
 474     ProperDihedral() = default;
 475     ProperDihedral(Radians phi, ForceConstant f, Multiplicity m) :
 476         phi_(phi), forceConstant_(f), multiplicity_(m)
 477     {
 478     }
 479     ProperDihedral(Degrees phi, ForceConstant f, Multiplicity m) :
 480         phi_(phi * DEG2RAD), forceConstant_(f), multiplicity_(m)
 481     {
 482     }
 483
 484     [[nodiscard]] const EquilConstant& equilDistance() const { return phi_; }
 485     [[nodiscard]] const ForceConstant& forceConstant() const { return forceConstant_; }
 486     [[nodiscard]] const Multiplicity&  multiplicity() const { return multiplicity_; }
 487
 488 private:
 489     EquilConstant phi_;
 490     ForceConstant forceConstant_;
 491     Multiplicity  multiplicity_;
 492 };
 493
 494 inline bool operator<(const ProperDihedral& a, const ProperDihedral& b)
 495 {
 496     return std::tie(a.equilDistance(), a.forceConstant(), a.multiplicity())
 497            < std::tie(b.equilDistance(), b.forceConstant(), b.multiplicity());
 498 }
 499
 500 inline bool operator==(const ProperDihedral& a, const ProperDihedral& b)
 501 {
 502     return std::tie(a.equilDistance(), a.forceConstant(), a.multiplicity())
 503            == std::tie(b.equilDistance(), b.forceConstant(), b.multiplicity());
 504 }
 505
 506
 507 /*! \brief Improper Dihedral Implementation
 508  */
 509 class ImproperDihedral : public TwoParameterInteraction<struct ImproperDihdedralParameter>
 510 {
 511 public:
 512     ImproperDihedral() = default;
 513     ImproperDihedral(Radians phi, ForceConstant f) :
 514         TwoParameterInteraction<struct ImproperDihdedralParameter>{ f, phi }
 515     {
 516     }
 517     ImproperDihedral(Degrees phi, ForceConstant f) :
 518         TwoParameterInteraction<struct ImproperDihdedralParameter>{ f, phi * DEG2RAD }
 519     {
 520     }
 521 };
 522
 523 /*! \brief Ryckaert-Belleman Dihedral Implementation
 524  */
 525 class RyckaertBellemanDihedral
 526 {
 527 public:
 528     RyckaertBellemanDihedral() = default;
 529     RyckaertBellemanDihedral(real p1, real p2, real p3, real p4, real p5, real p6) :
 530         parameters_{ p1, p2, p3, p4, p5, p6 }
 531     {
 532     }
 533
 534     const real& operator[](std::size_t i) const { return parameters_[i]; }
 535
 536     [[nodiscard]] const std::array<real, 6>& parameters() const { return parameters_; }
 537
 538     [[nodiscard]] std::size_t size() const { return parameters_.size(); }
 539
 540 private:
 541     std::array<real, 6> parameters_;
 542 };
 543
 544 inline bool operator<(const RyckaertBellemanDihedral& a, const RyckaertBellemanDihedral& b)
 545 {
 546     return a.parameters() < b.parameters();
 547 }
 548
 549 inline bool operator==(const RyckaertBellemanDihedral& a, const RyckaertBellemanDihedral& b)
 550 {
 551     return a.parameters() == b.parameters();
 552 }
 553
 554
 555 /*! \brief Type for 5-center interaction (C-MAP)
 556  *
 557  *  Note: no kernels currently implemented
 558  */
 559 class Default5Center
 560 {
 561 public:
 562     Default5Center() = default;
 563     Default5Center(Radians phi, Radians psi, ForceConstant fphi, ForceConstant fpsi) :
 564         phi_(phi), psi_(psi), fphi_(fphi), fpsi_(fpsi)
 565     {
 566     }
 567
 568     [[nodiscard]] const Radians&       phi() const { return phi_; }
 569     [[nodiscard]] const Radians&       psi() const { return psi_; }
 570     [[nodiscard]] const ForceConstant& fphi() const { return fphi_; }
 571     [[nodiscard]] const ForceConstant& fpsi() const { return fpsi_; }
 572
 573 private:
 574     Radians       phi_, psi_;
 575     ForceConstant fphi_, fpsi_;
 576 };
 577
 578 inline bool operator<(const Default5Center& a, const Default5Center& b)
 579 {
 580     return std::tie(a.phi(), a.psi(), a.fphi(), a.fpsi())
 581            < std::tie(b.phi(), b.psi(), b.fphi(), b.fpsi());
 582 }
 583
 584 inline bool operator==(const Default5Center& a, const Default5Center& b)
 585 {
 586     return std::tie(a.phi(), a.psi(), a.fphi(), a.fpsi())
 587            == std::tie(b.phi(), b.psi(), b.fphi(), b.fpsi());
 588 }
 589
 590
 591 } // namespace nblib
 592 #endif // NBLIB_LISTEDFORCES_BONDTYPES_H