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