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