src/gromacs/gmxpreprocess/ter_db.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,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 "ter_db.h"
  40
  41 #include <cctype>
  42 #include <cstring>
  43
  44 #include <algorithm>
  45 #include <string>
  46 #include <vector>
  47
  48 #include "gromacs/fileio/gmxfio.h"
  49 #include "gromacs/gmxpreprocess/fflibutil.h"
  50 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
  51 #include "gromacs/gmxpreprocess/grompp_impl.h"
  52 #include "gromacs/gmxpreprocess/h_db.h"
  53 #include "gromacs/gmxpreprocess/notset.h"
  54 #include "gromacs/gmxpreprocess/toputil.h"
  55 #include "gromacs/utility/cstringutil.h"
  56 #include "gromacs/utility/fatalerror.h"
  57 #include "gromacs/utility/futil.h"
  58 #include "gromacs/utility/smalloc.h"
  59 #include "gromacs/utility/strdb.h"
  60
  61 #include "hackblock.h"
  62 #include "resall.h"
  63
  64 /* use bonded types definitions in hackblock.h */
  65 #define ekwRepl (ebtsNR+1)
  66 #define ekwAdd  (ebtsNR+2)
  67 #define ekwDel  (ebtsNR+3)
  68 #define ekwNR   3
  69 static const char *kw_names[ekwNR] = {
  70     "replace", "add", "delete"
  71 };
  72
  73 static int find_kw(char *keyw)
  74 {
  75     int i;
  76
  77     for (i = 0; i < ebtsNR; i++)
  78     {
  79         if (gmx_strcasecmp(btsNames[i], keyw) == 0)
  80         {
  81             return i;
  82         }
  83     }
  84     for (i = 0; i < ekwNR; i++)
  85     {
  86         if (gmx_strcasecmp(kw_names[i], keyw) == 0)
  87         {
  88             return ebtsNR + 1 + i;
  89         }
  90     }
  91
  92     return NOTSET;
  93 }
  94
  95 #define FATAL() gmx_fatal(FARGS, "Reading Termini Database: not enough items on line\n%s", line)
  96
  97 static void read_atom(char *line, bool bAdd,
  98                       std::string *nname, t_atom *a, PreprocessingAtomTypes *atype, int *cgnr)
  99 {
 100     int    nr, i;
 101     char   buf[5][30];
 102     double m, q;
 103
 104     /* This code is messy, because of support for different formats:
 105      * for replace: [new name] <atom type> <m> <q> [cgnr (old format)]
 106      * for add:                <atom type> <m> <q> [cgnr]
 107      */
 108     nr = sscanf(line, "%s %s %s %s %s", buf[0], buf[1], buf[2], buf[3], buf[4]);
 109
 110     /* Here there an ambiguity due to the old replace format with cgnr,
 111      * which was read for years, but ignored in the rest of the code.
 112      * We have to assume that the atom type does not start with a digit
 113      * to make a line with 4 entries uniquely interpretable.
 114      */
 115     if (!bAdd && nr == 4 && isdigit(buf[1][0]))
 116     {
 117         nr = 3;
 118     }
 119
 120     if (nr < 3 || nr > 4)
 121     {
 122         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);
 123     }
 124     i = 0;
 125     if (!bAdd)
 126     {
 127         if (nr == 4)
 128         {
 129             *nname = buf[i++];
 130         }
 131         else
 132         {
 133             *nname = "";
 134         }
 135     }
 136     a->type = atype->atomTypeFromName(buf[i++]);
 137     sscanf(buf[i++], "%lf", &m);
 138     a->m = m;
 139     sscanf(buf[i++], "%lf", &q);
 140     a->q = q;
 141     if (bAdd && nr == 4)
 142     {
 143         sscanf(buf[i++], "%d", cgnr);
 144     }
 145     else
 146     {
 147         *cgnr = NOTSET;
 148     }
 149 }
 150
 151 static void print_atom(FILE *out, const t_atom &a, PreprocessingAtomTypes *atype)
 152 {
 153     fprintf(out, "\t%s\t%g\t%g\n",
 154             atype->atomNameFromAtomType(a.type), a.m, a.q);
 155 }
 156
 157 static void print_ter_db(const char *ff, char C, gmx::ArrayRef<const MoleculePatchDatabase> tb,
 158                          PreprocessingAtomTypes *atype)
 159 {
 160     std::string buf = gmx::formatString("%s-%c.tdb", ff, C);
 161     FILE       *out = gmx_fio_fopen(buf.c_str(), "w");
 162
 163     for (const auto &modification : tb)
 164     {
 165         fprintf(out, "[ %s ]\n", modification.name.c_str());
 166
 167         if (std::any_of(modification.hack.begin(), modification.hack.end(),
 168                         [](const auto &mod)
 169                         { return mod.type() == MoleculePatchType::Replace; }))
 170         {
 171             fprintf(out, "[ %s ]\n", kw_names[ekwRepl-ebtsNR-1]);
 172             for (const auto &hack : modification.hack)
 173             {
 174                 if (hack.type() == MoleculePatchType::Replace)
 175                 {
 176                     fprintf(out, "%s\t", hack.oname.c_str());
 177                     print_atom(out, hack.atom.back(), atype);
 178                 }
 179             }
 180         }
 181         if (std::any_of(modification.hack.begin(), modification.hack.end(),
 182                         [](const auto &mod)
 183                         { return mod.type() == MoleculePatchType::Add; }))
 184         {
 185             fprintf(out, "[ %s ]\n", kw_names[ekwAdd-ebtsNR-1]);
 186             for (const auto &hack : modification.hack)
 187             {
 188                 if (hack.type() == MoleculePatchType::Add)
 189                 {
 190                     print_ab(out, hack, hack.nname.c_str());
 191                     print_atom(out, hack.atom.back(), atype);
 192                 }
 193             }
 194         }
 195         if (std::any_of(modification.hack.begin(), modification.hack.end(),
 196                         [](const auto &mod)
 197                         { return mod.type() == MoleculePatchType::Delete; }))
 198         {
 199             fprintf(out, "[ %s ]\n", kw_names[ekwDel-ebtsNR-1]);
 200             for (const auto &hack : modification.hack)
 201             {
 202                 if (hack.type() == MoleculePatchType::Delete)
 203                 {
 204                     fprintf(out, "%s\n", hack.oname.c_str());
 205                 }
 206             }
 207         }
 208         for (int bt = 0; bt < ebtsNR; bt++)
 209         {
 210             if (!modification.rb[bt].b.empty())
 211             {
 212                 fprintf(out, "[ %s ]\n", btsNames[bt]);
 213                 for (const auto &b : modification.rb[bt].b)
 214                 {
 215                     for (int k = 0; k < btsNiatoms[bt]; k++)
 216                     {
 217                         fprintf(out, "%s%s", k ? "\t" : "", b.a[k].c_str());
 218                     }
 219                     if (!b.s.empty())
 220                     {
 221                         fprintf(out, "\t%s", b.s.c_str());
 222                     }
 223                     fprintf(out, "\n");
 224                 }
 225             }
 226         }
 227         fprintf(out, "\n");
 228     }
 229     gmx_fio_fclose(out);
 230 }
 231
 232 static void read_ter_db_file(const char                                         *fn,
 233                              std::vector<MoleculePatchDatabase>                 *tbptr,
 234                              PreprocessingAtomTypes                             *atype)
 235 {
 236     char         filebase[STRLEN];
 237     char         header[STRLEN], buf[STRLEN], line[STRLEN];
 238
 239     fflib_filename_base(fn, filebase, STRLEN);
 240     /* Remove the C/N termini extension */
 241     char *ptr = strrchr(filebase, '.');
 242     if (ptr != nullptr)
 243     {
 244         ptr[0] = '\0';
 245     }
 246
 247     FILE                                 *in = fflib_open(fn);
 248
 249     int                                   kwnr  = NOTSET;
 250     get_a_line(in, line, STRLEN);
 251     MoleculePatchDatabase                *block = nullptr;
 252     while (!feof(in))
 253     {
 254         if (get_header(line, header))
 255         {
 256             /* this is a new block, or a new keyword */
 257             kwnr = find_kw(header);
 258
 259             if (kwnr == NOTSET)
 260             {
 261                 tbptr->emplace_back(MoleculePatchDatabase());
 262                 block = &tbptr->back();
 263                 clearModificationBlock(block);
 264                 block->name     = header;
 265                 block->filebase = filebase;
 266             }
 267         }
 268         else
 269         {
 270             if (block == nullptr)
 271             {
 272                 gmx_fatal(FARGS, "reading termini database: "
 273                           "directive expected before line:\n%s\n"
 274                           "This might be a file in an old format.", line);
 275             }
 276             /* this is not a header, so it must be data */
 277             if (kwnr >= ebtsNR)
 278             {
 279                 /* this is a hack: add/rename/delete atoms */
 280                 /* make space for hacks */
 281                 block->hack.emplace_back(MoleculePatch());
 282                 MoleculePatch *hack = &block->hack.back();
 283
 284                 /* get data */
 285                 int n = 0;
 286                 if (kwnr == ekwRepl || kwnr == ekwDel)
 287                 {
 288                     if (sscanf(line, "%s%n", buf, &n) != 1)
 289                     {
 290                         gmx_fatal(FARGS, "Reading Termini Database '%s': "
 291                                   "expected atom name on line\n%s", fn, line);
 292                     }
 293                     hack->oname = buf;
 294                     /* we only replace or delete one atom at a time */
 295                     hack->nr = 1;
 296                 }
 297                 else if (kwnr == ekwAdd)
 298                 {
 299                     read_ab(line, fn, hack);
 300                     get_a_line(in, line, STRLEN);
 301                 }
 302                 else
 303                 {
 304                     gmx_fatal(FARGS, "unimplemented keyword number %d (%s:%d)",
 305                               kwnr, __FILE__, __LINE__);
 306                 }
 307                 if (kwnr == ekwRepl || kwnr == ekwAdd)
 308                 {
 309                     hack->atom.emplace_back();
 310                     read_atom(line+n, kwnr == ekwAdd,
 311                               &hack->nname, &hack->atom.back(), atype,
 312                               &hack->cgnr);
 313                     if (hack->nname.empty())
 314                     {
 315                         if (!hack->oname.empty())
 316                         {
 317                             hack->nname = hack->oname;
 318                         }
 319                         else
 320                         {
 321                             gmx_fatal(FARGS, "Reading Termini Database '%s': don't know which name the new atom should have on line\n%s", fn, line);
 322                         }
 323                     }
 324                 }
 325             }
 326             else if (kwnr >= 0 && kwnr < ebtsNR)
 327             {
 328                 /* this is bonded data: bonds, angles, dihedrals or impropers */
 329                 int                n = 0;
 330                 block->rb[kwnr].b.emplace_back();
 331                 BondedInteraction *newBond = &block->rb[kwnr].b.back();
 332                 for (int j = 0; j < btsNiatoms[kwnr]; j++)
 333                 {
 334                     int ni;
 335                     if (sscanf(line+n, "%s%n", buf, &ni) == 1)
 336                     {
 337                         newBond->a[j] = buf;
 338                     }
 339                     else
 340                     {
 341                         gmx_fatal(FARGS, "Reading Termini Database '%s': expected %d atom names (found %d) on line\n%s", fn, btsNiatoms[kwnr], j-1, line);
 342                     }
 343                     n += ni;
 344                 }
 345                 strcpy(buf, "");
 346                 sscanf(line+n, "%s", buf);
 347                 newBond->s = buf;
 348             }
 349             else
 350             {
 351                 gmx_fatal(FARGS, "Reading Termini Database: Expecting a header at line\n"
 352                           "%s", line);
 353             }
 354         }
 355         get_a_line(in, line, STRLEN);
 356     }
 357
 358     gmx_ffclose(in);
 359 }
 360
 361 int read_ter_db(const char *ffdir, char ter,
 362                 std::vector<MoleculePatchDatabase> *tbptr, PreprocessingAtomTypes *atype)
 363 {
 364     std::string ext = gmx::formatString(".%c.tdb", ter);
 365
 366     /* Search for termini database files.
 367      * Do not generate an error when none are found.
 368      */
 369     std::vector<std::string> tdbf  = fflib_search_file_end(ffdir, ext.c_str(), FALSE);
 370     tbptr->clear();
 371     for (const auto &filename : tdbf)
 372     {
 373         read_ter_db_file(filename.c_str(), tbptr, atype);
 374     }
 375
 376     if (debug)
 377     {
 378         print_ter_db("new", ter, *tbptr, atype);
 379     }
 380
 381     return tbptr->size();
 382 }
 383
 384 std::vector<MoleculePatchDatabase *>
 385 filter_ter(gmx::ArrayRef<MoleculePatchDatabase>                tb,
 386            const char                                         *resname)
 387 {
 388     // TODO Four years later, no force fields have ever used this, so decide status of this feature
 389     /* Since some force fields (e.g. OPLS) needs different
 390      * atomtypes for different residues there could be a lot
 391      * of entries in the databases for specific residues
 392      * (e.g. GLY-NH3+, SER-NH3+, ALA-NH3+).
 393      *
 394      * To reduce the database size, we assume that a terminus specifier liker
 395      *
 396      * [ GLY|SER|ALA-NH3+ ]
 397      *
 398      * would cover all of the three residue types above.
 399      * Wildcards (*,?) are not OK. Don't worry about e.g. GLU vs. GLUH since
 400      * pdb2gmx only uses the first 3 letters when calling this routine.
 401      *
 402      * To automate this, this routines scans a list of termini
 403      * for the residue name "resname" and returns an allocated list of
 404      * pointers to the termini that could be applied to the
 405      * residue in question. The variable pointed to by nret will
 406      * contain the number of valid pointers in the list.
 407      * Remember to free the list when you are done with it...
 408      */
 409
 410     auto none_idx = tb.end();
 411     std::vector<MoleculePatchDatabase *> list;
 412
 413     for (auto it = tb.begin(); it != tb.end(); it++)
 414     {
 415         const char *s     = it->name.c_str();
 416         bool        found = false;
 417         do
 418         {
 419             if (gmx_strncasecmp(resname, s, 3) == 0)
 420             {
 421                 found = true;
 422                 list.push_back(it);
 423             }
 424             else
 425             {
 426                 /* advance to next |-separated field */
 427                 s = strchr(s, '|');
 428                 if (s != nullptr)
 429                 {
 430                     s++;
 431                 }
 432             }
 433         }
 434         while (!found && s != nullptr);
 435     }
 436
 437     /* All residue-specific termini have been added. We might have to fall
 438      * back on generic termini, which are characterized by not having
 439      * '-' in the name prior to the last position (which indicates charge).
 440      * The [ None ] alternative is special since we don't want that
 441      * to be the default, so we put it last in the list we return.
 442      */
 443     for (auto it = tb.begin(); it != tb.end(); it++)
 444     {
 445         const char *s = it->name.c_str();
 446         if (gmx::equalCaseInsensitive("None", it->name))
 447         {
 448             none_idx = it;
 449         }
 450         else
 451         {
 452             /* Time to see if there's a generic terminus that matches.
 453                Is there a hyphen? */
 454             const char *c = strchr(s, '-');
 455
 456             /* A conjunction hyphen normally indicates a residue-specific
 457                terminus, which is named like "GLY-COOH". A generic terminus
 458                won't have a hyphen. */
 459             bool bFoundAnyHyphen = (c != nullptr);
 460             /* '-' as the last character indicates charge, so if that's
 461                the only one found e.g. "COO-", then it was not a conjunction
 462                hyphen, so this is a generic terminus */
 463             bool bOnlyFoundChargeHyphen = (bFoundAnyHyphen &&
 464                                            *(c+1) == '\0');
 465             /* Thus, "GLY-COO-" is not recognized as a generic terminus. */
 466             bool bFoundGenericTerminus = !bFoundAnyHyphen || bOnlyFoundChargeHyphen;
 467             if (bFoundGenericTerminus)
 468             {
 469                 /* Check that we haven't already added a residue-specific version
 470                  * of this terminus.
 471                  */
 472                 auto found = std::find_if(list.begin(), list.end(),
 473                                           [&s](const MoleculePatchDatabase *b)
 474                                           { return strstr(b->name.c_str(), s) != nullptr; });
 475                 if (found == list.end())
 476                 {
 477                     list.push_back(it);
 478                 }
 479             }
 480         }
 481     }
 482     if (none_idx != tb.end())
 483     {
 484         list.push_back(none_idx);
 485     }
 486
 487     return list;
 488 }
 489
 490
 491 MoleculePatchDatabase *choose_ter(gmx::ArrayRef<MoleculePatchDatabase *> tb, const char *title)
 492 {
 493     int sel, ret;
 494
 495     printf("%s\n", title);
 496     int i = 0;
 497     for (const auto &modification : tb)
 498     {
 499         bool bIsZwitterion = (0 == gmx_wcmatch("*ZWITTERION*", modification->name.c_str()));
 500         printf("%2d: %s%s\n", i, modification->name.c_str(),
 501                bIsZwitterion ? " (only use with zwitterions containing exactly one residue)" : "");
 502         i++;
 503     }
 504     do
 505     {
 506         ret = fscanf(stdin, "%d", &sel);
 507     }
 508     while ((ret != 1) || (sel < 0) || (sel >= tb.ssize()));
 509
 510     return tb[sel];
 511 }