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,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.
49 #include "gromacs/fileio/warninp.h"
50 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
51 #include "gromacs/gmxpreprocess/gpp_bond_atomtype.h"
52 #include "gromacs/gmxpreprocess/grompp_impl.h"
53 #include "gromacs/gmxpreprocess/notset.h"
54 #include "gromacs/gmxpreprocess/readir.h"
55 #include "gromacs/gmxpreprocess/topdirs.h"
56 #include "gromacs/gmxpreprocess/toputil.h"
57 #include "gromacs/math/functions.h"
58 #include "gromacs/mdtypes/md_enums.h"
59 #include "gromacs/topology/exclusionblocks.h"
60 #include "gromacs/topology/ifunc.h"
61 #include "gromacs/topology/symtab.h"
62 #include "gromacs/utility/arrayref.h"
63 #include "gromacs/utility/cstringutil.h"
64 #include "gromacs/utility/fatalerror.h"
65 #include "gromacs/utility/gmxassert.h"
66 #include "gromacs/utility/smalloc.h"
67 #include "gromacs/utility/stringutil.h"
69 void generate_nbparams(int comb,
71 InteractionsOfType* interactions,
72 PreprocessingAtomTypes* atypes,
76 real c, bi, bj, ci, cj, ci0, ci1, ci2, cj0, cj1, cj2;
78 /* Lean mean shortcuts */
81 interactions->interactionTypes.clear();
83 std::array<real, MAXFORCEPARAM> forceParam = { NOTSET };
84 /* Fill the matrix with force parameters */
92 for (int i = 0; (i < nr); i++)
94 for (int j = 0; (j < nr); j++)
96 for (int nf = 0; (nf < nrfp); nf++)
98 ci = atypes->atomNonBondedParamFromAtomType(i, nf);
99 cj = atypes->atomNonBondedParamFromAtomType(j, nf);
100 c = std::sqrt(ci * cj);
103 interactions->interactionTypes.emplace_back(InteractionOfType({}, forceParam));
108 case eCOMB_ARITHMETIC:
109 /* c0 and c1 are sigma and epsilon */
110 for (int i = 0; (i < nr); i++)
112 for (int j = 0; (j < nr); j++)
114 ci0 = atypes->atomNonBondedParamFromAtomType(i, 0);
115 cj0 = atypes->atomNonBondedParamFromAtomType(j, 0);
116 ci1 = atypes->atomNonBondedParamFromAtomType(i, 1);
117 cj1 = atypes->atomNonBondedParamFromAtomType(j, 1);
118 forceParam[0] = (fabs(ci0) + fabs(cj0)) * 0.5;
119 /* Negative sigma signals that c6 should be set to zero later,
120 * so we need to propagate that through the combination rules.
122 if (ci0 < 0 || cj0 < 0)
126 forceParam[1] = std::sqrt(ci1 * cj1);
127 interactions->interactionTypes.emplace_back(InteractionOfType({}, forceParam));
132 case eCOMB_GEOM_SIG_EPS:
133 /* c0 and c1 are sigma and epsilon */
134 for (int i = 0; (i < nr); i++)
136 for (int j = 0; (j < nr); j++)
138 ci0 = atypes->atomNonBondedParamFromAtomType(i, 0);
139 cj0 = atypes->atomNonBondedParamFromAtomType(j, 0);
140 ci1 = atypes->atomNonBondedParamFromAtomType(i, 1);
141 cj1 = atypes->atomNonBondedParamFromAtomType(j, 1);
142 forceParam[0] = std::sqrt(std::fabs(ci0 * cj0));
143 /* Negative sigma signals that c6 should be set to zero later,
144 * so we need to propagate that through the combination rules.
146 if (ci0 < 0 || cj0 < 0)
150 forceParam[1] = std::sqrt(ci1 * cj1);
151 interactions->interactionTypes.emplace_back(InteractionOfType({}, forceParam));
157 auto message = gmx::formatString("No such combination rule %d", comb);
158 warning_error_and_exit(wi, message, FARGS);
163 /* Buckingham rules */
164 for (int i = 0; (i < nr); i++)
166 for (int j = 0; (j < nr); j++)
168 ci0 = atypes->atomNonBondedParamFromAtomType(i, 0);
169 cj0 = atypes->atomNonBondedParamFromAtomType(j, 0);
170 ci2 = atypes->atomNonBondedParamFromAtomType(i, 2);
171 cj2 = atypes->atomNonBondedParamFromAtomType(j, 2);
172 bi = atypes->atomNonBondedParamFromAtomType(i, 1);
173 bj = atypes->atomNonBondedParamFromAtomType(j, 1);
174 forceParam[0] = std::sqrt(ci0 * cj0);
175 if ((bi == 0) || (bj == 0))
181 forceParam[1] = 2.0 / (1 / bi + 1 / bj);
183 forceParam[2] = std::sqrt(ci2 * cj2);
184 interactions->interactionTypes.emplace_back(InteractionOfType({}, forceParam));
190 auto message = gmx::formatString("Invalid nonbonded type %s",
191 interaction_function[ftype].longname);
192 warning_error(wi, message);
196 /*! \brief Used to temporarily store the explicit non-bonded parameter
197 * combinations, which will be copied to InteractionsOfType. */
200 //! Has this combination been set.
202 //! The non-bonded parameters
206 static void realloc_nb_params(PreprocessingAtomTypes* atypes, t_nbparam*** nbparam, t_nbparam*** pair)
208 /* Add space in the non-bonded parameters matrix */
209 int atnr = atypes->size();
210 srenew(*nbparam, atnr);
211 snew((*nbparam)[atnr - 1], atnr);
215 snew((*pair)[atnr - 1], atnr);
219 int copy_nbparams(t_nbparam** param, int ftype, InteractionsOfType* interactions, int nr)
226 for (int i = 0; i < nr; i++)
228 for (int j = 0; j <= i; j++)
230 GMX_RELEASE_ASSERT(param, "Must have valid parameters");
231 if (param[i][j].bSet)
233 for (int f = 0; f < nrfp; f++)
235 interactions->interactionTypes[nr * i + j].setForceParameter(f, param[i][j].c[f]);
236 interactions->interactionTypes[nr * j + i].setForceParameter(f, param[i][j].c[f]);
246 void free_nbparam(t_nbparam** param, int nr)
250 GMX_RELEASE_ASSERT(param, "Must have valid parameters");
251 for (i = 0; i < nr; i++)
253 GMX_RELEASE_ASSERT(param[i], "Must have valid parameters");
259 static void copy_B_from_A(int ftype, double* c)
263 nrfpA = NRFPA(ftype);
264 nrfpB = NRFPB(ftype);
266 /* Copy the B parameters from the first nrfpB A parameters */
267 for (i = 0; (i < nrfpB); i++)
273 void push_at(t_symtab* symtab,
274 PreprocessingAtomTypes* at,
275 PreprocessingBondAtomType* bondAtomType,
278 t_nbparam*** nbparam,
287 t_xlate xl[eptNR] = {
288 { "A", eptAtom }, { "N", eptNucleus }, { "S", eptShell },
289 { "B", eptBond }, { "V", eptVSite },
292 int nr, nfields, j, pt, nfp0 = -1;
293 int batype_nr, nread;
294 char type[STRLEN], btype[STRLEN], ptype[STRLEN];
296 double c[MAXFORCEPARAM];
297 char tmpfield[12][100]; /* Max 12 fields of width 100 */
300 bool have_atomic_number;
301 bool have_bonded_type;
305 /* First assign input line to temporary array */
306 nfields = sscanf(line, "%s%s%s%s%s%s%s%s%s%s%s%s", tmpfield[0], tmpfield[1], tmpfield[2],
307 tmpfield[3], tmpfield[4], tmpfield[5], tmpfield[6], tmpfield[7], tmpfield[8],
308 tmpfield[9], tmpfield[10], tmpfield[11]);
310 /* Comments on optional fields in the atomtypes section:
312 * The force field format is getting a bit old. For OPLS-AA we needed
313 * to add a special bonded atomtype, and for Gerrit Groenhofs QM/MM stuff
314 * we also needed the atomic numbers.
315 * To avoid making all old or user-generated force fields unusable we
316 * have introduced both these quantities as optional columns, and do some
317 * acrobatics to check whether they are present or not.
318 * This will all look much nicer when we switch to XML... sigh.
320 * Field 0 (mandatory) is the nonbonded type name. (string)
321 * Field 1 (optional) is the bonded type (string)
322 * Field 2 (optional) is the atomic number (int)
323 * Field 3 (mandatory) is the mass (numerical)
324 * Field 4 (mandatory) is the charge (numerical)
325 * Field 5 (mandatory) is the particle type (single character)
326 * This is followed by a number of nonbonded parameters.
328 * The safest way to identify the format is the particle type field.
330 * So, here is what we do:
332 * A. Read in the first six fields as strings
333 * B. If field 3 (starting from 0) is a single char, we have neither
334 * bonded_type or atomic numbers.
335 * C. If field 5 is a single char we have both.
336 * D. If field 4 is a single char we check field 1. If this begins with
337 * an alphabetical character we have bonded types, otherwise atomic numbers.
346 if ((strlen(tmpfield[5]) == 1) && isalpha(tmpfield[5][0]))
348 have_bonded_type = TRUE;
349 have_atomic_number = TRUE;
351 else if ((strlen(tmpfield[3]) == 1) && isalpha(tmpfield[3][0]))
353 have_bonded_type = FALSE;
354 have_atomic_number = FALSE;
358 have_bonded_type = (isalpha(tmpfield[1][0]) != 0);
359 have_atomic_number = !have_bonded_type;
362 /* optional fields */
371 if (have_atomic_number)
373 if (have_bonded_type)
375 nread = sscanf(line, "%s%s%d%lf%lf%s%lf%lf", type, btype, &atomnr, &m, &q,
376 ptype, &c[0], &c[1]);
385 /* have_atomic_number && !have_bonded_type */
386 nread = sscanf(line, "%s%d%lf%lf%s%lf%lf", type, &atomnr, &m, &q, ptype, &c[0], &c[1]);
396 if (have_bonded_type)
398 /* !have_atomic_number && have_bonded_type */
399 nread = sscanf(line, "%s%s%lf%lf%s%lf%lf", type, btype, &m, &q, ptype, &c[0], &c[1]);
408 /* !have_atomic_number && !have_bonded_type */
409 nread = sscanf(line, "%s%lf%lf%s%lf%lf", type, &m, &q, ptype, &c[0], &c[1]);
418 if (!have_bonded_type)
423 if (!have_atomic_number)
433 if (have_atomic_number)
435 if (have_bonded_type)
437 nread = sscanf(line, "%s%s%d%lf%lf%s%lf%lf%lf", type, btype, &atomnr, &m, &q,
438 ptype, &c[0], &c[1], &c[2]);
447 /* have_atomic_number && !have_bonded_type */
448 nread = sscanf(line, "%s%d%lf%lf%s%lf%lf%lf", type, &atomnr, &m, &q, ptype,
449 &c[0], &c[1], &c[2]);
459 if (have_bonded_type)
461 /* !have_atomic_number && have_bonded_type */
462 nread = sscanf(line, "%s%s%lf%lf%s%lf%lf%lf", type, btype, &m, &q, ptype, &c[0],
472 /* !have_atomic_number && !have_bonded_type */
473 nread = sscanf(line, "%s%lf%lf%s%lf%lf%lf", type, &m, &q, ptype, &c[0], &c[1], &c[2]);
482 if (!have_bonded_type)
487 if (!have_atomic_number)
495 auto message = gmx::formatString("Invalid function type %d in push_at", nb_funct);
496 warning_error_and_exit(wi, message, FARGS);
498 for (int j = nfp0; (j < MAXFORCEPARAM); j++)
502 std::array<real, MAXFORCEPARAM> forceParam;
504 if (strlen(type) == 1 && isdigit(type[0]))
506 warning_error_and_exit(wi, "Atom type names can't be single digits.", FARGS);
509 if (strlen(btype) == 1 && isdigit(btype[0]))
511 warning_error_and_exit(wi, "Bond atom type names can't be single digits.", FARGS);
514 /* Hack to read old topologies */
515 if (gmx_strcasecmp(ptype, "D") == 0)
519 for (j = 0; (j < eptNR); j++)
521 if (gmx_strcasecmp(ptype, xl[j].entry) == 0)
528 auto message = gmx::formatString("Invalid particle type %s", ptype);
529 warning_error_and_exit(wi, message, FARGS);
536 for (int i = 0; i < MAXFORCEPARAM; i++)
538 forceParam[i] = c[i];
541 InteractionOfType interactionType({}, forceParam, "");
543 batype_nr = bondAtomType->addBondAtomType(symtab, btype);
545 if ((nr = at->atomTypeFromName(type)) != NOTSET)
547 auto message = gmx::formatString(
548 "Atomtype %s was defined previously (e.g. in the forcefield files), "
549 "and has now been defined again. This could happen e.g. if you would "
550 "use a self-contained molecule .itp file that duplicates or replaces "
551 "the contents of the standard force-field files. You should check "
552 "the contents of your files and remove such repetition. If you know "
553 "you should override the previous definition, then you could choose "
554 "to suppress this warning with -maxwarn.",
556 warning(wi, message);
557 if ((nr = at->setType(nr, symtab, *atom, type, interactionType, batype_nr, atomnr)) == NOTSET)
559 auto message = gmx::formatString("Replacing atomtype %s failed", type);
560 warning_error_and_exit(wi, message, FARGS);
563 else if ((at->addType(symtab, *atom, type, interactionType, batype_nr, atomnr)) == NOTSET)
565 auto message = gmx::formatString("Adding atomtype %s failed", type);
566 warning_error_and_exit(wi, message, FARGS);
570 /* Add space in the non-bonded parameters matrix */
571 realloc_nb_params(at, nbparam, pair);
576 //! Return whether the contents of \c a and \c b are the same, considering also reversed order.
578 static bool equalEitherForwardOrBackward(gmx::ArrayRef<const T> a, gmx::ArrayRef<const T> b)
580 return (std::equal(a.begin(), a.end(), b.begin()) || std::equal(a.begin(), a.end(), b.rbegin()));
583 static void push_bondtype(InteractionsOfType* bt,
584 const InteractionOfType& b,
592 int nrfp = NRFP(ftype);
594 /* If bAllowRepeat is TRUE, we allow multiple entries as long as they
595 are on directly _adjacent_ lines.
598 /* First check if our atomtypes are _identical_ (not reversed) to the previous
599 entry. If they are not identical we search for earlier duplicates. If they are
600 we can skip it, since we already searched for the first line
604 bool isContinuationOfBlock = false;
605 if (bAllowRepeat && nr > 1)
607 isContinuationOfBlock = true;
608 gmx::ArrayRef<const int> newParAtom = b.atoms();
609 gmx::ArrayRef<const int> sysParAtom = bt->interactionTypes[nr - 2].atoms();
610 for (int j = 0; j < nral; j++)
612 if (newParAtom[j] != sysParAtom[j])
614 isContinuationOfBlock = false;
619 /* Search for earlier duplicates if this entry was not a continuation
620 from the previous line.
622 bool addBondType = true;
623 bool haveWarned = false;
624 bool haveErrored = false;
625 for (int i = 0; (i < nr); i++)
627 gmx::ArrayRef<const int> bParams = b.atoms();
628 gmx::ArrayRef<const int> testParams = bt->interactionTypes[i].atoms();
629 GMX_RELEASE_ASSERT(bParams.size() == testParams.size(),
630 "Number of atoms needs to be the same between parameters");
631 if (equalEitherForwardOrBackward(bParams, testParams))
633 GMX_ASSERT(nrfp <= MAXFORCEPARAM,
634 "This is ensured in other places, but we need this assert to keep the clang "
636 const bool identicalParameters = std::equal(
637 bt->interactionTypes[i].forceParam().begin(),
638 bt->interactionTypes[i].forceParam().begin() + nrfp, b.forceParam().begin());
640 if (!bAllowRepeat || identicalParameters)
645 if (!identicalParameters)
649 /* With dihedral type 9 we only allow for repeating
650 * of the same parameters with blocks with 1 entry.
651 * Allowing overriding is too complex to check.
653 if (!isContinuationOfBlock && !haveErrored)
656 "Encountered a second block of parameters for dihedral "
657 "type 9 for the same atoms, with either different parameters "
658 "and/or the first block has multiple lines. This is not "
663 else if (!haveWarned)
665 auto message = gmx::formatString(
666 "Bondtype %s was defined previously (e.g. in the forcefield files), "
667 "and has now been defined again. This could happen e.g. if you would "
668 "use a self-contained molecule .itp file that duplicates or replaces "
669 "the contents of the standard force-field files. You should check "
670 "the contents of your files and remove such repetition. If you know "
671 "you should override the previous definition, then you could choose "
672 "to suppress this warning with -maxwarn.%s",
673 interaction_function[ftype].longname,
674 (ftype == F_PDIHS) ? "\nUse dihedraltype 9 to allow several "
675 "multiplicity terms. Only consecutive "
676 "lines are combined. Non-consective lines "
677 "overwrite each other."
679 warning(wi, message);
681 fprintf(stderr, " old: ");
682 gmx::ArrayRef<const real> forceParam = bt->interactionTypes[i].forceParam();
683 for (int j = 0; j < nrfp; j++)
685 fprintf(stderr, " %g", forceParam[j]);
687 fprintf(stderr, " \n new: %s\n\n", line);
693 if (!identicalParameters && !bAllowRepeat)
695 /* Overwrite the parameters with the latest ones */
696 // TODO considering improving the following code by replacing with:
697 // std::copy(b->c, b->c + nrfp, bt->param[i].c);
698 gmx::ArrayRef<const real> forceParam = b.forceParam();
699 for (int j = 0; j < nrfp; j++)
701 bt->interactionTypes[i].setForceParameter(j, forceParam[j]);
709 /* fill the arrays up and down */
710 bt->interactionTypes.emplace_back(
711 InteractionOfType(b.atoms(), b.forceParam(), b.interactionTypeName()));
712 /* need to store force values because they might change below */
713 std::vector<real> forceParam(b.forceParam().begin(), b.forceParam().end());
715 /* The definitions of linear angles depend on the order of atoms,
716 * that means that for atoms i-j-k, with certain parameter a, the
717 * corresponding k-j-i angle will have parameter 1-a.
719 if (ftype == F_LINEAR_ANGLES)
721 forceParam[0] = 1 - forceParam[0];
722 forceParam[2] = 1 - forceParam[2];
724 std::vector<int> atoms;
725 gmx::ArrayRef<const int> oldAtoms = b.atoms();
726 for (auto oldAtom = oldAtoms.rbegin(); oldAtom != oldAtoms.rend(); oldAtom++)
728 atoms.emplace_back(*oldAtom);
730 bt->interactionTypes.emplace_back(InteractionOfType(atoms, forceParam, b.interactionTypeName()));
734 static std::vector<int> atomTypesFromAtomNames(const PreprocessingAtomTypes* atomTypes,
735 const PreprocessingBondAtomType* bondAtomTypes,
736 gmx::ArrayRef<const char[20]> atomNames,
740 GMX_RELEASE_ASSERT(!(!atomNames.empty() && !atomTypes && !bondAtomTypes),
741 "Need to have either valid atomtypes or bondatomtypes object");
743 std::vector<int> atomTypesFromAtomNames;
744 for (const auto& name : atomNames)
746 if (atomTypes != nullptr)
748 int atomType = atomTypes->atomTypeFromName(name);
749 if (atomType == NOTSET)
751 auto message = gmx::formatString("Unknown atomtype %s\n", name);
752 warning_error_and_exit(wi, message, FARGS);
754 atomTypesFromAtomNames.emplace_back(atomType);
756 else if (bondAtomTypes != nullptr)
758 int atomType = bondAtomTypes->bondAtomTypeFromName(name);
759 if (atomType == NOTSET)
761 auto message = gmx::formatString("Unknown bond_atomtype %s\n", name);
762 warning_error_and_exit(wi, message, FARGS);
764 atomTypesFromAtomNames.emplace_back(atomType);
767 return atomTypesFromAtomNames;
771 void push_bt(Directive d,
772 gmx::ArrayRef<InteractionsOfType> bt,
774 PreprocessingAtomTypes* at,
775 PreprocessingBondAtomType* bondAtomType,
779 const char* formal[MAXATOMLIST + 1] = {
780 "%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"
782 const char* formnl[MAXATOMLIST + 1] = { "%*s",
787 "%*s%*s%*s%*s%*s%*s",
788 "%*s%*s%*s%*s%*s%*s%*s" };
789 const char* formlf = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
790 int i, ft, ftype, nn, nrfp, nrfpA;
792 char alc[MAXATOMLIST + 1][20];
793 /* One force parameter more, so we can check if we read too many */
794 double c[MAXFORCEPARAM + 1];
796 if ((bondAtomType && at) || (!bondAtomType && !at))
798 gmx_incons("You should pass either bondAtomType or at to push_bt");
801 /* Make format string (nral ints+functype) */
802 if ((nn = sscanf(line, formal[nral], alc[0], alc[1], alc[2], alc[3], alc[4], alc[5])) != nral + 1)
804 auto message = gmx::formatString("Not enough atomtypes (%d instead of %d)", nn - 1, nral);
805 warning_error(wi, message);
809 ft = strtol(alc[nral], nullptr, 10);
810 ftype = ifunc_index(d, ft);
812 nrfpA = interaction_function[ftype].nrfpA;
813 strcpy(f1, formnl[nral]);
815 if ((nn = sscanf(line, f1, &c[0], &c[1], &c[2], &c[3], &c[4], &c[5], &c[6], &c[7], &c[8], &c[9],
816 &c[10], &c[11], &c[12]))
821 /* Copy the B-state from the A-state */
822 copy_B_from_A(ftype, c);
828 warning_error(wi, "Not enough parameters");
830 else if (nn > nrfpA && nn < nrfp)
832 warning_error(wi, "Too many parameters or not enough parameters for topology B");
836 warning_error(wi, "Too many parameters");
838 for (i = nn; (i < nrfp); i++)
844 std::vector<int> atomTypes =
845 atomTypesFromAtomNames(at, bondAtomType, gmx::arrayRefFromArray(alc, nral), wi);
846 std::array<real, MAXFORCEPARAM> forceParam;
847 for (int i = 0; (i < nrfp); i++)
849 forceParam[i] = c[i];
851 push_bondtype(&(bt[ftype]), InteractionOfType(atomTypes, forceParam), nral, ftype, FALSE, line, wi);
855 void push_dihedraltype(Directive d,
856 gmx::ArrayRef<InteractionsOfType> bt,
857 PreprocessingBondAtomType* bondAtomType,
861 const char* formal[MAXATOMLIST + 1] = {
862 "%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"
864 const char* formnl[MAXATOMLIST + 1] = { "%*s",
869 "%*s%*s%*s%*s%*s%*s",
870 "%*s%*s%*s%*s%*s%*s%*s" };
871 const char* formlf[MAXFORCEPARAM] = {
877 "%lf%lf%lf%lf%lf%lf",
878 "%lf%lf%lf%lf%lf%lf%lf",
879 "%lf%lf%lf%lf%lf%lf%lf%lf",
880 "%lf%lf%lf%lf%lf%lf%lf%lf%lf",
881 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
882 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
883 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
885 int i, ft, ftype, nn, nrfp, nrfpA, nral;
887 char alc[MAXATOMLIST + 1][20];
888 double c[MAXFORCEPARAM];
891 /* This routine accepts dihedraltypes defined from either 2 or 4 atoms.
893 * We first check for 2 atoms with the 3th column being an integer
894 * defining the type. If this isn't the case, we try it with 4 atoms
895 * and the 5th column defining the dihedral type.
897 nn = sscanf(line, formal[4], alc[0], alc[1], alc[2], alc[3], alc[4]);
898 if (nn >= 3 && strlen(alc[2]) == 1 && isdigit(alc[2][0]))
901 ft = strtol(alc[nral], nullptr, 10);
902 /* Move atom types around a bit and use 'X' for wildcard atoms
903 * to create a 4-atom dihedral definition with arbitrary atoms in
906 if (alc[2][0] == '2')
908 /* improper - the two atomtypes are 1,4. Use wildcards for 2,3 */
909 strcpy(alc[3], alc[1]);
910 sprintf(alc[2], "X");
911 sprintf(alc[1], "X");
912 /* alc[0] stays put */
916 /* proper - the two atomtypes are 2,3. Use wildcards for 1,4 */
917 sprintf(alc[3], "X");
918 strcpy(alc[2], alc[1]);
919 strcpy(alc[1], alc[0]);
920 sprintf(alc[0], "X");
923 else if (nn == 5 && strlen(alc[4]) == 1 && isdigit(alc[4][0]))
926 ft = strtol(alc[nral], nullptr, 10);
930 auto message = gmx::formatString(
931 "Incorrect number of atomtypes for dihedral (%d instead of 2 or 4)", nn - 1);
932 warning_error(wi, message);
938 /* Previously, we have always overwritten parameters if e.g. a torsion
939 with the same atomtypes occurs on multiple lines. However, CHARMM and
940 some other force fields specify multiple dihedrals over some bonds,
941 including cosines with multiplicity 6 and somethimes even higher.
942 Thus, they cannot be represented with Ryckaert-Bellemans terms.
943 To add support for these force fields, Dihedral type 9 is identical to
944 normal proper dihedrals, but repeated entries are allowed.
951 bAllowRepeat = FALSE;
955 ftype = ifunc_index(d, ft);
957 nrfpA = interaction_function[ftype].nrfpA;
959 strcpy(f1, formnl[nral]);
960 strcat(f1, formlf[nrfp - 1]);
962 /* Check number of parameters given */
963 if ((nn = sscanf(line, f1, &c[0], &c[1], &c[2], &c[3], &c[4], &c[5], &c[6], &c[7], &c[8], &c[9],
969 /* Copy the B-state from the A-state */
970 copy_B_from_A(ftype, c);
976 warning_error(wi, "Not enough parameters");
978 else if (nn > nrfpA && nn < nrfp)
980 warning_error(wi, "Too many parameters or not enough parameters for topology B");
984 warning_error(wi, "Too many parameters");
986 for (i = nn; (i < nrfp); i++)
993 std::vector<int> atoms;
994 std::array<real, MAXFORCEPARAM> forceParam;
995 for (int i = 0; (i < 4); i++)
997 if (!strcmp(alc[i], "X"))
999 atoms.emplace_back(-1);
1004 if ((atomNumber = bondAtomType->bondAtomTypeFromName(alc[i])) == NOTSET)
1006 auto message = gmx::formatString("Unknown bond_atomtype %s", alc[i]);
1007 warning_error_and_exit(wi, message, FARGS);
1009 atoms.emplace_back(atomNumber);
1012 for (int i = 0; (i < nrfp); i++)
1014 forceParam[i] = c[i];
1016 /* Always use 4 atoms here, since we created two wildcard atoms
1017 * if there wasn't of them 4 already.
1019 push_bondtype(&(bt[ftype]), InteractionOfType(atoms, forceParam), 4, ftype, bAllowRepeat, line, wi);
1023 void push_nbt(Directive d, t_nbparam** nbt, PreprocessingAtomTypes* atypes, char* pline, int nb_funct, warninp* wi)
1025 /* swap the atoms */
1026 const char* form3 = "%*s%*s%*s%lf%lf%lf";
1027 const char* form4 = "%*s%*s%*s%lf%lf%lf%lf";
1028 const char* form5 = "%*s%*s%*s%lf%lf%lf%lf%lf";
1029 char a0[80], a1[80];
1030 int i, f, n, ftype, nrfp;
1037 if (sscanf(pline, "%s%s%d", a0, a1, &f) != 3)
1043 ftype = ifunc_index(d, f);
1045 if (ftype != nb_funct)
1047 auto message = gmx::formatString("Trying to add %s while the default nonbond type is %s",
1048 interaction_function[ftype].longname,
1049 interaction_function[nb_funct].longname);
1050 warning_error(wi, message);
1054 /* Get the force parameters */
1056 if (ftype == F_LJ14)
1058 n = sscanf(pline, form4, &c[0], &c[1], &c[2], &c[3]);
1064 /* When the B topology parameters are not set,
1065 * copy them from topology A
1067 GMX_ASSERT(nrfp <= 4, "LJ-14 cannot have more than 4 parameters");
1068 for (i = n; i < nrfp; i++)
1073 else if (ftype == F_LJC14_Q)
1075 n = sscanf(pline, form5, &c[0], &c[1], &c[2], &c[3], &dum);
1078 incorrect_n_param(wi);
1084 if (sscanf(pline, form3, &c[0], &c[1], &dum) != 2)
1086 incorrect_n_param(wi);
1092 if (sscanf(pline, form4, &c[0], &c[1], &c[2], &dum) != 3)
1094 incorrect_n_param(wi);
1101 gmx::formatString("Number of force parameters for nonbonded interactions is %d", nrfp);
1102 warning_error_and_exit(wi, message, FARGS);
1104 for (i = 0; (i < nrfp); i++)
1109 /* Put the parameters in the matrix */
1110 if ((ai = atypes->atomTypeFromName(a0)) == NOTSET)
1112 auto message = gmx::formatString("Atomtype %s not found", a0);
1113 warning_error_and_exit(wi, message, FARGS);
1115 if ((aj = atypes->atomTypeFromName(a1)) == NOTSET)
1117 auto message = gmx::formatString("Atomtype %s not found", a1);
1118 warning_error_and_exit(wi, message, FARGS);
1120 nbp = &(nbt[std::max(ai, aj)][std::min(ai, aj)]);
1125 for (i = 0; i < nrfp; i++)
1127 bId = bId && (nbp->c[i] == cr[i]);
1131 auto message = gmx::formatString(
1132 "Non-bonded parameters were defined previously (e.g. in the forcefield files), "
1133 "and have now been defined again. This could happen e.g. if you would "
1134 "use a self-contained molecule .itp file that duplicates or replaces "
1135 "the contents of the standard force-field files. You should check "
1136 "the contents of your files and remove such repetition. If you know "
1137 "you should override the previous definitions, then you could choose "
1138 "to suppress this warning with -maxwarn.");
1139 warning(wi, message);
1140 fprintf(stderr, " old:");
1141 for (i = 0; i < nrfp; i++)
1143 fprintf(stderr, " %g", nbp->c[i]);
1145 fprintf(stderr, " new\n%s\n", pline);
1149 for (i = 0; i < nrfp; i++)
1155 void push_cmaptype(Directive d,
1156 gmx::ArrayRef<InteractionsOfType> bt,
1158 PreprocessingAtomTypes* atomtypes,
1159 PreprocessingBondAtomType* bondAtomType,
1163 const char* formal = "%s%s%s%s%s%s%s%s%n";
1165 int ft, ftype, nn, nrfp, nrfpA, nrfpB;
1166 int start, nchar_consumed;
1167 int nxcmap, nycmap, ncmap, read_cmap, sl, nct;
1168 char s[20], alc[MAXATOMLIST + 2][20];
1170 /* Keep the compiler happy */
1174 GMX_ASSERT(nral == 5, "CMAP requires 5 atoms per interaction");
1176 /* Here we can only check for < 8 */
1177 if ((nn = sscanf(line, formal, alc[0], alc[1], alc[2], alc[3], alc[4], alc[5], alc[6], alc[7], &nchar_consumed))
1181 gmx::formatString("Incorrect number of atomtypes for cmap (%d instead of 5)", nn - 1);
1182 warning_error(wi, message);
1185 start += nchar_consumed;
1187 ft = strtol(alc[nral], nullptr, 10);
1188 nxcmap = strtol(alc[nral + 1], nullptr, 10);
1189 nycmap = strtol(alc[nral + 2], nullptr, 10);
1191 /* Check for equal grid spacing in x and y dims */
1192 if (nxcmap != nycmap)
1194 auto message = gmx::formatString(
1195 "Not the same grid spacing in x and y for cmap grid: x=%d, y=%d", nxcmap, nycmap);
1196 warning_error(wi, message);
1199 ncmap = nxcmap * nycmap;
1200 ftype = ifunc_index(d, ft);
1201 nrfpA = strtol(alc[6], nullptr, 10) * strtol(alc[6], nullptr, 10);
1202 nrfpB = strtol(alc[7], nullptr, 10) * strtol(alc[7], nullptr, 10);
1203 nrfp = nrfpA + nrfpB;
1205 /* Read in CMAP parameters */
1207 for (int i = 0; i < ncmap; i++)
1209 while (isspace(*(line + start + sl)))
1213 nn = sscanf(line + start + sl, " %s ", s);
1215 bt[F_CMAP].cmap.emplace_back(strtod(s, nullptr));
1224 gmx::formatString("Error in reading cmap parameter for angle %s %s %s %s %s",
1225 alc[0], alc[1], alc[2], alc[3], alc[4]);
1226 warning_error(wi, message);
1230 /* Check do that we got the number of parameters we expected */
1231 if (read_cmap == nrfpA)
1233 for (int i = 0; i < ncmap; i++)
1235 bt[F_CMAP].cmap.emplace_back(bt[F_CMAP].cmap[i]);
1240 if (read_cmap < nrfpA)
1242 warning_error(wi, "Not enough cmap parameters");
1244 else if (read_cmap > nrfpA && read_cmap < nrfp)
1246 warning_error(wi, "Too many cmap parameters or not enough parameters for topology B");
1248 else if (read_cmap > nrfp)
1250 warning_error(wi, "Too many cmap parameters");
1255 /* Set grid spacing and the number of grids (we assume these numbers to be the same for all
1256 * grids so we can safely assign them each time
1258 bt[F_CMAP].cmakeGridSpacing = nxcmap; /* Or nycmap, they need to be equal */
1259 bt[F_CMAP].cmapAngles++; /* Since we are incrementing here, we need to subtract later, see (*****) */
1260 nct = (nral + 1) * bt[F_CMAP].cmapAngles;
1262 for (int i = 0; (i < nral); i++)
1264 /* Assign a grid number to each cmap_type */
1265 GMX_RELEASE_ASSERT(bondAtomType != nullptr, "Need valid PreprocessingBondAtomType object");
1266 bt[F_CMAP].cmapAtomTypes.emplace_back(bondAtomType->bondAtomTypeFromName(alc[i]));
1269 /* Assign a type number to this cmap */
1270 bt[F_CMAP].cmapAtomTypes.emplace_back(
1271 bt[F_CMAP].cmapAngles
1272 - 1); /* Since we inremented earlier, we need to subtrac here, to get the types right (****) */
1274 /* Check for the correct number of atoms (again) */
1275 if (bt[F_CMAP].nct() != nct)
1277 auto message = gmx::formatString("Incorrect number of atom types (%d) in cmap type %d\n",
1278 nct, bt[F_CMAP].cmapAngles);
1279 warning_error(wi, message);
1281 std::vector<int> atomTypes =
1282 atomTypesFromAtomNames(atomtypes, bondAtomType, gmx::constArrayRefFromArray(alc, nral), wi);
1283 std::array<real, MAXFORCEPARAM> forceParam = { NOTSET };
1285 /* Push the bond to the bondlist */
1286 push_bondtype(&(bt[ftype]), InteractionOfType(atomTypes, forceParam), nral, ftype, FALSE, line, wi);
1290 static void push_atom_now(t_symtab* symtab,
1308 int j, resind = 0, resnr;
1312 if (((nr == 0) && (atomnr != 1)) || (nr && (atomnr != at->nr + 1)))
1314 auto message = gmx::formatString(
1315 "Atoms in the .top are not numbered consecutively from 1 (rather, "
1316 "atomnr = %d, while at->nr = %d)",
1318 warning_error_and_exit(wi, message, FARGS);
1321 j = strlen(resnumberic) - 1;
1322 if (isdigit(resnumberic[j]))
1328 ric = resnumberic[j];
1329 if (j == 0 || !isdigit(resnumberic[j - 1]))
1332 gmx::formatString("Invalid residue number '%s' for atom %d", resnumberic, atomnr);
1333 warning_error_and_exit(wi, message, FARGS);
1336 resnr = strtol(resnumberic, nullptr, 10);
1340 resind = at->atom[nr - 1].resind;
1342 if (nr == 0 || strcmp(resname, *at->resinfo[resind].name) != 0
1343 || resnr != at->resinfo[resind].nr || ric != at->resinfo[resind].ic)
1353 at->nres = resind + 1;
1354 srenew(at->resinfo, at->nres);
1355 at->resinfo[resind].name = put_symtab(symtab, resname);
1356 at->resinfo[resind].nr = resnr;
1357 at->resinfo[resind].ic = ric;
1361 resind = at->atom[at->nr - 1].resind;
1364 /* New atom instance
1365 * get new space for arrays
1367 srenew(at->atom, nr + 1);
1368 srenew(at->atomname, nr + 1);
1369 srenew(at->atomtype, nr + 1);
1370 srenew(at->atomtypeB, nr + 1);
1373 at->atom[nr].type = type;
1374 at->atom[nr].ptype = ptype;
1375 at->atom[nr].q = q0;
1376 at->atom[nr].m = m0;
1377 at->atom[nr].typeB = typeB;
1378 at->atom[nr].qB = qB;
1379 at->atom[nr].mB = mB;
1381 at->atom[nr].resind = resind;
1382 at->atom[nr].atomnumber = atomicnumber;
1383 at->atomname[nr] = put_symtab(symtab, name);
1384 at->atomtype[nr] = put_symtab(symtab, ctype);
1385 at->atomtypeB[nr] = put_symtab(symtab, ctypeB);
1389 void push_atom(t_symtab* symtab, t_atoms* at, PreprocessingAtomTypes* atypes, char* line, warninp* wi)
1392 int cgnumber, atomnr, type, typeB, nscan;
1393 char id[STRLEN], ctype[STRLEN], ctypeB[STRLEN], resnumberic[STRLEN], resname[STRLEN],
1394 name[STRLEN], check[STRLEN];
1395 double m, q, mb, qb;
1396 real m0, q0, mB, qB;
1398 /* Fixed parameters */
1399 if (sscanf(line, "%s%s%s%s%s%d", id, ctype, resnumberic, resname, name, &cgnumber) != 6)
1404 sscanf(id, "%d", &atomnr);
1405 if ((type = atypes->atomTypeFromName(ctype)) == NOTSET)
1407 auto message = gmx::formatString("Atomtype %s not found", ctype);
1408 warning_error_and_exit(wi, message, FARGS);
1410 ptype = atypes->atomParticleTypeFromAtomType(type);
1412 /* Set default from type */
1413 q0 = atypes->atomChargeFromAtomType(type);
1414 m0 = atypes->atomMassFromAtomType(type);
1419 /* Optional parameters */
1420 nscan = sscanf(line, "%*s%*s%*s%*s%*s%*s%lf%lf%s%lf%lf%s", &q, &m, ctypeB, &qb, &mb, check);
1422 /* Nasty switch that falls thru all the way down! */
1431 if ((typeB = atypes->atomTypeFromName(ctypeB)) == NOTSET)
1433 auto message = gmx::formatString("Atomtype %s not found", ctypeB);
1434 warning_error_and_exit(wi, message, FARGS);
1436 qB = atypes->atomChargeFromAtomType(typeB);
1437 mB = atypes->atomMassFromAtomType(typeB);
1446 warning_error(wi, "Too many parameters");
1454 push_atom_now(symtab, at, atomnr, atypes->atomNumberFromAtomType(type), type, ctype, ptype,
1455 resnumberic, resname, name, m0, q0, typeB, typeB == type ? ctype : ctypeB, mB, qB, wi);
1458 void push_molt(t_symtab* symtab, std::vector<MoleculeInformation>* mol, char* line, warninp* wi)
1463 if ((sscanf(line, "%s%d", type, &nrexcl)) != 2)
1465 warning_error(wi, "Expected a molecule type name and nrexcl");
1468 /* Test if this moleculetype overwrites another */
1469 const auto found = std::find_if(mol->begin(), mol->end(),
1470 [&type](const auto& m) { return strcmp(*(m.name), type) == 0; });
1471 if (found != mol->end())
1473 auto message = gmx::formatString("moleculetype %s is redefined", type);
1474 warning_error_and_exit(wi, message, FARGS);
1477 mol->emplace_back();
1478 mol->back().initMolInfo();
1480 /* Fill in the values */
1481 mol->back().name = put_symtab(symtab, type);
1482 mol->back().nrexcl = nrexcl;
1483 mol->back().excl_set = false;
1486 static bool findIfAllNBAtomsMatch(gmx::ArrayRef<const int> atomsFromParameterArray,
1487 gmx::ArrayRef<const int> atomsFromCurrentParameter,
1491 if (atomsFromParameterArray.size() != atomsFromCurrentParameter.size())
1497 for (gmx::index i = 0; i < atomsFromCurrentParameter.ssize(); i++)
1499 if (at->atom[atomsFromCurrentParameter[i]].typeB != atomsFromParameterArray[i])
1508 for (gmx::index i = 0; i < atomsFromCurrentParameter.ssize(); i++)
1510 if (at->atom[atomsFromCurrentParameter[i]].type != atomsFromParameterArray[i])
1519 static bool default_nb_params(int ftype,
1520 gmx::ArrayRef<InteractionsOfType> bt,
1522 InteractionOfType* p,
1529 InteractionOfType* pi = nullptr;
1530 int nr = bt[ftype].size();
1531 int nral = NRAL(ftype);
1532 int nrfp = interaction_function[ftype].nrfpA;
1533 int nrfpB = interaction_function[ftype].nrfpB;
1535 if ((!bB && nrfp == 0) || (bB && nrfpB == 0))
1543 /* First test the generated-pair position to save
1544 * time when we have 1000*1000 entries for e.g. OPLS...
1546 ntype = static_cast<int>(std::sqrt(static_cast<double>(nr)));
1547 GMX_ASSERT(ntype * ntype == nr,
1548 "Number of pairs of generated non-bonded parameters should be a perfect square");
1551 ti = at->atom[p->ai()].typeB;
1552 tj = at->atom[p->aj()].typeB;
1556 ti = at->atom[p->ai()].type;
1557 tj = at->atom[p->aj()].type;
1559 pi = &(bt[ftype].interactionTypes[ntype * ti + tj]);
1560 if (pi->atoms().ssize() < nral)
1562 /* not initialized yet with atom names */
1567 bFound = ((ti == pi->ai()) && (tj == pi->aj()));
1571 gmx::ArrayRef<const int> paramAtoms = p->atoms();
1572 /* Search explicitly if we didnt find it */
1575 auto foundParameter =
1576 std::find_if(bt[ftype].interactionTypes.begin(), bt[ftype].interactionTypes.end(),
1577 [¶mAtoms, &at, &bB](const auto& param) {
1578 return findIfAllNBAtomsMatch(param.atoms(), paramAtoms, at, bB);
1580 if (foundParameter != bt[ftype].interactionTypes.end())
1583 pi = &(*foundParameter);
1589 gmx::ArrayRef<const real> forceParam = pi->forceParam();
1592 if (nrfp + nrfpB > MAXFORCEPARAM)
1594 gmx_incons("Too many force parameters");
1596 for (int j = c_start; j < nrfpB; j++)
1598 p->setForceParameter(nrfp + j, forceParam[j]);
1603 for (int j = c_start; j < nrfp; j++)
1605 p->setForceParameter(j, forceParam[j]);
1611 for (int j = c_start; j < nrfp; j++)
1613 p->setForceParameter(j, 0.0);
1619 static bool default_cmap_params(gmx::ArrayRef<InteractionsOfType> bondtype,
1621 PreprocessingAtomTypes* atypes,
1622 InteractionOfType* p,
1630 bool bFound = false;
1635 /* Match the current cmap angle against the list of cmap_types */
1636 for (int i = 0; i < bondtype[F_CMAP].nct() && !bFound; i += 6)
1641 if ((atypes->bondAtomTypeFromAtomType(at->atom[p->ai()].type)
1642 == bondtype[F_CMAP].cmapAtomTypes[i])
1643 && (atypes->bondAtomTypeFromAtomType(at->atom[p->aj()].type)
1644 == bondtype[F_CMAP].cmapAtomTypes[i + 1])
1645 && (atypes->bondAtomTypeFromAtomType(at->atom[p->ak()].type)
1646 == bondtype[F_CMAP].cmapAtomTypes[i + 2])
1647 && (atypes->bondAtomTypeFromAtomType(at->atom[p->al()].type)
1648 == bondtype[F_CMAP].cmapAtomTypes[i + 3])
1649 && (atypes->bondAtomTypeFromAtomType(at->atom[p->am()].type)
1650 == bondtype[F_CMAP].cmapAtomTypes[i + 4]))
1652 /* Found cmap torsion */
1654 ct = bondtype[F_CMAP].cmapAtomTypes[i + 5];
1660 /* If we did not find a matching type for this cmap torsion */
1663 auto message = gmx::formatString("Unknown cmap torsion between atoms %d %d %d %d %d", p->ai() + 1,
1664 p->aj() + 1, p->ak() + 1, p->al() + 1, p->am() + 1);
1665 warning_error_and_exit(wi, message, FARGS);
1668 *nparam_def = nparam_found;
1674 /* Returns the number of exact atom type matches, i.e. non wild-card matches,
1675 * returns -1 when there are no matches at all.
1677 static int natom_match(const InteractionOfType& pi,
1682 const PreprocessingAtomTypes* atypes)
1684 if ((pi.ai() == -1 || atypes->bondAtomTypeFromAtomType(type_i) == pi.ai())
1685 && (pi.aj() == -1 || atypes->bondAtomTypeFromAtomType(type_j) == pi.aj())
1686 && (pi.ak() == -1 || atypes->bondAtomTypeFromAtomType(type_k) == pi.ak())
1687 && (pi.al() == -1 || atypes->bondAtomTypeFromAtomType(type_l) == pi.al()))
1689 return (pi.ai() == -1 ? 0 : 1) + (pi.aj() == -1 ? 0 : 1) + (pi.ak() == -1 ? 0 : 1)
1690 + (pi.al() == -1 ? 0 : 1);
1698 static int findNumberOfDihedralAtomMatches(const InteractionOfType& currentParamFromParameterArray,
1699 const InteractionOfType& parameterToAdd,
1701 const PreprocessingAtomTypes* atypes,
1706 return natom_match(currentParamFromParameterArray, at->atom[parameterToAdd.ai()].typeB,
1707 at->atom[parameterToAdd.aj()].typeB, at->atom[parameterToAdd.ak()].typeB,
1708 at->atom[parameterToAdd.al()].typeB, atypes);
1712 return natom_match(currentParamFromParameterArray, at->atom[parameterToAdd.ai()].type,
1713 at->atom[parameterToAdd.aj()].type, at->atom[parameterToAdd.ak()].type,
1714 at->atom[parameterToAdd.al()].type, atypes);
1718 static bool findIfAllParameterAtomsMatch(gmx::ArrayRef<const int> atomsFromParameterArray,
1719 gmx::ArrayRef<const int> atomsFromCurrentParameter,
1721 const PreprocessingAtomTypes* atypes,
1724 if (atomsFromParameterArray.size() != atomsFromCurrentParameter.size())
1730 for (gmx::index i = 0; i < atomsFromCurrentParameter.ssize(); i++)
1732 if (atypes->bondAtomTypeFromAtomType(at->atom[atomsFromCurrentParameter[i]].typeB)
1733 != atomsFromParameterArray[i])
1742 for (gmx::index i = 0; i < atomsFromCurrentParameter.ssize(); i++)
1744 if (atypes->bondAtomTypeFromAtomType(at->atom[atomsFromCurrentParameter[i]].type)
1745 != atomsFromParameterArray[i])
1754 static std::vector<InteractionOfType>::iterator defaultInteractionsOfType(int ftype,
1755 gmx::ArrayRef<InteractionsOfType> bt,
1757 PreprocessingAtomTypes* atypes,
1758 const InteractionOfType& p,
1763 int nrfpA = interaction_function[ftype].nrfpA;
1764 int nrfpB = interaction_function[ftype].nrfpB;
1766 if ((!bB && nrfpA == 0) || (bB && nrfpB == 0))
1768 return bt[ftype].interactionTypes.end();
1773 if (ftype == F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS || ftype == F_PIDIHS)
1775 int nmatch_max = -1;
1777 /* For dihedrals we allow wildcards. We choose the first type
1778 * that has the most real matches, i.e. non-wildcard matches.
1780 auto prevPos = bt[ftype].interactionTypes.end();
1781 auto pos = bt[ftype].interactionTypes.begin();
1782 while (pos != bt[ftype].interactionTypes.end() && nmatch_max < 4)
1784 pos = std::find_if(bt[ftype].interactionTypes.begin(), bt[ftype].interactionTypes.end(),
1785 [&p, &at, &atypes, &bB, &nmatch_max](const auto& param) {
1786 return (findNumberOfDihedralAtomMatches(param, p, at, atypes, bB)
1789 if (pos != bt[ftype].interactionTypes.end())
1792 nmatch_max = findNumberOfDihedralAtomMatches(*pos, p, at, atypes, bB);
1796 if (prevPos != bt[ftype].interactionTypes.end())
1800 /* Find additional matches for this dihedral - necessary
1802 * The rule in that case is that additional matches
1803 * HAVE to be on adjacent lines!
1806 // Advance iterator (like std::advance) without incrementing past end (UB)
1807 const auto safeAdvance = [](auto& it, auto n, auto end) {
1808 it = end - it > n ? it + n : end;
1810 /* Continue from current iterator position */
1811 auto nextPos = prevPos;
1812 const auto endIter = bt[ftype].interactionTypes.end();
1813 safeAdvance(nextPos, 2, endIter);
1814 for (; nextPos < endIter && bSame; safeAdvance(nextPos, 2, endIter))
1816 bSame = (prevPos->ai() == nextPos->ai() && prevPos->aj() == nextPos->aj()
1817 && prevPos->ak() == nextPos->ak() && prevPos->al() == nextPos->al());
1822 /* nparam_found will be increased as long as the numbers match */
1825 *nparam_def = nparam_found;
1828 else /* Not a dihedral */
1830 gmx::ArrayRef<const int> atomParam = p.atoms();
1831 auto found = std::find_if(
1832 bt[ftype].interactionTypes.begin(), bt[ftype].interactionTypes.end(),
1833 [&atomParam, &at, &atypes, &bB](const auto& param) {
1834 return findIfAllParameterAtomsMatch(param.atoms(), atomParam, at, atypes, bB);
1836 if (found != bt[ftype].interactionTypes.end())
1840 *nparam_def = nparam_found;
1846 void push_bond(Directive d,
1847 gmx::ArrayRef<InteractionsOfType> bondtype,
1848 gmx::ArrayRef<InteractionsOfType> bond,
1850 PreprocessingAtomTypes* atypes,
1856 bool* bWarn_copy_A_B,
1859 const char* aaformat[MAXATOMLIST] = { "%d%d", "%d%d%d", "%d%d%d%d",
1860 "%d%d%d%d%d", "%d%d%d%d%d%d", "%d%d%d%d%d%d%d" };
1861 const char* asformat[MAXATOMLIST] = {
1862 "%*s%*s", "%*s%*s%*s", "%*s%*s%*s%*s",
1863 "%*s%*s%*s%*s%*s", "%*s%*s%*s%*s%*s%*s", "%*s%*s%*s%*s%*s%*s%*s"
1865 const char* ccformat = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
1866 int nral, nral_fmt, nread, ftype;
1867 char format[STRLEN];
1868 /* One force parameter more, so we can check if we read too many */
1869 double cc[MAXFORCEPARAM + 1];
1870 int aa[MAXATOMLIST + 1];
1871 bool bFoundA = FALSE, bFoundB = FALSE, bDef, bSwapParity = FALSE;
1872 int nparam_defA, nparam_defB;
1874 nparam_defA = nparam_defB = 0;
1876 ftype = ifunc_index(d, 1);
1878 for (int j = 0; j < nral; j++)
1882 bDef = (NRFP(ftype) > 0);
1884 if (ftype == F_SETTLE)
1886 /* SETTLE acts on 3 atoms, but the topology format only specifies
1887 * the first atom (for historical reasons).
1896 nread = sscanf(line, aaformat[nral_fmt - 1], &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
1898 if (ftype == F_SETTLE)
1905 if (nread < nral_fmt)
1910 else if (nread > nral_fmt)
1912 /* this is a hack to allow for virtual sites with swapped parity */
1913 bSwapParity = (aa[nral] < 0);
1916 aa[nral] = -aa[nral];
1918 ftype = ifunc_index(d, aa[nral]);
1924 case F_VSITE3OUT: break;
1927 gmx::formatString("Negative function types only allowed for %s and %s",
1928 interaction_function[F_VSITE3FAD].longname,
1929 interaction_function[F_VSITE3OUT].longname);
1930 warning_error_and_exit(wi, message, FARGS);
1936 /* Check for double atoms and atoms out of bounds */
1937 for (int i = 0; (i < nral); i++)
1939 if (aa[i] < 1 || aa[i] > at->nr)
1941 auto message = gmx::formatString(
1942 "Atom index (%d) in %s out of bounds (1-%d).\n"
1943 "This probably means that you have inserted topology section \"%s\"\n"
1944 "in a part belonging to a different molecule than you intended to.\n"
1945 "In that case move the \"%s\" section to the right molecule.",
1946 aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
1947 warning_error_and_exit(wi, message, FARGS);
1949 for (int j = i + 1; (j < nral); j++)
1951 GMX_ASSERT(j < MAXATOMLIST + 1,
1952 "Values from nral=NRAL() will satisfy this, we assert to keep gcc 4 happy");
1955 auto message = gmx::formatString("Duplicate atom index (%d) in %s", aa[i], dir2str(d));
1956 if (ftype == F_ANGRES)
1958 /* Since the angle restraints uses 2 pairs of atoms to
1959 * defines an angle between vectors, it can be useful
1960 * to use one atom twice, so we only issue a note here.
1962 warning_note(wi, message);
1966 warning_error(wi, message);
1972 /* default force parameters */
1973 std::vector<int> atoms;
1974 for (int j = 0; (j < nral); j++)
1976 atoms.emplace_back(aa[j] - 1);
1978 /* need to have an empty but initialized param array for some reason */
1979 std::array<real, MAXFORCEPARAM> forceParam = { 0.0 };
1981 /* Get force params for normal and free energy perturbation
1982 * studies, as determined by types!
1984 InteractionOfType param(atoms, forceParam, "");
1986 std::vector<InteractionOfType>::iterator foundAParameter = bondtype[ftype].interactionTypes.end();
1987 std::vector<InteractionOfType>::iterator foundBParameter = bondtype[ftype].interactionTypes.end();
1991 defaultInteractionsOfType(ftype, bondtype, at, atypes, param, FALSE, &nparam_defA);
1992 if (foundAParameter != bondtype[ftype].interactionTypes.end())
1994 /* Copy the A-state and B-state default parameters. */
1995 GMX_ASSERT(NRFPA(ftype) + NRFPB(ftype) <= MAXFORCEPARAM,
1996 "Bonded interactions may have at most 12 parameters");
1997 gmx::ArrayRef<const real> defaultParam = foundAParameter->forceParam();
1998 for (int j = 0; (j < NRFPA(ftype) + NRFPB(ftype)); j++)
2000 param.setForceParameter(j, defaultParam[j]);
2004 else if (NRFPA(ftype) == 0)
2009 defaultInteractionsOfType(ftype, bondtype, at, atypes, param, TRUE, &nparam_defB);
2010 if (foundBParameter != bondtype[ftype].interactionTypes.end())
2012 /* Copy only the B-state default parameters */
2013 gmx::ArrayRef<const real> defaultParam = foundBParameter->forceParam();
2014 for (int j = NRFPA(ftype); (j < NRFP(ftype)); j++)
2016 param.setForceParameter(j, defaultParam[j]);
2020 else if (NRFPB(ftype) == 0)
2025 else if (ftype == F_LJ14)
2027 bFoundA = default_nb_params(ftype, bondtype, at, ¶m, 0, FALSE, bGenPairs);
2028 bFoundB = default_nb_params(ftype, bondtype, at, ¶m, 0, TRUE, bGenPairs);
2030 else if (ftype == F_LJC14_Q)
2032 /* Fill in the A-state charges as default parameters */
2033 param.setForceParameter(0, fudgeQQ);
2034 param.setForceParameter(1, at->atom[param.ai()].q);
2035 param.setForceParameter(2, at->atom[param.aj()].q);
2036 /* The default LJ parameters are the standard 1-4 parameters */
2037 bFoundA = default_nb_params(F_LJ14, bondtype, at, ¶m, 3, FALSE, bGenPairs);
2040 else if (ftype == F_LJC_PAIRS_NB)
2042 /* Defaults are not supported here */
2048 gmx_incons("Unknown function type in push_bond");
2051 if (nread > nral_fmt)
2053 /* Manually specified parameters - in this case we discard multiple torsion info! */
2055 strcpy(format, asformat[nral_fmt - 1]);
2056 strcat(format, ccformat);
2058 nread = sscanf(line, format, &cc[0], &cc[1], &cc[2], &cc[3], &cc[4], &cc[5], &cc[6], &cc[7],
2059 &cc[8], &cc[9], &cc[10], &cc[11], &cc[12]);
2061 if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0))
2063 /* We only have to issue a warning if these atoms are perturbed! */
2065 gmx::ArrayRef<const int> paramAtoms = param.atoms();
2066 for (int j = 0; (j < nral); j++)
2068 bPert = bPert || PERTURBED(at->atom[paramAtoms[j]]);
2071 if (bPert && *bWarn_copy_A_B)
2073 auto message = gmx::formatString(
2074 "Some parameters for bonded interaction involving "
2075 "perturbed atoms are specified explicitly in "
2076 "state A, but not B - copying A to B");
2077 warning(wi, message);
2078 *bWarn_copy_A_B = FALSE;
2081 /* If only the A parameters were specified, copy them to the B state */
2082 /* The B-state parameters correspond to the first nrfpB
2083 * A-state parameters.
2085 for (int j = 0; (j < NRFPB(ftype)); j++)
2087 cc[nread++] = cc[j];
2091 /* If nread was 0 or EOF, no parameters were read => use defaults.
2092 * If nread was nrfpA we copied above so nread=nrfp.
2093 * If nread was nrfp we are cool.
2094 * For F_LJC14_Q we allow supplying fudgeQQ only.
2095 * Anything else is an error!
2097 if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) && !(ftype == F_LJC14_Q && nread == 1))
2099 auto message = gmx::formatString(
2100 "Incorrect number of parameters - found %d, expected %d "
2101 "or %d for %s (after the function type).",
2102 nread, NRFPA(ftype), NRFP(ftype), interaction_function[ftype].longname);
2103 warning_error_and_exit(wi, message, FARGS);
2106 for (int j = 0; (j < nread); j++)
2108 param.setForceParameter(j, cc[j]);
2110 /* Check whether we have to use the defaults */
2111 if (nread == NRFP(ftype))
2120 /* nread now holds the number of force parameters read! */
2125 /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */
2126 if (ftype == F_PDIHS)
2128 if ((nparam_defA != nparam_defB)
2129 || ((nparam_defA > 1 || nparam_defB > 1) && (foundAParameter != foundBParameter)))
2131 auto message = gmx::formatString(
2132 "Cannot automatically perturb a torsion with multiple terms to different "
2134 "Please specify perturbed parameters manually for this torsion in your "
2136 warning_error(wi, message);
2140 if (nread > 0 && nread < NRFPA(ftype))
2142 /* Issue an error, do not use defaults */
2143 auto message = gmx::formatString(
2144 "Not enough parameters, there should be at least %d (or 0 for defaults)", NRFPA(ftype));
2145 warning_error(wi, message);
2148 if (nread == 0 || nread == EOF)
2152 if (interaction_function[ftype].flags & IF_VSITE)
2154 for (int j = 0; j < MAXFORCEPARAM; j++)
2156 param.setForceParameter(j, NOTSET);
2160 /* flag to swap parity of vsi te construction */
2161 param.setForceParameter(1, -1);
2168 fprintf(stderr, "NOTE: No default %s types, using zeroes\n",
2169 interaction_function[ftype].longname);
2173 auto message = gmx::formatString("No default %s types",
2174 interaction_function[ftype].longname);
2175 warning_error(wi, message);
2185 case F_VSITE3FAD: param.setForceParameter(0, 360 - param.c0()); break;
2186 case F_VSITE3OUT: param.setForceParameter(2, -param.c2()); break;
2192 /* We only have to issue a warning if these atoms are perturbed! */
2194 gmx::ArrayRef<const int> paramAtoms = param.atoms();
2195 for (int j = 0; (j < nral); j++)
2197 bPert = bPert || PERTURBED(at->atom[paramAtoms[j]]);
2202 auto message = gmx::formatString(
2203 "No default %s types for perturbed atoms, "
2204 "using normal values",
2205 interaction_function[ftype].longname);
2206 warning(wi, message);
2212 gmx::ArrayRef<const real> paramValue = param.forceParam();
2213 if ((ftype == F_PDIHS || ftype == F_ANGRES || ftype == F_ANGRESZ) && paramValue[5] != paramValue[2])
2215 auto message = gmx::formatString("%s multiplicity can not be perturbed %f!=%f",
2216 interaction_function[ftype].longname, paramValue[2],
2218 warning_error_and_exit(wi, message, FARGS);
2221 if (IS_TABULATED(ftype) && param.c0() != param.c2())
2223 auto message = gmx::formatString("%s table number can not be perturbed %d!=%d",
2224 interaction_function[ftype].longname,
2225 gmx::roundToInt(param.c0()), gmx::roundToInt(param.c0()));
2226 warning_error_and_exit(wi, message, FARGS);
2229 /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
2230 if (ftype == F_RBDIHS)
2234 for (int i = 0; i < NRFP(ftype); i++)
2236 if (paramValue[i] != 0.0)
2247 /* Put the values in the appropriate arrays */
2248 add_param_to_list(&bond[ftype], param);
2250 /* Push additional torsions from FF for ftype==9 if we have them.
2251 * We have already checked that the A/B states do not differ in this case,
2252 * so we do not have to double-check that again, or the vsite stuff.
2253 * In addition, those torsions cannot be automatically perturbed.
2255 if (bDef && ftype == F_PDIHS)
2257 for (int i = 1; i < nparam_defA; i++)
2259 /* Advance pointer! */
2260 foundAParameter += 2;
2261 gmx::ArrayRef<const real> forceParam = foundAParameter->forceParam();
2262 for (int j = 0; j < (NRFPA(ftype) + NRFPB(ftype)); j++)
2264 param.setForceParameter(j, forceParam[j]);
2266 /* And push the next term for this torsion */
2267 add_param_to_list(&bond[ftype], param);
2272 void push_cmap(Directive d,
2273 gmx::ArrayRef<InteractionsOfType> bondtype,
2274 gmx::ArrayRef<InteractionsOfType> bond,
2276 PreprocessingAtomTypes* atypes,
2280 const char* aaformat[MAXATOMLIST + 1] = {
2281 "%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"
2284 int ftype, nral, nread, ncmap_params;
2286 int aa[MAXATOMLIST + 1];
2289 ftype = ifunc_index(d, 1);
2293 nread = sscanf(line, aaformat[nral - 1], &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
2300 else if (nread == nral)
2302 ftype = ifunc_index(d, 1);
2305 /* Check for double atoms and atoms out of bounds */
2306 for (int i = 0; i < nral; i++)
2308 if (aa[i] < 1 || aa[i] > at->nr)
2310 auto message = gmx::formatString(
2311 "Atom index (%d) in %s out of bounds (1-%d).\n"
2312 "This probably means that you have inserted topology section \"%s\"\n"
2313 "in a part belonging to a different molecule than you intended to.\n"
2314 "In that case move the \"%s\" section to the right molecule.",
2315 aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
2316 warning_error_and_exit(wi, message, FARGS);
2319 for (int j = i + 1; (j < nral); j++)
2323 auto message = gmx::formatString("Duplicate atom index (%d) in %s", aa[i], dir2str(d));
2324 warning_error(wi, message);
2329 /* default force parameters */
2330 std::vector<int> atoms;
2331 for (int j = 0; (j < nral); j++)
2333 atoms.emplace_back(aa[j] - 1);
2335 std::array<real, MAXFORCEPARAM> forceParam = { 0.0 };
2336 InteractionOfType param(atoms, forceParam, "");
2337 /* Get the cmap type for this cmap angle */
2338 bFound = default_cmap_params(bondtype, at, atypes, ¶m, FALSE, &cmap_type, &ncmap_params, wi);
2340 /* We want exactly one parameter (the cmap type in state A (currently no state B) back */
2341 if (bFound && ncmap_params == 1)
2343 /* Put the values in the appropriate arrays */
2344 param.setForceParameter(0, cmap_type);
2345 add_param_to_list(&bond[ftype], param);
2349 /* This is essentially the same check as in default_cmap_params() done one more time */
2350 auto message = gmx::formatString(
2351 "Unable to assign a cmap type to torsion %d %d %d %d and %d\n", param.ai() + 1,
2352 param.aj() + 1, param.ak() + 1, param.al() + 1, param.am() + 1);
2353 warning_error_and_exit(wi, message, FARGS);
2358 void push_vsitesn(Directive d, gmx::ArrayRef<InteractionsOfType> bond, t_atoms* at, char* line, warninp* wi)
2361 int type, ftype, n, ret, nj, a;
2363 double *weight = nullptr, weight_tot;
2365 std::array<real, MAXFORCEPARAM> forceParam = { 0.0 };
2367 ret = sscanf(ptr, "%d%n", &a, &n);
2371 auto message = gmx::formatString("Expected an atom index in section \"%s\"", dir2str(d));
2372 warning_error_and_exit(wi, message, FARGS);
2375 sscanf(ptr, "%d%n", &type, &n);
2377 ftype = ifunc_index(d, type);
2378 int firstAtom = a - 1;
2384 ret = sscanf(ptr, "%d%n", &a, &n);
2390 srenew(atc, nj + 20);
2391 srenew(weight, nj + 20);
2396 case 1: weight[nj] = 1; break;
2398 /* Here we use the A-state mass as a parameter.
2399 * Note that the B-state mass has no influence.
2401 weight[nj] = at->atom[atc[nj]].m;
2405 ret = sscanf(ptr, "%lf%n", &(weight[nj]), &n);
2409 auto message = gmx::formatString(
2410 "No weight or negative weight found for vsiten "
2411 "constructing atom %d (atom index %d)",
2412 nj + 1, atc[nj] + 1);
2413 warning_error_and_exit(wi, message, FARGS);
2417 auto message = gmx::formatString("Unknown vsiten type %d", type);
2418 warning_error_and_exit(wi, message, FARGS);
2420 weight_tot += weight[nj];
2428 gmx::formatString("Expected more than one atom index in section \"%s\"", dir2str(d));
2429 warning_error_and_exit(wi, message, FARGS);
2432 if (weight_tot == 0)
2434 warning_error_and_exit(wi, "The total mass of the construting atoms is zero", FARGS);
2437 for (int j = 0; j < nj; j++)
2439 std::vector<int> atoms = { firstAtom, atc[j] };
2441 forceParam[1] = weight[j] / weight_tot;
2442 /* Put the values in the appropriate arrays */
2443 add_param_to_list(&bond[ftype], InteractionOfType(atoms, forceParam));
2450 void push_mol(gmx::ArrayRef<MoleculeInformation> mols, char* pline, int* whichmol, int* nrcopies, warninp* wi)
2454 if (sscanf(pline, "%s%d", type, nrcopies) != 2)
2460 /* Search moleculename.
2461 * Here we originally only did case insensitive matching. But because
2462 * some PDB files can have many chains and use case to generate more
2463 * chain-identifiers, which in turn end up in our moleculetype name,
2464 * we added support for case-sensitivity.
2471 for (const auto& mol : mols)
2473 if (strcmp(type, *(mol.name)) == 0)
2478 if (gmx_strcasecmp(type, *(mol.name)) == 0)
2488 // select the case sensitive match
2489 *whichmol = matchcs;
2493 // avoid matching case-insensitive when we have multiple matches
2496 auto message = gmx::formatString(
2497 "For moleculetype '%s' in [ system ] %d case insensitive "
2498 "matches, but %d case sensitive matches were found. Check "
2499 "the case of the characters in the moleculetypes.",
2501 warning_error_and_exit(wi, message, FARGS);
2505 // select the unique case insensitive match
2506 *whichmol = matchci;
2510 auto message = gmx::formatString("No such moleculetype %s", type);
2511 warning_error_and_exit(wi, message, FARGS);
2516 void push_excl(char* line, gmx::ArrayRef<gmx::ExclusionBlock> b2, warninp* wi)
2520 char base[STRLEN], format[STRLEN];
2522 if (sscanf(line, "%d", &i) == 0)
2527 if ((1 <= i) && (i <= b2.ssize()))
2535 strcpy(base, "%*d");
2538 strcpy(format, base);
2539 strcat(format, "%d");
2540 n = sscanf(line, format, &j);
2543 if ((1 <= j) && (j <= b2.ssize()))
2546 b2[i].atomNumber.push_back(j);
2547 /* also add the reverse exclusion! */
2548 b2[j].atomNumber.push_back(i);
2549 strcat(base, "%*d");
2553 auto message = gmx::formatString("Invalid Atomnr j: %d, b2->nr: %zu\n", j, b2.size());
2554 warning_error_and_exit(wi, message, FARGS);
2560 int add_atomtype_decoupled(t_symtab* symtab, PreprocessingAtomTypes* at, t_nbparam*** nbparam, t_nbparam*** pair)
2565 /* Define an atom type with all parameters set to zero (no interactions) */
2568 /* Type for decoupled atoms could be anything,
2569 * this should be changed automatically later when required.
2571 atom.ptype = eptAtom;
2573 std::array<real, MAXFORCEPARAM> forceParam = { 0.0 };
2574 nr = at->addType(symtab, atom, "decoupled", InteractionOfType({}, forceParam, ""), -1, 0);
2576 /* Add space in the non-bonded parameters matrix */
2577 realloc_nb_params(at, nbparam, pair);
2582 static void convert_pairs_to_pairsQ(gmx::ArrayRef<InteractionsOfType> interactions, real fudgeQQ, t_atoms* atoms)
2584 /* Add the pair list to the pairQ list */
2585 std::vector<InteractionOfType> paramnew;
2587 gmx::ArrayRef<const InteractionOfType> paramp1 = interactions[F_LJ14].interactionTypes;
2588 gmx::ArrayRef<const InteractionOfType> paramp2 = interactions[F_LJC14_Q].interactionTypes;
2590 /* Fill in the new F_LJC14_Q array with the old one. NOTE:
2591 it may be possible to just ADD the converted F_LJ14 array
2592 to the old F_LJC14_Q array, but since we have to create
2593 a new sized memory structure, better just to deep copy it all.
2597 for (const auto& param : paramp2)
2599 paramnew.emplace_back(param);
2602 for (const auto& param : paramp1)
2604 std::vector<real> forceParam = { fudgeQQ, atoms->atom[param.ai()].q,
2605 atoms->atom[param.aj()].q, param.c0(), param.c1() };
2606 paramnew.emplace_back(InteractionOfType(param.atoms(), forceParam, ""));
2609 /* now assign the new data to the F_LJC14_Q structure */
2610 interactions[F_LJC14_Q].interactionTypes = paramnew;
2612 /* Empty the LJ14 pairlist */
2613 interactions[F_LJ14].interactionTypes.clear();
2616 static void generate_LJCpairsNB(MoleculeInformation* mol, int nb_funct, InteractionsOfType* nbp, warninp* wi)
2622 atom = mol->atoms.atom;
2624 ntype = static_cast<int>(std::sqrt(static_cast<double>(nbp->size())));
2625 GMX_ASSERT(ntype * ntype == gmx::ssize(*nbp),
2626 "Number of pairs of generated non-bonded parameters should be a perfect square");
2628 /* Add a pair interaction for all non-excluded atom pairs */
2629 const auto& excls = mol->excls;
2630 for (int i = 0; i < n; i++)
2632 for (int j = i + 1; j < n; j++)
2634 bool pairIsExcluded = false;
2635 for (const int atomK : excls[i])
2639 pairIsExcluded = true;
2642 if (!pairIsExcluded)
2644 if (nb_funct != F_LJ)
2646 auto message = gmx::formatString(
2647 "Can only generate non-bonded pair interactions "
2648 "for Van der Waals type Lennard-Jones");
2649 warning_error_and_exit(wi, message, FARGS);
2651 std::vector<int> atoms = { i, j };
2652 std::vector<real> forceParam = {
2653 atom[i].q, atom[j].q, nbp->interactionTypes[ntype * atom[i].type + atom[j].type].c0(),
2654 nbp->interactionTypes[ntype * atom[i].type + atom[j].type].c1()
2656 add_param_to_list(&mol->interactions[F_LJC_PAIRS_NB], InteractionOfType(atoms, forceParam));
2662 static void set_excl_all(gmx::ListOfLists<int>* excl)
2664 /* Get rid of the current exclusions and exclude all atom pairs */
2665 const int numAtoms = excl->ssize();
2666 std::vector<int> exclusionsForAtom(numAtoms);
2667 for (int i = 0; i < numAtoms; i++)
2669 exclusionsForAtom[i] = i;
2672 for (int i = 0; i < numAtoms; i++)
2674 excl->pushBack(exclusionsForAtom);
2678 static void decouple_atoms(t_atoms* atoms,
2679 int atomtype_decouple,
2682 const char* mol_name,
2687 for (i = 0; i < atoms->nr; i++)
2691 atom = &atoms->atom[i];
2693 if (atom->qB != atom->q || atom->typeB != atom->type)
2695 auto message = gmx::formatString(
2696 "Atom %d in molecule type '%s' has different A and B state "
2697 "charges and/or atom types set in the topology file as well "
2698 "as through the mdp option '%s'. You can not use both "
2699 "these methods simultaneously.",
2700 i + 1, mol_name, "couple-moltype");
2701 warning_error_and_exit(wi, message, FARGS);
2704 if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamVDW)
2708 if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamQ)
2710 atom->type = atomtype_decouple;
2712 if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamVDW)
2716 if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamQ)
2718 atom->typeB = atomtype_decouple;
2723 void convert_moltype_couple(MoleculeInformation* mol,
2724 int atomtype_decouple,
2730 InteractionsOfType* nbp,
2733 convert_pairs_to_pairsQ(mol->interactions, fudgeQQ, &mol->atoms);
2736 generate_LJCpairsNB(mol, nb_funct, nbp, wi);
2737 set_excl_all(&mol->excls);
2739 decouple_atoms(&mol->atoms, atomtype_decouple, couple_lam0, couple_lam1, *mol->name, wi);