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