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 * Green Red Orange Magenta Azure Cyan Skyblue
51 int nhelix(int nres,t_bb bb[])
55 for(i=n=0; (i<nres); i++)
61 real ellipticity(int nres,t_bb bb[])
67 static const t_ppwstr ppw[] = {
81 #define NPPW asize(ppw)
88 for(i=0; (i<nres); i++) {
91 for(j=0; (j<NPPW); j++) {
92 pp2=sqr(phi-ppw[j].phi)+sqr(psi-ppw[j].psi);
103 real ahx_len(int gnx,atom_id index[],rvec x[],matrix box)
104 /* Assume we have a list of Calpha atoms only! */
108 rvec_sub(x[index[0]],x[index[gnx-1]],dx);
113 real radius(FILE *fp,int nca,atom_id ca_index[],rvec x[])
114 /* Assume we have all the backbone */
120 for(i=0; (i<nca); i++) {
122 dl2=sqr(x[ai][XX])+sqr(x[ai][YY]);
125 fprintf(fp," %10g",dl2);
132 return sqrt(dlt/nca);
135 static real rot(rvec x1,rvec x2)
137 real phi1,dphi,cp,sp;
140 phi1=atan2(x1[YY],x1[XX]);
143 xx= cp*x2[XX]+sp*x2[YY];
144 yy=-sp*x2[XX]+cp*x2[YY];
146 dphi=RAD2DEG*atan2(yy,xx);
151 real twist(FILE *fp,int nca,atom_id caindex[],rvec x[])
158 for(i=1; (i<nca); i++) {
161 dphi=rot(x[a0],x[a1]);
171 real ca_phi(int gnx,atom_id index[],rvec x[],matrix box)
172 /* Assume we have a list of Calpha atoms only! */
175 int i,ai,aj,ak,al,t1,t2,t3;
176 rvec r_ij,r_kj,r_kl,m,n;
183 for(i=0; (i<gnx-4); i++) {
189 dih_angle(x[ai],x[aj],x[ak],x[al],NULL,
195 return (phitot/(gnx-4.0));
198 real dip(int nbb,atom_id bbind[],rvec x[],t_atom atom[])
205 for(i=0; (i<nbb); i++) {
208 for(m=0; (m<DIM); m++)
209 dipje[m]+=x[ai][m]*q;
214 real rise(int gnx,atom_id index[],rvec x[])
215 /* Assume we have a list of Calpha atoms only! */
223 for(i=1; (i<gnx); i++) {
230 return (ztot/(gnx-1.0));
233 void av_hblen(FILE *fp3,FILE *fp3a,
234 FILE *fp4,FILE *fp4a,
235 FILE *fp5,FILE *fp5a,
236 real t,int nres,t_bb bb[])
238 int i,n3=0,n4=0,n5=0;
241 for(i=0; (i<nres-3); i++)
243 fprintf(fp3a, "%10g",bb[i].d3);
247 fprintf(fp4a,"%10g",bb[i].d4);
252 fprintf(fp5a,"%10g",bb[i].d5);
257 fprintf(fp3,"%10g %10g\n",t,d3/n3);
258 fprintf(fp4,"%10g %10g\n",t,d4/n4);
259 fprintf(fp5,"%10g %10g\n",t,d5/n5);
266 void av_phipsi(FILE *fphi,FILE *fpsi,FILE *fphi2,FILE *fpsi2,
267 real t,int nres,t_bb bb[])
272 fprintf(fphi2,"%10g",t);
273 fprintf(fpsi2,"%10g",t);
274 for(i=0; (i<nres); i++)
278 fprintf(fphi2," %10g",bb[i].phi);
279 fprintf(fpsi2," %10g",bb[i].psi);
282 fprintf(fphi,"%10g %10g\n",t,(phi/n));
283 fprintf(fpsi,"%10g %10g\n",t,(psi/n));
288 static void set_ahcity(int nbb,t_bb bb[])
293 for(n=0; (n<nbb); n++) {
294 pp2=sqr(bb[n].phi-PHI_AHX)+sqr(bb[n].psi-PSI_AHX);
298 if ((bb[n].d4 < 0.36) || ((n > 0) && bb[n-1].bHelix))
304 t_bb *mkbbind(const char *fn,int *nres,int *nbb,int res0,
305 int *nall,atom_id **index,
306 char ***atomname,t_atom atom[],
309 static const char * bb_nm[] = { "N", "H", "CA", "C", "O" };
310 #define NBB asize(bb_nm)
313 int ai,i,i0,i1,j,k,ri,rnr,gnx,r0,r1;
315 fprintf(stderr,"Please select a group containing the entire backbone\n");
316 rd_index(fn,1,&gnx,index,&grpname);
318 fprintf(stderr,"Checking group %s\n",grpname);
319 r0=r1=atom[(*index)[0]].resind;
320 for(i=1; (i<gnx); i++) {
321 r0=min(r0,atom[(*index)[i]].resind);
322 r1=max(r1,atom[(*index)[i]].resind);
325 fprintf(stderr,"There are %d residues\n",rnr);
327 for(i=0; (i<rnr); i++)
328 bb[i].N=bb[i].H=bb[i].CA=bb[i].C=bb[i].O=-1,bb[i].resno=res0+i;
330 for(i=j=0; (i<gnx); i++) {
332 ri=atom[ai].resind-r0;
333 if (strcmp(*(resinfo[ri].name),"PRO") == 0) {
334 if (strcmp(*(atomname[ai]),"CD") == 0)
337 for(k=0; (k<NBB); k++)
338 if (strcmp(bb_nm[k],*(atomname[ai])) == 0) {
363 for(i0=0; (i0<rnr); i0++) {
364 if ((bb[i0].N != -1) && (bb[i0].H != -1) &&
366 (bb[i0].C != -1) && (bb[i0].O != -1))
369 for(i1=rnr-1; (i1>=0); i1--) {
370 if ((bb[i1].N != -1) && (bb[i1].H != -1) &&
372 (bb[i1].C != -1) && (bb[i1].O != -1))
380 for(i=i0; (i<i1); i++) {
381 bb[i].Cprev=bb[i-1].C;
382 bb[i].Nnext=bb[i+1].N;
385 fprintf(stderr,"There are %d complete backbone residues (from %d to %d)\n",
386 rnr,bb[i0].resno,bb[i1].resno);
388 gmx_fatal(FARGS,"rnr==0");
389 for(i=0; (i<rnr); i++,i0++)
393 for(i=0; (i<rnr); i++) {
394 ri=atom[bb[i].CA].resind;
395 sprintf(bb[i].label,"%s%d",*(resinfo[ri].name),ri+res0);
399 *nbb=rnr*asize(bb_nm);
404 real pprms(FILE *fp,int nbb,t_bb bb[])
410 for(i=n=0; (i<nbb); i++) {
412 rms=sqrt(bb[i].pprms2);
415 fprintf(fp,"%10g ",rms);
420 rms=sqrt(rms2/n-sqr(rmst/n));
425 void calc_hxprops(int nres,t_bb bb[],rvec x[],matrix box)
427 int i,ao,an,t1,t2,t3;
428 rvec dx,r_ij,r_kj,r_kl,m,n;
431 for(i=0; (i<nres); i++) {
433 bb[i].d4=bb[i].d3=bb[i].d5=0;
436 rvec_sub(x[ao],x[an],dx);
441 rvec_sub(x[ao],x[an],dx);
446 rvec_sub(x[ao],x[an],dx);
451 dih_angle(x[bb[i].Cprev],x[bb[i].N],x[bb[i].CA],x[bb[i].C],NULL,
455 dih_angle(x[bb[i].N],x[bb[i].CA],x[bb[i].C],x[bb[i].Nnext],NULL,
458 bb[i].pprms2=sqr(bb[i].phi-PHI_AHX)+sqr(bb[i].psi-PSI_AHX);
461 1.4*sin((bb[i].psi+138.0)*DEG2RAD) -
462 4.1*cos(2.0*DEG2RAD*(bb[i].psi+138.0)) +
463 2.0*cos(2.0*DEG2RAD*(bb[i].phi+30.0));
467 static void check_ahx(int nres,t_bb bb[],rvec x[],
468 int *hstart,int *hend)
470 int h0,h1,h0sav,h1sav;
475 for(; (!bb[h0].bHelix) && (h0<nres-4); h0++)
477 for(h1=h0; bb[h1+1].bHelix && (h1<nres-1); h1++)
480 /*fprintf(stderr,"Helix from %d to %d\n",h0,h1);*/
481 if (h1-h0 > h1sav-h0sav) {
487 } while (h1 < nres-1);
492 void do_start_end(int nres,t_bb bb[],rvec x[],int *nbb,atom_id bbindex[],
493 int *nca,atom_id caindex[],
494 gmx_bool bRange,int rStart,int rEnd)
496 int i,j,hstart=0,hend=0;
499 for(i=0; (i<nres); i++) {
500 if ((bb[i].resno >= rStart) && (bb[i].resno <= rEnd))
502 if (bb[i].resno == rStart)
504 if (bb[i].resno == rEnd)
509 /* Find start and end of longest helix fragment */
510 check_ahx(nres,bb,x,&hstart,&hend);
512 fprintf(stderr,"helix from: %d through %d\n",
513 bb[hstart].resno,bb[hend].resno);
515 for(j=0,i=hstart; (i<=hend); i++) {
516 bbindex[j++]=bb[i].N;
517 bbindex[j++]=bb[i].H;
518 bbindex[j++]=bb[i].CA;
519 bbindex[j++]=bb[i].C;
520 bbindex[j++]=bb[i].O;
521 caindex[i-hstart]=bb[i].CA;
524 *nca=(hend-hstart+1);
527 void pr_bb(FILE *fp,int nres,t_bb bb[])
532 fprintf(fp,"%3s %3s %3s %3s %3s %7s %7s %7s %7s %7s %3s\n",
533 "AA","N","Ca","C","O","Phi","Psi","D3","D4","D5","Hx?");
534 for(i=0; (i<nres); i++) {
535 fprintf(fp,"%3d %3d %3d %3d %3d %7.2f %7.2f %7.3f %7.3f %7.3f %3s\n",
536 bb[i].resno,bb[i].N,bb[i].CA,bb[i].C,bb[i].O,
537 bb[i].phi,bb[i].psi,bb[i].d3,bb[i].d4,bb[i].d5,
538 bb[i].bHelix ? "Yes" : "No");