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/math/functions.h"
57 #include "gromacs/mdtypes/md_enums.h"
58 #include "gromacs/topology/exclusionblocks.h"
59 #include "gromacs/topology/ifunc.h"
60 #include "gromacs/topology/symtab.h"
61 #include "gromacs/utility/arrayref.h"
62 #include "gromacs/utility/cstringutil.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/gmxassert.h"
65 #include "gromacs/utility/smalloc.h"
66 #include "gromacs/utility/stringutil.h"
68 void generate_nbparams(int comb
, int ftype
, t_params
*plist
, gpp_atomtype_t atype
,
73 real c
, bi
, bj
, ci
, cj
, ci0
, ci1
, ci2
, cj0
, cj1
, cj2
;
75 /* Lean mean shortcuts */
76 nr
= get_atomtype_ntypes(atype
);
78 snew(plist
->param
, nr
*nr
);
81 /* Fill the matrix with force parameters */
89 for (i
= k
= 0; (i
< nr
); i
++)
91 for (j
= 0; (j
< nr
); j
++, k
++)
93 for (nf
= 0; (nf
< nrfp
); nf
++)
95 ci
= get_atomtype_nbparam(i
, nf
, atype
);
96 cj
= get_atomtype_nbparam(j
, nf
, atype
);
97 c
= std::sqrt(ci
* cj
);
98 plist
->param
[k
].c
[nf
] = c
;
104 case eCOMB_ARITHMETIC
:
105 /* c0 and c1 are sigma and epsilon */
106 for (i
= k
= 0; (i
< nr
); i
++)
108 for (j
= 0; (j
< nr
); j
++, k
++)
110 ci0
= get_atomtype_nbparam(i
, 0, atype
);
111 cj0
= get_atomtype_nbparam(j
, 0, atype
);
112 ci1
= get_atomtype_nbparam(i
, 1, atype
);
113 cj1
= get_atomtype_nbparam(j
, 1, atype
);
114 plist
->param
[k
].c
[0] = (fabs(ci0
) + fabs(cj0
))*0.5;
115 /* Negative sigma signals that c6 should be set to zero later,
116 * so we need to propagate that through the combination rules.
118 if (ci0
< 0 || cj0
< 0)
120 plist
->param
[k
].c
[0] *= -1;
122 plist
->param
[k
].c
[1] = std::sqrt(ci1
*cj1
);
127 case eCOMB_GEOM_SIG_EPS
:
128 /* c0 and c1 are sigma and epsilon */
129 for (i
= k
= 0; (i
< nr
); i
++)
131 for (j
= 0; (j
< nr
); j
++, k
++)
133 ci0
= get_atomtype_nbparam(i
, 0, atype
);
134 cj0
= get_atomtype_nbparam(j
, 0, atype
);
135 ci1
= get_atomtype_nbparam(i
, 1, atype
);
136 cj1
= get_atomtype_nbparam(j
, 1, atype
);
137 plist
->param
[k
].c
[0] = std::sqrt(std::fabs(ci0
*cj0
));
138 /* Negative sigma signals that c6 should be set to zero later,
139 * so we need to propagate that through the combination rules.
141 if (ci0
< 0 || cj0
< 0)
143 plist
->param
[k
].c
[0] *= -1;
145 plist
->param
[k
].c
[1] = std::sqrt(ci1
*cj1
);
151 auto message
= gmx::formatString("No such combination rule %d", comb
);
152 warning_error_and_exit(wi
, message
, FARGS
);
156 gmx_incons("Topology processing, generate nb parameters");
161 /* Buckingham rules */
162 for (i
= k
= 0; (i
< nr
); i
++)
164 for (j
= 0; (j
< nr
); j
++, k
++)
166 ci0
= get_atomtype_nbparam(i
, 0, atype
);
167 cj0
= get_atomtype_nbparam(j
, 0, atype
);
168 ci2
= get_atomtype_nbparam(i
, 2, atype
);
169 cj2
= get_atomtype_nbparam(j
, 2, atype
);
170 bi
= get_atomtype_nbparam(i
, 1, atype
);
171 bj
= get_atomtype_nbparam(j
, 1, atype
);
172 plist
->param
[k
].c
[0] = std::sqrt(ci0
* cj0
);
173 if ((bi
== 0) || (bj
== 0))
175 plist
->param
[k
].c
[1] = 0;
179 plist
->param
[k
].c
[1] = 2.0/(1/bi
+1/bj
);
181 plist
->param
[k
].c
[2] = std::sqrt(ci2
* cj2
);
187 auto message
= gmx::formatString("Invalid nonbonded type %s",
188 interaction_function
[ftype
].longname
);
189 warning_error(wi
, message
);
193 static void realloc_nb_params(gpp_atomtype_t at
,
194 t_nbparam
***nbparam
, t_nbparam
***pair
)
196 /* Add space in the non-bonded parameters matrix */
197 int atnr
= get_atomtype_ntypes(at
);
198 srenew(*nbparam
, atnr
);
199 snew((*nbparam
)[atnr
-1], atnr
);
203 snew((*pair
)[atnr
-1], atnr
);
207 static void copy_B_from_A(int ftype
, double *c
)
211 nrfpA
= NRFPA(ftype
);
212 nrfpB
= NRFPB(ftype
);
214 /* Copy the B parameters from the first nrfpB A parameters */
215 for (i
= 0; (i
< nrfpB
); i
++)
221 void push_at (t_symtab
*symtab
, gpp_atomtype_t at
, t_bond_atomtype bat
,
222 char *line
, int nb_funct
,
223 t_nbparam
***nbparam
, t_nbparam
***pair
,
230 t_xlate xl
[eptNR
] = {
238 int nr
, i
, nfields
, j
, pt
, nfp0
= -1;
239 int batype_nr
, nread
;
240 char type
[STRLEN
], btype
[STRLEN
], ptype
[STRLEN
];
242 double c
[MAXFORCEPARAM
];
243 char tmpfield
[12][100]; /* Max 12 fields of width 100 */
247 bool have_atomic_number
;
248 bool have_bonded_type
;
253 /* First assign input line to temporary array */
254 nfields
= sscanf(line
, "%s%s%s%s%s%s%s%s%s%s%s%s",
255 tmpfield
[0], tmpfield
[1], tmpfield
[2], tmpfield
[3], tmpfield
[4], tmpfield
[5],
256 tmpfield
[6], tmpfield
[7], tmpfield
[8], tmpfield
[9], tmpfield
[10], tmpfield
[11]);
258 /* Comments on optional fields in the atomtypes section:
260 * The force field format is getting a bit old. For OPLS-AA we needed
261 * to add a special bonded atomtype, and for Gerrit Groenhofs QM/MM stuff
262 * we also needed the atomic numbers.
263 * To avoid making all old or user-generated force fields unusable we
264 * have introduced both these quantities as optional columns, and do some
265 * acrobatics to check whether they are present or not.
266 * This will all look much nicer when we switch to XML... sigh.
268 * Field 0 (mandatory) is the nonbonded type name. (string)
269 * Field 1 (optional) is the bonded type (string)
270 * Field 2 (optional) is the atomic number (int)
271 * Field 3 (mandatory) is the mass (numerical)
272 * Field 4 (mandatory) is the charge (numerical)
273 * Field 5 (mandatory) is the particle type (single character)
274 * This is followed by a number of nonbonded parameters.
276 * The safest way to identify the format is the particle type field.
278 * So, here is what we do:
280 * A. Read in the first six fields as strings
281 * B. If field 3 (starting from 0) is a single char, we have neither
282 * bonded_type or atomic numbers.
283 * C. If field 5 is a single char we have both.
284 * D. If field 4 is a single char we check field 1. If this begins with
285 * an alphabetical character we have bonded types, otherwise atomic numbers.
294 if ( (strlen(tmpfield
[5]) == 1) && isalpha(tmpfield
[5][0]) )
296 have_bonded_type
= TRUE
;
297 have_atomic_number
= TRUE
;
299 else if ( (strlen(tmpfield
[3]) == 1) && isalpha(tmpfield
[3][0]) )
301 have_bonded_type
= FALSE
;
302 have_atomic_number
= FALSE
;
306 have_bonded_type
= ( isalpha(tmpfield
[1][0]) != 0 );
307 have_atomic_number
= !have_bonded_type
;
310 /* optional fields */
319 if (have_atomic_number
)
321 if (have_bonded_type
)
323 nread
= sscanf(line
, "%s%s%d%lf%lf%s%lf%lf",
324 type
, btype
, &atomnr
, &m
, &q
, ptype
, &c
[0], &c
[1]);
333 /* have_atomic_number && !have_bonded_type */
334 nread
= sscanf(line
, "%s%d%lf%lf%s%lf%lf",
335 type
, &atomnr
, &m
, &q
, ptype
, &c
[0], &c
[1]);
345 if (have_bonded_type
)
347 /* !have_atomic_number && have_bonded_type */
348 nread
= sscanf(line
, "%s%s%lf%lf%s%lf%lf",
349 type
, btype
, &m
, &q
, ptype
, &c
[0], &c
[1]);
358 /* !have_atomic_number && !have_bonded_type */
359 nread
= sscanf(line
, "%s%lf%lf%s%lf%lf",
360 type
, &m
, &q
, ptype
, &c
[0], &c
[1]);
369 if (!have_bonded_type
)
374 if (!have_atomic_number
)
384 if (have_atomic_number
)
386 if (have_bonded_type
)
388 nread
= sscanf(line
, "%s%s%d%lf%lf%s%lf%lf%lf",
389 type
, btype
, &atomnr
, &m
, &q
, ptype
, &c
[0], &c
[1], &c
[2]);
398 /* have_atomic_number && !have_bonded_type */
399 nread
= sscanf(line
, "%s%d%lf%lf%s%lf%lf%lf",
400 type
, &atomnr
, &m
, &q
, ptype
, &c
[0], &c
[1], &c
[2]);
410 if (have_bonded_type
)
412 /* !have_atomic_number && have_bonded_type */
413 nread
= sscanf(line
, "%s%s%lf%lf%s%lf%lf%lf",
414 type
, btype
, &m
, &q
, ptype
, &c
[0], &c
[1], &c
[2]);
423 /* !have_atomic_number && !have_bonded_type */
424 nread
= sscanf(line
, "%s%lf%lf%s%lf%lf%lf",
425 type
, &m
, &q
, ptype
, &c
[0], &c
[1], &c
[2]);
434 if (!have_bonded_type
)
439 if (!have_atomic_number
)
447 auto message
= gmx::formatString("Invalid function type %d in push_at", nb_funct
);
448 warning_error_and_exit(wi
, message
, FARGS
);
450 for (j
= nfp0
; (j
< MAXFORCEPARAM
); j
++)
455 if (strlen(type
) == 1 && isdigit(type
[0]))
457 warning_error_and_exit(wi
, "Atom type names can't be single digits.", FARGS
);
460 if (strlen(btype
) == 1 && isdigit(btype
[0]))
462 warning_error_and_exit(wi
, "Bond atom type names can't be single digits.", FARGS
);
465 /* Hack to read old topologies */
466 if (gmx_strcasecmp(ptype
, "D") == 0)
470 for (j
= 0; (j
< eptNR
); j
++)
472 if (gmx_strcasecmp(ptype
, xl
[j
].entry
) == 0)
479 auto message
= gmx::formatString("Invalid particle type %s", ptype
);
480 warning_error_and_exit(wi
, message
, FARGS
);
487 for (i
= 0; (i
< MAXFORCEPARAM
); i
++)
492 if ((batype_nr
= get_bond_atomtype_type(btype
, bat
)) == NOTSET
)
494 add_bond_atomtype(bat
, symtab
, btype
);
496 batype_nr
= get_bond_atomtype_type(btype
, bat
);
498 if ((nr
= get_atomtype_type(type
, at
)) != NOTSET
)
500 auto message
= gmx::formatString("Overriding atomtype %s", type
);
501 warning(wi
, message
);
502 if ((nr
= set_atomtype(nr
, at
, symtab
, atom
, type
, param
, batype_nr
,
505 auto message
= gmx::formatString("Replacing atomtype %s failed", type
);
506 warning_error_and_exit(wi
, message
, FARGS
);
509 else if ((add_atomtype(at
, symtab
, atom
, type
, param
,
510 batype_nr
, atomnr
)) == NOTSET
)
512 auto message
= gmx::formatString("Adding atomtype %s failed", type
);
513 warning_error_and_exit(wi
, message
, FARGS
);
517 /* Add space in the non-bonded parameters matrix */
518 realloc_nb_params(at
, nbparam
, pair
);
524 //! Return whether the contents of \c a and \c b are the same, considering also reversed order.
525 template <typename T
>
526 static bool equalEitherForwardOrBackward(gmx::ArrayRef
<const T
> a
, gmx::ArrayRef
<const T
> b
)
528 return (std::equal(a
.begin(), a
.end(), b
.begin()) ||
529 std::equal(a
.begin(), a
.end(), b
.rbegin()));
532 static void push_bondtype(t_params
* bt
,
541 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 "
597 "type 9 for the same atoms, with either different parameters "
598 "and/or the first block has multiple lines. This is not supported.");
602 else if (!haveWarned
)
604 auto message
= gmx::formatString
605 ("Overriding %s parameters.%s",
606 interaction_function
[ftype
].longname
,
608 "\nUse dihedraltype 9 to allow several multiplicity terms. Only consecutive "
609 "lines are combined. Non-consective lines overwrite each other."
611 warning(wi
, message
);
613 fprintf(stderr
, " old: ");
614 for (int j
= 0; (j
< nrfp
); j
++)
616 fprintf(stderr
, " %g", bt
->param
[i
].c
[j
]);
618 fprintf(stderr
, " \n new: %s\n\n", line
);
624 if (!identicalParameters
&& !bAllowRepeat
)
626 /* Overwrite the parameters with the latest ones */
627 // TODO considering improving the following code by replacing with:
628 // std::copy(b->c, b->c + nrfp, bt->param[i].c);
629 for (int j
= 0; (j
< nrfp
); j
++)
631 bt
->param
[i
].c
[j
] = b
->c
[j
];
642 /* fill the arrays up and down */
643 memcpy(bt
->param
[bt
->nr
].c
, b
->c
, sizeof(b
->c
));
644 memcpy(bt
->param
[bt
->nr
].a
, b
->a
, sizeof(b
->a
));
645 memcpy(bt
->param
[bt
->nr
+1].c
, b
->c
, sizeof(b
->c
));
647 /* The definitions of linear angles depend on the order of atoms,
648 * that means that for atoms i-j-k, with certain parameter a, the
649 * corresponding k-j-i angle will have parameter 1-a.
651 if (ftype
== F_LINEAR_ANGLES
)
653 bt
->param
[bt
->nr
+1].c
[0] = 1-bt
->param
[bt
->nr
+1].c
[0];
654 bt
->param
[bt
->nr
+1].c
[2] = 1-bt
->param
[bt
->nr
+1].c
[2];
657 for (int j
= 0; (j
< nral
); j
++)
659 bt
->param
[bt
->nr
+1].a
[j
] = b
->a
[nral
-1-j
];
666 void push_bt(directive d
, t_params bt
[], int nral
,
668 t_bond_atomtype bat
, char *line
,
671 const char *formal
[MAXATOMLIST
+1] = {
680 const char *formnl
[MAXATOMLIST
+1] = {
686 "%*s%*s%*s%*s%*s%*s",
687 "%*s%*s%*s%*s%*s%*s%*s"
689 const char *formlf
= "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
690 int i
, ft
, ftype
, nn
, nrfp
, nrfpA
;
692 char alc
[MAXATOMLIST
+1][20];
693 /* One force parameter more, so we can check if we read too many */
694 double c
[MAXFORCEPARAM
+1];
697 if ((bat
&& at
) || (!bat
&& !at
))
699 gmx_incons("You should pass either bat or at to push_bt");
702 /* Make format string (nral ints+functype) */
703 if ((nn
= sscanf(line
, formal
[nral
],
704 alc
[0], alc
[1], alc
[2], alc
[3], alc
[4], alc
[5])) != nral
+1)
706 auto message
= gmx::formatString("Not enough atomtypes (%d instead of %d)", nn
-1, nral
);
707 warning_error(wi
, message
);
711 ft
= strtol(alc
[nral
], nullptr, 10);
712 ftype
= ifunc_index(d
, ft
);
714 nrfpA
= interaction_function
[ftype
].nrfpA
;
715 strcpy(f1
, formnl
[nral
]);
717 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]))
722 /* Copy the B-state from the A-state */
723 copy_B_from_A(ftype
, c
);
729 warning_error(wi
, "Not enough parameters");
731 else if (nn
> nrfpA
&& nn
< nrfp
)
733 warning_error(wi
, "Too many parameters or not enough parameters for topology B");
737 warning_error(wi
, "Too many parameters");
739 for (i
= nn
; (i
< nrfp
); i
++)
745 for (i
= 0; (i
< nral
); i
++)
747 if (at
&& ((p
.a
[i
] = get_atomtype_type(alc
[i
], at
)) == NOTSET
))
749 auto message
= gmx::formatString("Unknown atomtype %s\n", alc
[i
]);
750 warning_error_and_exit(wi
, message
, FARGS
);
752 else if (bat
&& ((p
.a
[i
] = get_bond_atomtype_type(alc
[i
], bat
)) == NOTSET
))
754 auto message
= gmx::formatString("Unknown bond_atomtype %s\n", alc
[i
]);
755 warning_error_and_exit(wi
, message
, FARGS
);
758 for (i
= 0; (i
< nrfp
); i
++)
762 push_bondtype (&(bt
[ftype
]), &p
, nral
, ftype
, FALSE
, line
, wi
);
766 void push_dihedraltype(directive d
, t_params bt
[],
767 t_bond_atomtype bat
, char *line
,
770 const char *formal
[MAXATOMLIST
+1] = {
779 const char *formnl
[MAXATOMLIST
+1] = {
785 "%*s%*s%*s%*s%*s%*s",
786 "%*s%*s%*s%*s%*s%*s%*s"
788 const char *formlf
[MAXFORCEPARAM
] = {
794 "%lf%lf%lf%lf%lf%lf",
795 "%lf%lf%lf%lf%lf%lf%lf",
796 "%lf%lf%lf%lf%lf%lf%lf%lf",
797 "%lf%lf%lf%lf%lf%lf%lf%lf%lf",
798 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
799 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
800 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
802 int i
, ft
, ftype
, nn
, nrfp
, nrfpA
, nral
;
804 char alc
[MAXATOMLIST
+1][20];
805 double c
[MAXFORCEPARAM
];
809 /* This routine accepts dihedraltypes defined from either 2 or 4 atoms.
811 * We first check for 2 atoms with the 3th column being an integer
812 * defining the type. If this isn't the case, we try it with 4 atoms
813 * and the 5th column defining the dihedral type.
815 nn
= sscanf(line
, formal
[4], alc
[0], alc
[1], alc
[2], alc
[3], alc
[4]);
816 if (nn
>= 3 && strlen(alc
[2]) == 1 && isdigit(alc
[2][0]))
819 ft
= strtol(alc
[nral
], nullptr, 10);
820 /* Move atom types around a bit and use 'X' for wildcard atoms
821 * to create a 4-atom dihedral definition with arbitrary atoms in
824 if (alc
[2][0] == '2')
826 /* improper - the two atomtypes are 1,4. Use wildcards for 2,3 */
827 strcpy(alc
[3], alc
[1]);
828 sprintf(alc
[2], "X");
829 sprintf(alc
[1], "X");
830 /* alc[0] stays put */
834 /* proper - the two atomtypes are 2,3. Use wildcards for 1,4 */
835 sprintf(alc
[3], "X");
836 strcpy(alc
[2], alc
[1]);
837 strcpy(alc
[1], alc
[0]);
838 sprintf(alc
[0], "X");
841 else if (nn
== 5 && strlen(alc
[4]) == 1 && isdigit(alc
[4][0]))
844 ft
= strtol(alc
[nral
], nullptr, 10);
848 auto message
= gmx::formatString("Incorrect number of atomtypes for dihedral (%d instead of 2 or 4)", nn
-1);
849 warning_error(wi
, message
);
855 /* Previously, we have always overwritten parameters if e.g. a torsion
856 with the same atomtypes occurs on multiple lines. However, CHARMM and
857 some other force fields specify multiple dihedrals over some bonds,
858 including cosines with multiplicity 6 and somethimes even higher.
859 Thus, they cannot be represented with Ryckaert-Bellemans terms.
860 To add support for these force fields, Dihedral type 9 is identical to
861 normal proper dihedrals, but repeated entries are allowed.
868 bAllowRepeat
= FALSE
;
872 ftype
= ifunc_index(d
, ft
);
874 nrfpA
= interaction_function
[ftype
].nrfpA
;
876 strcpy(f1
, formnl
[nral
]);
877 strcat(f1
, formlf
[nrfp
-1]);
879 /* Check number of parameters given */
880 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]))
885 /* Copy the B-state from the A-state */
886 copy_B_from_A(ftype
, c
);
892 warning_error(wi
, "Not enough parameters");
894 else if (nn
> nrfpA
&& nn
< nrfp
)
896 warning_error(wi
, "Too many parameters or not enough parameters for topology B");
900 warning_error(wi
, "Too many parameters");
902 for (i
= nn
; (i
< nrfp
); i
++)
909 for (i
= 0; (i
< 4); i
++)
911 if (!strcmp(alc
[i
], "X"))
917 if ((p
.a
[i
] = get_bond_atomtype_type(alc
[i
], bat
)) == NOTSET
)
919 auto message
= gmx::formatString("Unknown bond_atomtype %s", alc
[i
]);
920 warning_error_and_exit(wi
, message
, FARGS
);
924 for (i
= 0; (i
< nrfp
); i
++)
928 /* Always use 4 atoms here, since we created two wildcard atoms
929 * if there wasn't of them 4 already.
931 push_bondtype (&(bt
[ftype
]), &p
, 4, ftype
, bAllowRepeat
, line
, wi
);
935 void push_nbt(directive d
, t_nbparam
**nbt
, gpp_atomtype_t atype
,
936 char *pline
, int nb_funct
,
940 const char *form3
= "%*s%*s%*s%lf%lf%lf";
941 const char *form4
= "%*s%*s%*s%lf%lf%lf%lf";
942 const char *form5
= "%*s%*s%*s%lf%lf%lf%lf%lf";
944 int i
, f
, n
, ftype
, nrfp
;
951 if (sscanf (pline
, "%s%s%d", a0
, a1
, &f
) != 3)
957 ftype
= ifunc_index(d
, f
);
959 if (ftype
!= nb_funct
)
961 auto message
= gmx::formatString("Trying to add %s while the default nonbond type is %s",
962 interaction_function
[ftype
].longname
,
963 interaction_function
[nb_funct
].longname
);
964 warning_error(wi
, message
);
968 /* Get the force parameters */
972 n
= sscanf(pline
, form4
, &c
[0], &c
[1], &c
[2], &c
[3]);
978 /* When the B topology parameters are not set,
979 * copy them from topology A
981 GMX_ASSERT(nrfp
<= 4, "LJ-14 cannot have more than 4 parameters");
982 for (i
= n
; i
< nrfp
; i
++)
987 else if (ftype
== F_LJC14_Q
)
989 n
= sscanf(pline
, form5
, &c
[0], &c
[1], &c
[2], &c
[3], &dum
);
992 incorrect_n_param(wi
);
998 if (sscanf(pline
, form3
, &c
[0], &c
[1], &dum
) != 2)
1000 incorrect_n_param(wi
);
1006 if (sscanf(pline
, form4
, &c
[0], &c
[1], &c
[2], &dum
) != 3)
1008 incorrect_n_param(wi
);
1014 auto message
= gmx::formatString("Number of force parameters for nonbonded interactions is %d", nrfp
);
1015 warning_error_and_exit(wi
, message
, FARGS
);
1017 for (i
= 0; (i
< nrfp
); i
++)
1022 /* Put the parameters in the matrix */
1023 if ((ai
= get_atomtype_type (a0
, atype
)) == NOTSET
)
1025 auto message
= gmx::formatString("Atomtype %s not found", a0
);
1026 warning_error_and_exit(wi
, message
, FARGS
);
1028 if ((aj
= get_atomtype_type (a1
, atype
)) == NOTSET
)
1030 auto message
= gmx::formatString("Atomtype %s not found", a1
);
1031 warning_error_and_exit(wi
, message
, FARGS
);
1033 nbp
= &(nbt
[std::max(ai
, aj
)][std::min(ai
, aj
)]);
1038 for (i
= 0; i
< nrfp
; i
++)
1040 bId
= bId
&& (nbp
->c
[i
] == cr
[i
]);
1044 auto message
= gmx::formatString("Overriding non-bonded parameters,");
1045 warning(wi
, message
);
1046 fprintf(stderr
, " old:");
1047 for (i
= 0; i
< nrfp
; i
++)
1049 fprintf(stderr
, " %g", nbp
->c
[i
]);
1051 fprintf(stderr
, " new\n%s\n", pline
);
1055 for (i
= 0; i
< nrfp
; i
++)
1062 push_cmaptype(directive d
, t_params bt
[], int nral
, gpp_atomtype_t at
,
1063 t_bond_atomtype bat
, char *line
,
1066 const char *formal
= "%s%s%s%s%s%s%s%s%n";
1068 int i
, ft
, ftype
, nn
, nrfp
, nrfpA
, nrfpB
;
1069 int start
, nchar_consumed
;
1070 int nxcmap
, nycmap
, ncmap
, read_cmap
, sl
, nct
;
1071 char s
[20], alc
[MAXATOMLIST
+2][20];
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 auto message
= gmx::formatString("Incorrect number of atomtypes for cmap (%d instead of 5)", nn
-1);
1084 warning_error(wi
, message
);
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 auto message
= gmx::formatString("Not the same grid spacing in x and y for cmap grid: x=%d, y=%d", nxcmap
, nycmap
);
1097 warning_error(wi
, message
);
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 auto message
= gmx::formatString("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
, message
);
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 auto message
= gmx::formatString("Unknown atomtype %s\n", alc
[i
]);
1174 warning_error(wi
, message
);
1176 else if (bat
&& ((p
.a
[i
] = get_bond_atomtype_type(alc
[i
], bat
)) == NOTSET
))
1178 auto message
= gmx::formatString("Unknown bond_atomtype %s\n", alc
[i
]);
1179 warning_error(wi
, message
);
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 auto message
= gmx::formatString("Incorrect number of atom types (%d) in cmap type %d\n", nct
, bt
[F_CMAP
].nc
);
1193 warning_error(wi
, message
);
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
;
1219 if (((nr
== 0) && (atomnr
!= 1)) || (nr
&& (atomnr
!= at
->nr
+1)))
1221 auto message
= gmx::formatString
1222 ("Atoms in the .top are not numbered consecutively from 1 (rather, "
1223 "atomnr = %d, while at->nr = %d)", atomnr
, at
->nr
);
1224 warning_error_and_exit(wi
, message
, FARGS
);
1227 j
= strlen(resnumberic
) - 1;
1228 if (isdigit(resnumberic
[j
]))
1234 ric
= resnumberic
[j
];
1235 if (j
== 0 || !isdigit(resnumberic
[j
-1]))
1237 auto message
= gmx::formatString("Invalid residue number '%s' for atom %d",
1238 resnumberic
, atomnr
);
1239 warning_error_and_exit(wi
, message
, FARGS
);
1242 resnr
= strtol(resnumberic
, nullptr, 10);
1246 resind
= at
->atom
[nr
-1].resind
;
1248 if (nr
== 0 || strcmp(resname
, *at
->resinfo
[resind
].name
) != 0 ||
1249 resnr
!= at
->resinfo
[resind
].nr
||
1250 ric
!= at
->resinfo
[resind
].ic
)
1260 at
->nres
= resind
+ 1;
1261 srenew(at
->resinfo
, at
->nres
);
1262 at
->resinfo
[resind
].name
= put_symtab(symtab
, resname
);
1263 at
->resinfo
[resind
].nr
= resnr
;
1264 at
->resinfo
[resind
].ic
= ric
;
1268 resind
= at
->atom
[at
->nr
-1].resind
;
1271 /* New atom instance
1272 * get new space for arrays
1274 srenew(at
->atom
, nr
+1);
1275 srenew(at
->atomname
, nr
+1);
1276 srenew(at
->atomtype
, nr
+1);
1277 srenew(at
->atomtypeB
, nr
+1);
1280 at
->atom
[nr
].type
= type
;
1281 at
->atom
[nr
].ptype
= ptype
;
1282 at
->atom
[nr
].q
= q0
;
1283 at
->atom
[nr
].m
= m0
;
1284 at
->atom
[nr
].typeB
= typeB
;
1285 at
->atom
[nr
].qB
= qB
;
1286 at
->atom
[nr
].mB
= mB
;
1288 at
->atom
[nr
].resind
= resind
;
1289 at
->atom
[nr
].atomnumber
= atomicnumber
;
1290 at
->atomname
[nr
] = put_symtab(symtab
, name
);
1291 at
->atomtype
[nr
] = put_symtab(symtab
, ctype
);
1292 at
->atomtypeB
[nr
] = put_symtab(symtab
, ctypeB
);
1296 static void push_cg(t_block
*block
, int *lastindex
, int index
, int a
)
1298 if (((block
->nr
) && (*lastindex
!= index
)) || (!block
->nr
))
1300 /* add a new block */
1302 srenew(block
->index
, block
->nr
+1);
1304 block
->index
[block
->nr
] = a
+ 1;
1308 void push_atom(t_symtab
*symtab
, t_block
*cgs
,
1309 t_atoms
*at
, gpp_atomtype_t atype
, char *line
, int *lastcg
,
1313 int cgnumber
, atomnr
, type
, typeB
, nscan
;
1314 char id
[STRLEN
], ctype
[STRLEN
], ctypeB
[STRLEN
],
1315 resnumberic
[STRLEN
], resname
[STRLEN
], name
[STRLEN
], check
[STRLEN
];
1316 double m
, q
, mb
, qb
;
1317 real m0
, q0
, mB
, qB
;
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 auto message
= gmx::formatString("Atomtype %s not found", ctype
);
1333 warning_error_and_exit(wi
, message
, 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 auto message
= gmx::formatString("Atomtype %s not found", ctypeB
);
1360 warning_error_and_exit(wi
, message
, 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
,
1395 if ((sscanf(line
, "%s%d", type
, &nrexcl
)) != 2)
1397 warning_error(wi
, "Expected a molecule type name and nrexcl");
1400 /* Test if this moleculetype overwrites another */
1404 if (strcmp(*((*mol
)[i
].name
), type
) == 0)
1406 auto message
= gmx::formatString("moleculetype %s is redefined", type
);
1407 warning_error_and_exit(wi
, message
, FARGS
);
1413 srenew(*mol
, *nmol
);
1414 newmol
= &((*mol
)[*nmol
-1]);
1415 init_molinfo(newmol
);
1417 /* Fill in the values */
1418 newmol
->name
= put_symtab(symtab
, type
);
1419 newmol
->nrexcl
= nrexcl
;
1420 newmol
->excl_set
= FALSE
;
1423 static bool default_nb_params(int ftype
, t_params bt
[], t_atoms
*at
,
1424 t_param
*p
, int c_start
, bool bB
, bool bGenPairs
)
1426 int i
, j
, ti
, tj
, ntype
;
1428 t_param
*pi
= nullptr;
1429 int nr
= bt
[ftype
].nr
;
1430 int nral
= NRAL(ftype
);
1431 int nrfp
= interaction_function
[ftype
].nrfpA
;
1432 int nrfpB
= interaction_function
[ftype
].nrfpB
;
1434 if ((!bB
&& nrfp
== 0) || (bB
&& nrfpB
== 0))
1442 /* First test the generated-pair position to save
1443 * time when we have 1000*1000 entries for e.g. OPLS...
1445 ntype
= static_cast<int>(std::sqrt(static_cast<double>(nr
)));
1446 GMX_ASSERT(ntype
* ntype
== nr
, "Number of pairs of generated non-bonded parameters should be a perfect square");
1449 ti
= at
->atom
[p
->a
[0]].typeB
;
1450 tj
= at
->atom
[p
->a
[1]].typeB
;
1454 ti
= at
->atom
[p
->a
[0]].type
;
1455 tj
= at
->atom
[p
->a
[1]].type
;
1457 pi
= &(bt
[ftype
].param
[ntype
*ti
+tj
]);
1458 bFound
= ((ti
== pi
->a
[0]) && (tj
== pi
->a
[1]));
1461 /* Search explicitly if we didnt find it */
1464 for (i
= 0; ((i
< nr
) && !bFound
); i
++)
1466 pi
= &(bt
[ftype
].param
[i
]);
1469 for (j
= 0; ((j
< nral
) &&
1470 (at
->atom
[p
->a
[j
]].typeB
== pi
->a
[j
])); j
++)
1477 for (j
= 0; ((j
< nral
) &&
1478 (at
->atom
[p
->a
[j
]].type
== pi
->a
[j
])); j
++)
1483 bFound
= (j
== nral
);
1491 if (nrfp
+nrfpB
> MAXFORCEPARAM
)
1493 gmx_incons("Too many force parameters");
1495 for (j
= c_start
; (j
< nrfpB
); j
++)
1497 p
->c
[nrfp
+j
] = pi
->c
[j
];
1502 for (j
= c_start
; (j
< nrfp
); j
++)
1510 for (j
= c_start
; (j
< nrfp
); j
++)
1518 static bool default_cmap_params(t_params bondtype
[],
1519 t_atoms
*at
, gpp_atomtype_t atype
,
1520 t_param
*p
, bool bB
,
1521 int *cmap_type
, int *nparam_def
,
1524 int i
, nparam_found
;
1526 bool bFound
= FALSE
;
1531 /* Match the current cmap angle against the list of cmap_types */
1532 for (i
= 0; i
< bondtype
[F_CMAP
].nct
&& !bFound
; i
+= 6)
1541 (get_atomtype_batype(at
->atom
[p
->a
[0]].type
, atype
) == bondtype
[F_CMAP
].cmap_types
[i
]) &&
1542 (get_atomtype_batype(at
->atom
[p
->a
[1]].type
, atype
) == bondtype
[F_CMAP
].cmap_types
[i
+1]) &&
1543 (get_atomtype_batype(at
->atom
[p
->a
[2]].type
, atype
) == bondtype
[F_CMAP
].cmap_types
[i
+2]) &&
1544 (get_atomtype_batype(at
->atom
[p
->a
[3]].type
, atype
) == bondtype
[F_CMAP
].cmap_types
[i
+3]) &&
1545 (get_atomtype_batype(at
->atom
[p
->a
[4]].type
, atype
) == bondtype
[F_CMAP
].cmap_types
[i
+4]))
1547 /* Found cmap torsion */
1549 ct
= bondtype
[F_CMAP
].cmap_types
[i
+5];
1555 /* If we did not find a matching type for this cmap torsion */
1558 auto message
= gmx::formatString("Unknown cmap torsion between atoms %d %d %d %d %d",
1559 p
->a
[0]+1, p
->a
[1]+1, p
->a
[2]+1, p
->a
[3]+1, p
->a
[4]+1);
1560 warning_error_and_exit(wi
, message
, FARGS
);
1563 *nparam_def
= nparam_found
;
1569 /* Returns the number of exact atom type matches, i.e. non wild-card matches,
1570 * returns -1 when there are no matches at all.
1572 static int natom_match(t_param
*pi
,
1573 int type_i
, int type_j
, int type_k
, int type_l
,
1574 const gpp_atomtype
* atype
)
1576 if ((pi
->ai() == -1 || get_atomtype_batype(type_i
, atype
) == pi
->ai()) &&
1577 (pi
->aj() == -1 || get_atomtype_batype(type_j
, atype
) == pi
->aj()) &&
1578 (pi
->ak() == -1 || get_atomtype_batype(type_k
, atype
) == pi
->ak()) &&
1579 (pi
->al() == -1 || get_atomtype_batype(type_l
, atype
) == pi
->al()))
1582 (pi
->ai() == -1 ? 0 : 1) +
1583 (pi
->aj() == -1 ? 0 : 1) +
1584 (pi
->ak() == -1 ? 0 : 1) +
1585 (pi
->al() == -1 ? 0 : 1);
1593 static bool default_params(int ftype
, t_params bt
[],
1594 t_atoms
*at
, gpp_atomtype_t atype
,
1595 t_param
*p
, bool bB
,
1596 t_param
**param_def
,
1601 t_param
*pi
= nullptr;
1602 t_param
*pj
= nullptr;
1603 int nr
= bt
[ftype
].nr
;
1604 int nral
= NRAL(ftype
);
1605 int nrfpA
= interaction_function
[ftype
].nrfpA
;
1606 int nrfpB
= interaction_function
[ftype
].nrfpB
;
1608 if ((!bB
&& nrfpA
== 0) || (bB
&& nrfpB
== 0))
1616 if (ftype
== F_PDIHS
|| ftype
== F_RBDIHS
|| ftype
== F_IDIHS
|| ftype
== F_PIDIHS
)
1618 int nmatch_max
= -1;
1622 /* For dihedrals we allow wildcards. We choose the first type
1623 * that has the most real matches, i.e. non-wildcard matches.
1625 for (t
= 0; ((t
< nr
) && nmatch_max
< 4); t
++)
1630 pt
= &(bt
[ftype
].param
[t
]);
1633 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
);
1637 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
);
1639 if (nmatch
> nmatch_max
)
1641 nmatch_max
= nmatch
;
1651 pi
= &(bt
[ftype
].param
[i
]);
1654 /* Find additional matches for this dihedral - necessary
1656 * The rule in that case is that additional matches
1657 * HAVE to be on adjacent lines!
1660 /* Continue from current i value */
1661 for (j
= i
+ 2; j
< nr
&& bSame
; j
+= 2)
1663 pj
= &(bt
[ftype
].param
[j
]);
1664 bSame
= (pi
->ai() == pj
->ai() && pi
->aj() == pj
->aj() && pi
->ak() == pj
->ak() && pi
->al() == pj
->al());
1669 /* nparam_found will be increased as long as the numbers match */
1673 else /* Not a dihedral */
1677 for (i
= 0; ((i
< nr
) && !bFound
); i
++)
1679 pi
= &(bt
[ftype
].param
[i
]);
1682 for (j
= 0; ((j
< nral
) &&
1683 (get_atomtype_batype(at
->atom
[p
->a
[j
]].typeB
, atype
) == pi
->a
[j
])); j
++)
1690 for (j
= 0; ((j
< nral
) &&
1691 (get_atomtype_batype(at
->atom
[p
->a
[j
]].type
, atype
) == pi
->a
[j
])); j
++)
1696 bFound
= (j
== nral
);
1705 *nparam_def
= nparam_found
;
1712 void push_bond(directive d
, t_params bondtype
[], t_params bond
[],
1713 t_atoms
*at
, gpp_atomtype_t atype
, char *line
,
1714 bool bBonded
, bool bGenPairs
, real fudgeQQ
,
1715 bool bZero
, bool *bWarn_copy_A_B
,
1718 const char *aaformat
[MAXATOMLIST
] = {
1726 const char *asformat
[MAXATOMLIST
] = {
1731 "%*s%*s%*s%*s%*s%*s",
1732 "%*s%*s%*s%*s%*s%*s%*s"
1734 const char *ccformat
= "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
1735 int nr
, i
, j
, nral
, nral_fmt
, nread
, ftype
;
1736 char format
[STRLEN
];
1737 /* One force parameter more, so we can check if we read too many */
1738 double cc
[MAXFORCEPARAM
+1];
1739 int aa
[MAXATOMLIST
+1];
1740 t_param param
, *param_defA
, *param_defB
;
1741 bool bFoundA
= FALSE
, bFoundB
= FALSE
, bDef
, bPert
, bSwapParity
= FALSE
;
1742 int nparam_defA
, nparam_defB
;
1744 nparam_defA
= nparam_defB
= 0;
1746 ftype
= ifunc_index(d
, 1);
1748 for (j
= 0; j
< MAXATOMLIST
; j
++)
1752 bDef
= (NRFP(ftype
) > 0);
1754 if (ftype
== F_SETTLE
)
1756 /* SETTLE acts on 3 atoms, but the topology format only specifies
1757 * the first atom (for historical reasons).
1766 nread
= sscanf(line
, aaformat
[nral_fmt
-1],
1767 &aa
[0], &aa
[1], &aa
[2], &aa
[3], &aa
[4], &aa
[5]);
1769 if (ftype
== F_SETTLE
)
1776 if (nread
< nral_fmt
)
1781 else if (nread
> nral_fmt
)
1783 /* this is a hack to allow for virtual sites with swapped parity */
1784 bSwapParity
= (aa
[nral
] < 0);
1787 aa
[nral
] = -aa
[nral
];
1789 ftype
= ifunc_index(d
, aa
[nral
]);
1798 auto message
= gmx::formatString("Negative function types only allowed for %s and %s",
1799 interaction_function
[F_VSITE3FAD
].longname
,
1800 interaction_function
[F_VSITE3OUT
].longname
);
1801 warning_error_and_exit(wi
, message
, FARGS
);
1807 /* Check for double atoms and atoms out of bounds */
1808 for (i
= 0; (i
< nral
); i
++)
1810 if (aa
[i
] < 1 || aa
[i
] > at
->nr
)
1812 auto message
= gmx::formatString
1813 ("Atom index (%d) in %s out of bounds (1-%d).\n"
1814 "This probably means that you have inserted topology section \"%s\"\n"
1815 "in a part belonging to a different molecule than you intended to.\n"
1816 "In that case move the \"%s\" section to the right molecule.",
1817 aa
[i
], dir2str(d
), at
->nr
, dir2str(d
), dir2str(d
));
1818 warning_error_and_exit(wi
, message
, FARGS
);
1820 for (j
= i
+1; (j
< nral
); j
++)
1822 GMX_ASSERT(j
< MAXATOMLIST
+ 1, "Values from nral=NRAL() will satisfy this, we assert to keep gcc 4 happy");
1825 auto message
= gmx::formatString("Duplicate atom index (%d) in %s", aa
[i
], dir2str(d
));
1826 if (ftype
== F_ANGRES
)
1828 /* Since the angle restraints uses 2 pairs of atoms to
1829 * defines an angle between vectors, it can be useful
1830 * to use one atom twice, so we only issue a note here.
1832 warning_note(wi
, message
);
1836 warning_error(wi
, message
);
1842 /* default force parameters */
1843 for (j
= 0; (j
< MAXATOMLIST
); j
++)
1845 param
.a
[j
] = aa
[j
]-1;
1847 for (j
= 0; (j
< MAXFORCEPARAM
); j
++)
1852 /* Get force params for normal and free energy perturbation
1853 * studies, as determined by types!
1858 bFoundA
= default_params(ftype
, bondtype
, at
, atype
, ¶m
, FALSE
, ¶m_defA
, &nparam_defA
);
1861 /* Copy the A-state and B-state default parameters. */
1862 GMX_ASSERT(NRFPA(ftype
)+NRFPB(ftype
) <= MAXFORCEPARAM
, "Bonded interactions may have at most 12 parameters");
1863 for (j
= 0; (j
< NRFPA(ftype
)+NRFPB(ftype
)); j
++)
1865 param
.c
[j
] = param_defA
->c
[j
];
1868 bFoundB
= default_params(ftype
, bondtype
, at
, atype
, ¶m
, TRUE
, ¶m_defB
, &nparam_defB
);
1871 /* Copy only the B-state default parameters */
1872 for (j
= NRFPA(ftype
); (j
< NRFP(ftype
)); j
++)
1874 param
.c
[j
] = param_defB
->c
[j
];
1878 else if (ftype
== F_LJ14
)
1880 bFoundA
= default_nb_params(ftype
, bondtype
, at
, ¶m
, 0, FALSE
, bGenPairs
);
1881 bFoundB
= default_nb_params(ftype
, bondtype
, at
, ¶m
, 0, TRUE
, bGenPairs
);
1883 else if (ftype
== F_LJC14_Q
)
1885 param
.c
[0] = fudgeQQ
;
1886 /* Fill in the A-state charges as default parameters */
1887 param
.c
[1] = at
->atom
[param
.a
[0]].q
;
1888 param
.c
[2] = at
->atom
[param
.a
[1]].q
;
1889 /* The default LJ parameters are the standard 1-4 parameters */
1890 bFoundA
= default_nb_params(F_LJ14
, bondtype
, at
, ¶m
, 3, FALSE
, bGenPairs
);
1893 else if (ftype
== F_LJC_PAIRS_NB
)
1895 /* Defaults are not supported here */
1901 gmx_incons("Unknown function type in push_bond");
1904 if (nread
> nral_fmt
)
1906 /* Manually specified parameters - in this case we discard multiple torsion info! */
1908 strcpy(format
, asformat
[nral_fmt
-1]);
1909 strcat(format
, ccformat
);
1911 nread
= sscanf(line
, format
, &cc
[0], &cc
[1], &cc
[2], &cc
[3], &cc
[4], &cc
[5],
1912 &cc
[6], &cc
[7], &cc
[8], &cc
[9], &cc
[10], &cc
[11], &cc
[12]);
1914 if ((nread
== NRFPA(ftype
)) && (NRFPB(ftype
) != 0))
1916 /* We only have to issue a warning if these atoms are perturbed! */
1918 for (j
= 0; (j
< nral
); j
++)
1920 bPert
= bPert
|| PERTURBED(at
->atom
[param
.a
[j
]]);
1923 if (bPert
&& *bWarn_copy_A_B
)
1925 auto message
= gmx::formatString("Some parameters for bonded interaction involving "
1926 "perturbed atoms are specified explicitly in "
1927 "state A, but not B - copying A to B");
1928 warning(wi
, message
);
1929 *bWarn_copy_A_B
= FALSE
;
1932 /* If only the A parameters were specified, copy them to the B state */
1933 /* The B-state parameters correspond to the first nrfpB
1934 * A-state parameters.
1936 for (j
= 0; (j
< NRFPB(ftype
)); j
++)
1938 cc
[nread
++] = cc
[j
];
1942 /* If nread was 0 or EOF, no parameters were read => use defaults.
1943 * If nread was nrfpA we copied above so nread=nrfp.
1944 * If nread was nrfp we are cool.
1945 * For F_LJC14_Q we allow supplying fudgeQQ only.
1946 * Anything else is an error!
1948 if ((nread
!= 0) && (nread
!= EOF
) && (nread
!= NRFP(ftype
)) &&
1949 !(ftype
== F_LJC14_Q
&& nread
== 1))
1951 auto message
= gmx::formatString
1952 ("Incorrect number of parameters - found %d, expected %d "
1953 "or %d for %s (after the function type).",
1954 nread
, NRFPA(ftype
), NRFP(ftype
), interaction_function
[ftype
].longname
);
1955 warning_error_and_exit(wi
, message
, FARGS
);
1958 for (j
= 0; (j
< nread
); j
++)
1963 /* Check whether we have to use the defaults */
1964 if (nread
== NRFP(ftype
))
1973 /* nread now holds the number of force parameters read! */
1978 /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */
1979 if (ftype
== F_PDIHS
)
1981 if ((nparam_defA
!= nparam_defB
) || ((nparam_defA
> 1 || nparam_defB
> 1) && (param_defA
!= param_defB
)))
1984 gmx::formatString("Cannot automatically perturb a torsion with multiple terms to different form.\n"
1985 "Please specify perturbed parameters manually for this torsion in your topology!");
1986 warning_error(wi
, message
);
1990 if (nread
> 0 && nread
< NRFPA(ftype
))
1992 /* Issue an error, do not use defaults */
1993 auto message
= gmx::formatString("Not enough parameters, there should be at least %d (or 0 for defaults)", NRFPA(ftype
));
1994 warning_error(wi
, message
);
1997 if (nread
== 0 || nread
== EOF
)
2001 if (interaction_function
[ftype
].flags
& IF_VSITE
)
2003 /* set them to NOTSET, will be calculated later */
2004 for (j
= 0; (j
< MAXFORCEPARAM
); j
++)
2006 param
.c
[j
] = NOTSET
;
2011 param
.c1() = -1; /* flag to swap parity of vsite construction */
2018 fprintf(stderr
, "NOTE: No default %s types, using zeroes\n",
2019 interaction_function
[ftype
].longname
);
2023 auto message
= gmx::formatString("No default %s types", interaction_function
[ftype
].longname
);
2024 warning_error(wi
, message
);
2035 param
.c0() = 360 - param
.c0();
2038 param
.c2() = -param
.c2();
2045 /* We only have to issue a warning if these atoms are perturbed! */
2047 for (j
= 0; (j
< nral
); j
++)
2049 bPert
= bPert
|| PERTURBED(at
->atom
[param
.a
[j
]]);
2054 auto message
= gmx::formatString("No default %s types for perturbed atoms, "
2055 "using normal values", interaction_function
[ftype
].longname
);
2056 warning(wi
, message
);
2062 if ((ftype
== F_PDIHS
|| ftype
== F_ANGRES
|| ftype
== F_ANGRESZ
)
2063 && param
.c
[5] != param
.c
[2])
2065 auto message
= gmx::formatString("%s multiplicity can not be perturbed %f!=%f",
2066 interaction_function
[ftype
].longname
,
2067 param
.c
[2], param
.c
[5]);
2068 warning_error_and_exit(wi
, message
, FARGS
);
2071 if (IS_TABULATED(ftype
) && param
.c
[2] != param
.c
[0])
2073 auto message
= gmx::formatString("%s table number can not be perturbed %d!=%d",
2074 interaction_function
[ftype
].longname
,
2075 gmx::roundToInt(param
.c
[0]), gmx::roundToInt(param
.c
[2]));
2076 warning_error_and_exit(wi
, message
, FARGS
);
2079 /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
2080 if (ftype
== F_RBDIHS
)
2083 for (i
= 0; i
< NRFP(ftype
); i
++)
2085 if (param
.c
[i
] != 0)
2096 /* Put the values in the appropriate arrays */
2097 add_param_to_list (&bond
[ftype
], ¶m
);
2099 /* Push additional torsions from FF for ftype==9 if we have them.
2100 * We have already checked that the A/B states do not differ in this case,
2101 * so we do not have to double-check that again, or the vsite stuff.
2102 * In addition, those torsions cannot be automatically perturbed.
2104 if (bDef
&& ftype
== F_PDIHS
)
2106 for (i
= 1; i
< nparam_defA
; i
++)
2108 /* Advance pointer! */
2110 for (j
= 0; (j
< NRFPA(ftype
)+NRFPB(ftype
)); j
++)
2112 param
.c
[j
] = param_defA
->c
[j
];
2114 /* And push the next term for this torsion */
2115 add_param_to_list (&bond
[ftype
], ¶m
);
2120 void push_cmap(directive d
, t_params bondtype
[], t_params bond
[],
2121 t_atoms
*at
, gpp_atomtype_t atype
, char *line
,
2124 const char *aaformat
[MAXATOMLIST
+1] =
2135 int i
, j
, ftype
, nral
, nread
, ncmap_params
;
2137 int aa
[MAXATOMLIST
+1];
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 auto message
= gmx::formatString
2164 ("Atom index (%d) in %s out of bounds (1-%d).\n"
2165 "This probably means that you have inserted topology section \"%s\"\n"
2166 "in a part belonging to a different molecule than you intended to.\n"
2167 "In that case move the \"%s\" section to the right molecule.",
2168 aa
[i
], dir2str(d
), at
->nr
, dir2str(d
), dir2str(d
));
2169 warning_error_and_exit(wi
, message
, FARGS
);
2172 for (j
= i
+1; (j
< nral
); j
++)
2176 auto message
= gmx::formatString("Duplicate atom index (%d) in %s", aa
[i
], dir2str(d
));
2177 warning_error(wi
, message
);
2182 /* default force parameters */
2183 for (j
= 0; (j
< MAXATOMLIST
); j
++)
2185 param
.a
[j
] = aa
[j
]-1;
2187 for (j
= 0; (j
< MAXFORCEPARAM
); j
++)
2192 /* Get the cmap type for this cmap angle */
2193 bFound
= default_cmap_params(bondtype
, at
, atype
, ¶m
, FALSE
, &cmap_type
, &ncmap_params
, wi
);
2195 /* We want exactly one parameter (the cmap type in state A (currently no state B) back */
2196 if (bFound
&& ncmap_params
== 1)
2198 /* Put the values in the appropriate arrays */
2199 param
.c
[0] = cmap_type
;
2200 add_param_to_list(&bond
[ftype
], ¶m
);
2204 /* This is essentially the same check as in default_cmap_params() done one more time */
2205 auto message
= gmx::formatString("Unable to assign a cmap type to torsion %d %d %d %d and %d\n",
2206 param
.a
[0]+1, param
.a
[1]+1, param
.a
[2]+1, param
.a
[3]+1, param
.a
[4]+1);
2207 warning_error_and_exit(wi
, message
, FARGS
);
2213 void push_vsitesn(directive d
, t_params bond
[],
2214 t_atoms
*at
, char *line
,
2218 int type
, ftype
, j
, n
, ret
, nj
, a
;
2220 double *weight
= nullptr, weight_tot
;
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 auto message
= gmx::formatString("Expected an atom index in section \"%s\"",
2240 warning_error_and_exit(wi
, message
, 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 auto message
= gmx::formatString
2281 ("No weight or negative weight found for vsiten "
2282 "constructing atom %d (atom index %d)",
2284 warning_error_and_exit(wi
, message
, FARGS
);
2288 auto message
= gmx::formatString("Unknown vsiten type %d", type
);
2289 warning_error_and_exit(wi
, message
, FARGS
);
2291 weight_tot
+= weight
[nj
];
2299 auto message
= gmx::formatString("Expected more than one atom index in section \"%s\"",
2301 warning_error_and_exit(wi
, message
, FARGS
);
2304 if (weight_tot
== 0)
2306 warning_error_and_exit(wi
, "The total mass of the construting atoms is zero", FARGS
);
2309 for (j
= 0; j
< nj
; j
++)
2311 param
.a
[1] = atc
[j
];
2313 param
.c
[1] = weight
[j
]/weight_tot
;
2314 /* Put the values in the appropriate arrays */
2315 add_param_to_list (&bond
[ftype
], ¶m
);
2322 void push_mol(int nrmols
, t_molinfo mols
[], char *pline
, int *whichmol
,
2328 if (sscanf(pline
, "%s%d", type
, nrcopies
) != 2)
2334 /* Search moleculename.
2335 * Here we originally only did case insensitive matching. But because
2336 * some PDB files can have many chains and use case to generate more
2337 * chain-identifiers, which in turn end up in our moleculetype name,
2338 * we added support for case-sensitivity.
2344 for (int i
= 0; i
< nrmols
; i
++)
2346 if (strcmp(type
, *(mols
[i
].name
)) == 0)
2351 if (gmx_strcasecmp(type
, *(mols
[i
].name
)) == 0)
2360 // select the case sensitive match
2361 *whichmol
= matchcs
;
2365 // avoid matching case-insensitive when we have multiple matches
2368 auto message
= gmx::formatString
2369 ("For moleculetype '%s' in [ system ] %d case insensitive "
2370 "matches, but %d case sensitive matches were found. Check "
2371 "the case of the characters in the moleculetypes.",
2373 warning_error_and_exit(wi
, message
, FARGS
);
2377 // select the unique case insensitive match
2378 *whichmol
= matchci
;
2382 auto message
= gmx::formatString("No such moleculetype %s", type
);
2383 warning_error_and_exit(wi
, message
, FARGS
);
2388 void push_excl(char *line
, gmx::ExclusionBlocks
*b2
, warninp_t wi
)
2392 char base
[STRLEN
], format
[STRLEN
];
2394 if (sscanf(line
, "%d", &i
) == 0)
2399 if ((1 <= i
) && (i
<= b2
->nr
))
2407 strcpy(base
, "%*d");
2410 strcpy(format
, base
);
2411 strcat(format
, "%d");
2412 n
= sscanf(line
, format
, &j
);
2415 if ((1 <= j
) && (j
<= b2
->nr
))
2418 srenew(b2
->a
[i
], ++(b2
->nra
[i
]));
2419 b2
->a
[i
][b2
->nra
[i
]-1] = j
;
2420 /* also add the reverse exclusion! */
2421 srenew(b2
->a
[j
], ++(b2
->nra
[j
]));
2422 b2
->a
[j
][b2
->nra
[j
]-1] = i
;
2423 strcat(base
, "%*d");
2427 auto message
= gmx::formatString("Invalid Atomnr j: %d, b2->nr: %d\n", j
, b2
->nr
);
2428 warning_error_and_exit(wi
, message
, FARGS
);
2435 int add_atomtype_decoupled(t_symtab
*symtab
, gpp_atomtype_t at
,
2436 t_nbparam
***nbparam
, t_nbparam
***pair
)
2442 /* Define an atom type with all parameters set to zero (no interactions) */
2445 /* Type for decoupled atoms could be anything,
2446 * this should be changed automatically later when required.
2448 atom
.ptype
= eptAtom
;
2449 for (i
= 0; (i
< MAXFORCEPARAM
); i
++)
2454 nr
= add_atomtype(at
, symtab
, &atom
, "decoupled", ¶m
, -1, 0);
2456 /* Add space in the non-bonded parameters matrix */
2457 realloc_nb_params(at
, nbparam
, pair
);
2462 static void convert_pairs_to_pairsQ(t_params
*plist
,
2463 real fudgeQQ
, t_atoms
*atoms
)
2465 t_param
*paramp1
, *paramp2
, *paramnew
;
2466 int i
, j
, p1nr
, p2nr
, p2newnr
;
2468 /* Add the pair list to the pairQ list */
2469 p1nr
= plist
[F_LJ14
].nr
;
2470 p2nr
= plist
[F_LJC14_Q
].nr
;
2471 p2newnr
= p1nr
+ p2nr
;
2472 snew(paramnew
, p2newnr
);
2474 paramp1
= plist
[F_LJ14
].param
;
2475 paramp2
= plist
[F_LJC14_Q
].param
;
2477 /* Fill in the new F_LJC14_Q array with the old one. NOTE:
2478 it may be possible to just ADD the converted F_LJ14 array
2479 to the old F_LJC14_Q array, but since we have to create
2480 a new sized memory structure, better just to deep copy it all.
2483 for (i
= 0; i
< p2nr
; i
++)
2485 /* Copy over parameters */
2486 for (j
= 0; j
< 5; j
++) /* entries are 0-4 for F_LJC14_Q */
2488 paramnew
[i
].c
[j
] = paramp2
[i
].c
[j
];
2491 /* copy over atoms */
2492 for (j
= 0; j
< 2; j
++)
2494 paramnew
[i
].a
[j
] = paramp2
[i
].a
[j
];
2498 for (i
= p2nr
; i
< p2newnr
; i
++)
2502 /* Copy over parameters */
2503 paramnew
[i
].c
[0] = fudgeQQ
;
2504 paramnew
[i
].c
[1] = atoms
->atom
[paramp1
[j
].a
[0]].q
;
2505 paramnew
[i
].c
[2] = atoms
->atom
[paramp1
[j
].a
[1]].q
;
2506 paramnew
[i
].c
[3] = paramp1
[j
].c
[0];
2507 paramnew
[i
].c
[4] = paramp1
[j
].c
[1];
2509 /* copy over atoms */
2510 paramnew
[i
].a
[0] = paramp1
[j
].a
[0];
2511 paramnew
[i
].a
[1] = paramp1
[j
].a
[1];
2514 /* free the old pairlists */
2515 sfree(plist
[F_LJC14_Q
].param
);
2516 sfree(plist
[F_LJ14
].param
);
2518 /* now assign the new data to the F_LJC14_Q structure */
2519 plist
[F_LJC14_Q
].param
= paramnew
;
2520 plist
[F_LJC14_Q
].nr
= p2newnr
;
2522 /* Empty the LJ14 pairlist */
2523 plist
[F_LJ14
].nr
= 0;
2524 plist
[F_LJ14
].param
= nullptr;
2527 static void generate_LJCpairsNB(t_molinfo
*mol
, int nb_funct
, t_params
*nbp
, warninp_t wi
)
2529 int n
, ntype
, i
, j
, k
;
2536 atom
= mol
->atoms
.atom
;
2538 ntype
= static_cast<int>(std::sqrt(static_cast<double>(nbp
->nr
)));
2539 GMX_ASSERT(ntype
* ntype
== nbp
->nr
, "Number of pairs of generated non-bonded parameters should be a perfect square");
2541 for (i
= 0; i
< MAXATOMLIST
; i
++)
2543 param
.a
[i
] = NOTSET
;
2545 for (i
= 0; i
< MAXFORCEPARAM
; i
++)
2547 param
.c
[i
] = NOTSET
;
2550 /* Add a pair interaction for all non-excluded atom pairs */
2552 for (i
= 0; i
< n
; i
++)
2554 for (j
= i
+1; j
< n
; j
++)
2557 for (k
= excl
->index
[i
]; k
< excl
->index
[i
+1]; k
++)
2559 if (excl
->a
[k
] == j
)
2566 if (nb_funct
!= F_LJ
)
2568 auto message
= gmx::formatString
2569 ("Can only generate non-bonded pair interactions "
2570 "for Van der Waals type Lennard-Jones");
2571 warning_error_and_exit(wi
, message
, FARGS
);
2575 param
.c
[0] = atom
[i
].q
;
2576 param
.c
[1] = atom
[j
].q
;
2577 param
.c
[2] = nbp
->param
[ntype
*atom
[i
].type
+atom
[j
].type
].c
[0];
2578 param
.c
[3] = nbp
->param
[ntype
*atom
[i
].type
+atom
[j
].type
].c
[1];
2579 add_param_to_list(&mol
->plist
[F_LJC_PAIRS_NB
], ¶m
);
2585 static void set_excl_all(t_blocka
*excl
)
2589 /* Get rid of the current exclusions and exclude all atom pairs */
2591 excl
->nra
= nat
*nat
;
2592 srenew(excl
->a
, excl
->nra
);
2594 for (i
= 0; i
< nat
; i
++)
2597 for (j
= 0; j
< nat
; j
++)
2602 excl
->index
[nat
] = k
;
2605 static void decouple_atoms(t_atoms
*atoms
, int atomtype_decouple
,
2606 int couple_lam0
, int couple_lam1
,
2607 const char *mol_name
, warninp_t wi
)
2611 for (i
= 0; i
< atoms
->nr
; i
++)
2615 atom
= &atoms
->atom
[i
];
2617 if (atom
->qB
!= atom
->q
|| atom
->typeB
!= atom
->type
)
2619 auto message
= gmx::formatString
2620 ("Atom %d in molecule type '%s' has different A and B state "
2621 "charges and/or atom types set in the topology file as well "
2622 "as through the mdp option '%s'. You can not use both "
2623 "these methods simultaneously.",
2624 i
+ 1, mol_name
, "couple-moltype");
2625 warning_error_and_exit(wi
, message
, FARGS
);
2628 if (couple_lam0
== ecouplamNONE
|| couple_lam0
== ecouplamVDW
)
2632 if (couple_lam0
== ecouplamNONE
|| couple_lam0
== ecouplamQ
)
2634 atom
->type
= atomtype_decouple
;
2636 if (couple_lam1
== ecouplamNONE
|| couple_lam1
== ecouplamVDW
)
2640 if (couple_lam1
== ecouplamNONE
|| couple_lam1
== ecouplamQ
)
2642 atom
->typeB
= atomtype_decouple
;
2647 void convert_moltype_couple(t_molinfo
*mol
, int atomtype_decouple
, real fudgeQQ
,
2648 int couple_lam0
, int couple_lam1
,
2649 bool bCoupleIntra
, int nb_funct
, t_params
*nbp
,
2652 convert_pairs_to_pairsQ(mol
->plist
, fudgeQQ
, &mol
->atoms
);
2655 generate_LJCpairsNB(mol
, nb_funct
, nbp
, wi
);
2656 set_excl_all(&mol
->excls
);
2658 decouple_atoms(&mol
->atoms
, atomtype_decouple
, couple_lam0
, couple_lam1
,