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 *time,gmx_bool bRb,
117 const output_env_t oenv)
119 /* just a wrapper; declare extra args, then chuck away at end. */
123 int nlist = nangles ;
127 snew(multiplicity,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, multiplicity, time, bRb, 0.5,oenv);
140 void low_ana_dih_trans(gmx_bool bTrans, const char *fn_trans,
141 gmx_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 multiplicity[],
144 real *time, gmx_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);
160 /* Assumes the frames are equally spaced in time */
161 dt=(time[nframes-1]-time[0])/(nframes-1);
163 /* Analysis of dihedral transitions */
164 fprintf(stderr,"Now calculating transitions...\n");
171 for(k=0;k<NROT;k++) {
172 snew(rot_occ[k],nangles);
173 for (i=0; (i<nangles); i++)
179 /* dih[i][j] is the dihedral angle i in frame j */
181 for (i=0; (i<nangles); i++)
186 mind = maxd = prev = dih[i][0];
188 cur_bin = calc_bin(dih[i][0],multiplicity[i],core_frac);
189 rot_occ[cur_bin][i]++ ;
191 for (j=1; (j<nframes); j++)
193 new_bin = calc_bin(dih[i][j],multiplicity[i],core_frac);
194 rot_occ[new_bin][i]++ ;
198 else if ((new_bin != 0) && (cur_bin != new_bin)) {
205 /* why is all this md rubbish periodic? Remove 360 degree periodicity */
206 if ( (dih[i][j] - prev) > M_PI)
208 else if ( (dih[i][j] - prev) < -M_PI)
213 mind = min(mind, dih[i][j]);
214 maxd = max(maxd, dih[i][j]);
215 if ( (maxd - mind) > 2*M_PI/3) /* or 120 degrees, assuming */
216 { /* multiplicity 3. Not so general.*/
219 maxd = mind = dih[i][j]; /* get ready for next transition */
225 rot_occ[k][i] /= nframes ;
227 fprintf(stderr,"Total number of transitions: %10d\n",ntrans);
229 ttime = (dt*nframes*nangles)/ntrans;
230 fprintf(stderr,"Time between transitions: %10.3f ps\n",ttime);
233 /* new by grs - copy transitions from tr_h[] to dlist->ntr[]
234 * and rotamer populations from rot_occ to dlist->rot_occ[]
235 * based on fn histogramming in g_chi. diff roles for i and j here */
238 for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
239 for(i=0; (i<nlist); i++) {
240 if (((Dih < edOmega) ) ||
241 ((Dih == edOmega) && (has_dihedral(edOmega,&(dlist[i])))) ||
242 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1))) {
243 /* grs debug printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
244 dlist[i].ntr[Dih] = tr_h[j] ;
246 dlist[i].rot_occ[Dih][k] = rot_occ[k][j] ;
252 /* end addition by grs */
255 sprintf(title,"Number of transitions: %s",grpname);
256 fp=xvgropen(fn_trans,title,"Time (ps)","# transitions/timeframe",oenv);
257 for(j=0; (j<nframes); j++) {
258 fprintf(fp,"%10.3f %10d\n",time[j],tr_f[j]);
263 /* Compute histogram from # transitions per dihedral */
265 for(j=0; (j<nframes); j++)
267 for(i=0; (i<nangles); i++)
269 for(j=nframes; ((tr_f[j-1] == 0) && (j>0)); j--)
274 sprintf(title,"Transition time: %s",grpname);
275 fp=xvgropen(fn_histo,title,"Time (ps)","#",oenv);
276 for(i=j-1; (i>0); i--) {
278 fprintf(fp,"%10.3f %10d\n",ttime/i,tr_f[i]);
290 void mk_multiplicity_lookup (int *multiplicity, int maxchi, real **dih,
291 int nlist, t_dlist dlist[],int nangles)
293 /* new by grs - for dihedral j (as in dih[j]) get multiplicity from dlist
294 * and store in multiplicity[j]
301 for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
302 for(i=0; (i<nlist); i++) {
303 strncpy(name, dlist[i].name,3) ;
305 if (((Dih < edOmega) ) ||
306 ((Dih == edOmega) && (has_dihedral(edOmega,&(dlist[i])))) ||
307 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1))) {
308 /* default - we will correct the rest below */
309 multiplicity[j] = 3 ;
311 /* make omegas 2fold, though doesn't make much more sense than 3 */
312 if (Dih == edOmega && (has_dihedral(edOmega,&(dlist[i])))) {
313 multiplicity[j] = 2 ;
316 /* dihedrals to aromatic rings, COO, CONH2 or guanidinium are 2fold*/
317 if (Dih > edOmega && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)) {
318 if ( ((strstr(name,"PHE") != NULL) && (Dih == edChi2)) ||
319 ((strstr(name,"TYR") != NULL) && (Dih == edChi2)) ||
320 ((strstr(name,"PTR") != NULL) && (Dih == edChi2)) ||
321 ((strstr(name,"TRP") != NULL) && (Dih == edChi2)) ||
322 ((strstr(name,"HIS") != NULL) && (Dih == edChi2)) ||
323 ((strstr(name,"GLU") != NULL) && (Dih == edChi3)) ||
324 ((strstr(name,"ASP") != NULL) && (Dih == edChi2)) ||
325 ((strstr(name,"GLN") != NULL) && (Dih == edChi3)) ||
326 ((strstr(name,"ASN") != NULL) && (Dih == edChi2)) ||
327 ((strstr(name,"ARG") != NULL) && (Dih == edChi4)) ) {
336 fprintf(stderr,"WARNING: not all dihedrals found in topology (only %d out of %d)!\n",
338 /* Check for remaining dihedrals */
339 for(;(j < nangles); j++)
344 void mk_chi_lookup (int **lookup, int maxchi, real **dih,
345 int nlist, t_dlist dlist[])
348 /* by grs. should rewrite everything to use this. (but haven't,
349 * and at mmt only used in get_chi_product_traj
350 * returns the dihed number given the residue number (from-0)
351 * and chi (from-0) nr. -1 for chi undefined for that res (eg gly, ala..)*/
356 for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
357 for(i=0; (i<nlist); i++) {
359 if (((Dih < edOmega) ) ||
360 ((Dih == edOmega) && (has_dihedral(edOmega,&(dlist[i])))) ||
361 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1))) {
362 /* grs debug printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
363 if (Dih > edOmega ) {
368 lookup[i][Chi] = -1 ;
376 void get_chi_product_traj (real **dih,int nframes,int nangles, int nlist,
377 int maxchi, t_dlist dlist[], real time[],
378 int **lookup, int *multiplicity,gmx_bool bRb, gmx_bool bNormalize,
379 real core_frac, gmx_bool bAll, const char *fnall,
380 const output_env_t oenv)
383 gmx_bool bRotZero, bHaveChi=FALSE;
384 int accum=0, index, i,j,k,Xi,n,b ;
389 char hisfile[256],histitle[256], *namept;
391 int (*calc_bin)(real,int,real);
393 /* Analysis of dihedral transitions */
394 fprintf(stderr,"Now calculating Chi product trajectories...\n");
401 snew(chi_prtrj, nframes) ;
403 /* file for info on all residues */
405 fpall=xvgropen(fnall,"Cumulative Rotamers","Residue","Probability",oenv);
407 fpall=xvgropen(fnall,"Cumulative Rotamers","Residue","# Counts",oenv);
409 for(i=0; (i<nlist); i++) {
411 /* get nbin, the nr. of cumulative rotamers that need to be considered */
413 for (Xi = 0 ; Xi < maxchi ; Xi ++ ) {
414 index = lookup[i][Xi] ; /* chi_(Xi+1) of res i (-1 if off end) */
416 n = multiplicity[index];
420 nbin += 1 ; /* for the "zero rotamer", outside the core region */
422 for (j=0; (j<nframes); j++) {
426 index = lookup[i][0] ; /* index into dih of chi1 of res i */
432 b = calc_bin(dih[index][j],multiplicity[index],core_frac) ;
436 for (Xi = 1 ; Xi < maxchi ; Xi ++ ) {
437 index = lookup[i][Xi] ; /* chi_(Xi+1) of res i (-1 if off end) */
439 n = multiplicity[index];
440 b = calc_bin(dih[index][j],n,core_frac);
441 accum = n * accum + b - 1 ;
451 chi_prtrj[j] = accum ;
459 /* print cuml rotamer vs time */
460 print_one(oenv,"chiproduct", dlist[i].name, "chi product for",
461 "cumulative rotamer", nframes,time,chi_prtrj);
464 /* make a histogram pf culm. rotamer occupancy too */
465 snew(chi_prhist, nbin) ;
466 make_histo(NULL,nframes,chi_prtrj,nbin,chi_prhist,0,nbin);
468 sprintf(hisfile,"histo-chiprod%s.xvg",dlist[i].name);
469 sprintf(histitle,"cumulative rotamer distribution for %s",dlist[i].name);
470 fprintf(stderr," and %s ",hisfile);
471 fp=xvgropen(hisfile,histitle,"number","",oenv);
472 fprintf(fp,"@ xaxis tick on\n");
473 fprintf(fp,"@ xaxis tick major 1\n");
474 fprintf(fp,"@ type xy\n");
475 for(k=0; (k<nbin); k++) {
477 fprintf(fp,"%5d %10g\n",k,(1.0*chi_prhist[k])/nframes);
479 fprintf(fp,"%5d %10d\n",k,chi_prhist[k]);
485 /* and finally print out occupancies to a single file */
486 /* get the gmx from-1 res nr by setting a ptr to the number part
487 * of dlist[i].name - potential bug for 4-letter res names... */
488 namept = dlist[i].name + 3 ;
489 fprintf(fpall, "%5s ", namept);
490 for(k=0; (k<nbin); k++) {
492 fprintf(fpall," %10g",(1.0*chi_prhist[k])/nframes);
494 fprintf(fpall," %10d",chi_prhist[k]);
496 fprintf(fpall, "\n") ;
505 fprintf(stderr,"\n") ;
509 void calc_distribution_props(int nh,int histo[],real start,
510 int nkkk, t_karplus kkk[],
513 real d,dc,ds,c1,c2,tdc,tds;
514 real fac,ang,invth,Jc;
518 gmx_fatal(FARGS,"No points in histogram (%s, %d)",__FILE__,__LINE__);
521 /* Compute normalisation factor */
523 for(j=0; (j<nh); j++)
527 for(i=0; (i<nkkk); i++) {
532 for(j=0; (j<nh); j++) {
541 for(i=0; (i<nkkk); i++) {
542 c1 = cos(ang+kkk[i].offset);
544 Jc = (kkk[i].A*c2 + kkk[i].B*c1 + kkk[i].C);
545 kkk[i].Jc += histo[j]*Jc;
546 kkk[i].Jcsig += histo[j]*sqr(Jc);
549 for(i=0; (i<nkkk); i++) {
551 kkk[i].Jcsig = sqrt(kkk[i].Jcsig/th-sqr(kkk[i].Jc));
553 *S2 = tdc*tdc+tds*tds;
556 static void calc_angles(FILE *log,t_pbc *pbc,
557 int n3,atom_id index[],real ang[],rvec x_s[])
563 for(i=ix=0; (ix<n3); i++,ix+=3)
564 ang[i]=bond_angle(x_s[index[ix]],x_s[index[ix+1]],x_s[index[ix+2]],
565 pbc,r_ij,r_kj,&costh,&t1,&t2);
567 fprintf(debug,"Angle[0]=%g, costh=%g, index0 = %d, %d, %d\n",
568 ang[0],costh,index[0],index[1],index[2]);
569 pr_rvec(debug,0,"rij",r_ij,DIM,TRUE);
570 pr_rvec(debug,0,"rkj",r_kj,DIM,TRUE);
574 static real calc_fraction(real angles[], int nangles)
577 real trans = 0, gauche = 0;
580 for (i = 0; i < nangles; i++)
582 angle = angles[i] * RAD2DEG;
584 if (angle > 135 && angle < 225)
586 else if (angle > 270 && angle < 330)
588 else if (angle < 90 && angle > 30)
591 if (trans+gauche > 0)
592 return trans/(trans+gauche);
597 static void calc_dihs(FILE *log,t_pbc *pbc,
598 int n4,atom_id index[],real ang[],rvec x_s[])
601 rvec r_ij,r_kj,r_kl,m,n;
604 for(i=ix=0; (ix<n4); i++,ix+=4) {
605 aaa=dih_angle(x_s[index[ix]],x_s[index[ix+1]],x_s[index[ix+2]],
606 x_s[index[ix+3]],pbc,
610 ang[i]=aaa; /* not taking into account ryckaert bellemans yet */
614 void make_histo(FILE *log,
615 int ndata,real data[],int npoints,int histo[],
623 for(i=1; (i<ndata); i++) {
624 minx=min(minx,data[i]);
625 maxx=max(maxx,data[i]);
627 fprintf(log,"Min data: %10g Max data: %10g\n",minx,maxx);
629 dx=(double)npoints/(maxx-minx);
632 "Histogramming: ndata=%d, nhisto=%d, minx=%g,maxx=%g,dx=%g\n",
633 ndata,npoints,minx,maxx,dx);
634 for(i=0; (i<ndata); i++) {
635 ind=(data[i]-minx)*dx;
636 if ((ind >= 0) && (ind < npoints))
639 fprintf(log,"index = %d, data[%d] = %g\n",ind,i,data[i]);
643 void normalize_histo(int npoints,int histo[],real dx,real normhisto[])
649 for(i=0; (i<npoints); i++)
652 fprintf(stderr,"Empty histogram!\n");
656 for(i=0; (i<npoints); i++)
657 normhisto[i]=fac*histo[i];
660 void read_ang_dih(const char *trj_fn,
661 gmx_bool bAngles,gmx_bool bSaveAll,gmx_bool bRb,gmx_bool bPBC,
662 int maxangstat,int angstat[],
663 int *nframes,real **time,
664 int isize,atom_id index[],
668 const output_env_t oenv)
672 int i,angind,natoms,total,teller;
674 real t,fraction,pifac,aa,angle;
682 natoms = read_first_x(oenv,&status,trj_fn,&t,&x,box);
692 snew(angles[cur],nangles);
693 snew(angles[prev],nangles);
695 /* Start the loop over frames */
704 if (teller >= n_alloc) {
707 for (i=0; (i<nangles); i++)
708 srenew(dih[i],n_alloc);
709 srenew(*time,n_alloc);
710 srenew(*trans_frac,n_alloc);
711 srenew(*aver_angle,n_alloc);
720 calc_angles(stdout,pbc,isize,index,angles[cur],x);
723 calc_dihs(stdout,pbc,isize,index,angles[cur],x);
726 fraction = calc_fraction(angles[cur], nangles);
727 (*trans_frac)[teller] = fraction;
729 /* Change Ryckaert-Bellemans dihedrals to polymer convention
730 * Modified 990913 by Erik:
731 * We actually shouldn't change the convention, since it's
732 * calculated from polymer above, but we change the intervall
733 * from [-180,180] to [0,360].
736 for(i=0; (i<nangles); i++)
737 if (angles[cur][i] <= 0.0)
738 angles[cur][i] += 2*M_PI;
741 /* Periodicity in dihedral space... */
743 for(i=0; (i<nangles); i++) {
744 real dd = angles[cur][i];
745 angles[cur][i] = atan2(sin(dd),cos(dd));
750 for(i=0; (i<nangles); i++) {
751 while (angles[cur][i] <= angles[prev][i] - M_PI)
752 angles[cur][i]+=2*M_PI;
753 while (angles[cur][i] > angles[prev][i] + M_PI)
754 angles[cur][i]-=2*M_PI;
762 for(i=0; (i<nangles); i++) {
763 aa=aa+angles[cur][i];
765 /* angle in rad / 2Pi * max determines bin. bins go from 0 to maxangstat,
766 even though scale goes from -pi to pi (dihedral) or -pi/2 to pi/2
767 (angle) Basically: translate the x-axis by Pi. Translate it back by
771 angle = angles[cur][i];
773 while (angle < -M_PI)
775 while (angle >= M_PI)
781 /* Update the distribution histogram */
782 angind = (int) ((angle*maxangstat)/pifac + 0.5);
783 if (angind==maxangstat)
785 if ( (angind < 0) || (angind >= maxangstat) )
786 /* this will never happen */
787 gmx_fatal(FARGS,"angle (%f) index out of range (0..%d) : %d\n",
788 angle,maxangstat,angind);
791 if (angind==maxangstat)
792 fprintf(stderr,"angle %d fr %d = %g\n",i,cur,angle);
797 /* average over all angles */
798 (*aver_angle)[teller] = (aa/nangles);
800 /* this copies all current dih. angles to dih[i], teller is frame */
802 for (i = 0; i < nangles; i++)
803 dih[i][teller] = angles[cur][i];
808 /* Increment loop counter */
810 } while (read_next_x(oenv,status,&t,natoms,x,box));