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