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
53 void print_one(const output_env_t oenv
,const char *base
,const char *name
,
54 const char *title
, const char *ylabel
,int nf
,real time
[],
58 char buf
[256],t2
[256];
61 sprintf(buf
,"%s%s.xvg",base
,name
);
62 fprintf(stderr
,"\rPrinting %s ",buf
);
63 sprintf(t2
,"%s %s",title
,name
);
64 fp
=xvgropen(buf
,t2
,"Time (ps)",ylabel
,oenv
);
66 fprintf(fp
,"%10g %10g\n",time
[k
],data
[k
]);
70 static int calc_RBbin(real phi
, int multiplicity
, real core_frac
)
72 /* multiplicity and core_frac NOT used,
73 * just given to enable use of pt-to-fn in caller low_ana_dih_trans*/
74 static const real r30
= M_PI
/6.0;
75 static const real r90
= M_PI
/2.0;
76 static const real r150
= M_PI
*5.0/6.0;
78 if ((phi
< r30
) && (phi
> -r30
))
80 else if ((phi
> -r150
) && (phi
< -r90
))
82 else if ((phi
< r150
) && (phi
> r90
))
87 static int calc_Nbin(real phi
, int multiplicity
, real core_frac
)
89 static const real r360
= 360*DEG2RAD
;
90 real rot_width
, core_width
, core_offset
, low
, hi
;
92 /* with multiplicity 3 and core_frac 0.5
93 * 0<g(-)<120, 120<t<240, 240<g(+)<360
94 * 0< bin0 < 30, 30<bin1<90, 90<bin0<150, 150<bin2<210, 210<bin0<270, 270<bin3<330, 330<bin0<360
95 * so with multiplicity 3, bin1 is core g(-), bin2 is core t, bin3 is
96 core g(+), bin0 is between rotamers */
100 rot_width
= 360/multiplicity
;
101 core_width
= core_frac
* rot_width
;
102 core_offset
= (rot_width
- core_width
)/2.0 ;
103 for(bin
= 1 ; bin
<= multiplicity
; bin
++ ) {
104 low
= ((bin
- 1) * rot_width
) + core_offset
;
105 hi
= ((bin
- 1) * rot_width
) + core_offset
+ core_width
;
108 if ((phi
> low
) && (phi
< hi
))
114 void ana_dih_trans(const char *fn_trans
,const char *fn_histo
,
115 real
**dih
,int nframes
,int nangles
,
116 const char *grpname
,real t0
,real dt
,bool bRb
,
117 const output_env_t oenv
)
119 /* just a wrapper; declare extra args, then chuck away at end. */
123 int nlist
= nangles
;
128 for(k
=0; (k
<nangles
); k
++) {
132 low_ana_dih_trans(TRUE
, fn_trans
,TRUE
, fn_histo
, maxchi
,
133 dih
, nlist
, dlist
, nframes
,
134 nangles
, grpname
, xity
, t0
, dt
, bRb
, 0.5,oenv
);
140 void low_ana_dih_trans(bool bTrans
, const char *fn_trans
,
141 bool bHisto
, const char *fn_histo
, int maxchi
,
142 real
**dih
, int nlist
, t_dlist dlist
[], int nframes
,
143 int nangles
, const char *grpname
, int xity
[],
144 real t0
, real dt
, bool bRb
, real core_frac
,
145 const output_env_t oenv
)
150 int i
,j
,k
,Dih
,ntrans
;
153 real
*rot_occ
[NROT
] ;
154 int (*calc_bin
)(real
,int,real
);
156 /* Analysis of dihedral transitions */
157 fprintf(stderr
,"Now calculating transitions...\n");
164 for(k
=0;k
<NROT
;k
++) {
165 snew(rot_occ
[k
],nangles
);
166 for (i
=0; (i
<nangles
); i
++)
172 /* dih[i][j] is the dihedral angle i in frame j */
174 for (i
=0; (i
<nangles
); i
++)
179 mind
= maxd
= prev
= dih
[i
][0];
181 cur_bin
= calc_bin(dih
[i
][0],xity
[i
],core_frac
);
182 rot_occ
[cur_bin
][i
]++ ;
184 for (j
=1; (j
<nframes
); j
++)
186 new_bin
= calc_bin(dih
[i
][j
],xity
[i
],core_frac
);
187 rot_occ
[new_bin
][i
]++ ;
191 else if ((new_bin
!= 0) && (cur_bin
!= new_bin
)) {
198 /* why is all this md rubbish periodic? Remove 360 degree periodicity */
199 if ( (dih
[i
][j
] - prev
) > M_PI
)
201 else if ( (dih
[i
][j
] - prev
) < -M_PI
)
206 mind
= min(mind
, dih
[i
][j
]);
207 maxd
= max(maxd
, dih
[i
][j
]);
208 if ( (maxd
- mind
) > 2*M_PI
/3) /* or 120 degrees, assuming */
209 { /* multiplicity 3. Not so general.*/
212 maxd
= mind
= dih
[i
][j
]; /* get ready for next transition */
218 rot_occ
[k
][i
] /= nframes
;
220 fprintf(stderr
,"Total number of transitions: %10d\n",ntrans
);
222 ttime
= (dt
*nframes
*nangles
)/ntrans
;
223 fprintf(stderr
,"Time between transitions: %10.3f ps\n",ttime
);
226 /* new by grs - copy transitions from tr_h[] to dlist->ntr[]
227 * and rotamer populations from rot_occ to dlist->rot_occ[]
228 * based on fn histogramming in g_chi. diff roles for i and j here */
231 for (Dih
=0; (Dih
<NONCHI
+maxchi
); Dih
++) {
232 for(i
=0; (i
<nlist
); i
++) {
233 if (((Dih
< edOmega
) ) ||
234 ((Dih
== edOmega
) && (has_dihedral(edOmega
,&(dlist
[i
])))) ||
235 ((Dih
> edOmega
) && (dlist
[i
].atm
.Cn
[Dih
-NONCHI
+3] != -1))) {
236 /* grs debug printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
237 dlist
[i
].ntr
[Dih
] = tr_h
[j
] ;
239 dlist
[i
].rot_occ
[Dih
][k
] = rot_occ
[k
][j
] ;
245 /* end addition by grs */
248 sprintf(title
,"Number of transitions: %s",grpname
);
249 fp
=xvgropen(fn_trans
,title
,"Time (ps)","# transitions/timeframe",oenv
);
250 for(j
=0; (j
<nframes
); j
++) {
252 fprintf(fp
,"%10.3f %10d\n",tt
,tr_f
[j
]);
257 /* Compute histogram from # transitions per dihedral */
259 for(j
=0; (j
<nframes
); j
++)
261 for(i
=0; (i
<nangles
); i
++)
263 for(j
=nframes
; ((tr_f
[j
-1] == 0) && (j
>0)); j
--)
268 sprintf(title
,"Transition time: %s",grpname
);
269 fp
=xvgropen(fn_histo
,title
,"Time (ps)","#",oenv
);
270 for(i
=j
-1; (i
>0); i
--) {
272 fprintf(fp
,"%10.3f %10d\n",ttime
/i
,tr_f
[i
]);
284 void mk_multiplicity_lookup (int *xity
, int maxchi
, real
**dih
,
285 int nlist
, t_dlist dlist
[],int nangles
)
287 /* new by grs - for dihedral j (as in dih[j]) get multiplicity from dlist
288 * and store in xity[j]
295 for (Dih
=0; (Dih
<NONCHI
+maxchi
); Dih
++) {
296 for(i
=0; (i
<nlist
); i
++) {
297 strncpy(name
, dlist
[i
].name
,3) ;
299 if (((Dih
< edOmega
) ) ||
300 ((Dih
== edOmega
) && (has_dihedral(edOmega
,&(dlist
[i
])))) ||
301 ((Dih
> edOmega
) && (dlist
[i
].atm
.Cn
[Dih
-NONCHI
+3] != -1))) {
302 /* default - we will correct the rest below */
305 /* make omegas 2fold, though doesn't make much more sense than 3 */
306 if (Dih
== edOmega
&& (has_dihedral(edOmega
,&(dlist
[i
])))) {
310 /* dihedrals to aromatic rings, COO, CONH2 or guanidinium are 2fold*/
311 if (Dih
> edOmega
&& (dlist
[i
].atm
.Cn
[Dih
-NONCHI
+3] != -1)) {
312 if ( ((strstr(name
,"PHE") != NULL
) && (Dih
== edChi2
)) ||
313 ((strstr(name
,"TYR") != NULL
) && (Dih
== edChi2
)) ||
314 ((strstr(name
,"PTR") != NULL
) && (Dih
== edChi2
)) ||
315 ((strstr(name
,"TRP") != NULL
) && (Dih
== edChi2
)) ||
316 ((strstr(name
,"HIS") != NULL
) && (Dih
== edChi2
)) ||
317 ((strstr(name
,"GLU") != NULL
) && (Dih
== edChi3
)) ||
318 ((strstr(name
,"ASP") != NULL
) && (Dih
== edChi2
)) ||
319 ((strstr(name
,"GLN") != NULL
) && (Dih
== edChi3
)) ||
320 ((strstr(name
,"ASN") != NULL
) && (Dih
== edChi2
)) ||
321 ((strstr(name
,"ARG") != NULL
) && (Dih
== edChi4
)) ) {
330 fprintf(stderr
,"WARNING: not all dihedrals found in topology (only %d out of %d)!\n",
332 /* Check for remaining dihedrals */
333 for(;(j
< nangles
); j
++)
338 void mk_chi_lookup (int **lookup
, int maxchi
, real
**dih
,
339 int nlist
, t_dlist dlist
[])
342 /* by grs. should rewrite everything to use this. (but haven't,
343 * and at mmt only used in get_chi_product_traj
344 * returns the dihed number given the residue number (from-0)
345 * and chi (from-0) nr. -1 for chi undefined for that res (eg gly, ala..)*/
350 for (Dih
=0; (Dih
<NONCHI
+maxchi
); Dih
++) {
351 for(i
=0; (i
<nlist
); i
++) {
353 if (((Dih
< edOmega
) ) ||
354 ((Dih
== edOmega
) && (has_dihedral(edOmega
,&(dlist
[i
])))) ||
355 ((Dih
> edOmega
) && (dlist
[i
].atm
.Cn
[Dih
-NONCHI
+3] != -1))) {
356 /* grs debug printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
357 if (Dih
> edOmega
) {
362 lookup
[i
][Chi
] = -1 ;
370 void get_chi_product_traj (real
**dih
,int nframes
,int nangles
, int nlist
,
371 int maxchi
, t_dlist dlist
[], real time
[],
372 int **lookup
, int *xity
,bool bRb
, bool bNormalize
,
373 real core_frac
, bool bAll
, const char *fnall
,
374 const output_env_t oenv
)
377 bool bRotZero
, bHaveChi
=FALSE
;
378 int accum
=0, index
, i
,j
,k
,Xi
,n
,b
;
383 char hisfile
[256],histitle
[256], *namept
;
385 int (*calc_bin
)(real
,int,real
);
387 /* Analysis of dihedral transitions */
388 fprintf(stderr
,"Now calculating Chi product trajectories...\n");
395 snew(chi_prtrj
, nframes
) ;
397 /* file for info on all residues */
399 fpall
=xvgropen(fnall
,"Cumulative Rotamers","Residue","Probability",oenv
);
401 fpall
=xvgropen(fnall
,"Cumulative Rotamers","Residue","# Counts",oenv
);
403 for(i
=0; (i
<nlist
); i
++) {
405 /* get nbin, the nr. of cumulative rotamers that need to be considered */
407 for (Xi
= 0 ; Xi
< maxchi
; Xi
++ ) {
408 index
= lookup
[i
][Xi
] ; /* chi_(Xi+1) of res i (-1 if off end) */
414 nbin
+= 1 ; /* for the "zero rotamer", outside the core region */
416 for (j
=0; (j
<nframes
); j
++) {
420 index
= lookup
[i
][0] ; /* index into dih of chi1 of res i */
426 b
= calc_bin(dih
[index
][j
],xity
[index
],core_frac
) ;
430 for (Xi
= 1 ; Xi
< maxchi
; Xi
++ ) {
431 index
= lookup
[i
][Xi
] ; /* chi_(Xi+1) of res i (-1 if off end) */
434 b
= calc_bin(dih
[index
][j
],n
,core_frac
);
435 accum
= n
* accum
+ b
- 1 ;
445 chi_prtrj
[j
] = accum
;
453 /* print cuml rotamer vs time */
454 print_one(oenv
,"chiproduct", dlist
[i
].name
, "chi product for",
455 "cumulative rotamer", nframes
,time
,chi_prtrj
);
458 /* make a histogram pf culm. rotamer occupancy too */
459 snew(chi_prhist
, nbin
) ;
460 make_histo(NULL
,nframes
,chi_prtrj
,nbin
,chi_prhist
,0,nbin
);
462 sprintf(hisfile
,"histo-chiprod%s.xvg",dlist
[i
].name
);
463 sprintf(histitle
,"cumulative rotamer distribution for %s",dlist
[i
].name
);
464 fprintf(stderr
," and %s ",hisfile
);
465 fp
=xvgropen(hisfile
,histitle
,"number","",oenv
);
466 fprintf(fp
,"@ xaxis tick on\n");
467 fprintf(fp
,"@ xaxis tick major 1\n");
468 fprintf(fp
,"@ type xy\n");
469 for(k
=0; (k
<nbin
); k
++) {
471 fprintf(fp
,"%5d %10g\n",k
,(1.0*chi_prhist
[k
])/nframes
);
473 fprintf(fp
,"%5d %10d\n",k
,chi_prhist
[k
]);
479 /* and finally print out occupancies to a single file */
480 /* get the gmx from-1 res nr by setting a ptr to the number part
481 * of dlist[i].name - potential bug for 4-letter res names... */
482 namept
= dlist
[i
].name
+ 3 ;
483 fprintf(fpall
, "%5s ", namept
);
484 for(k
=0; (k
<nbin
); k
++) {
486 fprintf(fpall
," %10g",(1.0*chi_prhist
[k
])/nframes
);
488 fprintf(fpall
," %10d",chi_prhist
[k
]);
490 fprintf(fpall
, "\n") ;
499 fprintf(stderr
,"\n") ;
503 void calc_distribution_props(int nh
,int histo
[],real start
,
504 int nkkk
, t_karplus kkk
[],
507 real d
,dc
,ds
,c1
,c2
,tdc
,tds
;
508 real fac
,ang
,invth
,Jc
;
512 gmx_fatal(FARGS
,"No points in histogram (%s, %d)",__FILE__
,__LINE__
);
515 /* Compute normalisation factor */
517 for(j
=0; (j
<nh
); j
++)
521 for(i
=0; (i
<nkkk
); i
++) {
526 for(j
=0; (j
<nh
); j
++) {
535 for(i
=0; (i
<nkkk
); i
++) {
536 c1
= cos(ang
+kkk
[i
].offset
);
538 Jc
= (kkk
[i
].A
*c2
+ kkk
[i
].B
*c1
+ kkk
[i
].C
);
539 kkk
[i
].Jc
+= histo
[j
]*Jc
;
540 kkk
[i
].Jcsig
+= histo
[j
]*sqr(Jc
);
543 for(i
=0; (i
<nkkk
); i
++) {
545 kkk
[i
].Jcsig
= sqrt(kkk
[i
].Jcsig
/th
-sqr(kkk
[i
].Jc
));
547 *S2
= tdc
*tdc
+tds
*tds
;
550 static void calc_angles(FILE *log
,t_pbc
*pbc
,
551 int n3
,atom_id index
[],real ang
[],rvec x_s
[])
557 for(i
=ix
=0; (ix
<n3
); i
++,ix
+=3)
558 ang
[i
]=bond_angle(x_s
[index
[ix
]],x_s
[index
[ix
+1]],x_s
[index
[ix
+2]],
559 pbc
,r_ij
,r_kj
,&costh
,&t1
,&t2
);
561 fprintf(debug
,"Angle[0]=%g, costh=%g, index0 = %d, %d, %d\n",
562 ang
[0],costh
,index
[0],index
[1],index
[2]);
563 pr_rvec(debug
,0,"rij",r_ij
,DIM
,TRUE
);
564 pr_rvec(debug
,0,"rkj",r_kj
,DIM
,TRUE
);
568 static real
calc_fraction(real angles
[], int nangles
)
571 real trans
= 0, gauche
= 0;
574 for (i
= 0; i
< nangles
; i
++)
576 angle
= angles
[i
] * RAD2DEG
;
578 if (angle
> 135 && angle
< 225)
580 else if (angle
> 270 && angle
< 330)
582 else if (angle
< 90 && angle
> 30)
585 if (trans
+gauche
> 0)
586 return trans
/(trans
+gauche
);
591 static void calc_dihs(FILE *log
,t_pbc
*pbc
,
592 int n4
,atom_id index
[],real ang
[],rvec x_s
[])
595 rvec r_ij
,r_kj
,r_kl
,m
,n
;
596 real cos_phi
,sign
,aaa
;
598 for(i
=ix
=0; (ix
<n4
); i
++,ix
+=4) {
599 aaa
=dih_angle(x_s
[index
[ix
]],x_s
[index
[ix
+1]],x_s
[index
[ix
+2]],
600 x_s
[index
[ix
+3]],pbc
,
602 &cos_phi
,&sign
,&t1
,&t2
,&t3
);
604 ang
[i
]=aaa
; /* not taking into account ryckaert bellemans yet */
608 void make_histo(FILE *log
,
609 int ndata
,real data
[],int npoints
,int histo
[],
617 for(i
=1; (i
<ndata
); i
++) {
618 minx
=min(minx
,data
[i
]);
619 maxx
=max(maxx
,data
[i
]);
621 fprintf(log
,"Min data: %10g Max data: %10g\n",minx
,maxx
);
623 dx
=(double)npoints
/(maxx
-minx
);
626 "Histogramming: ndata=%d, nhisto=%d, minx=%g,maxx=%g,dx=%g\n",
627 ndata
,npoints
,minx
,maxx
,dx
);
628 for(i
=0; (i
<ndata
); i
++) {
629 ind
=(data
[i
]-minx
)*dx
;
630 if ((ind
>= 0) && (ind
< npoints
))
633 fprintf(log
,"index = %d, data[%d] = %g\n",ind
,i
,data
[i
]);
637 void normalize_histo(int npoints
,int histo
[],real dx
,real normhisto
[])
643 for(i
=0; (i
<npoints
); i
++)
646 fprintf(stderr
,"Empty histogram!\n");
650 for(i
=0; (i
<npoints
); i
++)
651 normhisto
[i
]=fac
*histo
[i
];
654 void read_ang_dih(const char *trj_fn
,
655 bool bAngles
,bool bSaveAll
,bool bRb
,bool bPBC
,
656 int maxangstat
,int angstat
[],
657 int *nframes
,real
**time
,
658 int isize
,atom_id index
[],
662 const output_env_t oenv
)
665 int i
,angind
,status
,natoms
,total
,teller
;
667 real t
,fraction
,pifac
,aa
,angle
;
675 natoms
= read_first_x(oenv
,&status
,trj_fn
,&t
,&x
,box
);
685 snew(angles
[cur
],nangles
);
686 snew(angles
[prev
],nangles
);
688 /* Start the loop over frames */
697 if (teller
>= n_alloc
) {
700 for (i
=0; (i
<nangles
); i
++)
701 srenew(dih
[i
],n_alloc
);
702 srenew(*time
,n_alloc
);
703 srenew(*trans_frac
,n_alloc
);
704 srenew(*aver_angle
,n_alloc
);
713 calc_angles(stdout
,pbc
,isize
,index
,angles
[cur
],x
);
716 calc_dihs(stdout
,pbc
,isize
,index
,angles
[cur
],x
);
719 fraction
= calc_fraction(angles
[cur
], nangles
);
720 (*trans_frac
)[teller
] = fraction
;
722 /* Change Ryckaert-Bellemans dihedrals to polymer convention
723 * Modified 990913 by Erik:
724 * We actually shouldn't change the convention, since it's
725 * calculated from polymer above, but we change the intervall
726 * from [-180,180] to [0,360].
729 for(i
=0; (i
<nangles
); i
++)
730 if (angles
[cur
][i
] <= 0.0)
731 angles
[cur
][i
] += 2*M_PI
;
734 /* Periodicity in dihedral space... */
736 for(i
=0; (i
<nangles
); i
++) {
737 real dd
= angles
[cur
][i
];
738 angles
[cur
][i
] = atan2(sin(dd
),cos(dd
));
743 for(i
=0; (i
<nangles
); i
++) {
744 while (angles
[cur
][i
] <= angles
[prev
][i
] - M_PI
)
745 angles
[cur
][i
]+=2*M_PI
;
746 while (angles
[cur
][i
] > angles
[prev
][i
] + M_PI
)
747 angles
[cur
][i
]-=2*M_PI
;
755 for(i
=0; (i
<nangles
); i
++) {
756 aa
=aa
+angles
[cur
][i
];
758 /* angle in rad / 2Pi * max determines bin. bins go from 0 to maxangstat,
759 even though scale goes from -pi to pi (dihedral) or -pi/2 to pi/2
760 (angle) Basically: translate the x-axis by Pi. Translate it back by
764 angle
= angles
[cur
][i
];
766 while (angle
< -M_PI
)
768 while (angle
>= M_PI
)
774 /* Update the distribution histogram */
775 angind
= (int) ((angle
*maxangstat
)/pifac
+ 0.5);
776 if (angind
==maxangstat
)
778 if ( (angind
< 0) || (angind
>= maxangstat
) )
779 /* this will never happen */
780 gmx_fatal(FARGS
,"angle (%f) index out of range (0..%d) : %d\n",
781 angle
,maxangstat
,angind
);
784 if (angind
==maxangstat
)
785 fprintf(stderr
,"angle %d fr %d = %g\n",i
,cur
,angle
);
790 /* average over all angles */
791 (*aver_angle
)[teller
] = (aa
/nangles
);
793 /* this copies all current dih. angles to dih[i], teller is frame */
795 for (i
= 0; i
< nangles
; i
++)
796 dih
[i
][teller
] = angles
[cur
][i
];
801 /* Increment loop counter */
803 } while (read_next_x(oenv
,status
,&t
,natoms
,x
,box
));