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.
54 real ellipticity(int nres,t_bb bb[])
60 static const t_ppwstr ppw[] = {
74 #define NPPW asize(ppw)
81 for(i=0; (i<nres); i++) {
84 for(j=0; (j<NPPW); j++) {
85 pp2=sqr(phi-ppw[j].phi)+sqr(psi-ppw[j].psi);
96 real ahx_len(int gnx,atom_id index[],rvec x[],matrix box)
97 /* Assume we have a list of Calpha atoms only! */
101 rvec_sub(x[index[0]],x[index[gnx-1]],dx);
106 real radius(FILE *fp,int nca,atom_id ca_index[],rvec x[])
107 /* Assume we have all the backbone */
113 for(i=0; (i<nca); i++) {
115 dl2=sqr(x[ai][XX])+sqr(x[ai][YY]);
118 fprintf(fp," %10g",dl2);
125 return sqrt(dlt/nca);
128 static real rot(rvec x1,rvec x2)
130 real phi1,dphi,cp,sp;
133 phi1=atan2(x1[YY],x1[XX]);
136 xx= cp*x2[XX]+sp*x2[YY];
137 yy=-sp*x2[XX]+cp*x2[YY];
139 dphi=RAD2DEG*atan2(yy,xx);
144 real twist(FILE *fp,int nca,atom_id caindex[],rvec x[])
151 for(i=1; (i<nca); i++) {
154 dphi=rot(x[a0],x[a1]);
164 real ca_phi(int gnx,atom_id index[],rvec x[],matrix box)
165 /* Assume we have a list of Calpha atoms only! */
168 int i,ai,aj,ak,al,t1,t2,t3;
169 rvec r_ij,r_kj,r_kl,m,n;
176 for(i=0; (i<gnx-4); i++) {
182 dih_angle(x[ai],x[aj],x[ak],x[al],NULL,
188 return (phitot/(gnx-4.0));
191 real dip(int nbb,atom_id bbind[],rvec x[],t_atom atom[])
198 for(i=0; (i<nbb); i++) {
201 for(m=0; (m<DIM); m++)
202 dipje[m]+=x[ai][m]*q;
207 real rise(int gnx,atom_id index[],rvec x[])
208 /* Assume we have a list of Calpha atoms only! */
216 for(i=1; (i<gnx); i++) {
223 return (ztot/(gnx-1.0));
226 void av_hblen(FILE *fp3,FILE *fp3a,
227 FILE *fp4,FILE *fp4a,
228 FILE *fp5,FILE *fp5a,
229 real t,int nres,t_bb bb[])
231 int i,n3=0,n4=0,n5=0;
234 for(i=0; (i<nres-3); i++)
236 fprintf(fp3a, "%10g",bb[i].d3);
240 fprintf(fp4a,"%10g",bb[i].d4);
245 fprintf(fp5a,"%10g",bb[i].d5);
250 fprintf(fp3,"%10g %10g\n",t,d3/n3);
251 fprintf(fp4,"%10g %10g\n",t,d4/n4);
252 fprintf(fp5,"%10g %10g\n",t,d5/n5);
259 void av_phipsi(FILE *fphi,FILE *fpsi,FILE *fphi2,FILE *fpsi2,
260 real t,int nres,t_bb bb[])
265 fprintf(fphi2,"%10g",t);
266 fprintf(fpsi2,"%10g",t);
267 for(i=0; (i<nres); i++)
271 fprintf(fphi2," %10g",bb[i].phi);
272 fprintf(fpsi2," %10g",bb[i].psi);
275 fprintf(fphi,"%10g %10g\n",t,(phi/n));
276 fprintf(fpsi,"%10g %10g\n",t,(psi/n));
281 static void set_ahcity(int nbb,t_bb bb[])
286 for(n=0; (n<nbb); n++) {
287 pp2=sqr(bb[n].phi-PHI_AHX)+sqr(bb[n].psi-PSI_AHX);
291 if ((bb[n].d4 < 0.36) || ((n > 0) && bb[n-1].bHelix))
297 t_bb *mkbbind(const char *fn,int *nres,int *nbb,int res0,
298 int *nall,atom_id **index,
299 char ***atomname,t_atom atom[],
302 static const char * bb_nm[] = { "N", "H", "CA", "C", "O", "HN" };
303 #define NBB asize(bb_nm)
306 int ai,i,i0,i1,j,k,ri,rnr,gnx,r0,r1;
308 fprintf(stderr,"Please select a group containing the entire backbone\n");
309 rd_index(fn,1,&gnx,index,&grpname);
311 fprintf(stderr,"Checking group %s\n",grpname);
312 r0=r1=atom[(*index)[0]].resind;
313 for(i=1; (i<gnx); i++) {
314 r0=min(r0,atom[(*index)[i]].resind);
315 r1=max(r1,atom[(*index)[i]].resind);
318 fprintf(stderr,"There are %d residues\n",rnr);
320 for(i=0; (i<rnr); i++)
321 bb[i].N=bb[i].H=bb[i].CA=bb[i].C=bb[i].O=-1,bb[i].resno=res0+i;
323 for(i=j=0; (i<gnx); i++) {
325 ri=atom[ai].resind-r0;
326 if (strcmp(*(resinfo[ri].name),"PRO") == 0) {
327 if (strcmp(*(atomname[ai]),"CD") == 0)
330 for(k=0; (k<NBB); k++)
331 if (strcmp(bb_nm[k],*(atomname[ai])) == 0) {
341 /* No attempt to address the case where some weird input has both H and HN atoms in the group */
358 for(i0=0; (i0<rnr); i0++) {
359 if ((bb[i0].N != -1) && (bb[i0].H != -1) &&
361 (bb[i0].C != -1) && (bb[i0].O != -1))
364 for(i1=rnr-1; (i1>=0); i1--) {
365 if ((bb[i1].N != -1) && (bb[i1].H != -1) &&
367 (bb[i1].C != -1) && (bb[i1].O != -1))
375 for(i=i0; (i<i1); i++) {
376 bb[i].Cprev=bb[i-1].C;
377 bb[i].Nnext=bb[i+1].N;
380 fprintf(stderr,"There are %d complete backbone residues (from %d to %d)\n",
381 rnr,bb[i0].resno,bb[i1].resno);
383 gmx_fatal(FARGS,"Zero complete backbone residues were found, cannot proceed");
384 for(i=0; (i<rnr); i++,i0++)
388 for(i=0; (i<rnr); i++) {
389 ri=atom[bb[i].CA].resind;
390 sprintf(bb[i].label,"%s%d",*(resinfo[ri].name),ri+res0);
394 *nbb=rnr*asize(bb_nm);
399 real pprms(FILE *fp,int nbb,t_bb bb[])
405 for(i=n=0; (i<nbb); i++) {
407 rms=sqrt(bb[i].pprms2);
410 fprintf(fp,"%10g ",rms);
415 rms=sqrt(rms2/n-sqr(rmst/n));
420 void calc_hxprops(int nres,t_bb bb[],rvec x[],matrix box)
422 int i,ao,an,t1,t2,t3;
423 rvec dx,r_ij,r_kj,r_kl,m,n;
426 for(i=0; (i<nres); i++) {
428 bb[i].d4=bb[i].d3=bb[i].d5=0;
431 rvec_sub(x[ao],x[an],dx);
436 rvec_sub(x[ao],x[an],dx);
441 rvec_sub(x[ao],x[an],dx);
446 dih_angle(x[bb[i].Cprev],x[bb[i].N],x[bb[i].CA],x[bb[i].C],NULL,
450 dih_angle(x[bb[i].N],x[bb[i].CA],x[bb[i].C],x[bb[i].Nnext],NULL,
453 bb[i].pprms2=sqr(bb[i].phi-PHI_AHX)+sqr(bb[i].psi-PSI_AHX);
456 1.4*sin((bb[i].psi+138.0)*DEG2RAD) -
457 4.1*cos(2.0*DEG2RAD*(bb[i].psi+138.0)) +
458 2.0*cos(2.0*DEG2RAD*(bb[i].phi+30.0));
462 static void check_ahx(int nres,t_bb bb[],rvec x[],
463 int *hstart,int *hend)
465 int h0,h1,h0sav,h1sav;
470 for(; (!bb[h0].bHelix) && (h0<nres-4); h0++)
472 for(h1=h0; bb[h1+1].bHelix && (h1<nres-1); h1++)
475 /*fprintf(stderr,"Helix from %d to %d\n",h0,h1);*/
476 if (h1-h0 > h1sav-h0sav) {
482 } while (h1 < nres-1);
487 void do_start_end(int nres,t_bb bb[],rvec x[],int *nbb,atom_id bbindex[],
488 int *nca,atom_id caindex[],
489 gmx_bool bRange,int rStart,int rEnd)
491 int i,j,hstart=0,hend=0;
494 for(i=0; (i<nres); i++) {
495 if ((bb[i].resno >= rStart) && (bb[i].resno <= rEnd))
497 if (bb[i].resno == rStart)
499 if (bb[i].resno == rEnd)
504 /* Find start and end of longest helix fragment */
505 check_ahx(nres,bb,x,&hstart,&hend);
507 fprintf(stderr,"helix from: %d through %d\n",
508 bb[hstart].resno,bb[hend].resno);
510 for(j=0,i=hstart; (i<=hend); i++) {
511 bbindex[j++]=bb[i].N;
512 bbindex[j++]=bb[i].H;
513 bbindex[j++]=bb[i].CA;
514 bbindex[j++]=bb[i].C;
515 bbindex[j++]=bb[i].O;
516 caindex[i-hstart]=bb[i].CA;
519 *nca=(hend-hstart+1);
522 void pr_bb(FILE *fp,int nres,t_bb bb[])
527 fprintf(fp,"%3s %3s %3s %3s %3s %7s %7s %7s %7s %7s %3s\n",
528 "AA","N","Ca","C","O","Phi","Psi","D3","D4","D5","Hx?");
529 for(i=0; (i<nres); i++) {
530 fprintf(fp,"%3d %3d %3d %3d %3d %7.2f %7.2f %7.3f %7.3f %7.3f %3s\n",
531 bb[i].resno,bb[i].N,bb[i].CA,bb[i].C,bb[i].O,
532 bb[i].phi,bb[i].psi,bb[i].d3,bb[i].d4,bb[i].d5,
533 bb[i].bHelix ? "Yes" : "No");