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, 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.
47 #include "gromacs/gmxpreprocess/fflibutil.h"
48 #include "gromacs/gmxpreprocess/notset.h"
49 #include "gromacs/gmxpreprocess/pgutil.h"
50 #include "gromacs/topology/symtab.h"
51 #include "gromacs/utility/cstringutil.h"
52 #include "gromacs/utility/fatalerror.h"
53 #include "gromacs/utility/futil.h"
54 #include "gromacs/utility/smalloc.h"
55 #include "gromacs/utility/strdb.h"
57 gpp_atomtype_t read_atype(const char *ffdir, t_symtab *tab)
62 char buf[STRLEN], name[STRLEN];
69 nfile = fflib_search_file_end(ffdir, ".atp", TRUE, &file);
74 for (f = 0; f < nfile; f++)
76 in = fflib_open(file[f]);
79 /* Skip blank or comment-only lines */
82 if (fgets2(buf, STRLEN, in) != NULL)
88 while (!feof(in) && strlen(buf) == 0);
90 if (sscanf(buf, "%s%lf", name, &m) == 2)
93 add_atomtype(at, tab, a, name, nb, 0, 0.0, 0.0, 0.0, 0, 0.0, 0.0 );
94 fprintf(stderr, "\rAtomtype %d", ++nratt);
99 fprintf(stderr, "\nInvalid format: %s\n", buf);
105 fprintf(stderr, "\n");
111 static void print_resatoms(FILE *out, gpp_atomtype_t atype, t_restp *rtp)
116 /* fprintf(out,"%5s\n",rtp->resname);
117 fprintf(out,"%5d\n",rtp->natom); */
118 fprintf(out, "[ %s ]\n", rtp->resname);
119 fprintf(out, " [ atoms ]\n");
121 for (j = 0; (j < rtp->natom); j++)
123 tp = rtp->atom[j].type;
124 tpnm = get_atomtype_name(tp, atype);
127 gmx_fatal(FARGS, "Incorrect atomtype (%d)", tp);
129 fprintf(out, "%6s %6s %8.3f %6d\n",
130 *(rtp->atomname[j]), tpnm, rtp->atom[j].q, rtp->cgnr[j]);
134 static gmx_bool read_atoms(FILE *in, char *line,
135 t_restp *r0, t_symtab *tab, gpp_atomtype_t atype)
137 int i, j, cg, maxentries;
138 char buf[256], buf1[256];
147 while (get_a_line(in, line, STRLEN) && (strchr(line, '[') == NULL))
149 if (sscanf(line, "%s%s%lf%d", buf, buf1, &q, &cg) != 4)
156 srenew(r0->atom, maxentries);
157 srenew(r0->atomname, maxentries);
158 srenew(r0->cgnr, maxentries);
160 r0->atomname[i] = put_symtab(tab, buf);
163 j = get_atomtype_type(buf1, atype);
166 gmx_fatal(FARGS, "Atom type %s (residue %s) not found in atomtype "
167 "database", buf1, r0->resname);
169 r0->atom[i].type = j;
170 r0->atom[i].m = get_atomtype_massA(j, atype);
175 srenew(r0->atomname, i);
181 gmx_bool read_bondeds(int bt, FILE *in, char *line, t_restp *rtp)
186 maxrb = rtp->rb[bt].nb;
187 while (get_a_line(in, line, STRLEN) && (strchr(line, '[') == NULL))
189 if (rtp->rb[bt].nb >= maxrb)
192 srenew(rtp->rb[bt].b, maxrb);
195 for (j = 0; j < btsNiatoms[bt]; j++)
197 if (sscanf(line+n, "%s%n", str, &ni) == 1)
199 rtp->rb[bt].b[rtp->rb[bt].nb].a[j] = gmx_strdup(str);
207 for (; j < MAXATOMLIST; j++)
209 rtp->rb[bt].b[rtp->rb[bt].nb].a[j] = NULL;
211 while (isspace(line[n]))
216 rtp->rb[bt].b[rtp->rb[bt].nb].s = gmx_strdup(line+n);
219 /* give back unused memory */
220 srenew(rtp->rb[bt].b, rtp->rb[bt].nb);
225 static void print_resbondeds(FILE *out, int bt, t_restp *rtp)
231 fprintf(out, " [ %s ]\n", btsNames[bt]);
233 for (i = 0; i < rtp->rb[bt].nb; i++)
235 for (j = 0; j < btsNiatoms[bt]; j++)
237 fprintf(out, "%6s ", rtp->rb[bt].b[i].a[j]);
239 if (rtp->rb[bt].b[i].s[0])
241 fprintf(out, " %s", rtp->rb[bt].b[i].s);
248 static void check_rtp(int nrtp, t_restp rtp[], char *libfn)
252 /* check for double entries, assuming list is already sorted */
253 for (i = 1; (i < nrtp); i++)
255 if (gmx_strcasecmp(rtp[i-1].resname, rtp[i].resname) == 0)
257 fprintf(stderr, "WARNING double entry %s in file %s\n",
258 rtp[i].resname, libfn);
263 int get_bt(char* header)
267 for (i = 0; i < ebtsNR; i++)
269 if (gmx_strcasecmp(btsNames[i], header) == 0)
277 void clear_t_restp(t_restp *rrtp)
279 memset((void *)rrtp, 0, sizeof(t_restp));
282 /* print all the ebtsNR type numbers */
283 void print_resall_header(FILE *out, t_restp rtp[])
285 fprintf(out, "[ bondedtypes ]\n");
286 fprintf(out, "; bonds angles dihedrals impropers all_dihedrals nr_exclusions HH14 remove_dih\n");
287 fprintf(out, " %5d %6d %9d %9d %14d %14d %14d %14d\n\n",
292 rtp[0].bKeepAllGeneratedDihedrals,
294 rtp[0].bGenerateHH14Interactions,
295 rtp[0].bRemoveDihedralIfWithImproper);
298 void print_resall(FILE *out, int nrtp, t_restp rtp[],
299 gpp_atomtype_t atype)
308 print_resall_header(out, rtp);
310 for (i = 0; i < nrtp; i++)
312 if (rtp[i].natom > 0)
314 print_resatoms(out, atype, &rtp[i]);
315 for (bt = 0; bt < ebtsNR; bt++)
317 print_resbondeds(out, bt, &rtp[i]);
323 void read_resall(char *rrdb, int *nrtpptr, t_restp **rtp,
324 gpp_atomtype_t atype, t_symtab *tab,
325 gmx_bool bAllowOverrideRTP)
328 char filebase[STRLEN], line[STRLEN], header[STRLEN];
329 int i, nrtp, maxrtp, bt, nparam;
330 int dum1, dum2, dum3;
331 t_restp *rrtp, *header_settings;
332 gmx_bool bNextResidue, bError;
335 fflib_filename_base(rrdb, filebase, STRLEN);
337 in = fflib_open(rrdb);
341 fprintf(debug, "%9s %5s", "Residue", "atoms");
342 for (i = 0; i < ebtsNR; i++)
344 fprintf(debug, " %10s", btsNames[i]);
346 fprintf(debug, "\n");
348 snew(header_settings, 1);
350 /* these bonded parameters will overwritten be when *
351 * there is a [ bondedtypes ] entry in the .rtp file */
352 header_settings->rb[ebtsBONDS].type = 1; /* normal bonds */
353 header_settings->rb[ebtsANGLES].type = 1; /* normal angles */
354 header_settings->rb[ebtsPDIHS].type = 1; /* normal dihedrals */
355 header_settings->rb[ebtsIDIHS].type = 2; /* normal impropers */
356 header_settings->rb[ebtsEXCLS].type = 1; /* normal exclusions */
357 header_settings->rb[ebtsCMAP].type = 1; /* normal cmap torsions */
359 header_settings->bKeepAllGeneratedDihedrals = FALSE;
360 header_settings->nrexcl = 3;
361 header_settings->bGenerateHH14Interactions = TRUE;
362 header_settings->bRemoveDihedralIfWithImproper = TRUE;
364 /* Column 5 & 6 aren't really bonded types, but we include
365 * them here to avoid introducing a new section:
366 * Column 5 : This controls the generation of dihedrals from the bonding.
367 * All possible dihedrals are generated automatically. A value of
368 * 1 here means that all these are retained. A value of
369 * 0 here requires generated dihedrals be removed if
370 * * there are any dihedrals on the same central atoms
371 * specified in the residue topology, or
372 * * there are other identical generated dihedrals
373 * sharing the same central atoms, or
374 * * there are other generated dihedrals sharing the
375 * same central bond that have fewer hydrogen atoms
376 * Column 6: Number of bonded neighbors to exclude.
377 * Column 7: Generate 1,4 interactions between two hydrogen atoms
378 * Column 8: Remove proper dihedrals if centered on the same bond
379 * as an improper dihedral
381 get_a_line(in, line, STRLEN);
382 if (!get_header(line, header))
384 gmx_fatal(FARGS, "in .rtp file at line:\n%s\n", line);
386 if (gmx_strncasecmp("bondedtypes", header, 5) == 0)
388 get_a_line(in, line, STRLEN);
389 if ((nparam = sscanf(line, "%d %d %d %d %d %d %d %d",
390 &header_settings->rb[ebtsBONDS].type, &header_settings->rb[ebtsANGLES].type,
391 &header_settings->rb[ebtsPDIHS].type, &header_settings->rb[ebtsIDIHS].type,
392 &dum1, &header_settings->nrexcl, &dum2, &dum3)) < 4)
394 gmx_fatal(FARGS, "need 4 to 8 parameters in the header of .rtp file %s at line:\n%s\n", rrdb, line);
396 header_settings->bKeepAllGeneratedDihedrals = (dum1 != 0);
397 header_settings->bGenerateHH14Interactions = (dum2 != 0);
398 header_settings->bRemoveDihedralIfWithImproper = (dum3 != 0);
399 get_a_line(in, line, STRLEN);
402 fprintf(stderr, "Using default: not generating all possible dihedrals\n");
403 header_settings->bKeepAllGeneratedDihedrals = FALSE;
407 fprintf(stderr, "Using default: excluding 3 bonded neighbors\n");
408 header_settings->nrexcl = 3;
412 fprintf(stderr, "Using default: generating 1,4 H--H interactions\n");
413 header_settings->bGenerateHH14Interactions = TRUE;
417 fprintf(stderr, "Using default: removing proper dihedrals found on the same bond as a proper dihedral\n");
418 header_settings->bRemoveDihedralIfWithImproper = TRUE;
424 "Reading .rtp file without '[ bondedtypes ]' directive,\n"
425 "Will proceed as if the entry was:\n");
426 print_resall_header(stderr, header_settings);
428 /* We don't know the current size of rrtp, but simply realloc immediately */
437 srenew(rrtp, maxrtp);
440 /* Initialise rtp entry structure */
441 rrtp[nrtp] = *header_settings;
442 if (!get_header(line, header))
444 gmx_fatal(FARGS, "in .rtp file at line:\n%s\n", line);
446 rrtp[nrtp].resname = gmx_strdup(header);
447 rrtp[nrtp].filebase = gmx_strdup(filebase);
449 get_a_line(in, line, STRLEN);
451 bNextResidue = FALSE;
454 if (!get_header(line, header))
463 /* header is an bonded directive */
464 bError = !read_bondeds(bt, in, line, &rrtp[nrtp]);
466 else if (gmx_strncasecmp("atoms", header, 5) == 0)
468 /* header is the atoms directive */
469 bError = !read_atoms(in, line, &(rrtp[nrtp]), tab, atype);
473 /* else header must be a residue name */
479 gmx_fatal(FARGS, "in .rtp file in residue %s at line:\n%s\n",
480 rrtp[nrtp].resname, line);
483 while (!feof(in) && !bNextResidue);
485 if (rrtp[nrtp].natom == 0)
487 gmx_fatal(FARGS, "No atoms found in .rtp file in residue %s\n",
492 fprintf(debug, "%3d %5s %5d",
493 nrtp+1, rrtp[nrtp].resname, rrtp[nrtp].natom);
494 for (i = 0; i < ebtsNR; i++)
496 fprintf(debug, " %10d", rrtp[nrtp].rb[i].nb);
498 fprintf(debug, "\n");
502 for (i = 0; i < nrtp; i++)
504 if (gmx_strcasecmp(rrtp[i].resname, rrtp[nrtp].resname) == 0)
513 fprintf(stderr, "\rResidue %d", nrtp);
518 if (firstrtp >= *nrtpptr)
520 gmx_fatal(FARGS, "Found a second entry for '%s' in '%s'",
521 rrtp[nrtp].resname, rrdb);
523 if (bAllowOverrideRTP)
525 fprintf(stderr, "\n");
526 fprintf(stderr, "Found another rtp entry for '%s' in '%s', ignoring this entry and keeping the one from '%s.rtp'\n",
527 rrtp[nrtp].resname, rrdb, rrtp[firstrtp].filebase);
528 /* We should free all the data for this entry.
529 * The current code gives a lot of dangling pointers.
531 clear_t_restp(&rrtp[nrtp]);
535 gmx_fatal(FARGS, "Found rtp entries for '%s' in both '%s' and '%s'. If you want the first definition to override the second one, set the -rtpo option of pdb2gmx.", rrtp[nrtp].resname, rrtp[firstrtp].filebase, rrdb);
540 /* give back unused memory */
543 fprintf(stderr, "\nSorting it all out...\n");
544 std::sort(rrtp, rrtp+nrtp, [](const t_restp &a, const t_restp &b) {return gmx_strcasecmp(a.resname, b.resname) < 0; });
546 check_rtp(nrtp, rrtp, rrdb);
552 /************************************************************
556 ***********************************************************/
557 static gmx_bool is_sign(char c)
559 return (c == '+' || c == '-');
562 /* Compares if the strins match except for a sign at the end
563 * of (only) one of the two.
565 static int neq_str_sign(const char *a1, const char *a2)
569 l1 = (int)strlen(a1);
570 l2 = (int)strlen(a2);
571 lm = std::min(l1, l2);
574 ((l1 == l2+1 && is_sign(a1[l1-1])) ||
575 (l2 == l1+1 && is_sign(a2[l2-1]))) &&
576 gmx_strncasecmp(a1, a2, lm) == 0)
586 char *search_rtp(const char *key, int nrtp, t_restp rtp[])
588 int i, n, nbest, best, besti;
589 char bestbuf[STRLEN];
593 /* We want to match at least one character */
595 for (i = 0; (i < nrtp); i++)
597 if (gmx_strcasecmp(key, rtp[i].resname) == 0)
605 /* Allow a mismatch of at most a sign character (with warning) */
606 n = neq_str_sign(key, rtp[i].resname);
608 n+1 >= (int)strlen(key) &&
609 n+1 >= (int)strlen(rtp[i].resname))
615 strcpy(bestbuf, rtp[besti].resname);
619 strcat(bestbuf, " or ");
620 strcat(bestbuf, rtp[i].resname);
635 gmx_fatal(FARGS, "Residue '%s' not found in residue topology database, looks a bit like %s", key, bestbuf);
637 else if (besti == -1)
639 gmx_fatal(FARGS, "Residue '%s' not found in residue topology database", key);
641 if (gmx_strcasecmp(rtp[besti].resname, key) != 0)
644 "\nWARNING: '%s' not found in residue topology database, "
645 "trying to use '%s'\n\n", key, rtp[besti].resname);
648 return rtp[besti].resname;
651 t_restp *get_restp(const char *rtpname, int nrtp, t_restp rtp[])
656 while (i < nrtp && gmx_strcasecmp(rtpname, rtp[i].resname) != 0)
662 /* This should never happen, since search_rtp should have been called
663 * before calling get_restp.
665 gmx_fatal(FARGS, "Residue type '%s' not found in residue topology database", rtpname);