src/gromacs/gmxpreprocess/ter_db.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include "sysstuff.h"
  42 #include "smalloc.h"
  43 #include "typedefs.h"
  44 #include "symtab.h"
  45 #include "gromacs/fileio/futil.h"
  46 #include "resall.h"
  47 #include "h_db.h"
  48 #include "string2.h"
  49 #include "gmx_fatal.h"
  50 #include "ter_db.h"
  51 #include "toputil.h"
  52 #include "gromacs/fileio/gmxfio.h"
  53 #include "fflibutil.h"
  54
  55 #include "gromacs/fileio/strdb.h"
  56
  57 /* use bonded types definitions in hackblock.h */
  58 #define ekwRepl ebtsNR+1
  59 #define ekwAdd  ebtsNR+2
  60 #define ekwDel  ebtsNR+3
  61 #define ekwNR   3
  62 const char *kw_names[ekwNR] = {
  63     "replace", "add", "delete"
  64 };
  65
  66 int find_kw(char *keyw)
  67 {
  68     int i;
  69
  70     for (i = 0; i < ebtsNR; i++)
  71     {
  72         if (gmx_strcasecmp(btsNames[i], keyw) == 0)
  73         {
  74             return i;
  75         }
  76     }
  77     for (i = 0; i < ekwNR; i++)
  78     {
  79         if (gmx_strcasecmp(kw_names[i], keyw) == 0)
  80         {
  81             return ebtsNR + 1 + i;
  82         }
  83     }
  84
  85     return NOTSET;
  86 }
  87
  88 #define FATAL() gmx_fatal(FARGS, "Reading Termini Database: not enough items on line\n%s", line)
  89
  90 static void read_atom(char *line, gmx_bool bAdd,
  91                       char **nname, t_atom *a, gpp_atomtype_t atype, int *cgnr)
  92 {
  93     int    nr, i;
  94     char   buf[5][30], type[12];
  95     double m, q;
  96
  97     /* This code is messy, because of support for different formats:
  98      * for replace: [new name] <atom type> <m> <q> [cgnr (old format)]
  99      * for add:                <atom type> <m> <q> [cgnr]
 100      */
 101     nr = sscanf(line, "%s %s %s %s %s", buf[0], buf[1], buf[2], buf[3], buf[4]);
 102
 103     /* Here there an ambiguity due to the old replace format with cgnr,
 104      * which was read for years, but ignored in the rest of the code.
 105      * We have to assume that the atom type does not start with a digit
 106      * to make a line with 4 entries uniquely interpretable.
 107      */
 108     if (!bAdd && nr == 4 && isdigit(buf[1][0]))
 109     {
 110         nr = 3;
 111     }
 112
 113     if (nr < 3 || nr > 4)
 114     {
 115         gmx_fatal(FARGS, "Reading Termini Database: expected %d or %d items of atom data in stead of %d on line\n%s", 3, 4, nr, line);
 116     }
 117     i = 0;
 118     if (!bAdd)
 119     {
 120         if (nr == 4)
 121         {
 122             *nname = strdup(buf[i++]);
 123         }
 124         else
 125         {
 126             *nname = NULL;
 127         }
 128     }
 129     a->type = get_atomtype_type(buf[i++], atype);
 130     sscanf(buf[i++], "%lf", &m);
 131     a->m = m;
 132     sscanf(buf[i++], "%lf", &q);
 133     a->q = q;
 134     if (bAdd && nr == 4)
 135     {
 136         sscanf(buf[i++], "%d", cgnr);
 137     }
 138     else
 139     {
 140         *cgnr = NOTSET;
 141     }
 142 }
 143
 144 static void print_atom(FILE *out, t_atom *a, gpp_atomtype_t atype)
 145 {
 146     fprintf(out, "\t%s\t%g\t%g\n",
 147             get_atomtype_name(a->type, atype), a->m, a->q);
 148 }
 149
 150 static void print_ter_db(const char *ff, char C, int nb, t_hackblock tb[],
 151                          gpp_atomtype_t atype)
 152 {
 153     FILE *out;
 154     int   i, j, k, bt, nrepl, nadd, ndel;
 155     char  buf[STRLEN], nname[STRLEN];
 156
 157     sprintf(buf, "%s-%c.tdb", ff, C);
 158     out = gmx_fio_fopen(buf, "w");
 159
 160     for (i = 0; (i < nb); i++)
 161     {
 162         fprintf(out, "[ %s ]\n", tb[i].name);
 163
 164         /* first count: */
 165         nrepl = 0;
 166         nadd  = 0;
 167         ndel  = 0;
 168         for (j = 0; j < tb[i].nhack; j++)
 169         {
 170             if (tb[i].hack[j].oname != NULL && tb[i].hack[j].nname != NULL)
 171             {
 172                 nrepl++;
 173             }
 174             else if (tb[i].hack[j].oname == NULL && tb[i].hack[j].nname != NULL)
 175             {
 176                 nadd++;
 177             }
 178             else if (tb[i].hack[j].oname != NULL && tb[i].hack[j].nname == NULL)
 179             {
 180                 ndel++;
 181             }
 182             else if (tb[i].hack[j].oname == NULL && tb[i].hack[j].nname == NULL)
 183             {
 184                 gmx_fatal(FARGS, "invalid hack (%s) in termini database", tb[i].name);
 185             }
 186         }
 187         if (nrepl)
 188         {
 189             fprintf(out, "[ %s ]\n", kw_names[ekwRepl-ebtsNR-1]);
 190             for (j = 0; j < tb[i].nhack; j++)
 191             {
 192                 if (tb[i].hack[j].oname != NULL && tb[i].hack[j].nname != NULL)
 193                 {
 194                     fprintf(out, "%s\t", tb[i].hack[j].oname);
 195                     print_atom(out, tb[i].hack[j].atom, atype);
 196                 }
 197             }
 198         }
 199         if (nadd)
 200         {
 201             fprintf(out, "[ %s ]\n", kw_names[ekwAdd-ebtsNR-1]);
 202             for (j = 0; j < tb[i].nhack; j++)
 203             {
 204                 if (tb[i].hack[j].oname == NULL && tb[i].hack[j].nname != NULL)
 205                 {
 206                     print_ab(out, &(tb[i].hack[j]), tb[i].hack[j].nname);
 207                     print_atom(out, tb[i].hack[j].atom, atype);
 208                 }
 209             }
 210         }
 211         if (ndel)
 212         {
 213             fprintf(out, "[ %s ]\n", kw_names[ekwDel-ebtsNR-1]);
 214             for (j = 0; j < tb[i].nhack; j++)
 215             {
 216                 if (tb[i].hack[j].oname != NULL && tb[i].hack[j].nname == NULL)
 217                 {
 218                     fprintf(out, "%s\n", tb[i].hack[j].oname);
 219                 }
 220             }
 221         }
 222         for (bt = 0; bt < ebtsNR; bt++)
 223         {
 224             if (tb[i].rb[bt].nb)
 225             {
 226                 fprintf(out, "[ %s ]\n", btsNames[bt]);
 227                 for (j = 0; j < tb[i].rb[bt].nb; j++)
 228                 {
 229                     for (k = 0; k < btsNiatoms[bt]; k++)
 230                     {
 231                         fprintf(out, "%s%s", k ? "\t" : "", tb[i].rb[bt].b[j].a[k]);
 232                     }
 233                     if (tb[i].rb[bt].b[j].s)
 234                     {
 235                         fprintf(out, "\t%s", tb[i].rb[bt].b[j].s);
 236                     }
 237                     fprintf(out, "\n");
 238                 }
 239             }
 240         }
 241         fprintf(out, "\n");
 242     }
 243     gmx_fio_fclose(out);
 244 }
 245
 246 static void read_ter_db_file(char *fn,
 247                              int *ntbptr, t_hackblock **tbptr,
 248                              gpp_atomtype_t atype)
 249 {
 250     char         filebase[STRLEN], *ptr;
 251     FILE        *in;
 252     char         header[STRLEN], buf[STRLEN], line[STRLEN];
 253     t_hackblock *tb;
 254     int          i, j, n, ni, kwnr, nb, maxnb, nh;
 255
 256     fflib_filename_base(fn, filebase, STRLEN);
 257     /* Remove the C/N termini extension */
 258     ptr = strrchr(filebase, '.');
 259     if (ptr != NULL)
 260     {
 261         ptr[0] = '\0';
 262     }
 263
 264     in = fflib_open(fn);
 265     if (debug)
 266     {
 267         fprintf(debug, "Opened %s\n", fn);
 268     }
 269
 270     tb    = *tbptr;
 271     nb    = *ntbptr - 1;
 272     maxnb = 0;
 273     kwnr  = NOTSET;
 274     get_a_line(in, line, STRLEN);
 275     while (!feof(in))
 276     {
 277         if (get_header(line, header))
 278         {
 279             /* this is a new block, or a new keyword */
 280             kwnr = find_kw(header);
 281
 282             if (kwnr == NOTSET)
 283             {
 284                 nb++;
 285                 /* here starts a new block */
 286                 if (nb >= maxnb)
 287                 {
 288                     maxnb = nb + 100;
 289                     srenew(tb, maxnb);
 290                 }
 291                 clear_t_hackblock(&tb[nb]);
 292                 tb[nb].name     = strdup(header);
 293                 tb[nb].filebase = strdup(filebase);
 294             }
 295         }
 296         else
 297         {
 298             if (nb < 0)
 299             {
 300                 gmx_fatal(FARGS, "reading termini database: "
 301                           "directive expected before line:\n%s\n"
 302                           "This might be a file in an old format.", line);
 303             }
 304             /* this is not a header, so it must be data */
 305             if (kwnr >= ebtsNR)
 306             {
 307                 /* this is a hack: add/rename/delete atoms */
 308                 /* make space for hacks */
 309                 if (tb[nb].nhack >= tb[nb].maxhack)
 310                 {
 311                     tb[nb].maxhack += 10;
 312                     srenew(tb[nb].hack, tb[nb].maxhack);
 313                 }
 314                 nh = tb[nb].nhack;
 315                 clear_t_hack(&(tb[nb].hack[nh]));
 316                 for (i = 0; i < 4; i++)
 317                 {
 318                     tb[nb].hack[nh].a[i] = NULL;
 319                 }
 320                 tb[nb].nhack++;
 321
 322                 /* get data */
 323                 n = 0;
 324                 if (kwnr == ekwRepl || kwnr == ekwDel)
 325                 {
 326                     if (sscanf(line, "%s%n", buf, &n) != 1)
 327                     {
 328                         gmx_fatal(FARGS, "Reading Termini Database '%s': "
 329                                   "expected atom name on line\n%s", fn, line);
 330                     }
 331                     tb[nb].hack[nh].oname = strdup(buf);
 332                     /* we only replace or delete one atom at a time */
 333                     tb[nb].hack[nh].nr = 1;
 334                 }
 335                 else if (kwnr == ekwAdd)
 336                 {
 337                     read_ab(line, fn, &(tb[nb].hack[nh]));
 338                     get_a_line(in, line, STRLEN);
 339                 }
 340                 else
 341                 {
 342                     gmx_fatal(FARGS, "unimplemented keyword number %d (%s:%d)",
 343                               kwnr, __FILE__, __LINE__);
 344                 }
 345                 if (kwnr == ekwRepl || kwnr == ekwAdd)
 346                 {
 347                     snew(tb[nb].hack[nh].atom, 1);
 348                     read_atom(line+n, kwnr == ekwAdd,
 349                               &tb[nb].hack[nh].nname, tb[nb].hack[nh].atom, atype,
 350                               &tb[nb].hack[nh].cgnr);
 351                     if (tb[nb].hack[nh].nname == NULL)
 352                     {
 353                         if (tb[nb].hack[nh].oname != NULL)
 354                         {
 355                             tb[nb].hack[nh].nname = strdup(tb[nb].hack[nh].oname);
 356                         }
 357                         else
 358                         {
 359                             gmx_fatal(FARGS, "Reading Termini Database '%s': don't know which name the new atom should have on line\n%s", fn, line);
 360                         }
 361                     }
 362                 }
 363             }
 364             else if (kwnr >= 0 && kwnr < ebtsNR)
 365             {
 366                 /* this is bonded data: bonds, angles, dihedrals or impropers */
 367                 srenew(tb[nb].rb[kwnr].b, tb[nb].rb[kwnr].nb+1);
 368                 n = 0;
 369                 for (j = 0; j < btsNiatoms[kwnr]; j++)
 370                 {
 371                     if (sscanf(line+n, "%s%n", buf, &ni) == 1)
 372                     {
 373                         tb[nb].rb[kwnr].b[tb[nb].rb[kwnr].nb].a[j] = strdup(buf);
 374                     }
 375                     else
 376                     {
 377                         gmx_fatal(FARGS, "Reading Termini Database '%s': expected %d atom names (found %d) on line\n%s", fn, btsNiatoms[kwnr], j-1, line);
 378                     }
 379                     n += ni;
 380                 }
 381                 for (; j < MAXATOMLIST; j++)
 382                 {
 383                     tb[nb].rb[kwnr].b[tb[nb].rb[kwnr].nb].a[j] = NULL;
 384                 }
 385                 strcpy(buf, "");
 386                 sscanf(line+n, "%s", buf);
 387                 tb[nb].rb[kwnr].b[tb[nb].rb[kwnr].nb].s = strdup(buf);
 388                 tb[nb].rb[kwnr].nb++;
 389             }
 390             else
 391             {
 392                 gmx_fatal(FARGS, "Reading Termini Database: Expecting a header at line\n"
 393                           "%s", line);
 394             }
 395         }
 396         get_a_line(in, line, STRLEN);
 397     }
 398     nb++;
 399     srenew(tb, nb);
 400
 401     gmx_ffclose(in);
 402
 403     *ntbptr = nb;
 404     *tbptr  = tb;
 405 }
 406
 407 int read_ter_db(const char *ffdir, char ter,
 408                 t_hackblock **tbptr, gpp_atomtype_t atype)
 409 {
 410     char   ext[STRLEN];
 411     int    ntdbf, f;
 412     char **tdbf;
 413     int    ntb;
 414
 415     sprintf(ext, ".%c.tdb", ter);
 416
 417     /* Search for termini database files.
 418      * Do not generate an error when none are found.
 419      */
 420     ntdbf  = fflib_search_file_end(ffdir, ext, FALSE, &tdbf);
 421     ntb    = 0;
 422     *tbptr = NULL;
 423     for (f = 0; f < ntdbf; f++)
 424     {
 425         read_ter_db_file(tdbf[f], &ntb, tbptr, atype);
 426         sfree(tdbf[f]);
 427     }
 428     sfree(tdbf);
 429
 430     if (debug)
 431     {
 432         print_ter_db("new", ter, ntb, *tbptr, atype);
 433     }
 434
 435     return ntb;
 436 }
 437
 438 t_hackblock **filter_ter(int nrtp, t_restp rtp[],
 439                          int nb, t_hackblock tb[],
 440                          const char *resname,
 441                          const char *rtpname,
 442                          int *nret)
 443 {
 444     /* Since some force fields (e.g. OPLS) needs different
 445      * atomtypes for different residues there could be a lot
 446      * of entries in the databases for specific residues
 447      * (e.g. GLY-NH3+, SER-NH3+, ALA-NH3+).
 448      *
 449      * To reduce the database size, we assume that a terminus specifier liker
 450      *
 451      * [ GLY|SER|ALA-NH3+ ]
 452      *
 453      * would cover all of the three residue types above.
 454      * Wildcards (*,?) are not OK. Don't worry about e.g. GLU vs. GLUH since
 455      * pdb2gmx only uses the first 3 letters when calling this routine.
 456      *
 457      * To automate this, this routines scans a list of termini
 458      * for the residue name "resname" and returns an allocated list of
 459      * pointers to the termini that could be applied to the
 460      * residue in question. The variable pointed to by nret will
 461      * contain the number of valid pointers in the list.
 462      * Remember to free the list when you are done with it...
 463      */
 464
 465     t_restp     *   restp;
 466     int             i, j, n, len, none_idx;
 467     gmx_bool        found;
 468     char           *rtpname_match, *s, *s2, *c;
 469     t_hackblock   **list;
 470
 471     rtpname_match = search_rtp(rtpname, nrtp, rtp);
 472     restp         = get_restp(rtpname_match, nrtp, rtp);
 473
 474     n    = 0;
 475     list = NULL;
 476
 477     for (i = 0; i < nb; i++)
 478     {
 479         s     = tb[i].name;
 480         found = FALSE;
 481         do
 482         {
 483             /* The residue name should appear in a tdb file with the same base name
 484              * as the file containing the rtp entry.
 485              * This makes termini selection for different molecule types
 486              * much cleaner.
 487              */
 488             if (gmx_strcasecmp(restp->filebase, tb[i].filebase) == 0 &&
 489                 gmx_strncasecmp(resname, s, 3) == 0)
 490             {
 491                 found = TRUE;
 492                 srenew(list, n+1);
 493                 list[n] = &(tb[i]);
 494                 n++;
 495             }
 496             else
 497             {
 498                 /* advance to next |-separated field */
 499                 s = strchr(s, '|');
 500                 if (s != NULL)
 501                 {
 502                     s++;
 503                 }
 504             }
 505         }
 506         while (!found && s != NULL);
 507     }
 508
 509     /* All residue-specific termini have been added. See if there
 510      * are some generic ones by searching for the occurence of
 511      * '-' in the name prior to the last position (which indicates charge).
 512      * The [ None ] alternative is special since we don't want that
 513      * to be the default, so we put it last in the list we return.
 514      */
 515     none_idx = -1;
 516     for (i = 0; i < nb; i++)
 517     {
 518         s = tb[i].name;
 519         /* The residue name should appear in a tdb file with the same base name
 520          * as the file containing the rtp entry.
 521          * This makes termini selection for different molecule types
 522          * much cleaner.
 523          */
 524         if (gmx_strcasecmp(restp->filebase, tb[i].filebase) == 0)
 525         {
 526             if (!gmx_strcasecmp("None", s))
 527             {
 528                 none_idx = i;
 529             }
 530             else
 531             {
 532                 c = strchr(s, '-');
 533                 if (c == NULL || ((c-s+1) == strlen(s)))
 534                 {
 535                     /* Check that we haven't already added a residue-specific version
 536                      * of this terminus.
 537                      */
 538                     for (j = 0; j < n && strstr((*list[j]).name, s) == NULL; j++)
 539                     {
 540                         ;
 541                     }
 542                     if (j == n)
 543                     {
 544                         srenew(list, n+1);
 545                         list[n] = &(tb[i]);
 546                         n++;
 547                     }
 548                 }
 549             }
 550         }
 551     }
 552     if (none_idx >= 0)
 553     {
 554         srenew(list, n+1);
 555         list[n] = &(tb[none_idx]);
 556         n++;
 557     }
 558
 559     *nret = n;
 560     return list;
 561 }
 562
 563
 564 t_hackblock *choose_ter(int nb, t_hackblock **tb, const char *title)
 565 {
 566     int i, sel, ret;
 567
 568     printf("%s\n", title);
 569     for (i = 0; (i < nb); i++)
 570     {
 571         char *advice_string = "";
 572         if (0 == gmx_wcmatch("*ZWITTERION*", (*tb[i]).name))
 573         {
 574             advice_string = " (only use with zwitterions containing exactly one residue)";
 575         }
 576         printf("%2d: %s%s\n", i, (*tb[i]).name, advice_string);
 577     }
 578     do
 579     {
 580         ret = fscanf(stdin, "%d", &sel);
 581     }
 582     while ((ret != 1) || (sel < 0) || (sel >= nb));
 583
 584     return tb[sel];
 585 }