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.
48 #include <unordered_set>
51 #include "gromacs/commandline/pargs.h"
52 #include "gromacs/commandline/viewit.h"
53 #include "gromacs/correlationfunctions/autocorr.h"
54 #include "gromacs/fileio/confio.h"
55 #include "gromacs/fileio/matio.h"
56 #include "gromacs/fileio/pdbio.h"
57 #include "gromacs/fileio/xvgr.h"
58 #include "gromacs/gmxana/gmx_ana.h"
59 #include "gromacs/gmxana/gstat.h"
60 #include "gromacs/math/units.h"
61 #include "gromacs/math/utilities.h"
62 #include "gromacs/math/vec.h"
63 #include "gromacs/topology/residuetypes.h"
64 #include "gromacs/topology/topology.h"
65 #include "gromacs/utility/arraysize.h"
66 #include "gromacs/utility/cstringutil.h"
67 #include "gromacs/utility/fatalerror.h"
68 #include "gromacs/utility/futil.h"
69 #include "gromacs/utility/smalloc.h"
70 #include "gromacs/utility/stringutil.h"
71 #include "gromacs/utility/unique_cptr.h"
73 static gmx_bool bAllowed(real phi, real psi)
75 static const char* map[] = { "1100000000000000001111111000000000001111111111111111111111111",
76 "1100000000000000001111110000000000011111111111111111111111111",
77 "1100000000000000001111110000000000011111111111111111111111111",
78 "1100000000000000001111100000000000111111111111111111111111111",
79 "1100000000000000001111100000000000111111111111111111111111111",
80 "1100000000000000001111100000000001111111111111111111111111111",
81 "1100000000000000001111100000000001111111111111111111111111111",
82 "1100000000000000001111100000000011111111111111111111111111111",
83 "1110000000000000001111110000000111111111111111111111111111111",
84 "1110000000000000001111110000001111111111111111111111111111111",
85 "1110000000000000001111111000011111111111111111111111111111111",
86 "1110000000000000001111111100111111111111111111111111111111111",
87 "1110000000000000001111111111111111111111111111111111111111111",
88 "1110000000000000001111111111111111111111111111111111111111111",
89 "1110000000000000001111111111111111111111111111111111111111111",
90 "1110000000000000001111111111111111111111111111111111111111111",
91 "1110000000000000001111111111111110011111111111111111111111111",
92 "1110000000000000001111111111111100000111111111111111111111111",
93 "1110000000000000001111111111111000000000001111111111111111111",
94 "1100000000000000001111111111110000000000000011111111111111111",
95 "1100000000000000001111111111100000000000000011111111111111111",
96 "1000000000000000001111111111000000000000000001111111111111110",
97 "0000000000000000001111111110000000000000000000111111111111100",
98 "0000000000000000000000000000000000000000000000000000000000000",
99 "0000000000000000000000000000000000000000000000000000000000000",
100 "0000000000000000000000000000000000000000000000000000000000000",
101 "0000000000000000000000000000000000000000000000000000000000000",
102 "0000000000000000000000000000000000000000000000000000000000000",
103 "0000000000000000000000000000000000000000000000000000000000000",
104 "0000000000000000000000000000000000000000000000000000000000000",
105 "0000000000000000000000000000000000000000000000000000000000000",
106 "0000000000000000000000000000000000000000000000000000000000000",
107 "0000000000000000000000000000000000000000000000000000000000000",
108 "0000000000000000000000000000000000000000000000000000000000000",
109 "0000000000000000000000000000000000000000000000000000000000000",
110 "0000000000000000000000000000000000000000000000000000000000000",
111 "0000000000000000000000000000000000000000000000000000000000000",
112 "0000000000000000000000000000000000000000000000000000000000000",
113 "0000000000000000000000000000000000111111111111000000000000000",
114 "1100000000000000000000000000000001111111111111100000000000111",
115 "1100000000000000000000000000000001111111111111110000000000111",
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 "0000000000000000000000000000000000000000000000000000000000000",
131 "0000000000000000000000000000000000000000000000000000000000000",
132 "0000000000000000000000000000000000000000000000000000000000000",
133 "0000000000000000000000000000000000000000000000000000000000000",
134 "0000000000000000000000000000000000000000000000000000000000000",
135 "0000000000000000000000000000000000000000000000000000000000000" };
138 #define INDEX(ppp) (((static_cast<int>(360 + (ppp)*gmx::c_rad2Deg)) % 360) / 6)
143 return map[x][y] == '1';
146 static std::vector<int> make_chi_ind(gmx::ArrayRef<t_dlist> dlist)
148 /* There are dlist.size() residues with max edMax dihedrals with 4
149 * atoms each. This may be an over-allocation, which is reduced
151 std::vector<int> id(dlist.size() * edMax * 4);
154 for (auto& dihedral : dlist)
156 /* Phi, fake the first one */
157 dihedral.j0[edPhi] = n / 4;
158 if (dihedral.atm.minC >= 0)
160 id[n++] = dihedral.atm.minC;
164 id[n++] = dihedral.atm.H;
166 id[n++] = dihedral.atm.N;
167 id[n++] = dihedral.atm.Cn[1];
168 id[n++] = dihedral.atm.C;
172 for (auto& dihedral : dlist)
174 /* Psi, fake the last one */
175 dihedral.j0[edPsi] = n / 4;
176 id[n++] = dihedral.atm.N;
177 id[n++] = dihedral.atm.Cn[1];
178 id[n++] = dihedral.atm.C;
179 if (i < (gmx::ssize(dlist) - 1))
181 id[n++] = dlist[i + 1].atm.N;
185 id[n++] = dihedral.atm.O;
190 for (auto& dihedral : dlist)
193 if (has_dihedral(edOmega, dihedral))
195 dihedral.j0[edOmega] = n / 4;
196 id[n++] = dihedral.atm.minCalpha;
197 id[n++] = dihedral.atm.minC;
198 id[n++] = dihedral.atm.N;
199 id[n++] = dihedral.atm.Cn[1];
202 for (int Xi = 0; (Xi < MAXCHI); Xi++)
205 for (auto& dihedral : dlist)
207 if (dihedral.atm.Cn[Xi + 3] != -1)
209 dihedral.j0[edChi1 + Xi] = n / 4;
210 id[n++] = dihedral.atm.Cn[Xi];
211 id[n++] = dihedral.atm.Cn[Xi + 1];
212 id[n++] = dihedral.atm.Cn[Xi + 2];
213 id[n++] = dihedral.atm.Cn[Xi + 3];
222 static void do_dihcorr(const char* fn,
227 gmx::ArrayRef<const t_dlist> dlist,
234 const gmx_output_env_t* oenv)
236 char name1[256], name2[256];
239 do_autocorr(fn, oenv, "Dihedral Autocorrelation Function", nf, ndih, dih, dt, eacCos, FALSE);
242 for (const auto& dihedral : dlist)
246 print_one(oenv, "corrphi", dihedral.name, "Phi ACF for", "C(t)", nf / 2, time, dih[j]);
250 for (const auto& dihedral : dlist)
254 print_one(oenv, "corrpsi", dihedral.name, "Psi ACF for", "C(t)", nf / 2, time, dih[j]);
258 for (const auto& dihedral : dlist)
260 if (has_dihedral(edOmega, dihedral))
264 print_one(oenv, "corromega", dihedral.name, "Omega ACF for", "C(t)", nf / 2, time, dih[j]);
269 for (Xi = 0; (Xi < maxchi); Xi++)
271 sprintf(name1, "corrchi%d", Xi + 1);
272 sprintf(name2, "Chi%d ACF for", Xi + 1);
273 for (const auto& dihedral : dlist)
275 if (dihedral.atm.Cn[Xi + 3] != -1)
279 print_one(oenv, name1, dihedral.name, name2, "C(t)", nf / 2, time, dih[j]);
285 fprintf(stderr, "\n");
288 static void copy_dih_data(const real in[], real out[], int nf, gmx_bool bLEAVE)
290 /* if bLEAVE, do nothing to data in copying to out
291 * otherwise multiply by 180/pi to convert rad to deg */
300 mult = (180.0 / M_PI);
302 for (i = 0; (i < nf); i++)
304 out[i] = in[i] * mult;
308 static void dump_em_all(gmx::ArrayRef<const t_dlist> dlist,
318 const gmx_output_env_t* oenv)
320 char name[256], titlestr[256], ystr[256];
324 gmx::sfree_guard dataGuard(data);
327 std::strcpy(ystr, "Angle (rad)");
331 std::strcpy(ystr, "Angle (degrees)");
336 for (const auto& dihedral : dlist)
340 copy_dih_data(dih[j], data, nf, bRAD);
341 print_one(oenv, "phi", dihedral.name, "\\xf\\f{}", ystr, nf, time, data);
345 for (const auto& dihedral : dlist)
349 copy_dih_data(dih[j], data, nf, bRAD);
350 print_one(oenv, "psi", dihedral.name, "\\xy\\f{}", ystr, nf, time, data);
354 for (const auto& dihedral : dlist)
356 if (has_dihedral(edOmega, dihedral))
360 copy_dih_data(dih[j], data, nf, bRAD);
361 print_one(oenv, "omega", dihedral.name, "\\xw\\f{}", ystr, nf, time, data);
367 for (int Xi = 0; (Xi < maxchi); Xi++)
369 for (const auto& dihedral : dlist)
371 if (dihedral.atm.Cn[Xi + 3] != -1)
375 sprintf(name, "chi%d", Xi + 1);
376 sprintf(titlestr, "\\xc\\f{}\\s%d\\N", Xi + 1);
377 copy_dih_data(dih[j], data, nf, bRAD);
378 print_one(oenv, name, dihedral.name, titlestr, ystr, nf, time, data);
384 fprintf(stderr, "\n");
387 static void reset_one(real dih[], int nf, real phase)
391 for (j = 0; (j < nf); j++)
394 while (dih[j] < -M_PI)
398 while (dih[j] >= M_PI)
405 static int reset_em_all(gmx::ArrayRef<const t_dlist> dlist, int nf, real** dih, int maxchi)
410 for (const auto& dihedral : dlist)
412 if (dihedral.atm.minC == -1)
414 reset_one(dih[j++], nf, M_PI);
418 reset_one(dih[j++], nf, 0);
422 for (size_t i = 0; i < dlist.size() - 1; ++i)
424 reset_one(dih[j++], nf, 0);
426 /* last Psi is faked from O */
427 reset_one(dih[j++], nf, M_PI);
430 for (const auto& dihedral : dlist)
432 if (has_dihedral(edOmega, dihedral))
434 reset_one(dih[j++], nf, 0);
437 /* Chi 1 thru maxchi */
438 for (int Xi = 0; (Xi < maxchi); Xi++)
440 for (const auto& dihedral : dlist)
442 if (dihedral.atm.Cn[Xi + 3] != -1)
444 reset_one(dih[j], nf, 0);
449 fprintf(stderr, "j after resetting (nr. active dihedrals) = %d\n", j);
453 static void histogramming(FILE* log,
458 gmx::ArrayRef<t_dlist> dlist,
459 gmx::ArrayRef<const int> index,
468 const t_atoms* atoms,
471 const gmx_output_env_t* oenv)
473 /* also gets 3J couplings and order parameters S2 */
474 // Avoid warnings about narrowing conversions from double to real
476 # pragma warning(disable : 4838)
478 t_karplus kkkphi[] = { { "J_NHa1", 6.51, -1.76, 1.6, -M_PI / 3, 0.0, 0.0 },
479 { "J_NHa2", 6.51, -1.76, 1.6, M_PI / 3, 0.0, 0.0 },
480 { "J_HaC'", 4.0, 1.1, 0.1, 0.0, 0.0, 0.0 },
481 { "J_NHCb", 4.7, -1.5, -0.2, M_PI / 3, 0.0, 0.0 },
482 { "J_Ci-1Hai", 4.5, -1.3, -1.2, 2 * M_PI / 3, 0.0, 0.0 } };
483 t_karplus kkkpsi[] = { { "J_HaN", -0.88, -0.61, -0.27, M_PI / 3, 0.0, 0.0 } };
484 t_karplus kkkchi1[] = { { "JHaHb2", 9.5, -1.6, 1.8, -M_PI / 3, 0, 0.0 },
485 { "JHaHb3", 9.5, -1.6, 1.8, 0, 0, 0.0 } };
487 # pragma warning(default : 4838)
489 #define NKKKPHI asize(kkkphi)
490 #define NKKKPSI asize(kkkpsi)
491 #define NKKKCHI asize(kkkchi1)
492 #define NJC (NKKKPHI + NKKKPSI + NKKKCHI)
494 FILE * fp, *ssfp[3] = { nullptr, nullptr, nullptr };
495 const char* sss[3] = { "sheet", "helix", "coil" };
499 int m, n, nn, nres, hindex, angle;
500 gmx_bool bBfac, bOccup;
501 char hisfile[256], hhisfile[256], title[256], *ss_str = nullptr;
506 fp = gmx_ffopen(ssdump, "r");
507 if (1 != fscanf(fp, "%d", &nres))
509 gmx_fatal(FARGS, "Error reading from file %s", ssdump);
512 snew(ss_str, nres + 1);
513 if (1 != fscanf(fp, "%s", ss_str))
515 gmx_fatal(FARGS, "Error reading from file %s", ssdump);
521 // Build a list of unique residue names found in the dihedral
522 // list, so we can loop over those unique names conveniently. The
523 // names are the same as the residue names found in rt in the
524 // caller, but ResidueType doesn't yet have a way to loop over its
526 std::unordered_set<std::string> uniqueResidueNames;
527 for (const auto& dihedral : dlist)
529 uniqueResidueNames.emplace(dihedral.residueName);
531 // Build the lookup tables for data relating to the all dihedrals
532 // from each unique residue name represented in the dihedral list.
533 std::array<std::map<std::string, std::vector<std::vector<int>>>, 3> his_aa_ss;
534 std::vector<std::map<std::string, std::vector<int>>> his_aa(edMax);
535 for (const auto& residueName : uniqueResidueNames)
539 for (auto& secondaryStructure : his_aa_ss)
541 secondaryStructure[residueName] =
542 std::vector<std::vector<int>>(edMax, std::vector<int>(nbin, 0));
545 for (auto& dihedraltype : his_aa)
547 dihedraltype[residueName] = std::vector<int>(nbin, 0);
552 snew(Jc, dlist.size());
553 snew(Jcsig, dlist.size());
554 for (size_t i = 0; i < dlist.size(); i++)
562 for (int Dih = 0; (Dih < NONCHI + maxchi); Dih++)
565 for (auto& dihedral : dlist)
567 if (((Dih < edOmega)) || ((Dih == edOmega) && (has_dihedral(edOmega, dihedral)))
568 || ((Dih > edOmega) && (dihedral.atm.Cn[Dih - NONCHI + 3] != -1)))
570 make_histo(log, nf, dih[j], nbin, histmp, -M_PI, M_PI);
574 /* Assume there is only one structure, the first.
575 * Compute index in histogram.
577 /* Check the atoms to see whether their B-factors are low enough
578 * Check atoms to see their occupancy is 1.
580 bBfac = bOccup = TRUE;
581 for (nn = 0; (nn < 4); nn++, n++)
583 bBfac = bBfac && (atoms->pdbinfo[index[n]].bfac <= bfac_max);
584 bOccup = bOccup && (atoms->pdbinfo[index[n]].occup == 1);
586 if (bOccup && ((bfac_max <= 0) || bBfac))
588 hindex = static_cast<int>(((dih[j][0] + M_PI) * nbin) / (2 * M_PI));
589 range_check(hindex, 0, nbin);
591 /* Assign dihedral to either of the structure determined
594 switch (ss_str[dihedral.resnr])
596 case 'E': his_aa_ss[0][dihedral.residueName][Dih][hindex]++; break;
597 case 'H': his_aa_ss[1][dihedral.residueName][Dih][hindex]++; break;
598 default: his_aa_ss[2][dihedral.residueName][Dih][hindex]++; break;
603 fprintf(debug, "Res. %d has imcomplete occupancy or bfacs > %g\n", dihedral.resnr, bfac_max);
614 calc_distribution_props(nbin, histmp, -M_PI, NKKKPHI, kkkphi, &S2);
616 for (m = 0; (m < NKKKPHI); m++)
618 Jc[i][m] = kkkphi[m].Jc;
619 Jcsig[i][m] = kkkphi[m].Jcsig;
623 calc_distribution_props(nbin, histmp, -M_PI, NKKKPSI, kkkpsi, &S2);
625 for (m = 0; (m < NKKKPSI); m++)
627 Jc[i][NKKKPHI + m] = kkkpsi[m].Jc;
628 Jcsig[i][NKKKPHI + m] = kkkpsi[m].Jcsig;
632 calc_distribution_props(nbin, histmp, -M_PI, NKKKCHI, kkkchi1, &S2);
633 for (m = 0; (m < NKKKCHI); m++)
635 Jc[i][NKKKPHI + NKKKPSI + m] = kkkchi1[m].Jc;
636 Jcsig[i][NKKKPHI + NKKKPSI + m] = kkkchi1[m].Jcsig;
639 default: /* covers edOmega and higher Chis than Chi1 */
640 calc_distribution_props(nbin, histmp, -M_PI, 0, nullptr, &S2);
643 dihedral.S2[Dih] = S2;
645 /* Sum distribution per amino acid type as well */
646 for (int k = 0; (k < nbin); k++)
648 his_aa[Dih][dihedral.residueName][k] += histmp[k];
653 else /* dihed not defined */
655 dihedral.S2[Dih] = 0.0;
662 /* Print out Jcouplings */
663 fprintf(log, "\n *** J-Couplings from simulation (plus std. dev.) ***\n\n");
664 fprintf(log, "Residue ");
665 for (int i = 0; (i < NKKKPHI); i++)
667 fprintf(log, "%7s SD", kkkphi[i].name);
669 for (int i = 0; (i < NKKKPSI); i++)
671 fprintf(log, "%7s SD", kkkpsi[i].name);
673 for (int i = 0; (i < NKKKCHI); i++)
675 fprintf(log, "%7s SD", kkkchi1[i].name);
678 for (int i = 0; (i < NJC + 1); i++)
680 fprintf(log, "------------");
685 for (const auto& dihedral : dlist)
687 fprintf(log, "%-10s", dihedral.name);
688 for (int j = 0; (j < NJC); j++)
690 fprintf(log, " %5.2f %4.2f", Jc[i][j], Jcsig[i][j]);
698 /* and to -jc file... */
701 fp = xvgropen(fn, "\\S3\\NJ-Couplings from Karplus Equation", "Residue", "Coupling", oenv);
703 for (int i = 0; (i < NKKKPHI); i++)
705 leg[i] = gmx_strdup(kkkphi[i].name);
707 for (int i = 0; (i < NKKKPSI); i++)
709 leg[i + NKKKPHI] = gmx_strdup(kkkpsi[i].name);
711 for (int i = 0; (i < NKKKCHI); i++)
713 leg[i + NKKKPHI + NKKKPSI] = gmx_strdup(kkkchi1[i].name);
715 xvgr_legend(fp, NJC, leg, oenv);
716 fprintf(fp, "%5s ", "#Res.");
717 for (int i = 0; (i < NJC); i++)
719 fprintf(fp, "%10s ", leg[i]);
724 for (const auto& dihedral : dlist)
726 fprintf(fp, "%5d ", dihedral.resnr);
727 for (int j = 0; (j < NJC); j++)
729 fprintf(fp, " %8.3f", Jc[i][j]);
736 for (int i = 0; (i < NJC); i++)
742 /* finished -jc stuff */
744 std::vector<real> normhisto(nbin);
745 for (const auto& residueName : uniqueResidueNames)
747 for (int Dih = 0; (Dih < edMax); Dih++)
749 /* First check whether something is in there */
751 for (j = 0; (j < nbin); j++)
753 if (his_aa[Dih][residueName][j] != 0)
759 && ((bPhi && (Dih == edPhi)) || (bPsi && (Dih == edPsi))
760 || (bOmega && (Dih == edOmega)) || (bChi && (Dih >= edChi1))))
764 normalize_histo(his_aa[Dih][residueName], (360.0 / nbin), normhisto);
770 sprintf(hisfile, "histo-phi%s", residueName.c_str());
771 sprintf(title, "\\xf\\f{} Distribution for %s", residueName.c_str());
774 sprintf(hisfile, "histo-psi%s", residueName.c_str());
775 sprintf(title, "\\xy\\f{} Distribution for %s", residueName.c_str());
778 sprintf(hisfile, "histo-omega%s", residueName.c_str());
779 sprintf(title, "\\xw\\f{} Distribution for %s", residueName.c_str());
782 sprintf(hisfile, "histo-chi%d%s", Dih - NONCHI + 1, residueName.c_str());
784 "\\xc\\f{}\\s%d\\N Distribution for %s",
786 residueName.c_str());
788 std::strcpy(hhisfile, hisfile);
789 std::strcat(hhisfile, ".xvg");
790 fp = xvgropen(hhisfile, title, "Degrees", "", oenv);
791 if (output_env_get_print_xvgr_codes(oenv))
793 fprintf(fp, "@ with g0\n");
795 xvgr_world(fp, -180, 0, 180, 0.1, oenv);
796 if (output_env_get_print_xvgr_codes(oenv))
799 "# this effort to set graph size fails unless you run with -autoscale "
800 "none or -autoscale y flags\n");
801 fprintf(fp, "@ xaxis tick on\n");
802 fprintf(fp, "@ xaxis tick major 90\n");
803 fprintf(fp, "@ xaxis tick minor 30\n");
804 fprintf(fp, "@ xaxis ticklabel prec 0\n");
805 fprintf(fp, "@ yaxis tick off\n");
806 fprintf(fp, "@ yaxis ticklabel off\n");
807 fprintf(fp, "@ type xy\n");
811 for (int k = 0; (k < 3); k++)
813 std::string sshisfile = gmx::formatString("%s-%s.xvg", hisfile, sss[k]);
814 ssfp[k] = gmx_ffopen(sshisfile, "w");
817 for (int j = 0; (j < nbin); j++)
819 angle = -180 + (360 / nbin) * j;
822 fprintf(fp, "%5d %10g\n", angle, normhisto[j]);
826 fprintf(fp, "%5d %10d\n", angle, his_aa[Dih][residueName][j]);
830 for (int k = 0; (k < 3); k++)
832 fprintf(ssfp[k], "%5d %10d\n", angle, his_aa_ss[k][residueName][Dih][j]);
836 fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
840 for (int k = 0; (k < 3); k++)
842 fprintf(ssfp[k], "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
843 gmx_ffclose(ssfp[k]);
854 for (size_t i = 0; i < dlist.size(); i++)
863 static FILE* rama_file(const char* fn, const char* title, const char* xaxis, const char* yaxis, const gmx_output_env_t* oenv)
867 fp = xvgropen(fn, title, xaxis, yaxis, oenv);
868 if (output_env_get_print_xvgr_codes(oenv))
870 fprintf(fp, "@ with g0\n");
872 xvgr_world(fp, -180, -180, 180, 180, oenv);
873 if (output_env_get_print_xvgr_codes(oenv))
875 fprintf(fp, "@ xaxis tick on\n");
876 fprintf(fp, "@ xaxis tick major 90\n");
877 fprintf(fp, "@ xaxis tick minor 30\n");
878 fprintf(fp, "@ xaxis ticklabel prec 0\n");
879 fprintf(fp, "@ yaxis tick on\n");
880 fprintf(fp, "@ yaxis tick major 90\n");
881 fprintf(fp, "@ yaxis tick minor 30\n");
882 fprintf(fp, "@ yaxis ticklabel prec 0\n");
883 fprintf(fp, "@ s0 type xy\n");
884 fprintf(fp, "@ s0 symbol 2\n");
885 fprintf(fp, "@ s0 symbol size 0.410000\n");
886 fprintf(fp, "@ s0 symbol fill 1\n");
887 fprintf(fp, "@ s0 symbol color 1\n");
888 fprintf(fp, "@ s0 symbol linewidth 1\n");
889 fprintf(fp, "@ s0 symbol linestyle 1\n");
890 fprintf(fp, "@ s0 symbol center false\n");
891 fprintf(fp, "@ s0 symbol char 0\n");
892 fprintf(fp, "@ s0 skip 0\n");
893 fprintf(fp, "@ s0 linestyle 0\n");
894 fprintf(fp, "@ s0 linewidth 1\n");
895 fprintf(fp, "@ type xy\n");
900 static void do_rama(int nf,
901 gmx::ArrayRef<const t_dlist> dlist,
905 const gmx_output_env_t* oenv)
907 FILE * fp, *gp = nullptr;
910 int Xi1, Xi2, Phi, Psi, Om = 0, nlevels;
911 constexpr int NMAT = 120;
912 real ** mat = nullptr, phi, psi, omega, axis[NMAT], lo, hi;
913 t_rgb rlo = { 1.0, 0.0, 0.0 };
914 t_rgb rmid = { 1.0, 1.0, 1.0 };
915 t_rgb rhi = { 0.0, 0.0, 1.0 };
917 for (const auto& dihedral : dlist)
919 if ((has_dihedral(edPhi, dihedral)) && has_dihedral(edPsi, dihedral))
921 sprintf(fn, "ramaPhiPsi%s.xvg", dihedral.name);
922 fp = rama_file(fn, "Ramachandran Plot", "\\8f\\4 (deg)", "\\8y\\4 (deg)", oenv);
923 bOm = bRamOmega && has_dihedral(edOmega, dihedral);
926 Om = dihedral.j0[edOmega];
928 for (int j = 0; (j < NMAT); j++)
931 axis[j] = -180 + gmx::exactDiv(360 * j, NMAT);
936 sprintf(fn, "violPhiPsi%s.xvg", dihedral.name);
937 gp = gmx_ffopen(fn, "w");
939 Phi = dihedral.j0[edPhi];
940 Psi = dihedral.j0[edPsi];
941 for (int j = 0; (j < nf); j++)
943 phi = gmx::c_rad2Deg * dih[Phi][j];
944 psi = gmx::c_rad2Deg * dih[Psi][j];
945 fprintf(fp, "%10g %10g\n", phi, psi);
950 static_cast<int>(!bAllowed(dih[Phi][j], gmx::c_rad2Deg * dih[Psi][j])));
954 omega = gmx::c_rad2Deg * dih[Om][j];
955 mat[static_cast<int>(((phi * NMAT) / 360) + gmx::exactDiv(NMAT, 2))]
956 [static_cast<int>(((psi * NMAT) / 360) + gmx::exactDiv(NMAT, 2))] += omega;
966 sprintf(fn, "ramomega%s.xpm", dihedral.name);
967 fp = gmx_ffopen(fn, "w");
969 for (int j = 0; (j < NMAT); j++)
971 for (int k = 0; (k < NMAT); k++)
974 lo = std::min(mat[j][k], lo);
975 hi = std::max(mat[j][k], hi);
979 if (std::abs(lo) > std::abs(hi))
988 for (int j = 0; (j < NMAT); j++)
990 for (int k = 0; (k < NMAT); k++)
1000 "Omega/Ramachandran Plot",
1017 for (int j = 0; (j < NMAT); j++)
1024 if (has_dihedral(edChi1, dihedral) && has_dihedral(edChi2, dihedral))
1026 sprintf(fn, "ramaX1X2%s.xvg", dihedral.name);
1028 "\\8c\\4\\s1\\N-\\8c\\4\\s2\\N Ramachandran Plot",
1029 "\\8c\\4\\s1\\N (deg)",
1030 "\\8c\\4\\s2\\N (deg)",
1032 Xi1 = dihedral.j0[edChi1];
1033 Xi2 = dihedral.j0[edChi2];
1034 for (int j = 0; (j < nf); j++)
1036 fprintf(fp, "%10g %10g\n", gmx::c_rad2Deg * dih[Xi1][j], gmx::c_rad2Deg * dih[Xi2][j]);
1042 fprintf(stderr, "No chi1 & chi2 angle for %s\n", dihedral.name);
1048 static void print_transitions(const char* fn,
1050 gmx::ArrayRef<const t_dlist> dlist,
1052 const gmx_output_env_t* oenv)
1054 /* based on order_params below */
1057 /* must correspond with enum in gstat.h */
1058 const char* leg[edMax] = {
1059 "Phi", "Psi", "Omega", "Chi1", "Chi2", "Chi3", "Chi4", "Chi5", "Chi6",
1061 /* Print order parameters */
1062 fp = xvgropen(fn, "Dihedral Rotamer Transitions", "Residue", "Transitions/ns", oenv);
1063 xvgr_legend(fp, NONCHI + maxchi, leg, oenv);
1065 fprintf(fp, "%5s ", "#Res.");
1066 fprintf(fp, "%10s %10s %10s ", leg[edPhi], leg[edPsi], leg[edOmega]);
1067 for (int Xi = 0; Xi < maxchi; Xi++)
1069 fprintf(fp, "%10s ", leg[NONCHI + Xi]);
1073 for (const auto& dihedral : dlist)
1075 fprintf(fp, "%5d ", dihedral.resnr);
1076 for (int Dih = 0; (Dih < NONCHI + maxchi); Dih++)
1078 fprintf(fp, "%10.3f ", dihedral.ntr[Dih] / dt);
1080 /* fprintf(fp,"%12s\n",dihedral.name); this confuses xmgrace */
1086 static void order_params(FILE* log,
1089 gmx::ArrayRef<const t_dlist> dlist,
1099 const gmx_output_env_t* oenv)
1105 /* except for S2Min/Max, must correspond with enum in pp2shift.h:38 */
1106 const char* const_leg[2 + edMax] = { "S2Min", "S2Max", "Phi", "Psi", "Omega", "Chi1",
1107 "Chi2", "Chi3", "Chi4", "Chi5", "Chi6" };
1108 #define NLEG asize(leg)
1110 char* leg[2 + edMax];
1112 for (int i = 0; i < NLEG; i++)
1114 leg[i] = gmx_strdup(const_leg[i]);
1117 /* Print order parameters */
1118 fp = xvgropen(fn, "Dihedral Order Parameters", "Residue", "S2", oenv);
1119 xvgr_legend(fp, 2 + NONCHI + maxchi, const_leg, oenv);
1121 for (int Dih = 0; (Dih < edMax); Dih++)
1126 fprintf(fp, "%5s ", "#Res.");
1127 fprintf(fp, "%10s %10s ", leg[0], leg[1]);
1128 fprintf(fp, "%10s %10s %10s ", leg[2 + edPhi], leg[2 + edPsi], leg[2 + edOmega]);
1129 for (int Xi = 0; Xi < maxchi; Xi++)
1131 fprintf(fp, "%10s ", leg[2 + NONCHI + Xi]);
1135 for (const auto& dihedral : dlist)
1139 for (int Dih = 0; (Dih < NONCHI + maxchi); Dih++)
1141 if (dihedral.S2[Dih] != 0)
1143 if (dihedral.S2[Dih] > S2Max)
1145 S2Max = dihedral.S2[Dih];
1147 if (dihedral.S2[Dih] < S2Min)
1149 S2Min = dihedral.S2[Dih];
1152 if (dihedral.S2[Dih] > 0.8)
1157 fprintf(fp, "%5d ", dihedral.resnr);
1158 fprintf(fp, "%10.3f %10.3f ", S2Min, S2Max);
1159 for (int Dih = 0; (Dih < NONCHI + maxchi); Dih++)
1161 fprintf(fp, "%10.3f ", dihedral.S2[Dih]);
1164 /* fprintf(fp,"%12s\n",dihedral.name); this confuses xmgrace */
1168 if (nullptr != pdbfn)
1172 atoms->havePdbInfo = TRUE;
1174 if (nullptr == atoms->pdbinfo)
1176 snew(atoms->pdbinfo, atoms->nr);
1178 for (int i = 0; (i < atoms->nr); i++)
1180 atoms->pdbinfo[i].bfac = bfac_init;
1183 for (const auto& dihedral : dlist)
1185 atoms->pdbinfo[dihedral.atm.N].bfac = -dihedral.S2[0]; /* Phi */
1186 atoms->pdbinfo[dihedral.atm.H].bfac = -dihedral.S2[0]; /* Phi */
1187 atoms->pdbinfo[dihedral.atm.C].bfac = -dihedral.S2[1]; /* Psi */
1188 atoms->pdbinfo[dihedral.atm.O].bfac = -dihedral.S2[1]; /* Psi */
1189 for (int Xi = 0; (Xi < maxchi); Xi++) /* Chi's */
1191 if (dihedral.atm.Cn[Xi + 3] != -1)
1193 atoms->pdbinfo[dihedral.atm.Cn[Xi + 1]].bfac = -dihedral.S2[NONCHI + Xi];
1198 fp = gmx_ffopen(pdbfn, "w");
1199 fprintf(fp, "REMARK generated by g_chi\n");
1202 "B-factor field contains negative of dihedral order parameters\n");
1203 write_pdbfile(fp, nullptr, atoms, x, pbcType, box, ' ', 0, nullptr);
1204 x0 = y0 = z0 = 1000.0;
1205 for (int i = 0; (i < atoms->nr); i++)
1207 x0 = std::min(x0, x[i][XX]);
1208 y0 = std::min(y0, x[i][YY]);
1209 z0 = std::min(z0, x[i][ZZ]);
1211 x0 *= 10.0; /* nm -> angstrom */
1212 y0 *= 10.0; /* nm -> angstrom */
1213 z0 *= 10.0; /* nm -> angstrom */
1214 for (int i = 0; (i < 10); i++)
1216 gmx_fprintf_pdb_atomline(fp,
1217 PdbRecordType::Atom,
1235 fprintf(log, "Dihedrals with S2 > 0.8\n");
1236 fprintf(log, "Dihedral: ");
1239 fprintf(log, " Phi ");
1243 fprintf(log, " Psi ");
1247 for (int Xi = 0; (Xi < maxchi); Xi++)
1249 fprintf(log, " %s ", leg[2 + NONCHI + Xi]);
1252 fprintf(log, "\nNumber: ");
1255 fprintf(log, "%4d ", nh[0]);
1259 fprintf(log, "%4d ", nh[1]);
1263 for (int Xi = 0; (Xi < maxchi); Xi++)
1265 fprintf(log, "%4d ", nh[NONCHI + Xi]);
1270 for (int i = 0; i < NLEG; i++)
1276 int gmx_chi(int argc, char* argv[])
1278 const char* desc[] = {
1279 "[THISMODULE] computes [GRK]phi[grk], [GRK]psi[grk], [GRK]omega[grk],",
1280 "and [GRK]chi[grk] dihedrals for all your",
1281 "amino acid backbone and sidechains.",
1282 "It can compute dihedral angle as a function of time, and as",
1283 "histogram distributions.",
1284 "The distributions [TT](histo-(dihedral)(RESIDUE).xvg[tt]) are cumulative over all ",
1285 "residues of each type.[PAR]",
1286 "If option [TT]-corr[tt] is given, the program will",
1287 "calculate dihedral autocorrelation functions. The function used",
1288 "is C(t) = [CHEVRON][COS][GRK]chi[grk]([GRK]tau[grk])[cos] ",
1289 "[COS][GRK]chi[grk]([GRK]tau[grk]+t)[cos][chevron]. The use of cosines",
1290 "rather than angles themselves, resolves the problem of periodicity.",
1291 "(Van der Spoel & Berendsen (1997), Biophys. J. 72, 2032-2041).",
1292 "Separate files for each dihedral of each residue",
1293 "[TT](corr(dihedral)(RESIDUE)(nresnr).xvg[tt]) are output, as well as a",
1294 "file containing the information for all residues (argument of [TT]-corr[tt]).[PAR]",
1295 "With option [TT]-all[tt], the angles themselves as a function of time for",
1296 "each residue are printed to separate files [TT](dihedral)(RESIDUE)(nresnr).xvg[tt].",
1297 "These can be in radians or degrees.[PAR]",
1298 "A log file (argument [TT]-g[tt]) is also written. This contains",
1300 " * information about the number of residues of each type.",
1301 " * The NMR ^3J coupling constants from the Karplus equation.",
1302 " * a table for each residue of the number of transitions between ",
1303 " rotamers per nanosecond, and the order parameter S^2 of each dihedral.",
1304 " * a table for each residue of the rotamer occupancy.",
1306 "All rotamers are taken as 3-fold, except for [GRK]omega[grk] and [GRK]chi[grk] dihedrals",
1307 "to planar groups (i.e. [GRK]chi[grk][SUB]2[sub] of aromatics, Asp and Asn; ",
1308 "[GRK]chi[grk][SUB]3[sub] of Glu",
1309 "and Gln; and [GRK]chi[grk][SUB]4[sub] of Arg), which are 2-fold. \"rotamer 0\" means ",
1310 "that the dihedral was not in the core region of each rotamer. ",
1311 "The width of the core region can be set with [TT]-core_rotamer[tt][PAR]",
1313 "The S^2 order parameters are also output to an [REF].xvg[ref] file",
1314 "(argument [TT]-o[tt] ) and optionally as a [REF].pdb[ref] file with",
1315 "the S^2 values as B-factor (argument [TT]-p[tt]). ",
1316 "The total number of rotamer transitions per timestep",
1317 "(argument [TT]-ot[tt]), the number of transitions per rotamer",
1318 "(argument [TT]-rt[tt]), and the ^3J couplings (argument [TT]-jc[tt]), ",
1319 "can also be written to [REF].xvg[ref] files. Note that the analysis",
1320 "of rotamer transitions assumes that the supplied trajectory frames",
1321 "are equally spaced in time.[PAR]",
1323 "If [TT]-chi_prod[tt] is set (and [TT]-maxchi[tt] > 0), cumulative rotamers, e.g.",
1324 "1+9([GRK]chi[grk][SUB]1[sub]-1)+3([GRK]chi[grk][SUB]2[sub]-1)+",
1325 "([GRK]chi[grk][SUB]3[sub]-1) (if the residue has three 3-fold ",
1326 "dihedrals and [TT]-maxchi[tt] >= 3)",
1327 "are calculated. As before, if any dihedral is not in the core region,",
1328 "the rotamer is taken to be 0. The occupancies of these cumulative ",
1329 "rotamers (starting with rotamer 0) are written to the file",
1330 "that is the argument of [TT]-cp[tt], and if the [TT]-all[tt] flag",
1331 "is given, the rotamers as functions of time",
1332 "are written to [TT]chiproduct(RESIDUE)(nresnr).xvg[tt] ",
1333 "and their occupancies to [TT]histo-chiproduct(RESIDUE)(nresnr).xvg[tt].[PAR]",
1335 "The option [TT]-r[tt] generates a contour plot of the average [GRK]omega[grk] angle",
1336 "as a function of the [GRK]phi[grk] and [GRK]psi[grk] angles, that is, in a Ramachandran ",
1337 "plot the average [GRK]omega[grk] angle is plotted using color coding.",
1341 const char* bugs[] = {
1342 "Produces MANY output files (up to about 4 times the number of residues in the "
1343 "protein, twice that if autocorrelation functions are calculated). Typically "
1344 "several hundred files are output.",
1345 "[GRK]phi[grk] and [GRK]psi[grk] dihedrals are calculated in a "
1346 "non-standard way, using H-N-CA-C for [GRK]phi[grk] instead of "
1347 "C(-)-N-CA-C, and N-CA-C-O for [GRK]psi[grk] instead of N-CA-C-N(+). "
1348 "This causes (usually small) discrepancies with the output of other "
1349 "tools like [gmx-rama].",
1350 "[TT]-r0[tt] option does not work properly",
1351 "Rotamers with multiplicity 2 are printed in [TT]chi.log[tt] as if they had ",
1352 "multiplicity 3, with the 3rd (g(+)) always having probability 0"
1356 static int r0 = 1, ndeg = 1, maxchi = 2;
1357 static gmx_bool bAll = FALSE;
1358 static gmx_bool bPhi = FALSE, bPsi = FALSE, bOmega = FALSE;
1359 static real bfac_init = -1.0, bfac_max = 0;
1360 static const char* maxchistr[] = { nullptr, "0", "1", "2", "3", "4", "5", "6", nullptr };
1361 static gmx_bool bRama = FALSE, bShift = FALSE, bViol = FALSE, bRamOmega = FALSE;
1362 static gmx_bool bNormHisto = TRUE, bChiProduct = FALSE, bHChi = FALSE, bRAD = FALSE, bPBC = TRUE;
1363 static real core_frac = 0.5;
1365 { "-r0", FALSE, etINT, { &r0 }, "starting residue" },
1366 { "-phi", FALSE, etBOOL, { &bPhi }, "Output for [GRK]phi[grk] dihedral angles" },
1367 { "-psi", FALSE, etBOOL, { &bPsi }, "Output for [GRK]psi[grk] dihedral angles" },
1372 "Output for [GRK]omega[grk] dihedrals (peptide bonds)" },
1377 "Generate [GRK]phi[grk]/[GRK]psi[grk] and "
1378 "[GRK]chi[grk][SUB]1[sub]/[GRK]chi[grk][SUB]2[sub] Ramachandran plots" },
1383 "Write a file that gives 0 or 1 for violated Ramachandran angles" },
1384 { "-periodic", FALSE, etBOOL, { &bPBC }, "Print dihedral angles modulo 360 degrees" },
1385 { "-all", FALSE, etBOOL, { &bAll }, "Output separate files for every dihedral." },
1390 "in angle vs time files, use radians rather than degrees." },
1395 "Compute chemical shifts from [GRK]phi[grk]/[GRK]psi[grk] angles" },
1396 { "-binwidth", FALSE, etINT, { &ndeg }, "bin width for histograms (degrees)" },
1401 "only the central [TT]-core_rotamer[tt]\\*(360/multiplicity) belongs to each rotamer "
1402 "(the rest is assigned to rotamer 0)" },
1403 { "-maxchi", FALSE, etENUM, { maxchistr }, "calculate first ndih [GRK]chi[grk] dihedrals" },
1404 { "-normhisto", FALSE, etBOOL, { &bNormHisto }, "Normalize histograms" },
1409 "compute average omega as a function of [GRK]phi[grk]/[GRK]psi[grk] and plot it in an "
1410 "[REF].xpm[ref] plot" },
1415 "B-factor value for [REF].pdb[ref] file for atoms with no calculated dihedral order "
1421 "compute a single cumulative rotamer for each residue" },
1422 { "-HChi", FALSE, etBOOL, { &bHChi }, "Include dihedrals to sidechain hydrogens" },
1427 "Maximum B-factor on any of the atoms that make up a dihedral, for the dihedral angle to "
1428 "be considere in the statistics. Applies to database work where a number of X-Ray "
1429 "structures is analyzed. [TT]-bmax[tt] <= 0 means no limit." }
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 real **dih, *trans_frac, *aver_angle, *time;
1445 int ** chi_lookup, *multiplicity;
1447 t_filenm fnm[] = { { efSTX, "-s", nullptr, ffREAD },
1448 { efTRX, "-f", nullptr, ffREAD },
1449 { efXVG, "-o", "order", ffWRITE },
1450 { efPDB, "-p", "order", ffOPTWR },
1451 { efDAT, "-ss", "ssdump", ffOPTRD },
1452 { efXVG, "-jc", "Jcoupling", ffWRITE },
1453 { efXVG, "-corr", "dihcorr", ffOPTWR },
1454 { efLOG, "-g", "chi", ffWRITE },
1455 /* add two more arguments copying from g_angle */
1456 { efXVG, "-ot", "dihtrans", ffOPTWR },
1457 { efXVG, "-oh", "trhisto", ffOPTWR },
1458 { efXVG, "-rt", "restrans", ffOPTWR },
1459 { efXVG, "-cp", "chiprodhisto", ffOPTWR } };
1460 #define NFILE asize(fnm)
1465 ppa = add_acf_pargs(&npargs, pa);
1466 gmx::sfree_guard ppaGuard(ppa);
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))
1472 gmx::unique_cptr<gmx_output_env_t, output_env_done> oenvGuard(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 */
1521 open_symtab(&symtab);
1522 gmx::unique_cptr<t_symtab, done_symtab> symtabGuard(&symtab);
1523 readConfAndAtoms(ftp2fn(efSTX, NFILE, fnm), &symtab, &name, &atoms, &pbcType, &x, nullptr, box);
1524 gmx::sfree_guard nameGuard(name);
1525 gmx::sfree_guard xGuard(x);
1526 gmx::unique_cptr<t_atoms, done_atom> atomsGuard(&atoms);
1527 if (atoms.pdbinfo == nullptr)
1529 snew(atoms.pdbinfo, atoms.nr);
1531 fprintf(log, "Title: %s\n", name);
1534 std::vector<t_dlist> dlist = mk_dlist(log, &atoms, bPhi, bPsi, bChi, bHChi, maxchi, r0, &rt);
1535 fprintf(stderr, "%zu residues with dihedrals found\n", dlist.size());
1539 gmx_fatal(FARGS, "No dihedrals in your structure!\n");
1542 /* Make a linear index for reading all dihedral atoms (4 per dihedral). */
1543 std::vector<int> index = make_chi_ind(dlist);
1544 int ndih = index.size() / 4; // 4 atoms per dihedral
1545 fprintf(stderr, "%d dihedrals found\n", ndih);
1549 /* COMPUTE ALL DIHEDRALS! */
1550 read_ang_dih(ftp2fn(efTRX, NFILE, fnm),
1565 gmx::sfree_guard timeGuard(time);
1566 gmx::sfree_guard transFracGuard(trans_frac);
1567 gmx::sfree_guard averAngleGuard(aver_angle);
1569 dt = (time[nf - 1] - time[0]) / (nf - 1); /* might want this for corr or n. transit*/
1574 gmx_fatal(FARGS, "Need at least 2 frames for correlation");
1578 /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi
1579 * pass nactdih instead of ndih to low_ana_dih_trans
1580 * to prevent accessing off end of arrays when maxchi < 5 or 6. */
1581 nactdih = reset_em_all(dlist, nf, dih, maxchi);
1585 dump_em_all(dlist, nf, time, dih, maxchi, bPhi, bPsi, bChi, bOmega, bRAD, oenv);
1588 /* Histogramming & J coupling constants & calc of S2 order params */
1602 ftp2fn(efDAT, NFILE, fnm),
1606 opt2fn("-jc", NFILE, fnm),
1611 * added multiplicity */
1613 snew(multiplicity, ndih);
1614 gmx::sfree_guard multiplicityGuard(multiplicity);
1615 mk_multiplicity_lookup(multiplicity, maxchi, dlist, ndih);
1617 std::strcpy(grpname, "All residues, ");
1620 std::strcat(grpname, "Phi ");
1624 std::strcat(grpname, "Psi ");
1628 std::strcat(grpname, "Omega ");
1632 std::strcat(grpname, "Chi 1-");
1633 sprintf(grpname + std::strlen(grpname), "%i", maxchi);
1637 low_ana_dih_trans(bDo_ot,
1638 opt2fn("-ot", NFILE, fnm),
1640 opt2fn("-oh", NFILE, fnm),
1653 /* Order parameters */
1655 opt2fn("-o", NFILE, fnm),
1658 ftp2fn_null(efPDB, NFILE, fnm),
1669 /* Print ramachandran maps! */
1672 do_rama(nf, dlist, dih, bViol, bRamOmega, oenv);
1677 do_pp2shifts(log, nf, dlist, dih);
1680 /* rprint S^2, transitions, and rotamer occupancies to log */
1681 traj_t_ns = 0.001 * (time[nf - 1] - time[0]);
1682 pr_dlist(log, dlist, traj_t_ns, edPrintST, bPhi, bPsi, bChi, bOmega, maxchi);
1683 pr_dlist(log, dlist, traj_t_ns, edPrintRO, bPhi, bPsi, bChi, bOmega, maxchi);
1685 /* transitions to xvg */
1688 print_transitions(opt2fn("-rt", NFILE, fnm), maxchi, dlist, traj_t_ns, oenv);
1691 /* chi_product trajectories (ie one "rotamer number" for each residue) */
1692 if (bChiProduct && bChi)
1694 snew(chi_lookup, dlist.size());
1695 for (size_t i = 0; i < dlist.size(); i++)
1697 snew(chi_lookup[i], maxchi);
1699 mk_chi_lookup(chi_lookup, maxchi, dlist);
1701 get_chi_product_traj(dih,
1712 opt2fn("-cp", NFILE, fnm),
1715 for (size_t i = 0; i < dlist.size(); i++)
1717 sfree(chi_lookup[i]);
1722 /* Correlation comes last because it messes up the angles */
1725 do_dihcorr(opt2fn("-corr", NFILE, fnm), nf, ndih, dih, dt, dlist, time, maxchi, bPhi, bPsi, bChi, bOmega, oenv);
1729 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1730 do_view(oenv, opt2fn("-jc", NFILE, fnm), "-nxy");
1733 do_view(oenv, opt2fn("-corr", NFILE, fnm), "-nxy");
1736 for (int i = 0; (i < ndih); i++)