api/nblib/interactions.cpp

   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 /*! \internal \file
  36  * \brief
  37  * Implements nblib particle-types interactions
  38  *
  39  * \author Victor Holanda <victor.holanda@cscs.ch>
  40  * \author Joe Jordan <ejjordan@kth.se>
  41  * \author Prashanth Kanduri <kanduri@cscs.ch>
  42  * \author Sebastian Keller <keller@cscs.ch>
  43  * \author Artem Zhmurov <zhmurov@gmail.com>
  44  */
  45 #include <cmath>
  46 #include <set>
  47
  48 #include "nblib/exception.h"
  49 #include "nblib/interactions.h"
  50
  51 namespace nblib
  52 {
  53
  54 C6 NonBondedInteractionMap::getC6(const ParticleTypeName& first, const ParticleTypeName& second) const
  55 {
  56     return std::get<0>(interactionMap_.at(std::make_tuple(first, second)));
  57 }
  58
  59 C12 NonBondedInteractionMap::getC12(const ParticleTypeName& first, const ParticleTypeName& second) const
  60 {
  61     return std::get<1>(interactionMap_.at(std::make_tuple(first, second)));
  62 }
  63
  64 void NonBondedInteractionMap::setInteractions(const ParticleTypeName& first,
  65                                               const ParticleTypeName& second,
  66                                               const C6                c6_combo,
  67                                               const C12               c12_combo)
  68 {
  69     auto interactionKey             = std::make_tuple(first, second);
  70     interactionMap_[interactionKey] = std::make_tuple(c6_combo, c12_combo);
  71 }
  72
  73 size_t NonBondedInteractionMap::count(const NonBondedInteractionMap::NamePairTuple& namePairTuple)
  74 {
  75     return interactionMap_.count(namePairTuple);
  76 }
  77
  78 namespace
  79 {
  80
  81 //! Combines the non-bonded parameters from two particles for pairwise interactions
  82 real combineNonbondedParameters(real v, real w, CombinationRule combinationRule)
  83 {
  84     if (combinationRule == CombinationRule::Geometric)
  85     {
  86         return std::sqrt(v * w);
  87     }
  88     else
  89     {
  90         throw InputException("unknown LJ Combination rule specified\n");
  91     }
  92 }
  93
  94 } // namespace
  95
  96 ParticleTypesInteractions::ParticleTypesInteractions(CombinationRule cr) : combinationRule_(cr) {}
  97
  98 ParticleTypesInteractions& ParticleTypesInteractions::add(const ParticleTypeName& particleTypeName,
  99                                                           C6                      c6,
 100                                                           C12                     c12)
 101 {
 102     auto insertLocation = singleParticleInteractionsMap_.insert(
 103             std::make_pair(particleTypeName, std::make_tuple(c6, c12)));
 104
 105     if (!insertLocation.second) // if particleTypeName already existed
 106     {
 107         if (std::get<0>(insertLocation.first->second) != c6
 108             || std::get<1>(insertLocation.first->second) != c12)
 109         {
 110             std::string message = formatString(
 111                     "Attempting to add nonbonded interaction parameters for particle "
 112                     "type {} twice",
 113                     particleTypeName.value());
 114             throw InputException(message);
 115         }
 116     }
 117     return *this;
 118 }
 119
 120 ParticleTypesInteractions& ParticleTypesInteractions::add(const ParticleTypeName& particleTypeName1,
 121                                                           const ParticleTypeName& particleTypeName2,
 122                                                           C6                      c6,
 123                                                           C12                     c12)
 124 {
 125     auto interactionKey         = std::make_tuple(particleTypeName1, particleTypeName2);
 126     auto possibleInteractionKey = std::make_tuple(particleTypeName2, particleTypeName1);
 127
 128     auto insertLocation = twoParticlesInteractionsMap_.insert(
 129             std::make_pair(interactionKey, std::make_tuple(c6, c12)));
 130     twoParticlesInteractionsMap_.insert(std::make_pair(possibleInteractionKey, std::make_tuple(c6, c12)));
 131
 132     if (!insertLocation.second) // if particleTypeName already existed
 133     {
 134         if (std::get<0>(insertLocation.first->second) != c6
 135             || std::get<1>(insertLocation.first->second) != c12)
 136         {
 137             std::string message = formatString(
 138                     "Attempting to add nonbonded interaction parameters between the particle types "
 139                     "{} {} twice",
 140                     particleTypeName1.value(),
 141                     particleTypeName2.value());
 142             throw InputException(message);
 143         }
 144     }
 145     return *this;
 146 }
 147
 148 NonBondedInteractionMap ParticleTypesInteractions::generateTable() const
 149 {
 150     NonBondedInteractionMap nonbondedParameters_;
 151
 152     // creating the combination rule based interaction matrix
 153     for (const auto& particleType1 : singleParticleInteractionsMap_)
 154     {
 155         C6  c6_1  = std::get<0>(particleType1.second);
 156         C12 c12_1 = std::get<1>(particleType1.second);
 157
 158         for (const auto& particleType2 : singleParticleInteractionsMap_)
 159         {
 160             C6  c6_2  = std::get<0>(particleType2.second);
 161             C12 c12_2 = std::get<1>(particleType2.second);
 162
 163             C6  c6_combo{ combineNonbondedParameters(c6_1, c6_2, combinationRule_) };
 164             C12 c12_combo{ combineNonbondedParameters(c12_1, c12_2, combinationRule_) };
 165
 166             nonbondedParameters_.setInteractions(
 167                     particleType1.first, particleType2.first, c6_combo, c12_combo);
 168         }
 169     }
 170
 171     // updating the interaction matrix based on the user fine tuned parameters
 172     for (const auto& particleTypeTuple : twoParticlesInteractionsMap_)
 173     {
 174         const auto& first  = std::get<0>(particleTypeTuple.first);
 175         const auto& second = std::get<1>(particleTypeTuple.first);
 176
 177         C6  c6_combo  = std::get<0>(particleTypeTuple.second);
 178         C12 c12_combo = std::get<1>(particleTypeTuple.second);
 179
 180         nonbondedParameters_.setInteractions(first, second, c6_combo, c12_combo);
 181     }
 182
 183     std::set<ParticleTypeName> particleTypes;
 184     for (auto const& typeKey : nonbondedParameters_)
 185     { // we don't need to get<1> because the list is guaranteed to be symmetric
 186         particleTypes.insert(std::get<0>(typeKey.first));
 187     }
 188
 189     // check whether there is any missing interaction
 190     for (const ParticleTypeName& particleTypeName1 : particleTypes)
 191     {
 192         for (const ParticleTypeName& particleTypeName2 : particleTypes)
 193         {
 194             auto interactionKey = std::make_tuple(particleTypeName1, particleTypeName2);
 195             if (nonbondedParameters_.count(interactionKey) == 0)
 196             {
 197                 std::string message = formatString("Missing interaction between {} {}",
 198                                                    particleTypeName1.value(),
 199                                                    particleTypeName2.value());
 200                 throw InputException(message);
 201             }
 202         }
 203     }
 204     return nonbondedParameters_;
 205 }
 206
 207 CombinationRule ParticleTypesInteractions::getCombinationRule() const
 208 {
 209     return combinationRule_;
 210 }
 211
 212 void ParticleTypesInteractions::merge(const ParticleTypesInteractions& other)
 213 {
 214     for (const auto& keyval : other.singleParticleInteractionsMap_)
 215     {
 216         add(keyval.first, std::get<0>(keyval.second), std::get<1>(keyval.second));
 217     }
 218
 219     for (const auto& keyval : other.twoParticlesInteractionsMap_)
 220     {
 221         add(std::get<0>(keyval.first),
 222             std::get<1>(keyval.first),
 223             std::get<0>(keyval.second),
 224             std::get<1>(keyval.second));
 225     }
 226 }
 227
 228 } // namespace nblib