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, char** name, const InteractionOfType& nb, const int bondAtomType, const int atomNumber) :
65 bondAtomType_(bondAtomType),
66 atomNumber_(atomNumber)
74 InteractionOfType nb_;
75 //! Bonded atomtype for the type.
77 //! Atom number for the atom type.
81 class PreprocessingAtomTypes::Impl
84 //! The number for currently loaded entries.
85 size_t size() const { return types.size(); }
86 //! The actual atom type data.
87 std::vector<AtomTypeData> types;
90 bool PreprocessingAtomTypes::isSet(int nt) const
92 return ((nt >= 0) && (nt < gmx::ssize(*this)));
95 int PreprocessingAtomTypes::atomTypeFromName(const std::string& str) const
97 /* Atom types are always case sensitive */
98 auto found = std::find_if(impl_->types.begin(), impl_->types.end(),
99 [&str](const auto& type) { return str == *type.name_; });
100 if (found == impl_->types.end())
106 return std::distance(impl_->types.begin(), found);
110 size_t PreprocessingAtomTypes::size() const
112 return impl_->size();
115 const char* PreprocessingAtomTypes::atomNameFromAtomType(int nt) const
117 return isSet(nt) ? *(impl_->types[nt].name_) : nullptr;
120 real PreprocessingAtomTypes::atomMassFromAtomType(int nt) const
122 return isSet(nt) ? impl_->types[nt].atom_.m : NOTSET;
125 real PreprocessingAtomTypes::atomChargeFromAtomType(int nt) const
127 return isSet(nt) ? impl_->types[nt].atom_.q : NOTSET;
130 int PreprocessingAtomTypes::atomParticleTypeFromAtomType(int nt) const
132 return isSet(nt) ? impl_->types[nt].atom_.ptype : NOTSET;
135 int PreprocessingAtomTypes::bondAtomTypeFromAtomType(int nt) const
137 return isSet(nt) ? impl_->types[nt].bondAtomType_ : NOTSET;
140 int PreprocessingAtomTypes::atomNumberFromAtomType(int nt) const
142 return isSet(nt) ? impl_->types[nt].atomNumber_ : NOTSET;
145 real PreprocessingAtomTypes::atomNonBondedParamFromAtomType(int nt, int param) const
151 gmx::ArrayRef<const real> forceParam = impl_->types[nt].nb_.forceParam();
152 if ((param < 0) || (param >= MAXFORCEPARAM))
156 return forceParam[param];
159 PreprocessingAtomTypes::PreprocessingAtomTypes() : impl_(new Impl) {}
161 PreprocessingAtomTypes::PreprocessingAtomTypes(PreprocessingAtomTypes&& old) noexcept :
162 impl_(std::move(old.impl_))
166 PreprocessingAtomTypes& PreprocessingAtomTypes::operator=(PreprocessingAtomTypes&& old) noexcept
168 impl_ = std::move(old.impl_);
172 PreprocessingAtomTypes::~PreprocessingAtomTypes() {}
174 int PreprocessingAtomTypes::addType(t_symtab* tab,
176 const std::string& name,
177 const InteractionOfType& nb,
181 int position = atomTypeFromName(name);
182 if (position == NOTSET)
184 impl_->types.emplace_back(a, put_symtab(tab, name.c_str()), nb, bondAtomType, atomNumber);
185 return atomTypeFromName(name);
193 int PreprocessingAtomTypes::setType(int nt,
196 const std::string& name,
197 const InteractionOfType& nb,
206 impl_->types[nt].atom_ = a;
207 impl_->types[nt].name_ = put_symtab(tab, name.c_str());
208 impl_->types[nt].nb_ = nb;
209 impl_->types[nt].bondAtomType_ = bondAtomType;
210 impl_->types[nt].atomNumber_ = atomNumber;
215 void PreprocessingAtomTypes::printTypes(FILE* out)
217 fprintf(out, "[ %s ]\n", dir2str(Directive::d_atomtypes));
218 fprintf(out, "; %6s %8s %8s %8s %12s %12s\n", "type", "mass", "charge", "particle", "c6",
220 for (auto& entry : impl_->types)
222 fprintf(out, "%8s %8.3f %8.3f %8s %12e %12e\n", *(entry.name_), entry.atom_.m,
223 entry.atom_.q, "A", entry.nb_.c0(), entry.nb_.c1());
229 static int search_atomtypes(const PreprocessingAtomTypes* ga,
231 gmx::ArrayRef<int> typelist,
233 gmx::ArrayRef<const InteractionOfType> interactionTypes,
237 int nrfp = NRFP(ftype);
238 int ntype = ga->size();
241 for (i = 0; (i < nn); i++)
243 if (typelist[i] == thistype)
245 /* This type number has already been added */
249 /* Otherwise, check if the parameters are identical to any previously added type */
252 for (int j = 0; j < ntype && bFound; j++)
254 /* Check nonbonded parameters */
255 gmx::ArrayRef<const real> forceParam1 =
256 interactionTypes[ntype * typelist[i] + j].forceParam();
257 gmx::ArrayRef<const real> forceParam2 = interactionTypes[ntype * thistype + j].forceParam();
258 for (int k = 0; (k < nrfp) && bFound; k++)
260 bFound = forceParam1[k] == forceParam2[k];
263 /* Check atomnumber */
264 int tli = typelist[i];
265 bFound = bFound && (ga->atomNumberFromAtomType(tli) == ga->atomNumberFromAtomType(thistype));
277 gmx_fatal(FARGS, "Atomtype horror n = %d, %s, %d", nn, __FILE__, __LINE__);
279 typelist[nn] = thistype;
287 void PreprocessingAtomTypes::renumberTypes(gmx::ArrayRef<InteractionsOfType> plist,
292 int nat, ftype, ntype;
295 std::vector<int> typelist(ntype);
299 fprintf(stderr, "renumbering atomtypes...\n");
302 /* Since the bonded interactions have been assigned now,
303 * we want to reduce the number of atom types by merging
304 * ones with identical nonbonded interactions, in addition
305 * to removing unused ones.
307 * With QM/MM we also check that the atom numbers match
310 /* Get nonbonded interaction type */
311 if (plist[F_LJ].size() > 0)
320 /* Renumber atomtypes by first making a list of which ones are actually used.
321 * We provide the list of nonbonded parameters so search_atomtypes
322 * can determine if two types should be merged.
325 for (const gmx_moltype_t& moltype : mtop->moltype)
327 const t_atoms* atoms = &moltype.atoms;
328 for (int i = 0; (i < atoms->nr); i++)
330 atoms->atom[i].type = search_atomtypes(this, &nat, typelist, atoms->atom[i].type,
331 plist[ftype].interactionTypes, ftype);
332 atoms->atom[i].typeB = search_atomtypes(this, &nat, typelist, atoms->atom[i].typeB,
333 plist[ftype].interactionTypes, ftype);
337 for (int i = 0; i < 2; i++)
339 if (wall_atomtype[i] >= 0)
341 wall_atomtype[i] = search_atomtypes(this, &nat, typelist, wall_atomtype[i],
342 plist[ftype].interactionTypes, ftype);
346 std::vector<AtomTypeData> new_types;
347 /* We now have a list of unique atomtypes in typelist */
350 std::vector<InteractionOfType> nbsnew;
352 for (int i = 0; (i < nat); i++)
354 int mi = typelist[i];
355 for (int j = 0; (j < nat); j++)
357 int mj = typelist[j];
358 const InteractionOfType& interactionType = plist[ftype].interactionTypes[ntype * mi + mj];
359 nbsnew.emplace_back(interactionType.atoms(), interactionType.forceParam(),
360 interactionType.interactionTypeName());
362 new_types.push_back(impl_->types[mi]);
365 mtop->ffparams.atnr = nat;
367 impl_->types = new_types;
368 plist[ftype].interactionTypes = nbsnew;
371 void PreprocessingAtomTypes::copyTot_atomtypes(t_atomtypes* atomtypes) const
373 /* Copy the atomtype data to the topology atomtype list */
375 atomtypes->nr = ntype;
376 snew(atomtypes->atomnumber, ntype);
378 for (int i = 0; i < ntype; i++)
380 atomtypes->atomnumber[i] = impl_->types[i].atomNumber_;