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, 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/notset.h"
53 #include "gromacs/gmxpreprocess/readir.h"
54 #include "gromacs/gmxpreprocess/topdirs.h"
55 #include "gromacs/gmxpreprocess/toputil.h"
56 #include "gromacs/mdtypes/md_enums.h"
57 #include "gromacs/topology/ifunc.h"
58 #include "gromacs/topology/symtab.h"
59 #include "gromacs/utility/arrayref.h"
60 #include "gromacs/utility/cstringutil.h"
61 #include "gromacs/utility/fatalerror.h"
62 #include "gromacs/utility/gmxassert.h"
63 #include "gromacs/utility/smalloc.h"
65 void generate_nbparams(int comb, int ftype, t_params *plist, gpp_atomtype_t atype,
70 real c, bi, bj, ci, cj, ci0, ci1, ci2, cj0, cj1, cj2;
73 /* Lean mean shortcuts */
74 nr = get_atomtype_ntypes(atype);
76 snew(plist->param, nr*nr);
79 /* Fill the matrix with force parameters */
87 for (i = k = 0; (i < nr); i++)
89 for (j = 0; (j < nr); j++, k++)
91 for (nf = 0; (nf < nrfp); nf++)
93 ci = get_atomtype_nbparam(i, nf, atype);
94 cj = get_atomtype_nbparam(j, nf, atype);
95 c = std::sqrt(ci * cj);
96 plist->param[k].c[nf] = c;
102 case eCOMB_ARITHMETIC:
103 /* c0 and c1 are sigma and epsilon */
104 for (i = k = 0; (i < nr); i++)
106 for (j = 0; (j < nr); j++, k++)
108 ci0 = get_atomtype_nbparam(i, 0, atype);
109 cj0 = get_atomtype_nbparam(j, 0, atype);
110 ci1 = get_atomtype_nbparam(i, 1, atype);
111 cj1 = get_atomtype_nbparam(j, 1, atype);
112 plist->param[k].c[0] = (fabs(ci0) + fabs(cj0))*0.5;
113 /* Negative sigma signals that c6 should be set to zero later,
114 * so we need to propagate that through the combination rules.
116 if (ci0 < 0 || cj0 < 0)
118 plist->param[k].c[0] *= -1;
120 plist->param[k].c[1] = std::sqrt(ci1*cj1);
125 case eCOMB_GEOM_SIG_EPS:
126 /* c0 and c1 are sigma and epsilon */
127 for (i = k = 0; (i < nr); i++)
129 for (j = 0; (j < nr); j++, k++)
131 ci0 = get_atomtype_nbparam(i, 0, atype);
132 cj0 = get_atomtype_nbparam(j, 0, atype);
133 ci1 = get_atomtype_nbparam(i, 1, atype);
134 cj1 = get_atomtype_nbparam(j, 1, atype);
135 plist->param[k].c[0] = std::sqrt(fabs(ci0*cj0));
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)
141 plist->param[k].c[0] *= -1;
143 plist->param[k].c[1] = std::sqrt(ci1*cj1);
149 sprintf(errbuf, "No such combination rule %d", comb);
150 warning_error_and_exit(wi, errbuf, FARGS);
154 gmx_incons("Topology processing, generate nb parameters");
159 /* Buckingham rules */
160 for (i = k = 0; (i < nr); i++)
162 for (j = 0; (j < nr); j++, k++)
164 ci0 = get_atomtype_nbparam(i, 0, atype);
165 cj0 = get_atomtype_nbparam(j, 0, atype);
166 ci2 = get_atomtype_nbparam(i, 2, atype);
167 cj2 = get_atomtype_nbparam(j, 2, atype);
168 bi = get_atomtype_nbparam(i, 1, atype);
169 bj = get_atomtype_nbparam(j, 1, atype);
170 plist->param[k].c[0] = std::sqrt(ci0 * cj0);
171 if ((bi == 0) || (bj == 0))
173 plist->param[k].c[1] = 0;
177 plist->param[k].c[1] = 2.0/(1/bi+1/bj);
179 plist->param[k].c[2] = std::sqrt(ci2 * cj2);
185 sprintf(errbuf, "Invalid nonbonded type %s",
186 interaction_function[ftype].longname);
187 warning_error(wi, errbuf);
191 static void realloc_nb_params(gpp_atomtype_t at,
192 t_nbparam ***nbparam, t_nbparam ***pair)
194 /* Add space in the non-bonded parameters matrix */
195 int atnr = get_atomtype_ntypes(at);
196 srenew(*nbparam, atnr);
197 snew((*nbparam)[atnr-1], atnr);
201 snew((*pair)[atnr-1], atnr);
205 static void copy_B_from_A(int ftype, double *c)
209 nrfpA = NRFPA(ftype);
210 nrfpB = NRFPB(ftype);
212 /* Copy the B parameters from the first nrfpB A parameters */
213 for (i = 0; (i < nrfpB); i++)
219 void push_at (t_symtab *symtab, gpp_atomtype_t at, t_bond_atomtype bat,
220 char *line, int nb_funct,
221 t_nbparam ***nbparam, t_nbparam ***pair,
228 t_xlate xl[eptNR] = {
236 int nr, i, nfields, j, pt, nfp0 = -1;
237 int batype_nr, nread;
238 char type[STRLEN], btype[STRLEN], ptype[STRLEN];
240 double c[MAXFORCEPARAM];
241 double radius, vol, surftens, gb_radius, S_hct;
242 char tmpfield[12][100]; /* Max 12 fields of width 100 */
247 gmx_bool have_atomic_number;
248 gmx_bool have_bonded_type;
253 /* First assign input line to temporary array */
254 nfields = sscanf(line, "%s%s%s%s%s%s%s%s%s%s%s%s",
255 tmpfield[0], tmpfield[1], tmpfield[2], tmpfield[3], tmpfield[4], tmpfield[5],
256 tmpfield[6], tmpfield[7], tmpfield[8], tmpfield[9], tmpfield[10], tmpfield[11]);
258 /* Comments on optional fields in the atomtypes section:
260 * The force field format is getting a bit old. For OPLS-AA we needed
261 * to add a special bonded atomtype, and for Gerrit Groenhofs QM/MM stuff
262 * we also needed the atomic numbers.
263 * To avoid making all old or user-generated force fields unusable we
264 * have introduced both these quantities as optional columns, and do some
265 * acrobatics to check whether they are present or not.
266 * This will all look much nicer when we switch to XML... sigh.
268 * Field 0 (mandatory) is the nonbonded type name. (string)
269 * Field 1 (optional) is the bonded type (string)
270 * Field 2 (optional) is the atomic number (int)
271 * Field 3 (mandatory) is the mass (numerical)
272 * Field 4 (mandatory) is the charge (numerical)
273 * Field 5 (mandatory) is the particle type (single character)
274 * This is followed by a number of nonbonded parameters.
276 * The safest way to identify the format is the particle type field.
278 * So, here is what we do:
280 * A. Read in the first six fields as strings
281 * B. If field 3 (starting from 0) is a single char, we have neither
282 * bonded_type or atomic numbers.
283 * C. If field 5 is a single char we have both.
284 * D. If field 4 is a single char we check field 1. If this begins with
285 * an alphabetical character we have bonded types, otherwise atomic numbers.
294 if ( (strlen(tmpfield[5]) == 1) && isalpha(tmpfield[5][0]) )
296 have_bonded_type = TRUE;
297 have_atomic_number = TRUE;
299 else if ( (strlen(tmpfield[3]) == 1) && isalpha(tmpfield[3][0]) )
301 have_bonded_type = FALSE;
302 have_atomic_number = FALSE;
306 have_bonded_type = ( isalpha(tmpfield[1][0]) != 0 );
307 have_atomic_number = !have_bonded_type;
310 /* optional fields */
324 if (have_atomic_number)
326 if (have_bonded_type)
328 nread = sscanf(line, "%s%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf",
329 type, btype, &atomnr, &m, &q, ptype, &c[0], &c[1],
330 &radius, &vol, &surftens, &gb_radius);
339 /* have_atomic_number && !have_bonded_type */
340 nread = sscanf(line, "%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf",
341 type, &atomnr, &m, &q, ptype, &c[0], &c[1],
342 &radius, &vol, &surftens, &gb_radius);
352 if (have_bonded_type)
354 /* !have_atomic_number && have_bonded_type */
355 nread = sscanf(line, "%s%s%lf%lf%s%lf%lf%lf%lf%lf%lf",
356 type, btype, &m, &q, ptype, &c[0], &c[1],
357 &radius, &vol, &surftens, &gb_radius);
366 /* !have_atomic_number && !have_bonded_type */
367 nread = sscanf(line, "%s%lf%lf%s%lf%lf%lf%lf%lf%lf",
368 type, &m, &q, ptype, &c[0], &c[1],
369 &radius, &vol, &surftens, &gb_radius);
378 if (!have_bonded_type)
383 if (!have_atomic_number)
393 if (have_atomic_number)
395 if (have_bonded_type)
397 nread = sscanf(line, "%s%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
398 type, btype, &atomnr, &m, &q, ptype, &c[0], &c[1], &c[2],
399 &radius, &vol, &surftens, &gb_radius);
408 /* have_atomic_number && !have_bonded_type */
409 nread = sscanf(line, "%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
410 type, &atomnr, &m, &q, ptype, &c[0], &c[1], &c[2],
411 &radius, &vol, &surftens, &gb_radius);
421 if (have_bonded_type)
423 /* !have_atomic_number && have_bonded_type */
424 nread = sscanf(line, "%s%s%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
425 type, btype, &m, &q, ptype, &c[0], &c[1], &c[2],
426 &radius, &vol, &surftens, &gb_radius);
435 /* !have_atomic_number && !have_bonded_type */
436 nread = sscanf(line, "%s%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
437 type, &m, &q, ptype, &c[0], &c[1], &c[2],
438 &radius, &vol, &surftens, &gb_radius);
447 if (!have_bonded_type)
452 if (!have_atomic_number)
460 sprintf(errbuf, "Invalid function type %d in push_at", nb_funct);
461 warning_error_and_exit(wi, errbuf, FARGS);
463 for (j = nfp0; (j < MAXFORCEPARAM); j++)
468 if (strlen(type) == 1 && isdigit(type[0]))
470 warning_error_and_exit(wi, "Atom type names can't be single digits.", FARGS);
473 if (strlen(btype) == 1 && isdigit(btype[0]))
475 warning_error_and_exit(wi, "Bond atom type names can't be single digits.", FARGS);
478 /* Hack to read old topologies */
479 if (gmx_strcasecmp(ptype, "D") == 0)
483 for (j = 0; (j < eptNR); j++)
485 if (gmx_strcasecmp(ptype, xl[j].entry) == 0)
492 sprintf(errbuf, "Invalid particle type %s", ptype);
493 warning_error_and_exit(wi, errbuf, FARGS);
498 fprintf(debug, "ptype: %s\n", ptype_str[pt]);
504 for (i = 0; (i < MAXFORCEPARAM); i++)
509 if ((batype_nr = get_bond_atomtype_type(btype, bat)) == NOTSET)
511 add_bond_atomtype(bat, symtab, btype);
513 batype_nr = get_bond_atomtype_type(btype, bat);
515 if ((nr = get_atomtype_type(type, at)) != NOTSET)
517 sprintf(errbuf, "Overriding atomtype %s", type);
519 if ((nr = set_atomtype(nr, at, symtab, atom, type, param, batype_nr,
520 radius, vol, surftens, atomnr, gb_radius, S_hct)) == NOTSET)
522 sprintf(errbuf, "Replacing atomtype %s failed", type);
523 warning_error_and_exit(wi, errbuf, FARGS);
526 else if ((add_atomtype(at, symtab, atom, type, param,
527 batype_nr, radius, vol,
528 surftens, atomnr, gb_radius, S_hct)) == NOTSET)
530 sprintf(errbuf, "Adding atomtype %s failed", type);
531 warning_error_and_exit(wi, errbuf, FARGS);
535 /* Add space in the non-bonded parameters matrix */
536 realloc_nb_params(at, nbparam, pair);
542 //! Return whether the contents of \c a and \c b are the same, considering also reversed order.
543 template <typename T>
544 static bool equalEitherForwardOrBackward(gmx::ConstArrayRef<T> a, gmx::ConstArrayRef<T> b)
546 return (std::equal(a.begin(), a.end(), b.begin()) ||
547 std::equal(a.begin(), a.end(), b.rbegin()));
550 static void push_bondtype(t_params * bt,
554 gmx_bool bAllowRepeat,
559 int nrfp = NRFP(ftype);
562 /* If bAllowRepeat is TRUE, we allow multiple entries as long as they
563 are on directly _adjacent_ lines.
566 /* First check if our atomtypes are _identical_ (not reversed) to the previous
567 entry. If they are not identical we search for earlier duplicates. If they are
568 we can skip it, since we already searched for the first line
572 bool isContinuationOfBlock = false;
573 if (bAllowRepeat && nr > 1)
575 isContinuationOfBlock = true;
576 for (int j = 0; j < nral; j++)
578 if (b->a[j] != bt->param[nr - 2].a[j])
580 isContinuationOfBlock = false;
585 /* Search for earlier duplicates if this entry was not a continuation
586 from the previous line.
588 bool addBondType = true;
589 bool haveWarned = false;
590 bool haveErrored = false;
591 for (int i = 0; (i < nr); i++)
593 gmx::ConstArrayRef<int> bParams(b->a, b->a + nral);
594 gmx::ConstArrayRef<int> testParams(bt->param[i].a, bt->param[i].a + nral);
595 if (equalEitherForwardOrBackward(bParams, testParams))
597 GMX_ASSERT(nrfp <= MAXFORCEPARAM, "This is ensured in other places, but we need this assert to keep the clang analyzer happy");
598 // TODO consider improving the following code by using:
599 // bool identicalParameters = std::equal(bt->param[i].c, bt->param[i].c + nrfp, b->c);
600 bool identicalParameters = true;
601 for (int j = 0; (j < nrfp); j++)
603 identicalParameters = identicalParameters && (bt->param[i].c[j] == b->c[j]);
606 if (!bAllowRepeat || identicalParameters)
611 if (!identicalParameters)
615 /* With dihedral type 9 we only allow for repeating
616 * of the same parameters with blocks with 1 entry.
617 * Allowing overriding is too complex to check.
619 if (!isContinuationOfBlock && !haveErrored)
621 warning_error(wi, "Encountered a second block of parameters for dihedral type 9 for the same atoms, with either different parameters and/or the first block has multiple lines. This is not supported.");
625 else if (!haveWarned)
627 sprintf(errbuf, "Overriding %s parameters.%s",
628 interaction_function[ftype].longname,
630 "\nUse dihedraltype 9 to allow several multiplicity terms. Only consecutive lines are combined. Non-consective lines overwrite each other."
634 fprintf(stderr, " old: ");
635 for (int j = 0; (j < nrfp); j++)
637 fprintf(stderr, " %g", bt->param[i].c[j]);
639 fprintf(stderr, " \n new: %s\n\n", line);
645 if (!identicalParameters && !bAllowRepeat)
647 /* Overwrite the parameters with the latest ones */
648 // TODO considering improving the following code by replacing with:
649 // std::copy(b->c, b->c + nrfp, bt->param[i].c);
650 for (int j = 0; (j < nrfp); j++)
652 bt->param[i].c[j] = b->c[j];
663 /* fill the arrays up and down */
664 memcpy(bt->param[bt->nr].c, b->c, sizeof(b->c));
665 memcpy(bt->param[bt->nr].a, b->a, sizeof(b->a));
666 memcpy(bt->param[bt->nr+1].c, b->c, sizeof(b->c));
668 /* The definitions of linear angles depend on the order of atoms,
669 * that means that for atoms i-j-k, with certain parameter a, the
670 * corresponding k-j-i angle will have parameter 1-a.
672 if (ftype == F_LINEAR_ANGLES)
674 bt->param[bt->nr+1].c[0] = 1-bt->param[bt->nr+1].c[0];
675 bt->param[bt->nr+1].c[2] = 1-bt->param[bt->nr+1].c[2];
678 for (int j = 0; (j < nral); j++)
680 bt->param[bt->nr+1].a[j] = b->a[nral-1-j];
687 void push_bt(directive d, t_params bt[], int nral,
689 t_bond_atomtype bat, char *line,
692 const char *formal[MAXATOMLIST+1] = {
701 const char *formnl[MAXATOMLIST+1] = {
707 "%*s%*s%*s%*s%*s%*s",
708 "%*s%*s%*s%*s%*s%*s%*s"
710 const char *formlf = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
711 int i, ft, ftype, nn, nrfp, nrfpA;
713 char alc[MAXATOMLIST+1][20];
714 /* One force parameter more, so we can check if we read too many */
715 double c[MAXFORCEPARAM+1];
719 if ((bat && at) || (!bat && !at))
721 gmx_incons("You should pass either bat or at to push_bt");
724 /* Make format string (nral ints+functype) */
725 if ((nn = sscanf(line, formal[nral],
726 alc[0], alc[1], alc[2], alc[3], alc[4], alc[5])) != nral+1)
728 sprintf(errbuf, "Not enough atomtypes (%d instead of %d)", nn-1, nral);
729 warning_error(wi, errbuf);
733 ft = strtol(alc[nral], nullptr, 10);
734 ftype = ifunc_index(d, ft);
736 nrfpA = interaction_function[ftype].nrfpA;
737 strcpy(f1, formnl[nral]);
739 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], &c[10], &c[11], &c[12]))
744 /* Copy the B-state from the A-state */
745 copy_B_from_A(ftype, c);
751 warning_error(wi, "Not enough parameters");
753 else if (nn > nrfpA && nn < nrfp)
755 warning_error(wi, "Too many parameters or not enough parameters for topology B");
759 warning_error(wi, "Too many parameters");
761 for (i = nn; (i < nrfp); i++)
767 for (i = 0; (i < nral); i++)
769 if (at && ((p.a[i] = get_atomtype_type(alc[i], at)) == NOTSET))
771 sprintf(errbuf, "Unknown atomtype %s\n", alc[i]);
772 warning_error_and_exit(wi, errbuf, FARGS);
774 else if (bat && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET))
776 sprintf(errbuf, "Unknown bond_atomtype %s\n", alc[i]);
777 warning_error_and_exit(wi, errbuf, FARGS);
780 for (i = 0; (i < nrfp); i++)
784 push_bondtype (&(bt[ftype]), &p, nral, ftype, FALSE, line, wi);
788 void push_dihedraltype(directive d, t_params bt[],
789 t_bond_atomtype bat, char *line,
792 const char *formal[MAXATOMLIST+1] = {
801 const char *formnl[MAXATOMLIST+1] = {
807 "%*s%*s%*s%*s%*s%*s",
808 "%*s%*s%*s%*s%*s%*s%*s"
810 const char *formlf[MAXFORCEPARAM] = {
816 "%lf%lf%lf%lf%lf%lf",
817 "%lf%lf%lf%lf%lf%lf%lf",
818 "%lf%lf%lf%lf%lf%lf%lf%lf",
819 "%lf%lf%lf%lf%lf%lf%lf%lf%lf",
820 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
821 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
822 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
824 int i, ft, ftype, nn, nrfp, nrfpA, nral;
826 char alc[MAXATOMLIST+1][20];
827 double c[MAXFORCEPARAM];
829 gmx_bool bAllowRepeat;
832 /* This routine accepts dihedraltypes defined from either 2 or 4 atoms.
834 * We first check for 2 atoms with the 3th column being an integer
835 * defining the type. If this isn't the case, we try it with 4 atoms
836 * and the 5th column defining the dihedral type.
838 nn = sscanf(line, formal[4], alc[0], alc[1], alc[2], alc[3], alc[4]);
839 if (nn >= 3 && strlen(alc[2]) == 1 && isdigit(alc[2][0]))
842 ft = strtol(alc[nral], nullptr, 10);
843 /* Move atom types around a bit and use 'X' for wildcard atoms
844 * to create a 4-atom dihedral definition with arbitrary atoms in
847 if (alc[2][0] == '2')
849 /* improper - the two atomtypes are 1,4. Use wildcards for 2,3 */
850 strcpy(alc[3], alc[1]);
851 sprintf(alc[2], "X");
852 sprintf(alc[1], "X");
853 /* alc[0] stays put */
857 /* proper - the two atomtypes are 2,3. Use wildcards for 1,4 */
858 sprintf(alc[3], "X");
859 strcpy(alc[2], alc[1]);
860 strcpy(alc[1], alc[0]);
861 sprintf(alc[0], "X");
864 else if (nn == 5 && strlen(alc[4]) == 1 && isdigit(alc[4][0]))
867 ft = strtol(alc[nral], nullptr, 10);
871 sprintf(errbuf, "Incorrect number of atomtypes for dihedral (%d instead of 2 or 4)", nn-1);
872 warning_error(wi, errbuf);
878 /* Previously, we have always overwritten parameters if e.g. a torsion
879 with the same atomtypes occurs on multiple lines. However, CHARMM and
880 some other force fields specify multiple dihedrals over some bonds,
881 including cosines with multiplicity 6 and somethimes even higher.
882 Thus, they cannot be represented with Ryckaert-Bellemans terms.
883 To add support for these force fields, Dihedral type 9 is identical to
884 normal proper dihedrals, but repeated entries are allowed.
891 bAllowRepeat = FALSE;
895 ftype = ifunc_index(d, ft);
897 nrfpA = interaction_function[ftype].nrfpA;
899 strcpy(f1, formnl[nral]);
900 strcat(f1, formlf[nrfp-1]);
902 /* Check number of parameters given */
903 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], &c[10], &c[11]))
908 /* Copy the B-state from the A-state */
909 copy_B_from_A(ftype, c);
915 warning_error(wi, "Not enough parameters");
917 else if (nn > nrfpA && nn < nrfp)
919 warning_error(wi, "Too many parameters or not enough parameters for topology B");
923 warning_error(wi, "Too many parameters");
925 for (i = nn; (i < nrfp); i++)
932 for (i = 0; (i < 4); i++)
934 if (!strcmp(alc[i], "X"))
940 if ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET)
942 sprintf(errbuf, "Unknown bond_atomtype %s", alc[i]);
943 warning_error_and_exit(wi, errbuf, FARGS);
947 for (i = 0; (i < nrfp); i++)
951 /* Always use 4 atoms here, since we created two wildcard atoms
952 * if there wasn't of them 4 already.
954 push_bondtype (&(bt[ftype]), &p, 4, ftype, bAllowRepeat, line, wi);
958 void push_nbt(directive d, t_nbparam **nbt, gpp_atomtype_t atype,
959 char *pline, int nb_funct,
963 const char *form3 = "%*s%*s%*s%lf%lf%lf";
964 const char *form4 = "%*s%*s%*s%lf%lf%lf%lf";
965 const char *form5 = "%*s%*s%*s%lf%lf%lf%lf%lf";
967 int i, f, n, ftype, nrfp;
975 if (sscanf (pline, "%s%s%d", a0, a1, &f) != 3)
981 ftype = ifunc_index(d, f);
983 if (ftype != nb_funct)
985 sprintf(errbuf, "Trying to add %s while the default nonbond type is %s",
986 interaction_function[ftype].longname,
987 interaction_function[nb_funct].longname);
988 warning_error(wi, errbuf);
992 /* Get the force parameters */
996 n = sscanf(pline, form4, &c[0], &c[1], &c[2], &c[3]);
1002 /* When the B topology parameters are not set,
1003 * copy them from topology A
1005 GMX_ASSERT(nrfp <= 4, "LJ-14 cannot have more than 4 parameters");
1006 for (i = n; i < nrfp; i++)
1011 else if (ftype == F_LJC14_Q)
1013 n = sscanf(pline, form5, &c[0], &c[1], &c[2], &c[3], &dum);
1016 incorrect_n_param(wi);
1022 if (sscanf(pline, form3, &c[0], &c[1], &dum) != 2)
1024 incorrect_n_param(wi);
1030 if (sscanf(pline, form4, &c[0], &c[1], &c[2], &dum) != 3)
1032 incorrect_n_param(wi);
1038 sprintf(errbuf, "Number of force parameters for nonbonded interactions is %d", nrfp);
1039 warning_error_and_exit(wi, errbuf, FARGS);
1041 for (i = 0; (i < nrfp); i++)
1046 /* Put the parameters in the matrix */
1047 if ((ai = get_atomtype_type (a0, atype)) == NOTSET)
1049 sprintf(errbuf, "Atomtype %s not found", a0);
1050 warning_error_and_exit(wi, errbuf, FARGS);
1052 if ((aj = get_atomtype_type (a1, atype)) == NOTSET)
1054 sprintf(errbuf, "Atomtype %s not found", a1);
1055 warning_error_and_exit(wi, errbuf, FARGS);
1057 nbp = &(nbt[std::max(ai, aj)][std::min(ai, aj)]);
1062 for (i = 0; i < nrfp; i++)
1064 bId = bId && (nbp->c[i] == cr[i]);
1068 sprintf(errbuf, "Overriding non-bonded parameters,");
1069 warning(wi, errbuf);
1070 fprintf(stderr, " old:");
1071 for (i = 0; i < nrfp; i++)
1073 fprintf(stderr, " %g", nbp->c[i]);
1075 fprintf(stderr, " new\n%s\n", pline);
1079 for (i = 0; i < nrfp; i++)
1086 push_gb_params (gpp_atomtype_t at, char *line,
1090 double radius, vol, surftens, gb_radius, S_hct;
1091 char atypename[STRLEN];
1092 char errbuf[STRLEN];
1094 if ( (sscanf(line, "%s%lf%lf%lf%lf%lf", atypename, &radius, &vol, &surftens, &gb_radius, &S_hct)) != 6)
1096 sprintf(errbuf, "Too few gb parameters for type %s\n", atypename);
1097 warning(wi, errbuf);
1100 /* Search for atomtype */
1101 atype = get_atomtype_type(atypename, at);
1103 if (atype == NOTSET)
1105 printf("Couldn't find topology match for atomtype %s\n", atypename);
1109 set_atomtype_gbparam(at, atype, radius, vol, surftens, gb_radius, S_hct);
1113 push_cmaptype(directive d, t_params bt[], int nral, gpp_atomtype_t at,
1114 t_bond_atomtype bat, char *line,
1117 const char *formal = "%s%s%s%s%s%s%s%s%n";
1119 int i, ft, ftype, nn, nrfp, nrfpA, nrfpB;
1120 int start, nchar_consumed;
1121 int nxcmap, nycmap, ncmap, read_cmap, sl, nct;
1122 char s[20], alc[MAXATOMLIST+2][20];
1124 char errbuf[STRLEN];
1126 /* Keep the compiler happy */
1130 GMX_ASSERT(nral == 5, "CMAP requires 5 atoms per interaction");
1132 /* Here we can only check for < 8 */
1133 if ((nn = sscanf(line, formal, alc[0], alc[1], alc[2], alc[3], alc[4], alc[5], alc[6], alc[7], &nchar_consumed)) < nral+3)
1135 sprintf(errbuf, "Incorrect number of atomtypes for cmap (%d instead of 5)", nn-1);
1136 warning_error(wi, errbuf);
1139 start += nchar_consumed;
1141 ft = strtol(alc[nral], nullptr, 10);
1142 nxcmap = strtol(alc[nral+1], nullptr, 10);
1143 nycmap = strtol(alc[nral+2], nullptr, 10);
1145 /* Check for equal grid spacing in x and y dims */
1146 if (nxcmap != nycmap)
1148 sprintf(errbuf, "Not the same grid spacing in x and y for cmap grid: x=%d, y=%d", nxcmap, nycmap);
1149 warning_error(wi, errbuf);
1152 ncmap = nxcmap*nycmap;
1153 ftype = ifunc_index(d, ft);
1154 nrfpA = strtol(alc[6], nullptr, 10)*strtol(alc[6], nullptr, 10);
1155 nrfpB = strtol(alc[7], nullptr, 10)*strtol(alc[7], nullptr, 10);
1158 /* Allocate memory for the CMAP grid */
1159 bt[F_CMAP].ncmap += nrfp;
1160 srenew(bt[F_CMAP].cmap, bt[F_CMAP].ncmap);
1162 /* Read in CMAP parameters */
1164 for (i = 0; i < ncmap; i++)
1166 while (isspace(*(line+start+sl)))
1170 nn = sscanf(line+start+sl, " %s ", s);
1172 bt[F_CMAP].cmap[i+(bt[F_CMAP].ncmap)-nrfp] = strtod(s, nullptr);
1180 sprintf(errbuf, "Error in reading cmap parameter for angle %s %s %s %s %s", alc[0], alc[1], alc[2], alc[3], alc[4]);
1181 warning_error(wi, errbuf);
1186 /* Check do that we got the number of parameters we expected */
1187 if (read_cmap == nrfpA)
1189 for (i = 0; i < ncmap; i++)
1191 bt[F_CMAP].cmap[i+ncmap] = bt[F_CMAP].cmap[i];
1196 if (read_cmap < nrfpA)
1198 warning_error(wi, "Not enough cmap parameters");
1200 else if (read_cmap > nrfpA && read_cmap < nrfp)
1202 warning_error(wi, "Too many cmap parameters or not enough parameters for topology B");
1204 else if (read_cmap > nrfp)
1206 warning_error(wi, "Too many cmap parameters");
1211 /* Set grid spacing and the number of grids (we assume these numbers to be the same for all grids
1212 * so we can safely assign them each time
1214 bt[F_CMAP].grid_spacing = nxcmap; /* Or nycmap, they need to be equal */
1215 bt[F_CMAP].nc = bt[F_CMAP].nc + 1; /* Since we are incrementing here, we need to subtract later, see (*****) */
1216 nct = (nral+1) * bt[F_CMAP].nc;
1218 /* Allocate memory for the cmap_types information */
1219 srenew(bt[F_CMAP].cmap_types, nct);
1221 for (i = 0; (i < nral); i++)
1223 if (at && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET))
1225 sprintf(errbuf, "Unknown atomtype %s\n", alc[i]);
1226 warning_error(wi, errbuf);
1228 else if (bat && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET))
1230 sprintf(errbuf, "Unknown bond_atomtype %s\n", alc[i]);
1231 warning_error(wi, errbuf);
1234 /* Assign a grid number to each cmap_type */
1235 bt[F_CMAP].cmap_types[bt[F_CMAP].nct++] = get_bond_atomtype_type(alc[i], bat);
1238 /* Assign a type number to this cmap */
1239 bt[F_CMAP].cmap_types[bt[F_CMAP].nct++] = bt[F_CMAP].nc-1; /* Since we inremented earlier, we need to subtrac here, to get the types right (****) */
1241 /* Check for the correct number of atoms (again) */
1242 if (bt[F_CMAP].nct != nct)
1244 sprintf(errbuf, "Incorrect number of atom types (%d) in cmap type %d\n", nct, bt[F_CMAP].nc);
1245 warning_error(wi, errbuf);
1248 /* Is this correct?? */
1249 for (i = 0; i < MAXFORCEPARAM; i++)
1254 /* Push the bond to the bondlist */
1255 push_bondtype (&(bt[ftype]), &p, nral, ftype, FALSE, line, wi);
1259 static void push_atom_now(t_symtab *symtab, t_atoms *at, int atomnr,
1261 int type, char *ctype, int ptype,
1263 char *resname, char *name, real m0, real q0,
1264 int typeB, char *ctypeB, real mB, real qB,
1267 int j, resind = 0, resnr;
1270 char errbuf[STRLEN];
1272 if (((nr == 0) && (atomnr != 1)) || (nr && (atomnr != at->nr+1)))
1274 sprintf(errbuf, "Atoms in the .top are not numbered consecutively from 1 (rather, atomnr = %d, while at->nr = %d)", atomnr, at->nr);
1275 warning_error_and_exit(wi, errbuf, FARGS);
1278 j = strlen(resnumberic) - 1;
1279 if (isdigit(resnumberic[j]))
1285 ric = resnumberic[j];
1286 if (j == 0 || !isdigit(resnumberic[j-1]))
1288 sprintf(errbuf, "Invalid residue number '%s' for atom %d",
1289 resnumberic, atomnr);
1290 warning_error_and_exit(wi, errbuf, FARGS);
1293 resnr = strtol(resnumberic, nullptr, 10);
1297 resind = at->atom[nr-1].resind;
1299 if (nr == 0 || strcmp(resname, *at->resinfo[resind].name) != 0 ||
1300 resnr != at->resinfo[resind].nr ||
1301 ric != at->resinfo[resind].ic)
1311 at->nres = resind + 1;
1312 srenew(at->resinfo, at->nres);
1313 at->resinfo[resind].name = put_symtab(symtab, resname);
1314 at->resinfo[resind].nr = resnr;
1315 at->resinfo[resind].ic = ric;
1319 resind = at->atom[at->nr-1].resind;
1322 /* New atom instance
1323 * get new space for arrays
1325 srenew(at->atom, nr+1);
1326 srenew(at->atomname, nr+1);
1327 srenew(at->atomtype, nr+1);
1328 srenew(at->atomtypeB, nr+1);
1331 at->atom[nr].type = type;
1332 at->atom[nr].ptype = ptype;
1333 at->atom[nr].q = q0;
1334 at->atom[nr].m = m0;
1335 at->atom[nr].typeB = typeB;
1336 at->atom[nr].qB = qB;
1337 at->atom[nr].mB = mB;
1339 at->atom[nr].resind = resind;
1340 at->atom[nr].atomnumber = atomicnumber;
1341 at->atomname[nr] = put_symtab(symtab, name);
1342 at->atomtype[nr] = put_symtab(symtab, ctype);
1343 at->atomtypeB[nr] = put_symtab(symtab, ctypeB);
1347 void push_cg(t_block *block, int *lastindex, int index, int a)
1351 fprintf (debug, "Index %d, Atom %d\n", index, a);
1354 if (((block->nr) && (*lastindex != index)) || (!block->nr))
1356 /* add a new block */
1358 srenew(block->index, block->nr+1);
1360 block->index[block->nr] = a + 1;
1364 void push_atom(t_symtab *symtab, t_block *cgs,
1365 t_atoms *at, gpp_atomtype_t atype, char *line, int *lastcg,
1369 int cgnumber, atomnr, type, typeB, nscan;
1370 char id[STRLEN], ctype[STRLEN], ctypeB[STRLEN],
1371 resnumberic[STRLEN], resname[STRLEN], name[STRLEN], check[STRLEN];
1372 double m, q, mb, qb;
1373 real m0, q0, mB, qB;
1374 char errbuf[STRLEN];
1376 /* Make a shortcut for writing in this molecule */
1379 /* Fixed parameters */
1380 if (sscanf(line, "%s%s%s%s%s%d",
1381 id, ctype, resnumberic, resname, name, &cgnumber) != 6)
1386 sscanf(id, "%d", &atomnr);
1387 if ((type = get_atomtype_type(ctype, atype)) == NOTSET)
1389 sprintf(errbuf, "Atomtype %s not found", ctype);
1390 warning_error_and_exit(wi, errbuf, FARGS);
1392 ptype = get_atomtype_ptype(type, atype);
1394 /* Set default from type */
1395 q0 = get_atomtype_qA(type, atype);
1396 m0 = get_atomtype_massA(type, atype);
1401 /* Optional parameters */
1402 nscan = sscanf(line, "%*s%*s%*s%*s%*s%*s%lf%lf%s%lf%lf%s",
1403 &q, &m, ctypeB, &qb, &mb, check);
1405 /* Nasty switch that falls thru all the way down! */
1414 if ((typeB = get_atomtype_type(ctypeB, atype)) == NOTSET)
1416 sprintf(errbuf, "Atomtype %s not found", ctypeB);
1417 warning_error_and_exit(wi, errbuf, FARGS);
1419 qB = get_atomtype_qA(typeB, atype);
1420 mB = get_atomtype_massA(typeB, atype);
1429 warning_error(wi, "Too many parameters");
1438 fprintf(debug, "mB=%g, qB=%g, typeB=%d\n", mB, qB, typeB);
1441 push_cg(cgs, lastcg, cgnumber, nr);
1443 push_atom_now(symtab, at, atomnr, get_atomtype_atomnumber(type, atype),
1444 type, ctype, ptype, resnumberic,
1445 resname, name, m0, q0, typeB,
1446 typeB == type ? ctype : ctypeB, mB, qB, wi);
1449 void push_molt(t_symtab *symtab, int *nmol, t_molinfo **mol, char *line,
1455 char errbuf[STRLEN];
1457 if ((sscanf(line, "%s%d", type, &nrexcl)) != 2)
1459 warning_error(wi, "Expected a molecule type name and nrexcl");
1462 /* Test if this moleculetype overwrites another */
1466 if (strcmp(*((*mol)[i].name), type) == 0)
1468 sprintf(errbuf, "moleculetype %s is redefined", type);
1469 warning_error_and_exit(wi, errbuf, FARGS);
1475 srenew(*mol, *nmol);
1476 newmol = &((*mol)[*nmol-1]);
1477 init_molinfo(newmol);
1479 /* Fill in the values */
1480 newmol->name = put_symtab(symtab, type);
1481 newmol->nrexcl = nrexcl;
1482 newmol->excl_set = FALSE;
1485 static gmx_bool default_nb_params(int ftype, t_params bt[], t_atoms *at,
1486 t_param *p, int c_start, gmx_bool bB, gmx_bool bGenPairs)
1488 int i, j, ti, tj, ntype;
1490 t_param *pi = nullptr;
1491 int nr = bt[ftype].nr;
1492 int nral = NRAL(ftype);
1493 int nrfp = interaction_function[ftype].nrfpA;
1494 int nrfpB = interaction_function[ftype].nrfpB;
1496 if ((!bB && nrfp == 0) || (bB && nrfpB == 0))
1504 /* First test the generated-pair position to save
1505 * time when we have 1000*1000 entries for e.g. OPLS...
1507 ntype = static_cast<int>(std::sqrt(static_cast<double>(nr)));
1508 GMX_ASSERT(ntype * ntype == nr, "Number of pairs of generated non-bonded parameters should be a perfect square");
1511 ti = at->atom[p->a[0]].typeB;
1512 tj = at->atom[p->a[1]].typeB;
1516 ti = at->atom[p->a[0]].type;
1517 tj = at->atom[p->a[1]].type;
1519 pi = &(bt[ftype].param[ntype*ti+tj]);
1520 bFound = ((ti == pi->a[0]) && (tj == pi->a[1]));
1523 /* Search explicitly if we didnt find it */
1526 for (i = 0; ((i < nr) && !bFound); i++)
1528 pi = &(bt[ftype].param[i]);
1531 for (j = 0; ((j < nral) &&
1532 (at->atom[p->a[j]].typeB == pi->a[j])); j++)
1539 for (j = 0; ((j < nral) &&
1540 (at->atom[p->a[j]].type == pi->a[j])); j++)
1545 bFound = (j == nral);
1553 if (nrfp+nrfpB > MAXFORCEPARAM)
1555 gmx_incons("Too many force parameters");
1557 for (j = c_start; (j < nrfpB); j++)
1559 p->c[nrfp+j] = pi->c[j];
1564 for (j = c_start; (j < nrfp); j++)
1572 for (j = c_start; (j < nrfp); j++)
1580 static gmx_bool default_cmap_params(t_params bondtype[],
1581 t_atoms *at, gpp_atomtype_t atype,
1582 t_param *p, gmx_bool bB,
1583 int *cmap_type, int *nparam_def,
1586 int i, nparam_found;
1588 gmx_bool bFound = FALSE;
1589 char errbuf[STRLEN];
1594 /* Match the current cmap angle against the list of cmap_types */
1595 for (i = 0; i < bondtype[F_CMAP].nct && !bFound; i += 6)
1604 (get_atomtype_batype(at->atom[p->a[0]].type, atype) == bondtype[F_CMAP].cmap_types[i]) &&
1605 (get_atomtype_batype(at->atom[p->a[1]].type, atype) == bondtype[F_CMAP].cmap_types[i+1]) &&
1606 (get_atomtype_batype(at->atom[p->a[2]].type, atype) == bondtype[F_CMAP].cmap_types[i+2]) &&
1607 (get_atomtype_batype(at->atom[p->a[3]].type, atype) == bondtype[F_CMAP].cmap_types[i+3]) &&
1608 (get_atomtype_batype(at->atom[p->a[4]].type, atype) == bondtype[F_CMAP].cmap_types[i+4]))
1610 /* Found cmap torsion */
1612 ct = bondtype[F_CMAP].cmap_types[i+5];
1618 /* If we did not find a matching type for this cmap torsion */
1621 sprintf(errbuf, "Unknown cmap torsion between atoms %d %d %d %d %d",
1622 p->a[0]+1, p->a[1]+1, p->a[2]+1, p->a[3]+1, p->a[4]+1);
1623 warning_error_and_exit(wi, errbuf, FARGS);
1626 *nparam_def = nparam_found;
1632 /* Returns the number of exact atom type matches, i.e. non wild-card matches,
1633 * returns -1 when there are no matches at all.
1635 static int natom_match(t_param *pi,
1636 int type_i, int type_j, int type_k, int type_l,
1637 const gpp_atomtype_t atype)
1639 if ((pi->ai() == -1 || get_atomtype_batype(type_i, atype) == pi->ai()) &&
1640 (pi->aj() == -1 || get_atomtype_batype(type_j, atype) == pi->aj()) &&
1641 (pi->ak() == -1 || get_atomtype_batype(type_k, atype) == pi->ak()) &&
1642 (pi->al() == -1 || get_atomtype_batype(type_l, atype) == pi->al()))
1645 (pi->ai() == -1 ? 0 : 1) +
1646 (pi->aj() == -1 ? 0 : 1) +
1647 (pi->ak() == -1 ? 0 : 1) +
1648 (pi->al() == -1 ? 0 : 1);
1656 static gmx_bool default_params(int ftype, t_params bt[],
1657 t_atoms *at, gpp_atomtype_t atype,
1658 t_param *p, gmx_bool bB,
1659 t_param **param_def,
1663 gmx_bool bFound, bSame;
1664 t_param *pi = nullptr;
1665 t_param *pj = nullptr;
1666 int nr = bt[ftype].nr;
1667 int nral = NRAL(ftype);
1668 int nrfpA = interaction_function[ftype].nrfpA;
1669 int nrfpB = interaction_function[ftype].nrfpB;
1671 if ((!bB && nrfpA == 0) || (bB && nrfpB == 0))
1679 if (ftype == F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS || ftype == F_PIDIHS)
1681 int nmatch_max = -1;
1685 /* For dihedrals we allow wildcards. We choose the first type
1686 * that has the most real matches, i.e. non-wildcard matches.
1688 for (t = 0; ((t < nr) && nmatch_max < 4); t++)
1693 pt = &(bt[ftype].param[t]);
1696 nmatch = natom_match(pt, at->atom[p->ai()].typeB, at->atom[p->aj()].typeB, at->atom[p->ak()].typeB, at->atom[p->al()].typeB, atype);
1700 nmatch = natom_match(pt, at->atom[p->ai()].type, at->atom[p->aj()].type, at->atom[p->ak()].type, at->atom[p->al()].type, atype);
1702 if (nmatch > nmatch_max)
1704 nmatch_max = nmatch;
1714 pi = &(bt[ftype].param[i]);
1717 /* Find additional matches for this dihedral - necessary
1719 * The rule in that case is that additional matches
1720 * HAVE to be on adjacent lines!
1723 /* Continue from current i value */
1724 for (j = i + 2; j < nr && bSame; j += 2)
1726 pj = &(bt[ftype].param[j]);
1727 bSame = (pi->ai() == pj->ai() && pi->aj() == pj->aj() && pi->ak() == pj->ak() && pi->al() == pj->al());
1732 /* nparam_found will be increased as long as the numbers match */
1736 else /* Not a dihedral */
1740 for (i = 0; ((i < nr) && !bFound); i++)
1742 pi = &(bt[ftype].param[i]);
1745 for (j = 0; ((j < nral) &&
1746 (get_atomtype_batype(at->atom[p->a[j]].typeB, atype) == pi->a[j])); j++)
1753 for (j = 0; ((j < nral) &&
1754 (get_atomtype_batype(at->atom[p->a[j]].type, atype) == pi->a[j])); j++)
1759 bFound = (j == nral);
1768 *nparam_def = nparam_found;
1775 void push_bond(directive d, t_params bondtype[], t_params bond[],
1776 t_atoms *at, gpp_atomtype_t atype, char *line,
1777 gmx_bool bBonded, gmx_bool bGenPairs, real fudgeQQ,
1778 gmx_bool bZero, gmx_bool *bWarn_copy_A_B,
1781 const char *aaformat[MAXATOMLIST] = {
1789 const char *asformat[MAXATOMLIST] = {
1794 "%*s%*s%*s%*s%*s%*s",
1795 "%*s%*s%*s%*s%*s%*s%*s"
1797 const char *ccformat = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
1798 int nr, i, j, nral, nral_fmt, nread, ftype;
1799 char format[STRLEN];
1800 /* One force parameter more, so we can check if we read too many */
1801 double cc[MAXFORCEPARAM+1];
1802 int aa[MAXATOMLIST+1];
1803 t_param param, *param_defA, *param_defB;
1804 gmx_bool bFoundA = FALSE, bFoundB = FALSE, bDef, bPert, bSwapParity = FALSE;
1805 int nparam_defA, nparam_defB;
1806 char errbuf[STRLEN];
1808 nparam_defA = nparam_defB = 0;
1810 ftype = ifunc_index(d, 1);
1812 for (j = 0; j < MAXATOMLIST; j++)
1816 bDef = (NRFP(ftype) > 0);
1818 if (ftype == F_SETTLE)
1820 /* SETTLE acts on 3 atoms, but the topology format only specifies
1821 * the first atom (for historical reasons).
1830 nread = sscanf(line, aaformat[nral_fmt-1],
1831 &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
1833 if (ftype == F_SETTLE)
1840 if (nread < nral_fmt)
1845 else if (nread > nral_fmt)
1847 /* this is a hack to allow for virtual sites with swapped parity */
1848 bSwapParity = (aa[nral] < 0);
1851 aa[nral] = -aa[nral];
1853 ftype = ifunc_index(d, aa[nral]);
1862 sprintf(errbuf, "Negative function types only allowed for %s and %s",
1863 interaction_function[F_VSITE3FAD].longname,
1864 interaction_function[F_VSITE3OUT].longname);
1865 warning_error_and_exit(wi, errbuf, FARGS);
1871 /* Check for double atoms and atoms out of bounds */
1872 for (i = 0; (i < nral); i++)
1874 if (aa[i] < 1 || aa[i] > at->nr)
1876 sprintf(errbuf, "Atom index (%d) in %s out of bounds (1-%d).\n"
1877 "This probably means that you have inserted topology section \"%s\"\n"
1878 "in a part belonging to a different molecule than you intended to.\n"
1879 "In that case move the \"%s\" section to the right molecule.",
1880 aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
1881 warning_error_and_exit(wi, errbuf, FARGS);
1883 for (j = i+1; (j < nral); j++)
1885 GMX_ASSERT(j < MAXATOMLIST + 1, "Values from nral=NRAL() will satisfy this, we assert to keep gcc 4 happy");
1888 sprintf(errbuf, "Duplicate atom index (%d) in %s", aa[i], dir2str(d));
1889 if (ftype == F_ANGRES)
1891 /* Since the angle restraints uses 2 pairs of atoms to
1892 * defines an angle between vectors, it can be useful
1893 * to use one atom twice, so we only issue a note here.
1895 warning_note(wi, errbuf);
1899 warning_error(wi, errbuf);
1905 /* default force parameters */
1906 for (j = 0; (j < MAXATOMLIST); j++)
1908 param.a[j] = aa[j]-1;
1910 for (j = 0; (j < MAXFORCEPARAM); j++)
1915 /* Get force params for normal and free energy perturbation
1916 * studies, as determined by types!
1921 bFoundA = default_params(ftype, bondtype, at, atype, ¶m, FALSE, ¶m_defA, &nparam_defA);
1924 /* Copy the A-state and B-state default parameters. */
1925 GMX_ASSERT(NRFPA(ftype)+NRFPB(ftype) <= MAXFORCEPARAM, "Bonded interactions may have at most 12 parameters");
1926 for (j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++)
1928 param.c[j] = param_defA->c[j];
1931 bFoundB = default_params(ftype, bondtype, at, atype, ¶m, TRUE, ¶m_defB, &nparam_defB);
1934 /* Copy only the B-state default parameters */
1935 for (j = NRFPA(ftype); (j < NRFP(ftype)); j++)
1937 param.c[j] = param_defB->c[j];
1941 else if (ftype == F_LJ14)
1943 bFoundA = default_nb_params(ftype, bondtype, at, ¶m, 0, FALSE, bGenPairs);
1944 bFoundB = default_nb_params(ftype, bondtype, at, ¶m, 0, TRUE, bGenPairs);
1946 else if (ftype == F_LJC14_Q)
1948 param.c[0] = fudgeQQ;
1949 /* Fill in the A-state charges as default parameters */
1950 param.c[1] = at->atom[param.a[0]].q;
1951 param.c[2] = at->atom[param.a[1]].q;
1952 /* The default LJ parameters are the standard 1-4 parameters */
1953 bFoundA = default_nb_params(F_LJ14, bondtype, at, ¶m, 3, FALSE, bGenPairs);
1956 else if (ftype == F_LJC_PAIRS_NB)
1958 /* Defaults are not supported here */
1964 gmx_incons("Unknown function type in push_bond");
1967 if (nread > nral_fmt)
1969 /* Manually specified parameters - in this case we discard multiple torsion info! */
1971 strcpy(format, asformat[nral_fmt-1]);
1972 strcat(format, ccformat);
1974 nread = sscanf(line, format, &cc[0], &cc[1], &cc[2], &cc[3], &cc[4], &cc[5],
1975 &cc[6], &cc[7], &cc[8], &cc[9], &cc[10], &cc[11], &cc[12]);
1977 if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0))
1979 /* We only have to issue a warning if these atoms are perturbed! */
1981 for (j = 0; (j < nral); j++)
1983 bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1986 if (bPert && *bWarn_copy_A_B)
1989 "Some parameters for bonded interaction involving perturbed atoms are specified explicitly in state A, but not B - copying A to B");
1990 warning(wi, errbuf);
1991 *bWarn_copy_A_B = FALSE;
1994 /* If only the A parameters were specified, copy them to the B state */
1995 /* The B-state parameters correspond to the first nrfpB
1996 * A-state parameters.
1998 for (j = 0; (j < NRFPB(ftype)); j++)
2000 cc[nread++] = cc[j];
2004 /* If nread was 0 or EOF, no parameters were read => use defaults.
2005 * If nread was nrfpA we copied above so nread=nrfp.
2006 * If nread was nrfp we are cool.
2007 * For F_LJC14_Q we allow supplying fudgeQQ only.
2008 * Anything else is an error!
2010 if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) &&
2011 !(ftype == F_LJC14_Q && nread == 1))
2013 sprintf(errbuf, "Incorrect number of parameters - found %d, expected %d or %d for %s.",
2014 nread, NRFPA(ftype), NRFP(ftype), interaction_function[ftype].longname);
2015 warning_error_and_exit(wi, errbuf, FARGS);
2018 for (j = 0; (j < nread); j++)
2023 /* Check whether we have to use the defaults */
2024 if (nread == NRFP(ftype))
2033 /* nread now holds the number of force parameters read! */
2038 /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */
2039 if (ftype == F_PDIHS)
2041 if ((nparam_defA != nparam_defB) || ((nparam_defA > 1 || nparam_defB > 1) && (param_defA != param_defB)))
2044 "Cannot automatically perturb a torsion with multiple terms to different form.\n"
2045 "Please specify perturbed parameters manually for this torsion in your topology!");
2046 warning_error(wi, errbuf);
2050 if (nread > 0 && nread < NRFPA(ftype))
2052 /* Issue an error, do not use defaults */
2053 sprintf(errbuf, "Not enough parameters, there should be at least %d (or 0 for defaults)", NRFPA(ftype));
2054 warning_error(wi, errbuf);
2057 if (nread == 0 || nread == EOF)
2061 if (interaction_function[ftype].flags & IF_VSITE)
2063 /* set them to NOTSET, will be calculated later */
2064 for (j = 0; (j < MAXFORCEPARAM); j++)
2066 param.c[j] = NOTSET;
2071 param.c1() = -1; /* flag to swap parity of vsite construction */
2078 fprintf(stderr, "NOTE: No default %s types, using zeroes\n",
2079 interaction_function[ftype].longname);
2083 sprintf(errbuf, "No default %s types", interaction_function[ftype].longname);
2084 warning_error(wi, errbuf);
2095 param.c0() = 360 - param.c0();
2098 param.c2() = -param.c2();
2105 /* We only have to issue a warning if these atoms are perturbed! */
2107 for (j = 0; (j < nral); j++)
2109 bPert = bPert || PERTURBED(at->atom[param.a[j]]);
2114 sprintf(errbuf, "No default %s types for perturbed atoms, "
2115 "using normal values", interaction_function[ftype].longname);
2116 warning(wi, errbuf);
2122 if ((ftype == F_PDIHS || ftype == F_ANGRES || ftype == F_ANGRESZ)
2123 && param.c[5] != param.c[2])
2125 sprintf(errbuf, "%s multiplicity can not be perturbed %f!=%f",
2126 interaction_function[ftype].longname,
2127 param.c[2], param.c[5]);
2128 warning_error_and_exit(wi, errbuf, FARGS);
2131 if (IS_TABULATED(ftype) && param.c[2] != param.c[0])
2133 sprintf(errbuf, "%s table number can not be perturbed %d!=%d",
2134 interaction_function[ftype].longname,
2135 (int)(param.c[0]+0.5), (int)(param.c[2]+0.5));
2136 warning_error_and_exit(wi, errbuf, FARGS);
2139 /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
2140 if (ftype == F_RBDIHS)
2143 for (i = 0; i < NRFP(ftype); i++)
2145 if (param.c[i] != 0)
2156 /* Put the values in the appropriate arrays */
2157 add_param_to_list (&bond[ftype], ¶m);
2159 /* Push additional torsions from FF for ftype==9 if we have them.
2160 * We have already checked that the A/B states do not differ in this case,
2161 * so we do not have to double-check that again, or the vsite stuff.
2162 * In addition, those torsions cannot be automatically perturbed.
2164 if (bDef && ftype == F_PDIHS)
2166 for (i = 1; i < nparam_defA; i++)
2168 /* Advance pointer! */
2170 for (j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++)
2172 param.c[j] = param_defA->c[j];
2174 /* And push the next term for this torsion */
2175 add_param_to_list (&bond[ftype], ¶m);
2180 void push_cmap(directive d, t_params bondtype[], t_params bond[],
2181 t_atoms *at, gpp_atomtype_t atype, char *line,
2184 const char *aaformat[MAXATOMLIST+1] =
2195 int i, j, ftype, nral, nread, ncmap_params;
2197 int aa[MAXATOMLIST+1];
2198 char errbuf[STRLEN];
2202 ftype = ifunc_index(d, 1);
2206 nread = sscanf(line, aaformat[nral-1],
2207 &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
2214 else if (nread == nral)
2216 ftype = ifunc_index(d, 1);
2219 /* Check for double atoms and atoms out of bounds */
2220 for (i = 0; i < nral; i++)
2222 if (aa[i] < 1 || aa[i] > at->nr)
2224 sprintf(errbuf, "Atom index (%d) in %s out of bounds (1-%d).\n"
2225 "This probably means that you have inserted topology section \"%s\"\n"
2226 "in a part belonging to a different molecule than you intended to.\n"
2227 "In that case move the \"%s\" section to the right molecule.",
2228 aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
2229 warning_error_and_exit(wi, errbuf, FARGS);
2232 for (j = i+1; (j < nral); j++)
2236 sprintf(errbuf, "Duplicate atom index (%d) in %s", aa[i], dir2str(d));
2237 warning_error(wi, errbuf);
2242 /* default force parameters */
2243 for (j = 0; (j < MAXATOMLIST); j++)
2245 param.a[j] = aa[j]-1;
2247 for (j = 0; (j < MAXFORCEPARAM); j++)
2252 /* Get the cmap type for this cmap angle */
2253 bFound = default_cmap_params(bondtype, at, atype, ¶m, FALSE, &cmap_type, &ncmap_params, wi);
2255 /* We want exactly one parameter (the cmap type in state A (currently no state B) back */
2256 if (bFound && ncmap_params == 1)
2258 /* Put the values in the appropriate arrays */
2259 param.c[0] = cmap_type;
2260 add_param_to_list(&bond[ftype], ¶m);
2264 /* This is essentially the same check as in default_cmap_params() done one more time */
2265 sprintf(errbuf, "Unable to assign a cmap type to torsion %d %d %d %d and %d\n",
2266 param.a[0]+1, param.a[1]+1, param.a[2]+1, param.a[3]+1, param.a[4]+1);
2267 warning_error_and_exit(wi, errbuf, FARGS);
2273 void push_vsitesn(directive d, t_params bond[],
2274 t_atoms *at, char *line,
2278 int type, ftype, j, n, ret, nj, a;
2280 double *weight = nullptr, weight_tot;
2282 char errbuf[STRLEN];
2284 /* default force parameters */
2285 for (j = 0; (j < MAXATOMLIST); j++)
2287 param.a[j] = NOTSET;
2289 for (j = 0; (j < MAXFORCEPARAM); j++)
2295 ret = sscanf(ptr, "%d%n", &a, &n);
2299 sprintf(errbuf, "Expected an atom index in section \"%s\"",
2301 warning_error_and_exit(wi, errbuf, FARGS);
2306 sscanf(ptr, "%d%n", &type, &n);
2308 ftype = ifunc_index(d, type);
2314 ret = sscanf(ptr, "%d%n", &a, &n);
2321 srenew(weight, nj+20);
2330 /* Here we use the A-state mass as a parameter.
2331 * Note that the B-state mass has no influence.
2333 weight[nj] = at->atom[atc[nj]].m;
2337 ret = sscanf(ptr, "%lf%n", &(weight[nj]), &n);
2341 sprintf(errbuf, "No weight or negative weight found for vsiten constructing atom %d (atom index %d)",
2343 warning_error_and_exit(wi, errbuf, FARGS);
2347 sprintf(errbuf, "Unknown vsiten type %d", type);
2348 warning_error_and_exit(wi, errbuf, FARGS);
2350 weight_tot += weight[nj];
2358 sprintf(errbuf, "Expected more than one atom index in section \"%s\"",
2360 warning_error_and_exit(wi, errbuf, FARGS);
2363 if (weight_tot == 0)
2365 warning_error_and_exit(wi, "The total mass of the construting atoms is zero", FARGS);
2368 for (j = 0; j < nj; j++)
2370 param.a[1] = atc[j];
2372 param.c[1] = weight[j]/weight_tot;
2373 /* Put the values in the appropriate arrays */
2374 add_param_to_list (&bond[ftype], ¶m);
2381 void push_mol(int nrmols, t_molinfo mols[], char *pline, int *whichmol,
2386 char errbuf[STRLEN];
2388 if (sscanf(pline, "%s%d", type, nrcopies) != 2)
2394 /* Search moleculename.
2395 * Here we originally only did case insensitive matching. But because
2396 * some PDB files can have many chains and use case to generate more
2397 * chain-identifiers, which in turn end up in our moleculetype name,
2398 * we added support for case-sensitivity.
2404 for (int i = 0; i < nrmols; i++)
2406 if (strcmp(type, *(mols[i].name)) == 0)
2411 if (gmx_strcasecmp(type, *(mols[i].name)) == 0)
2420 // select the case sensitive match
2421 *whichmol = matchcs;
2425 // avoid matching case-insensitive when we have multiple matches
2428 sprintf(errbuf, "For moleculetype '%s' in [ system ] %d case insensitive matches, but %d case sensitive matches were found. Check the case of the characters in the moleculetypes.", type, nrci, nrcs);
2429 warning_error_and_exit(wi, errbuf, FARGS);
2433 // select the unique case insensitive match
2434 *whichmol = matchci;
2438 sprintf(errbuf, "No such moleculetype %s", type);
2439 warning_error_and_exit(wi, errbuf, FARGS);
2444 void init_block2(t_block2 *b2, int natom)
2449 snew(b2->nra, b2->nr);
2450 snew(b2->a, b2->nr);
2451 for (i = 0; (i < b2->nr); i++)
2457 void done_block2(t_block2 *b2)
2463 for (i = 0; (i < b2->nr); i++)
2473 void push_excl(char *line, t_block2 *b2, warninp_t wi)
2477 char base[STRLEN], format[STRLEN];
2478 char errbuf[STRLEN];
2480 if (sscanf(line, "%d", &i) == 0)
2485 if ((1 <= i) && (i <= b2->nr))
2493 fprintf(debug, "Unbound atom %d\n", i-1);
2497 strcpy(base, "%*d");
2500 strcpy(format, base);
2501 strcat(format, "%d");
2502 n = sscanf(line, format, &j);
2505 if ((1 <= j) && (j <= b2->nr))
2508 srenew(b2->a[i], ++(b2->nra[i]));
2509 b2->a[i][b2->nra[i]-1] = j;
2510 /* also add the reverse exclusion! */
2511 srenew(b2->a[j], ++(b2->nra[j]));
2512 b2->a[j][b2->nra[j]-1] = i;
2513 strcat(base, "%*d");
2517 sprintf(errbuf, "Invalid Atomnr j: %d, b2->nr: %d\n", j, b2->nr);
2518 warning_error_and_exit(wi, errbuf, FARGS);
2525 void b_to_b2(t_blocka *b, t_block2 *b2)
2530 for (i = 0; (i < b->nr); i++)
2532 for (j = b->index[i]; (j < b->index[i+1]); j++)
2535 srenew(b2->a[i], ++b2->nra[i]);
2536 b2->a[i][b2->nra[i]-1] = a;
2541 void b2_to_b(t_block2 *b2, t_blocka *b)
2547 for (i = 0; (i < b2->nr); i++)
2550 for (j = 0; (j < b2->nra[i]); j++)
2552 b->a[nra+j] = b2->a[i][j];
2556 /* terminate list */
2560 static int icomp(const void *v1, const void *v2)
2562 return (*((int *) v1))-(*((int *) v2));
2565 void merge_excl(t_blocka *excl, t_block2 *b2, warninp_t wi)
2570 char errbuf[STRLEN];
2576 else if (b2->nr != excl->nr)
2578 sprintf(errbuf, "DEATH HORROR: b2->nr = %d, while excl->nr = %d",
2580 warning_error_and_exit(wi, errbuf, FARGS);
2584 fprintf(debug, "Entering merge_excl\n");
2587 /* First copy all entries from excl to b2 */
2590 /* Count and sort the exclusions */
2592 for (i = 0; (i < b2->nr); i++)
2596 /* remove double entries */
2597 qsort(b2->a[i], (size_t)b2->nra[i], (size_t)sizeof(b2->a[i][0]), icomp);
2599 for (j = 1; (j < b2->nra[i]); j++)
2601 if (b2->a[i][j] != b2->a[i][k-1])
2603 b2->a[i][k] = b2->a[i][j];
2612 srenew(excl->a, excl->nra);
2617 int add_atomtype_decoupled(t_symtab *symtab, gpp_atomtype_t at,
2618 t_nbparam ***nbparam, t_nbparam ***pair)
2624 /* Define an atom type with all parameters set to zero (no interactions) */
2627 /* Type for decoupled atoms could be anything,
2628 * this should be changed automatically later when required.
2630 atom.ptype = eptAtom;
2631 for (i = 0; (i < MAXFORCEPARAM); i++)
2636 nr = add_atomtype(at, symtab, &atom, "decoupled", ¶m, -1, 0.0, 0.0, 0.0, 0, 0, 0);
2638 /* Add space in the non-bonded parameters matrix */
2639 realloc_nb_params(at, nbparam, pair);
2644 static void convert_pairs_to_pairsQ(t_params *plist,
2645 real fudgeQQ, t_atoms *atoms)
2647 t_param *paramp1, *paramp2, *paramnew;
2648 int i, j, p1nr, p2nr, p2newnr;
2650 /* Add the pair list to the pairQ list */
2651 p1nr = plist[F_LJ14].nr;
2652 p2nr = plist[F_LJC14_Q].nr;
2653 p2newnr = p1nr + p2nr;
2654 snew(paramnew, p2newnr);
2656 paramp1 = plist[F_LJ14].param;
2657 paramp2 = plist[F_LJC14_Q].param;
2659 /* Fill in the new F_LJC14_Q array with the old one. NOTE:
2660 it may be possible to just ADD the converted F_LJ14 array
2661 to the old F_LJC14_Q array, but since we have to create
2662 a new sized memory structure, better just to deep copy it all.
2665 for (i = 0; i < p2nr; i++)
2667 /* Copy over parameters */
2668 for (j = 0; j < 5; j++) /* entries are 0-4 for F_LJC14_Q */
2670 paramnew[i].c[j] = paramp2[i].c[j];
2673 /* copy over atoms */
2674 for (j = 0; j < 2; j++)
2676 paramnew[i].a[j] = paramp2[i].a[j];
2680 for (i = p2nr; i < p2newnr; i++)
2684 /* Copy over parameters */
2685 paramnew[i].c[0] = fudgeQQ;
2686 paramnew[i].c[1] = atoms->atom[paramp1[j].a[0]].q;
2687 paramnew[i].c[2] = atoms->atom[paramp1[j].a[1]].q;
2688 paramnew[i].c[3] = paramp1[j].c[0];
2689 paramnew[i].c[4] = paramp1[j].c[1];
2691 /* copy over atoms */
2692 paramnew[i].a[0] = paramp1[j].a[0];
2693 paramnew[i].a[1] = paramp1[j].a[1];
2696 /* free the old pairlists */
2697 sfree(plist[F_LJC14_Q].param);
2698 sfree(plist[F_LJ14].param);
2700 /* now assign the new data to the F_LJC14_Q structure */
2701 plist[F_LJC14_Q].param = paramnew;
2702 plist[F_LJC14_Q].nr = p2newnr;
2704 /* Empty the LJ14 pairlist */
2705 plist[F_LJ14].nr = 0;
2706 plist[F_LJ14].param = nullptr;
2709 static void generate_LJCpairsNB(t_molinfo *mol, int nb_funct, t_params *nbp, warninp_t wi)
2711 int n, ntype, i, j, k;
2716 char errbuf[STRLEN];
2719 atom = mol->atoms.atom;
2721 ntype = static_cast<int>(std::sqrt(static_cast<double>(nbp->nr)));
2722 GMX_ASSERT(ntype * ntype == nbp->nr, "Number of pairs of generated non-bonded parameters should be a perfect square");
2724 for (i = 0; i < MAXATOMLIST; i++)
2726 param.a[i] = NOTSET;
2728 for (i = 0; i < MAXFORCEPARAM; i++)
2730 param.c[i] = NOTSET;
2733 /* Add a pair interaction for all non-excluded atom pairs */
2735 for (i = 0; i < n; i++)
2737 for (j = i+1; j < n; j++)
2740 for (k = excl->index[i]; k < excl->index[i+1]; k++)
2742 if (excl->a[k] == j)
2749 if (nb_funct != F_LJ)
2751 sprintf(errbuf, "Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones");
2752 warning_error_and_exit(wi, errbuf, FARGS);
2756 param.c[0] = atom[i].q;
2757 param.c[1] = atom[j].q;
2758 param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0];
2759 param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1];
2760 add_param_to_list(&mol->plist[F_LJC_PAIRS_NB], ¶m);
2766 static void set_excl_all(t_blocka *excl)
2770 /* Get rid of the current exclusions and exclude all atom pairs */
2772 excl->nra = nat*nat;
2773 srenew(excl->a, excl->nra);
2775 for (i = 0; i < nat; i++)
2778 for (j = 0; j < nat; j++)
2783 excl->index[nat] = k;
2786 static void decouple_atoms(t_atoms *atoms, int atomtype_decouple,
2787 int couple_lam0, int couple_lam1,
2788 const char *mol_name, warninp_t wi)
2791 char errbuf[STRLEN];
2793 for (i = 0; i < atoms->nr; i++)
2797 atom = &atoms->atom[i];
2799 if (atom->qB != atom->q || atom->typeB != atom->type)
2801 sprintf(errbuf, "Atom %d in molecule type '%s' has different A and B state charges and/or atom types set in the topology file as well as through the mdp option '%s'. You can not use both these methods simultaneously.",
2802 i + 1, mol_name, "couple-moltype");
2803 warning_error_and_exit(wi, errbuf, FARGS);
2806 if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamVDW)
2810 if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamQ)
2812 atom->type = atomtype_decouple;
2814 if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamVDW)
2818 if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamQ)
2820 atom->typeB = atomtype_decouple;
2825 void convert_moltype_couple(t_molinfo *mol, int atomtype_decouple, real fudgeQQ,
2826 int couple_lam0, int couple_lam1,
2827 gmx_bool bCoupleIntra, int nb_funct, t_params *nbp,
2830 convert_pairs_to_pairsQ(mol->plist, fudgeQQ, &mol->atoms);
2833 generate_LJCpairsNB(mol, nb_funct, nbp, wi);
2834 set_excl_all(&mol->excls);
2836 decouple_atoms(&mol->atoms, atomtype_decouple, couple_lam0, couple_lam1,