Merge branch release-2021
[alexxy/gromacs.git] / api / nblib / molecules.cpp
index 679d888cd160ac2a7026297bddada1db3522b07c..9b243bf6cdae5430f6a29b5b01a85e19a50ffbe5 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2020, by the GROMACS development team, led by
+ * Copyright (c) 2020,2021, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
 #include "nblib/exception.h"
 #include "nblib/molecules.h"
 #include "nblib/particletype.h"
-#include "nblib/util/internal.h"
 
 namespace nblib
 {
 
-
 Molecule::Molecule(MoleculeName moleculeName) : name_(std::move(moleculeName)) {}
 
 MoleculeName Molecule::name() const
@@ -115,139 +113,135 @@ Molecule& Molecule::addParticle(const ParticleName& particleName, const Particle
     return *this;
 }
 
-//! Two-particle interactions such as bonds and LJ1-4
-template<class Interaction>
-void Molecule::addInteraction(const ParticleName& particleNameI,
-                              const ResidueName&  residueNameI,
-                              const ParticleName& particleNameJ,
-                              const ResidueName&  residueNameJ,
-                              const Interaction&  interaction)
+ResidueName Molecule::residueName(const ParticleIdentifier& particleIdentifier)
+{
+    return (particleIdentifier.residueName() == ResidueName{}) ? ResidueName(name_)
+                                                               : particleIdentifier.residueName();
+}
+
+template<class ListedVariant, class... ParticleIdentifiers>
+void Molecule::addInteractionImpl(const ListedVariant& interaction, const ParticleIdentifiers&... particles)
+{
+    auto storeInteraction = [&](const auto& interaction_) {
+        using Interaction = std::decay_t<decltype(interaction_)>;
+
+        auto& interactionContainer = pickType<Interaction>(interactionData_);
+        interactionContainer.interactions_.emplace_back(particles...);
+        interactionContainer.interactionTypes_.push_back(interaction_);
+    };
+
+    // add the interaction to the correct location in interactionData_
+    std::visit(storeInteraction, interaction);
+}
+
+void Molecule::addInteraction(const ParticleIdentifier&   particleI,
+                              const ParticleIdentifier&   particleJ,
+                              const TwoCenterInteraction& interaction)
 {
-    if (particleNameI == particleNameJ and residueNameI == residueNameJ)
+    if (particleI == particleJ)
     {
         throw InputException(std::string("Cannot add interaction of particle ")
-                             + particleNameI.value() + " with itself in molecule " + name_.value());
+                             + particleI.particleName().value() + " with itself in molecule "
+                             + name_.value());
     }
 
-    auto& interactionContainer = pickType<Interaction>(interactionData_);
-    interactionContainer.interactions_.emplace_back(particleNameI, residueNameI, particleNameJ, residueNameJ);
-    interactionContainer.interactionTypes_.push_back(interaction);
+    addInteractionImpl(interaction,
+                       particleI.particleName(),
+                       residueName(particleI),
+                       particleJ.particleName(),
+                       residueName(particleJ));
 }
 
-//! \cond DO_NOT_DOCUMENT
-#define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x)                               \
-    template void Molecule::addInteraction(const ParticleName& particleNameI, \
-                                           const ResidueName&  residueNameI,  \
-                                           const ParticleName& particleNameJ, \
-                                           const ResidueName&  residueNameJ,  \
-                                           const x&            interaction);
-MAP(ADD_INTERACTION_INSTANTIATE_TEMPLATE, SUPPORTED_TWO_CENTER_TYPES)
-#undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
-//! \endcond
-
-// add interactions with default residue name
-template<class Interaction>
-void Molecule::addInteraction(const ParticleName& particleNameI,
-                              const ParticleName& particleNameJ,
-                              const Interaction&  interaction)
+void Molecule::addInteraction(const ParticleIdentifier&     particleI,
+                              const ParticleIdentifier&     particleJ,
+                              const ParticleIdentifier&     particleK,
+                              const ThreeCenterInteraction& interaction)
 {
-    addInteraction(particleNameI, ResidueName(name_), particleNameJ, ResidueName(name_), interaction);
+    if (particleI == particleJ and particleJ == particleK)
+    {
+        throw InputException(std::string("Cannot add interaction of particle ")
+                             + particleI.particleName().value() + " with itself in molecule "
+                             + name_.value());
+    }
+
+    addInteractionImpl(interaction,
+                       particleI.particleName(),
+                       residueName(particleI),
+                       particleJ.particleName(),
+                       residueName(particleJ),
+                       particleK.particleName(),
+                       residueName(particleK));
 }
 
-//! \cond DO_NOT_DOCUMENT
-#define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x) \
-    template void Molecule::addInteraction(     \
-            const ParticleName& particleNameI, const ParticleName& particleNameJ, const x& interaction);
-MAP(ADD_INTERACTION_INSTANTIATE_TEMPLATE, SUPPORTED_TWO_CENTER_TYPES)
-#undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
-
-//! 3-particle interactions such as angles
-template<class Interaction>
-void Molecule::addInteraction(const ParticleName& particleNameI,
-                              const ResidueName&  residueNameI,
-                              const ParticleName& particleNameJ,
-                              const ResidueName&  residueNameJ,
-                              const ParticleName& particleNameK,
-                              const ResidueName&  residueNameK,
-                              const Interaction&  interaction)
+void Molecule::addInteraction(const ParticleIdentifier&    particleI,
+                              const ParticleIdentifier&    particleJ,
+                              const ParticleIdentifier&    particleK,
+                              const ParticleIdentifier&    particleL,
+                              const FourCenterInteraction& interaction)
 {
-    if (particleNameI == particleNameJ and particleNameJ == particleNameK)
+    if (particleI == particleJ and particleJ == particleK and particleK == particleL)
     {
         throw InputException(std::string("Cannot add interaction of particle ")
-                             + particleNameI.value() + " with itself in molecule " + name_.value());
+                             + particleI.particleName().value() + " with itself in molecule "
+                             + name_.value());
     }
 
