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, 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"
48 #include "gromacs/gmxpreprocess/grompp_impl.h"
49 #include "gromacs/gmxpreprocess/notset.h"
50 #include "gromacs/gmxpreprocess/topdirs.h"
51 #include "gromacs/gmxpreprocess/toputil.h"
52 #include "gromacs/math/vecdump.h"
53 #include "gromacs/topology/ifunc.h"
54 #include "gromacs/topology/symtab.h"
55 #include "gromacs/topology/topology.h"
56 #include "gromacs/utility/fatalerror.h"
57 #include "gromacs/utility/smalloc.h"
61 //! Explicit constructor.
62 AtomTypeData(const t_atom& a, char** name, const InteractionOfType& nb, const int bondAtomType, const int atomNumber) :
66 bondAtomType_(bondAtomType),
67 atomNumber_(atomNumber)
75 InteractionOfType nb_;
76 //! Bonded atomtype for the type.
78 //! Atom number for the atom type.
82 class PreprocessingAtomTypes::Impl
85 //! The number for currently loaded entries.
86 size_t size() const { return types.size(); }
87 //! The actual atom type data.
88 std::vector<AtomTypeData> types;
91 bool PreprocessingAtomTypes::isSet(int nt) const
93 return ((nt >= 0) && (nt < gmx::ssize(*this)));
96 int PreprocessingAtomTypes::atomTypeFromName(const std::string& str) const
98 /* Atom types are always case sensitive */
99 auto found = std::find_if(impl_->types.begin(), impl_->types.end(), [&str](const auto& type) {
100 return str == *type.name_;
102 if (found == impl_->types.end())
108 return std::distance(impl_->types.begin(), found);
112 size_t PreprocessingAtomTypes::size() const
114 return impl_->size();
117 const char* PreprocessingAtomTypes::atomNameFromAtomType(int nt) const
119 return isSet(nt) ? *(impl_->types[nt].name_) : nullptr;
122 real PreprocessingAtomTypes::atomMassFromAtomType(int nt) const
124 return isSet(nt) ? impl_->types[nt].atom_.m : NOTSET;
127 real PreprocessingAtomTypes::atomChargeFromAtomType(int nt) const
129 return isSet(nt) ? impl_->types[nt].atom_.q : NOTSET;
132 int PreprocessingAtomTypes::atomParticleTypeFromAtomType(int nt) const
134 return isSet(nt) ? impl_->types[nt].atom_.ptype : NOTSET;
137 int PreprocessingAtomTypes::bondAtomTypeFromAtomType(int nt) const
139 return isSet(nt) ? impl_->types[nt].bondAtomType_ : NOTSET;
142 int PreprocessingAtomTypes::atomNumberFromAtomType(int nt) const
144 return isSet(nt) ? impl_->types[nt].atomNumber_ : NOTSET;
147 real PreprocessingAtomTypes::atomNonBondedParamFromAtomType(int nt, int param) const
153 gmx::ArrayRef<const real> forceParam = impl_->types[nt].nb_.forceParam();
154 if ((param < 0) || (param >= MAXFORCEPARAM))
158 return forceParam[param];
161 PreprocessingAtomTypes::PreprocessingAtomTypes() : impl_(new Impl) {}
163 PreprocessingAtomTypes::PreprocessingAtomTypes(PreprocessingAtomTypes&& old) noexcept :
164 impl_(std::move(old.impl_))
168 PreprocessingAtomTypes& PreprocessingAtomTypes::operator=(PreprocessingAtomTypes&& old) noexcept
170 impl_ = std::move(old.impl_);
174 PreprocessingAtomTypes::~PreprocessingAtomTypes() {}
176 int PreprocessingAtomTypes::addType(t_symtab* tab,
178 const std::string& name,
179 const InteractionOfType& nb,
183 int position = atomTypeFromName(name);
184 if (position == NOTSET)
186 impl_->types.emplace_back(a, put_symtab(tab, name.c_str()), nb, bondAtomType, atomNumber);
187 return atomTypeFromName(name);
195 int PreprocessingAtomTypes::setType(int nt,
198 const std::string& name,
199 const InteractionOfType& nb,
208 impl_->types[nt].atom_ = a;
209 impl_->types[nt].name_ = put_symtab(tab, name.c_str());
210 impl_->types[nt].nb_ = nb;
211 impl_->types[nt].bondAtomType_ = bondAtomType;
212 impl_->types[nt].atomNumber_ = atomNumber;
217 void PreprocessingAtomTypes::printTypes(FILE* out)
219 fprintf(out, "[ %s ]\n", dir2str(Directive::d_atomtypes));
221 "; %6s %8s %8s %8s %12s %12s\n",
228 for (auto& entry : impl_->types)
231 "%8s %8.3f %8.3f %8s %12e %12e\n",
243 static int search_atomtypes(const PreprocessingAtomTypes* ga,
245 gmx::ArrayRef<int> typelist,
247 gmx::ArrayRef<const InteractionOfType> interactionTypes,
251 int nrfp = NRFP(ftype);
252 int ntype = ga->size();
255 for (i = 0; (i < nn); i++)
257 if (typelist[i] == thistype)
259 /* This type number has already been added */
263 /* Otherwise, check if the parameters are identical to any previously added type */
266 for (int j = 0; j < ntype && bFound; j++)
268 /* Check nonbonded parameters */
269 gmx::ArrayRef<const real> forceParam1 =
270 interactionTypes[ntype * typelist[i] + j].forceParam();
271 gmx::ArrayRef<const real> forceParam2 = interactionTypes[ntype * thistype + j].forceParam();
272 for (int k = 0; (k < nrfp) && bFound; k++)
274 bFound = forceParam1[k] == forceParam2[k];
277 /* Check atomnumber */
278 int tli = typelist[i];
279 bFound = bFound && (ga->atomNumberFromAtomType(tli) == ga->atomNumberFromAtomType(thistype));
291 gmx_fatal(FARGS, "Atomtype horror n = %d, %s, %d", nn, __FILE__, __LINE__);
293 typelist[nn] = thistype;
301 void PreprocessingAtomTypes::renumberTypes(gmx::ArrayRef<InteractionsOfType> plist,
306 int nat, ftype, ntype;
309 std::vector<int> typelist(ntype);
313 fprintf(stderr, "renumbering atomtypes...\n");
316 /* Since the bonded interactions have been assigned now,
317 * we want to reduce the number of atom types by merging
318 * ones with identical nonbonded interactions, in addition
319 * to removing unused ones.
321 * With QM/MM we also check that the atom numbers match
324 /* Get nonbonded interaction type */
325 if (plist[F_LJ].size() > 0)
334 /* Renumber atomtypes by first making a list of which ones are actually used.
335 * We provide the list of nonbonded parameters so search_atomtypes
336 * can determine if two types should be merged.
339 for (const gmx_moltype_t& moltype : mtop->moltype)
341 const t_atoms* atoms = &moltype.atoms;
342 for (int i = 0; (i < atoms->nr); i++)
344 atoms->atom[i].type = search_atomtypes(
345 this, &nat, typelist, atoms->atom[i].type, plist[ftype].interactionTypes, ftype);
346 atoms->atom[i].typeB = search_atomtypes(
347 this, &nat, typelist, atoms->atom[i].typeB, plist[ftype].interactionTypes, ftype);
351 for (int i = 0; i < 2; i++)
353 if (wall_atomtype[i] >= 0)
355 wall_atomtype[i] = search_atomtypes(
356 this, &nat, typelist, wall_atomtype[i], plist[ftype].interactionTypes, ftype);
360 std::vector<AtomTypeData> new_types;
361 /* We now have a list of unique atomtypes in typelist */
364 std::vector<InteractionOfType> nbsnew;
366 for (int i = 0; (i < nat); i++)
368 int mi = typelist[i];
369 for (int j = 0; (j < nat); j++)
371 int mj = typelist[j];
372 const InteractionOfType& interactionType = plist[ftype].interactionTypes[ntype * mi + mj];
373 nbsnew.emplace_back(interactionType.atoms(),
374 interactionType.forceParam(),
375 interactionType.interactionTypeName());
377 new_types.push_back(impl_->types[mi]);
380 mtop->ffparams.atnr = nat;
382 impl_->types = new_types;
383 plist[ftype].interactionTypes = nbsnew;
386 void PreprocessingAtomTypes::copyTot_atomtypes(t_atomtypes* atomtypes) const
388 /* Copy the atomtype data to the topology atomtype list */
390 atomtypes->nr = ntype;
391 snew(atomtypes->atomnumber, ntype);
393 for (int i = 0; i < ntype; i++)
395 atomtypes->atomnumber[i] = impl_->types[i].atomNumber_;