src/gromacs/gmxpreprocess/genhydro.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-2019,2020, 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 "genhydro.h"
  40
  41 #include <cstring>
  42 #include <ctime>
  43
  44 #include "gromacs/fileio/confio.h"
  45 #include "gromacs/gmxlib/network.h"
  46 #include "gromacs/gmxpreprocess/calch.h"
  47 #include "gromacs/gmxpreprocess/h_db.h"
  48 #include "gromacs/gmxpreprocess/notset.h"
  49 #include "gromacs/gmxpreprocess/pgutil.h"
  50 #include "gromacs/gmxpreprocess/ter_db.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/topology/atoms.h"
  53 #include "gromacs/topology/symtab.h"
  54 #include "gromacs/utility/cstringutil.h"
  55 #include "gromacs/utility/exceptions.h"
  56 #include "gromacs/utility/fatalerror.h"
  57 #include "gromacs/utility/futil.h"
  58 #include "gromacs/utility/smalloc.h"
  59
  60 #include "hackblock.h"
  61 #include "resall.h"
  62
  63 static void copy_atom(const t_atoms* atoms1, int a1, t_atoms* atoms2, int a2, t_symtab* symtab)
  64 {
  65     atoms2->atom[a2]     = atoms1->atom[a1];
  66     atoms2->atomname[a2] = put_symtab(symtab, *atoms1->atomname[a1]);
  67 }
  68
  69 static int pdbasearch_atom(const char*              name,
  70                            int                      resind,
  71                            const t_atoms*           pdba,
  72                            const char*              searchtype,
  73                            bool                     bAllowMissing,
  74                            gmx::ArrayRef<const int> cyclicBondsIndex)
  75 {
  76     int i;
  77
  78     for (i = 0; (i < pdba->nr) && (pdba->atom[i].resind != resind); i++) {}
  79
  80     return search_atom(name, i, pdba, searchtype, bAllowMissing, cyclicBondsIndex);
  81 }
  82
  83 static void hacksearch_atom(int*                                            ii,
  84                             int*                                            jj,
  85                             const char*                                     name,
  86                             gmx::ArrayRef<const std::vector<MoleculePatch>> patches,
  87                             int                                             resind,
  88                             const t_atoms*                                  pdba)
  89 {
  90     int i;
  91
  92     *ii = -1;
  93     if (name[0] == '-')
  94     {
  95         name++;
  96         resind--;
  97     }
  98     for (i = 0; (i < pdba->nr) && (pdba->atom[i].resind != resind); i++) {}
  99     for (; (i < pdba->nr) && (pdba->atom[i].resind == resind) && (*ii < 0); i++)
 100     {
 101         int j = 0;
 102         for (const auto& patch : patches[i])
 103         {
 104             if (patch.nname == name)
 105             {
 106                 *ii = i;
 107                 *jj = j;
 108             }
 109             j++;
 110         }
 111     }
 112 }
 113
 114 static std::vector<MoleculePatchDatabase>
 115 getMoleculePatchDatabases(const t_atoms*                              pdba,
 116                           gmx::ArrayRef<const MoleculePatchDatabase>  globalPatches,
 117                           int                                         nterpairs,
 118                           gmx::ArrayRef<MoleculePatchDatabase* const> ntdb,
 119                           gmx::ArrayRef<MoleculePatchDatabase* const> ctdb,
 120                           gmx::ArrayRef<const int>                    rN,
 121                           gmx::ArrayRef<const int>                    rC)
 122 {
 123     std::vector<MoleculePatchDatabase> modBlock(pdba->nres);
 124     /* make space */
 125     /* first the termini */
 126     for (int i = 0; i < nterpairs; i++)
 127     {
 128         if (ntdb[i] != nullptr)
 129         {
 130             copyModificationBlocks(*ntdb[i], &modBlock[rN[i]]);
 131         }
 132         if (ctdb[i] != nullptr)
 133         {
 134             mergeAtomAndBondModifications(*ctdb[i], &modBlock[rC[i]]);
 135         }
 136     }
 137     /* then the whole hdb */
 138     for (int rnr = 0; rnr < pdba->nres; rnr++)
 139     {
 140         auto ahptr = search_h_db(globalPatches, *pdba->resinfo[rnr].rtp);
 141         if (ahptr != globalPatches.end())
 142         {
 143             if (modBlock[rnr].name.empty())
 144             {
 145                 modBlock[rnr].name = ahptr->name;
 146             }
 147             mergeAtomModifications(*ahptr, &modBlock[rnr]);
 148         }
 149     }
 150     return modBlock;
 151 }
 152
 153 static void expand_hackblocks_one(const MoleculePatchDatabase& newPatch,
 154                                   const std::string localAtomName, //NOLINT(performance-unnecessary-value-param)
 155                                   std::vector<MoleculePatch>* globalPatches,
 156                                   bool                        bN,
 157                                   bool                        bC)
 158 {
 159     /* we'll recursively add atoms to atoms */
 160     int pos = 0;
 161     for (auto& singlePatch : newPatch.hack)
 162     {
 163         /* first check if we're in the N- or C-terminus, then we should ignore
 164            all hacks involving atoms from resp. previous or next residue
 165            (i.e. which name begins with '-' (N) or '+' (C) */
 166         bool bIgnore = false;
 167         if (bN) /* N-terminus: ignore '-' */
 168         {
 169             for (int k = 0; k < 4 && !singlePatch.a[k].empty() && !bIgnore; k++)
 170             {
 171                 bIgnore = singlePatch.a[k][0] == '-';
 172             }
 173         }
 174         if (bC) /* C-terminus: ignore '+' */
 175         {
 176             for (int k = 0; k < 4 && !singlePatch.a[k].empty() && !bIgnore; k++)
 177             {
 178                 bIgnore = singlePatch.a[k][0] == '+';
 179             }
 180         }
 181         /* must be either hdb entry (tp>0) or add from tdb (oname==NULL)
 182            and first control aton (AI) matches this atom or
 183            delete/replace from tdb (oname!=NULL) and oname matches this atom */
 184
 185         if (!bIgnore
 186             && (((singlePatch.tp > 0 || singlePatch.oname.empty()) && singlePatch.a[0] == localAtomName)
 187                 || (singlePatch.oname == localAtomName)))
 188         {
 189             /* now expand all hacks for this atom */
 190             for (int k = 0; k < singlePatch.nr; k++)
 191             {
 192                 globalPatches->push_back(singlePatch);
 193                 MoleculePatch* patch = &globalPatches->back();
 194                 patch->bXSet         = false;
 195                 /* if we're adding (oname==NULL) and don't have a new name (nname)
 196                    yet, build it from localAtomName */
 197                 if (patch->nname.empty())
 198                 {
 199                     if (patch->oname.empty())
 200                     {
 201                         patch->nname    = localAtomName;
 202                         patch->nname[0] = 'H';
 203                     }
 204                 }
 205                 else
 206                 {
 207                     if (gmx_debug_at)
 208                     {
 209                         fprintf(debug,
 210                                 "Hack '%s' %d, replacing nname '%s' with '%s' (old name '%s')\n",
 211                                 localAtomName.c_str(), pos, patch->nname.c_str(),
 212                                 singlePatch.nname.c_str(),
 213                                 patch->oname.empty() ? "" : patch->oname.c_str());
 214                     }
 215                     patch->nname = singlePatch.nname;
 216                 }
 217
 218                 if (singlePatch.tp == 10 && k == 2)
 219                 {
 220                     /* This is a water virtual site, not a hydrogen */
 221                     /* Ugly hardcoded name hack */
 222                     patch->nname.assign("M");
 223                 }
 224                 else if (singlePatch.tp == 11 && k >= 2)
 225                 {
 226                     /* This is a water lone pair, not a hydrogen */
 227                     /* Ugly hardcoded name hack */
 228                     patch->nname.assign(gmx::formatString("LP%d", 1 + k - 2));
 229                 }
 230                 else if (singlePatch.nr > 1)
 231                 {
 232                     /* adding more than one atom, number them */
 233                     patch->nname.append(gmx::formatString("%d", 1 + k));
 234                 }
 235             }
 236
 237             /* add hacks to atoms we've just added */
 238             if (singlePatch.tp > 0 || singlePatch.oname.empty())
 239             {
 240                 for (int k = 0; k < singlePatch.nr; k++)
 241                 {
 242                     expand_hackblocks_one(
 243                             newPatch, globalPatches->at(globalPatches->size() - singlePatch.nr + k).nname,
 244                             globalPatches, bN, bC);
 245                 }
 246             }
 247         }
 248         pos++;
 249     }
 250 }
 251
 252 static void expand_hackblocks(const t_atoms*                             pdba,
 253                               gmx::ArrayRef<const MoleculePatchDatabase> hb,
 254                               gmx::ArrayRef<std::vector<MoleculePatch>>  patches,
 255                               int                                        nterpairs,
 256                               gmx::ArrayRef<const int>                   rN,
 257                               gmx::ArrayRef<const int>                   rC)
 258 {
 259     for (int i = 0; i < pdba->nr; i++)
 260     {
 261         bool bN = false;
 262         for (int j = 0; j < nterpairs && !bN; j++)
 263         {
 264             bN = pdba->atom[i].resind == rN[j];
 265         }
 266         bool bC = false;
 267         for (int j = 0; j < nterpairs && !bC; j++)
 268         {
 269             bC = pdba->atom[i].resind == rC[j];
 270         }
 271
 272         /* add hacks to this atom */
 273         expand_hackblocks_one(hb[pdba->atom[i].resind], *pdba->atomname[i], &patches[i], bN, bC);
 274     }
 275 }
 276
 277 static int check_atoms_present(const t_atoms*                            pdba,
 278                                gmx::ArrayRef<std::vector<MoleculePatch>> patches,
 279                                gmx::ArrayRef<const int>                  cyclicBondsIndex)
 280 {
 281     int nadd = 0;
 282     for (int i = 0; i < pdba->nr; i++)
 283     {
 284         int rnr = pdba->atom[i].resind;
 285         for (auto patch = patches[i].begin(); patch != patches[i].end(); patch++)
 286         {
 287             switch (patch->type())
 288             {
 289                 case MoleculePatchType::Add:
 290                 {
 291                     /* we're adding */
 292                     /* check if the atom is already present */
 293                     int k = pdbasearch_atom(patch->nname.c_str(), rnr, pdba, "check", TRUE, cyclicBondsIndex);
 294                     if (k != -1)
 295                     {
 296                         /* We found the added atom. */
 297                         patch->bAlreadyPresent = true;
 298                     }
 299                     else
 300                     {
 301                         patch->bAlreadyPresent = false;
 302                         /* count how many atoms we'll add */
 303                         nadd++;
 304                     }
 305                     break;
 306                 }
 307                 case MoleculePatchType::Delete:
 308                 {
 309                     /* we're deleting */
 310                     nadd--;
 311                     break;
 312                 }
 313                 case MoleculePatchType::Replace: { break;
 314                 }
 315                 default: { GMX_THROW(gmx::InternalError("Case not handled"));
 316                 }
 317             }
 318         }
 319     }
 320     return nadd;
 321 }
 322
 323 static void calc_all_pos(const t_atoms*                            pdba,
 324                          gmx::ArrayRef<const gmx::RVec>            x,
 325                          gmx::ArrayRef<std::vector<MoleculePatch>> patches,
 326                          bool                                      bCheckMissing,
 327                          gmx::ArrayRef<const int>                  cyclicBondsIndex)
 328 {
 329     int ii, l = 0;
 330 #define MAXH 4
 331     rvec xa[4];    /* control atoms for calc_h_pos */
 332     rvec xh[MAXH]; /* hydrogen positions from calc_h_pos */
 333
 334     int jj = 0;
 335
 336     for (int i = 0; i < pdba->nr; i++)
 337     {
 338         int rnr = pdba->atom[i].resind;
 339         for (auto patch = patches[i].begin(); patch != patches[i].end(); patch += patch->nr)
 340         {
 341             GMX_RELEASE_ASSERT(patch < patches[i].end(),
 342                                "The number of patches in the last patch can not exceed the total "
 343                                "number of patches");
 344             /* check if we're adding: */
 345             if (patch->type() == MoleculePatchType::Add && patch->tp > 0)
 346             {
 347                 bool bFoundAll = true;
 348                 for (int m = 0; (m < patch->nctl && bFoundAll); m++)
 349                 {
 350                     int ia = pdbasearch_atom(patch->a[m].c_str(), rnr, pdba, bCheckMissing ? "atom" : "check",
 351                                              !bCheckMissing, cyclicBondsIndex);
 352                     if (ia < 0)
 353                     {
 354                         /* not found in original atoms, might still be in
 355                          * the patch Instructions (patches) */
 356                         hacksearch_atom(&ii, &jj, patch->a[m].c_str(), patches, rnr, pdba);
 357                         if (ii >= 0)
 358                         {
 359                             copy_rvec(patches[ii][jj].newx, xa[m]);
 360                         }
 361                         else
 362                         {
 363                             bFoundAll = false;
 364                             if (bCheckMissing)
 365                             {
 366                                 gmx_fatal(FARGS,
 367                                           "Atom %s not found in residue %s %d"
 368                                           ", rtp entry %s"
 369                                           " while adding hydrogens",
 370                                           patch->a[m].c_str(), *pdba->resinfo[rnr].name,
 371                                           pdba->resinfo[rnr].nr, *pdba->resinfo[rnr].rtp);
 372                             }
 373                         }
 374                     }
 375                     else
 376                     {
 377                         copy_rvec(x[ia], xa[m]);
 378                     }
 379                 }
 380                 if (bFoundAll)
 381                 {
 382                     for (int m = 0; (m < MAXH); m++)
 383                     {
 384                         for (int d = 0; d < DIM; d++)
 385                         {
 386                             if (m < patch->nr)
 387                             {
 388                                 xh[m][d] = 0;
 389                             }
 390                             else
 391                             {
 392                                 xh[m][d] = NOTSET;
 393                             }
 394                         }
 395                     }
 396                     calc_h_pos(patch->tp, xa, xh, &l);
 397                     for (int m = 0; m < patch->nr; m++)
 398                     {
 399                         auto next = patch + m;
 400                         copy_rvec(xh[m], next->newx);
 401                         next->bXSet = true;
 402                     }
 403                 }
 404             }
 405         }
 406     }
 407 }
 408
 409 static int add_h_low(t_atoms**                                   initialAtoms,
 410                      t_atoms**                                   modifiedAtoms,
 411                      std::vector<gmx::RVec>*                     xptr,
 412                      gmx::ArrayRef<const MoleculePatchDatabase>  globalPatches,
 413                      t_symtab*                                   symtab,
 414                      const int                                   nterpairs,
 415                      gmx::ArrayRef<MoleculePatchDatabase* const> ntdb,
 416                      gmx::ArrayRef<MoleculePatchDatabase* const> ctdb,
 417                      gmx::ArrayRef<const int>                    rN,
 418                      gmx::ArrayRef<const int>                    rC,
 419                      const bool                                  bCheckMissing,
 420                      gmx::ArrayRef<const int>                    cyclicBondsIndex)
 421 {
 422     int                                     nadd;
 423     int                                     newi, natoms, nalreadypresent;
 424     std::vector<std::vector<MoleculePatch>> patches;
 425     std::vector<gmx::RVec>                  xn;
 426
 427     t_atoms* pdba = *initialAtoms;
 428
 429     /* set flags for adding hydrogens (according to hdb) */
 430     natoms = pdba->nr;
 431
 432     {
 433         /* We'll have to do all the hard work */
 434         /* first get all the hackblocks for each residue: */
 435         std::vector<MoleculePatchDatabase> hb =
 436                 getMoleculePatchDatabases(pdba, globalPatches, nterpairs, ntdb, ctdb, rN, rC);
 437
 438         /* expand the hackblocks to atom level */
 439         patches.resize(natoms);
 440         expand_hackblocks(pdba, hb, patches, nterpairs, rN, rC);
 441     }
 442
 443     /* Now calc the positions */
 444     calc_all_pos(pdba, *xptr, patches, bCheckMissing, cyclicBondsIndex);
 445
 446     /* we don't have to add atoms that are already present in initialAtoms,
 447        so we will remove them from the patches (MoleculePatch) */
 448     nadd = check_atoms_present(pdba, patches, cyclicBondsIndex);
 449
 450     /* Copy old atoms, making space for new ones */
 451     if (nadd > 0)
 452     {
 453         srenew(*modifiedAtoms, 1);
 454         init_t_atoms(*modifiedAtoms, natoms + nadd, FALSE);
 455         (*modifiedAtoms)->nres = pdba->nres;
 456         srenew((*modifiedAtoms)->resinfo, pdba->nres);
 457         std::copy(pdba->resinfo, pdba->resinfo + pdba->nres, (*modifiedAtoms)->resinfo);
 458     }
 459     if (nadd == 0)
 460     {
 461         return natoms;
 462     }
 463
 464     xn.resize(natoms + nadd);
 465     newi = 0;
 466     for (int i = 0; (i < natoms); i++)
 467     {
 468         /* check if this atom wasn't scheduled for deletion */
 469         if (patches[i].empty() || (!patches[i][0].nname.empty()))
 470         {
 471             if (newi >= natoms + nadd)
 472             {
 473                 /*gmx_fatal(FARGS,"Not enough space for adding atoms");*/
 474                 nadd += 10;
 475                 xn.resize(natoms + nadd);
 476                 srenew((*modifiedAtoms)->atom, natoms + nadd);
 477                 srenew((*modifiedAtoms)->atomname, natoms + nadd);
 478             }
 479             copy_atom(pdba, i, (*modifiedAtoms), newi, symtab);
 480             copy_rvec((*xptr)[i], xn[newi]);
 481             /* process the hacks for this atom */
 482             nalreadypresent = 0;
 483             for (auto patch = patches[i].begin(); patch != patches[i].end(); patch++)
 484             {
 485                 if (patch->type() == MoleculePatchType::Add) /* add */
 486                 {
 487                     newi++;
 488                     if (newi >= natoms + nadd)
 489                     {
 490                         /* gmx_fatal(FARGS,"Not enough space for adding atoms");*/
 491                         nadd += 10;
 492                         xn.resize(natoms + nadd);
 493                         srenew((*modifiedAtoms)->atom, natoms + nadd);
 494                         srenew((*modifiedAtoms)->atomname, natoms + nadd);
 495                     }
 496                     (*modifiedAtoms)->atom[newi].resind = pdba->atom[i].resind;
 497                 }
 498                 if (!patch->nname.empty()
 499                     && (patch->oname.empty() || patch->oname == *(*modifiedAtoms)->atomname[newi]))
 500                 {
 501                     /* add or replace */
 502                     if (patch->type() == MoleculePatchType::Add && patch->bAlreadyPresent)
 503                     {
 504                         /* This atom is already present, copy it from the input. */
 505                         nalreadypresent++;
 506                         copy_atom(pdba, i + nalreadypresent, (*modifiedAtoms), newi, symtab);
 507                         copy_rvec((*xptr)[i + nalreadypresent], xn[newi]);
 508                     }
 509                     else
 510                     {
 511                         if (gmx_debug_at)
 512                         {
 513                             fprintf(debug, "Replacing %d '%s' with (old name '%s') %s\n", newi,
 514                                     ((*modifiedAtoms)->atomname[newi] && *(*modifiedAtoms)->atomname[newi])
 515                                             ? *(*modifiedAtoms)->atomname[newi]
 516                                             : "",
 517                                     patch->oname.empty() ? "" : patch->oname.c_str(),
 518                                     patch->nname.c_str());
 519                         }
 520                         (*modifiedAtoms)->atomname[newi] = put_symtab(symtab, patch->nname.c_str());
 521                         if (patch->bXSet)
 522                         {
 523                             copy_rvec(patch->newx, xn[newi]);
 524                         }
 525                     }
 526                     if (debug)
 527                     {
 528                         fprintf(debug, " %s %g %g", *(*modifiedAtoms)->atomname[newi],
 529                                 (*modifiedAtoms)->atom[newi].m, (*modifiedAtoms)->atom[newi].q);
 530                     }
 531                 }
 532             }
 533             newi++;
 534             i += nalreadypresent;
 535         }
 536     }
 537     (*modifiedAtoms)->nr = newi;
 538
 539     done_atom(pdba);
 540     *initialAtoms = *modifiedAtoms;
 541
 542     *xptr = xn;
 543
 544     return newi;
 545 }
 546
 547 int add_h(t_atoms**                                   initialAtoms,
 548           t_atoms**                                   localAtoms,
 549           std::vector<gmx::RVec>*                     xptr,
 550           gmx::ArrayRef<const MoleculePatchDatabase>  globalPatches,
 551           t_symtab*                                   symtab,
 552           const int                                   nterpairs,
 553           gmx::ArrayRef<MoleculePatchDatabase* const> ntdb,
 554           gmx::ArrayRef<MoleculePatchDatabase* const> ctdb,
 555           gmx::ArrayRef<const int>                    rN,
 556           gmx::ArrayRef<const int>                    rC,
 557           const bool                                  bAllowMissing,
 558           gmx::ArrayRef<const int>                    cyclicBondsIndex)
 559 {
 560     int nold, nnew, niter;
 561
 562     /* Here we loop to be able to add atoms to added atoms.
 563      * We should not check for missing atoms here.
 564      */
 565     niter = 0;
 566     nnew  = 0;
 567     do
 568     {
 569         nold = nnew;
 570         nnew = add_h_low(initialAtoms, localAtoms, xptr, globalPatches, symtab, nterpairs, ntdb,
 571                          ctdb, rN, rC, FALSE, cyclicBondsIndex);
 572         niter++;
 573         if (niter > 100)
 574         {
 575             gmx_fatal(FARGS,
 576                       "More than 100 iterations of add_h. Maybe you are trying to replace an added "
 577                       "atom (this is not supported)?");
 578         }
 579     } while (nnew > nold);
 580
 581     if (!bAllowMissing)
 582     {
 583         /* Call add_h_low once more, now only for the missing atoms check */
 584         add_h_low(initialAtoms, localAtoms, xptr, globalPatches, symtab, nterpairs, ntdb, ctdb, rN,
 585                   rC, TRUE, cyclicBondsIndex);
 586     }
 587
 588     return nnew;
 589 }