src/gromacs/gmxpreprocess/genhydro.c

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