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