-    auto& interactionContainer = pickType<Interaction>(interactionData_);
-    interactionContainer.interactions_.emplace_back(
-            particleNameI, residueNameI, particleNameJ, residueNameJ, particleNameK, residueNameK);
-    interactionContainer.interactionTypes_.push_back(interaction);
+    addInteractionImpl(interaction,
+                       particleI.particleName(),
+                       residueName(particleI),
+                       particleJ.particleName(),
+                       residueName(particleJ),
+                       particleK.particleName(),
+                       residueName(particleK),
+                       particleL.particleName(),
+                       residueName(particleL));
 }
 
-#define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x)                               \
-    template void Molecule::addInteraction(const ParticleName& particleNameI, \
-                                           const ResidueName&  residueNameI,  \
-                                           const ParticleName& particleNameJ, \
-                                           const ResidueName&  residueNameJ,  \
-                                           const ParticleName& particleNameK, \
-                                           const ResidueName&  residueNameK,  \
-                                           const x&            interaction);
-MAP(ADD_INTERACTION_INSTANTIATE_TEMPLATE, SUPPORTED_THREE_CENTER_TYPES)
-#undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
-
-template<class Interaction>
-void Molecule::addInteraction(const ParticleName& particleNameI,
-                              const ParticleName& particleNameJ,
-                              const ParticleName& particleNameK,
-                              const Interaction&  interaction)
+void Molecule::addInteraction(const ParticleIdentifier&    particleI,
+                              const ParticleIdentifier&    particleJ,
+                              const ParticleIdentifier&    particleK,
+                              const ParticleIdentifier&    particleL,
+                              const ParticleIdentifier&    particleM,
+                              const FiveCenterInteraction& interaction)
 {
-    addInteraction(particleNameI,
-                   ResidueName(name_),
-                   particleNameJ,
-                   ResidueName(name_),
-                   particleNameK,
-                   ResidueName(name_),
-                   interaction);
+    if (particleI == particleJ and particleJ == particleK and particleK == particleL and particleL == particleM)
+    {
+        throw InputException(std::string("Cannot add interaction of particle ")
+                             + particleI.particleName().value() + " with itself in molecule "
+                             + name_.value());
+    }
+
+    addInteractionImpl(interaction,
+                       particleI.particleName(),
+                       residueName(particleI),
+                       particleJ.particleName(),
+                       residueName(particleJ),
+                       particleK.particleName(),
+                       residueName(particleK),
+                       particleL.particleName(),
+                       residueName(particleL),
+                       particleM.particleName(),
+                       residueName(particleM));
 }
 
-#define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x)                               \
-    template void Molecule::addInteraction(const ParticleName& particleNameI, \
-                                           const ParticleName& particleNameJ, \
-                                           const ParticleName& particleNameK, \
-                                           const x&            interaction);
-MAP(ADD_INTERACTION_INSTANTIATE_TEMPLATE, SUPPORTED_THREE_CENTER_TYPES)
-#undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
-//! \endcond
 
 int Molecule::numParticlesInMolecule() const
 {
     return particles_.size();
 }
 
-void Molecule::addExclusion(const int particleIndex, const int particleIndexToExclude)
+void Molecule::addExclusion(const ParticleIdentifier& particle, const ParticleIdentifier& particleToExclude)
 {
-    // We do not need to add exclusion in case the particle indexes are the same
-    // because self exclusion are added by addParticle
-    if (particleIndex != particleIndexToExclude)
+    if (particle == particleToExclude)
     {
-        exclusions_.emplace_back(particleIndex, particleIndexToExclude);
-        exclusions_.emplace_back(particleIndexToExclude, particleIndex);
+        return;
     }
-}
 
-void Molecule::addExclusion(std::tuple<ParticleName, ResidueName> particle,
-                            std::tuple<ParticleName, ResidueName> particleToExclude)
-{
     // duplication for the swapped pair happens in getExclusions()
-    exclusionsByName_.emplace_back(std::make_tuple(std::get<0>(particle),
-                                                   std::get<1>(particle),
-                                                   std::get<0>(particleToExclude),
-                                                   std::get<1>(particleToExclude)));
-}
-
-void Molecule::addExclusion(const ParticleName& particleName, const ParticleName& particleNameToExclude)
-{
-    addExclusion(std::make_tuple(particleName, ResidueName(name_)),
-                 std::make_tuple(particleNameToExclude, ResidueName(name_)));
+    exclusionsByName_.emplace_back(std::make_tuple(particle.particleName(),
+                                                   residueName(particle),
+                                                   particleToExclude.particleName(),
+                                                   residueName(particleToExclude)));
 }
 
 const Molecule::InteractionTuple& Molecule::interactionData() const