src/gromacs/gmxpreprocess/resall.cpp

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