2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
40 #include "gpp_atomtype.h"
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"
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)
72 InteractionOfType nb_;
73 //! Bonded atomtype for the type.
75 //! Atom number for the atom type.
79 class PreprocessingAtomTypes::Impl
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;
90 bool PreprocessingAtomTypes::isSet(int nt) const
92 return ((nt >= 0) && (nt < gmx::ssize(*this)));
95 std::optional<int> PreprocessingAtomTypes::atomTypeFromName(const std::string& str) const
97 /* Atom types are always case sensitive */
98 const auto found = impl_->nameToAtomType.find(str);
99 if (found == impl_->nameToAtomType.end())
105 GMX_ASSERT(str == *impl_->types[found->second].name_,
106 "Invalid data in atomTypeFromName lookup table");
107 return std::make_optional(found->second);
111 size_t PreprocessingAtomTypes::size() const
113 return impl_->size();
116 std::optional<const char*> PreprocessingAtomTypes::atomNameFromAtomType(int nt) const
118 return isSet(nt) ? std::make_optional(*(impl_->types[nt].name_)) : std::nullopt;
121 std::optional<real> PreprocessingAtomTypes::atomMassFromAtomType(int nt) const
123 return isSet(nt) ? std::make_optional(impl_->types[nt].atom_.m) : std::nullopt;
126 std::optional<real> PreprocessingAtomTypes::atomChargeFromAtomType(int nt) const
128 return isSet(nt) ? std::make_optional(impl_->types[nt].atom_.q) : std::nullopt;
131 std::optional<ParticleType> PreprocessingAtomTypes::atomParticleTypeFromAtomType(int nt) const
133 return isSet(nt) ? std::make_optional(impl_->types[nt].atom_.ptype) : std::nullopt;
136 std::optional<int> PreprocessingAtomTypes::bondAtomTypeFromAtomType(int nt) const
138 return isSet(nt) ? std::make_optional(impl_->types[nt].bondAtomType_) : std::nullopt;
141 std::optional<int> PreprocessingAtomTypes::atomNumberFromAtomType(int nt) const
143 return isSet(nt) ? std::make_optional(impl_->types[nt].atomNumber_) : std::nullopt;
146 std::optional<real> PreprocessingAtomTypes::atomNonBondedParamFromAtomType(int nt, int param) const
152 gmx::ArrayRef<const real> forceParam = impl_->types[nt].nb_.forceParam();
153 if ((param < 0) || (param >= MAXFORCEPARAM))
157 return std::make_optional(forceParam[param]);
160 PreprocessingAtomTypes::PreprocessingAtomTypes() : impl_(new Impl) {}
162 PreprocessingAtomTypes::PreprocessingAtomTypes(PreprocessingAtomTypes&& old) noexcept :
163 impl_(std::move(old.impl_))
167 PreprocessingAtomTypes& PreprocessingAtomTypes::operator=(PreprocessingAtomTypes&& old) noexcept
169 impl_ = std::move(old.impl_);
173 PreprocessingAtomTypes::~PreprocessingAtomTypes() {}
175 int PreprocessingAtomTypes::addType(t_symtab* tab,
177 const std::string& name,
178 const InteractionOfType& nb,
182 auto position = atomTypeFromName(name);
183 if (!position.has_value())
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;
196 std::optional<int> PreprocessingAtomTypes::setType(int nt,
199 const std::string& name,
200 const InteractionOfType& nb,
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;
215 return std::make_optional(nt);
218 void PreprocessingAtomTypes::printTypes(FILE* out)
220 fprintf(out, "[ %s ]\n", enumValueToString(Directive::d_atomtypes));
222 "; %6s %8s %8s %8s %12s %12s\n",
229 for (auto& entry : impl_->types)
232 "%8s %8.3f %8.3f %8s %12e %12e\n",
244 static int search_atomtypes(const PreprocessingAtomTypes* ga,
246 gmx::ArrayRef<int> typelist,
248 gmx::ArrayRef<const InteractionOfType> interactionTypes,
252 int nrfp = NRFP(ftype);
253 int ntype = ga->size();
256 for (i = 0; (i < nn); i++)
258 if (typelist[i] == thistype)
260 /* This type number has already been added */
264 /* Otherwise, check if the parameters are identical to any previously added type */
267 for (int j = 0; j < ntype && bFound; j++)
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++)
275 bFound = forceParam1[k] == forceParam2[k];
278 /* Check atomnumber */
279 int tli = typelist[i];
281 && (*ga->atomNumberFromAtomType(tli) == *ga->atomNumberFromAtomType(thistype));
293 gmx_fatal(FARGS, "Atomtype horror n = %d, %s, %d", nn, __FILE__, __LINE__);
295 typelist[nn] = thistype;
303 void PreprocessingAtomTypes::renumberTypes(gmx::ArrayRef<InteractionsOfType> plist,
308 int nat, ftype, ntype;
311 std::vector<int> typelist(ntype);
315 fprintf(stderr, "renumbering atomtypes...\n");
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.
323 * With QM/MM we also check that the atom numbers match
326 /* Get nonbonded interaction type */
327 if (plist[F_LJ].size() > 0)
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.
341 for (const gmx_moltype_t& moltype : mtop->moltype)
343 const t_atoms* atoms = &moltype.atoms;
344 for (int i = 0; (i < atoms->nr); i++)
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);
353 for (int i = 0; i < 2; i++)
355 if (wall_atomtype[i] >= 0)
357 wall_atomtype[i] = search_atomtypes(
358 this, &nat, typelist, wall_atomtype[i], plist[ftype].interactionTypes, ftype);
362 std::vector<AtomTypeData> new_types;
363 /* We now have a list of unique atomtypes in typelist */
366 std::vector<InteractionOfType> nbsnew;
368 // Reset the map used for fast lookups, and refill it below
369 impl_->nameToAtomType.clear();
371 for (int i = 0; (i < nat); i++)
373 int mi = typelist[i];
374 for (int j = 0; (j < nat); j++)
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());
382 new_types.push_back(impl_->types[mi]);
383 impl_->nameToAtomType[std::string(*impl_->types[mi].name_)] = new_types.size() - 1;
386 mtop->ffparams.atnr = nat;
388 impl_->types = new_types;
389 plist[ftype].interactionTypes = nbsnew;
392 void PreprocessingAtomTypes::copyTot_atomtypes(t_atomtypes* atomtypes) const
394 /* Copy the atomtype data to the topology atomtype list */
396 atomtypes->nr = ntype;
397 snew(atomtypes->atomnumber, ntype);
399 for (int i = 0; i < ntype; i++)
401 atomtypes->atomnumber[i] = impl_->types[i].atomNumber_;