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