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.
45 #include "gromacs/fileio/confio.h"
46 #include "gromacs/fileio/pdbio.h"
47 #include "gromacs/utility/fatalerror.h"
48 #include "gromacs/utility/futil.h"
50 #include "gromacs/legacyheaders/macros.h"
51 #include "gromacs/math/utilities.h"
52 #include "gromacs/math/units.h"
53 #include "gromacs/topology/residuetypes.h"
54 #include "gromacs/utility/smalloc.h"
55 #include "gromacs/commandline/pargs.h"
56 #include "gromacs/fileio/tpxio.h"
57 #include "gromacs/legacyheaders/txtdump.h"
58 #include "gromacs/legacyheaders/typedefs.h"
59 #include "gromacs/math/vec.h"
60 #include "gromacs/fileio/xvgr.h"
61 #include "gromacs/legacyheaders/viewit.h"
62 #include "gromacs/fileio/matio.h"
65 static gmx_bool bAllowed(real phi, real psi)
67 static const char *map[] = {
68 "1100000000000000001111111000000000001111111111111111111111111",
69 "1100000000000000001111110000000000011111111111111111111111111",
70 "1100000000000000001111110000000000011111111111111111111111111",
71 "1100000000000000001111100000000000111111111111111111111111111",
72 "1100000000000000001111100000000000111111111111111111111111111",
73 "1100000000000000001111100000000001111111111111111111111111111",
74 "1100000000000000001111100000000001111111111111111111111111111",
75 "1100000000000000001111100000000011111111111111111111111111111",
76 "1110000000000000001111110000000111111111111111111111111111111",
77 "1110000000000000001111110000001111111111111111111111111111111",
78 "1110000000000000001111111000011111111111111111111111111111111",
79 "1110000000000000001111111100111111111111111111111111111111111",
80 "1110000000000000001111111111111111111111111111111111111111111",
81 "1110000000000000001111111111111111111111111111111111111111111",
82 "1110000000000000001111111111111111111111111111111111111111111",
83 "1110000000000000001111111111111111111111111111111111111111111",
84 "1110000000000000001111111111111110011111111111111111111111111",
85 "1110000000000000001111111111111100000111111111111111111111111",
86 "1110000000000000001111111111111000000000001111111111111111111",
87 "1100000000000000001111111111110000000000000011111111111111111",
88 "1100000000000000001111111111100000000000000011111111111111111",
89 "1000000000000000001111111111000000000000000001111111111111110",
90 "0000000000000000001111111110000000000000000000111111111111100",
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 "0000000000000000000000000000000000000000000000000000000000000",
105 "0000000000000000000000000000000000000000000000000000000000000",
106 "0000000000000000000000000000000000111111111111000000000000000",
107 "1100000000000000000000000000000001111111111111100000000000111",
108 "1100000000000000000000000000000001111111111111110000000000111",
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",
127 "0000000000000000000000000000000000000000000000000000000000000",
128 "0000000000000000000000000000000000000000000000000000000000000"
130 #define NPP asize(map)
133 #define INDEX(ppp) ((((int) (360+ppp*RAD2DEG)) % 360)/6)
137 return (gmx_bool) map[x][y];
140 atom_id *make_chi_ind(int nl, t_dlist dl[], int *ndih)
145 /* There are nl residues with max edMax dihedrals with 4 atoms each */
146 snew(id, nl*edMax*4);
149 for (i = 0; (i < nl); i++)
151 /* Phi, fake the first one */
152 dl[i].j0[edPhi] = n/4;
153 if (dl[i].atm.minC >= 0)
155 id[n++] = dl[i].atm.minC;
159 id[n++] = dl[i].atm.H;
161 id[n++] = dl[i].atm.N;
162 id[n++] = dl[i].atm.Cn[1];
163 id[n++] = dl[i].atm.C;
165 for (i = 0; (i < nl); i++)
167 /* Psi, fake the last one */
168 dl[i].j0[edPsi] = n/4;
169 id[n++] = dl[i].atm.N;
170 id[n++] = dl[i].atm.Cn[1];
171 id[n++] = dl[i].atm.C;
174 id[n++] = dl[i+1].atm.N;
178 id[n++] = dl[i].atm.O;
181 for (i = 1; (i < nl); i++)
184 if (has_dihedral(edOmega, &(dl[i])))
186 dl[i].j0[edOmega] = n/4;
187 id[n++] = dl[i].atm.minCalpha;
188 id[n++] = dl[i].atm.minC;
189 id[n++] = dl[i].atm.N;
190 id[n++] = dl[i].atm.Cn[1];
193 for (Xi = 0; (Xi < MAXCHI); Xi++)
196 for (i = 0; (i < nl); i++)
198 if (dl[i].atm.Cn[Xi+3] != -1)
200 dl[i].j0[edChi1+Xi] = n/4;
201 id[n++] = dl[i].atm.Cn[Xi];
202 id[n++] = dl[i].atm.Cn[Xi+1];
203 id[n++] = dl[i].atm.Cn[Xi+2];
204 id[n++] = dl[i].atm.Cn[Xi+3];
213 int bin(real chi, int mult)
217 return (int) (chi*mult/360.0);
221 static void do_dihcorr(const char *fn, int nf, int ndih, real **dih, real dt,
222 int nlist, t_dlist dlist[], real time[], int maxchi,
223 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bOmega,
224 const output_env_t oenv)
226 char name1[256], name2[256];
229 do_autocorr(fn, oenv, "Dihedral Autocorrelation Function",
230 nf, ndih, dih, dt, eacCos, FALSE);
233 for (i = 0; (i < nlist); i++)
237 print_one(oenv, "corrphi", dlist[i].name, "Phi ACF for", "C(t)", nf/2, time,
242 for (i = 0; (i < nlist); i++)
246 print_one(oenv, "corrpsi", dlist[i].name, "Psi ACF for", "C(t)", nf/2, time,
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)",
263 for (Xi = 0; (Xi < maxchi); Xi++)
265 sprintf(name1, "corrchi%d", Xi+1);
266 sprintf(name2, "Chi%d ACF for", Xi+1);
267 for (i = 0; (i < nlist); i++)
269 if (dlist[i].atm.Cn[Xi+3] != -1)
273 print_one(oenv, name1, dlist[i].name, name2, "C(t)", nf/2, time, dih[j]);
279 fprintf(stderr, "\n");
282 static void copy_dih_data(real in[], real out[], int nf, gmx_bool bLEAVE)
284 /* if bLEAVE, do nothing to data in copying to out
285 * otherwise multiply by 180/pi to convert rad to deg */
296 for (i = 0; (i < nf); i++)
302 static void dump_em_all(int nlist, t_dlist dlist[], int nf, real time[],
303 real **dih, int maxchi,
304 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bOmega, gmx_bool bRAD,
305 const output_env_t oenv)
307 char name[256], titlestr[256], ystr[256];
314 strcpy(ystr, "Angle (rad)");
318 strcpy(ystr, "Angle (degrees)");
323 for (i = 0; (i < nlist); i++)
325 /* grs debug printf("OK i %d j %d\n", i, j) ; */
328 copy_dih_data(dih[j], data, nf, bRAD);
329 print_one(oenv, "phi", dlist[i].name, "\\xf\\f{}", ystr, nf, time, data);
333 for (i = 0; (i < nlist); i++)
337 copy_dih_data(dih[j], data, nf, bRAD);
338 print_one(oenv, "psi", dlist[i].name, "\\xy\\f{}", ystr, nf, time, data);
342 for (i = 0; (i < nlist); i++)
344 if (has_dihedral(edOmega, &(dlist[i])))
348 copy_dih_data(dih[j], data, nf, bRAD);
349 print_one(oenv, "omega", dlist[i].name, "\\xw\\f{}", ystr, nf, time, data);
355 for (Xi = 0; (Xi < maxchi); Xi++)
357 for (i = 0; (i < nlist); i++)
359 if (dlist[i].atm.Cn[Xi+3] != -1)
363 sprintf(name, "chi%d", Xi+1);
364 sprintf(titlestr, "\\xc\\f{}\\s%d\\N", Xi+1);
365 copy_dih_data(dih[j], data, nf, bRAD);
366 print_one(oenv, name, dlist[i].name, titlestr, ystr, nf, time, data);
372 fprintf(stderr, "\n");
375 static void reset_one(real dih[], int nf, real phase)
379 for (j = 0; (j < nf); j++)
382 while (dih[j] < -M_PI)
386 while (dih[j] >= M_PI)
393 static int reset_em_all(int nlist, t_dlist dlist[], int nf,
394 real **dih, int maxchi)
401 for (i = 0; (i < nlist); i++)
403 if (dlist[i].atm.minC == -1)
405 reset_one(dih[j++], nf, M_PI);
409 reset_one(dih[j++], nf, 0);
413 for (i = 0; (i < nlist-1); i++)
415 reset_one(dih[j++], nf, 0);
417 /* last Psi is faked from O */
418 reset_one(dih[j++], nf, M_PI);
421 for (i = 0; (i < nlist); i++)
423 if (has_dihedral(edOmega, &dlist[i]))
425 reset_one(dih[j++], nf, 0);
428 /* Chi 1 thru maxchi */
429 for (Xi = 0; (Xi < maxchi); Xi++)
431 for (i = 0; (i < nlist); i++)
433 if (dlist[i].atm.Cn[Xi+3] != -1)
435 reset_one(dih[j], nf, 0);
440 fprintf(stderr, "j after resetting (nr. active dihedrals) = %d\n", j);
444 static void histogramming(FILE *log, int nbin, gmx_residuetype_t *rt,
445 int nf, int maxchi, real **dih,
446 int nlist, t_dlist dlist[],
448 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bOmega, gmx_bool bChi,
449 gmx_bool bNormalize, gmx_bool bSSHisto, const char *ssdump,
450 real bfac_max, t_atoms *atoms,
451 gmx_bool bDo_jc, const char *fn,
452 const output_env_t oenv)
454 /* also gets 3J couplings and order parameters S2 */
455 t_karplus kkkphi[] = {
456 { "J_NHa1", 6.51, -1.76, 1.6, -M_PI/3, 0.0, 0.0 },
457 { "J_NHa2", 6.51, -1.76, 1.6, M_PI/3, 0.0, 0.0 },
458 { "J_HaC'", 4.0, 1.1, 0.1, 0.0, 0.0, 0.0 },
459 { "J_NHCb", 4.7, -1.5, -0.2, M_PI/3, 0.0, 0.0 },
460 { "J_Ci-1Hai", 4.5, -1.3, -1.2, 2*M_PI/3, 0.0, 0.0 }
462 t_karplus kkkpsi[] = {
463 { "J_HaN", -0.88, -0.61, -0.27, M_PI/3, 0.0, 0.0 }
465 t_karplus kkkchi1[] = {
466 { "JHaHb2", 9.5, -1.6, 1.8, -M_PI/3, 0, 0.0 },
467 { "JHaHb3", 9.5, -1.6, 1.8, 0, 0, 0.0 }
469 #define NKKKPHI asize(kkkphi)
470 #define NKKKPSI asize(kkkpsi)
471 #define NKKKCHI asize(kkkchi1)
472 #define NJC (NKKKPHI+NKKKPSI+NKKKCHI)
474 FILE *fp, *ssfp[3] = {NULL, NULL, NULL};
475 const char *sss[3] = { "sheet", "helix", "coil" };
479 int ****his_aa_ss = NULL;
480 int ***his_aa, **his_aa1, *histmp;
481 int i, j, k, m, n, nn, Dih, nres, hindex, angle;
482 gmx_bool bBfac, bOccup;
483 char hisfile[256], hhisfile[256], sshisfile[256], title[256], *ss_str = NULL;
485 const char *residue_name;
488 rt_size = gmx_residuetype_get_size(rt);
491 fp = gmx_ffopen(ssdump, "r");
492 if (1 != fscanf(fp, "%d", &nres))
494 gmx_fatal(FARGS, "Error reading from file %s", ssdump);
497 snew(ss_str, nres+1);
498 if (1 != fscanf(fp, "%s", ss_str))
500 gmx_fatal(FARGS, "Error reading from file %s", ssdump);
504 /* Four dimensional array... Very cool */
506 for (i = 0; (i < 3); i++)
508 snew(his_aa_ss[i], rt_size+1);
509 for (j = 0; (j <= rt_size); j++)
511 snew(his_aa_ss[i][j], edMax);
512 for (Dih = 0; (Dih < edMax); Dih++)
514 snew(his_aa_ss[i][j][Dih], nbin+1);
520 for (Dih = 0; (Dih < edMax); Dih++)
522 snew(his_aa[Dih], rt_size+1);
523 for (i = 0; (i <= rt_size); i++)
525 snew(his_aa[Dih][i], nbin+1);
532 for (i = 0; (i < nlist); i++)
540 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
542 for (i = 0; (i < nlist); i++)
544 if (((Dih < edOmega) ) ||
545 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
546 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)))
548 make_histo(log, nf, dih[j], nbin, histmp, -M_PI, M_PI);
552 /* Assume there is only one structure, the first.
553 * Compute index in histogram.
555 /* Check the atoms to see whether their B-factors are low enough
556 * Check atoms to see their occupancy is 1.
558 bBfac = bOccup = TRUE;
559 for (nn = 0; (nn < 4); nn++, n++)
561 bBfac = bBfac && (atoms->pdbinfo[index[n]].bfac <= bfac_max);
562 bOccup = bOccup && (atoms->pdbinfo[index[n]].occup == 1);
564 if (bOccup && ((bfac_max <= 0) || ((bfac_max > 0) && bBfac)))
566 hindex = ((dih[j][0]+M_PI)*nbin)/(2*M_PI);
567 range_check(hindex, 0, nbin);
569 /* Assign dihedral to either of the structure determined
572 switch (ss_str[dlist[i].resnr])
575 his_aa_ss[0][dlist[i].index][Dih][hindex]++;
578 his_aa_ss[1][dlist[i].index][Dih][hindex]++;
581 his_aa_ss[2][dlist[i].index][Dih][hindex]++;
587 fprintf(debug, "Res. %d has imcomplete occupancy or bfacs > %g\n",
588 dlist[i].resnr, bfac_max);
599 calc_distribution_props(nbin, histmp, -M_PI, NKKKPHI, kkkphi, &S2);
601 for (m = 0; (m < NKKKPHI); m++)
603 Jc[i][m] = kkkphi[m].Jc;
604 Jcsig[i][m] = kkkphi[m].Jcsig;
608 calc_distribution_props(nbin, histmp, -M_PI, NKKKPSI, kkkpsi, &S2);
610 for (m = 0; (m < NKKKPSI); m++)
612 Jc[i][NKKKPHI+m] = kkkpsi[m].Jc;
613 Jcsig[i][NKKKPHI+m] = kkkpsi[m].Jcsig;
617 calc_distribution_props(nbin, histmp, -M_PI, NKKKCHI, kkkchi1, &S2);
618 for (m = 0; (m < NKKKCHI); m++)
620 Jc[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jc;
621 Jcsig[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jcsig;
624 default: /* covers edOmega and higher Chis than Chi1 */
625 calc_distribution_props(nbin, histmp, -M_PI, 0, NULL, &S2);
628 dlist[i].S2[Dih] = S2;
630 /* Sum distribution per amino acid type as well */
631 for (k = 0; (k < nbin); k++)
633 his_aa[Dih][dlist[i].index][k] += histmp[k];
638 else /* dihed not defined */
640 dlist[i].S2[Dih] = 0.0;
646 /* Print out Jcouplings */
647 fprintf(log, "\n *** J-Couplings from simulation (plus std. dev.) ***\n\n");
648 fprintf(log, "Residue ");
649 for (i = 0; (i < NKKKPHI); i++)
651 fprintf(log, "%7s SD", kkkphi[i].name);
653 for (i = 0; (i < NKKKPSI); i++)
655 fprintf(log, "%7s SD", kkkpsi[i].name);
657 for (i = 0; (i < NKKKCHI); i++)
659 fprintf(log, "%7s SD", kkkchi1[i].name);
662 for (i = 0; (i < NJC+1); i++)
664 fprintf(log, "------------");
667 for (i = 0; (i < nlist); i++)
669 fprintf(log, "%-10s", dlist[i].name);
670 for (j = 0; (j < NJC); j++)
672 fprintf(log, " %5.2f %4.2f", Jc[i][j], Jcsig[i][j]);
678 /* and to -jc file... */
681 fp = xvgropen(fn, "\\S3\\NJ-Couplings from Karplus Equation", "Residue",
684 for (i = 0; (i < NKKKPHI); i++)
686 leg[i] = gmx_strdup(kkkphi[i].name);
688 for (i = 0; (i < NKKKPSI); i++)
690 leg[i+NKKKPHI] = gmx_strdup(kkkpsi[i].name);
692 for (i = 0; (i < NKKKCHI); i++)
694 leg[i+NKKKPHI+NKKKPSI] = gmx_strdup(kkkchi1[i].name);
696 xvgr_legend(fp, NJC, (const char**)leg, oenv);
697 fprintf(fp, "%5s ", "#Res.");
698 for (i = 0; (i < NJC); i++)
700 fprintf(fp, "%10s ", leg[i]);
703 for (i = 0; (i < nlist); i++)
705 fprintf(fp, "%5d ", dlist[i].resnr);
706 for (j = 0; (j < NJC); j++)
708 fprintf(fp, " %8.3f", Jc[i][j]);
713 for (i = 0; (i < NJC); i++)
718 /* finished -jc stuff */
720 snew(normhisto, nbin);
721 for (i = 0; (i < rt_size); i++)
723 for (Dih = 0; (Dih < edMax); Dih++)
725 /* First check whether something is in there */
726 for (j = 0; (j < nbin); j++)
728 if (his_aa[Dih][i][j] != 0)
734 ((bPhi && (Dih == edPhi)) ||
735 (bPsi && (Dih == edPsi)) ||
736 (bOmega && (Dih == edOmega)) ||
737 (bChi && (Dih >= edChi1))))
741 normalize_histo(nbin, his_aa[Dih][i], (360.0/nbin), normhisto);
744 residue_name = gmx_residuetype_get_name(rt, i);
748 sprintf(hisfile, "histo-phi%s", residue_name);
749 sprintf(title, "\\xf\\f{} Distribution for %s", residue_name);
752 sprintf(hisfile, "histo-psi%s", residue_name);
753 sprintf(title, "\\xy\\f{} Distribution for %s", residue_name);
756 sprintf(hisfile, "histo-omega%s", residue_name);
757 sprintf(title, "\\xw\\f{} Distribution for %s", residue_name);
760 sprintf(hisfile, "histo-chi%d%s", Dih-NONCHI+1, residue_name);
761 sprintf(title, "\\xc\\f{}\\s%d\\N Distribution for %s",
762 Dih-NONCHI+1, residue_name);
764 strcpy(hhisfile, hisfile);
765 strcat(hhisfile, ".xvg");
766 fp = xvgropen(hhisfile, title, "Degrees", "", oenv);
767 if (output_env_get_print_xvgr_codes(oenv))
769 fprintf(fp, "@ with g0\n");
771 xvgr_world(fp, -180, 0, 180, 0.1, oenv);
772 if (output_env_get_print_xvgr_codes(oenv))
774 fprintf(fp, "# this effort to set graph size fails unless you run with -autoscale none or -autoscale y flags\n");
775 fprintf(fp, "@ xaxis tick on\n");
776 fprintf(fp, "@ xaxis tick major 90\n");
777 fprintf(fp, "@ xaxis tick minor 30\n");
778 fprintf(fp, "@ xaxis ticklabel prec 0\n");
779 fprintf(fp, "@ yaxis tick off\n");
780 fprintf(fp, "@ yaxis ticklabel off\n");
781 fprintf(fp, "@ type xy\n");
785 for (k = 0; (k < 3); k++)
787 sprintf(sshisfile, "%s-%s.xvg", hisfile, sss[k]);
788 ssfp[k] = gmx_ffopen(sshisfile, "w");
791 for (j = 0; (j < nbin); j++)
793 angle = -180 + (360/nbin)*j;
796 fprintf(fp, "%5d %10g\n", angle, normhisto[j]);
800 fprintf(fp, "%5d %10d\n", angle, his_aa[Dih][i][j]);
804 for (k = 0; (k < 3); k++)
806 fprintf(ssfp[k], "%5d %10d\n", angle,
807 his_aa_ss[k][i][Dih][j]);
811 fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
815 for (k = 0; (k < 3); k++)
817 fprintf(ssfp[k], "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
818 gmx_ffclose(ssfp[k]);
828 /* Four dimensional array... Very cool */
829 for (i = 0; (i < 3); i++)
831 for (j = 0; (j <= rt_size); j++)
833 for (Dih = 0; (Dih < edMax); Dih++)
835 sfree(his_aa_ss[i][j][Dih]);
837 sfree(his_aa_ss[i][j]);
846 static FILE *rama_file(const char *fn, const char *title, const char *xaxis,
847 const char *yaxis, const output_env_t oenv)
851 fp = xvgropen(fn, title, xaxis, yaxis, oenv);
852 if (output_env_get_print_xvgr_codes(oenv))
854 fprintf(fp, "@ with g0\n");
856 xvgr_world(fp, -180, -180, 180, 180, oenv);
857 if (output_env_get_print_xvgr_codes(oenv))
859 fprintf(fp, "@ xaxis tick on\n");
860 fprintf(fp, "@ xaxis tick major 90\n");
861 fprintf(fp, "@ xaxis tick minor 30\n");
862 fprintf(fp, "@ xaxis ticklabel prec 0\n");
863 fprintf(fp, "@ yaxis tick on\n");
864 fprintf(fp, "@ yaxis tick major 90\n");
865 fprintf(fp, "@ yaxis tick minor 30\n");
866 fprintf(fp, "@ yaxis ticklabel prec 0\n");
867 fprintf(fp, "@ s0 type xy\n");
868 fprintf(fp, "@ s0 symbol 2\n");
869 fprintf(fp, "@ s0 symbol size 0.410000\n");
870 fprintf(fp, "@ s0 symbol fill 1\n");
871 fprintf(fp, "@ s0 symbol color 1\n");
872 fprintf(fp, "@ s0 symbol linewidth 1\n");
873 fprintf(fp, "@ s0 symbol linestyle 1\n");
874 fprintf(fp, "@ s0 symbol center false\n");
875 fprintf(fp, "@ s0 symbol char 0\n");
876 fprintf(fp, "@ s0 skip 0\n");
877 fprintf(fp, "@ s0 linestyle 0\n");
878 fprintf(fp, "@ s0 linewidth 1\n");
879 fprintf(fp, "@ type xy\n");
884 static void do_rama(int nf, int nlist, t_dlist dlist[], real **dih,
885 gmx_bool bViol, gmx_bool bRamOmega, const output_env_t oenv)
887 FILE *fp, *gp = NULL;
890 int i, j, k, Xi1, Xi2, Phi, Psi, Om = 0, nlevels;
892 real **mat = NULL, phi, psi, omega, axis[NMAT], lo, hi;
893 t_rgb rlo = { 1.0, 0.0, 0.0 };
894 t_rgb rmid = { 1.0, 1.0, 1.0 };
895 t_rgb rhi = { 0.0, 0.0, 1.0 };
897 for (i = 0; (i < nlist); i++)
899 if ((has_dihedral(edPhi, &(dlist[i]))) &&
900 (has_dihedral(edPsi, &(dlist[i]))))
902 sprintf(fn, "ramaPhiPsi%s.xvg", dlist[i].name);
903 fp = rama_file(fn, "Ramachandran Plot",
904 "\\8f\\4 (deg)", "\\8y\\4 (deg)", oenv);
905 bOm = bRamOmega && has_dihedral(edOmega, &(dlist[i]));
908 Om = dlist[i].j0[edOmega];
910 for (j = 0; (j < NMAT); j++)
913 axis[j] = -180+(360*j)/NMAT;
918 sprintf(fn, "violPhiPsi%s.xvg", dlist[i].name);
919 gp = gmx_ffopen(fn, "w");
921 Phi = dlist[i].j0[edPhi];
922 Psi = dlist[i].j0[edPsi];
923 for (j = 0; (j < nf); j++)
925 phi = RAD2DEG*dih[Phi][j];
926 psi = RAD2DEG*dih[Psi][j];
927 fprintf(fp, "%10g %10g\n", phi, psi);
930 fprintf(gp, "%d\n", !bAllowed(dih[Phi][j], RAD2DEG*dih[Psi][j]));
934 omega = RAD2DEG*dih[Om][j];
935 mat[(int)((phi*NMAT)/360)+NMAT/2][(int)((psi*NMAT)/360)+NMAT/2]
946 sprintf(fn, "ramomega%s.xpm", dlist[i].name);
947 fp = gmx_ffopen(fn, "w");
949 for (j = 0; (j < NMAT); j++)
951 for (k = 0; (k < NMAT); k++)
954 lo = min(mat[j][k], lo);
955 hi = max(mat[j][k], hi);
959 if (fabs(lo) > fabs(hi))
968 for (j = 0; (j < NMAT); j++)
970 for (k = 0; (k < NMAT); k++)
978 write_xpm3(fp, 0, "Omega/Ramachandran Plot", "Deg", "Phi", "Psi",
979 NMAT, NMAT, axis, axis, mat, lo, 180.0, hi, rlo, rmid, rhi, &nlevels);
981 for (j = 0; (j < NMAT); j++)
988 if ((has_dihedral(edChi1, &(dlist[i]))) &&
989 (has_dihedral(edChi2, &(dlist[i]))))
991 sprintf(fn, "ramaX1X2%s.xvg", dlist[i].name);
992 fp = rama_file(fn, "\\8c\\4\\s1\\N-\\8c\\4\\s2\\N Ramachandran Plot",
993 "\\8c\\4\\s1\\N (deg)", "\\8c\\4\\s2\\N (deg)", oenv);
994 Xi1 = dlist[i].j0[edChi1];
995 Xi2 = dlist[i].j0[edChi2];
996 for (j = 0; (j < nf); j++)
998 fprintf(fp, "%10g %10g\n", RAD2DEG*dih[Xi1][j], RAD2DEG*dih[Xi2][j]);
1004 fprintf(stderr, "No chi1 & chi2 angle for %s\n", dlist[i].name);
1010 static void print_transitions(const char *fn, int maxchi, int nlist,
1011 t_dlist dlist[], real dt,
1012 const output_env_t oenv)
1014 /* based on order_params below */
1019 /* must correspond with enum in pp2shift.h:38 */
1021 #define NLEG asize(leg)
1023 leg[0] = gmx_strdup("Phi");
1024 leg[1] = gmx_strdup("Psi");
1025 leg[2] = gmx_strdup("Omega");
1026 leg[3] = gmx_strdup("Chi1");
1027 leg[4] = gmx_strdup("Chi2");
1028 leg[5] = gmx_strdup("Chi3");
1029 leg[6] = gmx_strdup("Chi4");
1030 leg[7] = gmx_strdup("Chi5");
1031 leg[8] = gmx_strdup("Chi6");
1033 /* Print order parameters */
1034 fp = xvgropen(fn, "Dihedral Rotamer Transitions", "Residue", "Transitions/ns",
1036 xvgr_legend(fp, NONCHI+maxchi, (const char**)leg, oenv);
1038 for (Dih = 0; (Dih < edMax); Dih++)
1043 fprintf(fp, "%5s ", "#Res.");
1044 fprintf(fp, "%10s %10s %10s ", leg[edPhi], leg[edPsi], leg[edOmega]);
1045 for (Xi = 0; Xi < maxchi; Xi++)
1047 fprintf(fp, "%10s ", leg[NONCHI+Xi]);
1051 for (i = 0; (i < nlist); i++)
1053 fprintf(fp, "%5d ", dlist[i].resnr);
1054 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1056 fprintf(fp, "%10.3f ", dlist[i].ntr[Dih]/dt);
1058 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
1064 static void order_params(FILE *log,
1065 const char *fn, int maxchi, int nlist, t_dlist dlist[],
1066 const char *pdbfn, real bfac_init,
1067 t_atoms *atoms, rvec x[], int ePBC, matrix box,
1068 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, const output_env_t oenv)
1075 /* except for S2Min/Max, must correspond with enum in pp2shift.h:38 */
1076 const char *const_leg[2+edMax] = {
1077 "S2Min", "S2Max", "Phi", "Psi", "Omega",
1078 "Chi1", "Chi2", "Chi3", "Chi4", "Chi5",
1081 #define NLEG asize(leg)
1085 for (i = 0; i < NLEG; i++)
1087 leg[i] = gmx_strdup(const_leg[i]);
1090 /* Print order parameters */
1091 fp = xvgropen(fn, "Dihedral Order Parameters", "Residue", "S2", oenv);
1092 xvgr_legend(fp, 2+NONCHI+maxchi, const_leg, oenv);
1094 for (Dih = 0; (Dih < edMax); Dih++)
1099 fprintf(fp, "%5s ", "#Res.");
1100 fprintf(fp, "%10s %10s ", leg[0], leg[1]);
1101 fprintf(fp, "%10s %10s %10s ", leg[2+edPhi], leg[2+edPsi], leg[2+edOmega]);
1102 for (Xi = 0; Xi < maxchi; Xi++)
1104 fprintf(fp, "%10s ", leg[2+NONCHI+Xi]);
1108 for (i = 0; (i < nlist); i++)
1112 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1114 if (dlist[i].S2[Dih] != 0)
1116 if (dlist[i].S2[Dih] > S2Max)
1118 S2Max = dlist[i].S2[Dih];
1120 if (dlist[i].S2[Dih] < S2Min)
1122 S2Min = dlist[i].S2[Dih];
1125 if (dlist[i].S2[Dih] > 0.8)
1130 fprintf(fp, "%5d ", dlist[i].resnr);
1131 fprintf(fp, "%10.3f %10.3f ", S2Min, S2Max);
1132 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1134 fprintf(fp, "%10.3f ", dlist[i].S2[Dih]);
1137 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
1145 if (NULL == atoms->pdbinfo)
1147 snew(atoms->pdbinfo, atoms->nr);
1149 for (i = 0; (i < atoms->nr); i++)
1151 atoms->pdbinfo[i].bfac = bfac_init;
1154 for (i = 0; (i < nlist); i++)
1156 atoms->pdbinfo[dlist[i].atm.N].bfac = -dlist[i].S2[0]; /* Phi */
1157 atoms->pdbinfo[dlist[i].atm.H].bfac = -dlist[i].S2[0]; /* Phi */
1158 atoms->pdbinfo[dlist[i].atm.C].bfac = -dlist[i].S2[1]; /* Psi */
1159 atoms->pdbinfo[dlist[i].atm.O].bfac = -dlist[i].S2[1]; /* Psi */
1160 for (Xi = 0; (Xi < maxchi); Xi++) /* Chi's */
1162 if (dlist[i].atm.Cn[Xi+3] != -1)
1164 atoms->pdbinfo[dlist[i].atm.Cn[Xi+1]].bfac = -dlist[i].S2[NONCHI+Xi];
1169 fp = gmx_ffopen(pdbfn, "w");
1170 fprintf(fp, "REMARK generated by g_chi\n");
1171 fprintf(fp, "REMARK "
1172 "B-factor field contains negative of dihedral order parameters\n");
1173 write_pdbfile(fp, NULL, atoms, x, ePBC, box, ' ', 0, NULL, TRUE);
1174 x0 = y0 = z0 = 1000.0;
1175 for (i = 0; (i < atoms->nr); i++)
1177 x0 = min(x0, x[i][XX]);
1178 y0 = min(y0, x[i][YY]);
1179 z0 = min(z0, x[i][ZZ]);
1181 x0 *= 10.0; /* nm -> angstrom */
1182 y0 *= 10.0; /* nm -> angstrom */
1183 z0 *= 10.0; /* nm -> angstrom */
1184 for (i = 0; (i < 10); i++)
1186 gmx_fprintf_pdb_atomline(fp, epdbATOM, atoms->nr+1+i, "CA", ' ', "LEG", ' ', atoms->nres+1, ' ',
1187 x0, y0, z0+(1.2*i), 0.0, -0.1*i, "");
1192 fprintf(log, "Dihedrals with S2 > 0.8\n");
1193 fprintf(log, "Dihedral: ");
1196 fprintf(log, " Phi ");
1200 fprintf(log, " Psi ");
1204 for (Xi = 0; (Xi < maxchi); Xi++)
1206 fprintf(log, " %s ", leg[2+NONCHI+Xi]);
1209 fprintf(log, "\nNumber: ");
1212 fprintf(log, "%4d ", nh[0]);
1216 fprintf(log, "%4d ", nh[1]);
1220 for (Xi = 0; (Xi < maxchi); Xi++)
1222 fprintf(log, "%4d ", nh[NONCHI+Xi]);
1227 for (i = 0; i < NLEG; i++)
1234 int gmx_chi(int argc, char *argv[])
1236 const char *desc[] = {
1237 "[THISMODULE] computes [GRK]phi[grk], [GRK]psi[grk], [GRK]omega[grk],",
1238 "and [GRK]chi[grk] dihedrals for all your",
1239 "amino acid backbone and sidechains.",
1240 "It can compute dihedral angle as a function of time, and as",
1241 "histogram distributions.",
1242 "The distributions [TT](histo-(dihedral)(RESIDUE).xvg[tt]) are cumulative over all residues of each type.[PAR]",
1243 "If option [TT]-corr[tt] is given, the program will",
1244 "calculate dihedral autocorrelation functions. The function used",
1245 "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",
1246 "rather than angles themselves, resolves the problem of periodicity.",
1247 "(Van der Spoel & Berendsen (1997), Biophys. J. 72, 2032-2041).",
1248 "Separate files for each dihedral of each residue",
1249 "[TT](corr(dihedral)(RESIDUE)(nresnr).xvg[tt]) are output, as well as a",
1250 "file containing the information for all residues (argument of [TT]-corr[tt]).[PAR]",
1251 "With option [TT]-all[tt], the angles themselves as a function of time for",
1252 "each residue are printed to separate files [TT](dihedral)(RESIDUE)(nresnr).xvg[tt].",
1253 "These can be in radians or degrees.[PAR]",
1254 "A log file (argument [TT]-g[tt]) is also written. This contains [BR]",
1255 "(a) information about the number of residues of each type.[BR]",
1256 "(b) The NMR ^3J coupling constants from the Karplus equation.[BR]",
1257 "(c) a table for each residue of the number of transitions between ",
1258 "rotamers per nanosecond, and the order parameter S^2 of each dihedral.[BR]",
1259 "(d) a table for each residue of the rotamer occupancy.[PAR]",
1260 "All rotamers are taken as 3-fold, except for [GRK]omega[grk] and [GRK]chi[grk] dihedrals",
1261 "to planar groups (i.e. [GRK]chi[grk][SUB]2[sub] of aromatics, Asp and Asn; [GRK]chi[grk][SUB]3[sub] of Glu",
1262 "and Gln; and [GRK]chi[grk][SUB]4[sub] of Arg), which are 2-fold. \"rotamer 0\" means ",
1263 "that the dihedral was not in the core region of each rotamer. ",
1264 "The width of the core region can be set with [TT]-core_rotamer[tt][PAR]",
1266 "The S^2 order parameters are also output to an [TT].xvg[tt] file",
1267 "(argument [TT]-o[tt] ) and optionally as a [TT].pdb[tt] file with",
1268 "the S^2 values as B-factor (argument [TT]-p[tt]). ",
1269 "The total number of rotamer transitions per timestep",
1270 "(argument [TT]-ot[tt]), the number of transitions per rotamer",
1271 "(argument [TT]-rt[tt]), and the ^3J couplings (argument [TT]-jc[tt]), ",
1272 "can also be written to [TT].xvg[tt] files. Note that the analysis",
1273 "of rotamer transitions assumes that the supplied trajectory frames",
1274 "are equally spaced in time.[PAR]",
1276 "If [TT]-chi_prod[tt] is set (and [TT]-maxchi[tt] > 0), cumulative rotamers, e.g.",
1277 "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 ",
1278 "dihedrals and [TT]-maxchi[tt] >= 3)",
1279 "are calculated. As before, if any dihedral is not in the core region,",
1280 "the rotamer is taken to be 0. The occupancies of these cumulative ",
1281 "rotamers (starting with rotamer 0) are written to the file",
1282 "that is the argument of [TT]-cp[tt], and if the [TT]-all[tt] flag",
1283 "is given, the rotamers as functions of time",
1284 "are written to [TT]chiproduct(RESIDUE)(nresnr).xvg[tt] ",
1285 "and their occupancies to [TT]histo-chiproduct(RESIDUE)(nresnr).xvg[tt].[PAR]",
1287 "The option [TT]-r[tt] generates a contour plot of the average [GRK]omega[grk] angle",
1288 "as a function of the [GRK]phi[grk] and [GRK]psi[grk] angles, that is, in a Ramachandran plot",
1289 "the average [GRK]omega[grk] angle is plotted using color coding.",
1293 const char *bugs[] = {
1294 "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.",
1295 "[GRK]phi[grk] and [GRK]psi[grk] dihedrals are calculated in a "
1296 "non-standard way, using H-N-CA-C for [GRK]phi[grk] instead of "
1297 "C(-)-N-CA-C, and N-CA-C-O for [GRK]psi[grk] instead of N-CA-C-N(+). "
1298 "This causes (usually small) discrepancies with the output of other "
1299 "tools like [gmx-rama].",
1300 "[TT]-r0[tt] option does not work properly",
1301 "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"
1305 static int r0 = 1, ndeg = 1, maxchi = 2;
1306 static gmx_bool bAll = FALSE;
1307 static gmx_bool bPhi = FALSE, bPsi = FALSE, bOmega = FALSE;
1308 static real bfac_init = -1.0, bfac_max = 0;
1309 static const char *maxchistr[] = { NULL, "0", "1", "2", "3", "4", "5", "6", NULL };
1310 static gmx_bool bRama = FALSE, bShift = FALSE, bViol = FALSE, bRamOmega = FALSE;
1311 static gmx_bool bNormHisto = TRUE, bChiProduct = FALSE, bHChi = FALSE, bRAD = FALSE, bPBC = TRUE;
1312 static real core_frac = 0.5;
1314 { "-r0", FALSE, etINT, {&r0},
1315 "starting residue" },
1316 { "-phi", FALSE, etBOOL, {&bPhi},
1317 "Output for [GRK]phi[grk] dihedral angles" },
1318 { "-psi", FALSE, etBOOL, {&bPsi},
1319 "Output for [GRK]psi[grk] dihedral angles" },
1320 { "-omega", FALSE, etBOOL, {&bOmega},
1321 "Output for [GRK]omega[grk] dihedrals (peptide bonds)" },
1322 { "-rama", FALSE, etBOOL, {&bRama},
1323 "Generate [GRK]phi[grk]/[GRK]psi[grk] and [GRK]chi[grk][SUB]1[sub]/[GRK]chi[grk][SUB]2[sub] Ramachandran plots" },
1324 { "-viol", FALSE, etBOOL, {&bViol},
1325 "Write a file that gives 0 or 1 for violated Ramachandran angles" },
1326 { "-periodic", FALSE, etBOOL, {&bPBC},
1327 "Print dihedral angles modulo 360 degrees" },
1328 { "-all", FALSE, etBOOL, {&bAll},
1329 "Output separate files for every dihedral." },
1330 { "-rad", FALSE, etBOOL, {&bRAD},
1331 "in angle vs time files, use radians rather than degrees."},
1332 { "-shift", FALSE, etBOOL, {&bShift},
1333 "Compute chemical shifts from [GRK]phi[grk]/[GRK]psi[grk] angles" },
1334 { "-binwidth", FALSE, etINT, {&ndeg},
1335 "bin width for histograms (degrees)" },
1336 { "-core_rotamer", FALSE, etREAL, {&core_frac},
1337 "only the central [TT]-core_rotamer[tt]*(360/multiplicity) belongs to each rotamer (the rest is assigned to rotamer 0)" },
1338 { "-maxchi", FALSE, etENUM, {maxchistr},
1339 "calculate first ndih [GRK]chi[grk] dihedrals" },
1340 { "-normhisto", FALSE, etBOOL, {&bNormHisto},
1341 "Normalize histograms" },
1342 { "-ramomega", FALSE, etBOOL, {&bRamOmega},
1343 "compute average omega as a function of [GRK]phi[grk]/[GRK]psi[grk] and plot it in an [TT].xpm[tt] plot" },
1344 { "-bfact", FALSE, etREAL, {&bfac_init},
1345 "B-factor value for [TT].pdb[tt] file for atoms with no calculated dihedral order parameter"},
1346 { "-chi_prod", FALSE, etBOOL, {&bChiProduct},
1347 "compute a single cumulative rotamer for each residue"},
1348 { "-HChi", FALSE, etBOOL, {&bHChi},
1349 "Include dihedrals to sidechain hydrogens"},
1350 { "-bmax", FALSE, etREAL, {&bfac_max},
1351 "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." }
1355 int natoms, nlist, idum, nbin;
1360 char title[256], grpname[256];
1362 gmx_bool bChi, bCorr, bSSHisto;
1363 gmx_bool bDo_rt, bDo_oh, bDo_ot, bDo_jc;
1364 real dt = 0, traj_t_ns;
1366 gmx_residuetype_t *rt;
1368 atom_id isize, *index;
1369 int ndih, nactdih, nf;
1370 real **dih, *trans_frac, *aver_angle, *time;
1371 int i, j, **chi_lookup, *multiplicity;
1374 { efSTX, "-s", NULL, ffREAD },
1375 { efTRX, "-f", NULL, ffREAD },
1376 { efXVG, "-o", "order", ffWRITE },
1377 { efPDB, "-p", "order", ffOPTWR },
1378 { efDAT, "-ss", "ssdump", ffOPTRD },
1379 { efXVG, "-jc", "Jcoupling", ffWRITE },
1380 { efXVG, "-corr", "dihcorr", ffOPTWR },
1381 { efLOG, "-g", "chi", ffWRITE },
1382 /* add two more arguments copying from g_angle */
1383 { efXVG, "-ot", "dihtrans", ffOPTWR },
1384 { efXVG, "-oh", "trhisto", ffOPTWR },
1385 { efXVG, "-rt", "restrans", ffOPTWR },
1386 { efXVG, "-cp", "chiprodhisto", ffOPTWR }
1388 #define NFILE asize(fnm)
1393 ppa = add_acf_pargs(&npargs, pa);
1394 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
1395 NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs,
1401 /* Handle result from enumerated type */
1402 sscanf(maxchistr[0], "%d", &maxchi);
1403 bChi = (maxchi > 0);
1405 log = gmx_ffopen(ftp2fn(efLOG, NFILE, fnm), "w");
1414 /* set some options */
1415 bDo_rt = (opt2bSet("-rt", NFILE, fnm));
1416 bDo_oh = (opt2bSet("-oh", NFILE, fnm));
1417 bDo_ot = (opt2bSet("-ot", NFILE, fnm));
1418 bDo_jc = (opt2bSet("-jc", NFILE, fnm));
1419 bCorr = (opt2bSet("-corr", NFILE, fnm));
1422 fprintf(stderr, "Will calculate autocorrelation\n");
1425 if (core_frac > 1.0)
1427 fprintf(stderr, "core_rotamer fraction > 1.0 ; will use 1.0\n");
1430 if (core_frac < 0.0)
1432 fprintf(stderr, "core_rotamer fraction < 0.0 ; will use 0.0\n");
1436 if (maxchi > MAXCHI)
1439 "Will only calculate first %d Chi dihedrals in stead of %d.\n",
1443 bSSHisto = ftp2bSet(efDAT, NFILE, fnm);
1446 /* Find the chi angles using atoms struct and a list of amino acids */
1447 get_stx_coordnum(ftp2fn(efSTX, NFILE, fnm), &natoms);
1448 init_t_atoms(&atoms, natoms, TRUE);
1450 read_stx_conf(ftp2fn(efSTX, NFILE, fnm), title, &atoms, x, NULL, &ePBC, box);
1451 fprintf(log, "Title: %s\n", title);
1453 gmx_residuetype_init(&rt);
1454 dlist = mk_dlist(log, &atoms, &nlist, bPhi, bPsi, bChi, bHChi, maxchi, r0, rt);
1455 fprintf(stderr, "%d residues with dihedrals found\n", nlist);
1459 gmx_fatal(FARGS, "No dihedrals in your structure!\n");
1462 /* Make a linear index for reading all. */
1463 index = make_chi_ind(nlist, dlist, &ndih);
1465 fprintf(stderr, "%d dihedrals found\n", ndih);
1469 /* COMPUTE ALL DIHEDRALS! */
1470 read_ang_dih(ftp2fn(efTRX, NFILE, fnm), FALSE, TRUE, FALSE, bPBC, 1, &idum,
1471 &nf, &time, isize, index, &trans_frac, &aver_angle, dih, oenv);
1473 dt = (time[nf-1]-time[0])/(nf-1); /* might want this for corr or n. transit*/
1478 gmx_fatal(FARGS, "Need at least 2 frames for correlation");
1482 /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi
1483 * pass nactdih instead of ndih to low_ana_dih_trans and get_chi_product_traj
1484 * to prevent accessing off end of arrays when maxchi < 5 or 6. */
1485 nactdih = reset_em_all(nlist, dlist, nf, dih, maxchi);
1489 dump_em_all(nlist, dlist, nf, time, dih, maxchi, bPhi, bPsi, bChi, bOmega, bRAD, oenv);
1492 /* Histogramming & J coupling constants & calc of S2 order params */
1493 histogramming(log, nbin, rt, nf, maxchi, dih, nlist, dlist, index,
1494 bPhi, bPsi, bOmega, bChi,
1495 bNormHisto, bSSHisto, ftp2fn(efDAT, NFILE, fnm), bfac_max, &atoms,
1496 bDo_jc, opt2fn("-jc", NFILE, fnm), oenv);
1500 * added multiplicity */
1502 snew(multiplicity, ndih);
1503 mk_multiplicity_lookup(multiplicity, maxchi, nlist, dlist, ndih);
1505 strcpy(grpname, "All residues, ");
1508 strcat(grpname, "Phi ");
1512 strcat(grpname, "Psi ");
1516 strcat(grpname, "Omega ");
1520 strcat(grpname, "Chi 1-");
1521 sprintf(grpname + strlen(grpname), "%i", maxchi);
1525 low_ana_dih_trans(bDo_ot, opt2fn("-ot", NFILE, fnm),
1526 bDo_oh, opt2fn("-oh", NFILE, fnm), maxchi,
1527 dih, nlist, dlist, nf, nactdih, grpname, multiplicity,
1528 time, FALSE, core_frac, oenv);
1530 /* Order parameters */
1531 order_params(log, opt2fn("-o", NFILE, fnm), maxchi, nlist, dlist,
1532 ftp2fn_null(efPDB, NFILE, fnm), bfac_init,
1533 &atoms, x, ePBC, box, bPhi, bPsi, bChi, oenv);
1535 /* Print ramachandran maps! */
1538 do_rama(nf, nlist, dlist, dih, bViol, bRamOmega, oenv);
1543 do_pp2shifts(log, nf, nlist, dlist, dih);
1546 /* rprint S^2, transitions, and rotamer occupancies to log */
1547 traj_t_ns = 0.001 * (time[nf-1]-time[0]);
1548 pr_dlist(log, nlist, dlist, traj_t_ns, edPrintST, bPhi, bPsi, bChi, bOmega, maxchi);
1549 pr_dlist(log, nlist, dlist, traj_t_ns, edPrintRO, bPhi, bPsi, bChi, bOmega, maxchi);
1551 /* transitions to xvg */
1554 print_transitions(opt2fn("-rt", NFILE, fnm), maxchi, nlist, dlist, traj_t_ns, oenv);
1557 /* chi_product trajectories (ie one "rotamer number" for each residue) */
1558 if (bChiProduct && bChi)
1560 snew(chi_lookup, nlist);
1561 for (i = 0; i < nlist; i++)
1563 snew(chi_lookup[i], maxchi);
1565 mk_chi_lookup(chi_lookup, maxchi, nlist, dlist);
1567 get_chi_product_traj(dih, nf, nactdih,
1568 maxchi, dlist, time, chi_lookup, multiplicity,
1569 FALSE, bNormHisto, core_frac, bAll,
1570 opt2fn("-cp", NFILE, fnm), oenv);
1572 for (i = 0; i < nlist; i++)
1574 sfree(chi_lookup[i]);
1578 /* Correlation comes last because it messes up the angles */
1581 do_dihcorr(opt2fn("-corr", NFILE, fnm), nf, ndih, dih, dt, nlist, dlist, time,
1582 maxchi, bPhi, bPsi, bChi, bOmega, oenv);
1586 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1587 do_view(oenv, opt2fn("-jc", NFILE, fnm), "-nxy");
1590 do_view(oenv, opt2fn("-corr", NFILE, fnm), "-nxy");
1593 gmx_residuetype_destroy(rt);