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