2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
38 /* This file is completely threadsafe - keep it that way! */
55 static int in_strings(char *key,int nstr,const char **str)
59 for(j=0; (j<nstr); j++)
60 if (strcmp(str[j],key) == 0)
66 static gmx_bool hbond(rvec x[],int i,int j,real distance)
68 real tol = distance*distance;
71 rvec_sub(x[i],x[j],tmp);
73 return (iprod(tmp,tmp) < tol);
76 static void chk_allhb(t_atoms *pdba,rvec x[],t_blocka *hb,
77 gmx_bool donor[],gmx_bool accept[],real dist)
82 snew(hb->index,natom+1);
89 for(i=0; (i<natom); i++) {
91 for(j=i+1; (j<natom); j++)
92 if ((accept[j]) && (hbond(x,i,j,dist)))
96 for(j=i+1; (j<natom); j++)
97 if ((donor[j]) && (hbond(x,i,j,dist)))
105 static void pr_hbonds(FILE *fp,t_blocka *hb,t_atoms *pdba)
109 fprintf(fp,"Dumping all hydrogen bonds!\n");
110 for(i=0; (i<hb->nr); i++) {
113 for(j=j0; (j<j1); j++) {
115 fprintf(fp,"%5s%4d%5s - %5s%4d%5s\n",
116 *pdba->resinfo[pdba->atom[i].resind].name,
117 pdba->resinfo[pdba->atom[i].resind].nr,*pdba->atomname[i],
118 *pdba->resinfo[pdba->atom[k].resind].name,
119 pdba->resinfo[pdba->atom[k].resind].nr,*pdba->atomname[k]);
124 static gmx_bool chk_hbonds(int i,t_atoms *pdba, rvec x[],
125 gmx_bool ad[],gmx_bool hbond[],rvec xh,
126 real angle,real dist)
135 ri = pdba->atom[i].resind;
137 for(j=0; (j<natom); j++) {
138 /* Check whether the other atom is a donor/acceptor and not i */
139 if ((ad[j]) && (j != i)) {
140 /* Check whether the other atom is on the same ring as well */
141 if ((pdba->atom[j].resind != ri) ||
142 ((strcmp(*pdba->atomname[j],"ND1") != 0) &&
143 (strcmp(*pdba->atomname[j],"NE2") != 0))) {
145 d2 = distance2(x[i],x[j]);
146 rvec_sub(x[i],xh,nh);
147 rvec_sub(x[aj],xh,oh);
148 a = RAD2DEG * acos(cos_angle(nh,oh));
149 if ((d2 < dist2) && (a > angle)) {
152 "HBOND between %s%d-%s and %s%d-%s is %g nm, %g deg\n",
153 *pdba->resinfo[pdba->atom[i].resind].name,
154 pdba->resinfo[pdba->atom[i].resind].nr,*pdba->atomname[i],
155 *pdba->resinfo[pdba->atom[aj].resind].name,
156 pdba->resinfo[pdba->atom[aj].resind].nr,*pdba->atomname[aj],
167 static void calc_ringh(rvec xattach,rvec xb,rvec xc,rvec xh)
172 /* Add a proton on a ring to atom attach at distance 0.1 nm */
173 rvec_sub(xattach,xb,tab);
174 rvec_sub(xattach,xc,tac);
175 rvec_add(tab,tac,xh);
178 rvec_inc(xh,xattach);
181 void set_histp(t_atoms *pdba,rvec *x,real angle,real dist){
182 static const char *prot_acc[] = {
183 "O", "OD1", "OD2", "OE1", "OE2", "OG", "OG1", "OH", "OW"
185 #define NPA asize(prot_acc)
186 static const char *prot_don[] = {
187 "N", "NH1", "NH2", "NE", "ND1", "ND2", "NE2", "NZ", "OG", "OG1", "OH", "NE1", "OW"
189 #define NPD asize(prot_don)
191 gmx_bool *donor,*acceptor;
192 gmx_bool *hbond,bHaveH=FALSE;
196 int i,j,nd,na,aj,hisind,his0,type=-1;
197 int nd1,ne2,cg,cd2,ce1;
206 gmx_strcasecmp(*pdba->resinfo[pdba->atom[i].resind].name,"HIS") != 0)
215 /* A histidine residue exists that requires automated assignment, so
216 * doing the analysis of donors and acceptors is worthwhile. */
218 "Analysing hydrogen-bonding network for automated assigment of histidine\n"
222 snew(acceptor,natom);
227 for(j=0; (j<natom); j++) {
228 if (in_strings(*pdba->atomname[j],NPA,prot_acc) != -1) {
232 if (in_strings(*pdba->atomname[j],NPD,prot_don) != -1) {
237 fprintf(stderr," %d donors and %d acceptors were found.\n",nd,na);
238 chk_allhb(pdba,x,hb,donor,acceptor,dist);
240 pr_hbonds(debug,hb,pdba);
241 fprintf(stderr,"There are %d hydrogen bonds\n",hb->nra);
243 /* Now do the HIS stuff */
247 if (gmx_strcasecmp(*pdba->resinfo[pdba->atom[i].resind].name,"HIS") != 0)
253 if (pdba->atom[i].resind != hisind) {
254 hisind=pdba->atom[i].resind;
256 /* Find the atoms in the ring */
257 nd1=ne2=cg=cd2=ce1=-1;
258 while (i<natom && pdba->atom[i].resind==hisind) {
259 atomnm = *pdba->atomname[i];
260 if (strcmp(atomnm,"CD2") == 0)
262 else if (strcmp(atomnm,"CG") == 0)
264 else if (strcmp(atomnm,"CE1") == 0)
266 else if (strcmp(atomnm,"ND1") == 0)
268 else if (strcmp(atomnm,"NE2") == 0)
274 if (!((cg == -1 ) || (cd2 == -1) || (ce1 == -1) ||
275 (nd1 == -1) || (ne2 == -1))) {
276 calc_ringh(x[nd1],x[cg],x[ce1],xh1);
277 calc_ringh(x[ne2],x[ce1],x[cd2],xh2);
279 bHDd = chk_hbonds(nd1,pdba,x,acceptor,hbond,xh1,angle,dist);
280 chk_hbonds(nd1,pdba,x,donor,hbond,xh1,angle,dist);
281 bHEd = chk_hbonds(ne2,pdba,x,acceptor,hbond,xh2,angle,dist);
282 chk_hbonds(ne2,pdba,x,donor,hbond,xh2,angle,dist);
292 fprintf(stderr,"Will use %s for residue %d\n",
293 hh[type],pdba->resinfo[hisind].nr);
296 gmx_fatal(FARGS,"Incomplete ring in HIS%d",
297 pdba->resinfo[hisind].nr);
299 snew(pdba->resinfo[hisind].rtp,1);
300 *pdba->resinfo[hisind].rtp = strdup(hh[type]);