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_
;