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) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,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.
51 #include "gromacs/fileio/warninp.h"
52 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
53 #include "gromacs/gmxpreprocess/gpp_bond_atomtype.h"
54 #include "gromacs/gmxpreprocess/grompp_impl.h"
55 #include "gromacs/gmxpreprocess/notset.h"
56 #include "gromacs/gmxpreprocess/readir.h"
57 #include "gromacs/gmxpreprocess/topdirs.h"
58 #include "gromacs/gmxpreprocess/toputil.h"
59 #include "gromacs/math/functions.h"
60 #include "gromacs/mdtypes/md_enums.h"
61 #include "gromacs/topology/exclusionblocks.h"
62 #include "gromacs/topology/ifunc.h"
63 #include "gromacs/topology/symtab.h"
64 #include "gromacs/utility/arrayref.h"
65 #include "gromacs/utility/cstringutil.h"
66 #include "gromacs/utility/enumerationhelpers.h"
67 #include "gromacs/utility/fatalerror.h"
68 #include "gromacs/utility/gmxassert.h"
69 #include "gromacs/utility/smalloc.h"
70 #include "gromacs/utility/stringtoenumvalueconverter.h"
71 #include "gromacs/utility/stringutil.h"
73 void generate_nbparams(CombinationRule comb,
75 InteractionsOfType* interactions,
76 PreprocessingAtomTypes* atypes,
79 constexpr int c_nrfp2 = 2;
82 real c, bi, bj, ci, cj, ci0, ci1, ci2, cj0, cj1, cj2;
84 /* Lean mean shortcuts */
87 interactions->interactionTypes.clear();
89 std::array<real, MAXFORCEPARAM> forceParam = { NOTSET };
90 /* Fill the matrix with force parameters */
91 // Prefetch the parameters to improve cache hits and avoid dereference and call overhead
92 std::vector<std::pair<real, real>> cPrefetch;
93 cPrefetch.reserve(nr);
94 for (int i = 0; i < nr; i++)
96 cPrefetch.emplace_back(*atypes->atomNonBondedParamFromAtomType(i, 0),
97 *atypes->atomNonBondedParamFromAtomType(i, 1));
99 interactions->interactionTypes.reserve(nr * nr);
105 case CombinationRule::Geometric:
106 // Geometric combination rules, c6 and c12 are independent
107 GMX_RELEASE_ASSERT(nrfp == c_nrfp2, "nfrp should be 2");
108 for (int i = 0; (i < nr); i++)
110 for (int j = 0; (j < nr); j++)
112 for (int nf = 0; (nf < c_nrfp2); nf++)
114 ci = (nf == 0 ? cPrefetch[i].first : cPrefetch[i].second);
115 cj = (nf == 0 ? cPrefetch[j].first : cPrefetch[j].second);
116 c = std::sqrt(ci * cj);
119 interactions->interactionTypes.emplace_back(gmx::ArrayRef<const int>{},
125 case CombinationRule::Arithmetic:
126 /* c0 and c1 are sigma and epsilon */
127 for (int i = 0; (i < nr); i++)
129 for (int j = 0; (j < nr); j++)
131 ci0 = cPrefetch[i].first;
132 cj0 = cPrefetch[j].first;
133 ci1 = cPrefetch[i].second;
134 cj1 = cPrefetch[j].second;
135 forceParam[0] = (fabs(ci0) + fabs(cj0)) * 0.5;
136 /* Negative sigma signals that c6 should be set to zero later,
137 * so we need to propagate that through the combination rules.
139 if (ci0 < 0 || cj0 < 0)
143 forceParam[1] = std::sqrt(ci1 * cj1);
144 interactions->interactionTypes.emplace_back(gmx::ArrayRef<const int>{},
150 case CombinationRule::GeomSigEps:
151 /* c0 and c1 are sigma and epsilon */
152 for (int i = 0; (i < nr); i++)
154 for (int j = 0; (j < nr); j++)
156 ci0 = cPrefetch[i].first;
157 cj0 = cPrefetch[j].first;
158 ci1 = cPrefetch[i].second;
159 cj1 = cPrefetch[j].second;
160 forceParam[0] = std::sqrt(std::fabs(ci0 * cj0));
161 /* Negative sigma signals that c6 should be set to zero later,
162 * so we need to propagate that through the combination rules.
164 if (ci0 < 0 || cj0 < 0)
168 forceParam[1] = std::sqrt(ci1 * cj1);
169 interactions->interactionTypes.emplace_back(gmx::ArrayRef<const int>{},
177 gmx::formatString("No such combination rule %s", enumValueToString(comb));
178 warning_error_and_exit(wi, message, FARGS);
183 /* Buckingham rules */
184 for (int i = 0; (i < nr); i++)
186 for (int j = 0; (j < nr); j++)
188 ci0 = cPrefetch[i].first;
189 cj0 = cPrefetch[j].first;
190 ci2 = *atypes->atomNonBondedParamFromAtomType(i, 2);
191 cj2 = *atypes->atomNonBondedParamFromAtomType(j, 2);
192 bi = cPrefetch[i].second;
193 bj = cPrefetch[j].second;
194 forceParam[0] = std::sqrt(ci0 * cj0);
195 if ((bi == 0) || (bj == 0))
201 forceParam[1] = 2.0 / (1 / bi + 1 / bj);
203 forceParam[2] = std::sqrt(ci2 * cj2);
204 interactions->interactionTypes.emplace_back(gmx::ArrayRef<const int>{}, forceParam);
210 auto message = gmx::formatString("Invalid nonbonded type %s",
211 interaction_function[ftype].longname);
212 warning_error(wi, message);
216 /*! \brief Used to temporarily store the explicit non-bonded parameter
217 * combinations, which will be copied to InteractionsOfType. */
220 //! Has this combination been set.
222 //! The non-bonded parameters
226 static void realloc_nb_params(PreprocessingAtomTypes* atypes, t_nbparam*** nbparam, t_nbparam*** pair)
228 /* Add space in the non-bonded parameters matrix */
229 int atnr = atypes->size();
230 srenew(*nbparam, atnr);
231 snew((*nbparam)[atnr - 1], atnr);
235 snew((*pair)[atnr - 1], atnr);
239 int copy_nbparams(t_nbparam** param, int ftype, InteractionsOfType* interactions, int nr)
246 for (int i = 0; i < nr; i++)
248 for (int j = 0; j <= i; j++)
250 GMX_RELEASE_ASSERT(param, "Must have valid parameters");
251 if (param[i][j].bSet)
253 for (int f = 0; f < nrfp; f++)
255 interactions->interactionTypes[nr * i + j].setForceParameter(f, param[i][j].c[f]);
256 interactions->interactionTypes[nr * j + i].setForceParameter(f, param[i][j].c[f]);
266 void free_nbparam(t_nbparam** param, int nr)
270 GMX_RELEASE_ASSERT(param, "Must have valid parameters");
271 for (i = 0; i < nr; i++)
273 GMX_RELEASE_ASSERT(param[i], "Must have valid parameters");
279 static void copy_B_from_A(int ftype, double* c)
283 nrfpA = NRFPA(ftype);
284 nrfpB = NRFPB(ftype);
286 /* Copy the B parameters from the first nrfpB A parameters */
287 for (i = 0; (i < nrfpB); i++)
293 //! Local definition that supersedes the central one, as we only want the leading letter
294 static const char* enumValueToLetterAsString(ParticleType enumValue)
296 static constexpr gmx::EnumerationArray<ParticleType, const char*> particleTypeLetters = {
297 "A", "N", "S", "B", "V"
299 return particleTypeLetters[enumValue];
302 void push_at(t_symtab* symtab,
303 PreprocessingAtomTypes* at,
304 PreprocessingBondAtomType* bondAtomType,
307 t_nbparam*** nbparam,
311 int nfields, nfp0 = -1;
313 char type[STRLEN], btype[STRLEN], ptype[STRLEN];
315 double c[MAXFORCEPARAM];
316 char tmpfield[12][100]; /* Max 12 fields of width 100 */
319 bool have_atomic_number;
320 bool have_bonded_type;
324 /* First assign input line to temporary array */
325 nfields = sscanf(line,
326 "%s%s%s%s%s%s%s%s%s%s%s%s",
340 /* Comments on optional fields in the atomtypes section:
342 * The force field format is getting a bit old. For OPLS-AA we needed
343 * to add a special bonded atomtype, and for Gerrit Groenhofs QM/MM stuff
344 * we also needed the atomic numbers.
345 * To avoid making all old or user-generated force fields unusable we
346 * have introduced both these quantities as optional columns, and do some
347 * acrobatics to check whether they are present or not.
348 * This will all look much nicer when we switch to XML... sigh.
350 * Field 0 (mandatory) is the nonbonded type name. (string)
351 * Field 1 (optional) is the bonded type (string)
352 * Field 2 (optional) is the atomic number (int)
353 * Field 3 (mandatory) is the mass (numerical)
354 * Field 4 (mandatory) is the charge (numerical)
355 * Field 5 (mandatory) is the particle type (single character)
356 * This is followed by a number of nonbonded parameters.
358 * The safest way to identify the format is the particle type field.
360 * So, here is what we do:
362 * A. Read in the first six fields as strings
363 * B. If field 3 (starting from 0) is a single char, we have neither
364 * bonded_type or atomic numbers.
365 * C. If field 5 is a single char we have both.
366 * D. If field 4 is a single char we check field 1. If this begins with
367 * an alphabetical character we have bonded types, otherwise atomic numbers.
376 if ((strlen(tmpfield[5]) == 1) && isalpha(tmpfield[5][0]))
378 have_bonded_type = TRUE;
379 have_atomic_number = TRUE;
381 else if ((strlen(tmpfield[3]) == 1) && isalpha(tmpfield[3][0]))
383 have_bonded_type = FALSE;
384 have_atomic_number = FALSE;
388 have_bonded_type = (isalpha(tmpfield[1][0]) != 0);
389 have_atomic_number = !have_bonded_type;
392 /* optional fields */
401 if (have_atomic_number)
403 if (have_bonded_type)
406 line, "%s%s%d%lf%lf%s%lf%lf", type, btype, &atomnr, &m, &q, ptype, &c[0], &c[1]);
415 /* have_atomic_number && !have_bonded_type */
416 nread = sscanf(line, "%s%d%lf%lf%s%lf%lf", type, &atomnr, &m, &q, ptype, &c[0], &c[1]);
426 if (have_bonded_type)
428 /* !have_atomic_number && have_bonded_type */
429 nread = sscanf(line, "%s%s%lf%lf%s%lf%lf", type, btype, &m, &q, ptype, &c[0], &c[1]);
438 /* !have_atomic_number && !have_bonded_type */
439 nread = sscanf(line, "%s%lf%lf%s%lf%lf", type, &m, &q, ptype, &c[0], &c[1]);
448 if (!have_bonded_type)
453 if (!have_atomic_number)
463 if (have_atomic_number)
465 if (have_bonded_type)
468 line, "%s%s%d%lf%lf%s%lf%lf%lf", type, btype, &atomnr, &m, &q, ptype, &c[0], &c[1], &c[2]);
477 /* have_atomic_number && !have_bonded_type */
479 line, "%s%d%lf%lf%s%lf%lf%lf", type, &atomnr, &m, &q, ptype, &c[0], &c[1], &c[2]);
489 if (have_bonded_type)
491 /* !have_atomic_number && have_bonded_type */
493 line, "%s%s%lf%lf%s%lf%lf%lf", type, btype, &m, &q, ptype, &c[0], &c[1], &c[2]);
502 /* !have_atomic_number && !have_bonded_type */
503 nread = sscanf(line, "%s%lf%lf%s%lf%lf%lf", type, &m, &q, ptype, &c[0], &c[1], &c[2]);
512 if (!have_bonded_type)
517 if (!have_atomic_number)
525 auto message = gmx::formatString("Invalid function type %d in push_at", nb_funct);
526 warning_error_and_exit(wi, message, FARGS);
528 for (int j = nfp0; (j < MAXFORCEPARAM); j++)
532 std::array<real, MAXFORCEPARAM> forceParam;
534 if (strlen(type) == 1 && isdigit(type[0]))
536 warning_error_and_exit(wi, "Atom type names can't be single digits.", FARGS);
539 if (strlen(btype) == 1 && isdigit(btype[0]))
541 warning_error_and_exit(wi, "Bond atom type names can't be single digits.", FARGS);
544 /* Hack to read old topologies */
545 if (gmx_strcasecmp(ptype, "D") == 0)
549 static const gmx::StringToEnumValueConverter<ParticleType, enumValueToLetterAsString, gmx::StringCompareType::CaseInsensitive, gmx::StripStrings::No>
550 s_stringToParticleType;
551 std::optional<ParticleType> pt = s_stringToParticleType.valueFrom(ptype);
554 auto message = gmx::formatString("Invalid particle type %s", ptype);
555 warning_error_and_exit(wi, message, FARGS);
560 atom->ptype = pt.value();
561 for (int i = 0; i < MAXFORCEPARAM; i++)
563 forceParam[i] = c[i];
566 InteractionOfType interactionType({}, forceParam, "");
568 auto batype_nr = bondAtomType->addBondAtomType(symtab, btype);
570 auto atomType = at->atomTypeFromName(type);
571 if (atomType.has_value())
573 auto message = gmx::formatString(
574 "Atomtype %s was defined previously (e.g. in the forcefield files), "
575 "and has now been defined again. This could happen e.g. if you would "
576 "use a self-contained molecule .itp file that duplicates or replaces "
577 "the contents of the standard force-field files. You should check "
578 "the contents of your files and remove such repetition. If you know "
579 "you should override the previous definition, then you could choose "
580 "to suppress this warning with -maxwarn.",
582 warning(wi, message);
583 auto newAtomType = at->setType(*atomType, symtab, *atom, type, interactionType, batype_nr, atomnr);
584 if (!newAtomType.has_value())
586 auto message = gmx::formatString("Replacing atomtype %s failed", type);
587 warning_error_and_exit(wi, message, FARGS);
592 at->addType(symtab, *atom, type, interactionType, batype_nr, atomnr);
593 /* Add space in the non-bonded parameters matrix */
594 realloc_nb_params(at, nbparam, pair);
599 //! Return whether the contents of \c a and \c b are the same, considering also reversed order.
601 static bool equalEitherForwardOrBackward(gmx::ArrayRef<const T> a, gmx::ArrayRef<const T> b)
603 return (std::equal(a.begin(), a.end(), b.begin()) || std::equal(a.begin(), a.end(), b.rbegin()));
606 static void push_bondtype(InteractionsOfType* bt,
607 const InteractionOfType& b,
615 int nrfp = NRFP(ftype);
617 /* If bAllowRepeat is TRUE, we allow multiple entries as long as they
618 are on directly _adjacent_ lines.
621 /* First check if our atomtypes are _identical_ (not reversed) to the previous
622 entry. If they are not identical we search for earlier duplicates. If they are
623 we can skip it, since we already searched for the first line
627 bool isContinuationOfBlock = false;
628 if (bAllowRepeat && nr > 1)
630 isContinuationOfBlock = true;
631 gmx::ArrayRef<const int> newParAtom = b.atoms();
632 gmx::ArrayRef<const int> sysParAtom = bt->interactionTypes[nr - 2].atoms();
633 for (int j = 0; j < nral; j++)
635 if (newParAtom[j] != sysParAtom[j])
637 isContinuationOfBlock = false;
642 /* Search for earlier duplicates if this entry was not a continuation
643 from the previous line.
645 bool addBondType = true;
646 bool haveWarned = false;
647 bool haveErrored = false;
648 for (int i = 0; (i < nr); i++)
650 gmx::ArrayRef<const int> bParams = b.atoms();
651 gmx::ArrayRef<const int> testParams = bt->interactionTypes[i].atoms();
652 GMX_RELEASE_ASSERT(bParams.size() == testParams.size(),
653 "Number of atoms needs to be the same between parameters");
654 if (equalEitherForwardOrBackward(bParams, testParams))
656 GMX_ASSERT(nrfp <= MAXFORCEPARAM,
657 "This is ensured in other places, but we need this assert to keep the clang "
659 const bool identicalParameters = std::equal(bt->interactionTypes[i].forceParam().begin(),
660 bt->interactionTypes[i].forceParam().begin() + nrfp,
661 b.forceParam().begin());
663 if (!bAllowRepeat || identicalParameters)
668 if (!identicalParameters)
672 /* With dihedral type 9 we only allow for repeating
673 * of the same parameters with blocks with 1 entry.
674 * Allowing overriding is too complex to check.
676 if (!isContinuationOfBlock && !haveErrored)
679 "Encountered a second block of parameters for dihedral "
680 "type 9 for the same atoms, with either different parameters "
681 "and/or the first block has multiple lines. This is not "
686 else if (!haveWarned)
688 auto message = gmx::formatString(
689 "Bondtype %s was defined previously (e.g. in the forcefield files), "
690 "and has now been defined again. This could happen e.g. if you would "
691 "use a self-contained molecule .itp file that duplicates or replaces "
692 "the contents of the standard force-field files. You should check "
693 "the contents of your files and remove such repetition. If you know "
694 "you should override the previous definition, then you could choose "
695 "to suppress this warning with -maxwarn.%s",
696 interaction_function[ftype].longname,
697 (ftype == F_PDIHS) ? "\nUse dihedraltype 9 to allow several "
698 "multiplicity terms. Only consecutive "
699 "lines are combined. Non-consective lines "
700 "overwrite each other."
702 warning(wi, message);
704 fprintf(stderr, " old: ");
705 gmx::ArrayRef<const real> forceParam = bt->interactionTypes[i].forceParam();
706 for (int j = 0; j < nrfp; j++)
708 fprintf(stderr, " %g", forceParam[j]);
710 fprintf(stderr, " \n new: %s\n\n", line);
716 if (!identicalParameters && !bAllowRepeat)
718 /* Overwrite the parameters with the latest ones */
719 // TODO considering improving the following code by replacing with:
720 // std::copy(b->c, b->c + nrfp, bt->param[i].c);
721 gmx::ArrayRef<const real> forceParam = b.forceParam();
722 for (int j = 0; j < nrfp; j++)
724 bt->interactionTypes[i].setForceParameter(j, forceParam[j]);
732 /* fill the arrays up and down */
733 bt->interactionTypes.emplace_back(b.atoms(), b.forceParam(), b.interactionTypeName());
734 /* need to store force values because they might change below */
735 std::vector<real> forceParam(b.forceParam().begin(), b.forceParam().end());
737 /* The definitions of linear angles depend on the order of atoms,
738 * that means that for atoms i-j-k, with certain parameter a, the
739 * corresponding k-j-i angle will have parameter 1-a.
741 if (ftype == F_LINEAR_ANGLES)
743 forceParam[0] = 1 - forceParam[0];
744 forceParam[2] = 1 - forceParam[2];
746 std::vector<int> atoms;
747 gmx::ArrayRef<const int> oldAtoms = b.atoms();
748 for (auto oldAtom = oldAtoms.rbegin(); oldAtom != oldAtoms.rend(); oldAtom++)
750 atoms.emplace_back(*oldAtom);
752 bt->interactionTypes.emplace_back(atoms, forceParam, b.interactionTypeName());
756 static std::vector<int> atomTypesFromAtomNames(const PreprocessingAtomTypes* atomTypes,
757 const PreprocessingBondAtomType* bondAtomTypes,
758 gmx::ArrayRef<const char[20]> atomNames,
762 GMX_RELEASE_ASSERT(!(!atomNames.empty() && !atomTypes && !bondAtomTypes),
763 "Need to have either valid atomtypes or bondatomtypes object");
765 std::vector<int> atomTypesFromAtomNames;
766 for (const auto& name : atomNames)
768 if (atomTypes != nullptr)
770 auto atomType = atomTypes->atomTypeFromName(name);
771 if (!atomType.has_value())
773 auto message = gmx::formatString("Unknown atomtype %s\n", name);
774 warning_error_and_exit(wi, message, FARGS);
776 atomTypesFromAtomNames.emplace_back(*atomType);
778 else if (bondAtomTypes != nullptr)
780 auto bondAtomType = bondAtomTypes->bondAtomTypeFromName(name);
781 if (!bondAtomType.has_value())
783 auto message = gmx::formatString("Unknown bond_atomtype %s\n", name);
784 warning_error_and_exit(wi, message, FARGS);
786 atomTypesFromAtomNames.emplace_back(*bondAtomType);
789 return atomTypesFromAtomNames;
793 void push_bt(Directive d,
794 gmx::ArrayRef<InteractionsOfType> bt,
796 PreprocessingAtomTypes* at,
797 PreprocessingBondAtomType* bondAtomType,
801 const char* formal[MAXATOMLIST + 1] = {
802 "%s", "%s%s", "%s%s%s", "%s%s%s%s", "%s%s%s%s%s", "%s%s%s%s%s%s", "%s%s%s%s%s%s%s"
804 const char* formnl[MAXATOMLIST + 1] = { "%*s",
809 "%*s%*s%*s%*s%*s%*s",
810 "%*s%*s%*s%*s%*s%*s%*s" };
811 const char* formlf = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
812 int i, ft, ftype, nn, nrfp, nrfpA;
814 char alc[MAXATOMLIST + 1][20];
815 /* One force parameter more, so we can check if we read too many */
816 double c[MAXFORCEPARAM + 1];
818 if ((bondAtomType && at) || (!bondAtomType && !at))
820 gmx_incons("You should pass either bondAtomType or at to push_bt");
823 /* Make format string (nral ints+functype) */
824 if ((nn = sscanf(line, formal[nral], alc[0], alc[1], alc[2], alc[3], alc[4], alc[5])) != nral + 1)
826 auto message = gmx::formatString("Not enough atomtypes (%d instead of %d)", nn - 1, nral);
827 warning_error(wi, message);
831 ft = strtol(alc[nral], nullptr, 10);
832 ftype = ifunc_index(d, ft);
834 nrfpA = interaction_function[ftype].nrfpA;
835 strcpy(f1, formnl[nral]);
838 line, f1, &c[0], &c[1], &c[2], &c[3], &c[4], &c[5], &c[6], &c[7], &c[8], &c[9], &c[10], &c[11], &c[12]))
843 /* Copy the B-state from the A-state */
844 copy_B_from_A(ftype, c);
850 warning_error(wi, "Not enough parameters");
852 else if (nn > nrfpA && nn < nrfp)
854 warning_error(wi, "Too many parameters or not enough parameters for topology B");
858 warning_error(wi, "Too many parameters");
860 for (i = nn; (i < nrfp); i++)
866 std::vector<int> atomTypes =
867 atomTypesFromAtomNames(at, bondAtomType, gmx::arrayRefFromArray(alc, nral), wi);
868 std::array<real, MAXFORCEPARAM> forceParam;
869 for (int i = 0; (i < nrfp); i++)
871 forceParam[i] = c[i];
873 push_bondtype(&(bt[ftype]), InteractionOfType(atomTypes, forceParam), nral, ftype, FALSE, line, wi);
877 void push_dihedraltype(Directive d,
878 gmx::ArrayRef<InteractionsOfType> bt,
879 PreprocessingBondAtomType* bondAtomType,
883 const char* formal[MAXATOMLIST + 1] = {
884 "%s", "%s%s", "%s%s%s", "%s%s%s%s", "%s%s%s%s%s", "%s%s%s%s%s%s", "%s%s%s%s%s%s%s"
886 const char* formnl[MAXATOMLIST + 1] = { "%*s",
891 "%*s%*s%*s%*s%*s%*s",
892 "%*s%*s%*s%*s%*s%*s%*s" };
893 const char* formlf[MAXFORCEPARAM] = {
899 "%lf%lf%lf%lf%lf%lf",
900 "%lf%lf%lf%lf%lf%lf%lf",
901 "%lf%lf%lf%lf%lf%lf%lf%lf",
902 "%lf%lf%lf%lf%lf%lf%lf%lf%lf",
903 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
904 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
905 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
907 int i, ft, ftype, nn, nrfp, nrfpA, nral;
909 char alc[MAXATOMLIST + 1][20];
910 double c[MAXFORCEPARAM];
913 /* This routine accepts dihedraltypes defined from either 2 or 4 atoms.
915 * We first check for 2 atoms with the 3th column being an integer
916 * defining the type. If this isn't the case, we try it with 4 atoms
917 * and the 5th column defining the dihedral type.
919 nn = sscanf(line, formal[4], alc[0], alc[1], alc[2], alc[3], alc[4]);
920 if (nn >= 3 && strlen(alc[2]) == 1 && isdigit(alc[2][0]))
923 ft = strtol(alc[nral], nullptr, 10);
924 /* Move atom types around a bit and use 'X' for wildcard atoms
925 * to create a 4-atom dihedral definition with arbitrary atoms in
928 if (alc[2][0] == '2')
930 /* improper - the two atomtypes are 1,4. Use wildcards for 2,3 */
931 strcpy(alc[3], alc[1]);
932 sprintf(alc[2], "X");
933 sprintf(alc[1], "X");
934 /* alc[0] stays put */
938 /* proper - the two atomtypes are 2,3. Use wildcards for 1,4 */
939 sprintf(alc[3], "X");
940 strcpy(alc[2], alc[1]);
941 strcpy(alc[1], alc[0]);
942 sprintf(alc[0], "X");
945 else if (nn == 5 && strlen(alc[4]) == 1 && isdigit(alc[4][0]))
948 ft = strtol(alc[nral], nullptr, 10);
952 auto message = gmx::formatString(
953 "Incorrect number of atomtypes for dihedral (%d instead of 2 or 4)", nn - 1);
954 warning_error(wi, message);
960 /* Previously, we have always overwritten parameters if e.g. a torsion
961 with the same atomtypes occurs on multiple lines. However, CHARMM and
962 some other force fields specify multiple dihedrals over some bonds,
963 including cosines with multiplicity 6 and somethimes even higher.
964 Thus, they cannot be represented with Ryckaert-Bellemans terms.
965 To add support for these force fields, Dihedral type 9 is identical to
966 normal proper dihedrals, but repeated entries are allowed.
973 bAllowRepeat = FALSE;
977 ftype = ifunc_index(d, ft);
979 nrfpA = interaction_function[ftype].nrfpA;
981 strcpy(f1, formnl[nral]);
982 strcat(f1, formlf[nrfp - 1]);
984 /* Check number of parameters given */
986 line, f1, &c[0], &c[1], &c[2], &c[3], &c[4], &c[5], &c[6], &c[7], &c[8], &c[9], &c[10], &c[11]))
991 /* Copy the B-state from the A-state */
992 copy_B_from_A(ftype, c);
998 warning_error(wi, "Not enough parameters");
1000 else if (nn > nrfpA && nn < nrfp)
1002 warning_error(wi, "Too many parameters or not enough parameters for topology B");
1006 warning_error(wi, "Too many parameters");
1008 for (i = nn; (i < nrfp); i++)
1015 std::vector<int> atoms;
1016 std::array<real, MAXFORCEPARAM> forceParam;
1017 for (int i = 0; (i < 4); i++)
1019 if (!strcmp(alc[i], "X"))
1021 atoms.emplace_back(-1);
1025 auto atomNumber = bondAtomType->bondAtomTypeFromName(alc[i]);
1026 if (!atomNumber.has_value())
1028 auto message = gmx::formatString("Unknown bond_atomtype %s", alc[i]);
1029 warning_error_and_exit(wi, message, FARGS);
1031 atoms.emplace_back(*atomNumber);
1034 for (int i = 0; (i < nrfp); i++)
1036 forceParam[i] = c[i];
1038 /* Always use 4 atoms here, since we created two wildcard atoms
1039 * if there wasn't of them 4 already.
1041 push_bondtype(&(bt[ftype]), InteractionOfType(atoms, forceParam), 4, ftype, bAllowRepeat, line, wi);
1045 void push_nbt(Directive d, t_nbparam** nbt, PreprocessingAtomTypes* atypes, char* pline, int nb_funct, warninp* wi)
1047 /* swap the atoms */
1048 const char* form3 = "%*s%*s%*s%lf%lf%lf";
1049 const char* form4 = "%*s%*s%*s%lf%lf%lf%lf";
1050 const char* form5 = "%*s%*s%*s%lf%lf%lf%lf%lf";
1051 char a0[80], a1[80];
1052 int i, f, n, ftype, nrfp;
1058 if (sscanf(pline, "%s%s%d", a0, a1, &f) != 3)
1064 ftype = ifunc_index(d, f);
1066 if (ftype != nb_funct)
1068 auto message = gmx::formatString("Trying to add %s while the default nonbond type is %s",
1069 interaction_function[ftype].longname,
1070 interaction_function[nb_funct].longname);
1071 warning_error(wi, message);
1075 /* Get the force parameters */
1077 if (ftype == F_LJ14)
1079 n = sscanf(pline, form4, &c[0], &c[1], &c[2], &c[3]);
1085 /* When the B topology parameters are not set,
1086 * copy them from topology A
1088 GMX_ASSERT(nrfp <= 4, "LJ-14 cannot have more than 4 parameters");
1089 for (i = n; i < nrfp; i++)
1094 else if (ftype == F_LJC14_Q)
1096 n = sscanf(pline, form5, &c[0], &c[1], &c[2], &c[3], &dum);
1099 incorrect_n_param(wi);
1105 if (sscanf(pline, form3, &c[0], &c[1], &dum) != 2)
1107 incorrect_n_param(wi);
1113 if (sscanf(pline, form4, &c[0], &c[1], &c[2], &dum) != 3)
1115 incorrect_n_param(wi);
1122 gmx::formatString("Number of force parameters for nonbonded interactions is %d", nrfp);
1123 warning_error_and_exit(wi, message, FARGS);
1125 for (i = 0; (i < nrfp); i++)
1130 /* Put the parameters in the matrix */
1131 auto ai = atypes->atomTypeFromName(a0);
1132 if (!ai.has_value())
1134 auto message = gmx::formatString("Atomtype %s not found", a0);
1135 warning_error_and_exit(wi, message, FARGS);
1137 auto aj = atypes->atomTypeFromName(a1);
1138 if (!aj.has_value())
1140 auto message = gmx::formatString("Atomtype %s not found", a1);
1141 warning_error_and_exit(wi, message, FARGS);
1143 nbp = &(nbt[std::max(*ai, *aj)][std::min(*ai, *aj)]);
1148 for (i = 0; i < nrfp; i++)
1150 bId = bId && (nbp->c[i] == cr[i]);
1154 auto message = gmx::formatString(
1155 "Non-bonded parameters were defined previously (e.g. in the forcefield files), "
1156 "and have now been defined again. This could happen e.g. if you would "
1157 "use a self-contained molecule .itp file that duplicates or replaces "
1158 "the contents of the standard force-field files. You should check "
1159 "the contents of your files and remove such repetition. If you know "
1160 "you should override the previous definitions, then you could choose "
1161 "to suppress this warning with -maxwarn.");
1162 warning(wi, message);
1163 fprintf(stderr, " old:");
1164 for (i = 0; i < nrfp; i++)
1166 fprintf(stderr, " %g", nbp->c[i]);
1168 fprintf(stderr, " new\n%s\n", pline);
1172 for (i = 0; i < nrfp; i++)
1178 void push_cmaptype(Directive d,
1179 gmx::ArrayRef<InteractionsOfType> bt,
1181 PreprocessingAtomTypes* atomtypes,
1182 PreprocessingBondAtomType* bondAtomType,
1186 const char* formal = "%s%s%s%s%s%s%s%s%n";
1188 int ft, ftype, nn, nrfp, nrfpA, nrfpB;
1189 int start, nchar_consumed;
1190 int nxcmap, nycmap, ncmap, read_cmap, sl, nct;
1191 char s[20], alc[MAXATOMLIST + 2][20];
1193 /* Keep the compiler happy */
1197 GMX_ASSERT(nral == 5, "CMAP requires 5 atoms per interaction");
1199 /* Here we can only check for < 8 */
1200 if ((nn = sscanf(line, formal, alc[0], alc[1], alc[2], alc[3], alc[4], alc[5], alc[6], alc[7], &nchar_consumed))
1204 gmx::formatString("Incorrect number of atomtypes for cmap (%d instead of 5)", nn - 1);
1205 warning_error(wi, message);
1208 start += nchar_consumed;
1210 ft = strtol(alc[nral], nullptr, 10);
1211 nxcmap = strtol(alc[nral + 1], nullptr, 10);
1212 nycmap = strtol(alc[nral + 2], nullptr, 10);
1214 /* Check for equal grid spacing in x and y dims */
1215 if (nxcmap != nycmap)
1217 auto message = gmx::formatString(
1218 "Not the same grid spacing in x and y for cmap grid: x=%d, y=%d", nxcmap, nycmap);
1219 warning_error(wi, message);
1222 ncmap = nxcmap * nycmap;
1223 ftype = ifunc_index(d, ft);
1224 nrfpA = strtol(alc[6], nullptr, 10) * strtol(alc[6], nullptr, 10);
1225 nrfpB = strtol(alc[7], nullptr, 10) * strtol(alc[7], nullptr, 10);
1226 nrfp = nrfpA + nrfpB;
1228 /* Read in CMAP parameters */
1230 for (int i = 0; i < ncmap; i++)
1232 while (isspace(*(line + start + sl)))
1236 nn = sscanf(line + start + sl, " %s ", s);
1238 bt[F_CMAP].cmap.emplace_back(strtod(s, nullptr));
1247 gmx::formatString("Error in reading cmap parameter for angle %s %s %s %s %s",
1253 warning_error(wi, message);
1257 /* Check do that we got the number of parameters we expected */
1258 if (read_cmap == nrfpA)
1260 for (int i = 0; i < ncmap; i++)
1262 bt[F_CMAP].cmap.emplace_back(bt[F_CMAP].cmap[i]);
1267 if (read_cmap < nrfpA)
1269 warning_error(wi, "Not enough cmap parameters");
1271 else if (read_cmap > nrfpA && read_cmap < nrfp)
1273 warning_error(wi, "Too many cmap parameters or not enough parameters for topology B");
1275 else if (read_cmap > nrfp)
1277 warning_error(wi, "Too many cmap parameters");
1282 /* Set grid spacing and the number of grids (we assume these numbers to be the same for all
1283 * grids so we can safely assign them each time
1285 bt[F_CMAP].cmakeGridSpacing = nxcmap; /* Or nycmap, they need to be equal */
1286 bt[F_CMAP].cmapAngles++; /* Since we are incrementing here, we need to subtract later, see (*****) */
1287 nct = (nral + 1) * bt[F_CMAP].cmapAngles;
1289 for (int i = 0; (i < nral); i++)
1291 /* Assign a grid number to each cmap_type */
1292 GMX_RELEASE_ASSERT(bondAtomType != nullptr, "Need valid PreprocessingBondAtomType object");
1293 bt[F_CMAP].cmapAtomTypes.emplace_back(*bondAtomType->bondAtomTypeFromName(alc[i]));
1296 /* Assign a type number to this cmap */
1297 bt[F_CMAP].cmapAtomTypes.emplace_back(
1298 bt[F_CMAP].cmapAngles
1299 - 1); /* Since we inremented earlier, we need to subtrac here, to get the types right (****) */
1301 /* Check for the correct number of atoms (again) */
1302 if (bt[F_CMAP].nct() != nct)
1304 auto message = gmx::formatString(
1305 "Incorrect number of atom types (%d) in cmap type %d\n", nct, bt[F_CMAP].cmapAngles);
1306 warning_error(wi, message);
1308 std::vector<int> atomTypes =
1309 atomTypesFromAtomNames(atomtypes, bondAtomType, gmx::constArrayRefFromArray(alc, nral), wi);
1310 std::array<real, MAXFORCEPARAM> forceParam = { NOTSET };
1312 /* Push the bond to the bondlist */
1313 push_bondtype(&(bt[ftype]), InteractionOfType(atomTypes, forceParam), nral, ftype, FALSE, line, wi);
1317 static void push_atom_now(t_symtab* symtab,
1335 int j, resind = 0, resnr;
1339 if (((nr == 0) && (atomnr != 1)) || (nr && (atomnr != at->nr + 1)))
1341 auto message = gmx::formatString(
1342 "Atoms in the .top are not numbered consecutively from 1 (rather, "
1343 "atomnr = %d, while at->nr = %d)",
1346 warning_error_and_exit(wi, message, FARGS);
1349 j = strlen(resnumberic) - 1;
1350 if (isdigit(resnumberic[j]))
1356 ric = resnumberic[j];
1357 if (j == 0 || !isdigit(resnumberic[j - 1]))
1360 gmx::formatString("Invalid residue number '%s' for atom %d", resnumberic, atomnr);
1361 warning_error_and_exit(wi, message, FARGS);
1364 resnr = strtol(resnumberic, nullptr, 10);
1368 resind = at->atom[nr - 1].resind;
1370 if (nr == 0 || strcmp(resname, *at->resinfo[resind].name) != 0
1371 || resnr != at->resinfo[resind].nr || ric != at->resinfo[resind].ic)
1381 at->nres = resind + 1;
1382 srenew(at->resinfo, at->nres);
1383 at->resinfo[resind].name = put_symtab(symtab, resname);
1384 at->resinfo[resind].nr = resnr;
1385 at->resinfo[resind].ic = ric;
1389 resind = at->atom[at->nr - 1].resind;
1392 /* New atom instance
1393 * get new space for arrays
1395 srenew(at->atom, nr + 1);
1396 srenew(at->atomname, nr + 1);
1397 srenew(at->atomtype, nr + 1);
1398 srenew(at->atomtypeB, nr + 1);
1401 at->atom[nr].type = type;
1402 at->atom[nr].ptype = ptype;
1403 at->atom[nr].q = q0;
1404 at->atom[nr].m = m0;
1405 at->atom[nr].typeB = typeB;
1406 at->atom[nr].qB = qB;
1407 at->atom[nr].mB = mB;
1409 at->atom[nr].resind = resind;
1410 at->atom[nr].atomnumber = atomicnumber;
1411 at->atomname[nr] = put_symtab(symtab, name);
1412 at->atomtype[nr] = put_symtab(symtab, ctype);
1413 at->atomtypeB[nr] = put_symtab(symtab, ctypeB);
1417 void push_atom(t_symtab* symtab, t_atoms* at, PreprocessingAtomTypes* atypes, char* line, warninp* wi)
1419 int cgnumber, atomnr, nscan;
1420 char id[STRLEN], ctype[STRLEN], ctypeB[STRLEN], resnumberic[STRLEN], resname[STRLEN],
1421 name[STRLEN], check[STRLEN];
1422 double m, q, mb, qb;
1423 real m0, q0, mB, qB;
1425 /* Fixed parameters */
1426 if (sscanf(line, "%s%s%s%s%s%d", id, ctype, resnumberic, resname, name, &cgnumber) != 6)
1431 sscanf(id, "%d", &atomnr);
1432 auto type = atypes->atomTypeFromName(ctype);
1433 if (!type.has_value())
1435 auto message = gmx::formatString("Atomtype %s not found", ctype);
1436 warning_error_and_exit(wi, message, FARGS);
1438 ParticleType ptype = *atypes->atomParticleTypeFromAtomType(*type);
1440 /* Set default from type */
1441 q0 = *atypes->atomChargeFromAtomType(*type);
1442 m0 = *atypes->atomMassFromAtomType(*type);
1447 /* Optional parameters */
1448 nscan = sscanf(line, "%*s%*s%*s%*s%*s%*s%lf%lf%s%lf%lf%s", &q, &m, ctypeB, &qb, &mb, check);
1450 /* Nasty switch that falls thru all the way down! */
1459 typeB = atypes->atomTypeFromName(ctypeB);
1460 if (!typeB.has_value())
1462 auto message = gmx::formatString("Atomtype %s not found", ctypeB);
1463 warning_error_and_exit(wi, message, FARGS);
1465 qB = *atypes->atomChargeFromAtomType(*typeB);
1466 mB = *atypes->atomMassFromAtomType(*typeB);
1475 warning_error(wi, "Too many parameters");
1483 push_atom_now(symtab,
1486 *atypes->atomNumberFromAtomType(*type),
1496 typeB == type ? ctype : ctypeB,
1502 void push_molt(t_symtab* symtab, std::vector<MoleculeInformation>* mol, char* line, warninp* wi)
1507 if ((sscanf(line, "%s%d", type, &nrexcl)) != 2)
1509 warning_error(wi, "Expected a molecule type name and nrexcl");
1512 /* Test if this moleculetype overwrites another */
1513 const auto found = std::find_if(
1514 mol->begin(), mol->end(), [&type](const auto& m) { return strcmp(*(m.name), type) == 0; });
1515 if (found != mol->end())
1517 auto message = gmx::formatString("moleculetype %s is redefined", type);
1518 warning_error_and_exit(wi, message, FARGS);
1521 mol->emplace_back();
1522 mol->back().initMolInfo();
1524 /* Fill in the values */
1525 mol->back().name = put_symtab(symtab, type);
1526 mol->back().nrexcl = nrexcl;
1527 mol->back().excl_set = false;
1530 static bool findIfAllNBAtomsMatch(gmx::ArrayRef<const int> atomsFromParameterArray,
1531 gmx::ArrayRef<const int> atomsFromCurrentParameter,
1535 if (atomsFromParameterArray.size() != atomsFromCurrentParameter.size())
1541 for (gmx::index i = 0; i < atomsFromCurrentParameter.ssize(); i++)
1543 if (at->atom[atomsFromCurrentParameter[i]].typeB != atomsFromParameterArray[i])
1552 for (gmx::index i = 0; i < atomsFromCurrentParameter.ssize(); i++)
1554 if (at->atom[atomsFromCurrentParameter[i]].type != atomsFromParameterArray[i])
1563 static bool default_nb_params(int ftype,
1564 gmx::ArrayRef<InteractionsOfType> bt,
1566 InteractionOfType* p,
1573 InteractionOfType* pi = nullptr;
1574 int nr = bt[ftype].size();
1575 int nral = NRAL(ftype);
1576 int nrfp = interaction_function[ftype].nrfpA;
1577 int nrfpB = interaction_function[ftype].nrfpB;
1579 if ((!bB && nrfp == 0) || (bB && nrfpB == 0))
1587 /* First test the generated-pair position to save
1588 * time when we have 1000*1000 entries for e.g. OPLS...
1590 ntype = static_cast<int>(std::sqrt(static_cast<double>(nr)));
1591 GMX_ASSERT(ntype * ntype == nr,
1592 "Number of pairs of generated non-bonded parameters should be a perfect square");
1595 ti = at->atom[p->ai()].typeB;
1596 tj = at->atom[p->aj()].typeB;
1600 ti = at->atom[p->ai()].type;
1601 tj = at->atom[p->aj()].type;
1603 pi = &(bt[ftype].interactionTypes[ntype * ti + tj]);
1604 if (pi->atoms().ssize() < nral)
1606 /* not initialized yet with atom names */
1611 bFound = ((ti == pi->ai()) && (tj == pi->aj()));
1615 gmx::ArrayRef<const int> paramAtoms = p->atoms();
1616 /* Search explicitly if we didnt find it */
1619 auto foundParameter =
1620 std::find_if(bt[ftype].interactionTypes.begin(),
1621 bt[ftype].interactionTypes.end(),
1622 [¶mAtoms, &at, &bB](const auto& param) {
1623 return findIfAllNBAtomsMatch(param.atoms(), paramAtoms, at, bB);
1625 if (foundParameter != bt[ftype].interactionTypes.end())
1628 pi = &(*foundParameter);
1634 gmx::ArrayRef<const real> forceParam = pi->forceParam();
1637 if (nrfp + nrfpB > MAXFORCEPARAM)
1639 gmx_incons("Too many force parameters");
1641 for (int j = c_start; j < nrfpB; j++)
1643 p->setForceParameter(nrfp + j, forceParam[j]);
1648 for (int j = c_start; j < nrfp; j++)
1650 p->setForceParameter(j, forceParam[j]);
1656 for (int j = c_start; j < nrfp; j++)
1658 p->setForceParameter(j, 0.0);
1664 static bool default_cmap_params(gmx::ArrayRef<InteractionsOfType> bondtype,
1666 PreprocessingAtomTypes* atypes,
1667 InteractionOfType* p,
1675 bool bFound = false;
1680 /* Match the current cmap angle against the list of cmap_types */
1681 for (int i = 0; i < bondtype[F_CMAP].nct() && !bFound; i += 6)
1686 if ((atypes->bondAtomTypeFromAtomType(at->atom[p->ai()].type)
1687 == bondtype[F_CMAP].cmapAtomTypes[i])
1688 && (atypes->bondAtomTypeFromAtomType(at->atom[p->aj()].type)
1689 == bondtype[F_CMAP].cmapAtomTypes[i + 1])
1690 && (atypes->bondAtomTypeFromAtomType(at->atom[p->ak()].type)
1691 == bondtype[F_CMAP].cmapAtomTypes[i + 2])
1692 && (atypes->bondAtomTypeFromAtomType(at->atom[p->al()].type)
1693 == bondtype[F_CMAP].cmapAtomTypes[i + 3])
1694 && (atypes->bondAtomTypeFromAtomType(at->atom[p->am()].type)
1695 == bondtype[F_CMAP].cmapAtomTypes[i + 4]))
1697 /* Found cmap torsion */
1699 ct = bondtype[F_CMAP].cmapAtomTypes[i + 5];
1705 /* If we did not find a matching type for this cmap torsion */
1708 auto message = gmx::formatString("Unknown cmap torsion between atoms %d %d %d %d %d",
1714 warning_error_and_exit(wi, message, FARGS);
1717 *nparam_def = nparam_found;
1723 /* Returns the number of exact atom type matches, i.e. non wild-card matches,
1724 * returns -1 when there are no matches at all.
1726 static int findNumberOfDihedralAtomMatches(const InteractionOfType& bondType,
1727 const gmx::ArrayRef<const int> atomTypes)
1729 GMX_RELEASE_ASSERT(atomTypes.size() == 4, "Dihedrals have 4 atom types");
1730 const gmx::ArrayRef<const int> bondTypeAtomTypes = bondType.atoms();
1731 GMX_RELEASE_ASSERT(bondTypeAtomTypes.size() == 4, "Dihedral types have 4 atom types");
1732 int numExactMatches = 0;
1733 if (std::equal(bondTypeAtomTypes.begin(),
1734 bondTypeAtomTypes.end(),
1736 [&numExactMatches](int bondTypeAtomType, int atomType) {
1737 if (bondTypeAtomType == atomType)
1739 // Found an exact atom type match
1743 else if (bondTypeAtomType == -1)
1745 // Found a wildcard atom type match
1748 // Atom types do not match
1752 return numExactMatches;
1757 static std::vector<InteractionOfType>::iterator
1758 defaultInteractionsOfType(int ftype,
1759 gmx::ArrayRef<InteractionsOfType> bondType,
1760 const gmx::ArrayRef<const int> atomTypes,
1763 int nparam_found = 0;
1765 if (ftype == F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS || ftype == F_PIDIHS)
1767 int nmatch_max = -1;
1769 /* For dihedrals we allow wildcards. We choose the first type
1770 * that has the most exact matches, i.e. non-wildcard matches.
1772 auto prevPos = bondType[ftype].interactionTypes.end();
1773 auto pos = bondType[ftype].interactionTypes.begin();
1774 while (pos != bondType[ftype].interactionTypes.end() && nmatch_max < 4)
1776 pos = std::find_if(bondType[ftype].interactionTypes.begin(),
1777 bondType[ftype].interactionTypes.end(),
1778 [&atomTypes, &nmatch_max](const auto& bondType) {
1779 return (findNumberOfDihedralAtomMatches(bondType, atomTypes) > nmatch_max);
1781 if (pos != bondType[ftype].interactionTypes.end())
1784 nmatch_max = findNumberOfDihedralAtomMatches(*pos, atomTypes);
1788 if (prevPos != bondType[ftype].interactionTypes.end())
1792 /* Find additional matches for this dihedral - necessary
1794 * The rule in that case is that additional matches
1795 * HAVE to be on adjacent lines!
1798 // Advance iterator (like std::advance) without incrementing past end (UB)
1799 const auto safeAdvance = [](auto& it, auto n, auto end) {
1800 it = end - it > n ? it + n : end;
1802 /* Continue from current iterator position */
1803 auto nextPos = prevPos;
1804 const auto endIter = bondType[ftype].interactionTypes.end();
1805 safeAdvance(nextPos, 2, endIter);
1806 for (; nextPos < endIter && bSame; safeAdvance(nextPos, 2, endIter))
1808 bSame = (prevPos->ai() == nextPos->ai() && prevPos->aj() == nextPos->aj()
1809 && prevPos->ak() == nextPos->ak() && prevPos->al() == nextPos->al());
1814 /* nparam_found will be increased as long as the numbers match */
1817 *nparam_def = nparam_found;
1820 else /* Not a dihedral */
1822 auto found = std::find_if(
1823 bondType[ftype].interactionTypes.begin(),
1824 bondType[ftype].interactionTypes.end(),
1825 [&atomTypes](const auto& param) {
1826 return std::equal(param.atoms().begin(), param.atoms().end(), atomTypes.begin());
1828 if (found != bondType[ftype].interactionTypes.end())
1832 *nparam_def = nparam_found;
1838 void push_bond(Directive d,
1839 gmx::ArrayRef<InteractionsOfType> bondtype,
1840 gmx::ArrayRef<InteractionsOfType> bond,
1842 PreprocessingAtomTypes* atypes,
1848 bool* bWarn_copy_A_B,
1851 const char* aaformat[MAXATOMLIST] = { "%d%d", "%d%d%d", "%d%d%d%d",
1852 "%d%d%d%d%d", "%d%d%d%d%d%d", "%d%d%d%d%d%d%d" };
1853 const char* asformat[MAXATOMLIST] = {
1854 "%*s%*s", "%*s%*s%*s", "%*s%*s%*s%*s",
1855 "%*s%*s%*s%*s%*s", "%*s%*s%*s%*s%*s%*s", "%*s%*s%*s%*s%*s%*s%*s"
1857 const char* ccformat = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
1858 int nral, nral_fmt, nread, ftype;
1859 char format[STRLEN];
1860 /* One force parameter more, so we can check if we read too many */
1861 double cc[MAXFORCEPARAM + 1];
1862 std::array<int, MAXATOMLIST + 1> aa;
1863 bool bFoundA = FALSE, bFoundB = FALSE, bDef, bSwapParity = FALSE;
1864 int nparam_defA, nparam_defB;
1866 nparam_defA = nparam_defB = 0;
1868 ftype = ifunc_index(d, 1);
1870 for (int j = 0; j < nral; j++)
1874 bDef = (NRFP(ftype) > 0);
1876 if (ftype == F_SETTLE)
1878 /* SETTLE acts on 3 atoms, but the topology format only specifies
1879 * the first atom (for historical reasons).
1888 nread = sscanf(line, aaformat[nral_fmt - 1], &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
1890 if (ftype == F_SETTLE)
1897 if (nread < nral_fmt)
1902 else if (nread > nral_fmt)
1904 /* this is a hack to allow for virtual sites with swapped parity */
1905 bSwapParity = (aa[nral] < 0);
1908 aa[nral] = -aa[nral];
1910 ftype = ifunc_index(d, aa[nral]);
1916 case F_VSITE3OUT: break;
1919 gmx::formatString("Negative function types only allowed for %s and %s",
1920 interaction_function[F_VSITE3FAD].longname,
1921 interaction_function[F_VSITE3OUT].longname);
1922 warning_error_and_exit(wi, message, FARGS);
1928 /* Check for double atoms and atoms out of bounds, then convert to 0-based indexing */
1929 for (int i = 0; (i < nral); i++)
1931 if (aa[i] < 1 || aa[i] > at->nr)
1933 auto message = gmx::formatString(
1934 "Atom index (%d) in %s out of bounds (1-%d).\n"
1935 "This probably means that you have inserted topology section \"%s\"\n"
1936 "in a part belonging to a different molecule than you intended to.\n"
1937 "In that case move the \"%s\" section to the right molecule.",
1939 enumValueToString(d),
1941 enumValueToString(d),
1942 enumValueToString(d));
1943 warning_error_and_exit(wi, message, FARGS);
1945 for (int j = i + 1; (j < nral); j++)
1947 GMX_ASSERT(j < MAXATOMLIST + 1,
1948 "Values from nral=NRAL() will satisfy this, we assert to keep gcc 4 happy");
1951 auto message = gmx::formatString(
1952 "Duplicate atom index (%d) in %s", aa[i], enumValueToString(d));
1953 if (ftype == F_ANGRES)
1955 /* Since the angle restraints uses 2 pairs of atoms to
1956 * defines an angle between vectors, it can be useful
1957 * to use one atom twice, so we only issue a note here.
1959 warning_note(wi, message);
1963 warning_error(wi, message);
1968 // Convert to 0-based indexing
1972 // These are the atom indices for this interaction
1973 auto atomIndices = gmx::ArrayRef<const int>(aa).subArray(0, nral);
1975 // Look up the A-state atom types for this interaction
1976 std::vector<int> atomTypes(atomIndices.size());
1977 std::transform(atomIndices.begin(), atomIndices.end(), atomTypes.begin(), [at, atypes](const int atomIndex) {
1978 return atypes->bondAtomTypeFromAtomType(at->atom[atomIndex].type).value();
1980 // Look up the B-state atom types for this interaction
1981 std::vector<int> atomTypesB(atomIndices.size());
1982 std::transform(atomIndices.begin(), atomIndices.end(), atomTypesB.begin(), [at, atypes](const int atomIndex) {
1983 return atypes->bondAtomTypeFromAtomType(at->atom[atomIndex].typeB).value();
1986 /* default force parameters */
1987 /* need to have an empty but initialized param array for some reason */
1988 std::array<real, MAXFORCEPARAM> forceParam = { 0.0 };
1990 /* Get force params for normal and free energy perturbation
1991 * studies, as determined by types!
1993 InteractionOfType param(atomIndices, forceParam, "");
1995 std::vector<InteractionOfType>::iterator foundAParameter = bondtype[ftype].interactionTypes.end();
1996 std::vector<InteractionOfType>::iterator foundBParameter = bondtype[ftype].interactionTypes.end();
1999 if (NRFPA(ftype) == 0)
2005 foundAParameter = defaultInteractionsOfType(ftype, bondtype, atomTypes, &nparam_defA);
2006 if (foundAParameter != bondtype[ftype].interactionTypes.end())
2008 /* Copy the A-state and B-state default parameters. */
2009 GMX_ASSERT(NRFPA(ftype) + NRFPB(ftype) <= MAXFORCEPARAM,
2010 "Bonded interactions may have at most 12 parameters");
2011 gmx::ArrayRef<const real> defaultParam = foundAParameter->forceParam();
2012 for (int j = 0; (j < NRFPA(ftype) + NRFPB(ftype)); j++)
2014 param.setForceParameter(j, defaultParam[j]);
2020 if (NRFPB(ftype) == 0)
2026 foundBParameter = defaultInteractionsOfType(ftype, bondtype, atomTypesB, &nparam_defB);
2027 if (foundBParameter != bondtype[ftype].interactionTypes.end())
2029 /* Copy only the B-state default parameters */
2030 gmx::ArrayRef<const real> defaultParam = foundBParameter->forceParam();
2031 for (int j = NRFPA(ftype); (j < NRFP(ftype)); j++)
2033 param.setForceParameter(j, defaultParam[j]);
2039 else if (ftype == F_LJ14)
2041 bFoundA = default_nb_params(ftype, bondtype, at, ¶m, 0, FALSE, bGenPairs);
2042 bFoundB = default_nb_params(ftype, bondtype, at, ¶m, 0, TRUE, bGenPairs);
2044 else if (ftype == F_LJC14_Q)
2046 /* Fill in the A-state charges as default parameters */
2047 param.setForceParameter(0, fudgeQQ);
2048 param.setForceParameter(1, at->atom[param.ai()].q);
2049 param.setForceParameter(2, at->atom[param.aj()].q);
2050 /* The default LJ parameters are the standard 1-4 parameters */
2051 bFoundA = default_nb_params(F_LJ14, bondtype, at, ¶m, 3, FALSE, bGenPairs);
2054 else if (ftype == F_LJC_PAIRS_NB)
2056 /* Defaults are not supported here */
2062 gmx_incons("Unknown function type in push_bond");
2065 if (nread > nral_fmt)
2067 /* Manually specified parameters - in this case we discard multiple torsion info! */
2069 strcpy(format, asformat[nral_fmt - 1]);
2070 strcat(format, ccformat);
2072 nread = sscanf(line,
2088 if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0))
2090 /* We only have to issue a warning if these atoms are perturbed! */
2092 gmx::ArrayRef<const int> paramAtoms = param.atoms();
2093 for (int j = 0; (j < nral); j++)
2095 bPert = bPert || PERTURBED(at->atom[paramAtoms[j]]);
2098 if (bPert && *bWarn_copy_A_B)
2100 auto message = gmx::formatString(
2101 "Some parameters for bonded interaction involving "
2102 "perturbed atoms are specified explicitly in "
2103 "state A, but not B - copying A to B");
2104 warning(wi, message);
2105 *bWarn_copy_A_B = FALSE;
2108 /* If only the A parameters were specified, copy them to the B state */
2109 /* The B-state parameters correspond to the first nrfpB
2110 * A-state parameters.
2112 for (int j = 0; (j < NRFPB(ftype)); j++)
2114 cc[nread++] = cc[j];
2118 /* If nread was 0 or EOF, no parameters were read => use defaults.
2119 * If nread was nrfpA we copied above so nread=nrfp.
2120 * If nread was nrfp we are cool.
2121 * For F_LJC14_Q we allow supplying fudgeQQ only.
2122 * Anything else is an error!
2124 if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) && !(ftype == F_LJC14_Q && nread == 1))
2126 auto message = gmx::formatString(
2127 "Incorrect number of parameters - found %d, expected %d "
2128 "or %d for %s (after the function type).",
2132 interaction_function[ftype].longname);
2133 warning_error_and_exit(wi, message, FARGS);
2136 for (int j = 0; (j < nread); j++)
2138 param.setForceParameter(j, cc[j]);
2140 /* Check whether we have to use the defaults */
2141 if (nread == NRFP(ftype))
2150 /* nread now holds the number of force parameters read! */
2155 /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */
2156 if (ftype == F_PDIHS)
2158 if ((nparam_defA != nparam_defB)
2159 || ((nparam_defA > 1 || nparam_defB > 1) && (foundAParameter != foundBParameter)))
2161 auto message = gmx::formatString(
2162 "Cannot automatically perturb a torsion with multiple terms to different "
2164 "Please specify perturbed parameters manually for this torsion in your "
2166 warning_error(wi, message);
2170 if (nread > 0 && nread < NRFPA(ftype))
2172 /* Issue an error, do not use defaults */
2173 auto message = gmx::formatString(
2174 "Not enough parameters, there should be at least %d (or 0 for defaults)", NRFPA(ftype));
2175 warning_error(wi, message);
2178 if (nread == 0 || nread == EOF)
2182 if (interaction_function[ftype].flags & IF_VSITE)
2184 for (int j = 0; j < MAXFORCEPARAM; j++)
2186 param.setForceParameter(j, NOTSET);
2190 /* flag to swap parity of vsi te construction */
2191 param.setForceParameter(1, -1);
2199 "NOTE: No default %s types, using zeroes\n",
2200 interaction_function[ftype].longname);
2204 auto message = gmx::formatString("No default %s types",
2205 interaction_function[ftype].longname);
2206 warning_error(wi, message);
2216 case F_VSITE3FAD: param.setForceParameter(0, 360 - param.c0()); break;
2217 case F_VSITE3OUT: param.setForceParameter(2, -param.c2()); break;
2223 /* We only have to issue a warning if these atoms are perturbed! */
2225 gmx::ArrayRef<const int> paramAtoms = param.atoms();
2226 for (int j = 0; (j < nral); j++)
2228 bPert = bPert || PERTURBED(at->atom[paramAtoms[j]]);
2233 auto message = gmx::formatString(
2234 "No default %s types for perturbed atoms, "
2235 "using normal values",
2236 interaction_function[ftype].longname);
2237 warning(wi, message);
2243 gmx::ArrayRef<const real> paramValue = param.forceParam();
2244 if ((ftype == F_PDIHS || ftype == F_ANGRES || ftype == F_ANGRESZ) && paramValue[5] != paramValue[2])
2246 auto message = gmx::formatString("%s multiplicity can not be perturbed %f!=%f",
2247 interaction_function[ftype].longname,
2250 warning_error_and_exit(wi, message, FARGS);
2253 if (IS_TABULATED(ftype) && param.c0() != param.c2())
2255 auto message = gmx::formatString("%s table number can not be perturbed %d!=%d",
2256 interaction_function[ftype].longname,
2257 gmx::roundToInt(param.c0()),
2258 gmx::roundToInt(param.c0()));
2259 warning_error_and_exit(wi, message, FARGS);
2262 /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
2263 if (ftype == F_RBDIHS)
2267 for (int i = 0; i < NRFP(ftype); i++)
2269 if (paramValue[i] != 0.0)
2280 /* Put the values in the appropriate arrays */
2281 add_param_to_list(&bond[ftype], param);
2283 /* Push additional torsions from FF for ftype==9 if we have them.
2284 * We have already checked that the A/B states do not differ in this case,
2285 * so we do not have to double-check that again, or the vsite stuff.
2286 * In addition, those torsions cannot be automatically perturbed.
2288 if (bDef && ftype == F_PDIHS)
2290 for (int i = 1; i < nparam_defA; i++)
2292 /* Advance pointer! */
2293 foundAParameter += 2;
2294 gmx::ArrayRef<const real> forceParam = foundAParameter->forceParam();
2295 for (int j = 0; j < (NRFPA(ftype) + NRFPB(ftype)); j++)
2297 param.setForceParameter(j, forceParam[j]);
2299 /* And push the next term for this torsion */
2300 add_param_to_list(&bond[ftype], param);
2305 void push_cmap(Directive d,
2306 gmx::ArrayRef<InteractionsOfType> bondtype,
2307 gmx::ArrayRef<InteractionsOfType> bond,
2309 PreprocessingAtomTypes* atypes,
2313 const char* aaformat[MAXATOMLIST + 1] = {
2314 "%d", "%d%d", "%d%d%d", "%d%d%d%d", "%d%d%d%d%d", "%d%d%d%d%d%d", "%d%d%d%d%d%d%d"
2317 int ftype, nral, nread, ncmap_params;
2319 int aa[MAXATOMLIST + 1];
2322 ftype = ifunc_index(d, 1);
2326 nread = sscanf(line, aaformat[nral - 1], &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
2333 else if (nread == nral)
2335 ftype = ifunc_index(d, 1);
2338 /* Check for double atoms and atoms out of bounds */
2339 for (int i = 0; i < nral; i++)
2341 if (aa[i] < 1 || aa[i] > at->nr)
2343 auto message = gmx::formatString(
2344 "Atom index (%d) in %s out of bounds (1-%d).\n"
2345 "This probably means that you have inserted topology section \"%s\"\n"
2346 "in a part belonging to a different molecule than you intended to.\n"
2347 "In that case move the \"%s\" section to the right molecule.",
2349 enumValueToString(d),
2351 enumValueToString(d),
2352 enumValueToString(d));
2353 warning_error_and_exit(wi, message, FARGS);
2356 for (int j = i + 1; (j < nral); j++)
2360 auto message = gmx::formatString(
2361 "Duplicate atom index (%d) in %s", aa[i], enumValueToString(d));
2362 warning_error(wi, message);
2367 /* default force parameters */
2368 std::vector<int> atoms;
2369 for (int j = 0; (j < nral); j++)
2371 atoms.emplace_back(aa[j] - 1);
2373 std::array<real, MAXFORCEPARAM> forceParam = { 0.0 };
2374 InteractionOfType param(atoms, forceParam, "");
2375 /* Get the cmap type for this cmap angle */
2376 bFound = default_cmap_params(bondtype, at, atypes, ¶m, FALSE, &cmap_type, &ncmap_params, wi);
2378 /* We want exactly one parameter (the cmap type in state A (currently no state B) back */
2379 if (bFound && ncmap_params == 1)
2381 /* Put the values in the appropriate arrays */
2382 param.setForceParameter(0, cmap_type);
2383 add_param_to_list(&bond[ftype], param);
2387 /* This is essentially the same check as in default_cmap_params() done one more time */
2389 gmx::formatString("Unable to assign a cmap type to torsion %d %d %d %d and %d\n",
2395 warning_error_and_exit(wi, message, FARGS);
2400 void push_vsitesn(Directive d, gmx::ArrayRef<InteractionsOfType> bond, t_atoms* at, char* line, warninp* wi)
2403 int type, ftype, n, ret, nj, a;
2405 double *weight = nullptr, weight_tot;
2407 std::array<real, MAXFORCEPARAM> forceParam = { 0.0 };
2409 ret = sscanf(ptr, "%d%n", &a, &n);
2414 gmx::formatString("Expected an atom index in section \"%s\"", enumValueToString(d));
2415 warning_error_and_exit(wi, message, FARGS);
2418 sscanf(ptr, "%d%n", &type, &n);
2420 ftype = ifunc_index(d, type);
2421 int firstAtom = a - 1;
2427 ret = sscanf(ptr, "%d%n", &a, &n);
2433 srenew(atc, nj + 20);
2434 srenew(weight, nj + 20);
2439 case 1: weight[nj] = 1; break;
2441 /* Here we use the A-state mass as a parameter.
2442 * Note that the B-state mass has no influence.
2444 weight[nj] = at->atom[atc[nj]].m;
2448 ret = sscanf(ptr, "%lf%n", &(weight[nj]), &n);
2452 auto message = gmx::formatString(
2453 "No weight or negative weight found for vsiten "
2454 "constructing atom %d (atom index %d)",
2457 warning_error_and_exit(wi, message, FARGS);
2461 auto message = gmx::formatString("Unknown vsiten type %d", type);
2462 warning_error_and_exit(wi, message, FARGS);
2464 weight_tot += weight[nj];
2471 auto message = gmx::formatString("Expected more than one atom index in section \"%s\"",
2472 enumValueToString(d));
2473 warning_error_and_exit(wi, message, FARGS);
2476 if (weight_tot == 0)
2478 warning_error_and_exit(wi, "The total mass of the construting atoms is zero", FARGS);
2481 for (int j = 0; j < nj; j++)
2483 std::vector<int> atoms = { firstAtom, atc[j] };
2485 forceParam[1] = weight[j] / weight_tot;
2486 /* Put the values in the appropriate arrays */
2487 add_param_to_list(&bond[ftype], InteractionOfType(atoms, forceParam));
2494 void push_mol(gmx::ArrayRef<MoleculeInformation> mols, char* pline, int* whichmol, int* nrcopies, warninp* wi)
2498 if (sscanf(pline, "%s%d", type, nrcopies) != 2)
2504 /* Search moleculename.
2505 * Here we originally only did case insensitive matching. But because
2506 * some PDB files can have many chains and use case to generate more
2507 * chain-identifiers, which in turn end up in our moleculetype name,
2508 * we added support for case-sensitivity.
2515 for (const auto& mol : mols)
2517 if (strcmp(type, *(mol.name)) == 0)
2522 if (gmx_strcasecmp(type, *(mol.name)) == 0)
2532 // select the case sensitive match
2533 *whichmol = matchcs;
2537 // avoid matching case-insensitive when we have multiple matches
2540 auto message = gmx::formatString(
2541 "For moleculetype '%s' in [ system ] %d case insensitive "
2542 "matches, but %d case sensitive matches were found. Check "
2543 "the case of the characters in the moleculetypes.",
2547 warning_error_and_exit(wi, message, FARGS);
2551 // select the unique case insensitive match
2552 *whichmol = matchci;
2556 auto message = gmx::formatString("No such moleculetype %s", type);
2557 warning_error_and_exit(wi, message, FARGS);
2562 void push_excl(char* line, gmx::ArrayRef<gmx::ExclusionBlock> b2, warninp* wi)
2566 char base[STRLEN], format[STRLEN];
2568 if (sscanf(line, "%d", &i) == 0)
2573 if ((1 <= i) && (i <= b2.ssize()))
2581 strcpy(base, "%*d");
2584 strcpy(format, base);
2585 strcat(format, "%d");
2586 n = sscanf(line, format, &j);
2589 if ((1 <= j) && (j <= b2.ssize()))
2592 b2[i].atomNumber.push_back(j);
2593 /* also add the reverse exclusion! */
2594 b2[j].atomNumber.push_back(i);
2595 strcat(base, "%*d");
2599 auto message = gmx::formatString("Invalid Atomnr j: %d, b2->nr: %zu\n", j, b2.size());
2600 warning_error_and_exit(wi, message, FARGS);
2606 int add_atomtype_decoupled(t_symtab* symtab, PreprocessingAtomTypes* at, t_nbparam*** nbparam, t_nbparam*** pair)
2611 /* Define an atom type with all parameters set to zero (no interactions) */
2614 /* Type for decoupled atoms could be anything,
2615 * this should be changed automatically later when required.
2617 atom.ptype = ParticleType::Atom;
2619 std::array<real, MAXFORCEPARAM> forceParam = { 0.0 };
2620 nr = at->addType(symtab, atom, "decoupled", InteractionOfType({}, forceParam, ""), -1, 0);
2622 /* Add space in the non-bonded parameters matrix */
2623 realloc_nb_params(at, nbparam, pair);
2628 static void convert_pairs_to_pairsQ(gmx::ArrayRef<InteractionsOfType> interactions, real fudgeQQ, t_atoms* atoms)
2630 /* Add the pair list to the pairQ list */
2631 std::vector<InteractionOfType> paramnew;
2633 gmx::ArrayRef<const InteractionOfType> paramp1 = interactions[F_LJ14].interactionTypes;
2634 gmx::ArrayRef<const InteractionOfType> paramp2 = interactions[F_LJC14_Q].interactionTypes;
2636 /* Fill in the new F_LJC14_Q array with the old one. NOTE:
2637 it may be possible to just ADD the converted F_LJ14 array
2638 to the old F_LJC14_Q array, but since we have to create
2639 a new sized memory structure, better just to deep copy it all.
2643 for (const auto& param : paramp2)
2645 paramnew.emplace_back(param);
2648 for (const auto& param : paramp1)
2650 std::vector<real> forceParam = {
2651 fudgeQQ, atoms->atom[param.ai()].q, atoms->atom[param.aj()].q, param.c0(), param.c1()
2653 paramnew.emplace_back(param.atoms(), forceParam, "");
2656 /* now assign the new data to the F_LJC14_Q structure */
2657 interactions[F_LJC14_Q].interactionTypes = paramnew;
2659 /* Empty the LJ14 pairlist */
2660 interactions[F_LJ14].interactionTypes.clear();
2663 static void generate_LJCpairsNB(MoleculeInformation* mol, int nb_funct, InteractionsOfType* nbp, warninp* wi)
2669 atom = mol->atoms.atom;
2671 ntype = static_cast<int>(std::sqrt(static_cast<double>(nbp->size())));
2672 GMX_ASSERT(ntype * ntype == gmx::ssize(*nbp),
2673 "Number of pairs of generated non-bonded parameters should be a perfect square");
2675 /* Add a pair interaction for all non-excluded atom pairs */
2676 const auto& excls = mol->excls;
2677 for (int i = 0; i < n; i++)
2679 for (int j = i + 1; j < n; j++)
2681 bool pairIsExcluded = false;
2682 for (const int atomK : excls[i])
2686 pairIsExcluded = true;
2689 if (!pairIsExcluded)
2691 if (nb_funct != F_LJ)
2693 auto message = gmx::formatString(
2694 "Can only generate non-bonded pair interactions "
2695 "for Van der Waals type Lennard-Jones");
2696 warning_error_and_exit(wi, message, FARGS);
2698 std::vector<int> atoms = { i, j };
2699 std::vector<real> forceParam = {
2702 nbp->interactionTypes[ntype * atom[i].type + atom[j].type].c0(),
2703 nbp->interactionTypes[ntype * atom[i].type + atom[j].type].c1()
2705 add_param_to_list(&mol->interactions[F_LJC_PAIRS_NB], InteractionOfType(atoms, forceParam));
2711 static void set_excl_all(gmx::ListOfLists<int>* excl)
2713 /* Get rid of the current exclusions and exclude all atom pairs */
2714 const int numAtoms = excl->ssize();
2715 std::vector<int> exclusionsForAtom(numAtoms);
2716 for (int i = 0; i < numAtoms; i++)
2718 exclusionsForAtom[i] = i;
2721 for (int i = 0; i < numAtoms; i++)
2723 excl->pushBack(exclusionsForAtom);
2727 static void decouple_atoms(t_atoms* atoms,
2728 int atomtype_decouple,
2731 const char* mol_name,
2736 for (i = 0; i < atoms->nr; i++)
2740 atom = &atoms->atom[i];
2742 if (atom->qB != atom->q || atom->typeB != atom->type)
2744 auto message = gmx::formatString(
2745 "Atom %d in molecule type '%s' has different A and B state "
2746 "charges and/or atom types set in the topology file as well "
2747 "as through the mdp option '%s'. You can not use both "
2748 "these methods simultaneously.",
2752 warning_error_and_exit(wi, message, FARGS);
2755 if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamVDW)
2759 if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamQ)
2761 atom->type = atomtype_decouple;
2763 if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamVDW)
2767 if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamQ)
2769 atom->typeB = atomtype_decouple;
2774 void convert_moltype_couple(MoleculeInformation* mol,
2775 int atomtype_decouple,
2781 InteractionsOfType* nbp,
2784 convert_pairs_to_pairsQ(mol->interactions, fudgeQQ, &mol->atoms);
2787 generate_LJCpairsNB(mol, nb_funct, nbp, wi);
2788 set_excl_all(&mol->excls);
2790 decouple_atoms(&mol->atoms, atomtype_decouple, couple_lam0, couple_lam1, *mol->name, wi);