3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
52 #include "gmx_fatal.h"
63 static void copy_atom(t_atoms *atoms1,int a1,t_atoms *atoms2,int a2)
65 atoms2->atom[a2] = atoms1->atom[a1];
66 snew(atoms2->atomname[a2],1);
67 *atoms2->atomname[a2]=strdup(*atoms1->atomname[a1]);
70 static atom_id pdbasearch_atom(const char *name,int resind,t_atoms *pdba,
71 const char *searchtype,gmx_bool bAllowMissing)
75 for(i=0; (i<pdba->nr) && (pdba->atom[i].resind != resind); i++)
78 return search_atom(name,i,pdba,
79 searchtype,bAllowMissing);
82 static void hacksearch_atom(int *ii, int *jj, char *name,
83 int nab[], t_hack *ab[],
84 int resind, t_atoms *pdba)
93 for(i=0; (i<pdba->nr) && (pdba->atom[i].resind != resind); i++)
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) {
105 void dump_ab(FILE *out,int natom,int nab[], t_hack *ab[], gmx_bool bHeader)
109 #define SS(s) (s)?(s):"-"
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),
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]);
127 static t_hackblock *get_hackblocks(t_atoms *pdba, int nah, t_hackblock ah[],
129 t_hackblock **ntdb, t_hackblock **ctdb,
133 t_hackblock *hb,*ahptr;
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]]);
142 if (ctdb[i] != NULL) {
143 merge_t_hackblock(ctdb[i], &hb[rC[i]]);
146 /* then the whole hdb */
147 for(rnr=0; rnr < pdba->nres; rnr++) {
148 ahptr=search_h_db(nah,ah,*pdba->resinfo[rnr].rtp);
150 if (hb[rnr].name==NULL) {
151 hb[rnr].name=strdup(ahptr->name);
153 merge_hacks(ahptr, &hb[rnr]);
159 static void expand_hackblocks_one(t_hackblock *hbr, char *atomname,
160 int *nabi, t_hack **abi, gmx_bool bN, gmx_bool bC)
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) */
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]=='-';
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]=='+';
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 */
185 fprintf(debug," %s",hbr->hack[j].oname?hbr->hack[j].oname:hbr->hack[j].AI);
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 */
195 fprintf(debug," +%dh",hbr->hack[j].nr);
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';
210 fprintf(debug,"Hack '%s' %d, replacing nname '%s' with '%s' (old name '%s')\n",
212 (*abi)[*nabi + k].nname,hbr->hack[j].nname,
213 (*abi)[*nabi + k].oname ? (*abi)[*nabi + k].oname : "");
215 sfree((*abi)[*nabi + k].nname);
216 (*abi)[*nabi + k].nname=strdup(hbr->hack[j].nname);
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';
239 (*nabi) += hbr->hack[j].nr;
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,
252 static void expand_hackblocks(t_atoms *pdba, t_hackblock hb[],
253 int nab[], t_hack *ab[],
254 int nterpairs, int *rN, int *rC)
259 for(i=0; i < pdba->nr; i++) {
261 for(j=0; j<nterpairs && !bN; j++)
262 bN = pdba->atom[i].resind==rN[j];
264 for(j=0; j<nterpairs && !bC; j++)
265 bC = pdba->atom[i].resind==rC[j];
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);
271 if (debug) fprintf(debug,"\n");
274 static int check_atoms_present(t_atoms *pdba, int nab[], t_hack *ab[])
276 int i, j, k, d, rnr, nadd;
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 ) {
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);
289 /* We found the added atom. */
290 ab[i][j].bAlreadyPresent = TRUE;
292 fprintf(debug,"Atom '%s' in residue '%s' %d is already present\n",
294 *pdba->resinfo[rnr].name,pdba->resinfo[rnr].nr);
297 ab[i][j].bAlreadyPresent = FALSE;
298 /* count how many atoms we'll add */
301 } else if ( ab[i][j].nname==NULL ) {
311 static void calc_all_pos(t_atoms *pdba, rvec x[], int nab[], t_hack *ab[],
312 gmx_bool bCheckMissing)
314 int i, j, ii, jj, m, ia, d, rnr,l=0;
316 rvec xa[4]; /* control atoms for calc_h_pos */
317 rvec xh[MAXH]; /* hydrogen positions from calc_h_pos */
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) {
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",
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);
336 copy_rvec(ab[ii][jj].newx, xa[m]);
340 gmx_fatal(FARGS,"Atom %s not found in residue %s %d"
342 " while adding hydrogens",
344 *pdba->resinfo[rnr].name,
345 pdba->resinfo[rnr].nr,
346 *pdba->resinfo[rnr].rtp);
350 copy_rvec(x[ia], xa[m]);
354 for(m=0; (m<MAXH); m++)
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;
371 static void free_ab(int natoms,int *nab,t_hack **ab)
375 for(i=0; i<natoms; i++) {
376 free_t_hack(nab[i], &ab[i]);
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)
389 t_atoms *newpdba=NULL,*pdba=NULL;
391 int i,newi,j,d,natoms,nalreadypresent;
398 /* set flags for adding hydrogens (according to hdb) */
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
409 if (debug) fprintf(debug,"pointer to ab found\n");
415 /* WOW, everything was already figured out */
416 bUpdate_pdba = FALSE;
417 if (debug) fprintf(debug,"pointer to non-null ab found\n");
419 /* We'll have to do all the hard work */
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);
425 /* expand the hackblocks to atom level */
428 expand_hackblocks(pdba, hb, nab, ab, nterpairs, rN, rC);
429 free_t_hackblock(pdba->nres, &hb);
433 fprintf(debug,"before calc_all_pos\n");
434 dump_ab(debug, natoms, nab, ab, TRUE);
437 /* Now calc the positions */
438 calc_all_pos(pdba, *xptr, nab, ab, bCheckMissing);
441 fprintf(debug,"after calc_all_pos\n");
442 dump_ab(debug, natoms, nab, ab, TRUE);
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);
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);
455 /* Copy old atoms, making space for new ones */
457 init_t_atoms(newpdba,natoms+nadd,FALSE);
458 newpdba->nres = pdba->nres;
459 sfree(newpdba->resinfo);
460 newpdba->resinfo = pdba->resinfo;
464 if (debug) fprintf(debug,"snew xn for %d old + %d new atoms %d total)\n",
465 natoms, nadd, natoms+nadd);
468 /* There is nothing to do: return now */
470 free_ab(natoms,nab,ab);
476 snew(xn,natoms+nadd);
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");*/
484 srenew(xn,natoms+nadd);
486 srenew(newpdba->atom,natoms+nadd);
487 srenew(newpdba->atomname,natoms+nadd);
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]);
496 copy_atom(pdba,i, newpdba,newi);
497 copy_rvec((*xptr)[i],xn[newi]);
498 /* process the hacks for this atom */
500 for(j=0; j<nab[i]; j++) {
501 if ( ab[i][j].oname==NULL ) { /* add */
503 if (newi >= natoms+nadd) {
504 /* gmx_fatal(FARGS,"Not enough space for adding atoms");*/
506 srenew(xn,natoms+nadd);
508 srenew(newpdba->atom,natoms+nadd);
509 srenew(newpdba->atomname,natoms+nadd);
514 newpdba->atom[newi].resind=pdba->atom[i].resind;
516 if (debug) fprintf(debug," + %d",newi+1);
518 if (ab[i][j].nname != NULL &&
519 (ab[i][j].oname == NULL ||
520 strcmp(ab[i][j].oname,*newpdba->atomname[newi]) == 0)) {
522 if (ab[i][j].oname == NULL && ab[i][j].bAlreadyPresent) {
523 /* This atom is already present, copy it from the input. */
526 copy_atom(pdba,i+nalreadypresent, newpdba,newi);
528 copy_rvec((*xptr)[i+nalreadypresent],xn[newi]);
532 fprintf(debug,"Replacing %d '%s' with (old name '%s') %s\n",
534 (newpdba->atomname[newi] && *newpdba->atomname[newi]) ? *newpdba->atomname[newi] : "",
535 ab[i][j].oname ? ab[i][j].oname : "",
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; */
547 copy_rvec(ab[i][j].newx, xn[newi]);
549 if (bUpdate_pdba && debug)
550 fprintf(debug," %s %g %g",*newpdba->atomname[newi],
551 newpdba->atom[newi].m,newpdba->atom[newi].q);
555 i += nalreadypresent;
556 if (debug) fprintf(debug,"\n");
567 free_ab(natoms,nab,ab);
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]); */
577 sfree(pdba->atomname);
579 sfree(pdba->pdbinfo);
592 void deprotonate(t_atoms *atoms,rvec *x)
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]);
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)
617 /* Here we loop to be able to add atoms to added atoms.
618 * We should not check for missing atoms here.
624 nnew = add_h_low(pdbaptr,xptr,nah,ah,nterpairs,ntdb,ctdb,rN,rC,FALSE,
625 nabptr,abptr,bUpdate_pdba,bKeep_old_pdba);
628 gmx_fatal(FARGS,"More than 100 iterations of add_h. Maybe you are trying to replace an added atom (this is not supported)?");
630 } while (nnew > nold);
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);
641 int protonate(t_atoms **atomsptr,rvec **xptr,t_protonate *protdata)
645 gmx_bool bUpdate_pdba,bKeep_old_pdba;
646 int nntdb,nctdb,nt,ct;
650 if ( !protdata->bInit ) {
652 fprintf(debug,"protonate: Initializing protdata\n");
655 /* set forcefield to use: */
656 strcpy(protdata->FF,"gmx2.ff");
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);
664 gmx_fatal(FARGS,"no N-terminus database");
666 nctdb = read_ter_db(protdata->FF,'c',&protdata->ctdb,protdata->atype);
668 gmx_fatal(FARGS,"no C-terminus database");
671 /* set terminus types: -NH3+ (different for Proline) and -COO- */
673 snew(protdata->sel_ntdb, NTERPAIRS);
674 snew(protdata->sel_ctdb, NTERPAIRS);
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) {
688 protdata->sel_ntdb[0]=&(protdata->ntdb[nt]);
689 protdata->sel_ctdb[0]=&(protdata->ctdb[ct]);
691 /* set terminal residue numbers: */
692 snew(protdata->rN, NTERPAIRS);
693 snew(protdata->rC, NTERPAIRS);
695 protdata->rC[0]=atoms->atom[atoms->nr-1].resind;
697 /* keep unprotonated topology: */
698 protdata->upatoms = atoms;
699 /* we don't have these yet: */
700 protdata->patoms = NULL;
704 /* clear hackblocks: */
708 /* set flag to show we're initialized: */
709 protdata->bInit=TRUE;
712 fprintf(debug,"protonate: using available protdata\n");
714 /* add_h will need the unprotonated topology again: */
715 atoms = protdata->upatoms;
717 bKeep_old_pdba=FALSE;
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;
729 *atomsptr = protdata->patoms;
731 fprintf(debug,"natoms: %d -> %d (nadd=%d)\n",
732 protdata->upatoms->nr, protdata->patoms->nr, nadd);