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