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