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