2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
43 #include "gromacs/fileio/confio.h"
44 #include "gromacs/fileio/pdbio.h"
45 #include "gromacs/utility/fatalerror.h"
46 #include "gromacs/utility/futil.h"
49 #include "gromacs/math/utilities.h"
52 #include "gromacs/utility/smalloc.h"
53 #include "gromacs/commandline/pargs.h"
54 #include "gromacs/fileio/tpxio.h"
60 #include "gromacs/fileio/matio.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 = gmx_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] = gmx_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");
810 gmx_ffclose(ssfp[k]);
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 = gmx_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 = gmx_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[], real dt,
999 const output_env_t oenv)
1001 /* based on order_params below */
1006 /* must correspond with enum in pp2shift.h:38 */
1008 #define NLEG asize(leg)
1010 leg[0] = strdup("Phi");
1011 leg[1] = strdup("Psi");
1012 leg[2] = strdup("Omega");
1013 leg[3] = strdup("Chi1");
1014 leg[4] = strdup("Chi2");
1015 leg[5] = strdup("Chi3");
1016 leg[6] = strdup("Chi4");
1017 leg[7] = strdup("Chi5");
1018 leg[8] = strdup("Chi6");
1020 /* Print order parameters */
1021 fp = xvgropen(fn, "Dihedral Rotamer Transitions", "Residue", "Transitions/ns",
1023 xvgr_legend(fp, NONCHI+maxchi, (const char**)leg, oenv);
1025 for (Dih = 0; (Dih < edMax); Dih++)
1030 fprintf(fp, "%5s ", "#Res.");
1031 fprintf(fp, "%10s %10s %10s ", leg[edPhi], leg[edPsi], leg[edOmega]);
1032 for (Xi = 0; Xi < maxchi; Xi++)
1034 fprintf(fp, "%10s ", leg[NONCHI+Xi]);
1038 for (i = 0; (i < nlist); i++)
1040 fprintf(fp, "%5d ", dlist[i].resnr);
1041 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1043 fprintf(fp, "%10.3f ", dlist[i].ntr[Dih]/dt);
1045 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
1051 static void order_params(FILE *log,
1052 const char *fn, int maxchi, int nlist, t_dlist dlist[],
1053 const char *pdbfn, real bfac_init,
1054 t_atoms *atoms, rvec x[], int ePBC, matrix box,
1055 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, const output_env_t oenv)
1063 /* except for S2Min/Max, must correspond with enum in pp2shift.h:38 */
1064 const char *const_leg[2+edMax] = {
1065 "S2Min", "S2Max", "Phi", "Psi", "Omega",
1066 "Chi1", "Chi2", "Chi3", "Chi4", "Chi5",
1069 #define NLEG asize(leg)
1073 for (i = 0; i < NLEG; i++)
1075 leg[i] = strdup(const_leg[i]);
1078 /* Print order parameters */
1079 fp = xvgropen(fn, "Dihedral Order Parameters", "Residue", "S2", oenv);
1080 xvgr_legend(fp, 2+NONCHI+maxchi, const_leg, oenv);
1082 for (Dih = 0; (Dih < edMax); Dih++)
1087 fprintf(fp, "%5s ", "#Res.");
1088 fprintf(fp, "%10s %10s ", leg[0], leg[1]);
1089 fprintf(fp, "%10s %10s %10s ", leg[2+edPhi], leg[2+edPsi], leg[2+edOmega]);
1090 for (Xi = 0; Xi < maxchi; Xi++)
1092 fprintf(fp, "%10s ", leg[2+NONCHI+Xi]);
1096 for (i = 0; (i < nlist); i++)
1100 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1102 if (dlist[i].S2[Dih] != 0)
1104 if (dlist[i].S2[Dih] > S2Max)
1106 S2Max = dlist[i].S2[Dih];
1108 if (dlist[i].S2[Dih] < S2Min)
1110 S2Min = dlist[i].S2[Dih];
1113 if (dlist[i].S2[Dih] > 0.8)
1118 fprintf(fp, "%5d ", dlist[i].resnr);
1119 fprintf(fp, "%10.3f %10.3f ", S2Min, S2Max);
1120 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1122 fprintf(fp, "%10.3f ", dlist[i].S2[Dih]);
1125 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
1133 if (NULL == atoms->pdbinfo)
1135 snew(atoms->pdbinfo, atoms->nr);
1137 for (i = 0; (i < atoms->nr); i++)
1139 atoms->pdbinfo[i].bfac = bfac_init;
1142 for (i = 0; (i < nlist); i++)
1144 atoms->pdbinfo[dlist[i].atm.N].bfac = -dlist[i].S2[0]; /* Phi */
1145 atoms->pdbinfo[dlist[i].atm.H].bfac = -dlist[i].S2[0]; /* Phi */
1146 atoms->pdbinfo[dlist[i].atm.C].bfac = -dlist[i].S2[1]; /* Psi */
1147 atoms->pdbinfo[dlist[i].atm.O].bfac = -dlist[i].S2[1]; /* Psi */
1148 for (Xi = 0; (Xi < maxchi); Xi++) /* Chi's */
1150 if (dlist[i].atm.Cn[Xi+3] != -1)
1152 atoms->pdbinfo[dlist[i].atm.Cn[Xi+1]].bfac = -dlist[i].S2[NONCHI+Xi];
1157 fp = gmx_ffopen(pdbfn, "w");
1158 fprintf(fp, "REMARK generated by g_chi\n");
1159 fprintf(fp, "REMARK "
1160 "B-factor field contains negative of dihedral order parameters\n");
1161 write_pdbfile(fp, NULL, atoms, x, ePBC, box, ' ', 0, NULL, TRUE);
1162 x0 = y0 = z0 = 1000.0;
1163 for (i = 0; (i < atoms->nr); i++)
1165 x0 = min(x0, x[i][XX]);
1166 y0 = min(y0, x[i][YY]);
1167 z0 = min(z0, x[i][ZZ]);
1169 x0 *= 10.0; /* nm -> angstrom */
1170 y0 *= 10.0; /* nm -> angstrom */
1171 z0 *= 10.0; /* nm -> angstrom */
1172 sprintf(buf, "%s%%6.f%%6.2f\n", get_pdbformat());
1173 for (i = 0; (i < 10); i++)
1175 fprintf(fp, buf, "ATOM ", atoms->nr+1+i, "CA", "LEG", ' ',
1176 atoms->nres+1, ' ', x0, y0, z0+(1.2*i), 0.0, -0.1*i);
1181 fprintf(log, "Dihedrals with S2 > 0.8\n");
1182 fprintf(log, "Dihedral: ");
1185 fprintf(log, " Phi ");
1189 fprintf(log, " Psi ");
1193 for (Xi = 0; (Xi < maxchi); Xi++)
1195 fprintf(log, " %s ", leg[2+NONCHI+Xi]);
1198 fprintf(log, "\nNumber: ");
1201 fprintf(log, "%4d ", nh[0]);
1205 fprintf(log, "%4d ", nh[1]);
1209 for (Xi = 0; (Xi < maxchi); Xi++)
1211 fprintf(log, "%4d ", nh[NONCHI+Xi]);
1216 for (i = 0; i < NLEG; i++)
1223 int gmx_chi(int argc, char *argv[])
1225 const char *desc[] = {
1226 "[THISMODULE] computes [GRK]phi[grk], [GRK]psi[grk], [GRK]omega[grk],",
1227 "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 "
1285 "non-standard way, using H-N-CA-C for [GRK]phi[grk] instead of "
1286 "C(-)-N-CA-C, and N-CA-C-O for [GRK]psi[grk] instead of N-CA-C-N(+). "
1287 "This causes (usually small) discrepancies with the output of other "
1288 "tools like [gmx-rama].",
1289 "[TT]-r0[tt] option does not work properly",
1290 "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"
1294 static int r0 = 1, ndeg = 1, maxchi = 2;
1295 static gmx_bool bAll = FALSE;
1296 static gmx_bool bPhi = FALSE, bPsi = FALSE, bOmega = FALSE;
1297 static real bfac_init = -1.0, bfac_max = 0;
1298 static const char *maxchistr[] = { NULL, "0", "1", "2", "3", "4", "5", "6", NULL };
1299 static gmx_bool bRama = FALSE, bShift = FALSE, bViol = FALSE, bRamOmega = FALSE;
1300 static gmx_bool bNormHisto = TRUE, bChiProduct = FALSE, bHChi = FALSE, bRAD = FALSE, bPBC = TRUE;
1301 static real core_frac = 0.5;
1303 { "-r0", FALSE, etINT, {&r0},
1304 "starting residue" },
1305 { "-phi", FALSE, etBOOL, {&bPhi},
1306 "Output for [GRK]phi[grk] dihedral angles" },
1307 { "-psi", FALSE, etBOOL, {&bPsi},
1308 "Output for [GRK]psi[grk] dihedral angles" },
1309 { "-omega", FALSE, etBOOL, {&bOmega},
1310 "Output for [GRK]omega[grk] dihedrals (peptide bonds)" },
1311 { "-rama", FALSE, etBOOL, {&bRama},
1312 "Generate [GRK]phi[grk]/[GRK]psi[grk] and [GRK]chi[grk][SUB]1[sub]/[GRK]chi[grk][SUB]2[sub] Ramachandran plots" },
1313 { "-viol", FALSE, etBOOL, {&bViol},
1314 "Write a file that gives 0 or 1 for violated Ramachandran angles" },
1315 { "-periodic", FALSE, etBOOL, {&bPBC},
1316 "Print dihedral angles modulo 360 degrees" },
1317 { "-all", FALSE, etBOOL, {&bAll},
1318 "Output separate files for every dihedral." },
1319 { "-rad", FALSE, etBOOL, {&bRAD},
1320 "in angle vs time files, use radians rather than degrees."},
1321 { "-shift", FALSE, etBOOL, {&bShift},
1322 "Compute chemical shifts from [GRK]phi[grk]/[GRK]psi[grk] angles" },
1323 { "-binwidth", FALSE, etINT, {&ndeg},
1324 "bin width for histograms (degrees)" },
1325 { "-core_rotamer", FALSE, etREAL, {&core_frac},
1326 "only the central [TT]-core_rotamer[tt]*(360/multiplicity) belongs to each rotamer (the rest is assigned to rotamer 0)" },
1327 { "-maxchi", FALSE, etENUM, {maxchistr},
1328 "calculate first ndih [GRK]chi[grk] dihedrals" },
1329 { "-normhisto", FALSE, etBOOL, {&bNormHisto},
1330 "Normalize histograms" },
1331 { "-ramomega", FALSE, etBOOL, {&bRamOmega},
1332 "compute average omega as a function of [GRK]phi[grk]/[GRK]psi[grk] and plot it in an [TT].xpm[tt] plot" },
1333 { "-bfact", FALSE, etREAL, {&bfac_init},
1334 "B-factor value for [TT].pdb[tt] file for atoms with no calculated dihedral order parameter"},
1335 { "-chi_prod", FALSE, etBOOL, {&bChiProduct},
1336 "compute a single cumulative rotamer for each residue"},
1337 { "-HChi", FALSE, etBOOL, {&bHChi},
1338 "Include dihedrals to sidechain hydrogens"},
1339 { "-bmax", FALSE, etREAL, {&bfac_max},
1340 "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." }
1344 int natoms, nlist, idum, nbin;
1349 char title[256], grpname[256];
1351 gmx_bool bChi, bCorr, bSSHisto;
1352 gmx_bool bDo_rt, bDo_oh, bDo_ot, bDo_jc;
1353 real dt = 0, traj_t_ns;
1355 gmx_residuetype_t rt;
1357 atom_id isize, *index;
1358 int ndih, nactdih, nf;
1359 real **dih, *trans_frac, *aver_angle, *time;
1360 int i, j, **chi_lookup, *multiplicity;
1363 { efSTX, "-s", NULL, ffREAD },
1364 { efTRX, "-f", NULL, ffREAD },
1365 { efXVG, "-o", "order", ffWRITE },
1366 { efPDB, "-p", "order", ffOPTWR },
1367 { efDAT, "-ss", "ssdump", ffOPTRD },
1368 { efXVG, "-jc", "Jcoupling", ffWRITE },
1369 { efXVG, "-corr", "dihcorr", ffOPTWR },
1370 { efLOG, "-g", "chi", ffWRITE },
1371 /* add two more arguments copying from g_angle */
1372 { efXVG, "-ot", "dihtrans", ffOPTWR },
1373 { efXVG, "-oh", "trhisto", ffOPTWR },
1374 { efXVG, "-rt", "restrans", ffOPTWR },
1375 { efXVG, "-cp", "chiprodhisto", ffOPTWR }
1377 #define NFILE asize(fnm)
1382 ppa = add_acf_pargs(&npargs, pa);
1383 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
1384 NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs,
1390 /* Handle result from enumerated type */
1391 sscanf(maxchistr[0], "%d", &maxchi);
1392 bChi = (maxchi > 0);
1394 log = gmx_ffopen(ftp2fn(efLOG, NFILE, fnm), "w");
1403 /* set some options */
1404 bDo_rt = (opt2bSet("-rt", NFILE, fnm));
1405 bDo_oh = (opt2bSet("-oh", NFILE, fnm));
1406 bDo_ot = (opt2bSet("-ot", NFILE, fnm));
1407 bDo_jc = (opt2bSet("-jc", NFILE, fnm));
1408 bCorr = (opt2bSet("-corr", NFILE, fnm));
1411 fprintf(stderr, "Will calculate autocorrelation\n");
1414 if (core_frac > 1.0)
1416 fprintf(stderr, "core_rotamer fraction > 1.0 ; will use 1.0\n");
1419 if (core_frac < 0.0)
1421 fprintf(stderr, "core_rotamer fraction < 0.0 ; will use 0.0\n");
1425 if (maxchi > MAXCHI)
1428 "Will only calculate first %d Chi dihedrals in stead of %d.\n",
1432 bSSHisto = ftp2bSet(efDAT, NFILE, fnm);
1435 /* Find the chi angles using atoms struct and a list of amino acids */
1436 get_stx_coordnum(ftp2fn(efSTX, NFILE, fnm), &natoms);
1437 init_t_atoms(&atoms, natoms, TRUE);
1439 read_stx_conf(ftp2fn(efSTX, NFILE, fnm), title, &atoms, x, NULL, &ePBC, box);
1440 fprintf(log, "Title: %s\n", title);
1442 gmx_residuetype_init(&rt);
1443 dlist = mk_dlist(log, &atoms, &nlist, bPhi, bPsi, bChi, bHChi, maxchi, r0, rt);
1444 fprintf(stderr, "%d residues with dihedrals found\n", nlist);
1448 gmx_fatal(FARGS, "No dihedrals in your structure!\n");
1451 /* Make a linear index for reading all. */
1452 index = make_chi_ind(nlist, dlist, &ndih);
1454 fprintf(stderr, "%d dihedrals found\n", ndih);
1458 /* COMPUTE ALL DIHEDRALS! */
1459 read_ang_dih(ftp2fn(efTRX, NFILE, fnm), FALSE, TRUE, FALSE, bPBC, 1, &idum,
1460 &nf, &time, isize, index, &trans_frac, &aver_angle, dih, oenv);
1462 dt = (time[nf-1]-time[0])/(nf-1); /* might want this for corr or n. transit*/
1467 gmx_fatal(FARGS, "Need at least 2 frames for correlation");
1471 /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi
1472 * pass nactdih instead of ndih to low_ana_dih_trans and get_chi_product_traj
1473 * to prevent accessing off end of arrays when maxchi < 5 or 6. */
1474 nactdih = reset_em_all(nlist, dlist, nf, dih, maxchi);
1478 dump_em_all(nlist, dlist, nf, time, dih, maxchi, bPhi, bPsi, bChi, bOmega, bRAD, oenv);
1481 /* Histogramming & J coupling constants & calc of S2 order params */
1482 histogramming(log, nbin, rt, nf, maxchi, dih, nlist, dlist, index,
1483 bPhi, bPsi, bOmega, bChi,
1484 bNormHisto, bSSHisto, ftp2fn(efDAT, NFILE, fnm), bfac_max, &atoms,
1485 bDo_jc, opt2fn("-jc", NFILE, fnm), oenv);
1489 * added multiplicity */
1491 snew(multiplicity, ndih);
1492 mk_multiplicity_lookup(multiplicity, maxchi, nlist, dlist, ndih);
1494 strcpy(grpname, "All residues, ");
1497 strcat(grpname, "Phi ");
1501 strcat(grpname, "Psi ");
1505 strcat(grpname, "Omega ");
1509 strcat(grpname, "Chi 1-");
1510 sprintf(grpname + strlen(grpname), "%i", maxchi);
1514 low_ana_dih_trans(bDo_ot, opt2fn("-ot", NFILE, fnm),
1515 bDo_oh, opt2fn("-oh", NFILE, fnm), maxchi,
1516 dih, nlist, dlist, nf, nactdih, grpname, multiplicity,
1517 time, FALSE, core_frac, oenv);
1519 /* Order parameters */
1520 order_params(log, opt2fn("-o", NFILE, fnm), maxchi, nlist, dlist,
1521 ftp2fn_null(efPDB, NFILE, fnm), bfac_init,
1522 &atoms, x, ePBC, box, bPhi, bPsi, bChi, oenv);
1524 /* Print ramachandran maps! */
1527 do_rama(nf, nlist, dlist, dih, bViol, bRamOmega, oenv);
1532 do_pp2shifts(log, nf, nlist, dlist, dih);
1535 /* rprint S^2, transitions, and rotamer occupancies to log */
1536 traj_t_ns = 0.001 * (time[nf-1]-time[0]);
1537 pr_dlist(log, nlist, dlist, traj_t_ns, edPrintST, bPhi, bPsi, bChi, bOmega, maxchi);
1538 pr_dlist(log, nlist, dlist, traj_t_ns, edPrintRO, bPhi, bPsi, bChi, bOmega, maxchi);
1540 /* transitions to xvg */
1543 print_transitions(opt2fn("-rt", NFILE, fnm), maxchi, nlist, dlist, traj_t_ns, oenv);
1546 /* chi_product trajectories (ie one "rotamer number" for each residue) */
1547 if (bChiProduct && bChi)
1549 snew(chi_lookup, nlist);
1550 for (i = 0; i < nlist; i++)
1552 snew(chi_lookup[i], maxchi);
1554 mk_chi_lookup(chi_lookup, maxchi, nlist, dlist);
1556 get_chi_product_traj(dih, nf, nactdih,
1557 maxchi, dlist, time, chi_lookup, multiplicity,
1558 FALSE, bNormHisto, core_frac, bAll,
1559 opt2fn("-cp", NFILE, fnm), oenv);
1561 for (i = 0; i < nlist; i++)
1563 sfree(chi_lookup[i]);
1567 /* Correlation comes last because it fucks up the angles */
1570 do_dihcorr(opt2fn("-corr", NFILE, fnm), nf, ndih, dih, dt, nlist, dlist, time,
1571 maxchi, bPhi, bPsi, bChi, bOmega, oenv);
1575 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1576 do_view(oenv, opt2fn("-jc", NFILE, fnm), "-nxy");
1579 do_view(oenv, opt2fn("-corr", NFILE, fnm), "-nxy");
1582 gmx_residuetype_destroy(rt);