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"
62 static void copy_atom(t_atoms *atoms1,int a1,t_atoms *atoms2,int a2)
64 atoms2->atom[a2] = atoms1->atom[a1];
65 snew(atoms2->atomname[a2],1);
66 *atoms2->atomname[a2]=strdup(*atoms1->atomname[a1]);
69 static atom_id pdbasearch_atom(const char *name,int resind,t_atoms *pdba,
70 const char *searchtype,bool bAllowMissing)
74 for(i=0; (i<pdba->nr) && (pdba->atom[i].resind != resind); i++)
77 return search_atom(name,i,pdba->nr,pdba->atom,pdba->atomname,
78 searchtype,bAllowMissing);
81 static void hacksearch_atom(int *ii, int *jj, char *name,
82 int nab[], t_hack *ab[],
83 int resind, t_atoms *pdba)
92 for(i=0; (i<pdba->nr) && (pdba->atom[i].resind != resind); i++)
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) {
104 void dump_ab(FILE *out,int natom,int nab[], t_hack *ab[], bool bHeader)
108 #define SS(s) (s)?(s):"-"
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),
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]);
126 static t_hackblock *get_hackblocks(t_atoms *pdba, int nah, t_hackblock ah[],
128 t_hackblock **ntdb, t_hackblock **ctdb,
132 t_hackblock *hb,*ahptr;
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]]);
141 if (ctdb[i] != NULL) {
142 merge_t_hackblock(ctdb[i], &hb[rC[i]]);
145 /* then the whole hdb */
146 for(rnr=0; rnr < pdba->nres; rnr++) {
147 ahptr=search_h_db(nah,ah,*pdba->resinfo[rnr].rtp);
149 if (hb[rnr].name==NULL)
150 hb[rnr].name=strdup(ahptr->name);
151 merge_hacks(ahptr, &hb[rnr]);
157 static void expand_hackblocks_one(t_hackblock *hbr, char *atomname,
158 int *nabi, t_hack **abi, bool bN, bool bC)
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) */
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);
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';
201 fprintf(debug,"Hack '%s' %d, replacing nname '%s' with '%s' (old name '%s')\n",
203 (*abi)[*nabi + k].nname,hbr->hack[j].nname,
204 (*abi)[*nabi + k].oname ? (*abi)[*nabi + k].oname : "");
206 sfree((*abi)[*nabi + k].nname);
207 (*abi)[*nabi + k].nname=strdup(hbr->hack[j].nname);
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';
230 (*nabi) += hbr->hack[j].nr;
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,
241 static void expand_hackblocks(t_atoms *pdba, t_hackblock hb[],
242 int nab[], t_hack *ab[],
243 int nterpairs, int *rN, int *rC)
248 for(i=0; i < pdba->nr; i++) {
250 for(j=0; j<nterpairs && !bN; j++)
251 bN = pdba->atom[i].resind==rN[j];
253 for(j=0; j<nterpairs && !bC; j++)
254 bC = pdba->atom[i].resind==rC[j];
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);
260 if (debug) fprintf(debug,"\n");
263 static int check_atoms_present(t_atoms *pdba, int nab[], t_hack *ab[])
265 int i, j, k, d, rnr, nadd;
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 ) {
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);
278 /* We found the added atom. */
279 ab[i][j].bAlreadyPresent = TRUE;
281 fprintf(debug,"Atom '%s' in residue '%s' %d is already present\n",
283 *pdba->resinfo[rnr].name,pdba->resinfo[rnr].nr);
286 ab[i][j].bAlreadyPresent = FALSE;
287 /* count how many atoms we'll add */
290 } else if ( ab[i][j].nname==NULL ) {
300 static void calc_all_pos(t_atoms *pdba, rvec x[], int nab[], t_hack *ab[],
303 int i, j, ii, jj, m, ia, d, rnr,l=0;
305 rvec xa[4]; /* control atoms for calc_h_pos */
306 rvec xh[MAXH]; /* hydrogen positions from calc_h_pos */
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) {
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",
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);
325 copy_rvec(ab[ii][jj].newx, xa[m]);
329 gmx_fatal(FARGS,"Atom %s not found in residue %s %d"
331 " while adding hydrogens",
333 *pdba->resinfo[rnr].name,
334 pdba->resinfo[rnr].nr,
335 *pdba->resinfo[rnr].rtp);
339 copy_rvec(x[ia], xa[m]);
342 for(m=0; (m<MAXH); m++)
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;
359 static void free_ab(int natoms,int *nab,t_hack **ab)
363 for(i=0; i<natoms; i++) {
364 free_t_hack(nab[i], &ab[i]);
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, bool bCheckMissing,
374 int **nabptr, t_hack ***abptr,
375 bool bUpdate_pdba, bool bKeep_old_pdba)
377 t_atoms *newpdba=NULL,*pdba=NULL;
379 int i,newi,j,d,natoms,nalreadypresent;
386 /* set flags for adding hydrogens (according to hdb) */
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
397 if (debug) fprintf(debug,"pointer to ab found\n");
403 /* WOW, everything was already figured out */
404 bUpdate_pdba = FALSE;
405 if (debug) fprintf(debug,"pointer to non-null ab found\n");
407 /* We'll have to do all the hard work */
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);
413 /* expand the hackblocks to atom level */
416 expand_hackblocks(pdba, hb, nab, ab, nterpairs, rN, rC);
417 free_t_hackblock(pdba->nres, &hb);
421 fprintf(debug,"before calc_all_pos\n");
422 dump_ab(debug, natoms, nab, ab, TRUE);
425 /* Now calc the positions */
426 calc_all_pos(pdba, *xptr, nab, ab, bCheckMissing);
429 fprintf(debug,"after calc_all_pos\n");
430 dump_ab(debug, natoms, nab, ab, TRUE);
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);
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);
443 /* Copy old atoms, making space for new ones */
445 init_t_atoms(newpdba,natoms+nadd,FALSE);
446 newpdba->nres = pdba->nres;
447 sfree(newpdba->resinfo);
448 newpdba->resinfo = pdba->resinfo;
452 if (debug) fprintf(debug,"snew xn for %d old + %d new atoms %d total)\n",
453 natoms, nadd, natoms+nadd);
456 /* There is nothing to do: return now */
458 free_ab(natoms,nab,ab);
464 snew(xn,natoms+nadd);
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");*/
472 srenew(xn,natoms+nadd);
474 srenew(newpdba->atom,natoms+nadd);
475 srenew(newpdba->atomname,natoms+nadd);
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]);
484 copy_atom(pdba,i, newpdba,newi);
485 copy_rvec((*xptr)[i],xn[newi]);
486 /* process the hacks for this atom */
488 for(j=0; j<nab[i]; j++) {
489 if ( ab[i][j].oname==NULL ) { /* add */
491 if (newi >= natoms+nadd) {
492 /* gmx_fatal(FARGS,"Not enough space for adding atoms");*/
494 srenew(xn,natoms+nadd);
496 srenew(newpdba->atom,natoms+nadd);
497 srenew(newpdba->atomname,natoms+nadd);
502 newpdba->atom[newi].resind=pdba->atom[i].resind;
504 if (debug) fprintf(debug," + %d",newi+1);
506 if (ab[i][j].nname != NULL &&
507 (ab[i][j].oname == NULL ||
508 strcmp(ab[i][j].oname,*newpdba->atomname[newi]) == 0)) {
510 if (ab[i][j].oname == NULL && ab[i][j].bAlreadyPresent) {
511 /* This atom is already present, copy it from the input. */
514 copy_atom(pdba,i+nalreadypresent, newpdba,newi);
516 copy_rvec((*xptr)[i+nalreadypresent],xn[newi]);
520 fprintf(debug,"Replacing %d '%s' with (old name '%s') %s\n",
522 (newpdba->atomname[newi] && *newpdba->atomname[newi]) ? *newpdba->atomname[newi] : "",
523 ab[i][j].oname ? ab[i][j].oname : "",
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; */
535 copy_rvec(ab[i][j].newx, xn[newi]);
537 if (bUpdate_pdba && debug)
538 fprintf(debug," %s %g %g",*newpdba->atomname[newi],
539 newpdba->atom[newi].m,newpdba->atom[newi].q);
543 i += nalreadypresent;
544 if (debug) fprintf(debug,"\n");
555 free_ab(natoms,nab,ab);
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]); */
565 sfree(pdba->atomname);
567 sfree(pdba->pdbinfo);
580 void deprotonate(t_atoms *atoms,rvec *x)
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]);
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, bool bAllowMissing,
600 int **nabptr, t_hack ***abptr,
601 bool bUpdate_pdba, bool bKeep_old_pdba)
605 /* Here we loop to be able to add atoms to added atoms.
606 * We should not check for missing atoms here.
612 nnew = add_h_low(pdbaptr,xptr,nah,ah,nterpairs,ntdb,ctdb,rN,rC,FALSE,
613 nabptr,abptr,bUpdate_pdba,bKeep_old_pdba);
616 gmx_fatal(FARGS,"More than 100 iterations of add_h. Maybe you are trying to replace an added atom (this is not supported)?");
618 } while (nnew > nold);
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);
629 int protonate(t_atoms **atomsptr,rvec **xptr,t_protonate *protdata)
633 bool bUpdate_pdba,bKeep_old_pdba;
634 int nntdb,nctdb,nt,ct;
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);
648 gmx_fatal(FARGS,"no n-terminus db");
649 nctdb = read_ter_db(protdata->FF,'c',&protdata->ctdb,protdata->atype);
651 gmx_fatal(FARGS,"no c-terminus db");
653 /* set terminus types: -NH3+ (different for Proline) and -COO- */
655 snew(protdata->sel_ntdb, NTERPAIRS);
656 snew(protdata->sel_ctdb, NTERPAIRS);
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)
669 protdata->sel_ntdb[0]=&(protdata->ntdb[nt]);
670 protdata->sel_ctdb[0]=&(protdata->ctdb[ct]);
672 /* set terminal residue numbers: */
673 snew(protdata->rN, NTERPAIRS);
674 snew(protdata->rC, NTERPAIRS);
676 protdata->rC[0]=atoms->atom[atoms->nr-1].resind;
678 /* keep unprotonated topology: */
679 protdata->upatoms = atoms;
680 /* we don't have these yet: */
681 protdata->patoms = NULL;
685 /* clear hackblocks: */
689 /* set flag to show we're initialized: */
690 protdata->bInit=TRUE;
692 if (debug) fprintf(debug,"protonate: using available protdata\n");
693 /* add_h will need the unprotonated topoloy again: */
694 atoms = protdata->upatoms;
696 bKeep_old_pdba=FALSE;
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);