src/gromacs/gmxpreprocess/gpp_atomtype.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2011,2014,2015,2017,2018 by the GROMACS development team.
   7  * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "gpp_atomtype.h"
  41
  42 #include <climits>
  43 #include <cmath>
  44 #include <cstring>
  45
  46 #include <algorithm>
  47 #include <optional>
  48
  49 #include "gromacs/gmxpreprocess/grompp_impl.h"
  50 #include "gromacs/gmxpreprocess/notset.h"
  51 #include "gromacs/gmxpreprocess/topdirs.h"
  52 #include "gromacs/gmxpreprocess/toputil.h"
  53 #include "gromacs/math/vecdump.h"
  54 #include "gromacs/topology/ifunc.h"
  55 #include "gromacs/topology/symtab.h"
  56 #include "gromacs/topology/topology.h"
  57 #include "gromacs/utility/fatalerror.h"
  58 #include "gromacs/utility/smalloc.h"
  59
  60 struct AtomTypeData
  61 {
  62     //! Explicit constructor.
  63     AtomTypeData(const t_atom& a, char** name, const InteractionOfType& nb, const int bondAtomType, const int atomNumber) :
  64         atom_(a), name_(name), nb_(nb), bondAtomType_(bondAtomType), atomNumber_(atomNumber)
  65     {
  66     }
  67     //! Actual atom data.
  68     t_atom atom_;
  69     //! Atom name.
  70     char** name_;
  71     //! Nonbonded data.
  72     InteractionOfType nb_;
  73     //! Bonded atomtype for the type.
  74     int bondAtomType_;
  75     //! Atom number for the atom type.
  76     int atomNumber_;
  77 };
  78
  79 class PreprocessingAtomTypes::Impl
  80 {
  81 public:
  82     //! The number for currently loaded entries.
  83     size_t size() const { return types.size(); }
  84     //! The actual atom type data.
  85     std::vector<AtomTypeData> types;
  86     //! Map from \c types[i].name to \c i for quick look-up in \ref atomTypeFromName. Ref #3974.
  87     std::unordered_map<std::string, int> nameToAtomType;
  88 };
  89
  90 bool PreprocessingAtomTypes::isSet(int nt) const
  91 {
  92     return ((nt >= 0) && (nt < gmx::ssize(*this)));
  93 }
  94
  95 std::optional<int> PreprocessingAtomTypes::atomTypeFromName(const std::string& str) const
  96 {
  97     /* Atom types are always case sensitive */
  98     const auto found = impl_->nameToAtomType.find(str);
  99     if (found == impl_->nameToAtomType.end())
 100     {
 101         return std::nullopt;
 102     }
 103     else
 104     {
 105         GMX_ASSERT(str == *impl_->types[found->second].name_,
 106                    "Invalid data in atomTypeFromName lookup table");
 107         return std::make_optional(found->second);
 108     }
 109 }
 110
 111 size_t PreprocessingAtomTypes::size() const
 112 {
 113     return impl_->size();
 114 }
 115
 116 std::optional<const char*> PreprocessingAtomTypes::atomNameFromAtomType(int nt) const
 117 {
 118     return isSet(nt) ? std::make_optional(*(impl_->types[nt].name_)) : std::nullopt;
 119 }
 120
 121 std::optional<real> PreprocessingAtomTypes::atomMassFromAtomType(int nt) const
 122 {
 123     return isSet(nt) ? std::make_optional(impl_->types[nt].atom_.m) : std::nullopt;
 124 }
 125
 126 std::optional<real> PreprocessingAtomTypes::atomChargeFromAtomType(int nt) const
 127 {
 128     return isSet(nt) ? std::make_optional(impl_->types[nt].atom_.q) : std::nullopt;
 129 }
 130
 131 std::optional<ParticleType> PreprocessingAtomTypes::atomParticleTypeFromAtomType(int nt) const
 132 {
 133     return isSet(nt) ? std::make_optional(impl_->types[nt].atom_.ptype) : std::nullopt;
 134 }
 135
 136 std::optional<int> PreprocessingAtomTypes::bondAtomTypeFromAtomType(int nt) const
 137 {
 138     return isSet(nt) ? std::make_optional(impl_->types[nt].bondAtomType_) : std::nullopt;
 139 }
 140
 141 std::optional<int> PreprocessingAtomTypes::atomNumberFromAtomType(int nt) const
 142 {
 143     return isSet(nt) ? std::make_optional(impl_->types[nt].atomNumber_) : std::nullopt;
 144 }
 145
 146 std::optional<real> PreprocessingAtomTypes::atomNonBondedParamFromAtomType(int nt, int param) const
 147 {
 148     if (!isSet(nt))
 149     {
 150         return std::nullopt;
 151     }
 152     gmx::ArrayRef<const real> forceParam = impl_->types[nt].nb_.forceParam();
 153     if ((param < 0) || (param >= MAXFORCEPARAM))
 154     {
 155         return std::nullopt;
 156     }
 157     return std::make_optional(forceParam[param]);
 158 }
 159
 160 PreprocessingAtomTypes::PreprocessingAtomTypes() : impl_(new Impl) {}
 161
 162 PreprocessingAtomTypes::PreprocessingAtomTypes(PreprocessingAtomTypes&& old) noexcept :
 163     impl_(std::move(old.impl_))
 164 {
 165 }
 166
 167 PreprocessingAtomTypes& PreprocessingAtomTypes::operator=(PreprocessingAtomTypes&& old) noexcept
 168 {
 169     impl_ = std::move(old.impl_);
 170     return *this;
 171 }
 172
 173 PreprocessingAtomTypes::~PreprocessingAtomTypes() {}
 174
 175 int PreprocessingAtomTypes::addType(t_symtab*                tab,
 176                                     const t_atom&            a,
 177                                     const std::string&       name,
 178                                     const InteractionOfType& nb,
 179                                     int                      bondAtomType,
 180                                     int                      atomNumber)
 181 {
 182     auto position = atomTypeFromName(name);
 183     if (!position.has_value())
 184     {
 185         impl_->types.emplace_back(a, put_symtab(tab, name.c_str()), nb, bondAtomType, atomNumber);
 186         const int newType           = impl_->types.size() - 1;
 187         impl_->nameToAtomType[name] = newType;
 188         return newType;
 189     }
 190     else
 191     {
 192         return *position;
 193     }
 194 }
 195
 196 std::optional<int> PreprocessingAtomTypes::setType(int                      nt,
 197                                                    t_symtab*                tab,
 198                                                    const t_atom&            a,
 199                                                    const std::string&       name,
 200                                                    const InteractionOfType& nb,
 201                                                    int                      bondAtomType,
 202                                                    int                      atomNumber)
 203 {
 204     if (!isSet(nt))
 205     {
 206         return std::nullopt;
 207     }
 208
 209     impl_->types[nt].atom_         = a;
 210     impl_->types[nt].name_         = put_symtab(tab, name.c_str());
 211     impl_->types[nt].nb_           = nb;
 212     impl_->types[nt].bondAtomType_ = bondAtomType;
 213     impl_->types[nt].atomNumber_   = atomNumber;
 214
 215     return std::make_optional(nt);
 216 }
 217
 218 void PreprocessingAtomTypes::printTypes(FILE* out)
 219 {
 220     fprintf(out, "[ %s ]\n", enumValueToString(Directive::d_atomtypes));
 221     fprintf(out,
 222             "; %6s  %8s  %8s  %8s  %12s  %12s\n",
 223             "type",
 224             "mass",
 225             "charge",
 226             "particle",
 227             "c6",
 228             "c12");
 229     for (auto& entry : impl_->types)
 230     {
 231         fprintf(out,
 232                 "%8s  %8.3f  %8.3f  %8s  %12e  %12e\n",
 233                 *(entry.name_),
 234                 entry.atom_.m,
 235                 entry.atom_.q,
 236                 "A",
 237                 entry.nb_.c0(),
 238                 entry.nb_.c1());
 239     }
 240
 241     fprintf(out, "\n");
 242 }
 243
 244 static int search_atomtypes(const PreprocessingAtomTypes*          ga,
 245                             int*                                   n,
 246                             gmx::ArrayRef<int>                     typelist,
 247                             int                                    thistype,
 248                             gmx::ArrayRef<const InteractionOfType> interactionTypes,
 249                             int                                    ftype)
 250 {
 251     int nn    = *n;
 252     int nrfp  = NRFP(ftype);
 253     int ntype = ga->size();
 254
 255     int i;
 256     for (i = 0; (i < nn); i++)
 257     {
 258         if (typelist[i] == thistype)
 259         {
 260             /* This type number has already been added */
 261             break;
 262         }
 263
 264         /* Otherwise, check if the parameters are identical to any previously added type */
 265
 266         bool bFound = true;
 267         for (int j = 0; j < ntype && bFound; j++)
 268         {
 269             /* Check nonbonded parameters */
 270             gmx::ArrayRef<const real> forceParam1 =
 271                     interactionTypes[ntype * typelist[i] + j].forceParam();
 272             gmx::ArrayRef<const real> forceParam2 = interactionTypes[ntype * thistype + j].forceParam();
 273             for (int k = 0; (k < nrfp) && bFound; k++)
 274             {
 275                 bFound = forceParam1[k] == forceParam2[k];
 276             }
 277
 278             /* Check atomnumber */
 279             int tli = typelist[i];
 280             bFound  = bFound
 281                      && (*ga->atomNumberFromAtomType(tli) == *ga->atomNumberFromAtomType(thistype));
 282         }
 283         if (bFound)
 284         {
 285             break;
 286         }
 287     }
 288
 289     if (i == nn)
 290     {
 291         if (nn == ntype)
 292         {
 293             gmx_fatal(FARGS, "Atomtype horror n = %d, %s, %d", nn, __FILE__, __LINE__);
 294         }
 295         typelist[nn] = thistype;
 296         nn++;
 297     }
 298     *n = nn;
 299
 300     return i;
 301 }
 302
 303 void PreprocessingAtomTypes::renumberTypes(gmx::ArrayRef<InteractionsOfType> plist,
 304                                            gmx_mtop_t*                       mtop,
 305                                            int*                              wall_atomtype,
 306                                            bool                              bVerbose)
 307 {
 308     int nat, ftype, ntype;
 309
 310     ntype = size();
 311     std::vector<int> typelist(ntype);
 312
 313     if (bVerbose)
 314     {
 315         fprintf(stderr, "renumbering atomtypes...\n");
 316     }
 317
 318     /* Since the bonded interactions have been assigned now,
 319      * we want to reduce the number of atom types by merging
 320      * ones with identical nonbonded interactions, in addition
 321      * to removing unused ones.
 322      *
 323      * With QM/MM we also check that the atom numbers match
 324      */
 325
 326     /* Get nonbonded interaction type */
 327     if (plist[F_LJ].size() > 0)
 328     {
 329         ftype = F_LJ;
 330     }
 331     else
 332     {
 333         ftype = F_BHAM;
 334     }
 335
 336     /* Renumber atomtypes by first making a list of which ones are actually used.
 337      * We provide the list of nonbonded parameters so search_atomtypes
 338      * can determine if two types should be merged.
 339      */
 340     nat = 0;
 341     for (const gmx_moltype_t& moltype : mtop->moltype)
 342     {
 343         const t_atoms* atoms = &moltype.atoms;
 344         for (int i = 0; (i < atoms->nr); i++)
 345         {
 346             atoms->atom[i].type = search_atomtypes(
 347                     this, &nat, typelist, atoms->atom[i].type, plist[ftype].interactionTypes, ftype);
 348             atoms->atom[i].typeB = search_atomtypes(
 349                     this, &nat, typelist, atoms->atom[i].typeB, plist[ftype].interactionTypes, ftype);
 350         }
 351     }
 352
 353     for (int i = 0; i < 2; i++)
 354     {
 355         if (wall_atomtype[i] >= 0)
 356         {
 357             wall_atomtype[i] = search_atomtypes(
 358                     this, &nat, typelist, wall_atomtype[i], plist[ftype].interactionTypes, ftype);
 359         }
 360     }
 361
 362     std::vector<AtomTypeData> new_types;
 363     /* We now have a list of unique atomtypes in typelist */
 364
 365     /* Renumber nlist */
 366     std::vector<InteractionOfType> nbsnew;
 367
 368     // Reset the map used for fast lookups, and refill it below
 369     impl_->nameToAtomType.clear();
 370
 371     for (int i = 0; (i < nat); i++)
 372     {
 373         int mi = typelist[i];
 374         for (int j = 0; (j < nat); j++)
 375         {
 376             int                      mj              = typelist[j];
 377             const InteractionOfType& interactionType = plist[ftype].interactionTypes[ntype * mi + mj];
 378             nbsnew.emplace_back(interactionType.atoms(),
 379                                 interactionType.forceParam(),
 380                                 interactionType.interactionTypeName());
 381         }
 382         new_types.push_back(impl_->types[mi]);
 383         impl_->nameToAtomType[std::string(*impl_->types[mi].name_)] = new_types.size() - 1;
 384     }
 385
 386     mtop->ffparams.atnr = nat;
 387
 388     impl_->types                  = new_types;
 389     plist[ftype].interactionTypes = nbsnew;
 390 }
 391
 392 void PreprocessingAtomTypes::copyTot_atomtypes(t_atomtypes* atomtypes) const
 393 {
 394     /* Copy the atomtype data to the topology atomtype list */
 395     int ntype     = size();
 396     atomtypes->nr = ntype;
 397     snew(atomtypes->atomnumber, ntype);
 398
 399     for (int i = 0; i < ntype; i++)
 400     {
 401         atomtypes->atomnumber[i] = impl_->types[i].atomNumber_;
 402     }
 403 }