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,2019, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
39 #include "gpp_atomtype.h"
47 #include "gromacs/gmxpreprocess/grompp_impl.h"
48 #include "gromacs/gmxpreprocess/notset.h"
49 #include "gromacs/gmxpreprocess/topdirs.h"
50 #include "gromacs/gmxpreprocess/toputil.h"
51 #include "gromacs/math/vecdump.h"
52 #include "gromacs/topology/ifunc.h"
53 #include "gromacs/topology/symtab.h"
54 #include "gromacs/topology/topology.h"
55 #include "gromacs/utility/fatalerror.h"
56 #include "gromacs/utility/smalloc.h"
60 //! Explicit constructor.
61 AtomTypeData(const t_atom &a,
63 const InteractionType &nb,
64 const int bondAtomType,
65 const int atomNumber) :
66 atom_(a), name_(name), nb_(nb),
67 bondAtomType_(bondAtomType),
68 atomNumber_(atomNumber)
77 //! Bonded atomtype for the type.
79 //! Atom number for the atom type.
83 class PreprocessingAtomTypes::Impl
86 //! The number for currently loaded entries.
87 size_t size() const { return types.size(); }
88 //! The actual atom type data.
89 std::vector<AtomTypeData> types;
92 bool PreprocessingAtomTypes::isSet(int nt) const
94 return ((nt >= 0) || (nt < gmx::ssize(*this)));
97 int PreprocessingAtomTypes::atomTypeFromName(const std::string &str) const
99 /* Atom types are always case sensitive */
100 auto found = std::find_if(impl_->types.begin(), impl_->types.end(),
101 [&str](const auto &type)
102 { return str == *type.name_; });
103 if (found == impl_->types.end())
109 return std::distance(impl_->types.begin(), found);
113 size_t PreprocessingAtomTypes::size() const
115 return impl_->size();
118 const char *PreprocessingAtomTypes::atomNameFromAtomType(int nt) const
120 return isSet(nt) ? *(impl_->types[nt].name_) : nullptr;
123 real PreprocessingAtomTypes::atomMassAFromAtomType(int nt) const
125 return isSet(nt) ? impl_->types[nt].atom_.m : NOTSET;
128 real PreprocessingAtomTypes::atomMassBFromAtomType(int nt) const
130 return isSet(nt) ? impl_->types[nt].atom_.mB : NOTSET;
133 real PreprocessingAtomTypes::atomChargeAFromAtomType(int nt) const
135 return isSet(nt) ? impl_->types[nt].atom_.q : NOTSET;
138 real PreprocessingAtomTypes::atomChargeBFromAtomType(int nt) const
140 return isSet(nt) ? impl_->types[nt].atom_.qB : NOTSET;
143 int PreprocessingAtomTypes::atomParticleTypeFromAtomType(int nt) const
145 return isSet(nt) ? impl_->types[nt].atom_.ptype : NOTSET;
148 int PreprocessingAtomTypes::bondAtomTypeFromAtomType(int nt) const
150 return isSet(nt) ? impl_->types[nt].bondAtomType_ : NOTSET;
153 int PreprocessingAtomTypes::atomNumberFromAtomType(int nt) const
155 return isSet(nt) ? impl_->types[nt].atomNumber_ : NOTSET;
158 real PreprocessingAtomTypes::atomNonBondedParamFromAtomType(int nt, int param) const
164 gmx::ArrayRef<const real> forceParam = impl_->types[nt].nb_.forceParam();
165 if ((param < 0) || (param >= MAXFORCEPARAM))
169 return forceParam[param];
172 PreprocessingAtomTypes::PreprocessingAtomTypes()
176 PreprocessingAtomTypes::PreprocessingAtomTypes(PreprocessingAtomTypes &&old) noexcept
177 : impl_(std::move(old.impl_))
180 PreprocessingAtomTypes &PreprocessingAtomTypes::operator=(PreprocessingAtomTypes &&old) noexcept
182 impl_ = std::move(old.impl_);
186 PreprocessingAtomTypes::~PreprocessingAtomTypes()
189 int PreprocessingAtomTypes::addType(t_symtab *tab,
192 const InteractionType &nb,
196 auto found = std::find_if(impl_->types.begin(), impl_->types.end(),
197 [&name](const AtomTypeData &data)
198 { return strcmp(name, *data.name_) == 0; });
200 if (found == impl_->types.end())
202 impl_->types.emplace_back(a,
203 put_symtab(tab, name),
211 return std::distance(impl_->types.begin(), found);
215 int PreprocessingAtomTypes::setType(int nt,
219 const InteractionType &nb,
228 impl_->types[nt].atom_ = a;
229 impl_->types[nt].name_ = put_symtab(tab, name);
230 impl_->types[nt].nb_ = nb;
231 impl_->types[nt].bondAtomType_ = bondAtomType;
232 impl_->types[nt].atomNumber_ = atomNumber;
237 void PreprocessingAtomTypes::printTypes(FILE * out)
239 fprintf (out, "[ %s ]\n", dir2str(Directive::d_atomtypes));
240 fprintf (out, "; %6s %8s %8s %8s %12s %12s\n",
241 "type", "mass", "charge", "particle", "c6", "c12");
242 for (auto &entry : impl_->types)
244 fprintf(out, "%8s %8.3f %8.3f %8s %12e %12e\n",
245 *(entry.name_), entry.atom_.m, entry.atom_.q, "A",
246 entry.nb_.c0(), entry.nb_.c1());
252 static int search_atomtypes(const PreprocessingAtomTypes *ga,
254 gmx::ArrayRef<int> typelist,
256 gmx::ArrayRef<const InteractionType> interactionTypes,
260 int nrfp = NRFP(ftype);
261 int ntype = ga->size();
264 for (i = 0; (i < nn); i++)
266 if (typelist[i] == thistype)
268 /* This type number has already been added */
272 /* Otherwise, check if the parameters are identical to any previously added type */
275 for (int j = 0; j < ntype && bFound; j++)
277 /* Check nonbonded parameters */
278 gmx::ArrayRef<const real> forceParam1 = interactionTypes[ntype*typelist[i]+j].forceParam();
279 gmx::ArrayRef<const real> forceParam2 = interactionTypes[ntype*thistype+j].forceParam();
280 for (int k = 0; (k < nrfp) && bFound; k++)
282 bFound = forceParam1[k] == forceParam2[k];
285 /* Check atomnumber */
286 int tli = typelist[i];
288 (ga->atomNumberFromAtomType(tli) == ga->atomNumberFromAtomType(thistype));
300 gmx_fatal(FARGS, "Atomtype horror n = %d, %s, %d", nn, __FILE__, __LINE__);
302 typelist[nn] = thistype;
310 void PreprocessingAtomTypes::renumberTypes(gmx::ArrayRef<InteractionTypeParameters> plist,
315 int nat, ftype, ntype;
318 std::vector<int> typelist(ntype);
322 fprintf(stderr, "renumbering atomtypes...\n");
325 /* Since the bonded interactions have been assigned now,
326 * we want to reduce the number of atom types by merging
327 * ones with identical nonbonded interactions, in addition
328 * to removing unused ones.
330 * With QM/MM we also check that the atom numbers match
333 /* Get nonbonded interaction type */
334 if (plist[F_LJ].size() > 0)
343 /* Renumber atomtypes by first making a list of which ones are actually used.
344 * We provide the list of nonbonded parameters so search_atomtypes
345 * can determine if two types should be merged.
348 for (const gmx_moltype_t &moltype : mtop->moltype)
350 const t_atoms *atoms = &moltype.atoms;
351 for (int i = 0; (i < atoms->nr); i++)
353 atoms->atom[i].type =
354 search_atomtypes(this, &nat, typelist, atoms->atom[i].type,
355 plist[ftype].interactionTypes, ftype);
356 atoms->atom[i].typeB =
357 search_atomtypes(this, &nat, typelist, atoms->atom[i].typeB,
358 plist[ftype].interactionTypes, ftype);
362 for (int i = 0; i < 2; i++)
364 if (wall_atomtype[i] >= 0)
366 wall_atomtype[i] = search_atomtypes(this, &nat, typelist, wall_atomtype[i],
367 plist[ftype].interactionTypes, ftype);
371 std::vector<AtomTypeData> new_types;
372 /* We now have a list of unique atomtypes in typelist */
375 std::vector<InteractionType> nbsnew;
377 for (int i = 0; (i < nat); i++)
379 int mi = typelist[i];
380 for (int j = 0; (j < nat); j++)
382 int mj = typelist[j];
383 const InteractionType &interactionType = plist[ftype].interactionTypes[ntype*mi+mj];
384 nbsnew.emplace_back(interactionType.atoms(), interactionType.forceParam(), interactionType.interactionTypeName());
386 new_types.push_back(impl_->types[mi]);
389 mtop->ffparams.atnr = nat;
391 impl_->types = new_types;
392 plist[ftype].interactionTypes = nbsnew;
395 void PreprocessingAtomTypes::copyTot_atomtypes(t_atomtypes *atomtypes) const
397 /* Copy the atomtype data to the topology atomtype list */
399 atomtypes->nr = ntype;
400 snew(atomtypes->atomnumber, ntype);
402 for (int i = 0; i < ntype; i++)
404 atomtypes->atomnumber[i] = impl_->types[i].atomNumber_;