api/nblib/interactions.cpp

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