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