1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
46 #include "gmx_fatal.h"
51 #include "fflibutil.h"
53 gpp_atomtype_t read_atype(const char *ffdir, t_symtab *tab)
58 char buf[STRLEN], name[STRLEN];
65 nfile = fflib_search_file_end(ffdir, ".atp", TRUE, &file);
70 for (f = 0; f < nfile; f++)
72 in = fflib_open(file[f]);
75 /* Skip blank or comment-only lines */
78 fgets2(buf, STRLEN, in);
85 while (!feof(in) && NULL != buf && strlen(buf) == 0);
87 if ((buf != NULL) && (sscanf(buf, "%s%lf", name, &m) == 2))
90 add_atomtype(at, tab, a, name, nb, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 );
91 fprintf(stderr, "\rAtomtype %d", nratt+1);
97 fprintf(stderr, "\n");
103 static void print_resatoms(FILE *out, gpp_atomtype_t atype, t_restp *rtp)
108 /* fprintf(out,"%5s\n",rtp->resname);
109 fprintf(out,"%5d\n",rtp->natom); */
110 fprintf(out, "[ %s ]\n", rtp->resname);
111 fprintf(out, " [ atoms ]\n");
113 for (j = 0; (j < rtp->natom); j++)
115 tp = rtp->atom[j].type;
116 tpnm = get_atomtype_name(tp, atype);
119 gmx_fatal(FARGS, "Incorrect atomtype (%d)", tp);
121 fprintf(out, "%6s %6s %8.3f %6d\n",
122 *(rtp->atomname[j]), tpnm, rtp->atom[j].q, rtp->cgnr[j]);
126 static gmx_bool read_atoms(FILE *in, char *line,
127 t_restp *r0, t_symtab *tab, gpp_atomtype_t atype)
129 int i, j, cg, maxentries;
130 char buf[256], buf1[256];
139 while (get_a_line(in, line, STRLEN) && (strchr(line, '[') == NULL))
141 if (sscanf(line, "%s%s%lf%d", buf, buf1, &q, &cg) != 4)
148 srenew(r0->atom, maxentries);
149 srenew(r0->atomname, maxentries);
150 srenew(r0->cgnr, maxentries);
152 r0->atomname[i] = put_symtab(tab, buf);
155 j = get_atomtype_type(buf1, atype);
158 gmx_fatal(FARGS, "Atom type %s (residue %s) not found in atomtype "
159 "database", buf1, r0->resname);
161 r0->atom[i].type = j;
162 r0->atom[i].m = get_atomtype_massA(j, atype);
167 srenew(r0->atomname, i);
173 gmx_bool read_bondeds(int bt, FILE *in, char *line, t_restp *rtp)
178 maxrb = rtp->rb[bt].nb;
179 while (get_a_line(in, line, STRLEN) && (strchr(line, '[') == NULL))
181 if (rtp->rb[bt].nb >= maxrb)
184 srenew(rtp->rb[bt].b, maxrb);
187 for (j = 0; j < btsNiatoms[bt]; j++)
189 if (sscanf(line+n, "%s%n", str, &ni) == 1)
191 rtp->rb[bt].b[rtp->rb[bt].nb].a[j] = strdup(str);
199 for (; j < MAXATOMLIST; j++)
201 rtp->rb[bt].b[rtp->rb[bt].nb].a[j] = NULL;
203 while (isspace(line[n]))
208 rtp->rb[bt].b[rtp->rb[bt].nb].s = strdup(line+n);
211 /* give back unused memory */
212 srenew(rtp->rb[bt].b, rtp->rb[bt].nb);
217 static void print_resbondeds(FILE *out, int bt, t_restp *rtp)
223 fprintf(out, " [ %s ]\n", btsNames[bt]);
225 for (i = 0; i < rtp->rb[bt].nb; i++)
227 for (j = 0; j < btsNiatoms[bt]; j++)
229 fprintf(out, "%6s ", rtp->rb[bt].b[i].a[j]);
231 if (rtp->rb[bt].b[i].s[0])
233 fprintf(out, " %s", rtp->rb[bt].b[i].s);
240 static void check_rtp(int nrtp, t_restp rtp[], char *libfn)
244 /* check for double entries, assuming list is already sorted */
245 for (i = 1; (i < nrtp); i++)
247 if (gmx_strcasecmp(rtp[i-1].resname, rtp[i].resname) == 0)
249 fprintf(stderr, "WARNING double entry %s in file %s\n",
250 rtp[i].resname, libfn);
255 static int comprtp(const void *a, const void *b)
262 return gmx_strcasecmp(ra->resname, rb->resname);
265 int get_bt(char* header)
269 for (i = 0; i < ebtsNR; i++)
271 if (gmx_strcasecmp(btsNames[i], header) == 0)
279 void clear_t_restp(t_restp *rrtp)
281 memset((void *)rrtp, 0, sizeof(t_restp));
284 /* print all the ebtsNR type numbers */
285 void print_resall_header(FILE *out, t_restp rtp[])
287 fprintf(out, "[ bondedtypes ]\n");
288 fprintf(out, "; bonds angles dihedrals impropers all_dihedrals nr_exclusions HH14 remove_dih\n");
289 fprintf(out, " %5d %6d %9d %9d %14d %14d %14d %14d\n\n",
294 rtp[0].bKeepAllGeneratedDihedrals,
296 rtp[0].bGenerateHH14Interactions,
297 rtp[0].bRemoveDihedralIfWithImproper);
300 void print_resall(FILE *out, int nrtp, t_restp rtp[],
301 gpp_atomtype_t atype)
310 print_resall_header(out, rtp);
312 for (i = 0; i < nrtp; i++)
314 if (rtp[i].natom > 0)
316 print_resatoms(out, atype, &rtp[i]);
317 for (bt = 0; bt < ebtsNR; bt++)
319 print_resbondeds(out, bt, &rtp[i]);
325 void read_resall(char *rrdb, int *nrtpptr, t_restp **rtp,
326 gpp_atomtype_t atype, t_symtab *tab,
327 gmx_bool bAllowOverrideRTP)
330 char filebase[STRLEN], *ptr, line[STRLEN], header[STRLEN];
331 int i, nrtp, maxrtp, bt, nparam;
332 int dum1, dum2, dum3;
333 t_restp *rrtp, *header_settings;
334 gmx_bool bNextResidue, bError;
337 fflib_filename_base(rrdb, filebase, STRLEN);
339 in = fflib_open(rrdb);
343 fprintf(debug, "%9s %5s", "Residue", "atoms");
344 for (i = 0; i < ebtsNR; i++)
346 fprintf(debug, " %10s", btsNames[i]);
348 fprintf(debug, "\n");
350 snew(header_settings, 1);
352 /* these bonded parameters will overwritten be when *
353 * there is a [ bondedtypes ] entry in the .rtp file */
354 header_settings->rb[ebtsBONDS].type = 1; /* normal bonds */
355 header_settings->rb[ebtsANGLES].type = 1; /* normal angles */
356 header_settings->rb[ebtsPDIHS].type = 1; /* normal dihedrals */
357 header_settings->rb[ebtsIDIHS].type = 2; /* normal impropers */
358 header_settings->rb[ebtsEXCLS].type = 1; /* normal exclusions */
359 header_settings->rb[ebtsCMAP].type = 1; /* normal cmap torsions */
361 header_settings->bKeepAllGeneratedDihedrals = FALSE;
362 header_settings->nrexcl = 3;
363 header_settings->bGenerateHH14Interactions = TRUE;
364 header_settings->bRemoveDihedralIfWithImproper = TRUE;
366 /* Column 5 & 6 aren't really bonded types, but we include
367 * them here to avoid introducing a new section:
368 * Column 5 : This controls the generation of dihedrals from the bonding.
369 * All possible dihedrals are generated automatically. A value of
370 * 1 here means that all these are retained. A value of
371 * 0 here requires generated dihedrals be removed if
372 * * there are any dihedrals on the same central atoms
373 * specified in the residue topology, or
374 * * there are other identical generated dihedrals
375 * sharing the same central atoms, or
376 * * there are other generated dihedrals sharing the
377 * same central bond that have fewer hydrogen atoms
378 * Column 6: Number of bonded neighbors to exclude.
379 * Column 7: Generate 1,4 interactions between two hydrogen atoms
380 * Column 8: Remove proper dihedrals if centered on the same bond
381 * as an improper dihedral
383 get_a_line(in, line, STRLEN);
384 if (!get_header(line, header))
386 gmx_fatal(FARGS, "in .rtp file at line:\n%s\n", line);
388 if (gmx_strncasecmp("bondedtypes", header, 5) == 0)
390 get_a_line(in, line, STRLEN);
391 if ((nparam = sscanf(line, "%d %d %d %d %d %d %d %d",
392 &header_settings->rb[ebtsBONDS].type, &header_settings->rb[ebtsANGLES].type,
393 &header_settings->rb[ebtsPDIHS].type, &header_settings->rb[ebtsIDIHS].type,
394 &dum1, &header_settings->nrexcl, &dum2, &dum3)) < 4)
396 gmx_fatal(FARGS, "need 4 to 8 parameters in the header of .rtp file %s at line:\n%s\n", rrdb, line);
398 header_settings->bKeepAllGeneratedDihedrals = (dum1 != 0);
399 header_settings->bGenerateHH14Interactions = (dum2 != 0);
400 header_settings->bRemoveDihedralIfWithImproper = (dum3 != 0);
401 get_a_line(in, line, STRLEN);
404 fprintf(stderr, "Using default: not generating all possible dihedrals\n");
405 header_settings->bKeepAllGeneratedDihedrals = FALSE;
409 fprintf(stderr, "Using default: excluding 3 bonded neighbors\n");
410 header_settings->nrexcl = 3;
414 fprintf(stderr, "Using default: generating 1,4 H--H interactions\n");
415 header_settings->bGenerateHH14Interactions = TRUE;
419 fprintf(stderr, "Using default: removing proper dihedrals found on the same bond as a proper dihedral\n");
420 header_settings->bRemoveDihedralIfWithImproper = TRUE;
426 "Reading .rtp file without '[ bondedtypes ]' directive,\n"
427 "Will proceed as if the entry was:\n");
428 print_resall_header(stderr, header_settings);
430 /* We don't know the current size of rrtp, but simply realloc immediately */
439 srenew(rrtp, maxrtp);
442 /* Initialise rtp entry structure */
443 rrtp[nrtp] = *header_settings;
444 if (!get_header(line, header))
446 gmx_fatal(FARGS, "in .rtp file at line:\n%s\n", line);
448 rrtp[nrtp].resname = strdup(header);
449 rrtp[nrtp].filebase = strdup(filebase);
451 get_a_line(in, line, STRLEN);
453 bNextResidue = FALSE;
456 if (!get_header(line, header))
465 /* header is an bonded directive */
466 bError = !read_bondeds(bt, in, line, &rrtp[nrtp]);
468 else if (gmx_strncasecmp("atoms", header, 5) == 0)
470 /* header is the atoms directive */
471 bError = !read_atoms(in, line, &(rrtp[nrtp]), tab, atype);
475 /* else header must be a residue name */
481 gmx_fatal(FARGS, "in .rtp file in residue %s at line:\n%s\n",
482 rrtp[nrtp].resname, line);
485 while (!feof(in) && !bNextResidue);
487 if (rrtp[nrtp].natom == 0)
489 gmx_fatal(FARGS, "No atoms found in .rtp file in residue %s\n",
494 fprintf(debug, "%3d %5s %5d",
495 nrtp+1, rrtp[nrtp].resname, rrtp[nrtp].natom);
496 for (i = 0; i < ebtsNR; i++)
498 fprintf(debug, " %10d", rrtp[nrtp].rb[i].nb);
500 fprintf(debug, "\n");
504 for (i = 0; i < nrtp; i++)
506 if (gmx_strcasecmp(rrtp[i].resname, rrtp[nrtp].resname) == 0)
515 fprintf(stderr, "\rResidue %d", nrtp);
519 if (firstrtp >= *nrtpptr)
521 gmx_fatal(FARGS, "Found a second entry for '%s' in '%s'",
522 rrtp[nrtp].resname, rrdb);
524 if (bAllowOverrideRTP)
526 fprintf(stderr, "\n");
527 fprintf(stderr, "Found another rtp entry for '%s' in '%s', ignoring this entry and keeping the one from '%s.rtp'\n",
528 rrtp[nrtp].resname, rrdb, rrtp[firstrtp].filebase);
529 /* We should free all the data for this entry.
530 * The current code gives a lot of dangling pointers.
532 clear_t_restp(&rrtp[nrtp]);
536 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);
541 /* give back unused memory */
544 fprintf(stderr, "\nSorting it all out...\n");
545 qsort(rrtp, nrtp, (size_t)sizeof(rrtp[0]), comprtp);
547 check_rtp(nrtp, rrtp, rrdb);
553 /************************************************************
557 ***********************************************************/
558 static gmx_bool is_sign(char c)
560 return (c == '+' || c == '-');
563 /* Compares if the strins match except for a sign at the end
564 * of (only) one of the two.
566 static int neq_str_sign(const char *a1, const char *a2)
570 l1 = (int)strlen(a1);
571 l2 = (int)strlen(a2);
575 ((l1 == l2+1 && is_sign(a1[l1-1])) ||
576 (l2 == l1+1 && is_sign(a2[l2-1]))) &&
577 gmx_strncasecmp(a1, a2, lm) == 0)
587 char *search_rtp(const char *key, int nrtp, t_restp rtp[])
589 int i, n, nbest, best, besti;
590 char bestbuf[STRLEN];
594 /* We want to match at least one character */
596 for (i = 0; (i < nrtp); i++)
598 if (gmx_strcasecmp(key, rtp[i].resname) == 0)
606 /* Allow a mismatch of at most a sign character (with warning) */
607 n = neq_str_sign(key, rtp[i].resname);
609 n+1 >= (int)strlen(key) &&
610 n+1 >= (int)strlen(rtp[i].resname))
616 strcpy(bestbuf, rtp[besti].resname);
620 strcat(bestbuf, " or ");
621 strcat(bestbuf, rtp[i].resname);
636 gmx_fatal(FARGS, "Residue '%s' not found in residue topology database, looks a bit like %s", key, bestbuf);
638 else if (besti == -1)
640 gmx_fatal(FARGS, "Residue '%s' not found in residue topology database", key);
642 if (gmx_strcasecmp(rtp[besti].resname, key) != 0)
645 "\nWARNING: '%s' not found in residue topology database, "
646 "trying to use '%s'\n\n", key, rtp[besti].resname);
649 return rtp[besti].resname;
652 t_restp *get_restp(const char *rtpname, int nrtp, t_restp rtp[])
657 while (i < nrtp && gmx_strcasecmp(rtpname, rtp[i].resname) != 0)
663 /* This should never happen, since search_rtp should have been called
664 * before calling get_restp.
666 gmx_fatal(FARGS, "Residue type '%s' not found in residue topology database", rtpname);