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, 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.
44 #include "gromacs/legacyheaders/macros.h"
45 #include "gromacs/legacyheaders/names.h"
50 #include "gromacs/legacyheaders/warninp.h"
51 #include "gpp_atomtype.h"
52 #include "gpp_bond_atomtype.h"
54 #include "gromacs/topology/symtab.h"
55 #include "gromacs/utility/cstringutil.h"
56 #include "gromacs/utility/fatalerror.h"
57 #include "gromacs/utility/smalloc.h"
59 void generate_nbparams(int comb, int ftype, t_params *plist, gpp_atomtype_t atype,
64 real c, bi, bj, ci, cj, ci0, ci1, ci2, cj0, cj1, cj2;
67 /* Lean mean shortcuts */
68 nr = get_atomtype_ntypes(atype);
70 snew(plist->param, nr*nr);
73 /* Fill the matrix with force parameters */
81 for (i = k = 0; (i < nr); i++)
83 for (j = 0; (j < nr); j++, k++)
85 for (nf = 0; (nf < nrfp); nf++)
87 ci = get_atomtype_nbparam(i, nf, atype);
88 cj = get_atomtype_nbparam(j, nf, atype);
90 plist->param[k].c[nf] = c;
96 case eCOMB_ARITHMETIC:
97 /* c0 and c1 are sigma and epsilon */
98 for (i = k = 0; (i < nr); i++)
100 for (j = 0; (j < nr); j++, k++)
102 ci0 = get_atomtype_nbparam(i, 0, atype);
103 cj0 = get_atomtype_nbparam(j, 0, atype);
104 ci1 = get_atomtype_nbparam(i, 1, atype);
105 cj1 = get_atomtype_nbparam(j, 1, atype);
106 plist->param[k].c[0] = (fabs(ci0) + fabs(cj0))*0.5;
107 /* Negative sigma signals that c6 should be set to zero later,
108 * so we need to propagate that through the combination rules.
110 if (ci0 < 0 || cj0 < 0)
112 plist->param[k].c[0] *= -1;
114 plist->param[k].c[1] = sqrt(ci1*cj1);
119 case eCOMB_GEOM_SIG_EPS:
120 /* c0 and c1 are sigma and epsilon */
121 for (i = k = 0; (i < nr); i++)
123 for (j = 0; (j < nr); j++, k++)
125 ci0 = get_atomtype_nbparam(i, 0, atype);
126 cj0 = get_atomtype_nbparam(j, 0, atype);
127 ci1 = get_atomtype_nbparam(i, 1, atype);
128 cj1 = get_atomtype_nbparam(j, 1, atype);
129 plist->param[k].c[0] = sqrt(fabs(ci0*cj0));
130 /* Negative sigma signals that c6 should be set to zero later,
131 * so we need to propagate that through the combination rules.
133 if (ci0 < 0 || cj0 < 0)
135 plist->param[k].c[0] *= -1;
137 plist->param[k].c[1] = sqrt(ci1*cj1);
143 gmx_fatal(FARGS, "No such combination rule %d", comb);
147 gmx_incons("Topology processing, generate nb parameters");
152 /* Buckingham rules */
153 for (i = k = 0; (i < nr); i++)
155 for (j = 0; (j < nr); j++, k++)
157 ci0 = get_atomtype_nbparam(i, 0, atype);
158 cj0 = get_atomtype_nbparam(j, 0, atype);
159 ci2 = get_atomtype_nbparam(i, 2, atype);
160 cj2 = get_atomtype_nbparam(j, 2, atype);
161 bi = get_atomtype_nbparam(i, 1, atype);
162 bj = get_atomtype_nbparam(j, 1, atype);
163 plist->param[k].c[0] = sqrt(ci0 * cj0);
164 if ((bi == 0) || (bj == 0))
166 plist->param[k].c[1] = 0;
170 plist->param[k].c[1] = 2.0/(1/bi+1/bj);
172 plist->param[k].c[2] = sqrt(ci2 * cj2);
178 sprintf(errbuf, "Invalid nonbonded type %s",
179 interaction_function[ftype].longname);
180 warning_error(wi, errbuf);
184 static void realloc_nb_params(gpp_atomtype_t at,
185 t_nbparam ***nbparam, t_nbparam ***pair)
187 /* Add space in the non-bonded parameters matrix */
188 int atnr = get_atomtype_ntypes(at);
189 srenew(*nbparam, atnr);
190 snew((*nbparam)[atnr-1], atnr);
194 snew((*pair)[atnr-1], atnr);
198 static void copy_B_from_A(int ftype, double *c)
202 nrfpA = NRFPA(ftype);
203 nrfpB = NRFPB(ftype);
205 /* Copy the B parameters from the first nrfpB A parameters */
206 for (i = 0; (i < nrfpB); i++)
212 void push_at (t_symtab *symtab, gpp_atomtype_t at, t_bond_atomtype bat,
213 char *line, int nb_funct,
214 t_nbparam ***nbparam, t_nbparam ***pair,
221 t_xlate xl[eptNR] = {
229 int nr, i, nfields, j, pt, nfp0 = -1;
230 int batype_nr, nread;
231 char type[STRLEN], btype[STRLEN], ptype[STRLEN];
233 double c[MAXFORCEPARAM];
234 double radius, vol, surftens, gb_radius, S_hct;
235 char tmpfield[12][100]; /* Max 12 fields of width 100 */
240 gmx_bool have_atomic_number;
241 gmx_bool have_bonded_type;
246 /* First assign input line to temporary array */
247 nfields = sscanf(line, "%s%s%s%s%s%s%s%s%s%s%s%s",
248 tmpfield[0], tmpfield[1], tmpfield[2], tmpfield[3], tmpfield[4], tmpfield[5],
249 tmpfield[6], tmpfield[7], tmpfield[8], tmpfield[9], tmpfield[10], tmpfield[11]);
251 /* Comments on optional fields in the atomtypes section:
253 * The force field format is getting a bit old. For OPLS-AA we needed
254 * to add a special bonded atomtype, and for Gerrit Groenhofs QM/MM stuff
255 * we also needed the atomic numbers.
256 * To avoid making all old or user-generated force fields unusable we
257 * have introduced both these quantities as optional columns, and do some
258 * acrobatics to check whether they are present or not.
259 * This will all look much nicer when we switch to XML... sigh.
261 * Field 0 (mandatory) is the nonbonded type name. (string)
262 * Field 1 (optional) is the bonded type (string)
263 * Field 2 (optional) is the atomic number (int)
264 * Field 3 (mandatory) is the mass (numerical)
265 * Field 4 (mandatory) is the charge (numerical)
266 * Field 5 (mandatory) is the particle type (single character)
267 * This is followed by a number of nonbonded parameters.
269 * The safest way to identify the format is the particle type field.
271 * So, here is what we do:
273 * A. Read in the first six fields as strings
274 * B. If field 3 (starting from 0) is a single char, we have neither
275 * bonded_type or atomic numbers.
276 * C. If field 5 is a single char we have both.
277 * D. If field 4 is a single char we check field 1. If this begins with
278 * an alphabetical character we have bonded types, otherwise atomic numbers.
287 if ( (strlen(tmpfield[5]) == 1) && isalpha(tmpfield[5][0]) )
289 have_bonded_type = TRUE;
290 have_atomic_number = TRUE;
292 else if ( (strlen(tmpfield[3]) == 1) && isalpha(tmpfield[3][0]) )
294 have_bonded_type = FALSE;
295 have_atomic_number = FALSE;
299 have_bonded_type = ( isalpha(tmpfield[1][0]) != 0 );
300 have_atomic_number = !have_bonded_type;
303 /* optional fields */
317 if (have_atomic_number)
319 if (have_bonded_type)
321 nread = sscanf(line, "%s%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf",
322 type, btype, &atomnr, &m, &q, ptype, &c[0], &c[1],
323 &radius, &vol, &surftens, &gb_radius);
332 /* have_atomic_number && !have_bonded_type */
333 nread = sscanf(line, "%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf",
334 type, &atomnr, &m, &q, ptype, &c[0], &c[1],
335 &radius, &vol, &surftens, &gb_radius);
345 if (have_bonded_type)
347 /* !have_atomic_number && have_bonded_type */
348 nread = sscanf(line, "%s%s%lf%lf%s%lf%lf%lf%lf%lf%lf",
349 type, btype, &m, &q, ptype, &c[0], &c[1],
350 &radius, &vol, &surftens, &gb_radius);
359 /* !have_atomic_number && !have_bonded_type */
360 nread = sscanf(line, "%s%lf%lf%s%lf%lf%lf%lf%lf%lf",
361 type, &m, &q, ptype, &c[0], &c[1],
362 &radius, &vol, &surftens, &gb_radius);
371 if (!have_bonded_type)
376 if (!have_atomic_number)
386 if (have_atomic_number)
388 if (have_bonded_type)
390 nread = sscanf(line, "%s%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
391 type, btype, &atomnr, &m, &q, ptype, &c[0], &c[1], &c[2],
392 &radius, &vol, &surftens, &gb_radius);
401 /* have_atomic_number && !have_bonded_type */
402 nread = sscanf(line, "%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
403 type, &atomnr, &m, &q, ptype, &c[0], &c[1], &c[2],
404 &radius, &vol, &surftens, &gb_radius);
414 if (have_bonded_type)
416 /* !have_atomic_number && have_bonded_type */
417 nread = sscanf(line, "%s%s%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
418 type, btype, &m, &q, ptype, &c[0], &c[1], &c[2],
419 &radius, &vol, &surftens, &gb_radius);
428 /* !have_atomic_number && !have_bonded_type */
429 nread = sscanf(line, "%s%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
430 type, &m, &q, ptype, &c[0], &c[1], &c[2],
431 &radius, &vol, &surftens, &gb_radius);
440 if (!have_bonded_type)
445 if (!have_atomic_number)
453 gmx_fatal(FARGS, "Invalid function type %d in push_at %s %d", nb_funct,
456 for (j = nfp0; (j < MAXFORCEPARAM); j++)
461 if (strlen(type) == 1 && isdigit(type[0]))
463 gmx_fatal(FARGS, "Atom type names can't be single digits.");
466 if (strlen(btype) == 1 && isdigit(btype[0]))
468 gmx_fatal(FARGS, "Bond atom type names can't be single digits.");
471 /* Hack to read old topologies */
472 if (gmx_strcasecmp(ptype, "D") == 0)
476 for (j = 0; (j < eptNR); j++)
478 if (gmx_strcasecmp(ptype, xl[j].entry) == 0)
485 gmx_fatal(FARGS, "Invalid particle type %s on line %s",
491 fprintf(debug, "ptype: %s\n", ptype_str[pt]);
497 for (i = 0; (i < MAXFORCEPARAM); i++)
502 if ((batype_nr = get_bond_atomtype_type(btype, bat)) == NOTSET)
504 add_bond_atomtype(bat, symtab, btype);
506 batype_nr = get_bond_atomtype_type(btype, bat);
508 if ((nr = get_atomtype_type(type, at)) != NOTSET)
510 sprintf(errbuf, "Overriding atomtype %s", type);
512 if ((nr = set_atomtype(nr, at, symtab, atom, type, param, batype_nr,
513 radius, vol, surftens, atomnr, gb_radius, S_hct)) == NOTSET)
515 gmx_fatal(FARGS, "Replacing atomtype %s failed", type);
518 else if ((nr = add_atomtype(at, symtab, atom, type, param,
519 batype_nr, radius, vol,
520 surftens, atomnr, gb_radius, S_hct)) == NOTSET)
522 gmx_fatal(FARGS, "Adding atomtype %s failed", type);
526 /* Add space in the non-bonded parameters matrix */
527 realloc_nb_params(at, nbparam, pair);
533 static void push_bondtype(t_params * bt,
537 gmx_bool bAllowRepeat,
542 gmx_bool bTest, bFound, bCont, bId;
544 int nrfp = NRFP(ftype);
547 /* If bAllowRepeat is TRUE, we allow multiple entries as long as they
548 are on directly _adjacent_ lines.
551 /* First check if our atomtypes are _identical_ (not reversed) to the previous
552 entry. If they are not identical we search for earlier duplicates. If they are
553 we can skip it, since we already searched for the first line
560 if (bAllowRepeat && nr > 1)
562 for (j = 0, bCont = TRUE; (j < nral); j++)
564 bCont = bCont && (b->a[j] == bt->param[nr-2].a[j]);
568 /* Search for earlier duplicates if this entry was not a continuation
569 from the previous line.
574 for (i = 0; (i < nr); i++)
577 for (j = 0; (j < nral); j++)
579 bTest = (bTest && (b->a[j] == bt->param[i].a[j]));
584 for (j = 0; (j < nral); j++)
586 bTest = (bTest && (b->a[nral-1-j] == bt->param[i].a[j]));
594 for (j = 0; (j < nrfp); j++)
596 bId = bId && (bt->param[i].c[j] == b->c[j]);
600 sprintf(errbuf, "Overriding %s parameters.%s",
601 interaction_function[ftype].longname,
602 (ftype == F_PDIHS) ? "\nUse dihedraltype 4 to allow several multiplicity terms." : "");
604 fprintf(stderr, " old:");
605 for (j = 0; (j < nrfp); j++)
607 fprintf(stderr, " %g", bt->param[i].c[j]);
609 fprintf(stderr, " \n new: %s\n\n", line);
613 for (j = 0; (j < nrfp); j++)
615 bt->param[i].c[j] = b->c[j];
626 /* fill the arrays up and down */
627 memcpy(bt->param[bt->nr].c, b->c, sizeof(b->c));
628 memcpy(bt->param[bt->nr].a, b->a, sizeof(b->a));
629 memcpy(bt->param[bt->nr+1].c, b->c, sizeof(b->c));
631 /* The definitions of linear angles depend on the order of atoms,
632 * that means that for atoms i-j-k, with certain parameter a, the
633 * corresponding k-j-i angle will have parameter 1-a.
635 if (ftype == F_LINEAR_ANGLES)
637 bt->param[bt->nr+1].c[0] = 1-bt->param[bt->nr+1].c[0];
638 bt->param[bt->nr+1].c[2] = 1-bt->param[bt->nr+1].c[2];
641 for (j = 0; (j < nral); j++)
643 bt->param[bt->nr+1].a[j] = b->a[nral-1-j];
650 void push_bt(directive d, t_params bt[], int nral,
652 t_bond_atomtype bat, char *line,
655 const char *formal[MAXATOMLIST+1] = {
664 const char *formnl[MAXATOMLIST+1] = {
670 "%*s%*s%*s%*s%*s%*s",
671 "%*s%*s%*s%*s%*s%*s%*s"
673 const char *formlf = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
674 int i, ft, ftype, nn, nrfp, nrfpA, nrfpB;
676 char alc[MAXATOMLIST+1][20];
677 /* One force parameter more, so we can check if we read too many */
678 double c[MAXFORCEPARAM+1];
682 if ((bat && at) || (!bat && !at))
684 gmx_incons("You should pass either bat or at to push_bt");
687 /* Make format string (nral ints+functype) */
688 if ((nn = sscanf(line, formal[nral],
689 alc[0], alc[1], alc[2], alc[3], alc[4], alc[5])) != nral+1)
691 sprintf(errbuf, "Not enough atomtypes (%d instead of %d)", nn-1, nral);
692 warning_error(wi, errbuf);
696 ft = strtol(alc[nral], NULL, 10);
697 ftype = ifunc_index(d, ft);
699 nrfpA = interaction_function[ftype].nrfpA;
700 nrfpB = interaction_function[ftype].nrfpB;
701 strcpy(f1, formnl[nral]);
703 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]))
708 /* Copy the B-state from the A-state */
709 copy_B_from_A(ftype, c);
715 warning_error(wi, "Not enough parameters");
717 else if (nn > nrfpA && nn < nrfp)
719 warning_error(wi, "Too many parameters or not enough parameters for topology B");
723 warning_error(wi, "Too many parameters");
725 for (i = nn; (i < nrfp); i++)
731 for (i = 0; (i < nral); i++)
733 if (at && ((p.a[i] = get_atomtype_type(alc[i], at)) == NOTSET))
735 gmx_fatal(FARGS, "Unknown atomtype %s\n", alc[i]);
737 else if (bat && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET))
739 gmx_fatal(FARGS, "Unknown bond_atomtype %s\n", alc[i]);
742 for (i = 0; (i < nrfp); i++)
746 push_bondtype (&(bt[ftype]), &p, nral, ftype, FALSE, line, wi);
750 void push_dihedraltype(directive d, t_params bt[],
751 t_bond_atomtype bat, char *line,
754 const char *formal[MAXATOMLIST+1] = {
763 const char *formnl[MAXATOMLIST+1] = {
769 "%*s%*s%*s%*s%*s%*s",
770 "%*s%*s%*s%*s%*s%*s%*s"
772 const char *formlf[MAXFORCEPARAM] = {
778 "%lf%lf%lf%lf%lf%lf",
779 "%lf%lf%lf%lf%lf%lf%lf",
780 "%lf%lf%lf%lf%lf%lf%lf%lf",
781 "%lf%lf%lf%lf%lf%lf%lf%lf%lf",
782 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
783 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
784 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
786 int i, ft, ftype, nn, nrfp, nrfpA, nrfpB, nral;
788 char alc[MAXATOMLIST+1][20];
789 double c[MAXFORCEPARAM];
791 gmx_bool bAllowRepeat;
794 /* This routine accepts dihedraltypes defined from either 2 or 4 atoms.
796 * We first check for 2 atoms with the 3th column being an integer
797 * defining the type. If this isn't the case, we try it with 4 atoms
798 * and the 5th column defining the dihedral type.
800 nn = sscanf(line, formal[4], alc[0], alc[1], alc[2], alc[3], alc[4]);
801 if (nn >= 3 && strlen(alc[2]) == 1 && isdigit(alc[2][0]))
804 ft = strtol(alc[nral], NULL, 10);
805 /* Move atom types around a bit and use 'X' for wildcard atoms
806 * to create a 4-atom dihedral definition with arbitrary atoms in
809 if (alc[2][0] == '2')
811 /* improper - the two atomtypes are 1,4. Use wildcards for 2,3 */
812 strcpy(alc[3], alc[1]);
813 sprintf(alc[2], "X");
814 sprintf(alc[1], "X");
815 /* alc[0] stays put */
819 /* proper - the two atomtypes are 2,3. Use wildcards for 1,4 */
820 sprintf(alc[3], "X");
821 strcpy(alc[2], alc[1]);
822 strcpy(alc[1], alc[0]);
823 sprintf(alc[0], "X");
826 else if (nn == 5 && strlen(alc[4]) == 1 && isdigit(alc[4][0]))
829 ft = strtol(alc[nral], NULL, 10);
833 sprintf(errbuf, "Incorrect number of atomtypes for dihedral (%d instead of 2 or 4)", nn-1);
834 warning_error(wi, errbuf);
840 /* Previously, we have always overwritten parameters if e.g. a torsion
841 with the same atomtypes occurs on multiple lines. However, CHARMM and
842 some other force fields specify multiple dihedrals over some bonds,
843 including cosines with multiplicity 6 and somethimes even higher.
844 Thus, they cannot be represented with Ryckaert-Bellemans terms.
845 To add support for these force fields, Dihedral type 9 is identical to
846 normal proper dihedrals, but repeated entries are allowed.
853 bAllowRepeat = FALSE;
857 ftype = ifunc_index(d, ft);
859 nrfpA = interaction_function[ftype].nrfpA;
860 nrfpB = interaction_function[ftype].nrfpB;
862 strcpy(f1, formnl[nral]);
863 strcat(f1, formlf[nrfp-1]);
865 /* Check number of parameters given */
866 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]))
871 /* Copy the B-state from the A-state */
872 copy_B_from_A(ftype, c);
878 warning_error(wi, "Not enough parameters");
880 else if (nn > nrfpA && nn < nrfp)
882 warning_error(wi, "Too many parameters or not enough parameters for topology B");
886 warning_error(wi, "Too many parameters");
888 for (i = nn; (i < nrfp); i++)
895 for (i = 0; (i < 4); i++)
897 if (!strcmp(alc[i], "X"))
903 if ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET)
905 gmx_fatal(FARGS, "Unknown bond_atomtype %s", alc[i]);
909 for (i = 0; (i < nrfp); i++)
913 /* Always use 4 atoms here, since we created two wildcard atoms
914 * if there wasn't of them 4 already.
916 push_bondtype (&(bt[ftype]), &p, 4, ftype, bAllowRepeat, line, wi);
920 void push_nbt(directive d, t_nbparam **nbt, gpp_atomtype_t atype,
921 char *pline, int nb_funct,
925 const char *form2 = "%*s%*s%*s%lf%lf";
926 const char *form3 = "%*s%*s%*s%lf%lf%lf";
927 const char *form4 = "%*s%*s%*s%lf%lf%lf%lf";
928 const char *form5 = "%*s%*s%*s%lf%lf%lf%lf%lf";
930 int i, f, n, ftype, atnr, nrfp;
938 if (sscanf (pline, "%s%s%d", a0, a1, &f) != 3)
944 ftype = ifunc_index(d, f);
946 if (ftype != nb_funct)
948 sprintf(errbuf, "Trying to add %s while the default nonbond type is %s",
949 interaction_function[ftype].longname,
950 interaction_function[nb_funct].longname);
951 warning_error(wi, errbuf);
955 /* Get the force parameters */
959 n = sscanf(pline, form4, &c[0], &c[1], &c[2], &c[3]);
965 /* When the B topology parameters are not set,
966 * copy them from topology A
969 for (i = n; i < nrfp; i++)
974 else if (ftype == F_LJC14_Q)
976 n = sscanf(pline, form5, &c[0], &c[1], &c[2], &c[3], &dum);
979 incorrect_n_param(wi);
985 if (sscanf(pline, form3, &c[0], &c[1], &dum) != 2)
987 incorrect_n_param(wi);
993 if (sscanf(pline, form4, &c[0], &c[1], &c[2], &dum) != 3)
995 incorrect_n_param(wi);
1001 gmx_fatal(FARGS, "Number of force parameters for nonbonded interactions is %d"
1002 " in file %s, line %d", nrfp, __FILE__, __LINE__);
1004 for (i = 0; (i < nrfp); i++)
1009 /* Put the parameters in the matrix */
1010 if ((ai = get_atomtype_type (a0, atype)) == NOTSET)
1012 gmx_fatal(FARGS, "Atomtype %s not found", a0);
1014 if ((aj = get_atomtype_type (a1, atype)) == NOTSET)
1016 gmx_fatal(FARGS, "Atomtype %s not found", a1);
1018 nbp = &(nbt[max(ai, aj)][min(ai, aj)]);
1023 for (i = 0; i < nrfp; i++)
1025 bId = bId && (nbp->c[i] == cr[i]);
1029 sprintf(errbuf, "Overriding non-bonded parameters,");
1030 warning(wi, errbuf);
1031 fprintf(stderr, " old:");
1032 for (i = 0; i < nrfp; i++)
1034 fprintf(stderr, " %g", nbp->c[i]);
1036 fprintf(stderr, " new\n%s\n", pline);
1040 for (i = 0; i < nrfp; i++)
1047 push_gb_params (gpp_atomtype_t at, char *line,
1052 double radius, vol, surftens, gb_radius, S_hct;
1053 char atypename[STRLEN];
1054 char errbuf[STRLEN];
1056 if ( (nfield = sscanf(line, "%s%lf%lf%lf%lf%lf", atypename, &radius, &vol, &surftens, &gb_radius, &S_hct)) != 6)
1058 sprintf(errbuf, "Too few gb parameters for type %s\n", atypename);
1059 warning(wi, errbuf);
1062 /* Search for atomtype */
1063 atype = get_atomtype_type(atypename, at);
1065 if (atype == NOTSET)
1067 printf("Couldn't find topology match for atomtype %s\n", atypename);
1071 set_atomtype_gbparam(at, atype, radius, vol, surftens, gb_radius, S_hct);
1075 push_cmaptype(directive d, t_params bt[], int nral, gpp_atomtype_t at,
1076 t_bond_atomtype bat, char *line,
1079 const char *formal = "%s%s%s%s%s%s%s%s";
1081 int i, j, ft, ftype, nn, nrfp, nrfpA, nrfpB;
1083 int nxcmap, nycmap, ncmap, read_cmap, sl, nct;
1084 char s[20], alc[MAXATOMLIST+2][20];
1086 gmx_bool bAllowRepeat;
1089 /* Keep the compiler happy */
1093 if ((nn = sscanf(line, formal, alc[0], alc[1], alc[2], alc[3], alc[4], alc[5], alc[6], alc[7])) != nral+3)
1095 sprintf(errbuf, "Incorrect number of atomtypes for cmap (%d instead of 5)", nn-1);
1096 warning_error(wi, errbuf);
1100 /* Compute an offset for each line where the cmap parameters start
1101 * ie. where the atom types and grid spacing information ends
1103 for (i = 0; i < nn; i++)
1105 start += (int)strlen(alc[i]);
1108 /* There are nn-1 spaces between the atom types and the grid spacing info in the cmap.itp file */
1109 /* start is the position on the line where we start to read the actual cmap grid data from the itp file */
1110 start = start + nn -1;
1112 ft = strtol(alc[nral], NULL, 10);
1113 nxcmap = strtol(alc[nral+1], NULL, 10);
1114 nycmap = strtol(alc[nral+2], NULL, 10);
1116 /* Check for equal grid spacing in x and y dims */
1117 if (nxcmap != nycmap)
1119 gmx_fatal(FARGS, "Not the same grid spacing in x and y for cmap grid: x=%d, y=%d", nxcmap, nycmap);
1122 ncmap = nxcmap*nycmap;
1123 ftype = ifunc_index(d, ft);
1124 nrfpA = strtol(alc[6], NULL, 10)*strtol(alc[6], NULL, 10);
1125 nrfpB = strtol(alc[7], NULL, 10)*strtol(alc[7], NULL, 10);
1128 /* Allocate memory for the CMAP grid */
1129 bt[F_CMAP].ncmap += nrfp;
1130 srenew(bt[F_CMAP].cmap, bt[F_CMAP].ncmap);
1132 /* Read in CMAP parameters */
1134 for (i = 0; i < ncmap; i++)
1136 while (isspace(*(line+start+sl)))
1140 nn = sscanf(line+start+sl, " %s ", s);
1142 bt[F_CMAP].cmap[i+(bt[F_CMAP].ncmap)-nrfp] = strtod(s, NULL);
1150 gmx_fatal(FARGS, "Error in reading cmap parameter for angle %s %s %s %s %s", alc[0], alc[1], alc[2], alc[3], alc[4]);
1155 /* Check do that we got the number of parameters we expected */
1156 if (read_cmap == nrfpA)
1158 for (i = 0; i < ncmap; i++)
1160 bt[F_CMAP].cmap[i+ncmap] = bt[F_CMAP].cmap[i];
1165 if (read_cmap < nrfpA)
1167 warning_error(wi, "Not enough cmap parameters");
1169 else if (read_cmap > nrfpA && read_cmap < nrfp)
1171 warning_error(wi, "Too many cmap parameters or not enough parameters for topology B");
1173 else if (read_cmap > nrfp)
1175 warning_error(wi, "Too many cmap parameters");
1180 /* Set grid spacing and the number of grids (we assume these numbers to be the same for all grids
1181 * so we can safely assign them each time
1183 bt[F_CMAP].grid_spacing = nxcmap; /* Or nycmap, they need to be equal */
1184 bt[F_CMAP].nc = bt[F_CMAP].nc + 1; /* Since we are incrementing here, we need to subtract later, see (*****) */
1185 nct = (nral+1) * bt[F_CMAP].nc;
1187 /* Allocate memory for the cmap_types information */
1188 srenew(bt[F_CMAP].cmap_types, nct);
1190 for (i = 0; (i < nral); i++)
1192 if (at && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET))
1194 gmx_fatal(FARGS, "Unknown atomtype %s\n", alc[i]);
1196 else if (bat && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET))
1198 gmx_fatal(FARGS, "Unknown bond_atomtype %s\n", alc[i]);
1201 /* Assign a grid number to each cmap_type */
1202 bt[F_CMAP].cmap_types[bt[F_CMAP].nct++] = get_bond_atomtype_type(alc[i], bat);
1205 /* Assign a type number to this cmap */
1206 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 (****) */
1208 /* Check for the correct number of atoms (again) */
1209 if (bt[F_CMAP].nct != nct)
1211 gmx_fatal(FARGS, "Incorrect number of atom types (%d) in cmap type %d\n", nct, bt[F_CMAP].nc);
1214 /* Is this correct?? */
1215 for (i = 0; i < MAXFORCEPARAM; i++)
1220 /* Push the bond to the bondlist */
1221 push_bondtype (&(bt[ftype]), &p, nral, ftype, FALSE, line, wi);
1225 static void push_atom_now(t_symtab *symtab, t_atoms *at, int atomnr,
1227 int type, char *ctype, int ptype,
1229 char *resname, char *name, real m0, real q0,
1230 int typeB, char *ctypeB, real mB, real qB)
1232 int j, resind = 0, resnr;
1236 if (((nr == 0) && (atomnr != 1)) || (nr && (atomnr != at->nr+1)))
1238 gmx_fatal(FARGS, "Atoms in the .top are not numbered consecutively from 1 (rather, atomnr = %d, while at->nr = %d)", atomnr, at->nr);
1241 j = strlen(resnumberic) - 1;
1242 if (isdigit(resnumberic[j]))
1248 ric = resnumberic[j];
1249 if (j == 0 || !isdigit(resnumberic[j-1]))
1251 gmx_fatal(FARGS, "Invalid residue number '%s' for atom %d",
1252 resnumberic, atomnr);
1255 resnr = strtol(resnumberic, NULL, 10);
1259 resind = at->atom[nr-1].resind;
1261 if (nr == 0 || strcmp(resname, *at->resinfo[resind].name) != 0 ||
1262 resnr != at->resinfo[resind].nr ||
1263 ric != at->resinfo[resind].ic)
1273 at->nres = resind + 1;
1274 srenew(at->resinfo, at->nres);
1275 at->resinfo[resind].name = put_symtab(symtab, resname);
1276 at->resinfo[resind].nr = resnr;
1277 at->resinfo[resind].ic = ric;
1281 resind = at->atom[at->nr-1].resind;
1284 /* New atom instance
1285 * get new space for arrays
1287 srenew(at->atom, nr+1);
1288 srenew(at->atomname, nr+1);
1289 srenew(at->atomtype, nr+1);
1290 srenew(at->atomtypeB, nr+1);
1293 at->atom[nr].type = type;
1294 at->atom[nr].ptype = ptype;
1295 at->atom[nr].q = q0;
1296 at->atom[nr].m = m0;
1297 at->atom[nr].typeB = typeB;
1298 at->atom[nr].qB = qB;
1299 at->atom[nr].mB = mB;
1301 at->atom[nr].resind = resind;
1302 at->atom[nr].atomnumber = atomicnumber;
1303 at->atomname[nr] = put_symtab(symtab, name);
1304 at->atomtype[nr] = put_symtab(symtab, ctype);
1305 at->atomtypeB[nr] = put_symtab(symtab, ctypeB);
1309 void push_cg(t_block *block, int *lastindex, int index, int a)
1313 fprintf (debug, "Index %d, Atom %d\n", index, a);
1316 if (((block->nr) && (*lastindex != index)) || (!block->nr))
1318 /* add a new block */
1320 srenew(block->index, block->nr+1);
1322 block->index[block->nr] = a + 1;
1326 void push_atom(t_symtab *symtab, t_block *cgs,
1327 t_atoms *at, gpp_atomtype_t atype, char *line, int *lastcg,
1331 int cgnumber, atomnr, type, typeB, nscan;
1332 char id[STRLEN], ctype[STRLEN], ctypeB[STRLEN],
1333 resnumberic[STRLEN], resname[STRLEN], name[STRLEN], check[STRLEN];
1334 double m, q, mb, qb;
1335 real m0, q0, mB, qB;
1337 /* Make a shortcut for writing in this molecule */
1340 /* Fixed parameters */
1341 if (sscanf(line, "%s%s%s%s%s%d",
1342 id, ctype, resnumberic, resname, name, &cgnumber) != 6)
1347 sscanf(id, "%d", &atomnr);
1348 if ((type = get_atomtype_type(ctype, atype)) == NOTSET)
1350 gmx_fatal(FARGS, "Atomtype %s not found", ctype);
1352 ptype = get_atomtype_ptype(type, atype);
1354 /* Set default from type */
1355 q0 = get_atomtype_qA(type, atype);
1356 m0 = get_atomtype_massA(type, atype);
1361 /* Optional parameters */
1362 nscan = sscanf(line, "%*s%*s%*s%*s%*s%*s%lf%lf%s%lf%lf%s",
1363 &q, &m, ctypeB, &qb, &mb, check);
1365 /* Nasty switch that falls thru all the way down! */
1374 if ((typeB = get_atomtype_type(ctypeB, atype)) == NOTSET)
1376 gmx_fatal(FARGS, "Atomtype %s not found", ctypeB);
1378 qB = get_atomtype_qA(typeB, atype);
1379 mB = get_atomtype_massA(typeB, atype);
1388 warning_error(wi, "Too many parameters");
1397 fprintf(debug, "mB=%g, qB=%g, typeB=%d\n", mB, qB, typeB);
1400 push_cg(cgs, lastcg, cgnumber, nr);
1402 push_atom_now(symtab, at, atomnr, get_atomtype_atomnumber(type, atype),
1403 type, ctype, ptype, resnumberic,
1404 resname, name, m0, q0, typeB,
1405 typeB == type ? ctype : ctypeB, mB, qB);
1408 void push_molt(t_symtab *symtab, int *nmol, t_molinfo **mol, char *line,
1415 if ((sscanf(line, "%s%d", type, &nrexcl)) != 2)
1417 warning_error(wi, "Expected a molecule type name and nrexcl");
1420 /* Test if this atomtype overwrites another */
1424 if (gmx_strcasecmp(*((*mol)[i].name), type) == 0)
1426 gmx_fatal(FARGS, "moleculetype %s is redefined", type);
1432 srenew(*mol, *nmol);
1433 newmol = &((*mol)[*nmol-1]);
1434 init_molinfo(newmol);
1436 /* Fill in the values */
1437 newmol->name = put_symtab(symtab, type);
1438 newmol->nrexcl = nrexcl;
1439 newmol->excl_set = FALSE;
1442 static gmx_bool default_nb_params(int ftype, t_params bt[], t_atoms *at,
1443 t_param *p, int c_start, gmx_bool bB, gmx_bool bGenPairs)
1445 int i, j, ti, tj, ntype;
1448 int nr = bt[ftype].nr;
1449 int nral = NRAL(ftype);
1450 int nrfp = interaction_function[ftype].nrfpA;
1451 int nrfpB = interaction_function[ftype].nrfpB;
1453 if ((!bB && nrfp == 0) || (bB && nrfpB == 0))
1461 /* First test the generated-pair position to save
1462 * time when we have 1000*1000 entries for e.g. OPLS...
1467 ti = at->atom[p->a[0]].typeB;
1468 tj = at->atom[p->a[1]].typeB;
1472 ti = at->atom[p->a[0]].type;
1473 tj = at->atom[p->a[1]].type;
1475 pi = &(bt[ftype].param[ntype*ti+tj]);
1476 bFound = ((ti == pi->a[0]) && (tj == pi->a[1]));
1479 /* Search explicitly if we didnt find it */
1482 for (i = 0; ((i < nr) && !bFound); i++)
1484 pi = &(bt[ftype].param[i]);
1487 for (j = 0; ((j < nral) &&
1488 (at->atom[p->a[j]].typeB == pi->a[j])); j++)
1495 for (j = 0; ((j < nral) &&
1496 (at->atom[p->a[j]].type == pi->a[j])); j++)
1501 bFound = (j == nral);
1509 if (nrfp+nrfpB > MAXFORCEPARAM)
1511 gmx_incons("Too many force parameters");
1513 for (j = c_start; (j < nrfpB); j++)
1515 p->c[nrfp+j] = pi->c[j];
1520 for (j = c_start; (j < nrfp); j++)
1528 for (j = c_start; (j < nrfp); j++)
1536 static gmx_bool default_cmap_params(t_params bondtype[],
1537 t_atoms *at, gpp_atomtype_t atype,
1538 t_param *p, gmx_bool bB,
1539 int *cmap_type, int *nparam_def)
1541 int i, j, nparam_found;
1543 gmx_bool bFound = FALSE;
1548 /* Match the current cmap angle against the list of cmap_types */
1549 for (i = 0; i < bondtype[F_CMAP].nct && !bFound; i += 6)
1558 (get_atomtype_batype(at->atom[p->a[0]].type, atype) == bondtype[F_CMAP].cmap_types[i]) &&
1559 (get_atomtype_batype(at->atom[p->a[1]].type, atype) == bondtype[F_CMAP].cmap_types[i+1]) &&
1560 (get_atomtype_batype(at->atom[p->a[2]].type, atype) == bondtype[F_CMAP].cmap_types[i+2]) &&
1561 (get_atomtype_batype(at->atom[p->a[3]].type, atype) == bondtype[F_CMAP].cmap_types[i+3]) &&
1562 (get_atomtype_batype(at->atom[p->a[4]].type, atype) == bondtype[F_CMAP].cmap_types[i+4]))
1564 /* Found cmap torsion */
1566 ct = bondtype[F_CMAP].cmap_types[i+5];
1572 /* If we did not find a matching type for this cmap torsion */
1575 gmx_fatal(FARGS, "Unknown cmap torsion between atoms %d %d %d %d %d\n",
1576 p->a[0]+1, p->a[1]+1, p->a[2]+1, p->a[3]+1, p->a[4]+1);
1579 *nparam_def = nparam_found;
1585 static gmx_bool default_params(int ftype, t_params bt[],
1586 t_atoms *at, gpp_atomtype_t atype,
1587 t_param *p, gmx_bool bB,
1588 t_param **param_def,
1591 int i, j, nparam_found;
1592 gmx_bool bFound, bSame;
1595 int nr = bt[ftype].nr;
1596 int nral = NRAL(ftype);
1597 int nrfpA = interaction_function[ftype].nrfpA;
1598 int nrfpB = interaction_function[ftype].nrfpB;
1600 if ((!bB && nrfpA == 0) || (bB && nrfpB == 0))
1606 /* We allow wildcards now. The first type (with or without wildcards) that
1607 * fits is used, so you should probably put the wildcarded bondtypes
1608 * at the end of each section.
1612 /* OPLS uses 1000s of dihedraltypes, so in order to speed up the scanning we have a
1613 * special case for this. Check for B state outside loop to speed it up.
1615 if (ftype == F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS || ftype == F_PIDIHS)
1619 for (i = 0; ((i < nr) && !bFound); i++)
1621 pi = &(bt[ftype].param[i]);
1624 ((pi->AI == -1) || (get_atomtype_batype(at->atom[p->AI].typeB, atype) == pi->AI)) &&
1625 ((pi->AJ == -1) || (get_atomtype_batype(at->atom[p->AJ].typeB, atype) == pi->AJ)) &&
1626 ((pi->AK == -1) || (get_atomtype_batype(at->atom[p->AK].typeB, atype) == pi->AK)) &&
1627 ((pi->AL == -1) || (get_atomtype_batype(at->atom[p->AL].typeB, atype) == pi->AL))
1634 for (i = 0; ((i < nr) && !bFound); i++)
1636 pi = &(bt[ftype].param[i]);
1639 ((pi->AI == -1) || (get_atomtype_batype(at->atom[p->AI].type, atype) == pi->AI)) &&
1640 ((pi->AJ == -1) || (get_atomtype_batype(at->atom[p->AJ].type, atype) == pi->AJ)) &&
1641 ((pi->AK == -1) || (get_atomtype_batype(at->atom[p->AK].type, atype) == pi->AK)) &&
1642 ((pi->AL == -1) || (get_atomtype_batype(at->atom[p->AL].type, atype) == pi->AL))
1646 /* Find additional matches for this dihedral - necessary for ftype==9 which is used e.g. for charmm.
1647 * The rules in that case is that additional matches HAVE to be on adjacent lines!
1653 /* Continue from current i value */
1654 for (j = i+1; j < nr && bSame; j += 2)
1656 pj = &(bt[ftype].param[j]);
1657 bSame = (pi->AI == pj->AI && pi->AJ == pj->AJ && pi->AK == pj->AK && pi->AL == pj->AL);
1662 /* nparam_found will be increased as long as the numbers match */
1666 else /* Not a dihedral */
1668 for (i = 0; ((i < nr) && !bFound); i++)
1670 pi = &(bt[ftype].param[i]);
1673 for (j = 0; ((j < nral) &&
1674 (get_atomtype_batype(at->atom[p->a[j]].typeB, atype) == pi->a[j])); j++)
1681 for (j = 0; ((j < nral) &&
1682 (get_atomtype_batype(at->atom[p->a[j]].type, atype) == pi->a[j])); j++)
1687 bFound = (j == nral);
1696 *nparam_def = nparam_found;
1703 void push_bond(directive d, t_params bondtype[], t_params bond[],
1704 t_atoms *at, gpp_atomtype_t atype, char *line,
1705 gmx_bool bBonded, gmx_bool bGenPairs, real fudgeQQ,
1706 gmx_bool bZero, gmx_bool *bWarn_copy_A_B,
1709 const char *aaformat[MAXATOMLIST] = {
1717 const char *asformat[MAXATOMLIST] = {
1722 "%*s%*s%*s%*s%*s%*s",
1723 "%*s%*s%*s%*s%*s%*s%*s"
1725 const char *ccformat = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
1726 int nr, i, j, nral, nral_fmt, nread, ftype;
1727 char format[STRLEN];
1728 /* One force parameter more, so we can check if we read too many */
1729 double cc[MAXFORCEPARAM+1];
1730 int aa[MAXATOMLIST+1];
1731 t_param param, paramB, *param_defA, *param_defB;
1732 gmx_bool bFoundA = FALSE, bFoundB = FALSE, bDef, bPert, bSwapParity = FALSE;
1733 int nparam_defA, nparam_defB;
1736 nparam_defA = nparam_defB = 0;
1738 ftype = ifunc_index(d, 1);
1740 for (j = 0; j < MAXATOMLIST; j++)
1744 bDef = (NRFP(ftype) > 0);
1746 if (ftype == F_SETTLE)
1748 /* SETTLE acts on 3 atoms, but the topology format only specifies
1749 * the first atom (for historical reasons).
1758 nread = sscanf(line, aaformat[nral_fmt-1],
1759 &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
1761 if (ftype == F_SETTLE)
1768 if (nread < nral_fmt)
1773 else if (nread > nral_fmt)
1775 /* this is a hack to allow for virtual sites with swapped parity */
1776 bSwapParity = (aa[nral] < 0);
1779 aa[nral] = -aa[nral];
1781 ftype = ifunc_index(d, aa[nral]);
1790 gmx_fatal(FARGS, "Negative function types only allowed for %s and %s",
1791 interaction_function[F_VSITE3FAD].longname,
1792 interaction_function[F_VSITE3OUT].longname);
1798 /* Check for double atoms and atoms out of bounds */
1799 for (i = 0; (i < nral); i++)
1801 if (aa[i] < 1 || aa[i] > at->nr)
1803 gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
1804 "Atom index (%d) in %s out of bounds (1-%d).\n"
1805 "This probably means that you have inserted topology section \"%s\"\n"
1806 "in a part belonging to a different molecule than you intended to.\n"
1807 "In that case move the \"%s\" section to the right molecule.",
1808 get_warning_file(wi), get_warning_line(wi),
1809 aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
1811 for (j = i+1; (j < nral); j++)
1815 sprintf(errbuf, "Duplicate atom index (%d) in %s", aa[i], dir2str(d));
1816 warning(wi, errbuf);
1821 /* default force parameters */
1822 for (j = 0; (j < MAXATOMLIST); j++)
1824 param.a[j] = aa[j]-1;
1826 for (j = 0; (j < MAXFORCEPARAM); j++)
1831 /* Get force params for normal and free energy perturbation
1832 * studies, as determined by types!
1837 bFoundA = default_params(ftype, bondtype, at, atype, ¶m, FALSE, ¶m_defA, &nparam_defA);
1840 /* Copy the A-state and B-state default parameters. */
1841 assert(NRFPA(ftype)+NRFPB(ftype) <= MAXFORCEPARAM);
1842 for (j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++)
1844 param.c[j] = param_defA->c[j];
1847 bFoundB = default_params(ftype, bondtype, at, atype, ¶m, TRUE, ¶m_defB, &nparam_defB);
1850 /* Copy only the B-state default parameters */
1851 for (j = NRFPA(ftype); (j < NRFP(ftype)); j++)
1853 param.c[j] = param_defB->c[j];
1857 else if (ftype == F_LJ14)
1859 bFoundA = default_nb_params(ftype, bondtype, at, ¶m, 0, FALSE, bGenPairs);
1860 bFoundB = default_nb_params(ftype, bondtype, at, ¶m, 0, TRUE, bGenPairs);
1862 else if (ftype == F_LJC14_Q)
1864 param.c[0] = fudgeQQ;
1865 /* Fill in the A-state charges as default parameters */
1866 param.c[1] = at->atom[param.a[0]].q;
1867 param.c[2] = at->atom[param.a[1]].q;
1868 /* The default LJ parameters are the standard 1-4 parameters */
1869 bFoundA = default_nb_params(F_LJ14, bondtype, at, ¶m, 3, FALSE, bGenPairs);
1872 else if (ftype == F_LJC_PAIRS_NB)
1874 /* Defaults are not supported here */
1880 gmx_incons("Unknown function type in push_bond");
1883 if (nread > nral_fmt)
1885 /* Manually specified parameters - in this case we discard multiple torsion info! */
1887 strcpy(format, asformat[nral_fmt-1]);
1888 strcat(format, ccformat);
1890 nread = sscanf(line, format, &cc[0], &cc[1], &cc[2], &cc[3], &cc[4], &cc[5],
1891 &cc[6], &cc[7], &cc[8], &cc[9], &cc[10], &cc[11], &cc[12]);
1893 if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0))
1895 /* We only have to issue a warning if these atoms are perturbed! */
1897 for (j = 0; (j < nral); j++)
1899 bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1902 if (bPert && *bWarn_copy_A_B)
1905 "Some parameters for bonded interaction involving perturbed atoms are specified explicitly in state A, but not B - copying A to B");
1906 warning(wi, errbuf);
1907 *bWarn_copy_A_B = FALSE;
1910 /* If only the A parameters were specified, copy them to the B state */
1911 /* The B-state parameters correspond to the first nrfpB
1912 * A-state parameters.
1914 for (j = 0; (j < NRFPB(ftype)); j++)
1916 cc[nread++] = cc[j];
1920 /* If nread was 0 or EOF, no parameters were read => use defaults.
1921 * If nread was nrfpA we copied above so nread=nrfp.
1922 * If nread was nrfp we are cool.
1923 * For F_LJC14_Q we allow supplying fudgeQQ only.
1924 * Anything else is an error!
1926 if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) &&
1927 !(ftype == F_LJC14_Q && nread == 1))
1929 gmx_fatal(FARGS, "Incorrect number of parameters - found %d, expected %d or %d for %s.",
1930 nread, NRFPA(ftype), NRFP(ftype),
1931 interaction_function[ftype].longname);
1934 for (j = 0; (j < nread); j++)
1939 /* Check whether we have to use the defaults */
1940 if (nread == NRFP(ftype))
1949 /* nread now holds the number of force parameters read! */
1954 /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */
1955 if (ftype == F_PDIHS)
1957 if ((nparam_defA != nparam_defB) || ((nparam_defA > 1 || nparam_defB > 1) && (param_defA != param_defB)))
1960 "Cannot automatically perturb a torsion with multiple terms to different form.\n"
1961 "Please specify perturbed parameters manually for this torsion in your topology!");
1962 warning_error(wi, errbuf);
1966 if (nread > 0 && nread < NRFPA(ftype))
1968 /* Issue an error, do not use defaults */
1969 sprintf(errbuf, "Not enough parameters, there should be at least %d (or 0 for defaults)", NRFPA(ftype));
1970 warning_error(wi, errbuf);
1973 if (nread == 0 || nread == EOF)
1977 if (interaction_function[ftype].flags & IF_VSITE)
1979 /* set them to NOTSET, will be calculated later */
1980 for (j = 0; (j < MAXFORCEPARAM); j++)
1982 param.c[j] = NOTSET;
1987 param.C1 = -1; /* flag to swap parity of vsite construction */
1994 fprintf(stderr, "NOTE: No default %s types, using zeroes\n",
1995 interaction_function[ftype].longname);
1999 sprintf(errbuf, "No default %s types", interaction_function[ftype].longname);
2000 warning_error(wi, errbuf);
2011 param.C0 = 360 - param.C0;
2014 param.C2 = -param.C2;
2021 /* We only have to issue a warning if these atoms are perturbed! */
2023 for (j = 0; (j < nral); j++)
2025 bPert = bPert || PERTURBED(at->atom[param.a[j]]);
2030 sprintf(errbuf, "No default %s types for perturbed atoms, "
2031 "using normal values", interaction_function[ftype].longname);
2032 warning(wi, errbuf);
2038 if ((ftype == F_PDIHS || ftype == F_ANGRES || ftype == F_ANGRESZ)
2039 && param.c[5] != param.c[2])
2041 gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2042 " %s multiplicity can not be perturbed %f!=%f",
2043 get_warning_file(wi), get_warning_line(wi),
2044 interaction_function[ftype].longname,
2045 param.c[2], param.c[5]);
2048 if (IS_TABULATED(ftype) && param.c[2] != param.c[0])
2050 gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2051 " %s table number can not be perturbed %d!=%d",
2052 get_warning_file(wi), get_warning_line(wi),
2053 interaction_function[ftype].longname,
2054 (int)(param.c[0]+0.5), (int)(param.c[2]+0.5));
2057 /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
2058 if (ftype == F_RBDIHS)
2061 for (i = 0; i < NRFP(ftype); i++)
2063 if (param.c[i] != 0)
2074 /* Put the values in the appropriate arrays */
2075 add_param_to_list (&bond[ftype], ¶m);
2077 /* Push additional torsions from FF for ftype==9 if we have them.
2078 * We have already checked that the A/B states do not differ in this case,
2079 * so we do not have to double-check that again, or the vsite stuff.
2080 * In addition, those torsions cannot be automatically perturbed.
2082 if (bDef && ftype == F_PDIHS)
2084 for (i = 1; i < nparam_defA; i++)
2086 /* Advance pointer! */
2088 for (j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++)
2090 param.c[j] = param_defA->c[j];
2092 /* And push the next term for this torsion */
2093 add_param_to_list (&bond[ftype], ¶m);
2098 void push_cmap(directive d, t_params bondtype[], t_params bond[],
2099 t_atoms *at, gpp_atomtype_t atype, char *line,
2102 const char *aaformat[MAXATOMLIST+1] =
2113 int i, j, nr, ftype, nral, nread, ncmap_params;
2115 int aa[MAXATOMLIST+1];
2118 t_param param, paramB, *param_defA, *param_defB;
2120 ftype = ifunc_index(d, 1);
2122 nr = bondtype[ftype].nr;
2125 nread = sscanf(line, aaformat[nral-1],
2126 &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]);
2133 else if (nread == nral)
2135 ftype = ifunc_index(d, 1);
2138 /* Check for double atoms and atoms out of bounds */
2139 for (i = 0; i < nral; i++)
2141 if (aa[i] < 1 || aa[i] > at->nr)
2143 gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2144 "Atom index (%d) in %s out of bounds (1-%d).\n"
2145 "This probably means that you have inserted topology section \"%s\"\n"
2146 "in a part belonging to a different molecule than you intended to.\n"
2147 "In that case move the \"%s\" section to the right molecule.",
2148 get_warning_file(wi), get_warning_line(wi),
2149 aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d));
2152 for (j = i+1; (j < nral); j++)
2156 sprintf(errbuf, "Duplicate atom index (%d) in %s", aa[i], dir2str(d));
2157 warning(wi, errbuf);
2162 /* default force parameters */
2163 for (j = 0; (j < MAXATOMLIST); j++)
2165 param.a[j] = aa[j]-1;
2167 for (j = 0; (j < MAXFORCEPARAM); j++)
2172 /* Get the cmap type for this cmap angle */
2173 bFound = default_cmap_params(bondtype, at, atype, ¶m, FALSE, &cmap_type, &ncmap_params);
2175 /* We want exactly one parameter (the cmap type in state A (currently no state B) back */
2176 if (bFound && ncmap_params == 1)
2178 /* Put the values in the appropriate arrays */
2179 param.c[0] = cmap_type;
2180 add_param_to_list(&bond[ftype], ¶m);
2184 /* This is essentially the same check as in default_cmap_params() done one more time */
2185 gmx_fatal(FARGS, "Unable to assign a cmap type to torsion %d %d %d %d and %d\n",
2186 param.a[0]+1, param.a[1]+1, param.a[2]+1, param.a[3]+1, param.a[4]+1);
2192 void push_vsitesn(directive d, t_params bond[],
2193 t_atoms *at, char *line,
2197 int type, ftype, j, n, ret, nj, a;
2199 double *weight = NULL, weight_tot;
2202 /* default force parameters */
2203 for (j = 0; (j < MAXATOMLIST); j++)
2205 param.a[j] = NOTSET;
2207 for (j = 0; (j < MAXFORCEPARAM); j++)
2213 ret = sscanf(ptr, "%d%n", &a, &n);
2217 gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2218 " Expected an atom index in section \"%s\"",
2219 get_warning_file(wi), get_warning_line(wi),
2225 ret = sscanf(ptr, "%d%n", &type, &n);
2227 ftype = ifunc_index(d, type);
2233 ret = sscanf(ptr, "%d%n", &a, &n);
2240 srenew(weight, nj+20);
2249 /* Here we use the A-state mass as a parameter.
2250 * Note that the B-state mass has no influence.
2252 weight[nj] = at->atom[atc[nj]].m;
2256 ret = sscanf(ptr, "%lf%n", &(weight[nj]), &n);
2260 gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2261 " No weight or negative weight found for vsiten constructing atom %d (atom index %d)",
2262 get_warning_file(wi), get_warning_line(wi),
2267 gmx_fatal(FARGS, "Unknown vsiten type %d", type);
2269 weight_tot += weight[nj];
2277 gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2278 " Expected more than one atom index in section \"%s\"",
2279 get_warning_file(wi), get_warning_line(wi),
2283 if (weight_tot == 0)
2285 gmx_fatal(FARGS, "[ file %s, line %d ]:\n"
2286 " The total mass of the construting atoms is zero",
2287 get_warning_file(wi), get_warning_line(wi));
2290 for (j = 0; j < nj; j++)
2292 param.a[1] = atc[j];
2294 param.c[1] = weight[j]/weight_tot;
2295 /* Put the values in the appropriate arrays */
2296 add_param_to_list (&bond[ftype], ¶m);
2303 void push_mol(int nrmols, t_molinfo mols[], char *pline, int *whichmol,
2311 if (sscanf(pline, "%s%d", type, &copies) != 2)
2317 /* search moleculename */
2318 for (i = 0; ((i < nrmols) && gmx_strcasecmp(type, *(mols[i].name))); i++)
2330 gmx_fatal(FARGS, "No such moleculetype %s", type);
2334 void init_block2(t_block2 *b2, int natom)
2339 snew(b2->nra, b2->nr);
2340 snew(b2->a, b2->nr);
2341 for (i = 0; (i < b2->nr); i++)
2347 void done_block2(t_block2 *b2)
2353 for (i = 0; (i < b2->nr); i++)
2363 void push_excl(char *line, t_block2 *b2)
2367 char base[STRLEN], format[STRLEN];
2369 if (sscanf(line, "%d", &i) == 0)
2374 if ((1 <= i) && (i <= b2->nr))
2382 fprintf(debug, "Unbound atom %d\n", i-1);
2386 strcpy(base, "%*d");
2389 strcpy(format, base);
2390 strcat(format, "%d");
2391 n = sscanf(line, format, &j);
2394 if ((1 <= j) && (j <= b2->nr))
2397 srenew(b2->a[i], ++(b2->nra[i]));
2398 b2->a[i][b2->nra[i]-1] = j;
2399 /* also add the reverse exclusion! */
2400 srenew(b2->a[j], ++(b2->nra[j]));
2401 b2->a[j][b2->nra[j]-1] = i;
2402 strcat(base, "%*d");
2406 gmx_fatal(FARGS, "Invalid Atomnr j: %d, b2->nr: %d\n", j, b2->nr);
2413 void b_to_b2(t_blocka *b, t_block2 *b2)
2418 for (i = 0; (i < b->nr); i++)
2420 for (j = b->index[i]; (j < b->index[i+1]); j++)
2423 srenew(b2->a[i], ++b2->nra[i]);
2424 b2->a[i][b2->nra[i]-1] = a;
2429 void b2_to_b(t_block2 *b2, t_blocka *b)
2435 for (i = 0; (i < b2->nr); i++)
2438 for (j = 0; (j < b2->nra[i]); j++)
2440 b->a[nra+j] = b2->a[i][j];
2444 /* terminate list */
2448 static int icomp(const void *v1, const void *v2)
2450 return (*((atom_id *) v1))-(*((atom_id *) v2));
2453 void merge_excl(t_blocka *excl, t_block2 *b2)
2463 else if (b2->nr != excl->nr)
2465 gmx_fatal(FARGS, "DEATH HORROR: b2->nr = %d, while excl->nr = %d",
2470 fprintf(debug, "Entering merge_excl\n");
2473 /* First copy all entries from excl to b2 */
2476 /* Count and sort the exclusions */
2478 for (i = 0; (i < b2->nr); i++)
2482 /* remove double entries */
2483 qsort(b2->a[i], (size_t)b2->nra[i], (size_t)sizeof(b2->a[i][0]), icomp);
2485 for (j = 1; (j < b2->nra[i]); j++)
2487 if (b2->a[i][j] != b2->a[i][k-1])
2489 b2->a[i][k] = b2->a[i][j];
2498 srenew(excl->a, excl->nra);
2503 int add_atomtype_decoupled(t_symtab *symtab, gpp_atomtype_t at,
2504 t_nbparam ***nbparam, t_nbparam ***pair)
2510 /* Define an atom type with all parameters set to zero (no interactions) */
2513 /* Type for decoupled atoms could be anything,
2514 * this should be changed automatically later when required.
2516 atom.ptype = eptAtom;
2517 for (i = 0; (i < MAXFORCEPARAM); i++)
2522 nr = add_atomtype(at, symtab, &atom, "decoupled", ¶m, -1, 0.0, 0.0, 0.0, 0, 0, 0);
2524 /* Add space in the non-bonded parameters matrix */
2525 realloc_nb_params(at, nbparam, pair);
2530 static void convert_pairs_to_pairsQ(t_params *plist,
2531 real fudgeQQ, t_atoms *atoms)
2533 t_param *paramp1, *paramp2, *paramnew;
2534 int i, j, p1nr, p2nr, p2newnr;
2536 /* Add the pair list to the pairQ list */
2537 p1nr = plist[F_LJ14].nr;
2538 p2nr = plist[F_LJC14_Q].nr;
2539 p2newnr = p1nr + p2nr;
2540 snew(paramnew, p2newnr);
2542 paramp1 = plist[F_LJ14].param;
2543 paramp2 = plist[F_LJC14_Q].param;
2545 /* Fill in the new F_LJC14_Q array with the old one. NOTE:
2546 it may be possible to just ADD the converted F_LJ14 array
2547 to the old F_LJC14_Q array, but since we have to create
2548 a new sized memory structure, better just to deep copy it all.
2551 for (i = 0; i < p2nr; i++)
2553 /* Copy over parameters */
2554 for (j = 0; j < 5; j++) /* entries are 0-4 for F_LJC14_Q */
2556 paramnew[i].c[j] = paramp2[i].c[j];
2559 /* copy over atoms */
2560 for (j = 0; j < 2; j++)
2562 paramnew[i].a[j] = paramp2[i].a[j];
2566 for (i = p2nr; i < p2newnr; i++)
2570 /* Copy over parameters */
2571 paramnew[i].c[0] = fudgeQQ;
2572 paramnew[i].c[1] = atoms->atom[paramp1[j].a[0]].q;
2573 paramnew[i].c[2] = atoms->atom[paramp1[j].a[1]].q;
2574 paramnew[i].c[3] = paramp1[j].c[0];
2575 paramnew[i].c[4] = paramp1[j].c[1];
2577 /* copy over atoms */
2578 paramnew[i].a[0] = paramp1[j].a[0];
2579 paramnew[i].a[1] = paramp1[j].a[1];
2582 /* free the old pairlists */
2583 sfree(plist[F_LJC14_Q].param);
2584 sfree(plist[F_LJ14].param);
2586 /* now assign the new data to the F_LJC14_Q structure */
2587 plist[F_LJC14_Q].param = paramnew;
2588 plist[F_LJC14_Q].nr = p2newnr;
2590 /* Empty the LJ14 pairlist */
2591 plist[F_LJ14].nr = 0;
2592 plist[F_LJ14].param = NULL;
2595 static void generate_LJCpairsNB(t_molinfo *mol, int nb_funct, t_params *nbp)
2597 int n, ntype, i, j, k;
2604 atom = mol->atoms.atom;
2606 ntype = sqrt(nbp->nr);
2608 for (i = 0; i < MAXATOMLIST; i++)
2610 param.a[i] = NOTSET;
2612 for (i = 0; i < MAXFORCEPARAM; i++)
2614 param.c[i] = NOTSET;
2617 /* Add a pair interaction for all non-excluded atom pairs */
2619 for (i = 0; i < n; i++)
2621 for (j = i+1; j < n; j++)
2624 for (k = excl->index[i]; k < excl->index[i+1]; k++)
2626 if (excl->a[k] == j)
2633 if (nb_funct != F_LJ)
2635 gmx_fatal(FARGS, "Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones");
2639 param.c[0] = atom[i].q;
2640 param.c[1] = atom[j].q;
2641 param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0];
2642 param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1];
2643 add_param_to_list(&mol->plist[F_LJC_PAIRS_NB], ¶m);
2649 static void set_excl_all(t_blocka *excl)
2653 /* Get rid of the current exclusions and exclude all atom pairs */
2655 excl->nra = nat*nat;
2656 srenew(excl->a, excl->nra);
2658 for (i = 0; i < nat; i++)
2661 for (j = 0; j < nat; j++)
2666 excl->index[nat] = k;
2669 static void decouple_atoms(t_atoms *atoms, int atomtype_decouple,
2670 int couple_lam0, int couple_lam1)
2674 for (i = 0; i < atoms->nr; i++)
2676 if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamVDW)
2678 atoms->atom[i].q = 0.0;
2680 if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamQ)
2682 atoms->atom[i].type = atomtype_decouple;
2684 if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamVDW)
2686 atoms->atom[i].qB = 0.0;
2688 if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamQ)
2690 atoms->atom[i].typeB = atomtype_decouple;
2695 void convert_moltype_couple(t_molinfo *mol, int atomtype_decouple, real fudgeQQ,
2696 int couple_lam0, int couple_lam1,
2697 gmx_bool bCoupleIntra, int nb_funct, t_params *nbp)
2699 convert_pairs_to_pairsQ(mol->plist, fudgeQQ, &mol->atoms);
2702 generate_LJCpairsNB(mol, nb_funct, nbp);
2703 set_excl_all(&mol->excls);
2705 decouple_atoms(&mol->atoms, atomtype_decouple, couple_lam0, couple_lam1);