src/gromacs/gmxpreprocess/resall.c

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