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