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