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, 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(std::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 char tmpfield[12][100]; /* Max 12 fields of width 100 */
246 bool have_atomic_number;
247 bool have_bonded_type;
252 /* First assign input line to temporary array */
253 nfields = sscanf(line, "%s%s%s%s%s%s%s%s%s%s%s%s",
254 tmpfield[0], tmpfield[1], tmpfield[2], tmpfield[3], tmpfield[4], tmpfield[5],
255 tmpfield[6], tmpfield[7], tmpfield[8], tmpfield[9], tmpfield[10], tmpfield[11]);
257 /* Comments on optional fields in the atomtypes section:
259 * The force field format is getting a bit old. For OPLS-AA we needed
260 * to add a special bonded atomtype, and for Gerrit Groenhofs QM/MM stuff
261 * we also needed the atomic numbers.
262 * To avoid making all old or user-generated force fields unusable we
263 * have introduced both these quantities as optional columns, and do some
264 * acrobatics to check whether they are present or not.
265 * This will all look much nicer when we switch to XML... sigh.
267 * Field 0 (mandatory) is the nonbonded type name. (string)
268 * Field 1 (optional) is the bonded type (string)
269 * Field 2 (optional) is the atomic number (int)
270 * Field 3 (mandatory) is the mass (numerical)
271 * Field 4 (mandatory) is the charge (numerical)
272 * Field 5 (mandatory) is the particle type (single character)
273 * This is followed by a number of nonbonded parameters.
275 * The safest way to identify the format is the particle type field.
277 * So, here is what we do:
279 * A. Read in the first six fields as strings
280 * B. If field 3 (starting from 0) is a single char, we have neither
281 * bonded_type or atomic numbers.
282 * C. If field 5 is a single char we have both.
283 * D. If field 4 is a single char we check field 1. If this begins with
284 * an alphabetical character we have bonded types, otherwise atomic numbers.
293 if ( (strlen(tmpfield[5]) == 1) && isalpha(tmpfield[5][0]) )
295 have_bonded_type = TRUE;
296 have_atomic_number = TRUE;
298 else if ( (strlen(tmpfield[3]) == 1) && isalpha(tmpfield[3][0]) )
300 have_bonded_type = FALSE;
301 have_atomic_number = FALSE;
305 have_bonded_type = ( isalpha(tmpfield[1][0]) != 0 );
306 have_atomic_number = !have_bonded_type;
309 /* optional fields */
318 if (have_atomic_number)
320 if (have_bonded_type)
322 nread = sscanf(line, "%s%s%d%lf%lf%s%lf%lf",
323 type, btype, &atomnr, &m, &q, ptype, &c[0], &c[1]);
332 /* have_atomic_number && !have_bonded_type */
333 nread = sscanf(line, "%s%d%lf%lf%s%lf%lf",
334 type, &atomnr, &m, &q, ptype, &c[0], &c[1]);
344 if (have_bonded_type)
346 /* !have_atomic_number && have_bonded_type */
347 nread = sscanf(line, "%s%s%lf%lf%s%lf%lf",
348 type, btype, &m, &q, ptype, &c[0], &c[1]);
357 /* !have_atomic_number && !have_bonded_type */
358 nread = sscanf(line, "%s%lf%lf%s%lf%lf",
359 type, &m, &q, ptype, &c[0], &c[1]);
368 if (!have_bonded_type)
373 if (!have_atomic_number)
383 if (have_atomic_number)
385 if (have_bonded_type)
387 nread = sscanf(line, "%s%s%d%lf%lf%s%lf%lf%lf",
388 type, btype, &atomnr, &m, &q, ptype, &c[0], &c[1], &c[2]);
397 /* have_atomic_number && !have_bonded_type */
398 nread = sscanf(line, "%s%d%lf%lf%s%lf%lf%lf",
399 type, &atomnr, &m, &q, ptype, &c[0], &c[1], &c[2]);
409 if (have_bonded_type)
411 /* !have_atomic_number && have_bonded_type */
412 nread = sscanf(line, "%s%s%lf%lf%s%lf%lf%lf",
413 type, btype, &m, &q, ptype, &c[0], &c[1], &c[2]);
422 /* !have_atomic_number && !have_bonded_type */
423 nread = sscanf(line, "%s%lf%lf%s%lf%lf%lf",
424 type, &m, &q, ptype, &c[0], &c[1], &c[2]);
433 if (!have_bonded_type)
438 if (!have_atomic_number)
446 sprintf(errbuf, "Invalid function type %d in push_at", nb_funct);
447 warning_error_and_exit(wi, errbuf, FARGS);
449 for (j = nfp0; (j < MAXFORCEPARAM); j++)
454 if (strlen(type) == 1 && isdigit(type[0]))
456 warning_error_and_exit(wi, "Atom type names can't be single digits.", FARGS);
459 if (strlen(btype) == 1 && isdigit(btype[0]))
461 warning_error_and_exit(wi, "Bond atom type names can't be single digits.", FARGS);
464 /* Hack to read old topologies */
465 if (gmx_strcasecmp(ptype, "D") == 0)
469 for (j = 0; (j < eptNR); j++)
471 if (gmx_strcasecmp(ptype, xl[j].entry) == 0)
478 sprintf(errbuf, "Invalid particle type %s", ptype);
479 warning_error_and_exit(wi, errbuf, FARGS);
486 for (i = 0; (i < MAXFORCEPARAM); i++)
491 if ((batype_nr = get_bond_atomtype_type(btype, bat)) == NOTSET)
493 add_bond_atomtype(bat, symtab, btype);
495 batype_nr = get_bond_atomtype_type(btype, bat);
497 if ((nr = get_atomtype_type(type, at)) != NOTSET)
499 sprintf(errbuf, "Overriding atomtype %s", type);
501 if ((nr = set_atomtype(nr, at, symtab, atom, type, param, batype_nr,
504 sprintf(errbuf, "Replacing atomtype %s failed", type);
505 warning_error_and_exit(wi, errbuf, FARGS);
508 else if ((add_atomtype(at, symtab, atom, type, param,
509 batype_nr, atomnr)) == NOTSET)
511 sprintf(errbuf, "Adding atomtype %s failed", type);
512 warning_error_and_exit(wi, errbuf, FARGS);
516 /* Add space in the non-bonded parameters matrix */
517 realloc_nb_params(at, nbparam, pair);
523 //! Return whether the contents of \c a and \c b are the same, considering also reversed order.
524 template <typename T>
525 static bool equalEitherForwardOrBackward(gmx::ArrayRef<const T> a, gmx::ArrayRef<const T> b)
527 return (std::equal(a.begin(), a.end(), b.begin()) ||
528 std::equal(a.begin(), a.end(), b.rbegin()));
531 static void push_bondtype(t_params * bt,
540 int nrfp = NRFP(ftype);
543 /* If bAllowRepeat is TRUE, we allow multiple entries as long as they
544 are on directly _adjacent_ lines.
547 /* First check if our atomtypes are _identical_ (not reversed) to the previous
548 entry. If they are not identical we search for earlier duplicates. If they are
549 we can skip it, since we already searched for the first line
553 bool isContinuationOfBlock = false;
554 if (bAllowRepeat && nr > 1)
556 isContinuationOfBlock = true;
557 for (int j = 0; j < nral; j++)
559 if (b->a[j] != bt->param[nr - 2].a[j])
561 isContinuationOfBlock = false;
566 /* Search for earlier duplicates if this entry was not a continuation
567 from the previous line.
569 bool addBondType = true;
570 bool haveWarned = false;
571 bool haveErrored = false;
572 for (int i = 0; (i < nr); i++)
574 gmx::ArrayRef<const int> bParams(b->a, b->a + nral);
575 gmx::ArrayRef<const int> testParams(bt->param[i].a, bt->param[i].a + nral);
576 if (equalEitherForwardOrBackward(bParams, testParams))
578 GMX_ASSERT(nrfp <= MAXFORCEPARAM, "This is ensured in other places, but we need this assert to keep the clang analyzer happy");
579 const bool identicalParameters = std::equal(bt->param[i].c, bt->param[i].c + nrfp, b->c);
581 if (!bAllowRepeat || identicalParameters)
586 if (!identicalParameters)
590 /* With dihedral type 9 we only allow for repeating
591 * of the same parameters with blocks with 1 entry.
592 * Allowing overriding is too complex to check.
594 if (!isContinuationOfBlock && !haveErrored)
596 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.");
600 else if (!haveWarned)
602 sprintf(errbuf, "Overriding %s parameters.%s",
603 interaction_function[ftype].longname,
605 "\nUse dihedraltype 9 to allow several multiplicity terms. Only consecutive lines are combined. Non-consective lines overwrite each other."
609 fprintf(stderr, " old: ");
610 for (int j = 0; (j < nrfp); j++)
612 fprintf(stderr, " %g", bt->param[i].c[j]);
614 fprintf(stderr, " \n new: %s\n\n", line);
620 if (!identicalParameters && !bAllowRepeat)
622 /* Overwrite the parameters with the latest ones */
623 // TODO considering improving the following code by replacing with:
624 // std::copy(b->c, b->c + nrfp, bt->param[i].c);
625 for (int j = 0; (j < nrfp); j++)
627 bt->param[i].c[j] = b->c[j];
638 /* fill the arrays up and down */
639 memcpy(bt->param[bt->nr].c, b->c, sizeof(b->c));
640 memcpy(bt->param[bt->nr].a, b->a, sizeof(b->a));
641 memcpy(bt->param[bt->nr+1].c, b->c, sizeof(b->c));
643 /* The definitions of linear angles depend on the order of atoms,
644 * that means that for atoms i-j-k, with certain parameter a, the
645 * corresponding k-j-i angle will have parameter 1-a.
647 if (ftype == F_LINEAR_ANGLES)
649 bt->param[bt->nr+1].c[0] = 1-bt->param[bt->nr+1].c[0];
650 bt->param[bt->nr+1].c[2] = 1-bt->param[bt->nr+1].c[2];
653 for (int j = 0; (j < nral); j++)
655 bt->param[bt->nr+1].a[j] = b->a[nral-1-j];
662 void push_bt(directive d, t_params bt[], int nral,
664 t_bond_atomtype bat, char *line,
667 const char *formal[MAXATOMLIST+1] = {
676 const char *formnl[MAXATOMLIST+1] = {
682 "%*s%*s%*s%*s%*s%*s",
683 "%*s%*s%*s%*s%*s%*s%*s"
685 const char *formlf = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
686 int i, ft, ftype, nn, nrfp, nrfpA;
688 char alc[MAXATOMLIST+1][20];
689 /* One force parameter more, so we can check if we read too many */
690 double c[MAXFORCEPARAM+1];
694 if ((bat && at) || (!bat && !at))
696 gmx_incons("You should pass either bat or at to push_bt");
699 /* Make format string (nral ints+functype) */
700 if ((nn = sscanf(line, formal[nral],
701 alc[0], alc[1], alc[2], alc[3], alc[4], alc[5])) != nral+1)
703 sprintf(errbuf, "Not enough atomtypes (%d instead of %d)", nn-1, nral);
704 warning_error(wi, errbuf);
708 ft = strtol(alc[nral], nullptr, 10);
709 ftype = ifunc_index(d, ft);
711 nrfpA = interaction_function[ftype].nrfpA;
712 strcpy(f1, formnl[nral]);
714 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]))
719 /* Copy the B-state from the A-state */
720 copy_B_from_A(ftype, c);
726 warning_error(wi, "Not enough parameters");
728 else if (nn > nrfpA && nn < nrfp)
730 warning_error(wi, "Too many parameters or not enough parameters for topology B");
734 warning_error(wi, "Too many parameters");
736 for (i = nn; (i < nrfp); i++)
742 for (i = 0; (i < nral); i++)
744 if (at && ((p.a[i] = get_atomtype_type(alc[i], at)) == NOTSET))
746 sprintf(errbuf, "Unknown atomtype %s\n", alc[i]);
747 warning_error_and_exit(wi, errbuf, FARGS);
749 else if (bat && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET))
751 sprintf(errbuf, "Unknown bond_atomtype %s\n", alc[i]);
752 warning_error_and_exit(wi, errbuf, FARGS);
755 for (i = 0; (i < nrfp); i++)
759 push_bondtype (&(bt[ftype]), &p, nral, ftype, FALSE, line, wi);
763 void push_dihedraltype(directive d, t_params bt[],
764 t_bond_atomtype bat, char *line,
767 const char *formal[MAXATOMLIST+1] = {
776 const char *formnl[MAXATOMLIST+1] = {
782 "%*s%*s%*s%*s%*s%*s",
783 "%*s%*s%*s%*s%*s%*s%*s"
785 const char *formlf[MAXFORCEPARAM] = {
791 "%lf%lf%lf%lf%lf%lf",
792 "%lf%lf%lf%lf%lf%lf%lf",
793 "%lf%lf%lf%lf%lf%lf%lf%lf",
794 "%lf%lf%lf%lf%lf%lf%lf%lf%lf",
795 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
796 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
797 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
799 int i, ft, ftype, nn, nrfp, nrfpA, nral;
801 char alc[MAXATOMLIST+1][20];
802 double c[MAXFORCEPARAM];
807 /* This routine accepts dihedraltypes defined from either 2 or 4 atoms.
809 * We first check for 2 atoms with the 3th column being an integer
810 * defining the type. If this isn't the case, we try it with 4 atoms
811 * and the 5th column defining the dihedral type.
813 nn = sscanf(line, formal[4], alc[0], alc[1], alc[2], alc[3], alc[4]);
814 if (nn >= 3 && strlen(alc[2]) == 1 && isdigit(alc[2][0]))
817 ft = strtol(alc[nral], nullptr, 10);
818 /* Move atom types around a bit and use 'X' for wildcard atoms
819 * to create a 4-atom dihedral definition with arbitrary atoms in
822 if (alc[2][0] == '2')
824 /* improper - the two atomtypes are 1,4. Use wildcards for 2,3 */
825 strcpy(alc[3], alc[1]);
826 sprintf(alc[2], "X");
827 sprintf(alc[1], "X");
828 /* alc[0] stays put */
832 /* proper - the two atomtypes are 2,3. Use wildcards for 1,4 */
833 sprintf(alc[3], "X");
834 strcpy(alc[2], alc[1]);
835 strcpy(alc[1], alc[0]);
836 sprintf(alc[0], "X");
839 else if (nn == 5 && strlen(alc[4]) == 1 && isdigit(alc[4][0]))
842 ft = strtol(alc[nral], nullptr, 10);
846 sprintf(errbuf, "Incorrect number of atomtypes for dihedral (%d instead of 2 or 4)", nn-1);
847 warning_error(wi, errbuf);
853 /* Previously, we have always overwritten parameters if e.g. a torsion
854 with the same atomtypes occurs on multiple lines. However, CHARMM and
855 some other force fields specify multiple dihedrals over some bonds,
856 including cosines with multiplicity 6 and somethimes even higher.
857 Thus, they cannot be represented with Ryckaert-Bellemans terms.
858 To add support for these force fields, Dihedral type 9 is identical to
859 normal proper dihedrals, but repeated entries are allowed.
866 bAllowRepeat = FALSE;
870 ftype = ifunc_index(d, ft);
872 nrfpA = interaction_function[ftype].nrfpA;
874 strcpy(f1, formnl[nral]);
875 strcat(f1, formlf[nrfp-1]);
877 /* Check number of parameters given */
878 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]))
883 /* Copy the B-state from the A-state */
884 copy_B_from_A(ftype, c);
890 warning_error(wi, "Not enough parameters");
892 else if (nn > nrfpA && nn < nrfp)
894 warning_error(wi, "Too many parameters or not enough parameters for topology B");
898 warning_error(wi, "Too many parameters");
900 for (i = nn; (i < nrfp); i++)
907 for (i = 0; (i < 4); i++)
909 if (!strcmp(alc[i], "X"))
915 if ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET)
917 sprintf(errbuf, "Unknown bond_atomtype %s", alc[i]);
918 warning_error_and_exit(wi, errbuf, FARGS);
922 for (i = 0; (i < nrfp); i++)
926 /* Always use 4 atoms here, since we created two wildcard atoms
927 * if there wasn't of them 4 already.
929 push_bondtype (&(bt[ftype]), &p, 4, ftype, bAllowRepeat, line, wi);
933 void push_nbt(directive d, t_nbparam **nbt, gpp_atomtype_t atype,
934 char *pline, int nb_funct,
938 const char *form3 = "%*s%*s%*s%lf%lf%lf";
939 const char *form4 = "%*s%*s%*s%lf%lf%lf%lf";
940 const char *form5 = "%*s%*s%*s%lf%lf%lf%lf%lf";
942 int i, f, n, ftype, nrfp;
950 if (sscanf (pline, "%s%s%d", a0, a1, &f) != 3)
956 ftype = ifunc_index(d, f);
958 if (ftype != nb_funct)
960 sprintf(errbuf, "Trying to add %s while the default nonbond type is %s",
961 interaction_function[ftype].longname,
962 interaction_function[nb_funct].longname);
963 warning_error(wi, errbuf);
967 /* Get the force parameters */
971 n = sscanf(pline, form4, &c[0], &c[1], &c[2], &c[3]);
977 /* When the B topology parameters are not set,
978 * copy them from topology A
980 GMX_ASSERT(nrfp <= 4, "LJ-14 cannot have more than 4 parameters");
981 for (i = n; i < nrfp; i++)
986 else if (ftype == F_LJC14_Q)
988 n = sscanf(pline, form5, &c[0], &c[1], &c[2], &c[3], &dum);
991 incorrect_n_param(wi);
997 if (sscanf(pline, form3, &c[0], &c[1], &dum) != 2)
999 incorrect_n_param(wi);
1005 if (sscanf(pline, form4, &c[0], &c[1], &c[2], &dum) != 3)
1007 incorrect_n_param(wi);
1013 sprintf(errbuf, "Number of force parameters for nonbonded interactions is %d", nrfp);
1014 warning_error_and_exit(wi, errbuf, FARGS);
1016 for (i = 0; (i < nrfp); i++)
1021 /* Put the parameters in the matrix */
1022 if ((ai = get_atomtype_type (a0, atype)) == NOTSET)
1024 sprintf(errbuf, "Atomtype %s not found", a0);
1025 warning_error_and_exit(wi, errbuf, FARGS);
1027 if ((aj = get_atomtype_type (a1, atype)) == NOTSET)
1029 sprintf(errbuf, "Atomtype %s not found", a1);
1030 warning_error_and_exit(wi, errbuf, FARGS);
1032 nbp = &(nbt[std::max(ai, aj)][std::min(ai, aj)]);
1037 for (i = 0; i < nrfp; i++)
1039 bId = bId && (nbp->c[i] == cr[i]);
1043 sprintf(errbuf, "Overriding non-bonded parameters,");
1044 warning(wi, errbuf);
1045 fprintf(stderr, " old:");
1046 for (i = 0; i < nrfp; i++)
1048 fprintf(stderr, " %g", nbp->c[i]);
1050 fprintf(stderr, " new\n%s\n", pline);
1054 for (i = 0; i < nrfp; i++)
1061 push_cmaptype(directive d, t_params bt[], int nral, gpp_atomtype_t at,
1062 t_bond_atomtype bat, char *line,
1065 const char *formal = "%s%s%s%s%s%s%s%s%n";
1067 int i, ft, ftype, nn, nrfp, nrfpA, nrfpB;
1068 int start, nchar_consumed;
1069 int nxcmap, nycmap, ncmap, read_cmap, sl, nct;
1070 char s[20], alc[MAXATOMLIST+2][20];
1072 char errbuf[STRLEN];
1074 /* Keep the compiler happy */
1078 GMX_ASSERT(nral == 5, "CMAP requires 5 atoms per interaction");
1080 /* Here we can only check for < 8 */
1081 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)
1083 sprintf(errbuf, "Incorrect number of atomtypes for cmap (%d instead of 5)", nn-1);
1084 warning_error(wi, errbuf);
1087 start += nchar_consumed;
1089 ft = strtol(alc[nral], nullptr, 10);
1090 nxcmap = strtol(alc[nral+1], nullptr, 10);
1091 nycmap = strtol(alc[nral+2], nullptr, 10);
1093 /* Check for equal grid spacing in x and y dims */
1094 if (nxcmap != nycmap)
1096 sprintf(errbuf, "Not the same grid spacing in x and y for cmap grid: x=%d, y=%d", nxcmap, nycmap);
1097 warning_error(wi, errbuf);
1100 ncmap = nxcmap*nycmap;
1101 ftype = ifunc_index(d, ft);
1102 nrfpA = strtol(alc[6], nullptr, 10)*strtol(alc[6], nullptr, 10);
1103 nrfpB = strtol(alc[7], nullptr, 10)*strtol(alc[7], nullptr, 10);
1106 /* Allocate memory for the CMAP grid */
1107 bt[F_CMAP].ncmap += nrfp;
1108 srenew(bt[F_CMAP].cmap, bt[F_CMAP].ncmap);
1110 /* Read in CMAP parameters */
1112 for (i = 0; i < ncmap; i++)
1114 while (isspace(*(line+start+sl)))
1118 nn = sscanf(line+start+sl, " %s ", s);
1120 bt[F_CMAP].cmap[i+(bt[F_CMAP].ncmap)-nrfp] = strtod(s, nullptr);
1128 sprintf(errbuf, "Error in reading cmap parameter for angle %s %s %s %s %s", alc[0], alc[1], alc[2], alc[3], alc[4]);
1129 warning_error(wi, errbuf);
1134 /* Check do that we got the number of parameters we expected */
1135 if (read_cmap == nrfpA)
1137 for (i = 0; i < ncmap; i++)
1139 bt[F_CMAP].cmap[i+ncmap] = bt[F_CMAP].cmap[i];
1144 if (read_cmap < nrfpA)
1146 warning_error(wi, "Not enough cmap parameters");
1148 else if (read_cmap > nrfpA && read_cmap < nrfp)
1150 warning_error(wi, "Too many cmap parameters or not enough parameters for topology B");
1152 else if (read_cmap > nrfp)
1154 warning_error(wi, "Too many cmap parameters");
1159 /* Set grid spacing and the number of grids (we assume these numbers to be the same for all grids
1160 * so we can safely assign them each time
1162 bt[F_CMAP].grid_spacing = nxcmap; /* Or nycmap, they need to be equal */
1163 bt[F_CMAP].nc = bt[F_CMAP].nc + 1; /* Since we are incrementing here, we need to subtract later, see (*****) */
1164 nct = (nral+1) * bt[F_CMAP].nc;
1166 /* Allocate memory for the cmap_types information */
1167 srenew(bt[F_CMAP].cmap_types, nct);
1169 for (i = 0; (i < nral); i++)
1171 if (at && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET))
1173 sprintf(errbuf, "Unknown atomtype %s\n", alc[i]);
1174 warning_error(wi, errbuf);
1176 else if (bat && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET))
1178 sprintf(errbuf, "Unknown bond_atomtype %s\n", alc[i]);
1179 warning_error(wi, errbuf);
1182 /* Assign a grid number to each cmap_type */
1183 bt[F_CMAP].cmap_types[bt[F_CMAP].nct++] = get_bond_atomtype_type(alc[i], bat);
1186 /* Assign a type number to this cmap */
1187 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 (****) */
1189 /* Check for the correct number of atoms (again) */
1190 if (bt[F_CMAP].nct != nct)
1192 sprintf(errbuf, "Incorrect number of atom types (%d) in cmap type %d\n", nct, bt[F_CMAP].nc);
1193 warning_error(wi, errbuf);
1196 /* Is this correct?? */
1197 for (i = 0; i < MAXFORCEPARAM; i++)
1202 /* Push the bond to the bondlist */
1203 push_bondtype (&(bt[ftype]), &p, nral, ftype, FALSE, line, wi);
1207 static void push_atom_now(t_symtab *symtab, t_atoms *at, int atomnr,
1209 int type, char *ctype, int ptype,
1211 char *resname, char *name, real m0, real q0,
1212 int typeB, char *ctypeB, real mB, real qB,
1215 int j, resind = 0, resnr;
1218 char errbuf[STRLEN];
1220 if (((nr == 0) && (atomnr != 1)) || (nr && (atomnr != at->nr+1)))
1222 sprintf(errbuf, "Atoms in the .top are not numbered consecutively from 1 (rather, atomnr = %d, while at->nr = %d)", atomnr, at->nr);
1223 warning_error_and_exit(wi, errbuf, FARGS);
1226 j = strlen(resnumberic) - 1;
1227 if (isdigit(resnumberic[j]))
1233 ric = resnumberic[j];
1234 if (j == 0 || !isdigit(resnumberic[j-1]))
1236 sprintf(errbuf, "Invalid residue number '%s' for atom %d",
1237 resnumberic, atomnr);
1238 warning_error_and_exit(wi, errbuf, FARGS);
1241 resnr = strtol(resnumberic, nullptr, 10);
1245 resind = at->atom[nr-1].resind;
1247 if (nr == 0 || strcmp(resname, *at->resinfo[resind].name) != 0 ||
1248 resnr != at->resinfo[resind].nr ||
1249 ric != at->resinfo[resind].ic)
1259 at->nres = resind + 1;
1260 srenew(at->resinfo, at->nres);
1261 at->resinfo[resind].name = put_symtab(symtab, resname);
1262 at->resinfo[resind].nr = resnr;
1263 at->resinfo[resind].ic = ric;
1267 resind = at->atom[at->nr-1].resind;
1270 /* New atom instance
1271 * get new space for arrays
1273 srenew(at->atom, nr+1);
1274 srenew(at->atomname, nr+1);
1275 srenew(at->atomtype, nr+1);
1276 srenew(at->atomtypeB, nr+1);
1279 at->atom[nr].type = type;
1280 at->atom[nr].ptype = ptype;
1281 at->atom[nr].q = q0;
1282 at->atom[nr].m = m0;
1283 at->atom[nr].typeB = typeB;
1284 at->atom[nr].qB = qB;
1285 at->atom[nr].mB = mB;
1287 at->atom[nr].resind = resind;
1288 at->atom[nr].atomnumber = atomicnumber;
1289 at->atomname[nr] = put_symtab(symtab, name);
1290 at->atomtype[nr] = put_symtab(symtab, ctype);
1291 at->atomtypeB[nr] = put_symtab(symtab, ctypeB);
1295 static void push_cg(t_block *block, int *lastindex, int index, int a)
1297 if (((block->nr) && (*lastindex != index)) || (!block->nr))
1299 /* add a new block */
1301 srenew(block->index, block->nr+1);
1303 block->index[block->nr] = a + 1;
1307 void push_atom(t_symtab *symtab, t_block *cgs,
1308 t_atoms *at, gpp_atomtype_t atype, char *line, int *lastcg,
1312 int cgnumber, atomnr, type, typeB, nscan;
1313 char id[STRLEN], ctype[STRLEN], ctypeB[STRLEN],
1314 resnumberic[STRLEN], resname[STRLEN], name[STRLEN], check[STRLEN];
1315 double m, q, mb, qb;
1316 real m0, q0, mB, qB;
1317 char errbuf[STRLEN];
1319 /* Make a shortcut for writing in this molecule */
1322 /* Fixed parameters */
1323 if (sscanf(line, "%s%s%s%s%s%d",
1324 id, ctype, resnumberic, resname, name, &cgnumber) != 6)
1329 sscanf(id, "%d", &atomnr);
1330 if ((type = get_atomtype_type(ctype, atype)) == NOTSET)
1332 sprintf(errbuf, "Atomtype %s not found", ctype);
1333 warning_error_and_exit(wi, errbuf, FARGS);
1335 ptype = get_atomtype_ptype(type, atype);
1337 /* Set default from type */
1338 q0 = get_atomtype_qA(type, atype);
1339 m0 = get_atomtype_massA(type, atype);
1344 /* Optional parameters */
1345 nscan = sscanf(line, "%*s%*s%*s%*s%*s%*s%lf%lf%s%lf%lf%s",
1346 &q, &m, ctypeB, &qb, &mb, check);
1348 /* Nasty switch that falls thru all the way down! */
1357 if ((typeB = get_atomtype_type(ctypeB, atype)) == NOTSET)
1359 sprintf(errbuf, "Atomtype %s not found", ctypeB);
1360 warning_error_and_exit(wi, errbuf, FARGS);
1362 qB = get_atomtype_qA(typeB, atype);
1363 mB = get_atomtype_massA(typeB, atype);
1372 warning_error(wi, "Too many parameters");
1380 push_cg(cgs, lastcg, cgnumber, nr);
1382 push_atom_now(symtab, at, atomnr, get_atomtype_atomnumber(type, atype),
1383 type, ctype, ptype, resnumberic,
1384 resname, name, m0, q0, typeB,
1385 typeB == type ? ctype : ctypeB, mB, qB, wi);
1388 void push_molt(t_symtab *symtab, int *nmol, t_molinfo **mol, char *line,
1394 char errbuf[STRLEN];
1396 if ((sscanf(line, "%s%d", type, &nrexcl)) != 2)
1398 warning_error(wi, "Expected a molecule type name and nrexcl");
1401 /* Test if this moleculetype overwrites another */
1405 if (strcmp(*((*mol)[i].name), type) == 0)
1407 sprintf(errbuf, "moleculetype %s is redefined", type);
1408 warning_error_and_exit(wi, errbuf, FARGS);
1414 srenew(*mol, *nmol);
1415 newmol = &((*mol)[*nmol-1]);
1416 init_molinfo(newmol);
1418 /* Fill in the values */
1419 newmol->name = put_symtab(symtab, type);
1420 newmol->nrexcl = nrexcl;
1421 newmol->excl_set = FALSE;
1424 static bool default_nb_params(int ftype, t_params bt[], t_atoms *at,
1425 t_param *p, int c_start, bool bB, bool bGenPairs)
1427 int i, j, ti, tj, ntype;
1429 t_param *pi = nullptr;
1430 int nr = bt[ftype].nr;
1431 int nral = NRAL(ftype);
1432 int nrfp = interaction_function[ftype].nrfpA;
1433 int nrfpB = interaction_function[ftype].nrfpB;
1435 if ((!bB && nrfp == 0) || (bB && nrfpB == 0))
1443 /* First test the generated-pair position to save
1444 * time when we have 1000*1000 entries for e.g. OPLS...
1446 ntype = static_cast<int>(std::sqrt(static_cast<double>(nr)));
1447 GMX_ASSERT(ntype * ntype == nr, "Number of pairs of generated non-bonded parameters should be a perfect square");
1450 ti = at->atom[p->a[0]].typeB;
1451 tj = at->atom[p->a[1]].typeB;
1455 ti = at->atom[p->a[0]].type;
1456 tj = at->atom[p->a[1]].type;
1458 pi = &(bt[ftype].param[ntype*ti+tj]);
1459 bFound = ((ti == pi->a[0]) && (tj == pi->a[1]));
1462 /* Search explicitly if we didnt find it */
1465 for (i = 0; ((i < nr) && !bFound); i++)
1467 pi = &(bt[ftype].param[i]);
1470 for (j = 0; ((j < nral) &&
1471 (at->atom[p->a[j]].typeB == pi->a[j])); j++)
1478 for (j = 0; ((j < nral) &&
1479 (at->atom[p->a[j]].type == pi->a[j])); j++)
1484 bFound = (j == nral);
1492 if (nrfp+nrfpB > MAXFORCEPARAM)
1494 gmx_incons("Too many force parameters");
1496 for (j = c_start; (j < nrfpB); j++)
1498 p->c[nrfp+j] = pi->c[j];
1503 for (j = c_start; (j < nrfp); j++)
1511 for (j = c_start; (j < nrfp); j++)
1519 static bool default_cmap_params(t_params bondtype[],
1520 t_atoms *at, gpp_atomtype_t atype,
1521 t_param *p, bool bB,
1522 int *cmap_type, int *nparam_def,
1525 int i, nparam_found;
1527 bool bFound = FALSE;
1528 char errbuf[STRLEN];
1533 /* Match the current cmap angle against the list of cmap_types */
1534 for (i = 0; i < bondtype[F_CMAP].nct && !bFound; i += 6)
1543 (get_atomtype_batype(at->atom[p->a[0]].type, atype) == bondtype[F_CMAP].cmap_types[i]) &&
1544 (get_atomtype_batype(at->atom[p->a[1]].type, atype) == bondtype[F_CMAP].cmap_types[i+1]) &&
1545 (get_atomtype_batype(at->atom[p->a[2]].type, atype) == bondtype[F_CMAP].cmap_types[i+2]) &&
1546 (get_atomtype_batype(at->atom[p->a[3]].type, atype) == bondtype[F_CMAP].cmap_types[i+3]) &&
1547 (get_atomtype_batype(at->atom[p->a[4]].type, atype) == bondtype[F_CMAP].cmap_types[i+4]))
1549 /* Found cmap torsion */
1551 ct = bondtype[F_CMAP].cmap_types[i+5];
1557 /* If we did not find a matching type for this cmap torsion */
1560 sprintf(errbuf, "Unknown cmap torsion between atoms %d %d %d %d %d",
1561 p->a[0]+1, p->a[1]+1, p->a[2]+1, p->a[3]+1, p->a[4]+1);
1562 warning_error_and_exit(wi, errbuf, FARGS);
1565 *nparam_def = nparam_found;
1571 /* Returns the number of exact atom type matches, i.e. non wild-card matches,
1572 * returns -1 when there are no matches at all.
1574 static int natom_match(t_param *pi,
1575 int type_i, int type_j, int type_k, int type_l,
1576 const gpp_atomtype* atype)
1578 if ((pi->ai() == -1 || get_atomtype_batype(type_i, atype) == pi->ai()) &&
1579 (pi->aj() == -1 || get_atomtype_batype(type_j, atype) == pi->aj()) &&
1580 (pi->ak() == -1 || get_atomtype_batype(type_k, atype) == pi->ak()) &&
1581 (pi->al() == -1 || get_atomtype_batype(type_l, atype) == pi->al()))
1584 (pi->ai() == -1 ? 0 : 1) +
1585 (pi->aj() == -1 ? 0 : 1) +
1586 (pi->ak() == -1 ? 0 : 1) +
1587 (pi->al() == -1 ? 0 : 1);
1595 static bool default_params(int ftype, t_params bt[],
1596 t_atoms *at, gpp_atomtype_t atype,
1597 t_param *p, bool bB,
1598 t_param **param_def,
1603 t_param *pi = nullptr;
1604 t_param *pj = nullptr;
1605 int nr = bt[ftype].nr;
1606 int nral = NRAL(ftype);
1607 int nrfpA = interaction_function[ftype].nrfpA;
1608 int nrfpB = interaction_function[ftype].nrfpB;
1610 if ((!bB && nrfpA == 0) || (bB && nrfpB == 0))
1618 if (ftype == F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS || ftype == F_PIDIHS)
1620 int nmatch_max = -1;
1624 /* For dihedrals we allow wildcards. We choose the first type
1625 * that has the most real matches, i.e. non-wildcard matches.
1627 for (t = 0; ((t < nr) && nmatch_max < 4); t++)
1632 pt = &(bt[ftype].param[t]);
1635 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);
1639 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);
1641 if (nmatch > nmatch_max)
1643 nmatch_max = nmatch;
1653 pi = &(bt[ftype].param[i]);
1656 /* Find additional matches for this dihedral - necessary
1658 * The rule in that case is that additional matches
1659 * HAVE to be on adjacent lines!
1662 /* Continue from current i value */
1663 for (j = i + 2; j < nr && bSame; j += 2)
1665 pj = &(bt[ftype].param[j]);
1666 bSame = (pi->ai() == pj->ai() && pi->aj() == pj->aj() && pi->ak() == pj->ak() && pi->al() == pj->al());
1671 /* nparam_found will be increased as long as the numbers match */
1675 else /* Not a dihedral */
1679 for (i = 0; ((i < nr) && !bFound); i++)
1681 pi = &(bt[ftype].param[i]);
1684 for (j = 0; ((j < nral) &&
1685 (get_atomtype_batype(at->atom[p->a[j]].typeB, atype) == pi->a[j])); j++)
1692 for (j = 0; ((j < nral) &&
1693 (get_atomtype_batype(at->atom[p->a[j]].type, atype) == pi->a[j])); j++)
1698 bFound = (j == nral);
1707 *nparam_def = nparam_found;
1714 void push_bond(directive d, t_params bondtype[], t_params bond[],
1715 t_atoms *at, gpp_atomtype_t atype, char *line,
1716 bool bBonded, bool bGenPairs, real fudgeQQ,
1717 bool bZero, bool *bWarn_copy_A_B,
1720 const char *aaformat[MAXATOMLIST] = {
1728 const char *asformat[MAXATOMLIST] = {
1733 "%*s%*s%*s%*s%*s%*s",
1734 "%*s%*s%*s%*s%*s%*s%*s"
1736 const char *ccformat = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
1737 int nr, i, j, nral, nral_fmt, nread, ftype;
1738 char format[STRLEN];
1739 /* One force parameter more, so we can check if we read too many */
1740 double cc[MAXFORCEPARAM+1];
1741 int aa[MAXATOMLIST+1];
1742 t_param param, *param_defA, *param_defB;
1743 bool bFoundA = FALSE, bFoundB = FALSE, bDef, bPert, bSwapParity = FALSE;
1744 int nparam_defA, nparam_defB;
1745 char errbuf[STRLEN];
1747 nparam_defA = nparam_defB = 0;
1749 ftype = ifunc_index(d, 1);
1751 for (j = 0; j < MAXATOMLIST; j++)
1755 bDef = (NRFP(ftype) > 0);
1757 if (ftype == F_SETTLE)
1759 /* SETTLE acts on 3 atoms, but the topology format only specifies
1760 * the first atom (for historical reasons).
1769 nread = sscanf(line, aaformat[nral_fmt-1],
1770 &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
1772 if (ftype == F_SETTLE)
1779 if (nread < nral_fmt)
1784 else if (nread > nral_fmt)
1786 /* this is a hack to allow for virtual sites with swapped parity */
1787 bSwapParity = (aa[nral] < 0);
1790 aa[nral] = -aa[nral];
1792 ftype = ifunc_index(d, aa[nral]);
1801 sprintf(errbuf, "Negative function types only allowed for %s and %s",
1802 interaction_function[F_VSITE3FAD].longname,
1803 interaction_function[F_VSITE3OUT].longname);
1804 warning_error_and_exit(wi, errbuf, FARGS);
1810 /* Check for double atoms and atoms out of bounds */
1811 for (i = 0; (i < nral); i++)
1813 if (aa[i] < 1 || aa[i] > at->nr)
1815 sprintf(errbuf, "Atom index (%d) in %s out of bounds (1-%d).\n"
1816 "This probably means that you have inserted topology section \"%s\"\n"
1817 "in a part belonging to a different molecule than you intended to.\n"
1818 "In that case move the \"%s\" section to the right molecule.",
1819 aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
1820 warning_error_and_exit(wi, errbuf, FARGS);
1822 for (j = i+1; (j < nral); j++)
1824 GMX_ASSERT(j < MAXATOMLIST + 1, "Values from nral=NRAL() will satisfy this, we assert to keep gcc 4 happy");
1827 sprintf(errbuf, "Duplicate atom index (%d) in %s", aa[i], dir2str(d));
1828 if (ftype == F_ANGRES)
1830 /* Since the angle restraints uses 2 pairs of atoms to
1831 * defines an angle between vectors, it can be useful
1832 * to use one atom twice, so we only issue a note here.
1834 warning_note(wi, errbuf);
1838 warning_error(wi, errbuf);
1844 /* default force parameters */
1845 for (j = 0; (j < MAXATOMLIST); j++)
1847 param.a[j] = aa[j]-1;
1849 for (j = 0; (j < MAXFORCEPARAM); j++)
1854 /* Get force params for normal and free energy perturbation
1855 * studies, as determined by types!
1860 bFoundA = default_params(ftype, bondtype, at, atype, ¶m, FALSE, ¶m_defA, &nparam_defA);
1863 /* Copy the A-state and B-state default parameters. */
1864 GMX_ASSERT(NRFPA(ftype)+NRFPB(ftype) <= MAXFORCEPARAM, "Bonded interactions may have at most 12 parameters");
1865 for (j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++)
1867 param.c[j] = param_defA->c[j];
1870 bFoundB = default_params(ftype, bondtype, at, atype, ¶m, TRUE, ¶m_defB, &nparam_defB);
1873 /* Copy only the B-state default parameters */
1874 for (j = NRFPA(ftype); (j < NRFP(ftype)); j++)
1876 param.c[j] = param_defB->c[j];
1880 else if (ftype == F_LJ14)
1882 bFoundA = default_nb_params(ftype, bondtype, at, ¶m, 0, FALSE, bGenPairs);
1883 bFoundB = default_nb_params(ftype, bondtype, at, ¶m, 0, TRUE, bGenPairs);
1885 else if (ftype == F_LJC14_Q)
1887 param.c[0] = fudgeQQ;
1888 /* Fill in the A-state charges as default parameters */
1889 param.c[1] = at->atom[param.a[0]].q;
1890 param.c[2] = at->atom[param.a[1]].q;
1891 /* The default LJ parameters are the standard 1-4 parameters */
1892 bFoundA = default_nb_params(F_LJ14, bondtype, at, ¶m, 3, FALSE, bGenPairs);
1895 else if (ftype == F_LJC_PAIRS_NB)
1897 /* Defaults are not supported here */
1903 gmx_incons("Unknown function type in push_bond");
1906 if (nread > nral_fmt)
1908 /* Manually specified parameters - in this case we discard multiple torsion info! */
1910 strcpy(format, asformat[nral_fmt-1]);
1911 strcat(format, ccformat);
1913 nread = sscanf(line, format, &cc[0], &cc[1], &cc[2], &cc[3], &cc[4], &cc[5],
1914 &cc[6], &cc[7], &cc[8], &cc[9], &cc[10], &cc[11], &cc[12]);
1916 if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0))
1918 /* We only have to issue a warning if these atoms are perturbed! */
1920 for (j = 0; (j < nral); j++)
1922 bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1925 if (bPert && *bWarn_copy_A_B)
1928 "Some parameters for bonded interaction involving perturbed atoms are specified explicitly in state A, but not B - copying A to B");
1929 warning(wi, errbuf);
1930 *bWarn_copy_A_B = FALSE;
1933 /* If only the A parameters were specified, copy them to the B state */
1934 /* The B-state parameters correspond to the first nrfpB
1935 * A-state parameters.
1937 for (j = 0; (j < NRFPB(ftype)); j++)
1939 cc[nread++] = cc[j];
1943 /* If nread was 0 or EOF, no parameters were read => use defaults.
1944 * If nread was nrfpA we copied above so nread=nrfp.
1945 * If nread was nrfp we are cool.
1946 * For F_LJC14_Q we allow supplying fudgeQQ only.
1947 * Anything else is an error!
1949 if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) &&
1950 !(ftype == F_LJC14_Q && nread == 1))
1952 sprintf(errbuf, "Incorrect number of parameters - found %d, expected %d or %d for %s (after the function type).",
1953 nread, NRFPA(ftype), NRFP(ftype), interaction_function[ftype].longname);
1954 warning_error_and_exit(wi, errbuf, FARGS);
1957 for (j = 0; (j < nread); j++)
1962 /* Check whether we have to use the defaults */
1963 if (nread == NRFP(ftype))
1972 /* nread now holds the number of force parameters read! */
1977 /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */
1978 if (ftype == F_PDIHS)
1980 if ((nparam_defA != nparam_defB) || ((nparam_defA > 1 || nparam_defB > 1) && (param_defA != param_defB)))
1983 "Cannot automatically perturb a torsion with multiple terms to different form.\n"
1984 "Please specify perturbed parameters manually for this torsion in your topology!");
1985 warning_error(wi, errbuf);
1989 if (nread > 0 && nread < NRFPA(ftype))
1991 /* Issue an error, do not use defaults */
1992 sprintf(errbuf, "Not enough parameters, there should be at least %d (or 0 for defaults)", NRFPA(ftype));
1993 warning_error(wi, errbuf);
1996 if (nread == 0 || nread == EOF)
2000 if (interaction_function[ftype].flags & IF_VSITE)
2002 /* set them to NOTSET, will be calculated later */
2003 for (j = 0; (j < MAXFORCEPARAM); j++)
2005 param.c[j] = NOTSET;
2010 param.c1() = -1; /* flag to swap parity of vsite construction */
2017 fprintf(stderr, "NOTE: No default %s types, using zeroes\n",
2018 interaction_function[ftype].longname);
2022 sprintf(errbuf, "No default %s types", interaction_function[ftype].longname);
2023 warning_error(wi, errbuf);
2034 param.c0() = 360 - param.c0();
2037 param.c2() = -param.c2();
2044 /* We only have to issue a warning if these atoms are perturbed! */
2046 for (j = 0; (j < nral); j++)
2048 bPert = bPert || PERTURBED(at->atom[param.a[j]]);
2053 sprintf(errbuf, "No default %s types for perturbed atoms, "
2054 "using normal values", interaction_function[ftype].longname);
2055 warning(wi, errbuf);
2061 if ((ftype == F_PDIHS || ftype == F_ANGRES || ftype == F_ANGRESZ)
2062 && param.c[5] != param.c[2])
2064 sprintf(errbuf, "%s multiplicity can not be perturbed %f!=%f",
2065 interaction_function[ftype].longname,
2066 param.c[2], param.c[5]);
2067 warning_error_and_exit(wi, errbuf, FARGS);
2070 if (IS_TABULATED(ftype) && param.c[2] != param.c[0])
2072 sprintf(errbuf, "%s table number can not be perturbed %d!=%d",
2073 interaction_function[ftype].longname,
2074 static_cast<int>(param.c[0]+0.5), static_cast<int>(param.c[2]+0.5));
2075 warning_error_and_exit(wi, errbuf, FARGS);
2078 /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
2079 if (ftype == F_RBDIHS)
2082 for (i = 0; i < NRFP(ftype); i++)
2084 if (param.c[i] != 0)
2095 /* Put the values in the appropriate arrays */
2096 add_param_to_list (&bond[ftype], ¶m);
2098 /* Push additional torsions from FF for ftype==9 if we have them.
2099 * We have already checked that the A/B states do not differ in this case,
2100 * so we do not have to double-check that again, or the vsite stuff.
2101 * In addition, those torsions cannot be automatically perturbed.
2103 if (bDef && ftype == F_PDIHS)
2105 for (i = 1; i < nparam_defA; i++)
2107 /* Advance pointer! */
2109 for (j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++)
2111 param.c[j] = param_defA->c[j];
2113 /* And push the next term for this torsion */
2114 add_param_to_list (&bond[ftype], ¶m);
2119 void push_cmap(directive d, t_params bondtype[], t_params bond[],
2120 t_atoms *at, gpp_atomtype_t atype, char *line,
2123 const char *aaformat[MAXATOMLIST+1] =
2134 int i, j, ftype, nral, nread, ncmap_params;
2136 int aa[MAXATOMLIST+1];
2137 char errbuf[STRLEN];
2141 ftype = ifunc_index(d, 1);
2145 nread = sscanf(line, aaformat[nral-1],
2146 &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
2153 else if (nread == nral)
2155 ftype = ifunc_index(d, 1);
2158 /* Check for double atoms and atoms out of bounds */
2159 for (i = 0; i < nral; i++)
2161 if (aa[i] < 1 || aa[i] > at->nr)
2163 sprintf(errbuf, "Atom index (%d) in %s out of bounds (1-%d).\n"
2164 "This probably means that you have inserted topology section \"%s\"\n"
2165 "in a part belonging to a different molecule than you intended to.\n"
2166 "In that case move the \"%s\" section to the right molecule.",
2167 aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
2168 warning_error_and_exit(wi, errbuf, FARGS);
2171 for (j = i+1; (j < nral); j++)
2175 sprintf(errbuf, "Duplicate atom index (%d) in %s", aa[i], dir2str(d));
2176 warning_error(wi, errbuf);
2181 /* default force parameters */
2182 for (j = 0; (j < MAXATOMLIST); j++)
2184 param.a[j] = aa[j]-1;
2186 for (j = 0; (j < MAXFORCEPARAM); j++)
2191 /* Get the cmap type for this cmap angle */
2192 bFound = default_cmap_params(bondtype, at, atype, ¶m, FALSE, &cmap_type, &ncmap_params, wi);
2194 /* We want exactly one parameter (the cmap type in state A (currently no state B) back */
2195 if (bFound && ncmap_params == 1)
2197 /* Put the values in the appropriate arrays */
2198 param.c[0] = cmap_type;
2199 add_param_to_list(&bond[ftype], ¶m);
2203 /* This is essentially the same check as in default_cmap_params() done one more time */
2204 sprintf(errbuf, "Unable to assign a cmap type to torsion %d %d %d %d and %d\n",
2205 param.a[0]+1, param.a[1]+1, param.a[2]+1, param.a[3]+1, param.a[4]+1);
2206 warning_error_and_exit(wi, errbuf, FARGS);
2212 void push_vsitesn(directive d, t_params bond[],
2213 t_atoms *at, char *line,
2217 int type, ftype, j, n, ret, nj, a;
2219 double *weight = nullptr, weight_tot;
2221 char errbuf[STRLEN];
2223 /* default force parameters */
2224 for (j = 0; (j < MAXATOMLIST); j++)
2226 param.a[j] = NOTSET;
2228 for (j = 0; (j < MAXFORCEPARAM); j++)
2234 ret = sscanf(ptr, "%d%n", &a, &n);
2238 sprintf(errbuf, "Expected an atom index in section \"%s\"",
2240 warning_error_and_exit(wi, errbuf, FARGS);
2245 sscanf(ptr, "%d%n", &type, &n);
2247 ftype = ifunc_index(d, type);
2253 ret = sscanf(ptr, "%d%n", &a, &n);
2260 srenew(weight, nj+20);
2269 /* Here we use the A-state mass as a parameter.
2270 * Note that the B-state mass has no influence.
2272 weight[nj] = at->atom[atc[nj]].m;
2276 ret = sscanf(ptr, "%lf%n", &(weight[nj]), &n);
2280 sprintf(errbuf, "No weight or negative weight found for vsiten constructing atom %d (atom index %d)",
2282 warning_error_and_exit(wi, errbuf, FARGS);
2286 sprintf(errbuf, "Unknown vsiten type %d", type);
2287 warning_error_and_exit(wi, errbuf, FARGS);
2289 weight_tot += weight[nj];
2297 sprintf(errbuf, "Expected more than one atom index in section \"%s\"",
2299 warning_error_and_exit(wi, errbuf, FARGS);
2302 if (weight_tot == 0)
2304 warning_error_and_exit(wi, "The total mass of the construting atoms is zero", FARGS);
2307 for (j = 0; j < nj; j++)
2309 param.a[1] = atc[j];
2311 param.c[1] = weight[j]/weight_tot;
2312 /* Put the values in the appropriate arrays */
2313 add_param_to_list (&bond[ftype], ¶m);
2320 void push_mol(int nrmols, t_molinfo mols[], char *pline, int *whichmol,
2325 char errbuf[STRLEN];
2327 if (sscanf(pline, "%s%d", type, nrcopies) != 2)
2333 /* Search moleculename.
2334 * Here we originally only did case insensitive matching. But because
2335 * some PDB files can have many chains and use case to generate more
2336 * chain-identifiers, which in turn end up in our moleculetype name,
2337 * we added support for case-sensitivity.
2343 for (int i = 0; i < nrmols; i++)
2345 if (strcmp(type, *(mols[i].name)) == 0)
2350 if (gmx_strcasecmp(type, *(mols[i].name)) == 0)
2359 // select the case sensitive match
2360 *whichmol = matchcs;
2364 // avoid matching case-insensitive when we have multiple matches
2367 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);
2368 warning_error_and_exit(wi, errbuf, FARGS);
2372 // select the unique case insensitive match
2373 *whichmol = matchci;
2377 sprintf(errbuf, "No such moleculetype %s", type);
2378 warning_error_and_exit(wi, errbuf, FARGS);
2383 void init_block2(t_block2 *b2, int natom)
2388 snew(b2->nra, b2->nr);
2389 snew(b2->a, b2->nr);
2390 for (i = 0; (i < b2->nr); i++)
2396 void done_block2(t_block2 *b2)
2402 for (i = 0; (i < b2->nr); i++)
2412 void push_excl(char *line, t_block2 *b2, warninp_t wi)
2416 char base[STRLEN], format[STRLEN];
2417 char errbuf[STRLEN];
2419 if (sscanf(line, "%d", &i) == 0)
2424 if ((1 <= i) && (i <= b2->nr))
2432 strcpy(base, "%*d");
2435 strcpy(format, base);
2436 strcat(format, "%d");
2437 n = sscanf(line, format, &j);
2440 if ((1 <= j) && (j <= b2->nr))
2443 srenew(b2->a[i], ++(b2->nra[i]));
2444 b2->a[i][b2->nra[i]-1] = j;
2445 /* also add the reverse exclusion! */
2446 srenew(b2->a[j], ++(b2->nra[j]));
2447 b2->a[j][b2->nra[j]-1] = i;
2448 strcat(base, "%*d");
2452 sprintf(errbuf, "Invalid Atomnr j: %d, b2->nr: %d\n", j, b2->nr);
2453 warning_error_and_exit(wi, errbuf, FARGS);
2460 void b_to_b2(t_blocka *b, t_block2 *b2)
2465 for (i = 0; (i < b->nr); i++)
2467 for (j = b->index[i]; (j < b->index[i+1]); j++)
2470 srenew(b2->a[i], ++b2->nra[i]);
2471 b2->a[i][b2->nra[i]-1] = a;
2476 void b2_to_b(t_block2 *b2, t_blocka *b)
2482 for (i = 0; (i < b2->nr); i++)
2485 for (j = 0; (j < b2->nra[i]); j++)
2487 b->a[nra+j] = b2->a[i][j];
2491 /* terminate list */
2495 void merge_excl(t_blocka *excl, t_block2 *b2, warninp_t wi)
2500 char errbuf[STRLEN];
2506 else if (b2->nr != excl->nr)
2508 sprintf(errbuf, "DEATH HORROR: b2->nr = %d, while excl->nr = %d",
2510 warning_error_and_exit(wi, errbuf, FARGS);
2513 /* First copy all entries from excl to b2 */
2516 /* Count and sort the exclusions */
2518 for (i = 0; (i < b2->nr); i++)
2522 /* remove double entries */
2523 std::sort(b2->a[i], b2->a[i]+b2->nra[i]);
2525 for (j = 1; (j < b2->nra[i]); j++)
2527 if (b2->a[i][j] != b2->a[i][k-1])
2529 b2->a[i][k] = b2->a[i][j];
2538 srenew(excl->a, excl->nra);
2543 int add_atomtype_decoupled(t_symtab *symtab, gpp_atomtype_t at,
2544 t_nbparam ***nbparam, t_nbparam ***pair)
2550 /* Define an atom type with all parameters set to zero (no interactions) */
2553 /* Type for decoupled atoms could be anything,
2554 * this should be changed automatically later when required.
2556 atom.ptype = eptAtom;
2557 for (i = 0; (i < MAXFORCEPARAM); i++)
2562 nr = add_atomtype(at, symtab, &atom, "decoupled", ¶m, -1, 0);
2564 /* Add space in the non-bonded parameters matrix */
2565 realloc_nb_params(at, nbparam, pair);
2570 static void convert_pairs_to_pairsQ(t_params *plist,
2571 real fudgeQQ, t_atoms *atoms)
2573 t_param *paramp1, *paramp2, *paramnew;
2574 int i, j, p1nr, p2nr, p2newnr;
2576 /* Add the pair list to the pairQ list */
2577 p1nr = plist[F_LJ14].nr;
2578 p2nr = plist[F_LJC14_Q].nr;
2579 p2newnr = p1nr + p2nr;
2580 snew(paramnew, p2newnr);
2582 paramp1 = plist[F_LJ14].param;
2583 paramp2 = plist[F_LJC14_Q].param;
2585 /* Fill in the new F_LJC14_Q array with the old one. NOTE:
2586 it may be possible to just ADD the converted F_LJ14 array
2587 to the old F_LJC14_Q array, but since we have to create
2588 a new sized memory structure, better just to deep copy it all.
2591 for (i = 0; i < p2nr; i++)
2593 /* Copy over parameters */
2594 for (j = 0; j < 5; j++) /* entries are 0-4 for F_LJC14_Q */
2596 paramnew[i].c[j] = paramp2[i].c[j];
2599 /* copy over atoms */
2600 for (j = 0; j < 2; j++)
2602 paramnew[i].a[j] = paramp2[i].a[j];
2606 for (i = p2nr; i < p2newnr; i++)
2610 /* Copy over parameters */
2611 paramnew[i].c[0] = fudgeQQ;
2612 paramnew[i].c[1] = atoms->atom[paramp1[j].a[0]].q;
2613 paramnew[i].c[2] = atoms->atom[paramp1[j].a[1]].q;
2614 paramnew[i].c[3] = paramp1[j].c[0];
2615 paramnew[i].c[4] = paramp1[j].c[1];
2617 /* copy over atoms */
2618 paramnew[i].a[0] = paramp1[j].a[0];
2619 paramnew[i].a[1] = paramp1[j].a[1];
2622 /* free the old pairlists */
2623 sfree(plist[F_LJC14_Q].param);
2624 sfree(plist[F_LJ14].param);
2626 /* now assign the new data to the F_LJC14_Q structure */
2627 plist[F_LJC14_Q].param = paramnew;
2628 plist[F_LJC14_Q].nr = p2newnr;
2630 /* Empty the LJ14 pairlist */
2631 plist[F_LJ14].nr = 0;
2632 plist[F_LJ14].param = nullptr;
2635 static void generate_LJCpairsNB(t_molinfo *mol, int nb_funct, t_params *nbp, warninp_t wi)
2637 int n, ntype, i, j, k;
2642 char errbuf[STRLEN];
2645 atom = mol->atoms.atom;
2647 ntype = static_cast<int>(std::sqrt(static_cast<double>(nbp->nr)));
2648 GMX_ASSERT(ntype * ntype == nbp->nr, "Number of pairs of generated non-bonded parameters should be a perfect square");
2650 for (i = 0; i < MAXATOMLIST; i++)
2652 param.a[i] = NOTSET;
2654 for (i = 0; i < MAXFORCEPARAM; i++)
2656 param.c[i] = NOTSET;
2659 /* Add a pair interaction for all non-excluded atom pairs */
2661 for (i = 0; i < n; i++)
2663 for (j = i+1; j < n; j++)
2666 for (k = excl->index[i]; k < excl->index[i+1]; k++)
2668 if (excl->a[k] == j)
2675 if (nb_funct != F_LJ)
2677 sprintf(errbuf, "Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones");
2678 warning_error_and_exit(wi, errbuf, FARGS);
2682 param.c[0] = atom[i].q;
2683 param.c[1] = atom[j].q;
2684 param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0];
2685 param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1];
2686 add_param_to_list(&mol->plist[F_LJC_PAIRS_NB], ¶m);
2692 static void set_excl_all(t_blocka *excl)
2696 /* Get rid of the current exclusions and exclude all atom pairs */
2698 excl->nra = nat*nat;
2699 srenew(excl->a, excl->nra);
2701 for (i = 0; i < nat; i++)
2704 for (j = 0; j < nat; j++)
2709 excl->index[nat] = k;
2712 static void decouple_atoms(t_atoms *atoms, int atomtype_decouple,
2713 int couple_lam0, int couple_lam1,
2714 const char *mol_name, warninp_t wi)
2717 char errbuf[STRLEN];
2719 for (i = 0; i < atoms->nr; i++)
2723 atom = &atoms->atom[i];
2725 if (atom->qB != atom->q || atom->typeB != atom->type)
2727 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.",
2728 i + 1, mol_name, "couple-moltype");
2729 warning_error_and_exit(wi, errbuf, FARGS);
2732 if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamVDW)
2736 if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamQ)
2738 atom->type = atomtype_decouple;
2740 if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamVDW)
2744 if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamQ)
2746 atom->typeB = atomtype_decouple;
2751 void convert_moltype_couple(t_molinfo *mol, int atomtype_decouple, real fudgeQQ,
2752 int couple_lam0, int couple_lam1,
2753 bool bCoupleIntra, int nb_funct, t_params *nbp,
2756 convert_pairs_to_pairsQ(mol->plist, fudgeQQ, &mol->atoms);
2759 generate_LJCpairsNB(mol, nb_funct, nbp, wi);
2760 set_excl_all(&mol->excls);
2762 decouple_atoms(&mol->atoms, atomtype_decouple, couple_lam0, couple_lam1,