src/gromacs/gmxpreprocess/genhydro.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "config.h"
  40
  41 #include <string.h>
  42 #include <time.h>
  43
  44 #include "gromacs/legacyheaders/typedefs.h"
  45 #include "gromacs/fileio/confio.h"
  46 #include "gromacs/math/vec.h"
  47 #include "gromacs/utility/futil.h"
  48 #include "calch.h"
  49 #include "genhydro.h"
  50 #include "h_db.h"
  51 #include "ter_db.h"
  52 #include "resall.h"
  53 #include "pgutil.h"
  54 #include "gromacs/legacyheaders/network.h"
  55
  56 #include "gromacs/topology/symtab.h"
  57 #include "gromacs/utility/cstringutil.h"
  58 #include "gromacs/utility/fatalerror.h"
  59 #include "gromacs/utility/smalloc.h"
  60
  61 static void copy_atom(t_atoms *atoms1, int a1, t_atoms *atoms2, int a2)
  62 {
  63     atoms2->atom[a2] = atoms1->atom[a1];
  64     snew(atoms2->atomname[a2], 1);
  65     *atoms2->atomname[a2] = gmx_strdup(*atoms1->atomname[a1]);
  66 }
  67
  68 static atom_id pdbasearch_atom(const char *name, int resind, t_atoms *pdba,
  69                                const char *searchtype, gmx_bool bAllowMissing)
  70 {
  71     int  i;
  72
  73     for (i = 0; (i < pdba->nr) && (pdba->atom[i].resind != resind); i++)
  74     {
  75         ;
  76     }
  77
  78     return search_atom(name, i, pdba,
  79                        searchtype, bAllowMissing);
  80 }
  81
  82 static void hacksearch_atom(int *ii, int *jj, char *name,
  83                             int nab[], t_hack *ab[],
  84                             int resind, t_atoms *pdba)
  85 {
  86     int  i, j;
  87
  88     *ii = -1;
  89     if (name[0] == '-')
  90     {
  91         name++;
  92         resind--;
  93     }
  94     for (i = 0; (i < pdba->nr) && (pdba->atom[i].resind != resind); i++)
  95     {
  96         ;
  97     }
  98     for (; (i < pdba->nr) && (pdba->atom[i].resind == resind) && (*ii < 0); i++)
  99     {
 100         for (j = 0; (j < nab[i]) && (*ii < 0); j++)
 101         {
 102             if (ab[i][j].nname && strcmp(name, ab[i][j].nname) == 0)
 103             {
 104                 *ii = i;
 105                 *jj = j;
 106             }
 107         }
 108     }
 109
 110     return;
 111 }
 112
 113 void dump_ab(FILE *out, int natom, int nab[], t_hack *ab[], gmx_bool bHeader)
 114 {
 115     int i, j;
 116
 117 #define SS(s) (s) ? (s) : "-"
 118     /* dump ab */
 119     if (bHeader)
 120     {
 121         fprintf(out, "ADDBLOCK (t_hack) natom=%d\n"
 122                 "%4s %2s %-4s %-4s %2s %-4s %-4s %-4s %-4s %1s %s\n",
 123                 natom, "atom", "nr", "old", "new", "tp", "ai", "aj", "ak", "al", "a", "x");
 124     }
 125     for (i = 0; i < natom; i++)
 126     {
 127         for (j = 0; j < nab[i]; j++)
 128         {
 129             fprintf(out, "%4d %2d %-4s %-4s %2d %-4s %-4s %-4s %-4s %s %g %g %g\n",
 130                     i+1, ab[i][j].nr, SS(ab[i][j].oname), SS(ab[i][j].nname),
 131                     ab[i][j].tp,
 132                     SS(ab[i][j].AI), SS(ab[i][j].AJ),
 133                     SS(ab[i][j].AK), SS(ab[i][j].AL),
 134                     ab[i][j].atom ? "+" : "",
 135                     ab[i][j].newx[XX], ab[i][j].newx[YY], ab[i][j].newx[ZZ]);
 136         }
 137     }
 138 #undef SS
 139 }
 140
 141 static t_hackblock *get_hackblocks(t_atoms *pdba, int nah, t_hackblock ah[],
 142                                    int nterpairs,
 143                                    t_hackblock **ntdb, t_hackblock **ctdb,
 144                                    int *rN, int *rC)
 145 {
 146     int          i, rnr;
 147     t_hackblock *hb, *ahptr;
 148
 149     /* make space */
 150     snew(hb, pdba->nres);
 151     /* first the termini */
 152     for (i = 0; i < nterpairs; i++)
 153     {
 154         if (ntdb[i] != NULL)
 155         {
 156             copy_t_hackblock(ntdb[i], &hb[rN[i]]);
 157         }
 158         if (ctdb[i] != NULL)
 159         {
 160             merge_t_hackblock(ctdb[i], &hb[rC[i]]);
 161         }
 162     }
 163     /* then the whole hdb */
 164     for (rnr = 0; rnr < pdba->nres; rnr++)
 165     {
 166         ahptr = search_h_db(nah, ah, *pdba->resinfo[rnr].rtp);
 167         if (ahptr)
 168         {
 169             if (hb[rnr].name == NULL)
 170             {
 171                 hb[rnr].name = gmx_strdup(ahptr->name);
 172             }
 173             merge_hacks(ahptr, &hb[rnr]);
 174         }
 175     }
 176     return hb;
 177 }
 178
 179 static void expand_hackblocks_one(t_hackblock *hbr, char *atomname,
 180                                   int *nabi, t_hack **abi, gmx_bool bN, gmx_bool bC)
 181 {
 182     int      j, k, l, d;
 183     gmx_bool bIgnore;
 184
 185     /* we'll recursively add atoms to atoms */
 186     for (j = 0; j < hbr->nhack; j++)
 187     {
 188         /* first check if we're in the N- or C-terminus, then we should ignore
 189            all hacks involving atoms from resp. previous or next residue
 190            (i.e. which name begins with '-' (N) or '+' (C) */
 191         bIgnore = FALSE;
 192         if (bN) /* N-terminus: ignore '-' */
 193         {
 194             for (k = 0; k < 4 && hbr->hack[j].a[k] && !bIgnore; k++)
 195             {
 196                 bIgnore = hbr->hack[j].a[k][0] == '-';
 197             }
 198         }
 199         if (bC) /* C-terminus: ignore '+' */
 200         {
 201             for (k = 0; k < 4 && hbr->hack[j].a[k] && !bIgnore; k++)
 202             {
 203                 bIgnore = hbr->hack[j].a[k][0] == '+';
 204             }
 205         }
 206         /* must be either hdb entry (tp>0) or add from tdb (oname==NULL)
 207            and first control aton (AI) matches this atom or
 208            delete/replace from tdb (oname!=NULL) and oname matches this atom */
 209         if (debug)
 210         {
 211             fprintf(debug, " %s", hbr->hack[j].oname ? hbr->hack[j].oname : hbr->hack[j].AI);
 212         }
 213
 214         if (!bIgnore &&
 215             ( ( ( hbr->hack[j].tp > 0 || hbr->hack[j].oname == NULL ) &&
 216                 strcmp(atomname, hbr->hack[j].AI) == 0 ) ||
 217               ( hbr->hack[j].oname != NULL &&
 218                 strcmp(atomname, hbr->hack[j].oname) == 0) ) )
 219         {
 220             /* now expand all hacks for this atom */
 221             if (debug)
 222             {
 223                 fprintf(debug, " +%dh", hbr->hack[j].nr);
 224             }
 225             srenew(*abi, *nabi + hbr->hack[j].nr);
 226             for (k = 0; k < hbr->hack[j].nr; k++)
 227             {
 228                 copy_t_hack(&hbr->hack[j], &(*abi)[*nabi + k]);
 229                 (*abi)[*nabi + k].bXSet = FALSE;
 230                 /* if we're adding (oname==NULL) and don't have a new name (nname)
 231                    yet, build it from atomname */
 232                 if ( (*abi)[*nabi + k].nname == NULL)
 233                 {
 234                     if ( (*abi)[*nabi + k].oname == NULL)
 235                     {
 236                         (*abi)[*nabi + k].nname    = gmx_strdup(atomname);
 237                         (*abi)[*nabi + k].nname[0] = 'H';
 238                     }
 239                 }
 240                 else
 241                 {
 242                     if (gmx_debug_at)
 243                     {
 244                         fprintf(debug, "Hack '%s' %d, replacing nname '%s' with '%s' (old name '%s')\n",
 245                                 atomname, j,
 246                                 (*abi)[*nabi + k].nname, hbr->hack[j].nname,
 247                                 (*abi)[*nabi + k].oname ? (*abi)[*nabi + k].oname : "");
 248                     }
 249                     sfree((*abi)[*nabi + k].nname);
 250                     (*abi)[*nabi + k].nname = gmx_strdup(hbr->hack[j].nname);
 251                 }
 252
 253                 if (hbr->hack[j].tp == 10 && k == 2)
 254                 {
 255                     /* This is a water virtual site, not a hydrogen */
 256                     /* Ugly hardcoded name hack */
 257                     (*abi)[*nabi + k].nname[0] = 'M';
 258                 }
 259                 else if (hbr->hack[j].tp == 11 && k >= 2)
 260                 {
 261                     /* This is a water lone pair, not a hydrogen */
 262                     /* Ugly hardcoded name hack */
 263                     srenew((*abi)[*nabi + k].nname, 4);
 264                     (*abi)[*nabi + k].nname[0] = 'L';
 265                     (*abi)[*nabi + k].nname[1] = 'P';
 266                     (*abi)[*nabi + k].nname[2] = '1' + k - 2;
 267                     (*abi)[*nabi + k].nname[3] = '\0';
 268                 }
 269                 else if (hbr->hack[j].nr > 1)
 270                 {
 271                     /* adding more than one atom, number them */
 272                     l = strlen((*abi)[*nabi + k].nname);
 273                     srenew((*abi)[*nabi + k].nname, l+2);
 274                     (*abi)[*nabi + k].nname[l]   = '1' + k;
 275                     (*abi)[*nabi + k].nname[l+1] = '\0';
 276                 }
 277             }
 278             (*nabi) += hbr->hack[j].nr;
 279
 280             /* add hacks to atoms we've just added */
 281             if (hbr->hack[j].tp > 0 || hbr->hack[j].oname == NULL)
 282             {
 283                 for (k = 0; k < hbr->hack[j].nr; k++)
 284                 {
 285                     expand_hackblocks_one(hbr, (*abi)[*nabi-hbr->hack[j].nr+k].nname,
 286                                           nabi, abi, bN, bC);
 287                 }
 288             }
 289         }
 290     }
 291 }
 292
 293 static void expand_hackblocks(t_atoms *pdba, t_hackblock hb[],
 294                               int nab[], t_hack *ab[],
 295                               int nterpairs, int *rN, int *rC)
 296 {
 297     int      i, j;
 298     gmx_bool bN, bC;
 299
 300     for (i = 0; i < pdba->nr; i++)
 301     {
 302         bN = FALSE;
 303         for (j = 0; j < nterpairs && !bN; j++)
 304         {
 305             bN = pdba->atom[i].resind == rN[j];
 306         }
 307         bC = FALSE;
 308         for (j = 0; j < nterpairs && !bC; j++)
 309         {
 310             bC = pdba->atom[i].resind == rC[j];
 311         }
 312
 313         /* add hacks to this atom */
 314         expand_hackblocks_one(&hb[pdba->atom[i].resind], *pdba->atomname[i],
 315                               &nab[i], &ab[i], bN, bC);
 316     }
 317     if (debug)
 318     {
 319         fprintf(debug, "\n");
 320     }
 321 }
 322
 323 static int check_atoms_present(t_atoms *pdba, int nab[], t_hack *ab[])
 324 {
 325     int i, j, k, d, rnr, nadd;
 326
 327     nadd = 0;
 328     for (i = 0; i < pdba->nr; i++)
 329     {
 330         rnr = pdba->atom[i].resind;
 331         for (j = 0; j < nab[i]; j++)
 332         {
 333             if (ab[i][j].oname == NULL)
 334             {
 335                 /* we're adding */
 336                 if (ab[i][j].nname == NULL)
 337                 {
 338                     gmx_incons("ab[i][j].nname not allocated");
 339                 }
 340                 /* check if the atom is already present */
 341                 k = pdbasearch_atom(ab[i][j].nname, rnr, pdba, "check", TRUE);
 342                 if (k != -1)
 343                 {
 344                     /* We found the added atom. */
 345                     ab[i][j].bAlreadyPresent = TRUE;
 346                     if (debug)
 347                     {
 348                         fprintf(debug, "Atom '%s' in residue '%s' %d is already present\n",
 349                                 ab[i][j].nname,
 350                                 *pdba->resinfo[rnr].name, pdba->resinfo[rnr].nr);
 351                     }
 352                 }
 353                 else
 354                 {
 355                     ab[i][j].bAlreadyPresent = FALSE;
 356                     /* count how many atoms we'll add */
 357                     nadd++;
 358                 }
 359             }
 360             else if (ab[i][j].nname == NULL)
 361             {
 362                 /* we're deleting */
 363                 nadd--;
 364             }
 365         }
 366     }
 367
 368     return nadd;
 369 }
 370
 371 static void calc_all_pos(t_atoms *pdba, rvec x[], int nab[], t_hack *ab[],
 372                          gmx_bool bCheckMissing)
 373 {
 374     int      i, j, ii, jj, m, ia, d, rnr, l = 0;
 375 #define MAXH 4
 376     rvec     xa[4];    /* control atoms for calc_h_pos */
 377     rvec     xh[MAXH]; /* hydrogen positions from calc_h_pos */
 378     gmx_bool bFoundAll;
 379
 380     jj = 0;
 381
 382     for (i = 0; i < pdba->nr; i++)
 383     {
 384         rnr   = pdba->atom[i].resind;
 385         for (j = 0; j < nab[i]; j += ab[i][j].nr)
 386         {
 387             /* check if we're adding: */
 388             if (ab[i][j].oname == NULL && ab[i][j].tp > 0)
 389             {
 390                 bFoundAll = TRUE;
 391                 for (m = 0; (m < ab[i][j].nctl && bFoundAll); m++)
 392                 {
 393                     ia = pdbasearch_atom(ab[i][j].a[m], rnr, pdba,
 394                                          bCheckMissing ? "atom" : "check",
 395                                          !bCheckMissing);
 396                     if (ia < 0)
 397                     {
 398                         /* not found in original atoms, might still be in t_hack (ab) */
 399                         hacksearch_atom(&ii, &jj, ab[i][j].a[m], nab, ab, rnr, pdba);
 400                         if (ii >= 0)
 401                         {
 402                             copy_rvec(ab[ii][jj].newx, xa[m]);
 403                         }
 404                         else
 405                         {
 406                             bFoundAll = FALSE;
 407                             if (bCheckMissing)
 408                             {
 409                                 gmx_fatal(FARGS, "Atom %s not found in residue %s %d"
 410                                           ", rtp entry %s"
 411                                           " while adding hydrogens",
 412                                           ab[i][j].a[m],
 413                                           *pdba->resinfo[rnr].name,
 414                                           pdba->resinfo[rnr].nr,
 415                                           *pdba->resinfo[rnr].rtp);
 416                             }
 417                         }
 418                     }
 419                     else
 420                     {
 421                         copy_rvec(x[ia], xa[m]);
 422                     }
 423                 }
 424                 if (bFoundAll)
 425                 {
 426                     for (m = 0; (m < MAXH); m++)
 427                     {
 428                         for (d = 0; d < DIM; d++)
 429                         {
 430                             if (m < ab[i][j].nr)
 431                             {
 432                                 xh[m][d] = 0;
 433                             }
 434                             else
 435                             {
 436                                 xh[m][d] = NOTSET;
 437                             }
 438                         }
 439                     }
 440                     calc_h_pos(ab[i][j].tp, xa, xh, &l);
 441                     for (m = 0; m < ab[i][j].nr; m++)
 442                     {
 443                         copy_rvec(xh[m], ab[i][j+m].newx);
 444                         ab[i][j+m].bXSet = TRUE;
 445                     }
 446                 }
 447             }
 448         }
 449     }
 450 }
 451
 452 static void free_ab(int natoms, int *nab, t_hack **ab)
 453 {
 454     int i;
 455
 456     for (i = 0; i < natoms; i++)
 457     {
 458         free_t_hack(nab[i], &ab[i]);
 459     }
 460     sfree(nab);
 461     sfree(ab);
 462 }
 463
 464 static int add_h_low(t_atoms **pdbaptr, rvec *xptr[],
 465                      int nah, t_hackblock ah[],
 466                      int nterpairs, t_hackblock **ntdb, t_hackblock **ctdb,
 467                      int *rN, int *rC, gmx_bool bCheckMissing,
 468                      int **nabptr, t_hack ***abptr,
 469                      gmx_bool bUpdate_pdba, gmx_bool bKeep_old_pdba)
 470 {
 471     t_atoms        *newpdba = NULL, *pdba = NULL;
 472     int             nadd;
 473     int             i, newi, j, d, natoms, nalreadypresent;
 474     int            *nab = NULL;
 475     t_hack        **ab  = NULL;
 476     t_hackblock    *hb;
 477     rvec           *xn;
 478     gmx_bool        bKeep_ab;
 479
 480     /* set flags for adding hydrogens (according to hdb) */
 481     pdba   = *pdbaptr;
 482     natoms = pdba->nr;
 483
 484     if (nabptr && abptr)
 485     {
 486         /* the first time these will be pointers to NULL, but we will
 487            return in them the completed arrays, which we will get back
 488            the second time */
 489         nab      = *nabptr;
 490         ab       = *abptr;
 491         bKeep_ab = TRUE;
 492         if (debug)
 493         {
 494             fprintf(debug, "pointer to ab found\n");
 495         }
 496     }
 497     else
 498     {
 499         bKeep_ab = FALSE;
 500     }
 501
 502     if (nab && ab)
 503     {
 504         /* WOW, everything was already figured out */
 505         bUpdate_pdba = FALSE;
 506         if (debug)
 507         {
 508             fprintf(debug, "pointer to non-null ab found\n");
 509         }
 510     }
 511     else
 512     {
 513         /* We'll have to do all the hard work */
 514         bUpdate_pdba = TRUE;
 515         /* first get all the hackblocks for each residue: */
 516         hb = get_hackblocks(pdba, nah, ah, nterpairs, ntdb, ctdb, rN, rC);
 517         if (debug)
 518         {
 519             dump_hb(debug, pdba->nres, hb);
 520         }
 521
 522         /* expand the hackblocks to atom level */
 523         snew(nab, natoms);
 524         snew(ab, natoms);
 525         expand_hackblocks(pdba, hb, nab, ab, nterpairs, rN, rC);
 526         free_t_hackblock(pdba->nres, &hb);
 527     }
 528
 529     if (debug)
 530     {
 531         fprintf(debug, "before calc_all_pos\n");
 532         dump_ab(debug, natoms, nab, ab, TRUE);
 533     }
 534
 535     /* Now calc the positions */
 536     calc_all_pos(pdba, *xptr, nab, ab, bCheckMissing);
 537
 538     if (debug)
 539     {
 540         fprintf(debug, "after calc_all_pos\n");
 541         dump_ab(debug, natoms, nab, ab, TRUE);
 542     }
 543
 544     if (bUpdate_pdba)
 545     {
 546         /* we don't have to add atoms that are already present in pdba,
 547            so we will remove them from the ab (t_hack) */
 548         nadd = check_atoms_present(pdba, nab, ab);
 549         if (debug)
 550         {
 551             fprintf(debug, "removed add hacks that were already in pdba:\n");
 552             dump_ab(debug, natoms, nab, ab, TRUE);
 553             fprintf(debug, "will be adding %d atoms\n", nadd);
 554         }
 555
 556         /* Copy old atoms, making space for new ones */
 557         snew(newpdba, 1);
 558         init_t_atoms(newpdba, natoms+nadd, FALSE);
 559         newpdba->nres    = pdba->nres;
 560         sfree(newpdba->resinfo);
 561         newpdba->resinfo = pdba->resinfo;
 562     }
 563     else
 564     {
 565         nadd = 0;
 566     }
 567     if (debug)
 568     {
 569         fprintf(debug, "snew xn for %d old + %d new atoms %d total)\n",
 570                 natoms, nadd, natoms+nadd);
 571     }
 572
 573     if (nadd == 0)
 574     {
 575         /* There is nothing to do: return now */
 576         if (!bKeep_ab)
 577         {
 578             free_ab(natoms, nab, ab);
 579         }
 580
 581         return natoms;
 582     }
 583
 584     snew(xn, natoms+nadd);
 585     newi = 0;
 586     for (i = 0; (i < natoms); i++)
 587     {
 588         /* check if this atom wasn't scheduled for deletion */
 589         if (nab[i] == 0 || (ab[i][0].nname != NULL) )
 590         {
 591             if (newi >= natoms+nadd)
 592             {
 593                 /*gmx_fatal(FARGS,"Not enough space for adding atoms");*/
 594                 nadd += 10;
 595                 srenew(xn, natoms+nadd);
 596                 if (bUpdate_pdba)
 597                 {
 598                     srenew(newpdba->atom, natoms+nadd);
 599                     srenew(newpdba->atomname, natoms+nadd);
 600                 }
 601                 debug_gmx();
 602             }
 603             if (debug)
 604             {
 605                 fprintf(debug, "(%3d) %3d %4s %4s%3d %3d",
 606                         i+1, newi+1, *pdba->atomname[i],
 607                         *pdba->resinfo[pdba->atom[i].resind].name,
 608                         pdba->resinfo[pdba->atom[i].resind].nr, nab[i]);
 609             }
 610             if (bUpdate_pdba)
 611             {
 612                 copy_atom(pdba, i, newpdba, newi);
 613             }
 614             copy_rvec((*xptr)[i], xn[newi]);
 615             /* process the hacks for this atom */
 616             nalreadypresent = 0;
 617             for (j = 0; j < nab[i]; j++)
 618             {
 619                 if (ab[i][j].oname == NULL) /* add */
 620                 {
 621                     newi++;
 622                     if (newi >= natoms+nadd)
 623                     {
 624                         /* gmx_fatal(FARGS,"Not enough space for adding atoms");*/
 625                         nadd += 10;
 626                         srenew(xn, natoms+nadd);
 627                         if (bUpdate_pdba)
 628                         {
 629                             srenew(newpdba->atom, natoms+nadd);
 630                             srenew(newpdba->atomname, natoms+nadd);
 631                         }
 632                         debug_gmx();
 633                     }
 634                     if (bUpdate_pdba)
 635                     {
 636                         newpdba->atom[newi].resind = pdba->atom[i].resind;
 637                     }
 638                     if (debug)
 639                     {
 640                         fprintf(debug, " + %d", newi+1);
 641                     }
 642                 }
 643                 if (ab[i][j].nname != NULL &&
 644                     (ab[i][j].oname == NULL ||
 645                      strcmp(ab[i][j].oname, *newpdba->atomname[newi]) == 0))
 646                 {
 647                     /* add or replace */
 648                     if (ab[i][j].oname == NULL && ab[i][j].bAlreadyPresent)
 649                     {
 650                         /* This atom is already present, copy it from the input. */
 651                         nalreadypresent++;
 652                         if (bUpdate_pdba)
 653                         {
 654                             copy_atom(pdba, i+nalreadypresent, newpdba, newi);
 655                         }
 656                         copy_rvec((*xptr)[i+nalreadypresent], xn[newi]);
 657                     }
 658                     else
 659                     {
 660                         if (bUpdate_pdba)
 661                         {
 662                             if (gmx_debug_at)
 663                             {
 664                                 fprintf(debug, "Replacing %d '%s' with (old name '%s') %s\n",
 665                                         newi,
 666                                         (newpdba->atomname[newi] && *newpdba->atomname[newi]) ? *newpdba->atomname[newi] : "",
 667                                         ab[i][j].oname ? ab[i][j].oname : "",
 668                                         ab[i][j].nname);
 669                             }
 670                             snew(newpdba->atomname[newi], 1);
 671                             *newpdba->atomname[newi] = gmx_strdup(ab[i][j].nname);
 672                             if (ab[i][j].oname != NULL && ab[i][j].atom) /* replace */
 673                             {                                            /*          newpdba->atom[newi].m    = ab[i][j].atom->m; */
 674 /*        newpdba->atom[newi].q    = ab[i][j].atom->q; */
 675 /*        newpdba->atom[newi].type = ab[i][j].atom->type; */
 676                             }
 677                         }
 678                         if (ab[i][j].bXSet)
 679                         {
 680                             copy_rvec(ab[i][j].newx, xn[newi]);
 681                         }
 682                     }
 683                     if (bUpdate_pdba && debug)
 684                     {
 685                         fprintf(debug, " %s %g %g", *newpdba->atomname[newi],
 686                                 newpdba->atom[newi].m, newpdba->atom[newi].q);
 687                     }
 688                 }
 689             }
 690             newi++;
 691             i += nalreadypresent;
 692             if (debug)
 693             {
 694                 fprintf(debug, "\n");
 695             }
 696         }
 697     }
 698     if (bUpdate_pdba)
 699     {
 700         newpdba->nr = newi;
 701     }
 702
 703     if (bKeep_ab)
 704     {
 705         *nabptr = nab;
 706         *abptr  = ab;
 707     }
 708     else
 709     {
 710         /* Clean up */
 711         free_ab(natoms, nab, ab);
 712     }
 713
 714     if (bUpdate_pdba)
 715     {
 716         if (!bKeep_old_pdba)
 717         {
 718             for (i = 0; i < natoms; i++)
 719             {
 720                 /* Do not free the atomname string itself, it might be in symtab */
 721                 /* sfree(*(pdba->atomname[i])); */
 722                 /* sfree(pdba->atomname[i]); */
 723             }
 724             sfree(pdba->atomname);
 725             sfree(pdba->atom);
 726             sfree(pdba->pdbinfo);
 727             sfree(pdba);
 728         }
 729         *pdbaptr = newpdba;
 730     }
 731     else
 732     {
 733         nadd = newi-natoms;
 734     }
 735
 736     sfree(*xptr);
 737     *xptr = xn;
 738
 739     return newi;
 740 }
 741
 742 void deprotonate(t_atoms *atoms, rvec *x)
 743 {
 744     int  i, j;
 745
 746     j = 0;
 747     for (i = 0; i < atoms->nr; i++)
 748     {
 749         if ( (*atoms->atomname[i])[0] != 'H')
 750         {
 751             atoms->atomname[j] = atoms->atomname[i];
 752             atoms->atom[j]     = atoms->atom[i];
 753             copy_rvec(x[i], x[j]);
 754             j++;
 755         }
 756     }
 757     atoms->nr = j;
 758 }
 759
 760 int add_h(t_atoms **pdbaptr, rvec *xptr[],
 761           int nah, t_hackblock ah[],
 762           int nterpairs, t_hackblock **ntdb, t_hackblock **ctdb,
 763           int *rN, int *rC, gmx_bool bAllowMissing,
 764           int **nabptr, t_hack ***abptr,
 765           gmx_bool bUpdate_pdba, gmx_bool bKeep_old_pdba)
 766 {
 767     int nold, nnew, niter;
 768
 769     /* Here we loop to be able to add atoms to added atoms.
 770      * We should not check for missing atoms here.
 771      */
 772     niter = 0;
 773     nnew  = 0;
 774     do
 775     {
 776         nold = nnew;
 777         nnew = add_h_low(pdbaptr, xptr, nah, ah, nterpairs, ntdb, ctdb, rN, rC, FALSE,
 778                          nabptr, abptr, bUpdate_pdba, bKeep_old_pdba);
 779         niter++;
 780         if (niter > 100)
 781         {
 782             gmx_fatal(FARGS, "More than 100 iterations of add_h. Maybe you are trying to replace an added atom (this is not supported)?");
 783         }
 784     }
 785     while (nnew > nold);
 786
 787     if (!bAllowMissing)
 788     {
 789         /* Call add_h_low once more, now only for the missing atoms check */
 790         add_h_low(pdbaptr, xptr, nah, ah, nterpairs, ntdb, ctdb, rN, rC, TRUE,
 791                   nabptr, abptr, bUpdate_pdba, bKeep_old_pdba);
 792     }
 793
 794     return nnew;
 795 }
 796
 797 int protonate(t_atoms **atomsptr, rvec **xptr, t_protonate *protdata)
 798 {
 799 #define NTERPAIRS 1
 800     t_atoms    *atoms;
 801     gmx_bool    bUpdate_pdba, bKeep_old_pdba;
 802     int         nntdb, nctdb, nt, ct;
 803     int         nadd;
 804
 805     atoms = NULL;
 806     if (!protdata->bInit)
 807     {
 808         if (debug)
 809         {
 810             fprintf(debug, "protonate: Initializing protdata\n");
 811         }
 812
 813         /* set forcefield to use: */
 814         strcpy(protdata->FF, "oplsaa.ff");
 815
 816         /* get the databases: */
 817         protdata->nah = read_h_db(protdata->FF, &protdata->ah);
 818         open_symtab(&protdata->tab);
 819         protdata->atype = read_atype(protdata->FF, &protdata->tab);
 820         nntdb           = read_ter_db(protdata->FF, 'n', &protdata->ntdb, protdata->atype);
 821         if (nntdb < 1)
 822         {
 823             gmx_fatal(FARGS, "no N-terminus database");
 824         }
 825         nctdb = read_ter_db(protdata->FF, 'c', &protdata->ctdb, protdata->atype);
 826         if (nctdb < 1)
 827         {
 828             gmx_fatal(FARGS, "no C-terminus database");
 829         }
 830
 831         /* set terminus types: -NH3+ (different for Proline) and -COO- */
 832         atoms = *atomsptr;
 833         snew(protdata->sel_ntdb, NTERPAIRS);
 834         snew(protdata->sel_ctdb, NTERPAIRS);
 835
 836         if (nntdb >= 4 && nctdb >= 2)
 837         {
 838             /* Yuk, yuk, hardcoded default termini selections !!! */
 839             if (strncmp(*atoms->resinfo[atoms->atom[atoms->nr-1].resind].name, "PRO", 3) == 0)
 840             {
 841                 nt = 3;
 842             }
 843             else
 844             {
 845                 nt = 1;
 846             }
 847             ct = 1;
 848         }
 849         else
 850         {
 851             nt = 0;
 852             ct = 0;
 853         }
 854         protdata->sel_ntdb[0] = &(protdata->ntdb[nt]);
 855         protdata->sel_ctdb[0] = &(protdata->ctdb[ct]);
 856
 857         /* set terminal residue numbers: */
 858         snew(protdata->rN, NTERPAIRS);
 859         snew(protdata->rC, NTERPAIRS);
 860         protdata->rN[0] = 0;
 861         protdata->rC[0] = atoms->atom[atoms->nr-1].resind;
 862
 863         /* keep unprotonated topology: */
 864         protdata->upatoms = atoms;
 865         /* we don't have these yet: */
 866         protdata->patoms = NULL;
 867         bUpdate_pdba     = TRUE;
 868         bKeep_old_pdba   = TRUE;
 869
 870         /* clear hackblocks: */
 871         protdata->nab = NULL;
 872         protdata->ab  = NULL;
 873
 874         /* set flag to show we're initialized: */
 875         protdata->bInit = TRUE;
 876     }
 877     else
 878     {
 879         if (debug)
 880         {
 881             fprintf(debug, "protonate: using available protdata\n");
 882         }
 883         /* add_h will need the unprotonated topology again: */
 884         atoms          = protdata->upatoms;
 885         bUpdate_pdba   = FALSE;
 886         bKeep_old_pdba = FALSE;
 887     }
 888
 889     /* now protonate */
 890     nadd = add_h(&atoms, xptr, protdata->nah, protdata->ah,
 891                  NTERPAIRS, protdata->sel_ntdb, protdata->sel_ctdb,
 892                  protdata->rN, protdata->rC, TRUE,
 893                  &protdata->nab, &protdata->ab, bUpdate_pdba, bKeep_old_pdba);
 894     if (!protdata->patoms)
 895     {
 896         /* store protonated topology */
 897         protdata->patoms = atoms;
 898     }
 899     *atomsptr = protdata->patoms;
 900     if (debug)
 901     {
 902         fprintf(debug, "natoms: %d -> %d (nadd=%d)\n",
 903                 protdata->upatoms->nr, protdata->patoms->nr, nadd);
 904     }
 905     return nadd;
 906 #undef NTERPAIRS
 907 }