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
43 #include "gmx_fatal.h"
63 static gmx_bool bAllowed(real phi, real psi)
65 static const char *map[] = {
66 "1100000000000000001111111000000000001111111111111111111111111",
67 "1100000000000000001111110000000000011111111111111111111111111",
68 "1100000000000000001111110000000000011111111111111111111111111",
69 "1100000000000000001111100000000000111111111111111111111111111",
70 "1100000000000000001111100000000000111111111111111111111111111",
71 "1100000000000000001111100000000001111111111111111111111111111",
72 "1100000000000000001111100000000001111111111111111111111111111",
73 "1100000000000000001111100000000011111111111111111111111111111",
74 "1110000000000000001111110000000111111111111111111111111111111",
75 "1110000000000000001111110000001111111111111111111111111111111",
76 "1110000000000000001111111000011111111111111111111111111111111",
77 "1110000000000000001111111100111111111111111111111111111111111",
78 "1110000000000000001111111111111111111111111111111111111111111",
79 "1110000000000000001111111111111111111111111111111111111111111",
80 "1110000000000000001111111111111111111111111111111111111111111",
81 "1110000000000000001111111111111111111111111111111111111111111",
82 "1110000000000000001111111111111110011111111111111111111111111",
83 "1110000000000000001111111111111100000111111111111111111111111",
84 "1110000000000000001111111111111000000000001111111111111111111",
85 "1100000000000000001111111111110000000000000011111111111111111",
86 "1100000000000000001111111111100000000000000011111111111111111",
87 "1000000000000000001111111111000000000000000001111111111111110",
88 "0000000000000000001111111110000000000000000000111111111111100",
89 "0000000000000000000000000000000000000000000000000000000000000",
90 "0000000000000000000000000000000000000000000000000000000000000",
91 "0000000000000000000000000000000000000000000000000000000000000",
92 "0000000000000000000000000000000000000000000000000000000000000",
93 "0000000000000000000000000000000000000000000000000000000000000",
94 "0000000000000000000000000000000000000000000000000000000000000",
95 "0000000000000000000000000000000000000000000000000000000000000",
96 "0000000000000000000000000000000000000000000000000000000000000",
97 "0000000000000000000000000000000000000000000000000000000000000",
98 "0000000000000000000000000000000000000000000000000000000000000",
99 "0000000000000000000000000000000000000000000000000000000000000",
100 "0000000000000000000000000000000000000000000000000000000000000",
101 "0000000000000000000000000000000000000000000000000000000000000",
102 "0000000000000000000000000000000000000000000000000000000000000",
103 "0000000000000000000000000000000000000000000000000000000000000",
104 "0000000000000000000000000000000000111111111111000000000000000",
105 "1100000000000000000000000000000001111111111111100000000000111",
106 "1100000000000000000000000000000001111111111111110000000000111",
107 "0000000000000000000000000000000000000000000000000000000000000",
108 "0000000000000000000000000000000000000000000000000000000000000",
109 "0000000000000000000000000000000000000000000000000000000000000",
110 "0000000000000000000000000000000000000000000000000000000000000",
111 "0000000000000000000000000000000000000000000000000000000000000",
112 "0000000000000000000000000000000000000000000000000000000000000",
113 "0000000000000000000000000000000000000000000000000000000000000",
114 "0000000000000000000000000000000000000000000000000000000000000",
115 "0000000000000000000000000000000000000000000000000000000000000",
116 "0000000000000000000000000000000000000000000000000000000000000",
117 "0000000000000000000000000000000000000000000000000000000000000",
118 "0000000000000000000000000000000000000000000000000000000000000",
119 "0000000000000000000000000000000000000000000000000000000000000",
120 "0000000000000000000000000000000000000000000000000000000000000",
121 "0000000000000000000000000000000000000000000000000000000000000",
122 "0000000000000000000000000000000000000000000000000000000000000",
123 "0000000000000000000000000000000000000000000000000000000000000",
124 "0000000000000000000000000000000000000000000000000000000000000",
125 "0000000000000000000000000000000000000000000000000000000000000",
126 "0000000000000000000000000000000000000000000000000000000000000"
128 #define NPP asize(map)
131 #define INDEX(ppp) ((((int) (360+ppp*RAD2DEG)) % 360)/6)
135 return (gmx_bool) map[x][y];
138 atom_id *make_chi_ind(int nl, t_dlist dl[], int *ndih)
143 /* There are nl residues with max edMax dihedrals with 4 atoms each */
144 snew(id, nl*edMax*4);
147 for (i = 0; (i < nl); i++)
149 /* Phi, fake the first one */
150 dl[i].j0[edPhi] = n/4;
151 if (dl[i].atm.minC >= 0)
153 id[n++] = dl[i].atm.minC;
157 id[n++] = dl[i].atm.H;
159 id[n++] = dl[i].atm.N;
160 id[n++] = dl[i].atm.Cn[1];
161 id[n++] = dl[i].atm.C;
163 for (i = 0; (i < nl); i++)
165 /* Psi, fake the last one */
166 dl[i].j0[edPsi] = n/4;
167 id[n++] = dl[i].atm.N;
168 id[n++] = dl[i].atm.Cn[1];
169 id[n++] = dl[i].atm.C;
172 id[n++] = dl[i+1].atm.N;
176 id[n++] = dl[i].atm.O;
179 for (i = 1; (i < nl); i++)
182 if (has_dihedral(edOmega, &(dl[i])))
184 dl[i].j0[edOmega] = n/4;
185 id[n++] = dl[i].atm.minCalpha;
186 id[n++] = dl[i].atm.minC;
187 id[n++] = dl[i].atm.N;
188 id[n++] = dl[i].atm.Cn[1];
191 for (Xi = 0; (Xi < MAXCHI); Xi++)
194 for (i = 0; (i < nl); i++)
196 if (dl[i].atm.Cn[Xi+3] != -1)
198 dl[i].j0[edChi1+Xi] = n/4;
199 id[n++] = dl[i].atm.Cn[Xi];
200 id[n++] = dl[i].atm.Cn[Xi+1];
201 id[n++] = dl[i].atm.Cn[Xi+2];
202 id[n++] = dl[i].atm.Cn[Xi+3];
211 int bin(real chi, int mult)
215 return (int) (chi*mult/360.0);
219 static void do_dihcorr(const char *fn, int nf, int ndih, real **dih, real dt,
220 int nlist, t_dlist dlist[], real time[], int maxchi,
221 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bOmega,
222 const output_env_t oenv)
224 char name1[256], name2[256];
227 do_autocorr(fn, oenv, "Dihedral Autocorrelation Function",
228 nf, ndih, dih, dt, eacCos, FALSE);
231 for (i = 0; (i < nlist); i++)
235 print_one(oenv, "corrphi", dlist[i].name, "Phi ACF for", "C(t)", nf/2, time,
240 for (i = 0; (i < nlist); i++)
244 print_one(oenv, "corrpsi", dlist[i].name, "Psi ACF for", "C(t)", nf/2, time,
249 for (i = 0; (i < nlist); i++)
251 if (has_dihedral(edOmega, &dlist[i]))
255 print_one(oenv, "corromega", dlist[i].name, "Omega ACF for", "C(t)",
261 for (Xi = 0; (Xi < maxchi); Xi++)
263 sprintf(name1, "corrchi%d", Xi+1);
264 sprintf(name2, "Chi%d ACF for", Xi+1);
265 for (i = 0; (i < nlist); i++)
267 if (dlist[i].atm.Cn[Xi+3] != -1)
271 print_one(oenv, name1, dlist[i].name, name2, "C(t)", nf/2, time, dih[j]);
277 fprintf(stderr, "\n");
280 static void copy_dih_data(real in[], real out[], int nf, gmx_bool bLEAVE)
282 /* if bLEAVE, do nothing to data in copying to out
283 * otherwise multiply by 180/pi to convert rad to deg */
294 for (i = 0; (i < nf); i++)
300 static void dump_em_all(int nlist, t_dlist dlist[], int nf, real time[],
301 real **dih, int maxchi,
302 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bOmega, gmx_bool bRAD,
303 const output_env_t oenv)
305 char name[256], titlestr[256], ystr[256];
312 strcpy(ystr, "Angle (rad)");
316 strcpy(ystr, "Angle (degrees)");
321 for (i = 0; (i < nlist); i++)
323 /* grs debug printf("OK i %d j %d\n", i, j) ; */
326 copy_dih_data(dih[j], data, nf, bRAD);
327 print_one(oenv, "phi", dlist[i].name, "\\xf\\f{}", ystr, nf, time, data);
331 for (i = 0; (i < nlist); i++)
335 copy_dih_data(dih[j], data, nf, bRAD);
336 print_one(oenv, "psi", dlist[i].name, "\\xy\\f{}", ystr, nf, time, data);
340 for (i = 0; (i < nlist); i++)
342 if (has_dihedral(edOmega, &(dlist[i])))
346 copy_dih_data(dih[j], data, nf, bRAD);
347 print_one(oenv, "omega", dlist[i].name, "\\xw\\f{}", ystr, nf, time, data);
353 for (Xi = 0; (Xi < maxchi); Xi++)
355 for (i = 0; (i < nlist); i++)
357 if (dlist[i].atm.Cn[Xi+3] != -1)
361 sprintf(name, "chi%d", Xi+1);
362 sprintf(titlestr, "\\xc\\f{}\\s%d\\N", Xi+1);
363 copy_dih_data(dih[j], data, nf, bRAD);
364 print_one(oenv, name, dlist[i].name, titlestr, ystr, nf, time, data);
370 fprintf(stderr, "\n");
373 static void reset_one(real dih[], int nf, real phase)
377 for (j = 0; (j < nf); j++)
380 while (dih[j] < -M_PI)
384 while (dih[j] >= M_PI)
391 static int reset_em_all(int nlist, t_dlist dlist[], int nf,
392 real **dih, int maxchi)
399 for (i = 0; (i < nlist); i++)
401 if (dlist[i].atm.minC == -1)
403 reset_one(dih[j++], nf, M_PI);
407 reset_one(dih[j++], nf, 0);
411 for (i = 0; (i < nlist-1); i++)
413 reset_one(dih[j++], nf, 0);
415 /* last Psi is faked from O */
416 reset_one(dih[j++], nf, M_PI);
419 for (i = 0; (i < nlist); i++)
421 if (has_dihedral(edOmega, &dlist[i]))
423 reset_one(dih[j++], nf, 0);
426 /* Chi 1 thru maxchi */
427 for (Xi = 0; (Xi < maxchi); Xi++)
429 for (i = 0; (i < nlist); i++)
431 if (dlist[i].atm.Cn[Xi+3] != -1)
433 reset_one(dih[j], nf, 0);
438 fprintf(stderr, "j after resetting (nr. active dihedrals) = %d\n", j);
442 static void histogramming(FILE *log, int nbin, gmx_residuetype_t rt,
443 int nf, int maxchi, real **dih,
444 int nlist, t_dlist dlist[],
446 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bOmega, gmx_bool bChi,
447 gmx_bool bNormalize, gmx_bool bSSHisto, const char *ssdump,
448 real bfac_max, t_atoms *atoms,
449 gmx_bool bDo_jc, const char *fn,
450 const output_env_t oenv)
452 /* also gets 3J couplings and order parameters S2 */
453 t_karplus kkkphi[] = {
454 { "J_NHa1", 6.51, -1.76, 1.6, -M_PI/3, 0.0, 0.0 },
455 { "J_NHa2", 6.51, -1.76, 1.6, M_PI/3, 0.0, 0.0 },
456 { "J_HaC'", 4.0, 1.1, 0.1, 0.0, 0.0, 0.0 },
457 { "J_NHCb", 4.7, -1.5, -0.2, M_PI/3, 0.0, 0.0 },
458 { "J_Ci-1Hai", 4.5, -1.3, -1.2, 2*M_PI/3, 0.0, 0.0 }
460 t_karplus kkkpsi[] = {
461 { "J_HaN", -0.88, -0.61, -0.27, M_PI/3, 0.0, 0.0 }
463 t_karplus kkkchi1[] = {
464 { "JHaHb2", 9.5, -1.6, 1.8, -M_PI/3, 0, 0.0 },
465 { "JHaHb3", 9.5, -1.6, 1.8, 0, 0, 0.0 }
467 #define NKKKPHI asize(kkkphi)
468 #define NKKKPSI asize(kkkpsi)
469 #define NKKKCHI asize(kkkchi1)
470 #define NJC (NKKKPHI+NKKKPSI+NKKKCHI)
472 FILE *fp, *ssfp[3] = {NULL, NULL, NULL};
473 const char *sss[3] = { "sheet", "helix", "coil" };
477 int ****his_aa_ss = NULL;
478 int ***his_aa, **his_aa1, *histmp;
479 int i, j, k, m, n, nn, Dih, nres, hindex, angle;
480 gmx_bool bBfac, bOccup;
481 char hisfile[256], hhisfile[256], sshisfile[256], title[256], *ss_str = NULL;
483 const char *residue_name;
486 rt_size = gmx_residuetype_get_size(rt);
489 fp = ffopen(ssdump, "r");
490 if (1 != fscanf(fp, "%d", &nres))
492 gmx_fatal(FARGS, "Error reading from file %s", ssdump);
495 snew(ss_str, nres+1);
496 if (1 != fscanf(fp, "%s", ss_str))
498 gmx_fatal(FARGS, "Error reading from file %s", ssdump);
502 /* Four dimensional array... Very cool */
504 for (i = 0; (i < 3); i++)
506 snew(his_aa_ss[i], rt_size+1);
507 for (j = 0; (j <= rt_size); j++)
509 snew(his_aa_ss[i][j], edMax);
510 for (Dih = 0; (Dih < edMax); Dih++)
512 snew(his_aa_ss[i][j][Dih], nbin+1);
518 for (Dih = 0; (Dih < edMax); Dih++)
520 snew(his_aa[Dih], rt_size+1);
521 for (i = 0; (i <= rt_size); i++)
523 snew(his_aa[Dih][i], nbin+1);
530 for (i = 0; (i < nlist); i++)
538 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
540 for (i = 0; (i < nlist); i++)
542 if (((Dih < edOmega) ) ||
543 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
544 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)))
546 make_histo(log, nf, dih[j], nbin, histmp, -M_PI, M_PI);
550 /* Assume there is only one structure, the first.
551 * Compute index in histogram.
553 /* Check the atoms to see whether their B-factors are low enough
554 * Check atoms to see their occupancy is 1.
556 bBfac = bOccup = TRUE;
557 for (nn = 0; (nn < 4); nn++, n++)
559 bBfac = bBfac && (atoms->pdbinfo[index[n]].bfac <= bfac_max);
560 bOccup = bOccup && (atoms->pdbinfo[index[n]].occup == 1);
562 if (bOccup && ((bfac_max <= 0) || ((bfac_max > 0) && bBfac)))
564 hindex = ((dih[j][0]+M_PI)*nbin)/(2*M_PI);
565 range_check(hindex, 0, nbin);
567 /* Assign dihedral to either of the structure determined
570 switch (ss_str[dlist[i].resnr])
573 his_aa_ss[0][dlist[i].index][Dih][hindex]++;
576 his_aa_ss[1][dlist[i].index][Dih][hindex]++;
579 his_aa_ss[2][dlist[i].index][Dih][hindex]++;
585 fprintf(debug, "Res. %d has imcomplete occupancy or bfacs > %g\n",
586 dlist[i].resnr, bfac_max);
597 calc_distribution_props(nbin, histmp, -M_PI, NKKKPHI, kkkphi, &S2);
599 for (m = 0; (m < NKKKPHI); m++)
601 Jc[i][m] = kkkphi[m].Jc;
602 Jcsig[i][m] = kkkphi[m].Jcsig;
606 calc_distribution_props(nbin, histmp, -M_PI, NKKKPSI, kkkpsi, &S2);
608 for (m = 0; (m < NKKKPSI); m++)
610 Jc[i][NKKKPHI+m] = kkkpsi[m].Jc;
611 Jcsig[i][NKKKPHI+m] = kkkpsi[m].Jcsig;
615 calc_distribution_props(nbin, histmp, -M_PI, NKKKCHI, kkkchi1, &S2);
616 for (m = 0; (m < NKKKCHI); m++)
618 Jc[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jc;
619 Jcsig[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jcsig;
622 default: /* covers edOmega and higher Chis than Chi1 */
623 calc_distribution_props(nbin, histmp, -M_PI, 0, NULL, &S2);
626 dlist[i].S2[Dih] = S2;
628 /* Sum distribution per amino acid type as well */
629 for (k = 0; (k < nbin); k++)
631 his_aa[Dih][dlist[i].index][k] += histmp[k];
636 else /* dihed not defined */
638 dlist[i].S2[Dih] = 0.0;
644 /* Print out Jcouplings */
645 fprintf(log, "\n *** J-Couplings from simulation (plus std. dev.) ***\n\n");
646 fprintf(log, "Residue ");
647 for (i = 0; (i < NKKKPHI); i++)
649 fprintf(log, "%7s SD", kkkphi[i].name);
651 for (i = 0; (i < NKKKPSI); i++)
653 fprintf(log, "%7s SD", kkkpsi[i].name);
655 for (i = 0; (i < NKKKCHI); i++)
657 fprintf(log, "%7s SD", kkkchi1[i].name);
660 for (i = 0; (i < NJC+1); i++)
662 fprintf(log, "------------");
665 for (i = 0; (i < nlist); i++)
667 fprintf(log, "%-10s", dlist[i].name);
668 for (j = 0; (j < NJC); j++)
670 fprintf(log, " %5.2f %4.2f", Jc[i][j], Jcsig[i][j]);
676 /* and to -jc file... */
679 fp = xvgropen(fn, "\\S3\\NJ-Couplings from Karplus Equation", "Residue",
682 for (i = 0; (i < NKKKPHI); i++)
684 leg[i] = strdup(kkkphi[i].name);
686 for (i = 0; (i < NKKKPSI); i++)
688 leg[i+NKKKPHI] = strdup(kkkpsi[i].name);
690 for (i = 0; (i < NKKKCHI); i++)
692 leg[i+NKKKPHI+NKKKPSI] = strdup(kkkchi1[i].name);
694 xvgr_legend(fp, NJC, (const char**)leg, oenv);
695 fprintf(fp, "%5s ", "#Res.");
696 for (i = 0; (i < NJC); i++)
698 fprintf(fp, "%10s ", leg[i]);
701 for (i = 0; (i < nlist); i++)
703 fprintf(fp, "%5d ", dlist[i].resnr);
704 for (j = 0; (j < NJC); j++)
706 fprintf(fp, " %8.3f", Jc[i][j]);
711 for (i = 0; (i < NJC); i++)
716 /* finished -jc stuff */
718 snew(normhisto, nbin);
719 for (i = 0; (i < rt_size); i++)
721 for (Dih = 0; (Dih < edMax); Dih++)
723 /* First check whether something is in there */
724 for (j = 0; (j < nbin); j++)
726 if (his_aa[Dih][i][j] != 0)
732 ((bPhi && (Dih == edPhi)) ||
733 (bPsi && (Dih == edPsi)) ||
734 (bOmega && (Dih == edOmega)) ||
735 (bChi && (Dih >= edChi1))))
739 normalize_histo(nbin, his_aa[Dih][i], (360.0/nbin), normhisto);
742 residue_name = gmx_residuetype_get_name(rt, i);
746 sprintf(hisfile, "histo-phi%s", residue_name);
747 sprintf(title, "\\xf\\f{} Distribution for %s", residue_name);
750 sprintf(hisfile, "histo-psi%s", residue_name);
751 sprintf(title, "\\xy\\f{} Distribution for %s", residue_name);
754 sprintf(hisfile, "histo-omega%s", residue_name);
755 sprintf(title, "\\xw\\f{} Distribution for %s", residue_name);
758 sprintf(hisfile, "histo-chi%d%s", Dih-NONCHI+1, residue_name);
759 sprintf(title, "\\xc\\f{}\\s%d\\N Distribution for %s",
760 Dih-NONCHI+1, residue_name);
762 strcpy(hhisfile, hisfile);
763 strcat(hhisfile, ".xvg");
764 fp = xvgropen(hhisfile, title, "Degrees", "", oenv);
765 fprintf(fp, "@ with g0\n");
766 xvgr_world(fp, -180, 0, 180, 0.1, oenv);
767 fprintf(fp, "# this effort to set graph size fails unless you run with -autoscale none or -autoscale y flags\n");
768 fprintf(fp, "@ xaxis tick on\n");
769 fprintf(fp, "@ xaxis tick major 90\n");
770 fprintf(fp, "@ xaxis tick minor 30\n");
771 fprintf(fp, "@ xaxis ticklabel prec 0\n");
772 fprintf(fp, "@ yaxis tick off\n");
773 fprintf(fp, "@ yaxis ticklabel off\n");
774 fprintf(fp, "@ type xy\n");
777 for (k = 0; (k < 3); k++)
779 sprintf(sshisfile, "%s-%s.xvg", hisfile, sss[k]);
780 ssfp[k] = ffopen(sshisfile, "w");
783 for (j = 0; (j < nbin); j++)
785 angle = -180 + (360/nbin)*j;
788 fprintf(fp, "%5d %10g\n", angle, normhisto[j]);
792 fprintf(fp, "%5d %10d\n", angle, his_aa[Dih][i][j]);
796 for (k = 0; (k < 3); k++)
798 fprintf(ssfp[k], "%5d %10d\n", angle,
799 his_aa_ss[k][i][Dih][j]);
807 for (k = 0; (k < 3); k++)
809 fprintf(ssfp[k], "&\n");
820 /* Four dimensional array... Very cool */
821 for (i = 0; (i < 3); i++)
823 for (j = 0; (j <= rt_size); j++)
825 for (Dih = 0; (Dih < edMax); Dih++)
827 sfree(his_aa_ss[i][j][Dih]);
829 sfree(his_aa_ss[i][j]);
838 static FILE *rama_file(const char *fn, const char *title, const char *xaxis,
839 const char *yaxis, const output_env_t oenv)
843 fp = xvgropen(fn, title, xaxis, yaxis, oenv);
844 fprintf(fp, "@ with g0\n");
845 xvgr_world(fp, -180, -180, 180, 180, oenv);
846 fprintf(fp, "@ xaxis tick on\n");
847 fprintf(fp, "@ xaxis tick major 90\n");
848 fprintf(fp, "@ xaxis tick minor 30\n");
849 fprintf(fp, "@ xaxis ticklabel prec 0\n");
850 fprintf(fp, "@ yaxis tick on\n");
851 fprintf(fp, "@ yaxis tick major 90\n");
852 fprintf(fp, "@ yaxis tick minor 30\n");
853 fprintf(fp, "@ yaxis ticklabel prec 0\n");
854 fprintf(fp, "@ s0 type xy\n");
855 fprintf(fp, "@ s0 symbol 2\n");
856 fprintf(fp, "@ s0 symbol size 0.410000\n");
857 fprintf(fp, "@ s0 symbol fill 1\n");
858 fprintf(fp, "@ s0 symbol color 1\n");
859 fprintf(fp, "@ s0 symbol linewidth 1\n");
860 fprintf(fp, "@ s0 symbol linestyle 1\n");
861 fprintf(fp, "@ s0 symbol center false\n");
862 fprintf(fp, "@ s0 symbol char 0\n");
863 fprintf(fp, "@ s0 skip 0\n");
864 fprintf(fp, "@ s0 linestyle 0\n");
865 fprintf(fp, "@ s0 linewidth 1\n");
866 fprintf(fp, "@ type xy\n");
871 static void do_rama(int nf, int nlist, t_dlist dlist[], real **dih,
872 gmx_bool bViol, gmx_bool bRamOmega, const output_env_t oenv)
874 FILE *fp, *gp = NULL;
877 int i, j, k, Xi1, Xi2, Phi, Psi, Om = 0, nlevels;
879 real **mat = NULL, phi, psi, omega, axis[NMAT], lo, hi;
880 t_rgb rlo = { 1.0, 0.0, 0.0 };
881 t_rgb rmid = { 1.0, 1.0, 1.0 };
882 t_rgb rhi = { 0.0, 0.0, 1.0 };
884 for (i = 0; (i < nlist); i++)
886 if ((has_dihedral(edPhi, &(dlist[i]))) &&
887 (has_dihedral(edPsi, &(dlist[i]))))
889 sprintf(fn, "ramaPhiPsi%s.xvg", dlist[i].name);
890 fp = rama_file(fn, "Ramachandran Plot",
891 "\\8f\\4 (deg)", "\\8y\\4 (deg)", oenv);
892 bOm = bRamOmega && has_dihedral(edOmega, &(dlist[i]));
895 Om = dlist[i].j0[edOmega];
897 for (j = 0; (j < NMAT); j++)
900 axis[j] = -180+(360*j)/NMAT;
905 sprintf(fn, "violPhiPsi%s.xvg", dlist[i].name);
906 gp = ffopen(fn, "w");
908 Phi = dlist[i].j0[edPhi];
909 Psi = dlist[i].j0[edPsi];
910 for (j = 0; (j < nf); j++)
912 phi = RAD2DEG*dih[Phi][j];
913 psi = RAD2DEG*dih[Psi][j];
914 fprintf(fp, "%10g %10g\n", phi, psi);
917 fprintf(gp, "%d\n", !bAllowed(dih[Phi][j], RAD2DEG*dih[Psi][j]));
921 omega = RAD2DEG*dih[Om][j];
922 mat[(int)((phi*NMAT)/360)+NMAT/2][(int)((psi*NMAT)/360)+NMAT/2]
933 sprintf(fn, "ramomega%s.xpm", dlist[i].name);
934 fp = ffopen(fn, "w");
936 for (j = 0; (j < NMAT); j++)
938 for (k = 0; (k < NMAT); k++)
941 lo = min(mat[j][k], lo);
942 hi = max(mat[j][k], hi);
946 if (fabs(lo) > fabs(hi))
955 for (j = 0; (j < NMAT); j++)
957 for (k = 0; (k < NMAT); k++)
965 write_xpm3(fp, 0, "Omega/Ramachandran Plot", "Deg", "Phi", "Psi",
966 NMAT, NMAT, axis, axis, mat, lo, 180.0, hi, rlo, rmid, rhi, &nlevels);
968 for (j = 0; (j < NMAT); j++)
975 if ((has_dihedral(edChi1, &(dlist[i]))) &&
976 (has_dihedral(edChi2, &(dlist[i]))))
978 sprintf(fn, "ramaX1X2%s.xvg", dlist[i].name);
979 fp = rama_file(fn, "\\8c\\4\\s1\\N-\\8c\\4\\s2\\N Ramachandran Plot",
980 "\\8c\\4\\s1\\N (deg)", "\\8c\\4\\s2\\N (deg)", oenv);
981 Xi1 = dlist[i].j0[edChi1];
982 Xi2 = dlist[i].j0[edChi2];
983 for (j = 0; (j < nf); j++)
985 fprintf(fp, "%10g %10g\n", RAD2DEG*dih[Xi1][j], RAD2DEG*dih[Xi2][j]);
991 fprintf(stderr, "No chi1 & chi2 angle for %s\n", dlist[i].name);
997 static void print_transitions(const char *fn, int maxchi, int nlist,
998 t_dlist dlist[], t_atoms *atoms, rvec x[],
999 matrix box, gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, real dt,
1000 const output_env_t oenv)
1002 /* based on order_params below */
1007 /* must correspond with enum in pp2shift.h:38 */
1009 #define NLEG asize(leg)
1011 leg[0] = strdup("Phi");
1012 leg[1] = strdup("Psi");
1013 leg[2] = strdup("Omega");
1014 leg[3] = strdup("Chi1");
1015 leg[4] = strdup("Chi2");
1016 leg[5] = strdup("Chi3");
1017 leg[6] = strdup("Chi4");
1018 leg[7] = strdup("Chi5");
1019 leg[8] = strdup("Chi6");
1021 /* Print order parameters */
1022 fp = xvgropen(fn, "Dihedral Rotamer Transitions", "Residue", "Transitions/ns",
1024 xvgr_legend(fp, NONCHI+maxchi, (const char**)leg, oenv);
1026 for (Dih = 0; (Dih < edMax); Dih++)
1031 fprintf(fp, "%5s ", "#Res.");
1032 fprintf(fp, "%10s %10s %10s ", leg[edPhi], leg[edPsi], leg[edOmega]);
1033 for (Xi = 0; Xi < maxchi; Xi++)
1035 fprintf(fp, "%10s ", leg[NONCHI+Xi]);
1039 for (i = 0; (i < nlist); i++)
1041 fprintf(fp, "%5d ", dlist[i].resnr);
1042 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1044 fprintf(fp, "%10.3f ", dlist[i].ntr[Dih]/dt);
1046 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
1052 static void order_params(FILE *log,
1053 const char *fn, int maxchi, int nlist, t_dlist dlist[],
1054 const char *pdbfn, real bfac_init,
1055 t_atoms *atoms, rvec x[], int ePBC, matrix box,
1056 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, const output_env_t oenv)
1064 /* except for S2Min/Max, must correspond with enum in pp2shift.h:38 */
1065 const char *const_leg[2+edMax] = {
1066 "S2Min", "S2Max", "Phi", "Psi", "Omega",
1067 "Chi1", "Chi2", "Chi3", "Chi4", "Chi5",
1070 #define NLEG asize(leg)
1074 for (i = 0; i < NLEG; i++)
1076 leg[i] = strdup(const_leg[i]);
1079 /* Print order parameters */
1080 fp = xvgropen(fn, "Dihedral Order Parameters", "Residue", "S2", oenv);
1081 xvgr_legend(fp, 2+NONCHI+maxchi, const_leg, oenv);
1083 for (Dih = 0; (Dih < edMax); Dih++)
1088 fprintf(fp, "%5s ", "#Res.");
1089 fprintf(fp, "%10s %10s ", leg[0], leg[1]);
1090 fprintf(fp, "%10s %10s %10s ", leg[2+edPhi], leg[2+edPsi], leg[2+edOmega]);
1091 for (Xi = 0; Xi < maxchi; Xi++)
1093 fprintf(fp, "%10s ", leg[2+NONCHI+Xi]);
1097 for (i = 0; (i < nlist); i++)
1101 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1103 if (dlist[i].S2[Dih] != 0)
1105 if (dlist[i].S2[Dih] > S2Max)
1107 S2Max = dlist[i].S2[Dih];
1109 if (dlist[i].S2[Dih] < S2Min)
1111 S2Min = dlist[i].S2[Dih];
1114 if (dlist[i].S2[Dih] > 0.8)
1119 fprintf(fp, "%5d ", dlist[i].resnr);
1120 fprintf(fp, "%10.3f %10.3f ", S2Min, S2Max);
1121 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1123 fprintf(fp, "%10.3f ", dlist[i].S2[Dih]);
1126 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
1134 if (NULL == atoms->pdbinfo)
1136 snew(atoms->pdbinfo, atoms->nr);
1138 for (i = 0; (i < atoms->nr); i++)
1140 atoms->pdbinfo[i].bfac = bfac_init;
1143 for (i = 0; (i < nlist); i++)
1145 atoms->pdbinfo[dlist[i].atm.N].bfac = -dlist[i].S2[0]; /* Phi */
1146 atoms->pdbinfo[dlist[i].atm.H].bfac = -dlist[i].S2[0]; /* Phi */
1147 atoms->pdbinfo[dlist[i].atm.C].bfac = -dlist[i].S2[1]; /* Psi */
1148 atoms->pdbinfo[dlist[i].atm.O].bfac = -dlist[i].S2[1]; /* Psi */
1149 for (Xi = 0; (Xi < maxchi); Xi++) /* Chi's */
1151 if (dlist[i].atm.Cn[Xi+3] != -1)
1153 atoms->pdbinfo[dlist[i].atm.Cn[Xi+1]].bfac = -dlist[i].S2[NONCHI+Xi];
1158 fp = ffopen(pdbfn, "w");
1159 fprintf(fp, "REMARK generated by g_chi\n");
1160 fprintf(fp, "REMARK "
1161 "B-factor field contains negative of dihedral order parameters\n");
1162 write_pdbfile(fp, NULL, atoms, x, ePBC, box, ' ', 0, NULL, TRUE);
1163 x0 = y0 = z0 = 1000.0;
1164 for (i = 0; (i < atoms->nr); i++)
1166 x0 = min(x0, x[i][XX]);
1167 y0 = min(y0, x[i][YY]);
1168 z0 = min(z0, x[i][ZZ]);
1170 x0 *= 10.0; /* nm -> angstrom */
1171 y0 *= 10.0; /* nm -> angstrom */
1172 z0 *= 10.0; /* nm -> angstrom */
1173 sprintf(buf, "%s%%6.f%%6.2f\n", get_pdbformat());
1174 for (i = 0; (i < 10); i++)
1176 fprintf(fp, buf, "ATOM ", atoms->nr+1+i, "CA", "LEG", ' ',
1177 atoms->nres+1, ' ', x0, y0, z0+(1.2*i), 0.0, -0.1*i);
1182 fprintf(log, "Dihedrals with S2 > 0.8\n");
1183 fprintf(log, "Dihedral: ");
1186 fprintf(log, " Phi ");
1190 fprintf(log, " Psi ");
1194 for (Xi = 0; (Xi < maxchi); Xi++)
1196 fprintf(log, " %s ", leg[2+NONCHI+Xi]);
1199 fprintf(log, "\nNumber: ");
1202 fprintf(log, "%4d ", nh[0]);
1206 fprintf(log, "%4d ", nh[1]);
1210 for (Xi = 0; (Xi < maxchi); Xi++)
1212 fprintf(log, "%4d ", nh[NONCHI+Xi]);
1217 for (i = 0; i < NLEG; i++)
1224 int gmx_chi(int argc, char *argv[])
1226 const char *desc[] = {
1227 "[TT]g_chi[tt] computes [GRK]phi[grk], [GRK]psi[grk], [GRK]omega[grk], and [GRK]chi[grk] dihedrals for all your ",
1228 "amino acid backbone and sidechains.",
1229 "It can compute dihedral angle as a function of time, and as",
1230 "histogram distributions.",
1231 "The distributions [TT](histo-(dihedral)(RESIDUE).xvg[tt]) are cumulative over all residues of each type.[PAR]",
1232 "If option [TT]-corr[tt] is given, the program will",
1233 "calculate dihedral autocorrelation functions. The function used",
1234 "is C(t) = [CHEVRON][COS][GRK]chi[grk]([GRK]tau[grk])[cos] [COS][GRK]chi[grk]([GRK]tau[grk]+t)[cos][chevron]. The use of cosines",
1235 "rather than angles themselves, resolves the problem of periodicity.",
1236 "(Van der Spoel & Berendsen (1997), Biophys. J. 72, 2032-2041).",
1237 "Separate files for each dihedral of each residue",
1238 "[TT](corr(dihedral)(RESIDUE)(nresnr).xvg[tt]) are output, as well as a",
1239 "file containing the information for all residues (argument of [TT]-corr[tt]).[PAR]",
1240 "With option [TT]-all[tt], the angles themselves as a function of time for",
1241 "each residue are printed to separate files [TT](dihedral)(RESIDUE)(nresnr).xvg[tt].",
1242 "These can be in radians or degrees.[PAR]",
1243 "A log file (argument [TT]-g[tt]) is also written. This contains [BR]",
1244 "(a) information about the number of residues of each type.[BR]",
1245 "(b) The NMR ^3J coupling constants from the Karplus equation.[BR]",
1246 "(c) a table for each residue of the number of transitions between ",
1247 "rotamers per nanosecond, and the order parameter S^2 of each dihedral.[BR]",
1248 "(d) a table for each residue of the rotamer occupancy.[PAR]",
1249 "All rotamers are taken as 3-fold, except for [GRK]omega[grk] and [GRK]chi[grk] dihedrals",
1250 "to planar groups (i.e. [GRK]chi[grk][SUB]2[sub] of aromatics, Asp and Asn; [GRK]chi[grk][SUB]3[sub] of Glu",
1251 "and Gln; and [GRK]chi[grk][SUB]4[sub] of Arg), which are 2-fold. \"rotamer 0\" means ",
1252 "that the dihedral was not in the core region of each rotamer. ",
1253 "The width of the core region can be set with [TT]-core_rotamer[tt][PAR]",
1255 "The S^2 order parameters are also output to an [TT].xvg[tt] file",
1256 "(argument [TT]-o[tt] ) and optionally as a [TT].pdb[tt] file with",
1257 "the S^2 values as B-factor (argument [TT]-p[tt]). ",
1258 "The total number of rotamer transitions per timestep",
1259 "(argument [TT]-ot[tt]), the number of transitions per rotamer",
1260 "(argument [TT]-rt[tt]), and the ^3J couplings (argument [TT]-jc[tt]), ",
1261 "can also be written to [TT].xvg[tt] files. Note that the analysis",
1262 "of rotamer transitions assumes that the supplied trajectory frames",
1263 "are equally spaced in time.[PAR]",
1265 "If [TT]-chi_prod[tt] is set (and [TT]-maxchi[tt] > 0), cumulative rotamers, e.g.",
1266 "1+9([GRK]chi[grk][SUB]1[sub]-1)+3([GRK]chi[grk][SUB]2[sub]-1)+([GRK]chi[grk][SUB]3[sub]-1) (if the residue has three 3-fold ",
1267 "dihedrals and [TT]-maxchi[tt] >= 3)",
1268 "are calculated. As before, if any dihedral is not in the core region,",
1269 "the rotamer is taken to be 0. The occupancies of these cumulative ",
1270 "rotamers (starting with rotamer 0) are written to the file",
1271 "that is the argument of [TT]-cp[tt], and if the [TT]-all[tt] flag",
1272 "is given, the rotamers as functions of time",
1273 "are written to [TT]chiproduct(RESIDUE)(nresnr).xvg[tt] ",
1274 "and their occupancies to [TT]histo-chiproduct(RESIDUE)(nresnr).xvg[tt].[PAR]",
1276 "The option [TT]-r[tt] generates a contour plot of the average [GRK]omega[grk] angle",
1277 "as a function of the [GRK]phi[grk] and [GRK]psi[grk] angles, that is, in a Ramachandran plot",
1278 "the average [GRK]omega[grk] angle is plotted using color coding.",
1282 const char *bugs[] = {
1283 "Produces MANY output files (up to about 4 times the number of residues in the protein, twice that if autocorrelation functions are calculated). Typically several hundred files are output.",
1284 "[GRK]phi[grk] and [GRK]psi[grk] dihedrals are calculated in a non-standard way, using H-N-CA-C for [GRK]phi[grk] instead of C(-)-N-CA-C, and N-CA-C-O for [GRK]psi[grk] instead of N-CA-C-N(+). This causes (usually small) discrepancies with the output of other tools like [TT]g_rama[tt].",
1285 "[TT]-r0[tt] option does not work properly",
1286 "Rotamers with multiplicity 2 are printed in [TT]chi.log[tt] as if they had multiplicity 3, with the 3rd (g(+)) always having probability 0"
1290 static int r0 = 1, ndeg = 1, maxchi = 2;
1291 static gmx_bool bAll = FALSE;
1292 static gmx_bool bPhi = FALSE, bPsi = FALSE, bOmega = FALSE;
1293 static real bfac_init = -1.0, bfac_max = 0;
1294 static const char *maxchistr[] = { NULL, "0", "1", "2", "3", "4", "5", "6", NULL };
1295 static gmx_bool bRama = FALSE, bShift = FALSE, bViol = FALSE, bRamOmega = FALSE;
1296 static gmx_bool bNormHisto = TRUE, bChiProduct = FALSE, bHChi = FALSE, bRAD = FALSE, bPBC = TRUE;
1297 static real core_frac = 0.5;
1299 { "-r0", FALSE, etINT, {&r0},
1300 "starting residue" },
1301 { "-phi", FALSE, etBOOL, {&bPhi},
1302 "Output for [GRK]phi[grk] dihedral angles" },
1303 { "-psi", FALSE, etBOOL, {&bPsi},
1304 "Output for [GRK]psi[grk] dihedral angles" },
1305 { "-omega", FALSE, etBOOL, {&bOmega},
1306 "Output for [GRK]omega[grk] dihedrals (peptide bonds)" },
1307 { "-rama", FALSE, etBOOL, {&bRama},
1308 "Generate [GRK]phi[grk]/[GRK]psi[grk] and [GRK]chi[grk][SUB]1[sub]/[GRK]chi[grk][SUB]2[sub] Ramachandran plots" },
1309 { "-viol", FALSE, etBOOL, {&bViol},
1310 "Write a file that gives 0 or 1 for violated Ramachandran angles" },
1311 { "-periodic", FALSE, etBOOL, {&bPBC},
1312 "Print dihedral angles modulo 360 degrees" },
1313 { "-all", FALSE, etBOOL, {&bAll},
1314 "Output separate files for every dihedral." },
1315 { "-rad", FALSE, etBOOL, {&bRAD},
1316 "in angle vs time files, use radians rather than degrees."},
1317 { "-shift", FALSE, etBOOL, {&bShift},
1318 "Compute chemical shifts from [GRK]phi[grk]/[GRK]psi[grk] angles" },
1319 { "-binwidth", FALSE, etINT, {&ndeg},
1320 "bin width for histograms (degrees)" },
1321 { "-core_rotamer", FALSE, etREAL, {&core_frac},
1322 "only the central [TT]-core_rotamer[tt]*(360/multiplicity) belongs to each rotamer (the rest is assigned to rotamer 0)" },
1323 { "-maxchi", FALSE, etENUM, {maxchistr},
1324 "calculate first ndih [GRK]chi[grk] dihedrals" },
1325 { "-normhisto", FALSE, etBOOL, {&bNormHisto},
1326 "Normalize histograms" },
1327 { "-ramomega", FALSE, etBOOL, {&bRamOmega},
1328 "compute average omega as a function of [GRK]phi[grk]/[GRK]psi[grk] and plot it in an [TT].xpm[tt] plot" },
1329 { "-bfact", FALSE, etREAL, {&bfac_init},
1330 "B-factor value for [TT].pdb[tt] file for atoms with no calculated dihedral order parameter"},
1331 { "-chi_prod", FALSE, etBOOL, {&bChiProduct},
1332 "compute a single cumulative rotamer for each residue"},
1333 { "-HChi", FALSE, etBOOL, {&bHChi},
1334 "Include dihedrals to sidechain hydrogens"},
1335 { "-bmax", FALSE, etREAL, {&bfac_max},
1336 "Maximum B-factor on any of the atoms that make up a dihedral, for the dihedral angle to be considere in the statistics. Applies to database work where a number of X-Ray structures is analyzed. [TT]-bmax[tt] <= 0 means no limit." }
1340 int natoms, nlist, idum, nbin;
1345 char title[256], grpname[256];
1347 gmx_bool bChi, bCorr, bSSHisto;
1348 gmx_bool bDo_rt, bDo_oh, bDo_ot, bDo_jc;
1349 real dt = 0, traj_t_ns;
1351 gmx_residuetype_t rt;
1353 atom_id isize, *index;
1354 int ndih, nactdih, nf;
1355 real **dih, *trans_frac, *aver_angle, *time;
1356 int i, j, **chi_lookup, *multiplicity;
1359 { efSTX, "-s", NULL, ffREAD },
1360 { efTRX, "-f", NULL, ffREAD },
1361 { efXVG, "-o", "order", ffWRITE },
1362 { efPDB, "-p", "order", ffOPTWR },
1363 { efDAT, "-ss", "ssdump", ffOPTRD },
1364 { efXVG, "-jc", "Jcoupling", ffWRITE },
1365 { efXVG, "-corr", "dihcorr", ffOPTWR },
1366 { efLOG, "-g", "chi", ffWRITE },
1367 /* add two more arguments copying from g_angle */
1368 { efXVG, "-ot", "dihtrans", ffOPTWR },
1369 { efXVG, "-oh", "trhisto", ffOPTWR },
1370 { efXVG, "-rt", "restrans", ffOPTWR },
1371 { efXVG, "-cp", "chiprodhisto", ffOPTWR }
1373 #define NFILE asize(fnm)
1378 ppa = add_acf_pargs(&npargs, pa);
1379 parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
1380 NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs,
1383 /* Handle result from enumerated type */
1384 sscanf(maxchistr[0], "%d", &maxchi);
1385 bChi = (maxchi > 0);
1387 log = ffopen(ftp2fn(efLOG, NFILE, fnm), "w");
1396 /* set some options */
1397 bDo_rt = (opt2bSet("-rt", NFILE, fnm));
1398 bDo_oh = (opt2bSet("-oh", NFILE, fnm));
1399 bDo_ot = (opt2bSet("-ot", NFILE, fnm));
1400 bDo_jc = (opt2bSet("-jc", NFILE, fnm));
1401 bCorr = (opt2bSet("-corr", NFILE, fnm));
1404 fprintf(stderr, "Will calculate autocorrelation\n");
1407 if (core_frac > 1.0)
1409 fprintf(stderr, "core_rotamer fraction > 1.0 ; will use 1.0\n");
1412 if (core_frac < 0.0)
1414 fprintf(stderr, "core_rotamer fraction < 0.0 ; will use 0.0\n");
1418 if (maxchi > MAXCHI)
1421 "Will only calculate first %d Chi dihedrals in stead of %d.\n",
1425 bSSHisto = ftp2bSet(efDAT, NFILE, fnm);
1428 /* Find the chi angles using atoms struct and a list of amino acids */
1429 get_stx_coordnum(ftp2fn(efSTX, NFILE, fnm), &natoms);
1430 init_t_atoms(&atoms, natoms, TRUE);
1432 read_stx_conf(ftp2fn(efSTX, NFILE, fnm), title, &atoms, x, NULL, &ePBC, box);
1433 fprintf(log, "Title: %s\n", title);
1435 gmx_residuetype_init(&rt);
1436 dlist = mk_dlist(log, &atoms, &nlist, bPhi, bPsi, bChi, bHChi, maxchi, r0, rt);
1437 fprintf(stderr, "%d residues with dihedrals found\n", nlist);
1441 gmx_fatal(FARGS, "No dihedrals in your structure!\n");
1444 /* Make a linear index for reading all. */
1445 index = make_chi_ind(nlist, dlist, &ndih);
1447 fprintf(stderr, "%d dihedrals found\n", ndih);
1451 /* COMPUTE ALL DIHEDRALS! */
1452 read_ang_dih(ftp2fn(efTRX, NFILE, fnm), FALSE, TRUE, FALSE, bPBC, 1, &idum,
1453 &nf, &time, isize, index, &trans_frac, &aver_angle, dih, oenv);
1455 dt = (time[nf-1]-time[0])/(nf-1); /* might want this for corr or n. transit*/
1460 gmx_fatal(FARGS, "Need at least 2 frames for correlation");
1464 /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi
1465 * pass nactdih instead of ndih to low_ana_dih_trans and get_chi_product_traj
1466 * to prevent accessing off end of arrays when maxchi < 5 or 6. */
1467 nactdih = reset_em_all(nlist, dlist, nf, dih, maxchi);
1471 dump_em_all(nlist, dlist, nf, time, dih, maxchi, bPhi, bPsi, bChi, bOmega, bRAD, oenv);
1474 /* Histogramming & J coupling constants & calc of S2 order params */
1475 histogramming(log, nbin, rt, nf, maxchi, dih, nlist, dlist, index,
1476 bPhi, bPsi, bOmega, bChi,
1477 bNormHisto, bSSHisto, ftp2fn(efDAT, NFILE, fnm), bfac_max, &atoms,
1478 bDo_jc, opt2fn("-jc", NFILE, fnm), oenv);
1482 * added multiplicity */
1484 snew(multiplicity, ndih);
1485 mk_multiplicity_lookup(multiplicity, maxchi, dih, nlist, dlist, ndih);
1487 strcpy(grpname, "All residues, ");
1490 strcat(grpname, "Phi ");
1494 strcat(grpname, "Psi ");
1498 strcat(grpname, "Omega ");
1502 strcat(grpname, "Chi 1-");
1503 sprintf(grpname + strlen(grpname), "%i", maxchi);
1507 low_ana_dih_trans(bDo_ot, opt2fn("-ot", NFILE, fnm),
1508 bDo_oh, opt2fn("-oh", NFILE, fnm), maxchi,
1509 dih, nlist, dlist, nf, nactdih, grpname, multiplicity,
1510 time, FALSE, core_frac, oenv);
1512 /* Order parameters */
1513 order_params(log, opt2fn("-o", NFILE, fnm), maxchi, nlist, dlist,
1514 ftp2fn_null(efPDB, NFILE, fnm), bfac_init,
1515 &atoms, x, ePBC, box, bPhi, bPsi, bChi, oenv);
1517 /* Print ramachandran maps! */
1520 do_rama(nf, nlist, dlist, dih, bViol, bRamOmega, oenv);
1525 do_pp2shifts(log, nf, nlist, dlist, dih);
1528 /* rprint S^2, transitions, and rotamer occupancies to log */
1529 traj_t_ns = 0.001 * (time[nf-1]-time[0]);
1530 pr_dlist(log, nlist, dlist, traj_t_ns, edPrintST, bPhi, bPsi, bChi, bOmega, maxchi);
1531 pr_dlist(log, nlist, dlist, traj_t_ns, edPrintRO, bPhi, bPsi, bChi, bOmega, maxchi);
1533 /* transitions to xvg */
1536 print_transitions(opt2fn("-rt", NFILE, fnm), maxchi, nlist, dlist,
1537 &atoms, x, box, bPhi, bPsi, bChi, traj_t_ns, oenv);
1540 /* chi_product trajectories (ie one "rotamer number" for each residue) */
1541 if (bChiProduct && bChi)
1543 snew(chi_lookup, nlist);
1544 for (i = 0; i < nlist; i++)
1546 snew(chi_lookup[i], maxchi);
1548 mk_chi_lookup(chi_lookup, maxchi, dih, nlist, dlist);
1550 get_chi_product_traj(dih, nf, nactdih, nlist,
1551 maxchi, dlist, time, chi_lookup, multiplicity,
1552 FALSE, bNormHisto, core_frac, bAll,
1553 opt2fn("-cp", NFILE, fnm), oenv);
1555 for (i = 0; i < nlist; i++)
1557 sfree(chi_lookup[i]);
1561 /* Correlation comes last because it fucks up the angles */
1564 do_dihcorr(opt2fn("-corr", NFILE, fnm), nf, ndih, dih, dt, nlist, dlist, time,
1565 maxchi, bPhi, bPsi, bChi, bOmega, oenv);
1569 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1570 do_view(oenv, opt2fn("-jc", NFILE, fnm), "-nxy");
1573 do_view(oenv, opt2fn("-corr", NFILE, fnm), "-nxy");
1576 gmx_residuetype_destroy(rt);