From: Roland Schulz Date: Mon, 29 Sep 2014 14:47:53 +0000 (-0400) Subject: Merge release-4-6 into release-5-0 X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?p=alexxy%2Fgromacs.git;a=commitdiff_plain;h=3cf05dfc83395c34cdbcfed4ceed2ddb16330f08 Merge release-4-6 into release-5-0 Change-Id: Ie72eccf57febab7b7ac8092ce55988b9cb2737af --- 3cf05dfc83395c34cdbcfed4ceed2ddb16330f08 diff --cc src/gromacs/gmxpreprocess/toppush.c index fd2a8a1595,0000000000..fab21dda51 mode 100644,000000..100644 --- a/src/gromacs/gmxpreprocess/toppush.c +++ b/src/gromacs/gmxpreprocess/toppush.c @@@ -1,2707 -1,0 +1,2709 @@@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 1991-2000, University of Groningen, The Netherlands. + * Copyright (c) 2001-2004, The GROMACS development team. + * Copyright (c) 2013,2014, by the GROMACS development team, led by + * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, + * and including many others, as listed in the AUTHORS file in the + * top-level source directory and at http://www.gromacs.org. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * http://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at http://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out http://www.gromacs.org. + */ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "sysstuff.h" +#include "gromacs/utility/smalloc.h" +#include "macros.h" +#include "gromacs/utility/cstringutil.h" +#include "names.h" +#include "toputil.h" +#include "toppush.h" +#include "topdirs.h" +#include "readir.h" +#include "symtab.h" +#include "gmx_fatal.h" +#include "warninp.h" +#include "gpp_atomtype.h" +#include "gpp_bond_atomtype.h" + +void generate_nbparams(int comb, int ftype, t_params *plist, gpp_atomtype_t atype, + warninp_t wi) +{ + int i, j, k = -1, nf; + int nr, nrfp; + real c, bi, bj, ci, cj, ci0, ci1, ci2, cj0, cj1, cj2; + char errbuf[256]; + + /* Lean mean shortcuts */ + nr = get_atomtype_ntypes(atype); + nrfp = NRFP(ftype); + snew(plist->param, nr*nr); + plist->nr = nr*nr; + + /* Fill the matrix with force parameters */ + switch (ftype) + { + case F_LJ: + switch (comb) + { + case eCOMB_GEOMETRIC: + /* Gromos rules */ + for (i = k = 0; (i < nr); i++) + { + for (j = 0; (j < nr); j++, k++) + { + for (nf = 0; (nf < nrfp); nf++) + { + ci = get_atomtype_nbparam(i, nf, atype); + cj = get_atomtype_nbparam(j, nf, atype); + c = sqrt(ci * cj); + plist->param[k].c[nf] = c; + } + } + } + break; + + case eCOMB_ARITHMETIC: + /* c0 and c1 are sigma and epsilon */ + for (i = k = 0; (i < nr); i++) + { + for (j = 0; (j < nr); j++, k++) + { + ci0 = get_atomtype_nbparam(i, 0, atype); + cj0 = get_atomtype_nbparam(j, 0, atype); + ci1 = get_atomtype_nbparam(i, 1, atype); + cj1 = get_atomtype_nbparam(j, 1, atype); + plist->param[k].c[0] = (fabs(ci0) + fabs(cj0))*0.5; + /* Negative sigma signals that c6 should be set to zero later, + * so we need to propagate that through the combination rules. + */ + if (ci0 < 0 || cj0 < 0) + { + plist->param[k].c[0] *= -1; + } + plist->param[k].c[1] = sqrt(ci1*cj1); + } + } + + break; + case eCOMB_GEOM_SIG_EPS: + /* c0 and c1 are sigma and epsilon */ + for (i = k = 0; (i < nr); i++) + { + for (j = 0; (j < nr); j++, k++) + { + ci0 = get_atomtype_nbparam(i, 0, atype); + cj0 = get_atomtype_nbparam(j, 0, atype); + ci1 = get_atomtype_nbparam(i, 1, atype); + cj1 = get_atomtype_nbparam(j, 1, atype); + plist->param[k].c[0] = sqrt(fabs(ci0*cj0)); + /* Negative sigma signals that c6 should be set to zero later, + * so we need to propagate that through the combination rules. + */ + if (ci0 < 0 || cj0 < 0) + { + plist->param[k].c[0] *= -1; + } + plist->param[k].c[1] = sqrt(ci1*cj1); + } + } + + break; + default: + gmx_fatal(FARGS, "No such combination rule %d", comb); + } + if (plist->nr != k) + { + gmx_incons("Topology processing, generate nb parameters"); + } + break; + + case F_BHAM: + /* Buckingham rules */ + for (i = k = 0; (i < nr); i++) + { + for (j = 0; (j < nr); j++, k++) + { + ci0 = get_atomtype_nbparam(i, 0, atype); + cj0 = get_atomtype_nbparam(j, 0, atype); + ci2 = get_atomtype_nbparam(i, 2, atype); + cj2 = get_atomtype_nbparam(j, 2, atype); + bi = get_atomtype_nbparam(i, 1, atype); + bj = get_atomtype_nbparam(j, 1, atype); + plist->param[k].c[0] = sqrt(ci0 * cj0); + if ((bi == 0) || (bj == 0)) + { + plist->param[k].c[1] = 0; + } + else + { + plist->param[k].c[1] = 2.0/(1/bi+1/bj); + } + plist->param[k].c[2] = sqrt(ci2 * cj2); + } + } + + break; + default: + sprintf(errbuf, "Invalid nonbonded type %s", + interaction_function[ftype].longname); + warning_error(wi, errbuf); + } +} + +static void realloc_nb_params(gpp_atomtype_t at, + t_nbparam ***nbparam, t_nbparam ***pair) +{ + /* Add space in the non-bonded parameters matrix */ + int atnr = get_atomtype_ntypes(at); + srenew(*nbparam, atnr); + snew((*nbparam)[atnr-1], atnr); + if (pair) + { + srenew(*pair, atnr); + snew((*pair)[atnr-1], atnr); + } +} + +static void copy_B_from_A(int ftype, double *c) +{ + int nrfpA, nrfpB, i; + + nrfpA = NRFPA(ftype); + nrfpB = NRFPB(ftype); + + /* Copy the B parameters from the first nrfpB A parameters */ + for (i = 0; (i < nrfpB); i++) + { + c[nrfpA+i] = c[i]; + } +} + +void push_at (t_symtab *symtab, gpp_atomtype_t at, t_bond_atomtype bat, + char *line, int nb_funct, + t_nbparam ***nbparam, t_nbparam ***pair, + warninp_t wi) +{ + typedef struct { + const char *entry; + int ptype; + } t_xlate; + t_xlate xl[eptNR] = { + { "A", eptAtom }, + { "N", eptNucleus }, + { "S", eptShell }, + { "B", eptBond }, + { "V", eptVSite }, + }; + + int nr, i, nfields, j, pt, nfp0 = -1; + int batype_nr, nread; + char type[STRLEN], btype[STRLEN], ptype[STRLEN]; + double m, q; + double c[MAXFORCEPARAM]; + double radius, vol, surftens, gb_radius, S_hct; + char tmpfield[12][100]; /* Max 12 fields of width 100 */ + char errbuf[256]; + t_atom *atom; + t_param *param; + int atomnr; + gmx_bool have_atomic_number; + gmx_bool have_bonded_type; + + snew(atom, 1); + snew(param, 1); + + /* First assign input line to temporary array */ + nfields = sscanf(line, "%s%s%s%s%s%s%s%s%s%s%s%s", + tmpfield[0], tmpfield[1], tmpfield[2], tmpfield[3], tmpfield[4], tmpfield[5], + tmpfield[6], tmpfield[7], tmpfield[8], tmpfield[9], tmpfield[10], tmpfield[11]); + + /* Comments on optional fields in the atomtypes section: + * + * The force field format is getting a bit old. For OPLS-AA we needed + * to add a special bonded atomtype, and for Gerrit Groenhofs QM/MM stuff + * we also needed the atomic numbers. + * To avoid making all old or user-generated force fields unusable we + * have introduced both these quantities as optional columns, and do some + * acrobatics to check whether they are present or not. + * This will all look much nicer when we switch to XML... sigh. + * + * Field 0 (mandatory) is the nonbonded type name. (string) + * Field 1 (optional) is the bonded type (string) + * Field 2 (optional) is the atomic number (int) + * Field 3 (mandatory) is the mass (numerical) + * Field 4 (mandatory) is the charge (numerical) + * Field 5 (mandatory) is the particle type (single character) + * This is followed by a number of nonbonded parameters. + * + * The safest way to identify the format is the particle type field. + * + * So, here is what we do: + * + * A. Read in the first six fields as strings + * B. If field 3 (starting from 0) is a single char, we have neither + * bonded_type or atomic numbers. + * C. If field 5 is a single char we have both. + * D. If field 4 is a single char we check field 1. If this begins with + * an alphabetical character we have bonded types, otherwise atomic numbers. + */ + + if (nfields < 6) + { + too_few(wi); + return; + } + + if ( (strlen(tmpfield[5]) == 1) && isalpha(tmpfield[5][0]) ) + { + have_bonded_type = TRUE; + have_atomic_number = TRUE; + } + else if ( (strlen(tmpfield[3]) == 1) && isalpha(tmpfield[3][0]) ) + { + have_bonded_type = FALSE; + have_atomic_number = FALSE; + } + else + { + have_bonded_type = ( isalpha(tmpfield[1][0]) != 0 ); + have_atomic_number = !have_bonded_type; + } + + /* optional fields */ + surftens = -1; + vol = -1; + radius = -1; + gb_radius = -1; + atomnr = -1; + S_hct = -1; + + switch (nb_funct) + { + + case F_LJ: + nfp0 = 2; + + if (have_atomic_number) + { + if (have_bonded_type) + { + nread = sscanf(line, "%s%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf", + type, btype, &atomnr, &m, &q, ptype, &c[0], &c[1], + &radius, &vol, &surftens, &gb_radius); + if (nread < 8) + { + too_few(wi); + return; + } + } + else + { + /* have_atomic_number && !have_bonded_type */ + nread = sscanf(line, "%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf", + type, &atomnr, &m, &q, ptype, &c[0], &c[1], + &radius, &vol, &surftens, &gb_radius); + if (nread < 7) + { + too_few(wi); + return; + } + } + } + else + { + if (have_bonded_type) + { + /* !have_atomic_number && have_bonded_type */ + nread = sscanf(line, "%s%s%lf%lf%s%lf%lf%lf%lf%lf%lf", + type, btype, &m, &q, ptype, &c[0], &c[1], + &radius, &vol, &surftens, &gb_radius); + if (nread < 7) + { + too_few(wi); + return; + } + } + else + { + /* !have_atomic_number && !have_bonded_type */ + nread = sscanf(line, "%s%lf%lf%s%lf%lf%lf%lf%lf%lf", + type, &m, &q, ptype, &c[0], &c[1], + &radius, &vol, &surftens, &gb_radius); + if (nread < 6) + { + too_few(wi); + return; + } + } + } + + if (!have_bonded_type) + { + strcpy(btype, type); + } + + if (!have_atomic_number) + { + atomnr = -1; + } + + break; + + case F_BHAM: + nfp0 = 3; + + if (have_atomic_number) + { + if (have_bonded_type) + { + nread = sscanf(line, "%s%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf%lf", + type, btype, &atomnr, &m, &q, ptype, &c[0], &c[1], &c[2], + &radius, &vol, &surftens, &gb_radius); + if (nread < 9) + { + too_few(wi); + return; + } + } + else + { + /* have_atomic_number && !have_bonded_type */ + nread = sscanf(line, "%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf%lf", + type, &atomnr, &m, &q, ptype, &c[0], &c[1], &c[2], + &radius, &vol, &surftens, &gb_radius); + if (nread < 8) + { + too_few(wi); + return; + } + } + } + else + { + if (have_bonded_type) + { + /* !have_atomic_number && have_bonded_type */ + nread = sscanf(line, "%s%s%lf%lf%s%lf%lf%lf%lf%lf%lf%lf", + type, btype, &m, &q, ptype, &c[0], &c[1], &c[2], + &radius, &vol, &surftens, &gb_radius); + if (nread < 8) + { + too_few(wi); + return; + } + } + else + { + /* !have_atomic_number && !have_bonded_type */ + nread = sscanf(line, "%s%lf%lf%s%lf%lf%lf%lf%lf%lf%lf", + type, &m, &q, ptype, &c[0], &c[1], &c[2], + &radius, &vol, &surftens, &gb_radius); + if (nread < 7) + { + too_few(wi); + return; + } + } + } + + if (!have_bonded_type) + { + strcpy(btype, type); + } + + if (!have_atomic_number) + { + atomnr = -1; + } + + break; + + default: + gmx_fatal(FARGS, "Invalid function type %d in push_at %s %d", nb_funct, + __FILE__, __LINE__); + } + for (j = nfp0; (j < MAXFORCEPARAM); j++) + { + c[j] = 0.0; + } + + if (strlen(type) == 1 && isdigit(type[0])) + { + gmx_fatal(FARGS, "Atom type names can't be single digits."); + } + + if (strlen(btype) == 1 && isdigit(btype[0])) + { + gmx_fatal(FARGS, "Bond atom type names can't be single digits."); + } + + /* Hack to read old topologies */ + if (gmx_strcasecmp(ptype, "D") == 0) + { + sprintf(ptype, "V"); + } + for (j = 0; (j < eptNR); j++) + { + if (gmx_strcasecmp(ptype, xl[j].entry) == 0) + { + break; + } + } + if (j == eptNR) + { + gmx_fatal(FARGS, "Invalid particle type %s on line %s", + ptype, line); + } + pt = xl[j].ptype; + if (debug) + { + fprintf(debug, "ptype: %s\n", ptype_str[pt]); + } + + atom->q = q; + atom->m = m; + atom->ptype = pt; + for (i = 0; (i < MAXFORCEPARAM); i++) + { + param->c[i] = c[i]; + } + + if ((batype_nr = get_bond_atomtype_type(btype, bat)) == NOTSET) + { + add_bond_atomtype(bat, symtab, btype); + } + batype_nr = get_bond_atomtype_type(btype, bat); + + if ((nr = get_atomtype_type(type, at)) != NOTSET) + { + sprintf(errbuf, "Overriding atomtype %s", type); + warning(wi, errbuf); + if ((nr = set_atomtype(nr, at, symtab, atom, type, param, batype_nr, + radius, vol, surftens, atomnr, gb_radius, S_hct)) == NOTSET) + { + gmx_fatal(FARGS, "Replacing atomtype %s failed", type); + } + } + else if ((nr = add_atomtype(at, symtab, atom, type, param, + batype_nr, radius, vol, + surftens, atomnr, gb_radius, S_hct)) == NOTSET) + { + gmx_fatal(FARGS, "Adding atomtype %s failed", type); + } + else + { + /* Add space in the non-bonded parameters matrix */ + realloc_nb_params(at, nbparam, pair); + } + sfree(atom); + sfree(param); +} + +static void push_bondtype(t_params * bt, + t_param * b, + int nral, + int ftype, + gmx_bool bAllowRepeat, + char * line, + warninp_t wi) +{ + int i, j; + gmx_bool bTest, bFound, bCont, bId; + int nr = bt->nr; + int nrfp = NRFP(ftype); + char errbuf[256]; + + /* If bAllowRepeat is TRUE, we allow multiple entries as long as they + are on directly _adjacent_ lines. + */ + + /* First check if our atomtypes are _identical_ (not reversed) to the previous + entry. If they are not identical we search for earlier duplicates. If they are + we can skip it, since we already searched for the first line + in this group. + */ + + bFound = FALSE; + bCont = FALSE; + + if (bAllowRepeat && nr > 1) + { + for (j = 0, bCont = TRUE; (j < nral); j++) + { + bCont = bCont && (b->a[j] == bt->param[nr-2].a[j]); + } + } + + /* Search for earlier duplicates if this entry was not a continuation + from the previous line. + */ + if (!bCont) + { + bFound = FALSE; + for (i = 0; (i < nr); i++) + { + bTest = TRUE; + for (j = 0; (j < nral); j++) + { + bTest = (bTest && (b->a[j] == bt->param[i].a[j])); + } + if (!bTest) + { + bTest = TRUE; + for (j = 0; (j < nral); j++) + { + bTest = (bTest && (b->a[nral-1-j] == bt->param[i].a[j])); + } + } + if (bTest) + { + if (!bFound) + { + bId = TRUE; + for (j = 0; (j < nrfp); j++) + { + bId = bId && (bt->param[i].c[j] == b->c[j]); + } + if (!bId) + { + sprintf(errbuf, "Overriding %s parameters.%s", + interaction_function[ftype].longname, - (ftype == F_PDIHS) ? "\nUse dihedraltype 4 to allow several multiplicity terms." : ""); ++ (ftype == F_PDIHS) ? ++ "\nUse dihedraltype 9 to allow several multiplicity terms. Only consecutive lines are combined. Non-consective lines overwrite each other." ++ : ""); + warning(wi, errbuf); - fprintf(stderr, " old:"); ++ fprintf(stderr, " old: "); + for (j = 0; (j < nrfp); j++) + { + fprintf(stderr, " %g", bt->param[i].c[j]); + } + fprintf(stderr, " \n new: %s\n\n", line); + } + } + /* Overwrite it! */ + for (j = 0; (j < nrfp); j++) + { + bt->param[i].c[j] = b->c[j]; + } + bFound = TRUE; + } + } + } + if (!bFound) + { + /* alloc */ + pr_alloc (2, bt); + + /* fill the arrays up and down */ + memcpy(bt->param[bt->nr].c, b->c, sizeof(b->c)); + memcpy(bt->param[bt->nr].a, b->a, sizeof(b->a)); + memcpy(bt->param[bt->nr+1].c, b->c, sizeof(b->c)); + + /* The definitions of linear angles depend on the order of atoms, + * that means that for atoms i-j-k, with certain parameter a, the + * corresponding k-j-i angle will have parameter 1-a. + */ + if (ftype == F_LINEAR_ANGLES) + { + bt->param[bt->nr+1].c[0] = 1-bt->param[bt->nr+1].c[0]; + bt->param[bt->nr+1].c[2] = 1-bt->param[bt->nr+1].c[2]; + } + + for (j = 0; (j < nral); j++) + { + bt->param[bt->nr+1].a[j] = b->a[nral-1-j]; + } + + bt->nr += 2; + } +} + +void push_bt(directive d, t_params bt[], int nral, + gpp_atomtype_t at, + t_bond_atomtype bat, char *line, + warninp_t wi) +{ + const char *formal[MAXATOMLIST+1] = { + "%s", + "%s%s", + "%s%s%s", + "%s%s%s%s", + "%s%s%s%s%s", + "%s%s%s%s%s%s", + "%s%s%s%s%s%s%s" + }; + const char *formnl[MAXATOMLIST+1] = { + "%*s", + "%*s%*s", + "%*s%*s%*s", + "%*s%*s%*s%*s", + "%*s%*s%*s%*s%*s", + "%*s%*s%*s%*s%*s%*s", + "%*s%*s%*s%*s%*s%*s%*s" + }; + const char *formlf = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf"; + int i, ft, ftype, nn, nrfp, nrfpA, nrfpB; + char f1[STRLEN]; + char alc[MAXATOMLIST+1][20]; + /* One force parameter more, so we can check if we read too many */ + double c[MAXFORCEPARAM+1]; + t_param p; + char errbuf[256]; + + if ((bat && at) || (!bat && !at)) + { + gmx_incons("You should pass either bat or at to push_bt"); + } + + /* Make format string (nral ints+functype) */ + if ((nn = sscanf(line, formal[nral], + alc[0], alc[1], alc[2], alc[3], alc[4], alc[5])) != nral+1) + { + sprintf(errbuf, "Not enough atomtypes (%d instead of %d)", nn-1, nral); + warning_error(wi, errbuf); + return; + } + + ft = strtol(alc[nral], NULL, 10); + ftype = ifunc_index(d, ft); + nrfp = NRFP(ftype); + nrfpA = interaction_function[ftype].nrfpA; + nrfpB = interaction_function[ftype].nrfpB; + strcpy(f1, formnl[nral]); + strcat(f1, formlf); + 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])) + != nrfp) + { + if (nn == nrfpA) + { + /* Copy the B-state from the A-state */ + copy_B_from_A(ftype, c); + } + else + { + if (nn < nrfpA) + { + warning_error(wi, "Not enough parameters"); + } + else if (nn > nrfpA && nn < nrfp) + { + warning_error(wi, "Too many parameters or not enough parameters for topology B"); + } + else if (nn > nrfp) + { + warning_error(wi, "Too many parameters"); + } + for (i = nn; (i < nrfp); i++) + { + c[i] = 0.0; + } + } + } + for (i = 0; (i < nral); i++) + { + if (at && ((p.a[i] = get_atomtype_type(alc[i], at)) == NOTSET)) + { + gmx_fatal(FARGS, "Unknown atomtype %s\n", alc[i]); + } + else if (bat && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET)) + { + gmx_fatal(FARGS, "Unknown bond_atomtype %s\n", alc[i]); + } + } + for (i = 0; (i < nrfp); i++) + { + p.c[i] = c[i]; + } + push_bondtype (&(bt[ftype]), &p, nral, ftype, FALSE, line, wi); +} + + +void push_dihedraltype(directive d, t_params bt[], + t_bond_atomtype bat, char *line, + warninp_t wi) +{ + const char *formal[MAXATOMLIST+1] = { + "%s", + "%s%s", + "%s%s%s", + "%s%s%s%s", + "%s%s%s%s%s", + "%s%s%s%s%s%s", + "%s%s%s%s%s%s%s" + }; + const char *formnl[MAXATOMLIST+1] = { + "%*s", + "%*s%*s", + "%*s%*s%*s", + "%*s%*s%*s%*s", + "%*s%*s%*s%*s%*s", + "%*s%*s%*s%*s%*s%*s", + "%*s%*s%*s%*s%*s%*s%*s" + }; + const char *formlf[MAXFORCEPARAM] = { + "%lf", + "%lf%lf", + "%lf%lf%lf", + "%lf%lf%lf%lf", + "%lf%lf%lf%lf%lf", + "%lf%lf%lf%lf%lf%lf", + "%lf%lf%lf%lf%lf%lf%lf", + "%lf%lf%lf%lf%lf%lf%lf%lf", + "%lf%lf%lf%lf%lf%lf%lf%lf%lf", + "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf", + "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf", + "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf", + }; + int i, ft, ftype, nn, nrfp, nrfpA, nrfpB, nral; + char f1[STRLEN]; + char alc[MAXATOMLIST+1][20]; + double c[MAXFORCEPARAM]; + t_param p; + gmx_bool bAllowRepeat; + char errbuf[256]; + + /* This routine accepts dihedraltypes defined from either 2 or 4 atoms. + * + * We first check for 2 atoms with the 3th column being an integer + * defining the type. If this isn't the case, we try it with 4 atoms + * and the 5th column defining the dihedral type. + */ + nn = sscanf(line, formal[4], alc[0], alc[1], alc[2], alc[3], alc[4]); + if (nn >= 3 && strlen(alc[2]) == 1 && isdigit(alc[2][0])) + { + nral = 2; + ft = strtol(alc[nral], NULL, 10); + /* Move atom types around a bit and use 'X' for wildcard atoms + * to create a 4-atom dihedral definition with arbitrary atoms in + * position 1 and 4. + */ + if (alc[2][0] == '2') + { + /* improper - the two atomtypes are 1,4. Use wildcards for 2,3 */ + strcpy(alc[3], alc[1]); + sprintf(alc[2], "X"); + sprintf(alc[1], "X"); + /* alc[0] stays put */ + } + else + { + /* proper - the two atomtypes are 2,3. Use wildcards for 1,4 */ + sprintf(alc[3], "X"); + strcpy(alc[2], alc[1]); + strcpy(alc[1], alc[0]); + sprintf(alc[0], "X"); + } + } + else if (nn == 5 && strlen(alc[4]) == 1 && isdigit(alc[4][0])) + { + nral = 4; + ft = strtol(alc[nral], NULL, 10); + } + else + { + sprintf(errbuf, "Incorrect number of atomtypes for dihedral (%d instead of 2 or 4)", nn-1); + warning_error(wi, errbuf); + return; + } + + if (ft == 9) + { + /* Previously, we have always overwritten parameters if e.g. a torsion + with the same atomtypes occurs on multiple lines. However, CHARMM and + some other force fields specify multiple dihedrals over some bonds, + including cosines with multiplicity 6 and somethimes even higher. + Thus, they cannot be represented with Ryckaert-Bellemans terms. + To add support for these force fields, Dihedral type 9 is identical to + normal proper dihedrals, but repeated entries are allowed. + */ + bAllowRepeat = TRUE; + ft = 1; + } + else + { + bAllowRepeat = FALSE; + } + + + ftype = ifunc_index(d, ft); + nrfp = NRFP(ftype); + nrfpA = interaction_function[ftype].nrfpA; + nrfpB = interaction_function[ftype].nrfpB; + + strcpy(f1, formnl[nral]); + strcat(f1, formlf[nrfp-1]); + + /* Check number of parameters given */ + 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])) + != nrfp) + { + if (nn == nrfpA) + { + /* Copy the B-state from the A-state */ + copy_B_from_A(ftype, c); + } + else + { + if (nn < nrfpA) + { + warning_error(wi, "Not enough parameters"); + } + else if (nn > nrfpA && nn < nrfp) + { + warning_error(wi, "Too many parameters or not enough parameters for topology B"); + } + else if (nn > nrfp) + { + warning_error(wi, "Too many parameters"); + } + for (i = nn; (i < nrfp); i++) + { + c[i] = 0.0; + } + } + } + + for (i = 0; (i < 4); i++) + { + if (!strcmp(alc[i], "X")) + { + p.a[i] = -1; + } + else + { + if ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET) + { + gmx_fatal(FARGS, "Unknown bond_atomtype %s", alc[i]); + } + } + } + for (i = 0; (i < nrfp); i++) + { + p.c[i] = c[i]; + } + /* Always use 4 atoms here, since we created two wildcard atoms + * if there wasn't of them 4 already. + */ + push_bondtype (&(bt[ftype]), &p, 4, ftype, bAllowRepeat, line, wi); +} + + +void push_nbt(directive d, t_nbparam **nbt, gpp_atomtype_t atype, + char *pline, int nb_funct, + warninp_t wi) +{ + /* swap the atoms */ + const char *form2 = "%*s%*s%*s%lf%lf"; + const char *form3 = "%*s%*s%*s%lf%lf%lf"; + const char *form4 = "%*s%*s%*s%lf%lf%lf%lf"; + const char *form5 = "%*s%*s%*s%lf%lf%lf%lf%lf"; + char a0[80], a1[80]; + int i, f, n, ftype, atnr, nrfp; + double c[4], dum; + real cr[4], sig6; + atom_id ai, aj; + t_nbparam *nbp; + gmx_bool bId; + char errbuf[256]; + + if (sscanf (pline, "%s%s%d", a0, a1, &f) != 3) + { + too_few(wi); + return; + } + + ftype = ifunc_index(d, f); + + if (ftype != nb_funct) + { + sprintf(errbuf, "Trying to add %s while the default nonbond type is %s", + interaction_function[ftype].longname, + interaction_function[nb_funct].longname); + warning_error(wi, errbuf); + return; + } + + /* Get the force parameters */ + nrfp = NRFP(ftype); + if (ftype == F_LJ14) + { + n = sscanf(pline, form4, &c[0], &c[1], &c[2], &c[3]); + if (n < 2) + { + too_few(wi); + return; + } + /* When the B topology parameters are not set, + * copy them from topology A + */ + assert(nrfp <= 4); + for (i = n; i < nrfp; i++) + { + c[i] = c[i-2]; + } + } + else if (ftype == F_LJC14_Q) + { + n = sscanf(pline, form5, &c[0], &c[1], &c[2], &c[3], &dum); + if (n != 4) + { + incorrect_n_param(wi); + return; + } + } + else if (nrfp == 2) + { + if (sscanf(pline, form3, &c[0], &c[1], &dum) != 2) + { + incorrect_n_param(wi); + return; + } + } + else if (nrfp == 3) + { + if (sscanf(pline, form4, &c[0], &c[1], &c[2], &dum) != 3) + { + incorrect_n_param(wi); + return; + } + } + else + { + gmx_fatal(FARGS, "Number of force parameters for nonbonded interactions is %d" + " in file %s, line %d", nrfp, __FILE__, __LINE__); + } + for (i = 0; (i < nrfp); i++) + { + cr[i] = c[i]; + } + + /* Put the parameters in the matrix */ + if ((ai = get_atomtype_type (a0, atype)) == NOTSET) + { + gmx_fatal(FARGS, "Atomtype %s not found", a0); + } + if ((aj = get_atomtype_type (a1, atype)) == NOTSET) + { + gmx_fatal(FARGS, "Atomtype %s not found", a1); + } + nbp = &(nbt[max(ai, aj)][min(ai, aj)]); + + if (nbp->bSet) + { + bId = TRUE; + for (i = 0; i < nrfp; i++) + { + bId = bId && (nbp->c[i] == cr[i]); + } + if (!bId) + { + sprintf(errbuf, "Overriding non-bonded parameters,"); + warning(wi, errbuf); + fprintf(stderr, " old:"); + for (i = 0; i < nrfp; i++) + { + fprintf(stderr, " %g", nbp->c[i]); + } + fprintf(stderr, " new\n%s\n", pline); + } + } + nbp->bSet = TRUE; + for (i = 0; i < nrfp; i++) + { + nbp->c[i] = cr[i]; + } +} + +void +push_gb_params (gpp_atomtype_t at, char *line, + warninp_t wi) +{ + int nfield; + int atype; + double radius, vol, surftens, gb_radius, S_hct; + char atypename[STRLEN]; + char errbuf[STRLEN]; + + if ( (nfield = sscanf(line, "%s%lf%lf%lf%lf%lf", atypename, &radius, &vol, &surftens, &gb_radius, &S_hct)) != 6) + { + sprintf(errbuf, "Too few gb parameters for type %s\n", atypename); + warning(wi, errbuf); + } + + /* Search for atomtype */ + atype = get_atomtype_type(atypename, at); + + if (atype == NOTSET) + { + printf("Couldn't find topology match for atomtype %s\n", atypename); + abort(); + } + + set_atomtype_gbparam(at, atype, radius, vol, surftens, gb_radius, S_hct); +} + +void +push_cmaptype(directive d, t_params bt[], int nral, gpp_atomtype_t at, + t_bond_atomtype bat, char *line, + warninp_t wi) +{ + const char *formal = "%s%s%s%s%s%s%s%s"; + + int i, j, ft, ftype, nn, nrfp, nrfpA, nrfpB; + int start; + int nxcmap, nycmap, ncmap, read_cmap, sl, nct; + char s[20], alc[MAXATOMLIST+2][20]; + t_param p; + gmx_bool bAllowRepeat; + char errbuf[256]; + + /* Keep the compiler happy */ + read_cmap = 0; + start = 0; + + if ((nn = sscanf(line, formal, alc[0], alc[1], alc[2], alc[3], alc[4], alc[5], alc[6], alc[7])) != nral+3) + { + sprintf(errbuf, "Incorrect number of atomtypes for cmap (%d instead of 5)", nn-1); + warning_error(wi, errbuf); + return; + } + + /* Compute an offset for each line where the cmap parameters start + * ie. where the atom types and grid spacing information ends + */ + for (i = 0; i < nn; i++) + { + start += (int)strlen(alc[i]); + } + + /* There are nn-1 spaces between the atom types and the grid spacing info in the cmap.itp file */ + /* start is the position on the line where we start to read the actual cmap grid data from the itp file */ + start = start + nn -1; + + ft = strtol(alc[nral], NULL, 10); + nxcmap = strtol(alc[nral+1], NULL, 10); + nycmap = strtol(alc[nral+2], NULL, 10); + + /* Check for equal grid spacing in x and y dims */ + if (nxcmap != nycmap) + { + gmx_fatal(FARGS, "Not the same grid spacing in x and y for cmap grid: x=%d, y=%d", nxcmap, nycmap); + } + + ncmap = nxcmap*nycmap; + ftype = ifunc_index(d, ft); + nrfpA = strtol(alc[6], NULL, 10)*strtol(alc[6], NULL, 10); + nrfpB = strtol(alc[7], NULL, 10)*strtol(alc[7], NULL, 10); + nrfp = nrfpA+nrfpB; + + /* Allocate memory for the CMAP grid */ + bt[F_CMAP].ncmap += nrfp; + srenew(bt[F_CMAP].cmap, bt[F_CMAP].ncmap); + + /* Read in CMAP parameters */ + sl = 0; + for (i = 0; i < ncmap; i++) + { + while (isspace(*(line+start+sl))) + { + sl++; + } + nn = sscanf(line+start+sl, " %s ", s); + sl += strlen(s); + bt[F_CMAP].cmap[i+(bt[F_CMAP].ncmap)-nrfp] = strtod(s, NULL); + + if (nn == 1) + { + read_cmap++; + } + else + { + 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]); + } + + } + + /* Check do that we got the number of parameters we expected */ + if (read_cmap == nrfpA) + { + for (i = 0; i < ncmap; i++) + { + bt[F_CMAP].cmap[i+ncmap] = bt[F_CMAP].cmap[i]; + } + } + else + { + if (read_cmap < nrfpA) + { + warning_error(wi, "Not enough cmap parameters"); + } + else if (read_cmap > nrfpA && read_cmap < nrfp) + { + warning_error(wi, "Too many cmap parameters or not enough parameters for topology B"); + } + else if (read_cmap > nrfp) + { + warning_error(wi, "Too many cmap parameters"); + } + } + + + /* Set grid spacing and the number of grids (we assume these numbers to be the same for all grids + * so we can safely assign them each time + */ + bt[F_CMAP].grid_spacing = nxcmap; /* Or nycmap, they need to be equal */ + bt[F_CMAP].nc = bt[F_CMAP].nc + 1; /* Since we are incrementing here, we need to subtract later, see (*****) */ + nct = (nral+1) * bt[F_CMAP].nc; + + /* Allocate memory for the cmap_types information */ + srenew(bt[F_CMAP].cmap_types, nct); + + for (i = 0; (i < nral); i++) + { + if (at && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET)) + { + gmx_fatal(FARGS, "Unknown atomtype %s\n", alc[i]); + } + else if (bat && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET)) + { + gmx_fatal(FARGS, "Unknown bond_atomtype %s\n", alc[i]); + } + + /* Assign a grid number to each cmap_type */ + bt[F_CMAP].cmap_types[bt[F_CMAP].nct++] = get_bond_atomtype_type(alc[i], bat); + } + + /* Assign a type number to this cmap */ + 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 (****) */ + + /* Check for the correct number of atoms (again) */ + if (bt[F_CMAP].nct != nct) + { + gmx_fatal(FARGS, "Incorrect number of atom types (%d) in cmap type %d\n", nct, bt[F_CMAP].nc); + } + + /* Is this correct?? */ + for (i = 0; i < MAXFORCEPARAM; i++) + { + p.c[i] = NOTSET; + } + + /* Push the bond to the bondlist */ + push_bondtype (&(bt[ftype]), &p, nral, ftype, FALSE, line, wi); +} + + +static void push_atom_now(t_symtab *symtab, t_atoms *at, int atomnr, + int atomicnumber, + int type, char *ctype, int ptype, + char *resnumberic, + char *resname, char *name, real m0, real q0, + int typeB, char *ctypeB, real mB, real qB) +{ + int j, resind = 0, resnr; + unsigned char ric; + int nr = at->nr; + + if (((nr == 0) && (atomnr != 1)) || (nr && (atomnr != at->nr+1))) + { + gmx_fatal(FARGS, "Atoms in the .top are not numbered consecutively from 1 (rather, atomnr = %d, while at->nr = %d)", atomnr, at->nr); + } + + j = strlen(resnumberic) - 1; + if (isdigit(resnumberic[j])) + { + ric = ' '; + } + else + { + ric = resnumberic[j]; + if (j == 0 || !isdigit(resnumberic[j-1])) + { + gmx_fatal(FARGS, "Invalid residue number '%s' for atom %d", + resnumberic, atomnr); + } + } + resnr = strtol(resnumberic, NULL, 10); + + if (nr > 0) + { + resind = at->atom[nr-1].resind; + } + if (nr == 0 || strcmp(resname, *at->resinfo[resind].name) != 0 || + resnr != at->resinfo[resind].nr || + ric != at->resinfo[resind].ic) + { + if (nr == 0) + { + resind = 0; + } + else + { + resind++; + } + at->nres = resind + 1; + srenew(at->resinfo, at->nres); + at->resinfo[resind].name = put_symtab(symtab, resname); + at->resinfo[resind].nr = resnr; + at->resinfo[resind].ic = ric; + } + else + { + resind = at->atom[at->nr-1].resind; + } + + /* New atom instance + * get new space for arrays + */ + srenew(at->atom, nr+1); + srenew(at->atomname, nr+1); + srenew(at->atomtype, nr+1); + srenew(at->atomtypeB, nr+1); + + /* fill the list */ + at->atom[nr].type = type; + at->atom[nr].ptype = ptype; + at->atom[nr].q = q0; + at->atom[nr].m = m0; + at->atom[nr].typeB = typeB; + at->atom[nr].qB = qB; + at->atom[nr].mB = mB; + + at->atom[nr].resind = resind; + at->atom[nr].atomnumber = atomicnumber; + at->atomname[nr] = put_symtab(symtab, name); + at->atomtype[nr] = put_symtab(symtab, ctype); + at->atomtypeB[nr] = put_symtab(symtab, ctypeB); + at->nr++; +} + +void push_cg(t_block *block, int *lastindex, int index, int a) +{ + if (debug) + { + fprintf (debug, "Index %d, Atom %d\n", index, a); + } + + if (((block->nr) && (*lastindex != index)) || (!block->nr)) + { + /* add a new block */ + block->nr++; + srenew(block->index, block->nr+1); + } + block->index[block->nr] = a + 1; + *lastindex = index; +} + +void push_atom(t_symtab *symtab, t_block *cgs, + t_atoms *at, gpp_atomtype_t atype, char *line, int *lastcg, + warninp_t wi) +{ + int nr, ptype; + int cgnumber, atomnr, type, typeB, nscan; + char id[STRLEN], ctype[STRLEN], ctypeB[STRLEN], + resnumberic[STRLEN], resname[STRLEN], name[STRLEN], check[STRLEN]; + double m, q, mb, qb; + real m0, q0, mB, qB; + + /* Make a shortcut for writing in this molecule */ + nr = at->nr; + + /* Fixed parameters */ + if (sscanf(line, "%s%s%s%s%s%d", + id, ctype, resnumberic, resname, name, &cgnumber) != 6) + { + too_few(wi); + return; + } + sscanf(id, "%d", &atomnr); + if ((type = get_atomtype_type(ctype, atype)) == NOTSET) + { + gmx_fatal(FARGS, "Atomtype %s not found", ctype); + } + ptype = get_atomtype_ptype(type, atype); + + /* Set default from type */ + q0 = get_atomtype_qA(type, atype); + m0 = get_atomtype_massA(type, atype); + typeB = type; + qB = q0; + mB = m0; + + /* Optional parameters */ + nscan = sscanf(line, "%*s%*s%*s%*s%*s%*s%lf%lf%s%lf%lf%s", + &q, &m, ctypeB, &qb, &mb, check); + + /* Nasty switch that falls thru all the way down! */ + if (nscan > 0) + { + q0 = qB = q; + if (nscan > 1) + { + m0 = mB = m; + if (nscan > 2) + { + if ((typeB = get_atomtype_type(ctypeB, atype)) == NOTSET) + { + gmx_fatal(FARGS, "Atomtype %s not found", ctypeB); + } + qB = get_atomtype_qA(typeB, atype); + mB = get_atomtype_massA(typeB, atype); + if (nscan > 3) + { + qB = qb; + if (nscan > 4) + { + mB = mb; + if (nscan > 5) + { + warning_error(wi, "Too many parameters"); + } + } + } + } + } + } + if (debug) + { + fprintf(debug, "mB=%g, qB=%g, typeB=%d\n", mB, qB, typeB); + } + + push_cg(cgs, lastcg, cgnumber, nr); + + push_atom_now(symtab, at, atomnr, get_atomtype_atomnumber(type, atype), + type, ctype, ptype, resnumberic, + resname, name, m0, q0, typeB, + typeB == type ? ctype : ctypeB, mB, qB); +} + +void push_molt(t_symtab *symtab, int *nmol, t_molinfo **mol, char *line, + warninp_t wi) +{ + char type[STRLEN]; + int nrexcl, i; + t_molinfo *newmol; + + if ((sscanf(line, "%s%d", type, &nrexcl)) != 2) + { + warning_error(wi, "Expected a molecule type name and nrexcl"); + } + + /* Test if this atomtype overwrites another */ + i = 0; + while (i < *nmol) + { + if (gmx_strcasecmp(*((*mol)[i].name), type) == 0) + { + gmx_fatal(FARGS, "moleculetype %s is redefined", type); + } + i++; + } + + (*nmol)++; + srenew(*mol, *nmol); + newmol = &((*mol)[*nmol-1]); + init_molinfo(newmol); + + /* Fill in the values */ + newmol->name = put_symtab(symtab, type); + newmol->nrexcl = nrexcl; + newmol->excl_set = FALSE; +} + +static gmx_bool default_nb_params(int ftype, t_params bt[], t_atoms *at, + t_param *p, int c_start, gmx_bool bB, gmx_bool bGenPairs) +{ + int i, j, ti, tj, ntype; + gmx_bool bFound; + t_param *pi = NULL; + int nr = bt[ftype].nr; + int nral = NRAL(ftype); + int nrfp = interaction_function[ftype].nrfpA; + int nrfpB = interaction_function[ftype].nrfpB; + + if ((!bB && nrfp == 0) || (bB && nrfpB == 0)) + { + return TRUE; + } + + bFound = FALSE; + if (bGenPairs) + { + /* First test the generated-pair position to save + * time when we have 1000*1000 entries for e.g. OPLS... + */ + ntype = sqrt(nr); + if (bB) + { + ti = at->atom[p->a[0]].typeB; + tj = at->atom[p->a[1]].typeB; + } + else + { + ti = at->atom[p->a[0]].type; + tj = at->atom[p->a[1]].type; + } + pi = &(bt[ftype].param[ntype*ti+tj]); + bFound = ((ti == pi->a[0]) && (tj == pi->a[1])); + } + + /* Search explicitly if we didnt find it */ + if (!bFound) + { + for (i = 0; ((i < nr) && !bFound); i++) + { + pi = &(bt[ftype].param[i]); + if (bB) + { + for (j = 0; ((j < nral) && + (at->atom[p->a[j]].typeB == pi->a[j])); j++) + { + ; + } + } + else + { + for (j = 0; ((j < nral) && + (at->atom[p->a[j]].type == pi->a[j])); j++) + { + ; + } + } + bFound = (j == nral); + } + } + + if (bFound) + { + if (bB) + { + if (nrfp+nrfpB > MAXFORCEPARAM) + { + gmx_incons("Too many force parameters"); + } + for (j = c_start; (j < nrfpB); j++) + { + p->c[nrfp+j] = pi->c[j]; + } + } + else + { + for (j = c_start; (j < nrfp); j++) + { + p->c[j] = pi->c[j]; + } + } + } + else + { + for (j = c_start; (j < nrfp); j++) + { + p->c[j] = 0.0; + } + } + return bFound; +} + +static gmx_bool default_cmap_params(t_params bondtype[], + t_atoms *at, gpp_atomtype_t atype, + t_param *p, gmx_bool bB, + int *cmap_type, int *nparam_def) +{ + int i, j, nparam_found; + int ct; + gmx_bool bFound = FALSE; + + nparam_found = 0; + ct = 0; + + /* Match the current cmap angle against the list of cmap_types */ + for (i = 0; i < bondtype[F_CMAP].nct && !bFound; i += 6) + { + if (bB) + { + + } + else + { + if ( + (get_atomtype_batype(at->atom[p->a[0]].type, atype) == bondtype[F_CMAP].cmap_types[i]) && + (get_atomtype_batype(at->atom[p->a[1]].type, atype) == bondtype[F_CMAP].cmap_types[i+1]) && + (get_atomtype_batype(at->atom[p->a[2]].type, atype) == bondtype[F_CMAP].cmap_types[i+2]) && + (get_atomtype_batype(at->atom[p->a[3]].type, atype) == bondtype[F_CMAP].cmap_types[i+3]) && + (get_atomtype_batype(at->atom[p->a[4]].type, atype) == bondtype[F_CMAP].cmap_types[i+4])) + { + /* Found cmap torsion */ + bFound = TRUE; + ct = bondtype[F_CMAP].cmap_types[i+5]; + nparam_found = 1; + } + } + } + + /* If we did not find a matching type for this cmap torsion */ + if (!bFound) + { + gmx_fatal(FARGS, "Unknown cmap torsion between atoms %d %d %d %d %d\n", + p->a[0]+1, p->a[1]+1, p->a[2]+1, p->a[3]+1, p->a[4]+1); + } + + *nparam_def = nparam_found; + *cmap_type = ct; + + return bFound; +} + +static gmx_bool default_params(int ftype, t_params bt[], + t_atoms *at, gpp_atomtype_t atype, + t_param *p, gmx_bool bB, + t_param **param_def, + int *nparam_def) +{ + int i, j, nparam_found; + gmx_bool bFound, bSame; + t_param *pi = NULL; + t_param *pj = NULL; + int nr = bt[ftype].nr; + int nral = NRAL(ftype); + int nrfpA = interaction_function[ftype].nrfpA; + int nrfpB = interaction_function[ftype].nrfpB; + + if ((!bB && nrfpA == 0) || (bB && nrfpB == 0)) + { + return TRUE; + } + + + /* We allow wildcards now. The first type (with or without wildcards) that + * fits is used, so you should probably put the wildcarded bondtypes + * at the end of each section. + */ + bFound = FALSE; + nparam_found = 0; + /* OPLS uses 1000s of dihedraltypes, so in order to speed up the scanning we have a + * special case for this. Check for B state outside loop to speed it up. + */ + if (ftype == F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS || ftype == F_PIDIHS) + { + if (bB) + { + for (i = 0; ((i < nr) && !bFound); i++) + { + pi = &(bt[ftype].param[i]); + bFound = + ( + ((pi->AI == -1) || (get_atomtype_batype(at->atom[p->AI].typeB, atype) == pi->AI)) && + ((pi->AJ == -1) || (get_atomtype_batype(at->atom[p->AJ].typeB, atype) == pi->AJ)) && + ((pi->AK == -1) || (get_atomtype_batype(at->atom[p->AK].typeB, atype) == pi->AK)) && + ((pi->AL == -1) || (get_atomtype_batype(at->atom[p->AL].typeB, atype) == pi->AL)) + ); + } + } + else + { + /* State A */ + for (i = 0; ((i < nr) && !bFound); i++) + { + pi = &(bt[ftype].param[i]); + bFound = + ( + ((pi->AI == -1) || (get_atomtype_batype(at->atom[p->AI].type, atype) == pi->AI)) && + ((pi->AJ == -1) || (get_atomtype_batype(at->atom[p->AJ].type, atype) == pi->AJ)) && + ((pi->AK == -1) || (get_atomtype_batype(at->atom[p->AK].type, atype) == pi->AK)) && + ((pi->AL == -1) || (get_atomtype_batype(at->atom[p->AL].type, atype) == pi->AL)) + ); + } + } + /* Find additional matches for this dihedral - necessary for ftype==9 which is used e.g. for charmm. + * The rules in that case is that additional matches HAVE to be on adjacent lines! + */ + if (bFound == TRUE) + { + nparam_found++; + bSame = TRUE; + /* Continue from current i value */ + for (j = i+1; j < nr && bSame; j += 2) + { + pj = &(bt[ftype].param[j]); + bSame = (pi->AI == pj->AI && pi->AJ == pj->AJ && pi->AK == pj->AK && pi->AL == pj->AL); + if (bSame) + { + nparam_found++; + } + /* nparam_found will be increased as long as the numbers match */ + } + } + } + else /* Not a dihedral */ + { + for (i = 0; ((i < nr) && !bFound); i++) + { + pi = &(bt[ftype].param[i]); + if (bB) + { + for (j = 0; ((j < nral) && + (get_atomtype_batype(at->atom[p->a[j]].typeB, atype) == pi->a[j])); j++) + { + ; + } + } + else + { + for (j = 0; ((j < nral) && + (get_atomtype_batype(at->atom[p->a[j]].type, atype) == pi->a[j])); j++) + { + ; + } + } + bFound = (j == nral); + } + if (bFound) + { + nparam_found = 1; + } + } + + *param_def = pi; + *nparam_def = nparam_found; + + return bFound; +} + + + +void push_bond(directive d, t_params bondtype[], t_params bond[], + t_atoms *at, gpp_atomtype_t atype, char *line, + gmx_bool bBonded, gmx_bool bGenPairs, real fudgeQQ, + gmx_bool bZero, gmx_bool *bWarn_copy_A_B, + warninp_t wi) +{ + const char *aaformat[MAXATOMLIST] = { + "%d%d", + "%d%d%d", + "%d%d%d%d", + "%d%d%d%d%d", + "%d%d%d%d%d%d", + "%d%d%d%d%d%d%d" + }; + const char *asformat[MAXATOMLIST] = { + "%*s%*s", + "%*s%*s%*s", + "%*s%*s%*s%*s", + "%*s%*s%*s%*s%*s", + "%*s%*s%*s%*s%*s%*s", + "%*s%*s%*s%*s%*s%*s%*s" + }; + const char *ccformat = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf"; + int nr, i, j, nral, nral_fmt, nread, ftype; + char format[STRLEN]; + /* One force parameter more, so we can check if we read too many */ + double cc[MAXFORCEPARAM+1]; + int aa[MAXATOMLIST+1]; + t_param param, paramB, *param_defA, *param_defB; + gmx_bool bFoundA = FALSE, bFoundB = FALSE, bDef, bPert, bSwapParity = FALSE; + int nparam_defA, nparam_defB; + char errbuf[256]; + + nparam_defA = nparam_defB = 0; + + ftype = ifunc_index(d, 1); + nral = NRAL(ftype); + for (j = 0; j < MAXATOMLIST; j++) + { + aa[j] = NOTSET; + } + bDef = (NRFP(ftype) > 0); + + if (ftype == F_SETTLE) + { + /* SETTLE acts on 3 atoms, but the topology format only specifies + * the first atom (for historical reasons). + */ + nral_fmt = 1; + } + else + { + nral_fmt = nral; + } + + nread = sscanf(line, aaformat[nral_fmt-1], + &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]); + + if (ftype == F_SETTLE) + { + aa[3] = aa[1]; + aa[1] = aa[0] + 1; + aa[2] = aa[0] + 2; + } + + if (nread < nral_fmt) + { + too_few(wi); + return; + } + else if (nread > nral_fmt) + { + /* this is a hack to allow for virtual sites with swapped parity */ + bSwapParity = (aa[nral] < 0); + if (bSwapParity) + { + aa[nral] = -aa[nral]; + } + ftype = ifunc_index(d, aa[nral]); + if (bSwapParity) + { + switch (ftype) + { + case F_VSITE3FAD: + case F_VSITE3OUT: + break; + default: + gmx_fatal(FARGS, "Negative function types only allowed for %s and %s", + interaction_function[F_VSITE3FAD].longname, + interaction_function[F_VSITE3OUT].longname); + } + } + } + + + /* Check for double atoms and atoms out of bounds */ + for (i = 0; (i < nral); i++) + { + if (aa[i] < 1 || aa[i] > at->nr) + { + gmx_fatal(FARGS, "[ file %s, line %d ]:\n" + "Atom index (%d) in %s out of bounds (1-%d).\n" + "This probably means that you have inserted topology section \"%s\"\n" + "in a part belonging to a different molecule than you intended to.\n" + "In that case move the \"%s\" section to the right molecule.", + get_warning_file(wi), get_warning_line(wi), + aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d)); + } + for (j = i+1; (j < nral); j++) + { + if (aa[i] == aa[j]) + { + sprintf(errbuf, "Duplicate atom index (%d) in %s", aa[i], dir2str(d)); + warning(wi, errbuf); + } + } + } + + /* default force parameters */ + for (j = 0; (j < MAXATOMLIST); j++) + { + param.a[j] = aa[j]-1; + } + for (j = 0; (j < MAXFORCEPARAM); j++) + { + param.c[j] = 0.0; + } + + /* Get force params for normal and free energy perturbation + * studies, as determined by types! + */ + + if (bBonded) + { + bFoundA = default_params(ftype, bondtype, at, atype, ¶m, FALSE, ¶m_defA, &nparam_defA); + if (bFoundA) + { + /* Copy the A-state and B-state default parameters. */ + assert(NRFPA(ftype)+NRFPB(ftype) <= MAXFORCEPARAM); + for (j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++) + { + param.c[j] = param_defA->c[j]; + } + } + bFoundB = default_params(ftype, bondtype, at, atype, ¶m, TRUE, ¶m_defB, &nparam_defB); + if (bFoundB) + { + /* Copy only the B-state default parameters */ + for (j = NRFPA(ftype); (j < NRFP(ftype)); j++) + { + param.c[j] = param_defB->c[j]; + } + } + } + else if (ftype == F_LJ14) + { + bFoundA = default_nb_params(ftype, bondtype, at, ¶m, 0, FALSE, bGenPairs); + bFoundB = default_nb_params(ftype, bondtype, at, ¶m, 0, TRUE, bGenPairs); + } + else if (ftype == F_LJC14_Q) + { + param.c[0] = fudgeQQ; + /* Fill in the A-state charges as default parameters */ + param.c[1] = at->atom[param.a[0]].q; + param.c[2] = at->atom[param.a[1]].q; + /* The default LJ parameters are the standard 1-4 parameters */ + bFoundA = default_nb_params(F_LJ14, bondtype, at, ¶m, 3, FALSE, bGenPairs); + bFoundB = TRUE; + } + else if (ftype == F_LJC_PAIRS_NB) + { + /* Defaults are not supported here */ + bFoundA = FALSE; + bFoundB = TRUE; + } + else + { + gmx_incons("Unknown function type in push_bond"); + } + + if (nread > nral_fmt) + { + /* Manually specified parameters - in this case we discard multiple torsion info! */ + + strcpy(format, asformat[nral_fmt-1]); + strcat(format, ccformat); + + nread = sscanf(line, format, &cc[0], &cc[1], &cc[2], &cc[3], &cc[4], &cc[5], + &cc[6], &cc[7], &cc[8], &cc[9], &cc[10], &cc[11], &cc[12]); + + if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0)) + { + /* We only have to issue a warning if these atoms are perturbed! */ + bPert = FALSE; + for (j = 0; (j < nral); j++) + { + bPert = bPert || PERTURBED(at->atom[param.a[j]]); + } + + if (bPert && *bWarn_copy_A_B) + { + sprintf(errbuf, + "Some parameters for bonded interaction involving perturbed atoms are specified explicitly in state A, but not B - copying A to B"); + warning(wi, errbuf); + *bWarn_copy_A_B = FALSE; + } + + /* If only the A parameters were specified, copy them to the B state */ + /* The B-state parameters correspond to the first nrfpB + * A-state parameters. + */ + for (j = 0; (j < NRFPB(ftype)); j++) + { + cc[nread++] = cc[j]; + } + } + + /* If nread was 0 or EOF, no parameters were read => use defaults. + * If nread was nrfpA we copied above so nread=nrfp. + * If nread was nrfp we are cool. + * For F_LJC14_Q we allow supplying fudgeQQ only. + * Anything else is an error! + */ + if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) && + !(ftype == F_LJC14_Q && nread == 1)) + { + gmx_fatal(FARGS, "Incorrect number of parameters - found %d, expected %d or %d for %s.", + nread, NRFPA(ftype), NRFP(ftype), + interaction_function[ftype].longname); + } + + for (j = 0; (j < nread); j++) + { + param.c[j] = cc[j]; + } + + /* Check whether we have to use the defaults */ + if (nread == NRFP(ftype)) + { + bDef = FALSE; + } + } + else + { + nread = 0; + } + /* nread now holds the number of force parameters read! */ + + if (bDef) + { + /* Use defaults */ + /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */ + if (ftype == F_PDIHS) + { + if ((nparam_defA != nparam_defB) || ((nparam_defA > 1 || nparam_defB > 1) && (param_defA != param_defB))) + { + sprintf(errbuf, + "Cannot automatically perturb a torsion with multiple terms to different form.\n" + "Please specify perturbed parameters manually for this torsion in your topology!"); + warning_error(wi, errbuf); + } + } + + if (nread > 0 && nread < NRFPA(ftype)) + { + /* Issue an error, do not use defaults */ + sprintf(errbuf, "Not enough parameters, there should be at least %d (or 0 for defaults)", NRFPA(ftype)); + warning_error(wi, errbuf); + } + + if (nread == 0 || nread == EOF) + { + if (!bFoundA) + { + if (interaction_function[ftype].flags & IF_VSITE) + { + /* set them to NOTSET, will be calculated later */ + for (j = 0; (j < MAXFORCEPARAM); j++) + { + param.c[j] = NOTSET; + } + + if (bSwapParity) + { + param.C1 = -1; /* flag to swap parity of vsite construction */ + } + } + else + { + if (bZero) + { + fprintf(stderr, "NOTE: No default %s types, using zeroes\n", + interaction_function[ftype].longname); + } + else + { + sprintf(errbuf, "No default %s types", interaction_function[ftype].longname); + warning_error(wi, errbuf); + } + } + } + else + { + if (bSwapParity) + { + switch (ftype) + { + case F_VSITE3FAD: + param.C0 = 360 - param.C0; + break; + case F_VSITE3OUT: + param.C2 = -param.C2; + break; + } + } + } + if (!bFoundB) + { + /* We only have to issue a warning if these atoms are perturbed! */ + bPert = FALSE; + for (j = 0; (j < nral); j++) + { + bPert = bPert || PERTURBED(at->atom[param.a[j]]); + } + + if (bPert) + { + sprintf(errbuf, "No default %s types for perturbed atoms, " + "using normal values", interaction_function[ftype].longname); + warning(wi, errbuf); + } + } + } + } + + if ((ftype == F_PDIHS || ftype == F_ANGRES || ftype == F_ANGRESZ) + && param.c[5] != param.c[2]) + { + gmx_fatal(FARGS, "[ file %s, line %d ]:\n" + " %s multiplicity can not be perturbed %f!=%f", + get_warning_file(wi), get_warning_line(wi), + interaction_function[ftype].longname, + param.c[2], param.c[5]); + } + + if (IS_TABULATED(ftype) && param.c[2] != param.c[0]) + { + gmx_fatal(FARGS, "[ file %s, line %d ]:\n" + " %s table number can not be perturbed %d!=%d", + get_warning_file(wi), get_warning_line(wi), + interaction_function[ftype].longname, + (int)(param.c[0]+0.5), (int)(param.c[2]+0.5)); + } + + /* Dont add R-B dihedrals where all parameters are zero (no interaction) */ + if (ftype == F_RBDIHS) + { + nr = 0; + for (i = 0; i < NRFP(ftype); i++) + { + if (param.c[i] != 0) + { + nr++; + } + } + if (nr == 0) + { + return; + } + } + + /* Put the values in the appropriate arrays */ + add_param_to_list (&bond[ftype], ¶m); + + /* Push additional torsions from FF for ftype==9 if we have them. + * We have already checked that the A/B states do not differ in this case, + * so we do not have to double-check that again, or the vsite stuff. + * In addition, those torsions cannot be automatically perturbed. + */ + if (bDef && ftype == F_PDIHS) + { + for (i = 1; i < nparam_defA; i++) + { + /* Advance pointer! */ + param_defA += 2; + for (j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++) + { + param.c[j] = param_defA->c[j]; + } + /* And push the next term for this torsion */ + add_param_to_list (&bond[ftype], ¶m); + } + } +} + +void push_cmap(directive d, t_params bondtype[], t_params bond[], + t_atoms *at, gpp_atomtype_t atype, char *line, + warninp_t wi) +{ + const char *aaformat[MAXATOMLIST+1] = + { + "%d", + "%d%d", + "%d%d%d", + "%d%d%d%d", + "%d%d%d%d%d", + "%d%d%d%d%d%d", + "%d%d%d%d%d%d%d" + }; + + int i, j, nr, ftype, nral, nread, ncmap_params; + int cmap_type; + int aa[MAXATOMLIST+1]; + char errbuf[256]; + gmx_bool bFound; + t_param param, paramB, *param_defA, *param_defB; + + ftype = ifunc_index(d, 1); + nral = NRAL(ftype); + nr = bondtype[ftype].nr; + ncmap_params = 0; + + nread = sscanf(line, aaformat[nral-1], + &aa[0], &aa[1], &aa[2], &aa[3], &aa[4], &aa[5]); + + if (nread < nral) + { + too_few(wi); + return; + } + else if (nread == nral) + { + ftype = ifunc_index(d, 1); + } + + /* Check for double atoms and atoms out of bounds */ + for (i = 0; i < nral; i++) + { + if (aa[i] < 1 || aa[i] > at->nr) + { + gmx_fatal(FARGS, "[ file %s, line %d ]:\n" + "Atom index (%d) in %s out of bounds (1-%d).\n" + "This probably means that you have inserted topology section \"%s\"\n" + "in a part belonging to a different molecule than you intended to.\n" + "In that case move the \"%s\" section to the right molecule.", + get_warning_file(wi), get_warning_line(wi), + aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d)); + } + + for (j = i+1; (j < nral); j++) + { + if (aa[i] == aa[j]) + { + sprintf(errbuf, "Duplicate atom index (%d) in %s", aa[i], dir2str(d)); + warning(wi, errbuf); + } + } + } + + /* default force parameters */ + for (j = 0; (j < MAXATOMLIST); j++) + { + param.a[j] = aa[j]-1; + } + for (j = 0; (j < MAXFORCEPARAM); j++) + { + param.c[j] = 0.0; + } + + /* Get the cmap type for this cmap angle */ + bFound = default_cmap_params(bondtype, at, atype, ¶m, FALSE, &cmap_type, &ncmap_params); + + /* We want exactly one parameter (the cmap type in state A (currently no state B) back */ + if (bFound && ncmap_params == 1) + { + /* Put the values in the appropriate arrays */ + param.c[0] = cmap_type; + add_param_to_list(&bond[ftype], ¶m); + } + else + { + /* This is essentially the same check as in default_cmap_params() done one more time */ + gmx_fatal(FARGS, "Unable to assign a cmap type to torsion %d %d %d %d and %d\n", + param.a[0]+1, param.a[1]+1, param.a[2]+1, param.a[3]+1, param.a[4]+1); + } +} + + + +void push_vsitesn(directive d, t_params bond[], + t_atoms *at, char *line, + warninp_t wi) +{ + char *ptr; + int type, ftype, j, n, ret, nj, a; + int *atc = NULL; + double *weight = NULL, weight_tot; + t_param param; + + /* default force parameters */ + for (j = 0; (j < MAXATOMLIST); j++) + { + param.a[j] = NOTSET; + } + for (j = 0; (j < MAXFORCEPARAM); j++) + { + param.c[j] = 0.0; + } + + ptr = line; + ret = sscanf(ptr, "%d%n", &a, &n); + ptr += n; + if (ret == 0) + { + gmx_fatal(FARGS, "[ file %s, line %d ]:\n" + " Expected an atom index in section \"%s\"", + get_warning_file(wi), get_warning_line(wi), + dir2str(d)); + } + + param.a[0] = a - 1; + + ret = sscanf(ptr, "%d%n", &type, &n); + ptr += n; + ftype = ifunc_index(d, type); + + weight_tot = 0; + nj = 0; + do + { + ret = sscanf(ptr, "%d%n", &a, &n); + ptr += n; + if (ret > 0) + { + if (nj % 20 == 0) + { + srenew(atc, nj+20); + srenew(weight, nj+20); + } + atc[nj] = a - 1; + switch (type) + { + case 1: + weight[nj] = 1; + break; + case 2: + /* Here we use the A-state mass as a parameter. + * Note that the B-state mass has no influence. + */ + weight[nj] = at->atom[atc[nj]].m; + break; + case 3: + weight[nj] = -1; + ret = sscanf(ptr, "%lf%n", &(weight[nj]), &n); + ptr += n; + if (weight[nj] < 0) + { + gmx_fatal(FARGS, "[ file %s, line %d ]:\n" + " No weight or negative weight found for vsiten constructing atom %d (atom index %d)", + get_warning_file(wi), get_warning_line(wi), + nj+1, atc[nj]+1); + } + break; + default: + gmx_fatal(FARGS, "Unknown vsiten type %d", type); + } + weight_tot += weight[nj]; + nj++; + } + } + while (ret > 0); + + if (nj == 0) + { + gmx_fatal(FARGS, "[ file %s, line %d ]:\n" + " Expected more than one atom index in section \"%s\"", + get_warning_file(wi), get_warning_line(wi), + dir2str(d)); + } + + if (weight_tot == 0) + { + gmx_fatal(FARGS, "[ file %s, line %d ]:\n" + " The total mass of the construting atoms is zero", + get_warning_file(wi), get_warning_line(wi)); + } + + for (j = 0; j < nj; j++) + { + param.a[1] = atc[j]; + param.c[0] = nj; + param.c[1] = weight[j]/weight_tot; + /* Put the values in the appropriate arrays */ + add_param_to_list (&bond[ftype], ¶m); + } + + sfree(atc); + sfree(weight); +} + +void push_mol(int nrmols, t_molinfo mols[], char *pline, int *whichmol, + int *nrcopies, + warninp_t wi) +{ + int i, copies; + char type[STRLEN]; + + *nrcopies = 0; + if (sscanf(pline, "%s%d", type, &copies) != 2) + { + too_few(wi); + return; + } + + /* search moleculename */ + for (i = 0; ((i < nrmols) && gmx_strcasecmp(type, *(mols[i].name))); i++) + { + ; + } + + if (i < nrmols) + { + *nrcopies = copies; + *whichmol = i; + } + else + { + gmx_fatal(FARGS, "No such moleculetype %s", type); + } +} + +void init_block2(t_block2 *b2, int natom) +{ + int i; + + b2->nr = natom; + snew(b2->nra, b2->nr); + snew(b2->a, b2->nr); + for (i = 0; (i < b2->nr); i++) + { + b2->a[i] = NULL; + } +} + +void done_block2(t_block2 *b2) +{ + int i; + + if (b2->nr) + { + for (i = 0; (i < b2->nr); i++) + { + sfree(b2->a[i]); + } + sfree(b2->a); + sfree(b2->nra); + b2->nr = 0; + } +} + +void push_excl(char *line, t_block2 *b2) +{ + int i, j; + int n; + char base[STRLEN], format[STRLEN]; + + if (sscanf(line, "%d", &i) == 0) + { + return; + } + + if ((1 <= i) && (i <= b2->nr)) + { + i--; + } + else + { + if (debug) + { + fprintf(debug, "Unbound atom %d\n", i-1); + } + return; + } + strcpy(base, "%*d"); + do + { + strcpy(format, base); + strcat(format, "%d"); + n = sscanf(line, format, &j); + if (n == 1) + { + if ((1 <= j) && (j <= b2->nr)) + { + j--; + srenew(b2->a[i], ++(b2->nra[i])); + b2->a[i][b2->nra[i]-1] = j; + /* also add the reverse exclusion! */ + srenew(b2->a[j], ++(b2->nra[j])); + b2->a[j][b2->nra[j]-1] = i; + strcat(base, "%*d"); + } + else + { + gmx_fatal(FARGS, "Invalid Atomnr j: %d, b2->nr: %d\n", j, b2->nr); + } + } + } + while (n == 1); +} + +void b_to_b2(t_blocka *b, t_block2 *b2) +{ + int i; + atom_id j, a; + + for (i = 0; (i < b->nr); i++) + { + for (j = b->index[i]; (j < b->index[i+1]); j++) + { + a = b->a[j]; + srenew(b2->a[i], ++b2->nra[i]); + b2->a[i][b2->nra[i]-1] = a; + } + } +} + +void b2_to_b(t_block2 *b2, t_blocka *b) +{ + int i, nra; + atom_id j; + + nra = 0; + for (i = 0; (i < b2->nr); i++) + { + b->index[i] = nra; + for (j = 0; (j < b2->nra[i]); j++) + { + b->a[nra+j] = b2->a[i][j]; + } + nra += b2->nra[i]; + } + /* terminate list */ + b->index[i] = nra; +} + +static int icomp(const void *v1, const void *v2) +{ + return (*((atom_id *) v1))-(*((atom_id *) v2)); +} + +void merge_excl(t_blocka *excl, t_block2 *b2) +{ + int i, k; + atom_id j; + int nra; + + if (!b2->nr) + { + return; + } + else if (b2->nr != excl->nr) + { + gmx_fatal(FARGS, "DEATH HORROR: b2->nr = %d, while excl->nr = %d", + b2->nr, excl->nr); + } + else if (debug) + { + fprintf(debug, "Entering merge_excl\n"); + } + + /* First copy all entries from excl to b2 */ + b_to_b2(excl, b2); + + /* Count and sort the exclusions */ + nra = 0; + for (i = 0; (i < b2->nr); i++) + { + if (b2->nra[i] > 0) + { + /* remove double entries */ + qsort(b2->a[i], (size_t)b2->nra[i], (size_t)sizeof(b2->a[i][0]), icomp); + k = 1; + for (j = 1; (j < b2->nra[i]); j++) + { + if (b2->a[i][j] != b2->a[i][k-1]) + { + b2->a[i][k] = b2->a[i][j]; + k++; + } + } + b2->nra[i] = k; + nra += b2->nra[i]; + } + } + excl->nra = nra; + srenew(excl->a, excl->nra); + + b2_to_b(b2, excl); +} + +int add_atomtype_decoupled(t_symtab *symtab, gpp_atomtype_t at, + t_nbparam ***nbparam, t_nbparam ***pair) +{ + t_atom atom; + t_param param; + int i, nr; + + /* Define an atom type with all parameters set to zero (no interactions) */ + atom.q = 0.0; + atom.m = 0.0; + /* Type for decoupled atoms could be anything, + * this should be changed automatically later when required. + */ + atom.ptype = eptAtom; + for (i = 0; (i < MAXFORCEPARAM); i++) + { + param.c[i] = 0.0; + } + + nr = add_atomtype(at, symtab, &atom, "decoupled", ¶m, -1, 0.0, 0.0, 0.0, 0, 0, 0); + + /* Add space in the non-bonded parameters matrix */ + realloc_nb_params(at, nbparam, pair); + + return nr; +} + +static void convert_pairs_to_pairsQ(t_params *plist, + real fudgeQQ, t_atoms *atoms) +{ + t_param *paramp1, *paramp2, *paramnew; + int i, j, p1nr, p2nr, p2newnr; + + /* Add the pair list to the pairQ list */ + p1nr = plist[F_LJ14].nr; + p2nr = plist[F_LJC14_Q].nr; + p2newnr = p1nr + p2nr; + snew(paramnew, p2newnr); + + paramp1 = plist[F_LJ14].param; + paramp2 = plist[F_LJC14_Q].param; + + /* Fill in the new F_LJC14_Q array with the old one. NOTE: + it may be possible to just ADD the converted F_LJ14 array + to the old F_LJC14_Q array, but since we have to create + a new sized memory structure, better just to deep copy it all. + */ + + for (i = 0; i < p2nr; i++) + { + /* Copy over parameters */ + for (j = 0; j < 5; j++) /* entries are 0-4 for F_LJC14_Q */ + { + paramnew[i].c[j] = paramp2[i].c[j]; + } + + /* copy over atoms */ + for (j = 0; j < 2; j++) + { + paramnew[i].a[j] = paramp2[i].a[j]; + } + } + + for (i = p2nr; i < p2newnr; i++) + { + j = i-p2nr; + + /* Copy over parameters */ + paramnew[i].c[0] = fudgeQQ; + paramnew[i].c[1] = atoms->atom[paramp1[j].a[0]].q; + paramnew[i].c[2] = atoms->atom[paramp1[j].a[1]].q; + paramnew[i].c[3] = paramp1[j].c[0]; + paramnew[i].c[4] = paramp1[j].c[1]; + + /* copy over atoms */ + paramnew[i].a[0] = paramp1[j].a[0]; + paramnew[i].a[1] = paramp1[j].a[1]; + } + + /* free the old pairlists */ + sfree(plist[F_LJC14_Q].param); + sfree(plist[F_LJ14].param); + + /* now assign the new data to the F_LJC14_Q structure */ + plist[F_LJC14_Q].param = paramnew; + plist[F_LJC14_Q].nr = p2newnr; + + /* Empty the LJ14 pairlist */ + plist[F_LJ14].nr = 0; + plist[F_LJ14].param = NULL; +} + +static void generate_LJCpairsNB(t_molinfo *mol, int nb_funct, t_params *nbp) +{ + int n, ntype, i, j, k; + t_atom *atom; + t_blocka *excl; + gmx_bool bExcl; + t_param param; + + n = mol->atoms.nr; + atom = mol->atoms.atom; + + ntype = sqrt(nbp->nr); + + for (i = 0; i < MAXATOMLIST; i++) + { + param.a[i] = NOTSET; + } + for (i = 0; i < MAXFORCEPARAM; i++) + { + param.c[i] = NOTSET; + } + + /* Add a pair interaction for all non-excluded atom pairs */ + excl = &mol->excls; + for (i = 0; i < n; i++) + { + for (j = i+1; j < n; j++) + { + bExcl = FALSE; + for (k = excl->index[i]; k < excl->index[i+1]; k++) + { + if (excl->a[k] == j) + { + bExcl = TRUE; + } + } + if (!bExcl) + { + if (nb_funct != F_LJ) + { + gmx_fatal(FARGS, "Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones"); + } + param.a[0] = i; + param.a[1] = j; + param.c[0] = atom[i].q; + param.c[1] = atom[j].q; + param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0]; + param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1]; + add_param_to_list(&mol->plist[F_LJC_PAIRS_NB], ¶m); + } + } + } +} + +static void set_excl_all(t_blocka *excl) +{ + int nat, i, j, k; + + /* Get rid of the current exclusions and exclude all atom pairs */ + nat = excl->nr; + excl->nra = nat*nat; + srenew(excl->a, excl->nra); + k = 0; + for (i = 0; i < nat; i++) + { + excl->index[i] = k; + for (j = 0; j < nat; j++) + { + excl->a[k++] = j; + } + } + excl->index[nat] = k; +} + +static void decouple_atoms(t_atoms *atoms, int atomtype_decouple, + int couple_lam0, int couple_lam1) +{ + int i; + + for (i = 0; i < atoms->nr; i++) + { + if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamVDW) + { + atoms->atom[i].q = 0.0; + } + if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamQ) + { + atoms->atom[i].type = atomtype_decouple; + } + if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamVDW) + { + atoms->atom[i].qB = 0.0; + } + if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamQ) + { + atoms->atom[i].typeB = atomtype_decouple; + } + } +} + +void convert_moltype_couple(t_molinfo *mol, int atomtype_decouple, real fudgeQQ, + int couple_lam0, int couple_lam1, + gmx_bool bCoupleIntra, int nb_funct, t_params *nbp) +{ + convert_pairs_to_pairsQ(mol->plist, fudgeQQ, &mol->atoms); + if (!bCoupleIntra) + { + generate_LJCpairsNB(mol, nb_funct, nbp); + set_excl_all(&mol->excls); + } + decouple_atoms(&mol->atoms, atomtype_decouple, couple_lam0, couple_lam1); +}