src/gromacs/gmxpreprocess/resall.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  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.
  15  *
  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.
  20  *
  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.
  25  *
  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.
  33  *
  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.
  36  */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include "sysstuff.h"
  42 #include <ctype.h>
  43 #include "string2.h"
  44 #include "gromacs/fileio/futil.h"
  45 #include "smalloc.h"
  46 #include "gmx_fatal.h"
  47 #include "symtab.h"
  48 #include "macros.h"
  49 #include "resall.h"
  50 #include "pgutil.h"
  51 #include "fflibutil.h"
  52
  53 #include "gromacs/fileio/strdb.h"
  54
  55 gpp_atomtype_t read_atype(const char *ffdir, t_symtab *tab)
  56 {
  57     int            nfile, f;
  58     char         **file;
  59     FILE          *in;
  60     char           buf[STRLEN], name[STRLEN];
  61     double         m;
  62     int            nratt = 0;
  63     gpp_atomtype_t at;
  64     t_atom        *a;
  65     t_param       *nb;
  66
  67     nfile = fflib_search_file_end(ffdir, ".atp", TRUE, &file);
  68     at    = init_atomtype();
  69     snew(a, 1);
  70     snew(nb, 1);
  71
  72     for (f = 0; f < nfile; f++)
  73     {
  74         in = fflib_open(file[f]);
  75         while (!feof(in))
  76         {
  77             /* Skip blank or comment-only lines */
  78             do
  79             {
  80                 fgets2(buf, STRLEN, in);
  81                 if (NULL != buf)
  82                 {
  83                     strip_comment(buf);
  84                     trim(buf);
  85                 }
  86             }
  87             while (!feof(in) && NULL != buf && strlen(buf) == 0);
  88
  89             if ((buf != NULL) && (sscanf(buf, "%s%lf", name, &m) == 2))
  90             {
  91                 a->m = m;
  92                 add_atomtype(at, tab, a, name, nb, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 );
  93                 fprintf(stderr, "\rAtomtype %d", nratt+1);
  94             }
  95         }
  96         gmx_ffclose(in);
  97         sfree(file[f]);
  98     }
  99     fprintf(stderr, "\n");
 100     sfree(file);
 101
 102     return at;
 103 }
 104
 105 static void print_resatoms(FILE *out, gpp_atomtype_t atype, t_restp *rtp)
 106 {
 107     int   j, tp;
 108     char *tpnm;
 109
 110     /* fprintf(out,"%5s\n",rtp->resname);
 111        fprintf(out,"%5d\n",rtp->natom); */
 112     fprintf(out, "[ %s ]\n", rtp->resname);
 113     fprintf(out, " [ atoms ]\n");
 114
 115     for (j = 0; (j < rtp->natom); j++)
 116     {
 117         tp   = rtp->atom[j].type;
 118         tpnm = get_atomtype_name(tp, atype);
 119         if (tpnm == NULL)
 120         {
 121             gmx_fatal(FARGS, "Incorrect atomtype (%d)", tp);
 122         }
 123         fprintf(out, "%6s  %6s  %8.3f  %6d\n",
 124                 *(rtp->atomname[j]), tpnm, rtp->atom[j].q, rtp->cgnr[j]);
 125     }
 126 }
 127
 128 static gmx_bool read_atoms(FILE *in, char *line,
 129                            t_restp *r0, t_symtab *tab, gpp_atomtype_t atype)
 130 {
 131     int    i, j, cg, maxentries;
 132     char   buf[256], buf1[256];
 133     double q;
 134
 135     /* Read Atoms */
 136     maxentries   = 0;
 137     r0->atom     =     NULL;
 138     r0->atomname = NULL;
 139     r0->cgnr     =     NULL;
 140     i            = 0;
 141     while (get_a_line(in, line, STRLEN) && (strchr(line, '[') == NULL))
 142     {
 143         if (sscanf(line, "%s%s%lf%d", buf, buf1, &q, &cg) != 4)
 144         {
 145             return FALSE;
 146         }
 147         if (i >= maxentries)
 148         {
 149             maxentries += 100;
 150             srenew(r0->atom,     maxentries);
 151             srenew(r0->atomname, maxentries);
 152             srenew(r0->cgnr,     maxentries);
 153         }
 154         r0->atomname[i] = put_symtab(tab, buf);
 155         r0->atom[i].q   = q;
 156         r0->cgnr[i]     = cg;
 157         j               = get_atomtype_type(buf1, atype);
 158         if (j == NOTSET)
 159         {
 160             gmx_fatal(FARGS, "Atom type %s (residue %s) not found in atomtype "
 161                       "database", buf1, r0->resname);
 162         }
 163         r0->atom[i].type = j;
 164         r0->atom[i].m    = get_atomtype_massA(j, atype);
 165         i++;
 166     }
 167     r0->natom = i;
 168     srenew(r0->atom, i);
 169     srenew(r0->atomname, i);
 170     srenew(r0->cgnr, i);
 171
 172     return TRUE;
 173 }
 174
 175 gmx_bool read_bondeds(int bt, FILE *in, char *line, t_restp *rtp)
 176 {
 177     char str[STRLEN];
 178     int  j, n, ni, maxrb;
 179
 180     maxrb = rtp->rb[bt].nb;
 181     while (get_a_line(in, line, STRLEN) && (strchr(line, '[') == NULL))
 182     {
 183         if (rtp->rb[bt].nb >= maxrb)
 184         {
 185             maxrb += 100;
 186             srenew(rtp->rb[bt].b, maxrb);
 187         }
 188         n = 0;
 189         for (j = 0; j < btsNiatoms[bt]; j++)
 190         {
 191             if (sscanf(line+n, "%s%n", str, &ni) == 1)
 192             {
 193                 rtp->rb[bt].b[rtp->rb[bt].nb].a[j] = strdup(str);
 194             }
 195             else
 196             {
 197                 return FALSE;
 198             }
 199             n += ni;
 200         }
 201         for (; j < MAXATOMLIST; j++)
 202         {
 203             rtp->rb[bt].b[rtp->rb[bt].nb].a[j] = NULL;
 204         }
 205         while (isspace(line[n]))
 206         {
 207             n++;
 208         }
 209         rtrim(line+n);
 210         rtp->rb[bt].b[rtp->rb[bt].nb].s = strdup(line+n);
 211         rtp->rb[bt].nb++;
 212     }
 213     /* give back unused memory */
 214     srenew(rtp->rb[bt].b, rtp->rb[bt].nb);
 215
 216     return TRUE;
 217 }
 218
 219 static void print_resbondeds(FILE *out, int bt, t_restp *rtp)
 220 {
 221     int i, j;
 222
 223     if (rtp->rb[bt].nb)
 224     {
 225         fprintf(out, " [ %s ]\n", btsNames[bt]);
 226
 227         for (i = 0; i < rtp->rb[bt].nb; i++)
 228         {
 229             for (j = 0; j < btsNiatoms[bt]; j++)
 230             {
 231                 fprintf(out, "%6s ", rtp->rb[bt].b[i].a[j]);
 232             }
 233             if (rtp->rb[bt].b[i].s[0])
 234             {
 235                 fprintf(out, "    %s", rtp->rb[bt].b[i].s);
 236             }
 237             fprintf(out, "\n");
 238         }
 239     }
 240 }
 241
 242 static void check_rtp(int nrtp, t_restp rtp[], char *libfn)
 243 {
 244     int i;
 245
 246     /* check for double entries, assuming list is already sorted */
 247     for (i = 1; (i < nrtp); i++)
 248     {
 249         if (gmx_strcasecmp(rtp[i-1].resname, rtp[i].resname) == 0)
 250         {
 251             fprintf(stderr, "WARNING double entry %s in file %s\n",
 252                     rtp[i].resname, libfn);
 253         }
 254     }
 255 }
 256
 257 static int comprtp(const void *a, const void *b)
 258 {
 259     t_restp *ra, *rb;
 260
 261     ra = (t_restp *)a;
 262     rb = (t_restp *)b;
 263
 264     return gmx_strcasecmp(ra->resname, rb->resname);
 265 }
 266
 267 int get_bt(char* header)
 268 {
 269     int i;
 270
 271     for (i = 0; i < ebtsNR; i++)
 272     {
 273         if (gmx_strcasecmp(btsNames[i], header) == 0)
 274         {
 275             return i;
 276         }
 277     }
 278     return NOTSET;
 279 }
 280
 281 void clear_t_restp(t_restp *rrtp)
 282 {
 283     memset((void *)rrtp, 0, sizeof(t_restp));
 284 }
 285
 286 /* print all the ebtsNR type numbers */
 287 void print_resall_header(FILE *out, t_restp rtp[])
 288 {
 289     fprintf(out, "[ bondedtypes ]\n");
 290     fprintf(out, "; bonds  angles  dihedrals  impropers all_dihedrals nr_exclusions  HH14  remove_dih\n");
 291     fprintf(out, " %5d  %6d  %9d  %9d  %14d  %14d %14d %14d\n\n",
 292             rtp[0].rb[0].type,
 293             rtp[0].rb[1].type,
 294             rtp[0].rb[2].type,
 295             rtp[0].rb[3].type,
 296             rtp[0].bKeepAllGeneratedDihedrals,
 297             rtp[0].nrexcl,
 298             rtp[0].bGenerateHH14Interactions,
 299             rtp[0].bRemoveDihedralIfWithImproper);
 300 }
 301
 302 void print_resall(FILE *out, int nrtp, t_restp rtp[],
 303                   gpp_atomtype_t atype)
 304 {
 305     int i, bt;
 306
 307     if (nrtp == 0)
 308     {
 309         return;
 310     }
 311
 312     print_resall_header(out, rtp);
 313
 314     for (i = 0; i < nrtp; i++)
 315     {
 316         if (rtp[i].natom > 0)
 317         {
 318             print_resatoms(out, atype, &rtp[i]);
 319             for (bt = 0; bt < ebtsNR; bt++)
 320             {
 321                 print_resbondeds(out, bt, &rtp[i]);
 322             }
 323         }
 324     }
 325 }
 326
 327 void read_resall(char *rrdb, int *nrtpptr, t_restp **rtp,
 328                  gpp_atomtype_t atype, t_symtab *tab,
 329                  gmx_bool bAllowOverrideRTP)
 330 {
 331     FILE         *in;
 332     char          filebase[STRLEN], *ptr, line[STRLEN], header[STRLEN];
 333     int           i, nrtp, maxrtp, bt, nparam;
 334     int           dum1, dum2, dum3;
 335     t_restp      *rrtp, *header_settings;
 336     gmx_bool      bNextResidue, bError;
 337     int           firstrtp;
 338
 339     fflib_filename_base(rrdb, filebase, STRLEN);
 340
 341     in = fflib_open(rrdb);
 342
 343     if (debug)
 344     {
 345         fprintf(debug, "%9s %5s", "Residue", "atoms");
 346         for (i = 0; i < ebtsNR; i++)
 347         {
 348             fprintf(debug, " %10s", btsNames[i]);
 349         }
 350         fprintf(debug, "\n");
 351     }
 352     snew(header_settings, 1);
 353
 354     /* these bonded parameters will overwritten be when  *
 355      * there is a [ bondedtypes ] entry in the .rtp file */
 356     header_settings->rb[ebtsBONDS].type  = 1; /* normal bonds     */
 357     header_settings->rb[ebtsANGLES].type = 1; /* normal angles    */
 358     header_settings->rb[ebtsPDIHS].type  = 1; /* normal dihedrals */
 359     header_settings->rb[ebtsIDIHS].type  = 2; /* normal impropers */
 360     header_settings->rb[ebtsEXCLS].type  = 1; /* normal exclusions */
 361     header_settings->rb[ebtsCMAP].type   = 1; /* normal cmap torsions */
 362
 363     header_settings->bKeepAllGeneratedDihedrals    = FALSE;
 364     header_settings->nrexcl                        = 3;
 365     header_settings->bGenerateHH14Interactions     = TRUE;
 366     header_settings->bRemoveDihedralIfWithImproper = TRUE;
 367
 368     /* Column 5 & 6 aren't really bonded types, but we include
 369      * them here to avoid introducing a new section:
 370      * Column 5 : This controls the generation of dihedrals from the bonding.
 371      *            All possible dihedrals are generated automatically. A value of
 372      *            1 here means that all these are retained. A value of
 373      *            0 here requires generated dihedrals be removed if
 374      *              * there are any dihedrals on the same central atoms
 375      *                specified in the residue topology, or
 376      *              * there are other identical generated dihedrals
 377      *                sharing the same central atoms, or
 378      *              * there are other generated dihedrals sharing the
 379      *                same central bond that have fewer hydrogen atoms
 380      * Column 6: Number of bonded neighbors to exclude.
 381      * Column 7: Generate 1,4 interactions between two hydrogen atoms
 382      * Column 8: Remove proper dihedrals if centered on the same bond
 383      *           as an improper dihedral
 384      */
 385     get_a_line(in, line, STRLEN);
 386     if (!get_header(line, header))
 387     {
 388         gmx_fatal(FARGS, "in .rtp file at line:\n%s\n", line);
 389     }
 390     if (gmx_strncasecmp("bondedtypes", header, 5) == 0)
 391     {
 392         get_a_line(in, line, STRLEN);
 393         if ((nparam = sscanf(line, "%d %d %d %d %d %d %d %d",
 394                              &header_settings->rb[ebtsBONDS].type, &header_settings->rb[ebtsANGLES].type,
 395                              &header_settings->rb[ebtsPDIHS].type, &header_settings->rb[ebtsIDIHS].type,
 396                              &dum1, &header_settings->nrexcl, &dum2, &dum3)) < 4)
 397         {
 398             gmx_fatal(FARGS, "need 4 to 8 parameters in the header of .rtp file %s at line:\n%s\n", rrdb, line);
 399         }
 400         header_settings->bKeepAllGeneratedDihedrals    = (dum1 != 0);
 401         header_settings->bGenerateHH14Interactions     = (dum2 != 0);
 402         header_settings->bRemoveDihedralIfWithImproper = (dum3 != 0);
 403         get_a_line(in, line, STRLEN);
 404         if (nparam < 5)
 405         {
 406             fprintf(stderr, "Using default: not generating all possible dihedrals\n");
 407             header_settings->bKeepAllGeneratedDihedrals = FALSE;
 408         }
 409         if (nparam < 6)
 410         {
 411             fprintf(stderr, "Using default: excluding 3 bonded neighbors\n");
 412             header_settings->nrexcl = 3;
 413         }
 414         if (nparam < 7)
 415         {
 416             fprintf(stderr, "Using default: generating 1,4 H--H interactions\n");
 417             header_settings->bGenerateHH14Interactions = TRUE;
 418         }
 419         if (nparam < 8)
 420         {
 421             fprintf(stderr, "Using default: removing proper dihedrals found on the same bond as a proper dihedral\n");
 422             header_settings->bRemoveDihedralIfWithImproper = TRUE;
 423         }
 424     }
 425     else
 426     {
 427         fprintf(stderr,
 428                 "Reading .rtp file without '[ bondedtypes ]' directive,\n"
 429                 "Will proceed as if the entry was:\n");
 430         print_resall_header(stderr, header_settings);
 431     }
 432     /* We don't know the current size of rrtp, but simply realloc immediately */
 433     nrtp   = *nrtpptr;
 434     rrtp   = *rtp;
 435     maxrtp = nrtp;
 436     while (!feof(in))
 437     {
 438         if (nrtp >= maxrtp)
 439         {
 440             maxrtp += 100;
 441             srenew(rrtp, maxrtp);
 442         }
 443
 444         /* Initialise rtp entry structure */
 445         rrtp[nrtp] = *header_settings;
 446         if (!get_header(line, header))
 447         {
 448             gmx_fatal(FARGS, "in .rtp file at line:\n%s\n", line);
 449         }
 450         rrtp[nrtp].resname  = strdup(header);
 451         rrtp[nrtp].filebase = strdup(filebase);
 452
 453         get_a_line(in, line, STRLEN);
 454         bError       = FALSE;
 455         bNextResidue = FALSE;
 456         do
 457         {
 458             if (!get_header(line, header))
 459             {
 460                 bError = TRUE;
 461             }
 462             else
 463             {
 464                 bt = get_bt(header);
 465                 if (bt != NOTSET)
 466                 {
 467                     /* header is an bonded directive */
 468                     bError = !read_bondeds(bt, in, line, &rrtp[nrtp]);
 469                 }
 470                 else if (gmx_strncasecmp("atoms", header, 5) == 0)
 471                 {
 472                     /* header is the atoms directive */
 473                     bError = !read_atoms(in, line, &(rrtp[nrtp]), tab, atype);
 474                 }
 475                 else
 476                 {
 477                     /* else header must be a residue name */
 478                     bNextResidue = TRUE;
 479                 }
 480             }
 481             if (bError)
 482             {
 483                 gmx_fatal(FARGS, "in .rtp file in residue %s at line:\n%s\n",
 484                           rrtp[nrtp].resname, line);
 485             }
 486         }
 487         while (!feof(in) && !bNextResidue);
 488
 489         if (rrtp[nrtp].natom == 0)
 490         {
 491             gmx_fatal(FARGS, "No atoms found in .rtp file in residue %s\n",
 492                       rrtp[nrtp].resname);
 493         }
 494         if (debug)
 495         {
 496             fprintf(debug, "%3d %5s %5d",
 497                     nrtp+1, rrtp[nrtp].resname, rrtp[nrtp].natom);
 498             for (i = 0; i < ebtsNR; i++)
 499             {
 500                 fprintf(debug, " %10d", rrtp[nrtp].rb[i].nb);
 501             }
 502             fprintf(debug, "\n");
 503         }
 504
 505         firstrtp = -1;
 506         for (i = 0; i < nrtp; i++)
 507         {
 508             if (gmx_strcasecmp(rrtp[i].resname, rrtp[nrtp].resname) == 0)
 509             {
 510                 firstrtp = i;
 511             }
 512         }
 513
 514         if (firstrtp == -1)
 515         {
 516             nrtp++;
 517             fprintf(stderr, "\rResidue %d", nrtp);
 518         }
 519         else
 520         {
 521             if (firstrtp >= *nrtpptr)
 522             {
 523                 gmx_fatal(FARGS, "Found a second entry for '%s' in '%s'",
 524                           rrtp[nrtp].resname, rrdb);
 525             }
 526             if (bAllowOverrideRTP)
 527             {
 528                 fprintf(stderr, "\n");
 529                 fprintf(stderr, "Found another rtp entry for '%s' in '%s', ignoring this entry and keeping the one from '%s.rtp'\n",
 530                         rrtp[nrtp].resname, rrdb, rrtp[firstrtp].filebase);
 531                 /* We should free all the data for this entry.
 532                  * The current code gives a lot of dangling pointers.
 533                  */
 534                 clear_t_restp(&rrtp[nrtp]);
 535             }
 536             else
 537             {
 538                 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);
 539             }
 540         }
 541     }
 542     gmx_ffclose(in);
 543     /* give back unused memory */
 544     srenew(rrtp, nrtp);
 545
 546     fprintf(stderr, "\nSorting it all out...\n");
 547     qsort(rrtp, nrtp, (size_t)sizeof(rrtp[0]), comprtp);
 548
 549     check_rtp(nrtp, rrtp, rrdb);
 550
 551     *nrtpptr = nrtp;
 552     *rtp     = rrtp;
 553 }
 554
 555 /************************************************************
 556  *
 557  *                  SEARCH   ROUTINES
 558  *
 559  ***********************************************************/
 560 static gmx_bool is_sign(char c)
 561 {
 562     return (c == '+' || c == '-');
 563 }
 564
 565 /* Compares if the strins match except for a sign at the end
 566  * of (only) one of the two.
 567  */
 568 static int neq_str_sign(const char *a1, const char *a2)
 569 {
 570     int l1, l2, lm;
 571
 572     l1 = (int)strlen(a1);
 573     l2 = (int)strlen(a2);
 574     lm = min(l1, l2);
 575
 576     if (lm >= 1 &&
 577         ((l1 == l2+1 && is_sign(a1[l1-1])) ||
 578          (l2 == l1+1 && is_sign(a2[l2-1]))) &&
 579         gmx_strncasecmp(a1, a2, lm) == 0)
 580     {
 581         return lm;
 582     }
 583     else
 584     {
 585         return 0;
 586     }
 587 }
 588
 589 char *search_rtp(const char *key, int nrtp, t_restp rtp[])
 590 {
 591     int  i, n, nbest, best, besti;
 592     char bestbuf[STRLEN];
 593
 594     nbest =  0;
 595     besti = -1;
 596     /* We want to match at least one character */
 597     best  =  1;
 598     for (i = 0; (i < nrtp); i++)
 599     {
 600         if (gmx_strcasecmp(key, rtp[i].resname) == 0)
 601         {
 602             besti = i;
 603             nbest = 1;
 604             break;
 605         }
 606         else
 607         {
 608             /* Allow a mismatch of at most a sign character (with warning) */
 609             n = neq_str_sign(key, rtp[i].resname);
 610             if (n >= best &&
 611                 n+1 >= (int)strlen(key) &&
 612                 n+1 >= (int)strlen(rtp[i].resname))
 613             {
 614                 if (n == best)
 615                 {
 616                     if (nbest == 1)
 617                     {
 618                         strcpy(bestbuf, rtp[besti].resname);
 619                     }
 620                     if (nbest >= 1)
 621                     {
 622                         strcat(bestbuf, " or ");
 623                         strcat(bestbuf, rtp[i].resname);
 624                     }
 625                 }
 626                 else
 627                 {
 628                     nbest = 0;
 629                 }
 630                 besti = i;
 631                 best  = n;
 632                 nbest++;
 633             }
 634         }
 635     }
 636     if (nbest > 1)
 637     {
 638         gmx_fatal(FARGS, "Residue '%s' not found in residue topology database, looks a bit like %s", key, bestbuf);
 639     }
 640     else if (besti == -1)
 641     {
 642         gmx_fatal(FARGS, "Residue '%s' not found in residue topology database", key);
 643     }
 644     if (gmx_strcasecmp(rtp[besti].resname, key) != 0)
 645     {
 646         fprintf(stderr,
 647                 "\nWARNING: '%s' not found in residue topology database, "
 648                 "trying to use '%s'\n\n", key, rtp[besti].resname);
 649     }
 650
 651     return rtp[besti].resname;
 652 }
 653
 654 t_restp *get_restp(const char *rtpname, int nrtp, t_restp rtp[])
 655 {
 656     int  i;
 657
 658     i = 0;
 659     while (i < nrtp && gmx_strcasecmp(rtpname, rtp[i].resname) != 0)
 660     {
 661         i++;
 662     }
 663     if (i >= nrtp)
 664     {
 665         /* This should never happen, since search_rtp should have been called
 666          * before calling get_restp.
 667          */
 668         gmx_fatal(FARGS, "Residue type '%s' not found in residue topology database", rtpname);
 669     }
 670
 671     return &rtp[i];
 672 }