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,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
46 #include "gromacs/commandline/pargs.h"
47 #include "gromacs/commandline/viewit.h"
48 #include "gromacs/correlationfunctions/autocorr.h"
49 #include "gromacs/fileio/confio.h"
50 #include "gromacs/fileio/matio.h"
51 #include "gromacs/fileio/pdbio.h"
52 #include "gromacs/fileio/xvgr.h"
53 #include "gromacs/gmxana/gmx_ana.h"
54 #include "gromacs/gmxana/gstat.h"
55 #include "gromacs/math/units.h"
56 #include "gromacs/math/utilities.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/topology/residuetypes.h"
59 #include "gromacs/topology/topology.h"
60 #include "gromacs/utility/arraysize.h"
61 #include "gromacs/utility/cstringutil.h"
62 #include "gromacs/utility/fatalerror.h"
63 #include "gromacs/utility/futil.h"
64 #include "gromacs/utility/smalloc.h"
65 #include "gromacs/utility/stringutil.h"
67 static gmx_bool bAllowed(real phi, real psi)
69 static const char* map[] = { "1100000000000000001111111000000000001111111111111111111111111",
70 "1100000000000000001111110000000000011111111111111111111111111",
71 "1100000000000000001111110000000000011111111111111111111111111",
72 "1100000000000000001111100000000000111111111111111111111111111",
73 "1100000000000000001111100000000000111111111111111111111111111",
74 "1100000000000000001111100000000001111111111111111111111111111",
75 "1100000000000000001111100000000001111111111111111111111111111",
76 "1100000000000000001111100000000011111111111111111111111111111",
77 "1110000000000000001111110000000111111111111111111111111111111",
78 "1110000000000000001111110000001111111111111111111111111111111",
79 "1110000000000000001111111000011111111111111111111111111111111",
80 "1110000000000000001111111100111111111111111111111111111111111",
81 "1110000000000000001111111111111111111111111111111111111111111",
82 "1110000000000000001111111111111111111111111111111111111111111",
83 "1110000000000000001111111111111111111111111111111111111111111",
84 "1110000000000000001111111111111111111111111111111111111111111",
85 "1110000000000000001111111111111110011111111111111111111111111",
86 "1110000000000000001111111111111100000111111111111111111111111",
87 "1110000000000000001111111111111000000000001111111111111111111",
88 "1100000000000000001111111111110000000000000011111111111111111",
89 "1100000000000000001111111111100000000000000011111111111111111",
90 "1000000000000000001111111111000000000000000001111111111111110",
91 "0000000000000000001111111110000000000000000000111111111111100",
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 "0000000000000000000000000000000000000000000000000000000000000",
105 "0000000000000000000000000000000000000000000000000000000000000",
106 "0000000000000000000000000000000000000000000000000000000000000",
107 "0000000000000000000000000000000000111111111111000000000000000",
108 "1100000000000000000000000000000001111111111111100000000000111",
109 "1100000000000000000000000000000001111111111111110000000000111",
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",
127 "0000000000000000000000000000000000000000000000000000000000000",
128 "0000000000000000000000000000000000000000000000000000000000000",
129 "0000000000000000000000000000000000000000000000000000000000000" };
130 #define NPP asize(map)
133 #define INDEX(ppp) (((static_cast<int>(360 + (ppp)*RAD2DEG)) % 360) / 6)
138 return map[x][y] == '1';
141 static int* make_chi_ind(int nl, t_dlist dl[], int* ndih)
146 /* There are nl residues with max edMax dihedrals with 4 atoms each */
147 snew(id, nl * edMax * 4);
150 for (i = 0; (i < nl); i++)
152 /* Phi, fake the first one */
153 dl[i].j0[edPhi] = n / 4;
154 if (dl[i].atm.minC >= 0)
156 id[n++] = dl[i].atm.minC;
160 id[n++] = dl[i].atm.H;
162 id[n++] = dl[i].atm.N;
163 id[n++] = dl[i].atm.Cn[1];
164 id[n++] = dl[i].atm.C;
166 for (i = 0; (i < nl); i++)
168 /* Psi, fake the last one */
169 dl[i].j0[edPsi] = n / 4;
170 id[n++] = dl[i].atm.N;
171 id[n++] = dl[i].atm.Cn[1];
172 id[n++] = dl[i].atm.C;
175 id[n++] = dl[i + 1].atm.N;
179 id[n++] = dl[i].atm.O;
182 for (i = 0; (i < nl); i++)
185 if (has_dihedral(edOmega, &(dl[i])))
187 dl[i].j0[edOmega] = n / 4;
188 id[n++] = dl[i].atm.minCalpha;
189 id[n++] = dl[i].atm.minC;
190 id[n++] = dl[i].atm.N;
191 id[n++] = dl[i].atm.Cn[1];
194 for (Xi = 0; (Xi < MAXCHI); Xi++)
197 for (i = 0; (i < nl); i++)
199 if (dl[i].atm.Cn[Xi + 3] != -1)
201 dl[i].j0[edChi1 + Xi] = n / 4;
202 id[n++] = dl[i].atm.Cn[Xi];
203 id[n++] = dl[i].atm.Cn[Xi + 1];
204 id[n++] = dl[i].atm.Cn[Xi + 2];
205 id[n++] = dl[i].atm.Cn[Xi + 3];
214 static void do_dihcorr(const char* fn,
227 const gmx_output_env_t* oenv)
229 char name1[256], name2[256];
232 do_autocorr(fn, oenv, "Dihedral Autocorrelation Function", nf, ndih, dih, dt, eacCos, FALSE);
235 for (i = 0; (i < nlist); i++)
239 print_one(oenv, "corrphi", dlist[i].name, "Phi ACF for", "C(t)", nf / 2, time, dih[j]);
243 for (i = 0; (i < nlist); i++)
247 print_one(oenv, "corrpsi", dlist[i].name, "Psi ACF for", "C(t)", nf / 2, time, dih[j]);
251 for (i = 0; (i < nlist); i++)
253 if (has_dihedral(edOmega, &dlist[i]))
257 print_one(oenv, "corromega", dlist[i].name, "Omega ACF for", "C(t)", nf / 2, time, dih[j]);
262 for (Xi = 0; (Xi < maxchi); Xi++)
264 sprintf(name1, "corrchi%d", Xi + 1);
265 sprintf(name2, "Chi%d ACF for", Xi + 1);
266 for (i = 0; (i < nlist); i++)
268 if (dlist[i].atm.Cn[Xi + 3] != -1)
272 print_one(oenv, name1, dlist[i].name, name2, "C(t)", nf / 2, time, dih[j]);
278 fprintf(stderr, "\n");
281 static void copy_dih_data(const real in[], real out[], int nf, gmx_bool bLEAVE)
283 /* if bLEAVE, do nothing to data in copying to out
284 * otherwise multiply by 180/pi to convert rad to deg */
293 mult = (180.0 / M_PI);
295 for (i = 0; (i < nf); i++)
297 out[i] = in[i] * mult;
301 static void dump_em_all(int nlist,
312 const gmx_output_env_t* oenv)
314 char name[256], titlestr[256], ystr[256];
321 std::strcpy(ystr, "Angle (rad)");
325 std::strcpy(ystr, "Angle (degrees)");
330 for (i = 0; (i < nlist); i++)
332 /* grs debug printf("OK i %d j %d\n", i, j) ; */
335 copy_dih_data(dih[j], data, nf, bRAD);
336 print_one(oenv, "phi", dlist[i].name, "\\xf\\f{}", ystr, nf, time, data);
340 for (i = 0; (i < nlist); i++)
344 copy_dih_data(dih[j], data, nf, bRAD);
345 print_one(oenv, "psi", dlist[i].name, "\\xy\\f{}", ystr, nf, time, data);
349 for (i = 0; (i < nlist); i++)
351 if (has_dihedral(edOmega, &(dlist[i])))
355 copy_dih_data(dih[j], data, nf, bRAD);
356 print_one(oenv, "omega", dlist[i].name, "\\xw\\f{}", ystr, nf, time, data);
362 for (Xi = 0; (Xi < maxchi); Xi++)
364 for (i = 0; (i < nlist); i++)
366 if (dlist[i].atm.Cn[Xi + 3] != -1)
370 sprintf(name, "chi%d", Xi + 1);
371 sprintf(titlestr, "\\xc\\f{}\\s%d\\N", Xi + 1);
372 copy_dih_data(dih[j], data, nf, bRAD);
373 print_one(oenv, name, dlist[i].name, titlestr, ystr, nf, time, data);
379 fprintf(stderr, "\n");
382 static void reset_one(real dih[], int nf, real phase)
386 for (j = 0; (j < nf); j++)
389 while (dih[j] < -M_PI)
393 while (dih[j] >= M_PI)
400 static int reset_em_all(int nlist, t_dlist dlist[], int nf, real** dih, int maxchi)
407 for (i = 0; (i < nlist); i++)
409 if (dlist[i].atm.minC == -1)
411 reset_one(dih[j++], nf, M_PI);
415 reset_one(dih[j++], nf, 0);
419 for (i = 0; (i < nlist - 1); i++)
421 reset_one(dih[j++], nf, 0);
423 /* last Psi is faked from O */
424 reset_one(dih[j++], nf, M_PI);
427 for (i = 0; (i < nlist); i++)
429 if (has_dihedral(edOmega, &dlist[i]))
431 reset_one(dih[j++], nf, 0);
434 /* Chi 1 thru maxchi */
435 for (Xi = 0; (Xi < maxchi); Xi++)
437 for (i = 0; (i < nlist); i++)
439 if (dlist[i].atm.Cn[Xi + 3] != -1)
441 reset_one(dih[j], nf, 0);
446 fprintf(stderr, "j after resetting (nr. active dihedrals) = %d\n", j);
450 static void histogramming(FILE* log,
467 const t_atoms* atoms,
470 const gmx_output_env_t* oenv)
472 /* also gets 3J couplings and order parameters S2 */
473 // Avoid warnings about narrowing conversions from double to real
475 # pragma warning(disable : 4838)
477 t_karplus kkkphi[] = { { "J_NHa1", 6.51, -1.76, 1.6, -M_PI / 3, 0.0, 0.0 },
478 { "J_NHa2", 6.51, -1.76, 1.6, M_PI / 3, 0.0, 0.0 },
479 { "J_HaC'", 4.0, 1.1, 0.1, 0.0, 0.0, 0.0 },
480 { "J_NHCb", 4.7, -1.5, -0.2, M_PI / 3, 0.0, 0.0 },
481 { "J_Ci-1Hai", 4.5, -1.3, -1.2, 2 * M_PI / 3, 0.0, 0.0 } };
482 t_karplus kkkpsi[] = { { "J_HaN", -0.88, -0.61, -0.27, M_PI / 3, 0.0, 0.0 } };
483 t_karplus kkkchi1[] = { { "JHaHb2", 9.5, -1.6, 1.8, -M_PI / 3, 0, 0.0 },
484 { "JHaHb3", 9.5, -1.6, 1.8, 0, 0, 0.0 } };
486 # pragma warning(default : 4838)
488 #define NKKKPHI asize(kkkphi)
489 #define NKKKPSI asize(kkkpsi)
490 #define NKKKCHI asize(kkkchi1)
491 #define NJC (NKKKPHI + NKKKPSI + NKKKCHI)
493 FILE * fp, *ssfp[3] = { nullptr, nullptr, nullptr };
494 const char* sss[3] = { "sheet", "helix", "coil" };
498 int**** his_aa_ss = nullptr;
499 int *** his_aa, *histmp;
500 int i, j, k, m, n, nn, Dih, nres, hindex, angle;
501 gmx_bool bBfac, bOccup;
502 char hisfile[256], hhisfile[256], title[256], *ss_str = nullptr;
504 const char* residue_name;
507 rt_size = rt->numberOfEntries();
510 fp = gmx_ffopen(ssdump, "r");
511 if (1 != fscanf(fp, "%d", &nres))
513 gmx_fatal(FARGS, "Error reading from file %s", ssdump);
516 snew(ss_str, nres + 1);
517 if (1 != fscanf(fp, "%s", ss_str))
519 gmx_fatal(FARGS, "Error reading from file %s", ssdump);
523 /* Four dimensional array... Very cool */
525 for (i = 0; (i < 3); i++)
527 snew(his_aa_ss[i], rt_size + 1);
528 for (j = 0; (j <= rt_size); j++)
530 snew(his_aa_ss[i][j], edMax);
531 for (Dih = 0; (Dih < edMax); Dih++)
533 snew(his_aa_ss[i][j][Dih], nbin + 1);
539 for (Dih = 0; (Dih < edMax); Dih++)
541 snew(his_aa[Dih], rt_size + 1);
542 for (i = 0; (i <= rt_size); i++)
544 snew(his_aa[Dih][i], nbin + 1);
551 for (i = 0; (i < nlist); i++)
559 for (Dih = 0; (Dih < NONCHI + maxchi); Dih++)
561 for (i = 0; (i < nlist); i++)
563 if (((Dih < edOmega)) || ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i]))))
564 || ((Dih > edOmega) && (dlist[i].atm.Cn[Dih - NONCHI + 3] != -1)))
566 make_histo(log, nf, dih[j], nbin, histmp, -M_PI, M_PI);
570 /* Assume there is only one structure, the first.
571 * Compute index in histogram.
573 /* Check the atoms to see whether their B-factors are low enough
574 * Check atoms to see their occupancy is 1.
576 bBfac = bOccup = TRUE;
577 for (nn = 0; (nn < 4); nn++, n++)
579 bBfac = bBfac && (atoms->pdbinfo[index[n]].bfac <= bfac_max);
580 bOccup = bOccup && (atoms->pdbinfo[index[n]].occup == 1);
582 if (bOccup && ((bfac_max <= 0) || bBfac))
584 hindex = static_cast<int>(((dih[j][0] + M_PI) * nbin) / (2 * M_PI));
585 range_check(hindex, 0, nbin);
587 /* Assign dihedral to either of the structure determined
590 switch (ss_str[dlist[i].resnr])
592 case 'E': his_aa_ss[0][dlist[i].index][Dih][hindex]++; break;
593 case 'H': his_aa_ss[1][dlist[i].index][Dih][hindex]++; break;
594 default: his_aa_ss[2][dlist[i].index][Dih][hindex]++; break;
599 fprintf(debug, "Res. %d has imcomplete occupancy or bfacs > %g\n", dlist[i].resnr, bfac_max);
610 calc_distribution_props(nbin, histmp, -M_PI, NKKKPHI, kkkphi, &S2);
612 for (m = 0; (m < NKKKPHI); m++)
614 Jc[i][m] = kkkphi[m].Jc;
615 Jcsig[i][m] = kkkphi[m].Jcsig;
619 calc_distribution_props(nbin, histmp, -M_PI, NKKKPSI, kkkpsi, &S2);
621 for (m = 0; (m < NKKKPSI); m++)
623 Jc[i][NKKKPHI + m] = kkkpsi[m].Jc;
624 Jcsig[i][NKKKPHI + m] = kkkpsi[m].Jcsig;
628 calc_distribution_props(nbin, histmp, -M_PI, NKKKCHI, kkkchi1, &S2);
629 for (m = 0; (m < NKKKCHI); m++)
631 Jc[i][NKKKPHI + NKKKPSI + m] = kkkchi1[m].Jc;
632 Jcsig[i][NKKKPHI + NKKKPSI + m] = kkkchi1[m].Jcsig;
635 default: /* covers edOmega and higher Chis than Chi1 */
636 calc_distribution_props(nbin, histmp, -M_PI, 0, nullptr, &S2);
639 dlist[i].S2[Dih] = S2;
641 /* Sum distribution per amino acid type as well */
642 for (k = 0; (k < nbin); k++)
644 his_aa[Dih][dlist[i].index][k] += histmp[k];
649 else /* dihed not defined */
651 dlist[i].S2[Dih] = 0.0;
657 /* Print out Jcouplings */
658 fprintf(log, "\n *** J-Couplings from simulation (plus std. dev.) ***\n\n");
659 fprintf(log, "Residue ");
660 for (i = 0; (i < NKKKPHI); i++)
662 fprintf(log, "%7s SD", kkkphi[i].name);
664 for (i = 0; (i < NKKKPSI); i++)
666 fprintf(log, "%7s SD", kkkpsi[i].name);
668 for (i = 0; (i < NKKKCHI); i++)
670 fprintf(log, "%7s SD", kkkchi1[i].name);
673 for (i = 0; (i < NJC + 1); i++)
675 fprintf(log, "------------");
678 for (i = 0; (i < nlist); i++)
680 fprintf(log, "%-10s", dlist[i].name);
681 for (j = 0; (j < NJC); j++)
683 fprintf(log, " %5.2f %4.2f", Jc[i][j], Jcsig[i][j]);
689 /* and to -jc file... */
692 fp = xvgropen(fn, "\\S3\\NJ-Couplings from Karplus Equation", "Residue", "Coupling", oenv);
694 for (i = 0; (i < NKKKPHI); i++)
696 leg[i] = gmx_strdup(kkkphi[i].name);
698 for (i = 0; (i < NKKKPSI); i++)
700 leg[i + NKKKPHI] = gmx_strdup(kkkpsi[i].name);
702 for (i = 0; (i < NKKKCHI); i++)
704 leg[i + NKKKPHI + NKKKPSI] = gmx_strdup(kkkchi1[i].name);
706 xvgr_legend(fp, NJC, leg, oenv);
707 fprintf(fp, "%5s ", "#Res.");
708 for (i = 0; (i < NJC); i++)
710 fprintf(fp, "%10s ", leg[i]);
713 for (i = 0; (i < nlist); i++)
715 fprintf(fp, "%5d ", dlist[i].resnr);
716 for (j = 0; (j < NJC); j++)
718 fprintf(fp, " %8.3f", Jc[i][j]);
723 for (i = 0; (i < NJC); i++)
728 /* finished -jc stuff */
730 snew(normhisto, nbin);
731 for (i = 0; (i < rt_size); i++)
733 for (Dih = 0; (Dih < edMax); Dih++)
735 /* First check whether something is in there */
736 for (j = 0; (j < nbin); j++)
738 if (his_aa[Dih][i][j] != 0)
744 && ((bPhi && (Dih == edPhi)) || (bPsi && (Dih == edPsi))
745 || (bOmega && (Dih == edOmega)) || (bChi && (Dih >= edChi1))))
749 normalize_histo(nbin, his_aa[Dih][i], (360.0 / nbin), normhisto);
752 residue_name = rt->nameFromResidueIndex(i).c_str();
756 sprintf(hisfile, "histo-phi%s", residue_name);
757 sprintf(title, "\\xf\\f{} Distribution for %s", residue_name);
760 sprintf(hisfile, "histo-psi%s", residue_name);
761 sprintf(title, "\\xy\\f{} Distribution for %s", residue_name);
764 sprintf(hisfile, "histo-omega%s", residue_name);
765 sprintf(title, "\\xw\\f{} Distribution for %s", residue_name);
768 sprintf(hisfile, "histo-chi%d%s", Dih - NONCHI + 1, residue_name);
769 sprintf(title, "\\xc\\f{}\\s%d\\N Distribution for %s", Dih - NONCHI + 1, residue_name);
771 std::strcpy(hhisfile, hisfile);
772 std::strcat(hhisfile, ".xvg");
773 fp = xvgropen(hhisfile, title, "Degrees", "", oenv);
774 if (output_env_get_print_xvgr_codes(oenv))
776 fprintf(fp, "@ with g0\n");
778 xvgr_world(fp, -180, 0, 180, 0.1, oenv);
779 if (output_env_get_print_xvgr_codes(oenv))
782 "# this effort to set graph size fails unless you run with -autoscale "
783 "none or -autoscale y flags\n");
784 fprintf(fp, "@ xaxis tick on\n");
785 fprintf(fp, "@ xaxis tick major 90\n");
786 fprintf(fp, "@ xaxis tick minor 30\n");
787 fprintf(fp, "@ xaxis ticklabel prec 0\n");
788 fprintf(fp, "@ yaxis tick off\n");
789 fprintf(fp, "@ yaxis ticklabel off\n");
790 fprintf(fp, "@ type xy\n");
794 for (k = 0; (k < 3); k++)
796 std::string sshisfile = gmx::formatString("%s-%s.xvg", hisfile, sss[k]);
797 ssfp[k] = gmx_ffopen(sshisfile, "w");
800 for (j = 0; (j < nbin); j++)
802 angle = -180 + (360 / nbin) * j;
805 fprintf(fp, "%5d %10g\n", angle, normhisto[j]);
809 fprintf(fp, "%5d %10d\n", angle, his_aa[Dih][i][j]);
813 for (k = 0; (k < 3); k++)
815 fprintf(ssfp[k], "%5d %10d\n", angle, his_aa_ss[k][i][Dih][j]);
819 fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
823 for (k = 0; (k < 3); k++)
825 fprintf(ssfp[k], "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
826 gmx_ffclose(ssfp[k]);
836 /* Four dimensional array... Very cool */
837 for (i = 0; (i < 3); i++)
839 for (j = 0; (j <= rt_size); j++)
841 for (Dih = 0; (Dih < edMax); Dih++)
843 sfree(his_aa_ss[i][j][Dih]);
845 sfree(his_aa_ss[i][j]);
854 static FILE* rama_file(const char* fn, const char* title, const char* xaxis, const char* yaxis, const gmx_output_env_t* oenv)
858 fp = xvgropen(fn, title, xaxis, yaxis, oenv);
859 if (output_env_get_print_xvgr_codes(oenv))
861 fprintf(fp, "@ with g0\n");
863 xvgr_world(fp, -180, -180, 180, 180, oenv);
864 if (output_env_get_print_xvgr_codes(oenv))
866 fprintf(fp, "@ xaxis tick on\n");
867 fprintf(fp, "@ xaxis tick major 90\n");
868 fprintf(fp, "@ xaxis tick minor 30\n");
869 fprintf(fp, "@ xaxis ticklabel prec 0\n");
870 fprintf(fp, "@ yaxis tick on\n");
871 fprintf(fp, "@ yaxis tick major 90\n");
872 fprintf(fp, "@ yaxis tick minor 30\n");
873 fprintf(fp, "@ yaxis ticklabel prec 0\n");
874 fprintf(fp, "@ s0 type xy\n");
875 fprintf(fp, "@ s0 symbol 2\n");
876 fprintf(fp, "@ s0 symbol size 0.410000\n");
877 fprintf(fp, "@ s0 symbol fill 1\n");
878 fprintf(fp, "@ s0 symbol color 1\n");
879 fprintf(fp, "@ s0 symbol linewidth 1\n");
880 fprintf(fp, "@ s0 symbol linestyle 1\n");
881 fprintf(fp, "@ s0 symbol center false\n");
882 fprintf(fp, "@ s0 symbol char 0\n");
883 fprintf(fp, "@ s0 skip 0\n");
884 fprintf(fp, "@ s0 linestyle 0\n");
885 fprintf(fp, "@ s0 linewidth 1\n");
886 fprintf(fp, "@ type xy\n");
891 static void do_rama(int nf,
897 const gmx_output_env_t* oenv)
899 FILE * fp, *gp = nullptr;
902 int i, j, k, Xi1, Xi2, Phi, Psi, Om = 0, nlevels;
903 constexpr int NMAT = 120;
904 real ** mat = nullptr, phi, psi, omega, axis[NMAT], lo, hi;
905 t_rgb rlo = { 1.0, 0.0, 0.0 };
906 t_rgb rmid = { 1.0, 1.0, 1.0 };
907 t_rgb rhi = { 0.0, 0.0, 1.0 };
909 for (i = 0; (i < nlist); i++)
911 if ((has_dihedral(edPhi, &(dlist[i]))) && (has_dihedral(edPsi, &(dlist[i]))))
913 sprintf(fn, "ramaPhiPsi%s.xvg", dlist[i].name);
914 fp = rama_file(fn, "Ramachandran Plot", "\\8f\\4 (deg)", "\\8y\\4 (deg)", oenv);
915 bOm = bRamOmega && has_dihedral(edOmega, &(dlist[i]));
918 Om = dlist[i].j0[edOmega];
920 for (j = 0; (j < NMAT); j++)
923 axis[j] = -180 + gmx::exactDiv(360 * j, NMAT);
928 sprintf(fn, "violPhiPsi%s.xvg", dlist[i].name);
929 gp = gmx_ffopen(fn, "w");
931 Phi = dlist[i].j0[edPhi];
932 Psi = dlist[i].j0[edPsi];
933 for (j = 0; (j < nf); j++)
935 phi = RAD2DEG * dih[Phi][j];
936 psi = RAD2DEG * dih[Psi][j];
937 fprintf(fp, "%10g %10g\n", phi, psi);
940 fprintf(gp, "%d\n", static_cast<int>(!bAllowed(dih[Phi][j], RAD2DEG * dih[Psi][j])));
944 omega = RAD2DEG * dih[Om][j];
945 mat[static_cast<int>(((phi * NMAT) / 360) + gmx::exactDiv(NMAT, 2))]
946 [static_cast<int>(((psi * NMAT) / 360) + gmx::exactDiv(NMAT, 2))] += omega;
956 sprintf(fn, "ramomega%s.xpm", dlist[i].name);
957 fp = gmx_ffopen(fn, "w");
959 for (j = 0; (j < NMAT); j++)
961 for (k = 0; (k < NMAT); k++)
964 lo = std::min(mat[j][k], lo);
965 hi = std::max(mat[j][k], hi);
969 if (std::abs(lo) > std::abs(hi))
978 for (j = 0; (j < NMAT); j++)
980 for (k = 0; (k < NMAT); k++)
990 "Omega/Ramachandran Plot",
1007 for (j = 0; (j < NMAT); j++)
1014 if ((has_dihedral(edChi1, &(dlist[i]))) && (has_dihedral(edChi2, &(dlist[i]))))
1016 sprintf(fn, "ramaX1X2%s.xvg", dlist[i].name);
1018 "\\8c\\4\\s1\\N-\\8c\\4\\s2\\N Ramachandran Plot",
1019 "\\8c\\4\\s1\\N (deg)",
1020 "\\8c\\4\\s2\\N (deg)",
1022 Xi1 = dlist[i].j0[edChi1];
1023 Xi2 = dlist[i].j0[edChi2];
1024 for (j = 0; (j < nf); j++)
1026 fprintf(fp, "%10g %10g\n", RAD2DEG * dih[Xi1][j], RAD2DEG * dih[Xi2][j]);
1032 fprintf(stderr, "No chi1 & chi2 angle for %s\n", dlist[i].name);
1038 static void print_transitions(const char* fn, int maxchi, int nlist, t_dlist dlist[], real dt, const gmx_output_env_t* oenv)
1040 /* based on order_params below */
1044 /* must correspond with enum in pp2shift.h:38 */
1046 #define NLEG asize(leg)
1048 leg[0] = gmx_strdup("Phi");
1049 leg[1] = gmx_strdup("Psi");
1050 leg[2] = gmx_strdup("Omega");
1051 leg[3] = gmx_strdup("Chi1");
1052 leg[4] = gmx_strdup("Chi2");
1053 leg[5] = gmx_strdup("Chi3");
1054 leg[6] = gmx_strdup("Chi4");
1055 leg[7] = gmx_strdup("Chi5");
1056 leg[8] = gmx_strdup("Chi6");
1058 /* Print order parameters */
1059 fp = xvgropen(fn, "Dihedral Rotamer Transitions", "Residue", "Transitions/ns", oenv);
1060 xvgr_legend(fp, NONCHI + maxchi, leg, oenv);
1062 fprintf(fp, "%5s ", "#Res.");
1063 fprintf(fp, "%10s %10s %10s ", leg[edPhi], leg[edPsi], leg[edOmega]);
1064 for (Xi = 0; Xi < maxchi; Xi++)
1066 fprintf(fp, "%10s ", leg[NONCHI + Xi]);
1070 for (i = 0; (i < nlist); i++)
1072 fprintf(fp, "%5d ", dlist[i].resnr);
1073 for (Dih = 0; (Dih < NONCHI + maxchi); Dih++)
1075 fprintf(fp, "%10.3f ", dlist[i].ntr[Dih] / dt);
1077 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
1083 static void order_params(FILE* log,
1097 const gmx_output_env_t* oenv)
1104 /* except for S2Min/Max, must correspond with enum in pp2shift.h:38 */
1105 const char* const_leg[2 + edMax] = { "S2Min", "S2Max", "Phi", "Psi", "Omega", "Chi1",
1106 "Chi2", "Chi3", "Chi4", "Chi5", "Chi6" };
1107 #define NLEG asize(leg)
1109 char* leg[2 + edMax];
1111 for (i = 0; i < NLEG; i++)
1113 leg[i] = gmx_strdup(const_leg[i]);
1116 /* Print order parameters */
1117 fp = xvgropen(fn, "Dihedral Order Parameters", "Residue", "S2", oenv);
1118 xvgr_legend(fp, 2 + NONCHI + maxchi, const_leg, oenv);
1120 for (Dih = 0; (Dih < edMax); Dih++)
1125 fprintf(fp, "%5s ", "#Res.");
1126 fprintf(fp, "%10s %10s ", leg[0], leg[1]);
1127 fprintf(fp, "%10s %10s %10s ", leg[2 + edPhi], leg[2 + edPsi], leg[2 + edOmega]);
1128 for (Xi = 0; Xi < maxchi; Xi++)
1130 fprintf(fp, "%10s ", leg[2 + NONCHI + Xi]);
1134 for (i = 0; (i < nlist); i++)
1138 for (Dih = 0; (Dih < NONCHI + maxchi); Dih++)
1140 if (dlist[i].S2[Dih] != 0)
1142 if (dlist[i].S2[Dih] > S2Max)
1144 S2Max = dlist[i].S2[Dih];
1146 if (dlist[i].S2[Dih] < S2Min)
1148 S2Min = dlist[i].S2[Dih];
1151 if (dlist[i].S2[Dih] > 0.8)
1156 fprintf(fp, "%5d ", dlist[i].resnr);
1157 fprintf(fp, "%10.3f %10.3f ", S2Min, S2Max);
1158 for (Dih = 0; (Dih < NONCHI + maxchi); Dih++)
1160 fprintf(fp, "%10.3f ", dlist[i].S2[Dih]);
1163 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
1167 if (nullptr != pdbfn)
1171 atoms->havePdbInfo = TRUE;
1173 if (nullptr == atoms->pdbinfo)
1175 snew(atoms->pdbinfo, atoms->nr);
1177 for (i = 0; (i < atoms->nr); i++)
1179 atoms->pdbinfo[i].bfac = bfac_init;
1182 for (i = 0; (i < nlist); i++)
1184 atoms->pdbinfo[dlist[i].atm.N].bfac = -dlist[i].S2[0]; /* Phi */
1185 atoms->pdbinfo[dlist[i].atm.H].bfac = -dlist[i].S2[0]; /* Phi */
1186 atoms->pdbinfo[dlist[i].atm.C].bfac = -dlist[i].S2[1]; /* Psi */
1187 atoms->pdbinfo[dlist[i].atm.O].bfac = -dlist[i].S2[1]; /* Psi */
1188 for (Xi = 0; (Xi < maxchi); Xi++) /* Chi's */
1190 if (dlist[i].atm.Cn[Xi + 3] != -1)
1192 atoms->pdbinfo[dlist[i].atm.Cn[Xi + 1]].bfac = -dlist[i].S2[NONCHI + Xi];
1197 fp = gmx_ffopen(pdbfn, "w");
1198 fprintf(fp, "REMARK generated by g_chi\n");
1201 "B-factor field contains negative of dihedral order parameters\n");
1202 write_pdbfile(fp, nullptr, atoms, x, pbcType, box, ' ', 0, nullptr);
1203 x0 = y0 = z0 = 1000.0;
1204 for (i = 0; (i < atoms->nr); i++)
1206 x0 = std::min(x0, x[i][XX]);
1207 y0 = std::min(y0, x[i][YY]);
1208 z0 = std::min(z0, x[i][ZZ]);
1210 x0 *= 10.0; /* nm -> angstrom */
1211 y0 *= 10.0; /* nm -> angstrom */
1212 z0 *= 10.0; /* nm -> angstrom */
1213 for (i = 0; (i < 10); i++)
1215 gmx_fprintf_pdb_atomline(fp,
1216 PdbRecordType::Atom,
1234 fprintf(log, "Dihedrals with S2 > 0.8\n");
1235 fprintf(log, "Dihedral: ");
1238 fprintf(log, " Phi ");
1242 fprintf(log, " Psi ");
1246 for (Xi = 0; (Xi < maxchi); Xi++)
1248 fprintf(log, " %s ", leg[2 + NONCHI + Xi]);
1251 fprintf(log, "\nNumber: ");
1254 fprintf(log, "%4d ", nh[0]);
1258 fprintf(log, "%4d ", nh[1]);
1262 for (Xi = 0; (Xi < maxchi); Xi++)
1264 fprintf(log, "%4d ", nh[NONCHI + Xi]);
1269 for (i = 0; i < NLEG; i++)
1275 int gmx_chi(int argc, char* argv[])
1277 const char* desc[] = {
1278 "[THISMODULE] computes [GRK]phi[grk], [GRK]psi[grk], [GRK]omega[grk],",
1279 "and [GRK]chi[grk] dihedrals for all your",
1280 "amino acid backbone and sidechains.",
1281 "It can compute dihedral angle as a function of time, and as",
1282 "histogram distributions.",
1283 "The distributions [TT](histo-(dihedral)(RESIDUE).xvg[tt]) are cumulative over all ",
1284 "residues of each type.[PAR]",
1285 "If option [TT]-corr[tt] is given, the program will",
1286 "calculate dihedral autocorrelation functions. The function used",
1287 "is C(t) = [CHEVRON][COS][GRK]chi[grk]([GRK]tau[grk])[cos] ",
1288 "[COS][GRK]chi[grk]([GRK]tau[grk]+t)[cos][chevron]. The use of cosines",
1289 "rather than angles themselves, resolves the problem of periodicity.",
1290 "(Van der Spoel & Berendsen (1997), Biophys. J. 72, 2032-2041).",
1291 "Separate files for each dihedral of each residue",
1292 "[TT](corr(dihedral)(RESIDUE)(nresnr).xvg[tt]) are output, as well as a",
1293 "file containing the information for all residues (argument of [TT]-corr[tt]).[PAR]",
1294 "With option [TT]-all[tt], the angles themselves as a function of time for",
1295 "each residue are printed to separate files [TT](dihedral)(RESIDUE)(nresnr).xvg[tt].",
1296 "These can be in radians or degrees.[PAR]",
1297 "A log file (argument [TT]-g[tt]) is also written. This contains",
1299 " * information about the number of residues of each type.",
1300 " * The NMR ^3J coupling constants from the Karplus equation.",
1301 " * a table for each residue of the number of transitions between ",
1302 " rotamers per nanosecond, and the order parameter S^2 of each dihedral.",
1303 " * a table for each residue of the rotamer occupancy.",
1305 "All rotamers are taken as 3-fold, except for [GRK]omega[grk] and [GRK]chi[grk] dihedrals",
1306 "to planar groups (i.e. [GRK]chi[grk][SUB]2[sub] of aromatics, Asp and Asn; ",
1307 "[GRK]chi[grk][SUB]3[sub] of Glu",
1308 "and Gln; and [GRK]chi[grk][SUB]4[sub] of Arg), which are 2-fold. \"rotamer 0\" means ",
1309 "that the dihedral was not in the core region of each rotamer. ",
1310 "The width of the core region can be set with [TT]-core_rotamer[tt][PAR]",
1312 "The S^2 order parameters are also output to an [REF].xvg[ref] file",
1313 "(argument [TT]-o[tt] ) and optionally as a [REF].pdb[ref] file with",
1314 "the S^2 values as B-factor (argument [TT]-p[tt]). ",
1315 "The total number of rotamer transitions per timestep",
1316 "(argument [TT]-ot[tt]), the number of transitions per rotamer",
1317 "(argument [TT]-rt[tt]), and the ^3J couplings (argument [TT]-jc[tt]), ",
1318 "can also be written to [REF].xvg[ref] files. Note that the analysis",
1319 "of rotamer transitions assumes that the supplied trajectory frames",
1320 "are equally spaced in time.[PAR]",
1322 "If [TT]-chi_prod[tt] is set (and [TT]-maxchi[tt] > 0), cumulative rotamers, e.g.",
1323 "1+9([GRK]chi[grk][SUB]1[sub]-1)+3([GRK]chi[grk][SUB]2[sub]-1)+",
1324 "([GRK]chi[grk][SUB]3[sub]-1) (if the residue has three 3-fold ",
1325 "dihedrals and [TT]-maxchi[tt] >= 3)",
1326 "are calculated. As before, if any dihedral is not in the core region,",
1327 "the rotamer is taken to be 0. The occupancies of these cumulative ",
1328 "rotamers (starting with rotamer 0) are written to the file",
1329 "that is the argument of [TT]-cp[tt], and if the [TT]-all[tt] flag",
1330 "is given, the rotamers as functions of time",
1331 "are written to [TT]chiproduct(RESIDUE)(nresnr).xvg[tt] ",
1332 "and their occupancies to [TT]histo-chiproduct(RESIDUE)(nresnr).xvg[tt].[PAR]",
1334 "The option [TT]-r[tt] generates a contour plot of the average [GRK]omega[grk] angle",
1335 "as a function of the [GRK]phi[grk] and [GRK]psi[grk] angles, that is, in a Ramachandran ",
1336 "plot the average [GRK]omega[grk] angle is plotted using color coding.",
1340 const char* bugs[] = {
1341 "Produces MANY output files (up to about 4 times the number of residues in the "
1342 "protein, twice that if autocorrelation functions are calculated). Typically "
1343 "several hundred files are output.",
1344 "[GRK]phi[grk] and [GRK]psi[grk] dihedrals are calculated in a "
1345 "non-standard way, using H-N-CA-C for [GRK]phi[grk] instead of "
1346 "C(-)-N-CA-C, and N-CA-C-O for [GRK]psi[grk] instead of N-CA-C-N(+). "
1347 "This causes (usually small) discrepancies with the output of other "
1348 "tools like [gmx-rama].",
1349 "[TT]-r0[tt] option does not work properly",
1350 "Rotamers with multiplicity 2 are printed in [TT]chi.log[tt] as if they had ",
1351 "multiplicity 3, with the 3rd (g(+)) always having probability 0"
1355 static int r0 = 1, ndeg = 1, maxchi = 2;
1356 static gmx_bool bAll = FALSE;
1357 static gmx_bool bPhi = FALSE, bPsi = FALSE, bOmega = FALSE;
1358 static real bfac_init = -1.0, bfac_max = 0;
1359 static const char* maxchistr[] = { nullptr, "0", "1", "2", "3", "4", "5", "6", nullptr };
1360 static gmx_bool bRama = FALSE, bShift = FALSE, bViol = FALSE, bRamOmega = FALSE;
1361 static gmx_bool bNormHisto = TRUE, bChiProduct = FALSE, bHChi = FALSE, bRAD = FALSE, bPBC = TRUE;
1362 static real core_frac = 0.5;
1364 { "-r0", FALSE, etINT, { &r0 }, "starting residue" },
1365 { "-phi", FALSE, etBOOL, { &bPhi }, "Output for [GRK]phi[grk] dihedral angles" },
1366 { "-psi", FALSE, etBOOL, { &bPsi }, "Output for [GRK]psi[grk] dihedral angles" },
1371 "Output for [GRK]omega[grk] dihedrals (peptide bonds)" },
1376 "Generate [GRK]phi[grk]/[GRK]psi[grk] and "
1377 "[GRK]chi[grk][SUB]1[sub]/[GRK]chi[grk][SUB]2[sub] Ramachandran plots" },
1382 "Write a file that gives 0 or 1 for violated Ramachandran angles" },
1383 { "-periodic", FALSE, etBOOL, { &bPBC }, "Print dihedral angles modulo 360 degrees" },
1384 { "-all", FALSE, etBOOL, { &bAll }, "Output separate files for every dihedral." },
1389 "in angle vs time files, use radians rather than degrees." },
1394 "Compute chemical shifts from [GRK]phi[grk]/[GRK]psi[grk] angles" },
1395 { "-binwidth", FALSE, etINT, { &ndeg }, "bin width for histograms (degrees)" },
1400 "only the central [TT]-core_rotamer[tt]\\*(360/multiplicity) belongs to each rotamer "
1401 "(the rest is assigned to rotamer 0)" },
1402 { "-maxchi", FALSE, etENUM, { maxchistr }, "calculate first ndih [GRK]chi[grk] dihedrals" },
1403 { "-normhisto", FALSE, etBOOL, { &bNormHisto }, "Normalize histograms" },
1408 "compute average omega as a function of [GRK]phi[grk]/[GRK]psi[grk] and plot it in an "
1409 "[REF].xpm[ref] plot" },
1414 "B-factor value for [REF].pdb[ref] file for atoms with no calculated dihedral order "
1420 "compute a single cumulative rotamer for each residue" },
1421 { "-HChi", FALSE, etBOOL, { &bHChi }, "Include dihedrals to sidechain hydrogens" },
1426 "Maximum B-factor on any of the atoms that make up a dihedral, for the dihedral angle to "
1427 "be considere in the statistics. Applies to database work where a number of X-Ray "
1428 "structures is analyzed. [TT]-bmax[tt] <= 0 means no limit." }
1432 int nlist, idum, nbin;
1438 gmx_bool bChi, bCorr, bSSHisto;
1439 gmx_bool bDo_rt, bDo_oh, bDo_ot, bDo_jc;
1440 real dt = 0, traj_t_ns;
1441 gmx_output_env_t* oenv;
1444 int ndih, nactdih, nf;
1445 real **dih, *trans_frac, *aver_angle, *time;
1446 int i, **chi_lookup, *multiplicity;
1448 t_filenm fnm[] = { { efSTX, "-s", nullptr, ffREAD },
1449 { efTRX, "-f", nullptr, ffREAD },
1450 { efXVG, "-o", "order", ffWRITE },
1451 { efPDB, "-p", "order", ffOPTWR },
1452 { efDAT, "-ss", "ssdump", ffOPTRD },
1453 { efXVG, "-jc", "Jcoupling", ffWRITE },
1454 { efXVG, "-corr", "dihcorr", ffOPTWR },
1455 { efLOG, "-g", "chi", ffWRITE },
1456 /* add two more arguments copying from g_angle */
1457 { efXVG, "-ot", "dihtrans", ffOPTWR },
1458 { efXVG, "-oh", "trhisto", ffOPTWR },
1459 { efXVG, "-rt", "restrans", ffOPTWR },
1460 { efXVG, "-cp", "chiprodhisto", ffOPTWR } };
1461 #define NFILE asize(fnm)
1466 ppa = add_acf_pargs(&npargs, pa);
1467 if (!parse_common_args(
1468 &argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs, &oenv))
1474 /* Handle result from enumerated type */
1475 sscanf(maxchistr[0], "%d", &maxchi);
1476 bChi = (maxchi > 0);
1478 log = gmx_ffopen(ftp2fn(efLOG, NFILE, fnm), "w");
1487 /* set some options */
1488 bDo_rt = (opt2bSet("-rt", NFILE, fnm));
1489 bDo_oh = (opt2bSet("-oh", NFILE, fnm));
1490 bDo_ot = (opt2bSet("-ot", NFILE, fnm));
1491 bDo_jc = (opt2bSet("-jc", NFILE, fnm));
1492 bCorr = (opt2bSet("-corr", NFILE, fnm));
1495 fprintf(stderr, "Will calculate autocorrelation\n");
1498 if (core_frac > 1.0)
1500 fprintf(stderr, "core_rotamer fraction > 1.0 ; will use 1.0\n");
1503 if (core_frac < 0.0)
1505 fprintf(stderr, "core_rotamer fraction < 0.0 ; will use 0.0\n");
1509 if (maxchi > MAXCHI)
1511 fprintf(stderr, "Will only calculate first %d Chi dihedrals instead of %d.\n", MAXCHI, maxchi);
1514 bSSHisto = ftp2bSet(efDAT, NFILE, fnm);
1517 /* Find the chi angles using atoms struct and a list of amino acids */
1520 read_tps_conf(ftp2fn(efSTX, NFILE, fnm), top, &pbcType, &x, nullptr, box, FALSE);
1521 t_atoms& atoms = top->atoms;
1522 if (atoms.pdbinfo == nullptr)
1524 snew(atoms.pdbinfo, atoms.nr);
1526 fprintf(log, "Title: %s\n", *top->name);
1529 dlist = mk_dlist(log, &atoms, &nlist, bPhi, bPsi, bChi, bHChi, maxchi, r0, &rt);
1530 fprintf(stderr, "%d residues with dihedrals found\n", nlist);
1534 gmx_fatal(FARGS, "No dihedrals in your structure!\n");
1537 /* Make a linear index for reading all. */
1538 index = make_chi_ind(nlist, dlist, &ndih);
1540 fprintf(stderr, "%d dihedrals found\n", ndih);
1544 /* COMPUTE ALL DIHEDRALS! */
1546 ftp2fn(efTRX, NFILE, fnm), FALSE, TRUE, FALSE, bPBC, 1, &idum, &nf, &time, isize, index, &trans_frac, &aver_angle, dih, oenv);
1548 dt = (time[nf - 1] - time[0]) / (nf - 1); /* might want this for corr or n. transit*/
1553 gmx_fatal(FARGS, "Need at least 2 frames for correlation");
1557 /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi
1558 * pass nactdih instead of ndih to low_ana_dih_trans and get_chi_product_traj
1559 * to prevent accessing off end of arrays when maxchi < 5 or 6. */
1560 nactdih = reset_em_all(nlist, dlist, nf, dih, maxchi);
1564 dump_em_all(nlist, dlist, nf, time, dih, maxchi, bPhi, bPsi, bChi, bOmega, bRAD, oenv);
1567 /* Histogramming & J coupling constants & calc of S2 order params */
1583 ftp2fn(efDAT, NFILE, fnm),
1587 opt2fn("-jc", NFILE, fnm),
1592 * added multiplicity */
1594 snew(multiplicity, ndih);
1595 mk_multiplicity_lookup(multiplicity, maxchi, nlist, dlist, ndih);
1597 std::strcpy(grpname, "All residues, ");
1600 std::strcat(grpname, "Phi ");
1604 std::strcat(grpname, "Psi ");
1608 std::strcat(grpname, "Omega ");
1612 std::strcat(grpname, "Chi 1-");
1613 sprintf(grpname + std::strlen(grpname), "%i", maxchi);
1617 low_ana_dih_trans(bDo_ot,
1618 opt2fn("-ot", NFILE, fnm),
1620 opt2fn("-oh", NFILE, fnm),
1634 /* Order parameters */
1636 opt2fn("-o", NFILE, fnm),
1640 ftp2fn_null(efPDB, NFILE, fnm),
1651 /* Print ramachandran maps! */
1654 do_rama(nf, nlist, dlist, dih, bViol, bRamOmega, oenv);
1659 do_pp2shifts(log, nf, nlist, dlist, dih);
1662 /* rprint S^2, transitions, and rotamer occupancies to log */
1663 traj_t_ns = 0.001 * (time[nf - 1] - time[0]);
1664 pr_dlist(log, nlist, dlist, traj_t_ns, edPrintST, bPhi, bPsi, bChi, bOmega, maxchi);
1665 pr_dlist(log, nlist, dlist, traj_t_ns, edPrintRO, bPhi, bPsi, bChi, bOmega, maxchi);
1667 /* transitions to xvg */
1670 print_transitions(opt2fn("-rt", NFILE, fnm), maxchi, nlist, dlist, traj_t_ns, oenv);
1673 /* chi_product trajectories (ie one "rotamer number" for each residue) */
1674 if (bChiProduct && bChi)
1676 snew(chi_lookup, nlist);
1677 for (i = 0; i < nlist; i++)
1679 snew(chi_lookup[i], maxchi);
1681 mk_chi_lookup(chi_lookup, maxchi, nlist, dlist);
1683 get_chi_product_traj(dih,
1695 opt2fn("-cp", NFILE, fnm),
1698 for (i = 0; i < nlist; i++)
1700 sfree(chi_lookup[i]);
1704 /* Correlation comes last because it messes up the angles */
1707 do_dihcorr(opt2fn("-corr", NFILE, fnm), nf, ndih, dih, dt, nlist, dlist, time, maxchi, bPhi, bPsi, bChi, bOmega, oenv);
1711 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1712 do_view(oenv, opt2fn("-jc", NFILE, fnm), "-nxy");
1715 do_view(oenv, opt2fn("-corr", NFILE, fnm), "-nxy");