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