4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Green Red Orange Magenta Azure Cyan Skyblue
50 void print_one(char *base
,char *name
,char *title
, char *ylabel
,
51 int nf
,real time
[],real data
[])
54 char buf
[256],t2
[256];
57 sprintf(buf
,"%s%s.xvg",base
,name
);
58 fprintf(stderr
,"\rPrinting %s ",buf
);
59 sprintf(t2
,"%s %s",title
,name
);
60 fp
=xvgropen(buf
,t2
,"Time (ps)",ylabel
);
62 fprintf(fp
,"%10g %10g\n",time
[k
],data
[k
]);
66 static int calc_RBbin(real phi
, int multiplicity
, real core_frac
)
68 /* multiplicity and core_frac NOT used,
69 * just given to enable use of pt-to-fn in caller low_ana_dih_trans*/
70 static const real r30
= M_PI
/6.0;
71 static const real r90
= M_PI
/2.0;
72 static const real r150
= M_PI
*5.0/6.0;
74 if ((phi
< r30
) && (phi
> -r30
))
76 else if ((phi
> -r150
) && (phi
< -r90
))
78 else if ((phi
< r150
) && (phi
> r90
))
83 static int calc_Nbin(real phi
, int multiplicity
, real core_frac
)
85 static const real r360
= 360*DEG2RAD
;
86 real rot_width
, core_width
, core_offset
, low
, hi
;
88 /* with multiplicity 3 and core_frac 0.5
89 * 0<g(-)<120, 120<t<240, 240<g(+)<360
90 * 0< bin0 < 30, 30<bin1<90, 90<bin0<150, 150<bin2<210, 210<bin0<270, 270<bin3<330, 330<bin0<360
91 * so with multiplicity 3, bin1 is core g(-), bin2 is core t, bin3 is
92 core g(+), bin0 is between rotamers */
96 rot_width
= 360/multiplicity
;
97 core_width
= core_frac
* rot_width
;
98 core_offset
= (rot_width
- core_width
)/2.0 ;
99 for(bin
= 1 ; bin
<= multiplicity
; bin
++ ) {
100 low
= ((bin
- 1) * rot_width
) + core_offset
;
101 hi
= ((bin
- 1) * rot_width
) + core_offset
+ core_width
;
104 if ((phi
> low
) && (phi
< hi
))
110 void ana_dih_trans(char *fn_trans
,char *fn_histo
,
111 real
**dih
,int nframes
,int nangles
,
112 char *grpname
,real t0
,real dt
,bool bRb
){
114 /* just a wrapper; declare extra args, then chuck away at end. */
118 int nlist
= nangles
;
123 for(k
=0; (k
<nangles
); k
++) {
127 low_ana_dih_trans(TRUE
, fn_trans
,TRUE
, fn_histo
, maxchi
,
128 dih
, nlist
, dlist
, nframes
,
129 nangles
, grpname
, xity
, t0
, dt
, bRb
, 0.5);
135 void low_ana_dih_trans(bool bTrans
, char *fn_trans
,
136 bool bHisto
, char *fn_histo
, int maxchi
,
137 real
**dih
, int nlist
, t_dlist dlist
[], int nframes
,
138 int nangles
, char *grpname
, int xity
[],
139 real t0
, real dt
, bool bRb
, real core_frac
)
144 int i
,j
,k
,Dih
,ntrans
;
146 real ttime
,tt
,mind
, maxd
, prev
;
147 real
*rot_occ
[NROT
] ;
148 int (*calc_bin
)(real
,int,real
);
150 /* Analysis of dihedral transitions */
151 fprintf(stderr
,"Now calculating transitions...\n");
158 for(k
=0;k
<NROT
;k
++) {
159 snew(rot_occ
[k
],nangles
);
160 for (i
=0; (i
<nangles
); i
++)
166 /* dih[i][j] is the dihedral angle i in frame j */
168 for (i
=0; (i
<nangles
); i
++)
173 mind
= maxd
= prev
= dih
[i
][0];
175 cur_bin
= calc_bin(dih
[i
][0],xity
[i
],core_frac
);
176 rot_occ
[cur_bin
][i
]++ ;
178 for (j
=1; (j
<nframes
); j
++)
180 new_bin
= calc_bin(dih
[i
][j
],xity
[i
],core_frac
);
181 rot_occ
[new_bin
][i
]++ ;
185 else if ((new_bin
!= 0) && (cur_bin
!= new_bin
)) {
192 /* why is all this md rubbish periodic? Remove 360 degree periodicity */
193 if ( (dih
[i
][j
] - prev
) > M_PI
)
195 else if ( (dih
[i
][j
] - prev
) < -M_PI
)
200 mind
= min(mind
, dih
[i
][j
]);
201 maxd
= max(maxd
, dih
[i
][j
]);
202 if ( (maxd
- mind
) > 2*M_PI
/3) /* or 120 degrees, assuming */
203 { /* multiplicity 3. Not so general.*/
206 maxd
= mind
= dih
[i
][j
]; /* get ready for next transition */
212 rot_occ
[k
][i
] /= nframes
;
214 fprintf(stderr
,"Total number of transitions: %10d\n",ntrans
);
216 ttime
= (dt
*nframes
*nangles
)/ntrans
;
217 fprintf(stderr
,"Time between transitions: %10.3f ps\n",ttime
);
220 /* new by grs - copy transitions from tr_h[] to dlist->ntr[]
221 * and rotamer populations from rot_occ to dlist->rot_occ[]
222 * based on fn histogramming in g_chi. diff roles for i and j here */
225 for (Dih
=0; (Dih
<NONCHI
+maxchi
); Dih
++) {
226 for(i
=0; (i
<nlist
); i
++) {
227 if (((Dih
< edOmega
) ) ||
228 ((Dih
== edOmega
) && (has_dihedral(edOmega
,&(dlist
[i
])))) ||
229 ((Dih
> edOmega
) && (dlist
[i
].atm
.Cn
[Dih
-NONCHI
+3] != -1))) {
230 /* grs debug printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
231 dlist
[i
].ntr
[Dih
] = tr_h
[j
] ;
233 dlist
[i
].rot_occ
[Dih
][k
] = rot_occ
[k
][j
] ;
239 /* end addition by grs */
242 sprintf(title
,"Number of transitions: %s",grpname
);
243 fp
=xvgropen(fn_trans
,title
,"Time (ps)","# transitions/timeframe");
244 for(j
=0; (j
<nframes
); j
++) {
246 fprintf(fp
,"%10.3f %10d\n",tt
,tr_f
[j
]);
251 /* Compute histogram from # transitions per dihedral */
253 for(j
=0; (j
<nframes
); j
++)
255 for(i
=0; (i
<nangles
); i
++)
257 for(j
=nframes
; ((tr_f
[j
-1] == 0) && (j
>0)); j
--)
262 sprintf(title
,"Transition time: %s",grpname
);
263 fp
=xvgropen(fn_histo
,title
,"Time (ps)","#");
264 for(i
=j
-1; (i
>0); i
--) {
266 fprintf(fp
,"%10.3f %10d\n",ttime
/i
,tr_f
[i
]);
278 void mk_multiplicity_lookup (int *xity
, int maxchi
, real
**dih
,
279 int nlist
, t_dlist dlist
[],int nangles
)
281 /* new by grs - for dihedral j (as in dih[j]) get multiplicity from dlist
282 * and store in xity[j]
289 for (Dih
=0; (Dih
<NONCHI
+maxchi
); Dih
++) {
290 for(i
=0; (i
<nlist
); i
++) {
291 strncpy(name
, dlist
[i
].name
,3) ;
293 if (((Dih
< edOmega
) ) ||
294 ((Dih
== edOmega
) && (has_dihedral(edOmega
,&(dlist
[i
])))) ||
295 ((Dih
> edOmega
) && (dlist
[i
].atm
.Cn
[Dih
-NONCHI
+3] != -1))) {
296 /* default - we will correct the rest below */
299 /* make omegas 2fold, though doesn't make much more sense than 3 */
300 if (Dih
== edOmega
&& (has_dihedral(edOmega
,&(dlist
[i
])))) {
304 /* dihedrals to aromatic rings, COO, CONH2 or guanidinium are 2fold*/
305 if (Dih
> edOmega
&& (dlist
[i
].atm
.Cn
[Dih
-NONCHI
+3] != -1)) {
306 if ( ((strstr(name
,"PHE") != NULL
) && (Dih
== edChi2
)) ||
307 ((strstr(name
,"TYR") != NULL
) && (Dih
== edChi2
)) ||
308 ((strstr(name
,"PTR") != NULL
) && (Dih
== edChi2
)) ||
309 ((strstr(name
,"TRP") != NULL
) && (Dih
== edChi2
)) ||
310 ((strstr(name
,"HIS") != NULL
) && (Dih
== edChi2
)) ||
311 ((strstr(name
,"GLU") != NULL
) && (Dih
== edChi3
)) ||
312 ((strstr(name
,"ASP") != NULL
) && (Dih
== edChi2
)) ||
313 ((strstr(name
,"GLN") != NULL
) && (Dih
== edChi3
)) ||
314 ((strstr(name
,"ASN") != NULL
) && (Dih
== edChi2
)) ||
315 ((strstr(name
,"ARG") != NULL
) && (Dih
== edChi4
)) ) {
324 fprintf(stderr
,"WARNING: not all dihedrals found in topology (only %d out of %d)!\n",
326 /* Check for remaining dihedrals */
327 for(;(j
< nangles
); j
++)
332 void mk_chi_lookup (int **lookup
, int maxchi
, real
**dih
,
333 int nlist
, t_dlist dlist
[])
336 /* by grs. should rewrite everything to use this. (but haven't,
337 * and at mmt only used in get_chi_product_traj
338 * returns the dihed number given the residue number (from-0)
339 * and chi (from-0) nr. -1 for chi undefined for that res (eg gly, ala..)*/
344 for (Dih
=0; (Dih
<NONCHI
+maxchi
); Dih
++) {
345 for(i
=0; (i
<nlist
); i
++) {
347 if (((Dih
< edOmega
) ) ||
348 ((Dih
== edOmega
) && (has_dihedral(edOmega
,&(dlist
[i
])))) ||
349 ((Dih
> edOmega
) && (dlist
[i
].atm
.Cn
[Dih
-NONCHI
+3] != -1))) {
350 /* grs debug printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
351 if (Dih
> edOmega
) {
356 lookup
[i
][Chi
] = -1 ;
364 void get_chi_product_traj (real
**dih
,int nframes
,int nangles
, int nlist
,
365 int maxchi
, t_dlist dlist
[], real time
[],
366 int **lookup
, int *xity
,bool bRb
, bool bNormalize
,
367 real core_frac
, bool bAll
, char *fnall
)
370 bool bRotZero
, bHaveChi
=FALSE
;
371 int accum
=0, index
, i
,j
,k
,Xi
,n
,b
;
376 char hisfile
[256],histitle
[256], *namept
;
378 int (*calc_bin
)(real
,int,real
);
380 /* Analysis of dihedral transitions */
381 fprintf(stderr
,"Now calculating Chi product trajectories...\n");
388 snew(chi_prtrj
, nframes
) ;
390 /* file for info on all residues */
392 fpall
=xvgropen(fnall
,"Cumulative Rotamers","Residue","Probability");
394 fpall
=xvgropen(fnall
,"Cumulative Rotamers","Residue","# Counts");
396 for(i
=0; (i
<nlist
); i
++) {
398 /* get nbin, the nr. of cumulative rotamers that need to be considered */
400 for (Xi
= 0 ; Xi
< maxchi
; Xi
++ ) {
401 index
= lookup
[i
][Xi
] ; /* chi_(Xi+1) of res i (-1 if off end) */
407 nbin
+= 1 ; /* for the "zero rotamer", outside the core region */
409 for (j
=0; (j
<nframes
); j
++) {
413 index
= lookup
[i
][0] ; /* index into dih of chi1 of res i */
419 b
= calc_bin(dih
[index
][j
],xity
[index
],core_frac
) ;
423 for (Xi
= 1 ; Xi
< maxchi
; Xi
++ ) {
424 index
= lookup
[i
][Xi
] ; /* chi_(Xi+1) of res i (-1 if off end) */
427 b
= calc_bin(dih
[index
][j
],n
,core_frac
);
428 accum
= n
* accum
+ b
- 1 ;
438 chi_prtrj
[j
] = accum
;
446 /* print cuml rotamer vs time */
447 print_one("chiproduct", dlist
[i
].name
, "chi product for",
448 "cumulative rotamer", nframes
,time
,chi_prtrj
);
451 /* make a histogram pf culm. rotamer occupancy too */
452 snew(chi_prhist
, nbin
) ;
453 make_histo(NULL
,nframes
,chi_prtrj
,nbin
,chi_prhist
,0,nbin
);
455 sprintf(hisfile
,"histo-chiprod%s.xvg",dlist
[i
].name
);
456 sprintf(histitle
,"cumulative rotamer distribution for %s",dlist
[i
].name
);
457 fprintf(stderr
," and %s ",hisfile
);
458 fp
=xvgropen(hisfile
,histitle
,"number","");
459 fprintf(fp
,"@ xaxis tick on\n");
460 fprintf(fp
,"@ xaxis tick major 1\n");
461 fprintf(fp
,"@ type xy\n");
462 for(k
=0; (k
<nbin
); k
++) {
464 fprintf(fp
,"%5d %10g\n",k
,(1.0*chi_prhist
[k
])/nframes
);
466 fprintf(fp
,"%5d %10d\n",k
,chi_prhist
[k
]);
472 /* and finally print out occupancies to a single file */
473 /* get the gmx from-1 res nr by setting a ptr to the number part
474 * of dlist[i].name - potential bug for 4-letter res names... */
475 namept
= dlist
[i
].name
+ 3 ;
476 fprintf(fpall
, "%5s ", namept
);
477 for(k
=0; (k
<nbin
); k
++) {
479 fprintf(fpall
," %10g",(1.0*chi_prhist
[k
])/nframes
);
481 fprintf(fpall
," %10d",chi_prhist
[k
]);
483 fprintf(fpall
, "\n") ;
492 fprintf(stderr
,"\n") ;
496 void calc_distribution_props(int nh
,int histo
[],real start
,
497 int nkkk
, t_karplus kkk
[],
500 real d
,dc
,ds
,c1
,c2
,tdc
,tds
;
501 real fac
,ang
,invth
,Jc
;
505 fatal_error(0,"No points in histogram (%s, %d)",__FILE__
,__LINE__
);
508 /* Compute normalisation factor */
510 for(j
=0; (j
<nh
); j
++)
514 for(i
=0; (i
<nkkk
); i
++) {
519 for(j
=0; (j
<nh
); j
++) {
528 for(i
=0; (i
<nkkk
); i
++) {
529 c1
= cos(ang
+kkk
[i
].offset
);
531 Jc
= (kkk
[i
].A
*c2
+ kkk
[i
].B
*c1
+ kkk
[i
].C
);
532 kkk
[i
].Jc
+= histo
[j
]*Jc
;
533 kkk
[i
].Jcsig
+= histo
[j
]*sqr(Jc
);
536 for(i
=0; (i
<nkkk
); i
++) {
538 kkk
[i
].Jcsig
= sqrt(kkk
[i
].Jcsig
/th
-sqr(kkk
[i
].Jc
));
540 *S2
= tdc
*tdc
+tds
*tds
;
543 static void calc_angles(FILE *log
,matrix box
,
544 int n3
,atom_id index
[],real ang
[],rvec x_s
[])
550 for(i
=ix
=0; (ix
<n3
); i
++,ix
+=3)
551 ang
[i
]=bond_angle(box
,x_s
[index
[ix
]],x_s
[index
[ix
+1]],x_s
[index
[ix
+2]],
552 r_ij
,r_kj
,&costh
,&t1
,&t2
);
554 fprintf(debug
,"Angle[0]=%g, costh=%g, index0 = %d, %d, %d\n",
555 ang
[0],costh
,index
[0],index
[1],index
[2]);
556 pr_rvec(debug
,0,"rij",r_ij
,DIM
,TRUE
);
557 pr_rvec(debug
,0,"rkj",r_kj
,DIM
,TRUE
);
558 pr_rvecs(debug
,0,"box",box
,DIM
);
562 static real
calc_fraction(real angles
[], int nangles
)
565 real trans
= 0, gauche
= 0;
568 for (i
= 0; i
< nangles
; i
++)
570 angle
= angles
[i
] * RAD2DEG
;
572 if (angle
> 135 && angle
< 225)
574 else if (angle
> 270 && angle
< 330)
576 else if (angle
< 90 && angle
> 30)
579 if (trans
+gauche
> 0)
580 return trans
/(trans
+gauche
);
585 static void calc_dihs(FILE *log
,matrix box
,
586 int n4
,atom_id index
[],real ang
[],rvec x_s
[])
589 rvec r_ij
,r_kj
,r_kl
,m
,n
;
590 real cos_phi
,sign
,aaa
;
592 for(i
=ix
=0; (ix
<n4
); i
++,ix
+=4) {
594 x_s
[index
[ix
]],x_s
[index
[ix
+1]],x_s
[index
[ix
+2]],
597 &cos_phi
,&sign
,&t1
,&t2
,&t3
);
598 ang
[i
]=aaa
; /* not taking into account ryckaert bellemans yet */
602 void make_histo(FILE *log
,
603 int ndata
,real data
[],int npoints
,int histo
[],
611 for(i
=1; (i
<ndata
); i
++) {
612 minx
=min(minx
,data
[i
]);
613 maxx
=max(maxx
,data
[i
]);
615 fprintf(log
,"Min data: %10g Max data: %10g\n",minx
,maxx
);
617 dx
=(double)npoints
/(maxx
-minx
);
620 "Histogramming: ndata=%d, nhisto=%d, minx=%g,maxx=%g,dx=%g\n",
621 ndata
,npoints
,minx
,maxx
,dx
);
622 for(i
=0; (i
<ndata
); i
++) {
623 ind
=(data
[i
]-minx
)*dx
;
624 if ((ind
>= 0) && (ind
< npoints
))
627 fprintf(log
,"index = %d, data[%d] = %g\n",ind
,i
,data
[i
]);
631 void normalize_histo(int npoints
,int histo
[],real dx
,real normhisto
[])
637 for(i
=0; (i
<npoints
); i
++)
640 fprintf(stderr
,"Empty histogram!\n");
644 for(i
=0; (i
<npoints
); i
++)
645 normhisto
[i
]=fac
*histo
[i
];
648 void read_ang_dih(char *trj_fn
,char *stx_fn
,
649 bool bAngles
,bool bSaveAll
,bool bRb
,
650 int maxangstat
,int angstat
[],
651 int *nframes
,real
**time
,
652 int isize
,atom_id index
[],
658 int ftp
,i
,angind
,status
,natoms
,nat
,total
,teller
;
659 int nangles
,nat_trj
,n_alloc
;
660 real t
,fraction
,pifac
,aa
,angle
;
669 ftp
= fn2ftp(stx_fn
);
670 if ((ftp
== efTPR
) || (ftp
== efTPB
) || (ftp
== efTPA
)) {
671 top
= read_top(stx_fn
);
672 natoms
= top
->atoms
.nr
;
676 get_stx_coordnum(stx_fn
,&natoms
);
677 fprintf(stderr
,"Can not remove periodicity since %s is not a GROMACS topology\n",stx_fn
);
679 nat_trj
= read_first_x(&status
,trj_fn
,&t
,&x
,box
);
681 /* Check for consistency of topology and trajectory */
682 if (natoms
< nat_trj
)
683 fprintf(stderr
,"WARNING! Topology has fewer atoms than trajectory\n");
698 snew(angles
[cur
],nangles
);
699 snew(angles
[prev
],nangles
);
701 /* Start the loop over frames */
710 if (teller
>= n_alloc
) {
713 for (i
=0; (i
<nangles
); i
++)
714 srenew(dih
[i
],n_alloc
);
715 srenew(*time
,n_alloc
);
716 srenew(*trans_frac
,n_alloc
);
717 srenew(*aver_angle
,n_alloc
);
723 rm_pbc(&(top
->idef
),nat_trj
,box
,x
,x_s
);
726 calc_angles(stdout
,box
,isize
,index
,angles
[cur
],x_s
);
728 calc_dihs(stdout
,box
,isize
,index
,angles
[cur
],x_s
);
731 fraction
= calc_fraction(angles
[cur
], nangles
);
732 (*trans_frac
)[teller
] = fraction
;
734 /* Change Ryckaert-Bellemans dihedrals to polymer convention
735 * Modified 990913 by Erik:
736 * We actually shouldn't change the convention, since it's
737 * calculated from polymer above, but we change the intervall
738 * from [-180,180] to [0,360].
741 for(i
=0; (i
<nangles
); i
++)
742 if (angles
[cur
][i
] <= 0.0)
743 angles
[cur
][i
] += 2*M_PI
;
746 /* Periodicity in dihedral space... */
748 for(i
=0; (i
<nangles
); i
++) {
749 while (angles
[cur
][i
] <= angles
[prev
][i
] - M_PI
)
750 angles
[cur
][i
]+=2*M_PI
;
751 while (angles
[cur
][i
] > angles
[prev
][i
] + M_PI
)
752 angles
[cur
][i
]-=2*M_PI
;
759 for(i
=0; (i
<nangles
); i
++) {
760 aa
=aa
+angles
[cur
][i
];
762 /* angle in rad / 2Pi * max determines bin. bins go from 0 to maxangstat,
763 even though scale goes from -pi to pi (dihedral) or -pi/2 to pi/2
764 (angle) Basically: translate the x-axis by Pi. Translate it back by
768 angle
= angles
[cur
][i
];
770 while (angle
< -M_PI
)
772 while (angle
>= M_PI
)
778 /* Update the distribution histogram */
779 angind
= (int) ((angle
*maxangstat
)/pifac
+ 0.5);
780 if (angind
==maxangstat
)
782 if ( (angind
< 0) || (angind
>= maxangstat
) )
783 /* this will never happen */
784 fatal_error(0,"angle (%f) index out of range (0..%d) : %d\n",
785 angle
,maxangstat
,angind
);
788 if (angind
==maxangstat
)
789 fprintf(stderr
,"angle %d fr %d = %g\n",i
,cur
,angle
);
794 /* average over all angles */
795 (*aver_angle
)[teller
] = (aa
/nangles
);
797 /* this copies all current dih. angles to dih[i], teller is frame */
799 for (i
= 0; i
< nangles
; i
++)
800 dih
[i
][teller
] = angles
[cur
][i
];
805 /* Increment loop counter */
807 } while (read_next_x(status
,&t
,nat_trj
,x
,box
));