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,2014,2015,2017,2018, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "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/resall.h"
  51 #include "gromacs/gmxpreprocess/ter_db.h"
  52 #include "gromacs/math/vec.h"
  53 #include "gromacs/topology/symtab.h"
  54 #include "gromacs/utility/cstringutil.h"
  55 #include "gromacs/utility/fatalerror.h"
  56 #include "gromacs/utility/futil.h"
  57 #include "gromacs/utility/smalloc.h"
  58
  59 static void copy_atom(t_atoms *atoms1, int a1, t_atoms *atoms2, int a2)
  60 {
  61     atoms2->atom[a2] = atoms1->atom[a1];
  62     snew(atoms2->atomname[a2], 1);
  63     *atoms2->atomname[a2] = gmx_strdup(*atoms1->atomname[a1]);
  64 }
  65
  66 static int pdbasearch_atom(const char *name, int resind, t_atoms *pdba,
  67                            const char *searchtype, bool bAllowMissing)
  68 {
  69     int  i;
  70
  71     for (i = 0; (i < pdba->nr) && (pdba->atom[i].resind != resind); i++)
  72     {
  73         ;
  74     }
  75
  76     return search_atom(name, i, pdba,
  77                        searchtype, bAllowMissing);
  78 }
  79
  80 static void hacksearch_atom(int *ii, int *jj, char *name,
  81                             const int nab[], t_hack *ab[],
  82                             int resind, t_atoms *pdba)
  83 {
  84     int  i, j;
  85
  86     *ii = -1;
  87     if (name[0] == '-')
  88     {
  89         name++;
  90         resind--;
  91     }
  92     for (i = 0; (i < pdba->nr) && (pdba->atom[i].resind != resind); i++)
  93     {
  94         ;
  95     }
  96     for (; (i < pdba->nr) && (pdba->atom[i].resind == resind) && (*ii < 0); i++)
  97     {
  98         for (j = 0; (j < nab[i]) && (*ii < 0); j++)
  99         {
 100             if (ab[i][j].nname && strcmp(name, ab[i][j].nname) == 0)
 101             {
 102                 *ii = i;
 103                 *jj = j;
 104             }
 105         }
 106     }
 107
 108 }
 109
 110 static t_hackblock *get_hackblocks(t_atoms *pdba, int nah, t_hackblock ah[],
 111                                    int nterpairs,
 112                                    t_hackblock **ntdb, t_hackblock **ctdb,
 113                                    const int *rN, const int *rC)
 114 {
 115     int          i, rnr;
 116     t_hackblock *hb, *ahptr;
 117
 118     /* make space */
 119     snew(hb, pdba->nres);
 120     /* first the termini */
 121     for (i = 0; i < nterpairs; i++)
 122     {
 123         if (ntdb[i] != nullptr)
 124         {
 125             copy_t_hackblock(ntdb[i], &hb[rN[i]]);
 126         }
 127         if (ctdb[i] != nullptr)
 128         {
 129             merge_t_hackblock(ctdb[i], &hb[rC[i]]);
 130         }
 131     }
 132     /* then the whole hdb */
 133     for (rnr = 0; rnr < pdba->nres; rnr++)
 134     {
 135         ahptr = search_h_db(nah, ah, *pdba->resinfo[rnr].rtp);
 136         if (ahptr)
 137         {
 138             if (hb[rnr].name == nullptr)
 139             {
 140                 hb[rnr].name = gmx_strdup(ahptr->name);
 141             }
 142             merge_hacks(ahptr, &hb[rnr]);
 143         }
 144     }
 145     return hb;
 146 }
 147
 148 static void expand_hackblocks_one(t_hackblock *hbr, char *atomname,
 149                                   int *nabi, t_hack **abi, bool bN, bool bC)
 150 {
 151     int      j, k, l;
 152     bool     bIgnore;
 153
 154     /* we'll recursively add atoms to atoms */
 155     for (j = 0; j < hbr->nhack; j++)
 156     {
 157         /* first check if we're in the N- or C-terminus, then we should ignore
 158            all hacks involving atoms from resp. previous or next residue
 159            (i.e. which name begins with '-' (N) or '+' (C) */
 160         bIgnore = FALSE;
 161         if (bN) /* N-terminus: ignore '-' */
 162         {
 163             for (k = 0; k < 4 && hbr->hack[j].a[k] && !bIgnore; k++)
 164             {
 165                 bIgnore = hbr->hack[j].a[k][0] == '-';
 166             }
 167         }
 168         if (bC) /* C-terminus: ignore '+' */
 169         {
 170             for (k = 0; k < 4 && hbr->hack[j].a[k] && !bIgnore; k++)
 171             {
 172                 bIgnore = hbr->hack[j].a[k][0] == '+';
 173             }
 174         }
 175         /* must be either hdb entry (tp>0) or add from tdb (oname==NULL)
 176            and first control aton (AI) matches this atom or
 177            delete/replace from tdb (oname!=NULL) and oname matches this atom */
 178
 179         if (!bIgnore &&
 180             ( ( ( hbr->hack[j].tp > 0 || hbr->hack[j].oname == nullptr ) &&
 181                 strcmp(atomname, hbr->hack[j].ai()) == 0 ) ||
 182               ( hbr->hack[j].oname != nullptr &&
 183                 strcmp(atomname, hbr->hack[j].oname) == 0) ) )
 184         {
 185             /* now expand all hacks for this atom */
 186             srenew(*abi, *nabi + hbr->hack[j].nr);
 187             for (k = 0; k < hbr->hack[j].nr; k++)
 188             {
 189                 copy_t_hack(&hbr->hack[j], &(*abi)[*nabi + k]);
 190                 (*abi)[*nabi + k].bXSet = FALSE;
 191                 /* if we're adding (oname==NULL) and don't have a new name (nname)
 192                    yet, build it from atomname */
 193                 if ( (*abi)[*nabi + k].nname == nullptr)
 194                 {
 195                     if ( (*abi)[*nabi + k].oname == nullptr)
 196                     {
 197                         (*abi)[*nabi + k].nname    = gmx_strdup(atomname);
 198                         (*abi)[*nabi + k].nname[0] = 'H';
 199                     }
 200                 }
 201                 else
 202                 {
 203                     if (gmx_debug_at)
 204                     {
 205                         fprintf(debug, "Hack '%s' %d, replacing nname '%s' with '%s' (old name '%s')\n",
 206                                 atomname, j,
 207                                 (*abi)[*nabi + k].nname, hbr->hack[j].nname,
 208                                 (*abi)[*nabi + k].oname ? (*abi)[*nabi + k].oname : "");
 209                     }
 210                     sfree((*abi)[*nabi + k].nname);
 211                     (*abi)[*nabi + k].nname = gmx_strdup(hbr->hack[j].nname);
 212                 }
 213
 214                 if (hbr->hack[j].tp == 10 && k == 2)
 215                 {
 216                     /* This is a water virtual site, not a hydrogen */
 217                     /* Ugly hardcoded name hack */
 218                     (*abi)[*nabi + k].nname[0] = 'M';
 219                 }
 220                 else if (hbr->hack[j].tp == 11 && k >= 2)
 221                 {
 222                     /* This is a water lone pair, not a hydrogen */
 223                     /* Ugly hardcoded name hack */
 224                     srenew((*abi)[*nabi + k].nname, 4);
 225                     (*abi)[*nabi + k].nname[0] = 'L';
 226                     (*abi)[*nabi + k].nname[1] = 'P';
 227                     (*abi)[*nabi + k].nname[2] = '1' + k - 2;
 228                     (*abi)[*nabi + k].nname[3] = '\0';
 229                 }
 230                 else if (hbr->hack[j].nr > 1)
 231                 {
 232                     /* adding more than one atom, number them */
 233                     l = strlen((*abi)[*nabi + k].nname);
 234                     srenew((*abi)[*nabi + k].nname, l+2);
 235                     (*abi)[*nabi + k].nname[l]   = '1' + k;
 236                     (*abi)[*nabi + k].nname[l+1] = '\0';
 237                 }
 238             }
 239             (*nabi) += hbr->hack[j].nr;
 240
 241             /* add hacks to atoms we've just added */
 242             if (hbr->hack[j].tp > 0 || hbr->hack[j].oname == nullptr)
 243             {
 244                 for (k = 0; k < hbr->hack[j].nr; k++)
 245                 {
 246                     expand_hackblocks_one(hbr, (*abi)[*nabi-hbr->hack[j].nr+k].nname,
 247                                           nabi, abi, bN, bC);
 248                 }
 249             }
 250         }
 251     }
 252 }
 253
 254 static void expand_hackblocks(t_atoms *pdba, t_hackblock hb[],
 255                               int nab[], t_hack *ab[],
 256                               int nterpairs, const int *rN, const int *rC)
 257 {
 258     int      i, j;
 259     bool     bN, bC;
 260
 261     for (i = 0; i < pdba->nr; i++)
 262     {
 263         bN = FALSE;
 264         for (j = 0; j < nterpairs && !bN; j++)
 265         {
 266             bN = pdba->atom[i].resind == rN[j];
 267         }
 268         bC = FALSE;
 269         for (j = 0; j < nterpairs && !bC; j++)
 270         {
 271             bC = pdba->atom[i].resind == rC[j];
 272         }
 273
 274         /* add hacks to this atom */
 275         expand_hackblocks_one(&hb[pdba->atom[i].resind], *pdba->atomname[i],
 276                               &nab[i], &ab[i], bN, bC);
 277     }
 278 }
 279
 280 static int check_atoms_present(t_atoms *pdba, const int nab[], t_hack *ab[])
 281 {
 282     int i, j, k, rnr, nadd;
 283
 284     nadd = 0;
 285     for (i = 0; i < pdba->nr; i++)
 286     {
 287         rnr = pdba->atom[i].resind;
 288         for (j = 0; j < nab[i]; j++)
 289         {
 290             if (ab[i][j].oname == nullptr)
 291             {
 292                 /* we're adding */
 293                 if (ab[i][j].nname == nullptr)
 294                 {
 295                     gmx_incons("ab[i][j].nname not allocated");
 296                 }
 297                 /* check if the atom is already present */
 298                 k = pdbasearch_atom(ab[i][j].nname, rnr, pdba, "check", TRUE);
 299                 if (k != -1)
 300                 {
 301                     /* We found the added atom. */
 302                     ab[i][j].bAlreadyPresent = TRUE;
 303                 }
 304                 else
 305                 {
 306                     ab[i][j].bAlreadyPresent = FALSE;
 307                     /* count how many atoms we'll add */
 308                     nadd++;
 309                 }
 310             }
 311             else if (ab[i][j].nname == nullptr)
 312             {
 313                 /* we're deleting */
 314                 nadd--;
 315             }
 316         }
 317     }
 318
 319     return nadd;
 320 }
 321
 322 static void calc_all_pos(t_atoms *pdba, rvec x[], int nab[], t_hack *ab[],
 323                          bool bCheckMissing)
 324 {
 325     int      i, j, ii, jj, m, ia, d, rnr, l = 0;
 326 #define MAXH 4
 327     rvec     xa[4];    /* control atoms for calc_h_pos */
 328     rvec     xh[MAXH]; /* hydrogen positions from calc_h_pos */
 329     bool     bFoundAll;
 330
 331     jj = 0;
 332
 333     for (i = 0; i < pdba->nr; i++)
 334     {
 335         rnr   = pdba->atom[i].resind;
 336         for (j = 0; j < nab[i]; j += ab[i][j].nr)
 337         {
 338             /* check if we're adding: */
 339             if (ab[i][j].oname == nullptr && ab[i][j].tp > 0)
 340             {
 341                 bFoundAll = TRUE;
 342                 for (m = 0; (m < ab[i][j].nctl && bFoundAll); m++)
 343                 {
 344                     ia = pdbasearch_atom(ab[i][j].a[m], rnr, pdba,
 345                                          bCheckMissing ? "atom" : "check",
 346                                          !bCheckMissing);
 347                     if (ia < 0)
 348                     {
 349                         /* not found in original atoms, might still be in t_hack (ab) */
 350                         hacksearch_atom(&ii, &jj, ab[i][j].a[m], nab, ab, rnr, pdba);
 351                         if (ii >= 0)
 352                         {
 353                             copy_rvec(ab[ii][jj].newx, xa[m]);
 354                         }
 355                         else
 356                         {
 357                             bFoundAll = FALSE;
 358                             if (bCheckMissing)
 359                             {
 360                                 gmx_fatal(FARGS, "Atom %s not found in residue %s %d"
 361                                           ", rtp entry %s"
 362                                           " while adding hydrogens",
 363                                           ab[i][j].a[m],
 364                                           *pdba->resinfo[rnr].name,
 365                                           pdba->resinfo[rnr].nr,
 366                                           *pdba->resinfo[rnr].rtp);
 367                             }
 368                         }
 369                     }
 370                     else
 371                     {
 372                         copy_rvec(x[ia], xa[m]);
 373                     }
 374                 }
 375                 if (bFoundAll)
 376                 {
 377                     for (m = 0; (m < MAXH); m++)
 378                     {
 379                         for (d = 0; d < DIM; d++)
 380                         {
 381                             if (m < ab[i][j].nr)
 382                             {
 383                                 xh[m][d] = 0;
 384                             }
 385                             else
 386                             {
 387                                 xh[m][d] = NOTSET;
 388                             }
 389                         }
 390                     }
 391                     calc_h_pos(ab[i][j].tp, xa, xh, &l);
 392                     for (m = 0; m < ab[i][j].nr; m++)
 393                     {
 394                         copy_rvec(xh[m], ab[i][j+m].newx);
 395                         ab[i][j+m].bXSet = TRUE;
 396                     }
 397                 }
 398             }
 399         }
 400     }
 401 }
 402
 403 static void free_ab(int natoms, int *nab, t_hack **ab)
 404 {
 405     int i;
 406
 407     for (i = 0; i < natoms; i++)
 408     {
 409         free_t_hack(nab[i], &ab[i]);
 410     }
 411     sfree(nab);
 412     sfree(ab);
 413 }
 414
 415 static int add_h_low(t_atoms **pdbaptr, rvec *xptr[],
 416                      int nah, t_hackblock ah[],
 417                      int nterpairs, t_hackblock **ntdb, t_hackblock **ctdb,
 418                      int *rN, int *rC, bool bCheckMissing,
 419                      int **nabptr, t_hack ***abptr,
 420                      bool bUpdate_pdba, bool bKeep_old_pdba)
 421 {
 422     t_atoms        *newpdba = nullptr, *pdba = nullptr;
 423     int             nadd;
 424     int             i, newi, j, natoms, nalreadypresent;
 425     int            *nab = nullptr;
 426     t_hack        **ab  = nullptr;
 427     t_hackblock    *hb;
 428     rvec           *xn;
 429     bool            bKeep_ab;
 430
 431     /* set flags for adding hydrogens (according to hdb) */
 432     pdba   = *pdbaptr;
 433     natoms = pdba->nr;
 434
 435     if (nabptr && abptr)
 436     {
 437         /* the first time these will be pointers to NULL, but we will
 438            return in them the completed arrays, which we will get back
 439            the second time */
 440         nab      = *nabptr;
 441         ab       = *abptr;
 442         bKeep_ab = TRUE;
 443     }
 444     else
 445     {
 446         bKeep_ab = FALSE;
 447     }
 448
 449     if (nab && ab)
 450     {
 451         /* WOW, everything was already figured out */
 452         bUpdate_pdba = FALSE;
 453     }
 454     else
 455     {
 456         /* We'll have to do all the hard work */
 457         bUpdate_pdba = TRUE;
 458         /* first get all the hackblocks for each residue: */
 459         hb = get_hackblocks(pdba, nah, ah, nterpairs, ntdb, ctdb, rN, rC);
 460
 461         /* expand the hackblocks to atom level */
 462         snew(nab, natoms);
 463         snew(ab, natoms);
 464         expand_hackblocks(pdba, hb, nab, ab, nterpairs, rN, rC);
 465         free_t_hackblock(pdba->nres, &hb);
 466     }
 467
 468     /* Now calc the positions */
 469     calc_all_pos(pdba, *xptr, nab, ab, bCheckMissing);
 470
 471     if (bUpdate_pdba)
 472     {
 473         /* we don't have to add atoms that are already present in pdba,
 474            so we will remove them from the ab (t_hack) */
 475         nadd = check_atoms_present(pdba, nab, ab);
 476
 477         /* Copy old atoms, making space for new ones */
 478         snew(newpdba, 1);
 479         init_t_atoms(newpdba, natoms+nadd, FALSE);
 480         newpdba->nres    = pdba->nres;
 481         sfree(newpdba->resinfo);
 482         newpdba->resinfo = pdba->resinfo;
 483     }
 484     else
 485     {
 486         nadd = 0;
 487     }
 488
 489     if (nadd == 0)
 490     {
 491         /* There is nothing to do: return now */
 492         if (!bKeep_ab)
 493         {
 494             free_ab(natoms, nab, ab);
 495         }
 496
 497         return natoms;
 498     }
 499
 500     snew(xn, natoms+nadd);
 501     newi = 0;
 502     for (i = 0; (i < natoms); i++)
 503     {
 504         /* check if this atom wasn't scheduled for deletion */
 505         if (nab[i] == 0 || (ab[i][0].nname != nullptr) )
 506         {
 507             if (newi >= natoms+nadd)
 508             {
 509                 /*gmx_fatal(FARGS,"Not enough space for adding atoms");*/
 510                 nadd += 10;
 511                 srenew(xn, natoms+nadd);
 512                 if (bUpdate_pdba)
 513                 {
 514                     srenew(newpdba->atom, natoms+nadd);
 515                     srenew(newpdba->atomname, natoms+nadd);
 516                 }
 517             }
 518             if (bUpdate_pdba)
 519             {
 520                 copy_atom(pdba, i, newpdba, newi);
 521             }
 522             copy_rvec((*xptr)[i], xn[newi]);
 523             /* process the hacks for this atom */
 524             nalreadypresent = 0;
 525             for (j = 0; j < nab[i]; j++)
 526             {
 527                 if (ab[i][j].oname == nullptr) /* add */
 528                 {
 529                     newi++;
 530                     if (newi >= natoms+nadd)
 531                     {
 532                         /* gmx_fatal(FARGS,"Not enough space for adding atoms");*/
 533                         nadd += 10;
 534                         srenew(xn, natoms+nadd);
 535                         if (bUpdate_pdba)
 536                         {
 537                             srenew(newpdba->atom, natoms+nadd);
 538                             srenew(newpdba->atomname, natoms+nadd);
 539                         }
 540                     }
 541                     if (bUpdate_pdba)
 542                     {
 543                         newpdba->atom[newi].resind = pdba->atom[i].resind;
 544                     }
 545                 }
 546                 if (ab[i][j].nname != nullptr &&
 547                     (ab[i][j].oname == nullptr ||
 548                      strcmp(ab[i][j].oname, *newpdba->atomname[newi]) == 0))
 549                 {
 550                     /* add or replace */
 551                     if (ab[i][j].oname == nullptr && ab[i][j].bAlreadyPresent)
 552                     {
 553                         /* This atom is already present, copy it from the input. */
 554                         nalreadypresent++;
 555                         if (bUpdate_pdba)
 556                         {
 557                             copy_atom(pdba, i+nalreadypresent, newpdba, newi);
 558                         }
 559                         copy_rvec((*xptr)[i+nalreadypresent], xn[newi]);
 560                     }
 561                     else
 562                     {
 563                         if (bUpdate_pdba)
 564                         {
 565                             if (gmx_debug_at)
 566                             {
 567                                 fprintf(debug, "Replacing %d '%s' with (old name '%s') %s\n",
 568                                         newi,
 569                                         (newpdba->atomname[newi] && *newpdba->atomname[newi]) ? *newpdba->atomname[newi] : "",
 570                                         ab[i][j].oname ? ab[i][j].oname : "",
 571                                         ab[i][j].nname);
 572                             }
 573                             snew(newpdba->atomname[newi], 1);
 574                             *newpdba->atomname[newi] = gmx_strdup(ab[i][j].nname);
 575                             if (ab[i][j].oname != nullptr && ab[i][j].atom) /* replace */
 576                             {                                               /*          newpdba->atom[newi].m    = ab[i][j].atom->m; */
 577 /*        newpdba->atom[newi].q    = ab[i][j].atom->q; */
 578 /*        newpdba->atom[newi].type = ab[i][j].atom->type; */
 579                             }
 580                         }
 581                         if (ab[i][j].bXSet)
 582                         {
 583                             copy_rvec(ab[i][j].newx, xn[newi]);
 584                         }
 585                     }
 586                     if (bUpdate_pdba && debug)
 587                     {
 588                         fprintf(debug, " %s %g %g", *newpdba->atomname[newi],
 589                                 newpdba->atom[newi].m, newpdba->atom[newi].q);
 590                     }
 591                 }
 592             }
 593             newi++;
 594             i += nalreadypresent;
 595         }
 596     }
 597     if (bUpdate_pdba)
 598     {
 599         newpdba->nr = newi;
 600     }
 601
 602     if (bKeep_ab)
 603     {
 604         *nabptr = nab;
 605         *abptr  = ab;
 606     }
 607     else
 608     {
 609         /* Clean up */
 610         free_ab(natoms, nab, ab);
 611     }
 612
 613     if (bUpdate_pdba)
 614     {
 615         if (!bKeep_old_pdba)
 616         {
 617             for (i = 0; i < natoms; i++)
 618             {
 619                 /* Do not free the atomname string itself, it might be in symtab */
 620                 /* sfree(*(pdba->atomname[i])); */
 621                 /* sfree(pdba->atomname[i]); */
 622             }
 623             sfree(pdba->atomname);
 624             sfree(pdba->atom);
 625             sfree(pdba->pdbinfo);
 626             sfree(pdba);
 627         }
 628         *pdbaptr = newpdba;
 629     }
 630
 631     sfree(*xptr);
 632     *xptr = xn;
 633
 634     return newi;
 635 }
 636
 637 int add_h(t_atoms **pdbaptr, rvec *xptr[],
 638           int nah, t_hackblock ah[],
 639           int nterpairs, t_hackblock **ntdb, t_hackblock **ctdb,
 640           int *rN, int *rC, bool bAllowMissing,
 641           int **nabptr, t_hack ***abptr,
 642           bool bUpdate_pdba, bool bKeep_old_pdba)
 643 {
 644     int nold, nnew, niter;
 645
 646     /* Here we loop to be able to add atoms to added atoms.
 647      * We should not check for missing atoms here.
 648      */
 649     niter = 0;
 650     nnew  = 0;
 651     do
 652     {
 653         nold = nnew;
 654         nnew = add_h_low(pdbaptr, xptr, nah, ah, nterpairs, ntdb, ctdb, rN, rC, FALSE,
 655                          nabptr, abptr, bUpdate_pdba, bKeep_old_pdba);
 656         niter++;
 657         if (niter > 100)
 658         {
 659             gmx_fatal(FARGS, "More than 100 iterations of add_h. Maybe you are trying to replace an added atom (this is not supported)?");
 660         }
 661     }
 662     while (nnew > nold);
 663
 664     if (!bAllowMissing)
 665     {
 666         /* Call add_h_low once more, now only for the missing atoms check */
 667         add_h_low(pdbaptr, xptr, nah, ah, nterpairs, ntdb, ctdb, rN, rC, TRUE,
 668                   nabptr, abptr, bUpdate_pdba, bKeep_old_pdba);
 669     }
 670
 671     return nnew;
 672 }