2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
43 #include "gromacs/fileio/confio.h"
44 #include "gromacs/fileio/pdbio.h"
45 #include "gromacs/utility/fatalerror.h"
46 #include "gromacs/utility/futil.h"
49 #include "gromacs/math/utilities.h"
50 #include "gromacs/math/units.h"
51 #include "gromacs/topology/residuetypes.h"
52 #include "gromacs/utility/smalloc.h"
53 #include "gromacs/commandline/pargs.h"
54 #include "gromacs/fileio/tpxio.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/fileio/xvgr.h"
60 #include "gromacs/fileio/matio.h"
63 static gmx_bool bAllowed(real phi, real psi)
65 static const char *map[] = {
66 "1100000000000000001111111000000000001111111111111111111111111",
67 "1100000000000000001111110000000000011111111111111111111111111",
68 "1100000000000000001111110000000000011111111111111111111111111",
69 "1100000000000000001111100000000000111111111111111111111111111",
70 "1100000000000000001111100000000000111111111111111111111111111",
71 "1100000000000000001111100000000001111111111111111111111111111",
72 "1100000000000000001111100000000001111111111111111111111111111",
73 "1100000000000000001111100000000011111111111111111111111111111",
74 "1110000000000000001111110000000111111111111111111111111111111",
75 "1110000000000000001111110000001111111111111111111111111111111",
76 "1110000000000000001111111000011111111111111111111111111111111",
77 "1110000000000000001111111100111111111111111111111111111111111",
78 "1110000000000000001111111111111111111111111111111111111111111",
79 "1110000000000000001111111111111111111111111111111111111111111",
80 "1110000000000000001111111111111111111111111111111111111111111",
81 "1110000000000000001111111111111111111111111111111111111111111",
82 "1110000000000000001111111111111110011111111111111111111111111",
83 "1110000000000000001111111111111100000111111111111111111111111",
84 "1110000000000000001111111111111000000000001111111111111111111",
85 "1100000000000000001111111111110000000000000011111111111111111",
86 "1100000000000000001111111111100000000000000011111111111111111",
87 "1000000000000000001111111111000000000000000001111111111111110",
88 "0000000000000000001111111110000000000000000000111111111111100",
89 "0000000000000000000000000000000000000000000000000000000000000",
90 "0000000000000000000000000000000000000000000000000000000000000",
91 "0000000000000000000000000000000000000000000000000000000000000",
92 "0000000000000000000000000000000000000000000000000000000000000",
93 "0000000000000000000000000000000000000000000000000000000000000",
94 "0000000000000000000000000000000000000000000000000000000000000",
95 "0000000000000000000000000000000000000000000000000000000000000",
96 "0000000000000000000000000000000000000000000000000000000000000",
97 "0000000000000000000000000000000000000000000000000000000000000",
98 "0000000000000000000000000000000000000000000000000000000000000",
99 "0000000000000000000000000000000000000000000000000000000000000",
100 "0000000000000000000000000000000000000000000000000000000000000",
101 "0000000000000000000000000000000000000000000000000000000000000",
102 "0000000000000000000000000000000000000000000000000000000000000",
103 "0000000000000000000000000000000000000000000000000000000000000",
104 "0000000000000000000000000000000000111111111111000000000000000",
105 "1100000000000000000000000000000001111111111111100000000000111",
106 "1100000000000000000000000000000001111111111111110000000000111",
107 "0000000000000000000000000000000000000000000000000000000000000",
108 "0000000000000000000000000000000000000000000000000000000000000",
109 "0000000000000000000000000000000000000000000000000000000000000",
110 "0000000000000000000000000000000000000000000000000000000000000",
111 "0000000000000000000000000000000000000000000000000000000000000",
112 "0000000000000000000000000000000000000000000000000000000000000",
113 "0000000000000000000000000000000000000000000000000000000000000",
114 "0000000000000000000000000000000000000000000000000000000000000",
115 "0000000000000000000000000000000000000000000000000000000000000",
116 "0000000000000000000000000000000000000000000000000000000000000",
117 "0000000000000000000000000000000000000000000000000000000000000",
118 "0000000000000000000000000000000000000000000000000000000000000",
119 "0000000000000000000000000000000000000000000000000000000000000",
120 "0000000000000000000000000000000000000000000000000000000000000",
121 "0000000000000000000000000000000000000000000000000000000000000",
122 "0000000000000000000000000000000000000000000000000000000000000",
123 "0000000000000000000000000000000000000000000000000000000000000",
124 "0000000000000000000000000000000000000000000000000000000000000",
125 "0000000000000000000000000000000000000000000000000000000000000",
126 "0000000000000000000000000000000000000000000000000000000000000"
128 #define NPP asize(map)
131 #define INDEX(ppp) ((((int) (360+ppp*RAD2DEG)) % 360)/6)
135 return (gmx_bool) map[x][y];
138 atom_id *make_chi_ind(int nl, t_dlist dl[], int *ndih)
143 /* There are nl residues with max edMax dihedrals with 4 atoms each */
144 snew(id, nl*edMax*4);
147 for (i = 0; (i < nl); i++)
149 /* Phi, fake the first one */
150 dl[i].j0[edPhi] = n/4;
151 if (dl[i].atm.minC >= 0)
153 id[n++] = dl[i].atm.minC;
157 id[n++] = dl[i].atm.H;
159 id[n++] = dl[i].atm.N;
160 id[n++] = dl[i].atm.Cn[1];
161 id[n++] = dl[i].atm.C;
163 for (i = 0; (i < nl); i++)
165 /* Psi, fake the last one */
166 dl[i].j0[edPsi] = n/4;
167 id[n++] = dl[i].atm.N;
168 id[n++] = dl[i].atm.Cn[1];
169 id[n++] = dl[i].atm.C;
172 id[n++] = dl[i+1].atm.N;
176 id[n++] = dl[i].atm.O;
179 for (i = 1; (i < nl); i++)
182 if (has_dihedral(edOmega, &(dl[i])))
184 dl[i].j0[edOmega] = n/4;
185 id[n++] = dl[i].atm.minCalpha;
186 id[n++] = dl[i].atm.minC;
187 id[n++] = dl[i].atm.N;
188 id[n++] = dl[i].atm.Cn[1];
191 for (Xi = 0; (Xi < MAXCHI); Xi++)
194 for (i = 0; (i < nl); i++)
196 if (dl[i].atm.Cn[Xi+3] != -1)
198 dl[i].j0[edChi1+Xi] = n/4;
199 id[n++] = dl[i].atm.Cn[Xi];
200 id[n++] = dl[i].atm.Cn[Xi+1];
201 id[n++] = dl[i].atm.Cn[Xi+2];
202 id[n++] = dl[i].atm.Cn[Xi+3];
211 int bin(real chi, int mult)
215 return (int) (chi*mult/360.0);
219 static void do_dihcorr(const char *fn, int nf, int ndih, real **dih, real dt,
220 int nlist, t_dlist dlist[], real time[], int maxchi,
221 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bOmega,
222 const output_env_t oenv)
224 char name1[256], name2[256];
227 do_autocorr(fn, oenv, "Dihedral Autocorrelation Function",
228 nf, ndih, dih, dt, eacCos, FALSE);
231 for (i = 0; (i < nlist); i++)
235 print_one(oenv, "corrphi", dlist[i].name, "Phi ACF for", "C(t)", nf/2, time,
240 for (i = 0; (i < nlist); i++)
244 print_one(oenv, "corrpsi", dlist[i].name, "Psi ACF for", "C(t)", nf/2, time,
249 for (i = 0; (i < nlist); i++)
251 if (has_dihedral(edOmega, &dlist[i]))
255 print_one(oenv, "corromega", dlist[i].name, "Omega ACF for", "C(t)",
261 for (Xi = 0; (Xi < maxchi); Xi++)
263 sprintf(name1, "corrchi%d", Xi+1);
264 sprintf(name2, "Chi%d ACF for", Xi+1);
265 for (i = 0; (i < nlist); i++)
267 if (dlist[i].atm.Cn[Xi+3] != -1)
271 print_one(oenv, name1, dlist[i].name, name2, "C(t)", nf/2, time, dih[j]);
277 fprintf(stderr, "\n");
280 static void copy_dih_data(real in[], real out[], int nf, gmx_bool bLEAVE)
282 /* if bLEAVE, do nothing to data in copying to out
283 * otherwise multiply by 180/pi to convert rad to deg */
294 for (i = 0; (i < nf); i++)
300 static void dump_em_all(int nlist, t_dlist dlist[], int nf, real time[],
301 real **dih, int maxchi,
302 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bOmega, gmx_bool bRAD,
303 const output_env_t oenv)
305 char name[256], titlestr[256], ystr[256];
312 strcpy(ystr, "Angle (rad)");
316 strcpy(ystr, "Angle (degrees)");
321 for (i = 0; (i < nlist); i++)
323 /* grs debug printf("OK i %d j %d\n", i, j) ; */
326 copy_dih_data(dih[j], data, nf, bRAD);
327 print_one(oenv, "phi", dlist[i].name, "\\xf\\f{}", ystr, nf, time, data);
331 for (i = 0; (i < nlist); i++)
335 copy_dih_data(dih[j], data, nf, bRAD);
336 print_one(oenv, "psi", dlist[i].name, "\\xy\\f{}", ystr, nf, time, data);
340 for (i = 0; (i < nlist); i++)
342 if (has_dihedral(edOmega, &(dlist[i])))
346 copy_dih_data(dih[j], data, nf, bRAD);
347 print_one(oenv, "omega", dlist[i].name, "\\xw\\f{}", ystr, nf, time, data);
353 for (Xi = 0; (Xi < maxchi); Xi++)
355 for (i = 0; (i < nlist); i++)
357 if (dlist[i].atm.Cn[Xi+3] != -1)
361 sprintf(name, "chi%d", Xi+1);
362 sprintf(titlestr, "\\xc\\f{}\\s%d\\N", Xi+1);
363 copy_dih_data(dih[j], data, nf, bRAD);
364 print_one(oenv, name, dlist[i].name, titlestr, ystr, nf, time, data);
370 fprintf(stderr, "\n");
373 static void reset_one(real dih[], int nf, real phase)
377 for (j = 0; (j < nf); j++)
380 while (dih[j] < -M_PI)
384 while (dih[j] >= M_PI)
391 static int reset_em_all(int nlist, t_dlist dlist[], int nf,
392 real **dih, int maxchi)
399 for (i = 0; (i < nlist); i++)
401 if (dlist[i].atm.minC == -1)
403 reset_one(dih[j++], nf, M_PI);
407 reset_one(dih[j++], nf, 0);
411 for (i = 0; (i < nlist-1); i++)
413 reset_one(dih[j++], nf, 0);
415 /* last Psi is faked from O */
416 reset_one(dih[j++], nf, M_PI);
419 for (i = 0; (i < nlist); i++)
421 if (has_dihedral(edOmega, &dlist[i]))
423 reset_one(dih[j++], nf, 0);
426 /* Chi 1 thru maxchi */
427 for (Xi = 0; (Xi < maxchi); Xi++)
429 for (i = 0; (i < nlist); i++)
431 if (dlist[i].atm.Cn[Xi+3] != -1)
433 reset_one(dih[j], nf, 0);
438 fprintf(stderr, "j after resetting (nr. active dihedrals) = %d\n", j);
442 static void histogramming(FILE *log, int nbin, gmx_residuetype_t *rt,
443 int nf, int maxchi, real **dih,
444 int nlist, t_dlist dlist[],
446 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bOmega, gmx_bool bChi,
447 gmx_bool bNormalize, gmx_bool bSSHisto, const char *ssdump,
448 real bfac_max, t_atoms *atoms,
449 gmx_bool bDo_jc, const char *fn,
450 const output_env_t oenv)
452 /* also gets 3J couplings and order parameters S2 */
453 t_karplus kkkphi[] = {
454 { "J_NHa1", 6.51, -1.76, 1.6, -M_PI/3, 0.0, 0.0 },
455 { "J_NHa2", 6.51, -1.76, 1.6, M_PI/3, 0.0, 0.0 },
456 { "J_HaC'", 4.0, 1.1, 0.1, 0.0, 0.0, 0.0 },
457 { "J_NHCb", 4.7, -1.5, -0.2, M_PI/3, 0.0, 0.0 },
458 { "J_Ci-1Hai", 4.5, -1.3, -1.2, 2*M_PI/3, 0.0, 0.0 }
460 t_karplus kkkpsi[] = {
461 { "J_HaN", -0.88, -0.61, -0.27, M_PI/3, 0.0, 0.0 }
463 t_karplus kkkchi1[] = {
464 { "JHaHb2", 9.5, -1.6, 1.8, -M_PI/3, 0, 0.0 },
465 { "JHaHb3", 9.5, -1.6, 1.8, 0, 0, 0.0 }
467 #define NKKKPHI asize(kkkphi)
468 #define NKKKPSI asize(kkkpsi)
469 #define NKKKCHI asize(kkkchi1)
470 #define NJC (NKKKPHI+NKKKPSI+NKKKCHI)
472 FILE *fp, *ssfp[3] = {NULL, NULL, NULL};
473 const char *sss[3] = { "sheet", "helix", "coil" };
477 int ****his_aa_ss = NULL;
478 int ***his_aa, **his_aa1, *histmp;
479 int i, j, k, m, n, nn, Dih, nres, hindex, angle;
480 gmx_bool bBfac, bOccup;
481 char hisfile[256], hhisfile[256], sshisfile[256], title[256], *ss_str = NULL;
483 const char *residue_name;
486 rt_size = gmx_residuetype_get_size(rt);
489 fp = gmx_ffopen(ssdump, "r");
490 if (1 != fscanf(fp, "%d", &nres))
492 gmx_fatal(FARGS, "Error reading from file %s", ssdump);
495 snew(ss_str, nres+1);
496 if (1 != fscanf(fp, "%s", ss_str))
498 gmx_fatal(FARGS, "Error reading from file %s", ssdump);
502 /* Four dimensional array... Very cool */
504 for (i = 0; (i < 3); i++)
506 snew(his_aa_ss[i], rt_size+1);
507 for (j = 0; (j <= rt_size); j++)
509 snew(his_aa_ss[i][j], edMax);
510 for (Dih = 0; (Dih < edMax); Dih++)
512 snew(his_aa_ss[i][j][Dih], nbin+1);
518 for (Dih = 0; (Dih < edMax); Dih++)
520 snew(his_aa[Dih], rt_size+1);
521 for (i = 0; (i <= rt_size); i++)
523 snew(his_aa[Dih][i], nbin+1);
530 for (i = 0; (i < nlist); i++)
538 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
540 for (i = 0; (i < nlist); i++)
542 if (((Dih < edOmega) ) ||
543 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
544 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)))
546 make_histo(log, nf, dih[j], nbin, histmp, -M_PI, M_PI);
550 /* Assume there is only one structure, the first.
551 * Compute index in histogram.
553 /* Check the atoms to see whether their B-factors are low enough
554 * Check atoms to see their occupancy is 1.
556 bBfac = bOccup = TRUE;
557 for (nn = 0; (nn < 4); nn++, n++)
559 bBfac = bBfac && (atoms->pdbinfo[index[n]].bfac <= bfac_max);
560 bOccup = bOccup && (atoms->pdbinfo[index[n]].occup == 1);
562 if (bOccup && ((bfac_max <= 0) || ((bfac_max > 0) && bBfac)))
564 hindex = ((dih[j][0]+M_PI)*nbin)/(2*M_PI);
565 range_check(hindex, 0, nbin);
567 /* Assign dihedral to either of the structure determined
570 switch (ss_str[dlist[i].resnr])
573 his_aa_ss[0][dlist[i].index][Dih][hindex]++;
576 his_aa_ss[1][dlist[i].index][Dih][hindex]++;
579 his_aa_ss[2][dlist[i].index][Dih][hindex]++;
585 fprintf(debug, "Res. %d has imcomplete occupancy or bfacs > %g\n",
586 dlist[i].resnr, bfac_max);
597 calc_distribution_props(nbin, histmp, -M_PI, NKKKPHI, kkkphi, &S2);
599 for (m = 0; (m < NKKKPHI); m++)
601 Jc[i][m] = kkkphi[m].Jc;
602 Jcsig[i][m] = kkkphi[m].Jcsig;
606 calc_distribution_props(nbin, histmp, -M_PI, NKKKPSI, kkkpsi, &S2);
608 for (m = 0; (m < NKKKPSI); m++)
610 Jc[i][NKKKPHI+m] = kkkpsi[m].Jc;
611 Jcsig[i][NKKKPHI+m] = kkkpsi[m].Jcsig;
615 calc_distribution_props(nbin, histmp, -M_PI, NKKKCHI, kkkchi1, &S2);
616 for (m = 0; (m < NKKKCHI); m++)
618 Jc[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jc;
619 Jcsig[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jcsig;
622 default: /* covers edOmega and higher Chis than Chi1 */
623 calc_distribution_props(nbin, histmp, -M_PI, 0, NULL, &S2);
626 dlist[i].S2[Dih] = S2;
628 /* Sum distribution per amino acid type as well */
629 for (k = 0; (k < nbin); k++)
631 his_aa[Dih][dlist[i].index][k] += histmp[k];
636 else /* dihed not defined */
638 dlist[i].S2[Dih] = 0.0;
644 /* Print out Jcouplings */
645 fprintf(log, "\n *** J-Couplings from simulation (plus std. dev.) ***\n\n");
646 fprintf(log, "Residue ");
647 for (i = 0; (i < NKKKPHI); i++)
649 fprintf(log, "%7s SD", kkkphi[i].name);
651 for (i = 0; (i < NKKKPSI); i++)
653 fprintf(log, "%7s SD", kkkpsi[i].name);
655 for (i = 0; (i < NKKKCHI); i++)
657 fprintf(log, "%7s SD", kkkchi1[i].name);
660 for (i = 0; (i < NJC+1); i++)
662 fprintf(log, "------------");
665 for (i = 0; (i < nlist); i++)
667 fprintf(log, "%-10s", dlist[i].name);
668 for (j = 0; (j < NJC); j++)
670 fprintf(log, " %5.2f %4.2f", Jc[i][j], Jcsig[i][j]);
676 /* and to -jc file... */
679 fp = xvgropen(fn, "\\S3\\NJ-Couplings from Karplus Equation", "Residue",
682 for (i = 0; (i < NKKKPHI); i++)
684 leg[i] = strdup(kkkphi[i].name);
686 for (i = 0; (i < NKKKPSI); i++)
688 leg[i+NKKKPHI] = strdup(kkkpsi[i].name);
690 for (i = 0; (i < NKKKCHI); i++)
692 leg[i+NKKKPHI+NKKKPSI] = strdup(kkkchi1[i].name);
694 xvgr_legend(fp, NJC, (const char**)leg, oenv);
695 fprintf(fp, "%5s ", "#Res.");
696 for (i = 0; (i < NJC); i++)
698 fprintf(fp, "%10s ", leg[i]);
701 for (i = 0; (i < nlist); i++)
703 fprintf(fp, "%5d ", dlist[i].resnr);
704 for (j = 0; (j < NJC); j++)
706 fprintf(fp, " %8.3f", Jc[i][j]);
711 for (i = 0; (i < NJC); i++)
716 /* finished -jc stuff */
718 snew(normhisto, nbin);
719 for (i = 0; (i < rt_size); i++)
721 for (Dih = 0; (Dih < edMax); Dih++)
723 /* First check whether something is in there */
724 for (j = 0; (j < nbin); j++)
726 if (his_aa[Dih][i][j] != 0)
732 ((bPhi && (Dih == edPhi)) ||
733 (bPsi && (Dih == edPsi)) ||
734 (bOmega && (Dih == edOmega)) ||
735 (bChi && (Dih >= edChi1))))
739 normalize_histo(nbin, his_aa[Dih][i], (360.0/nbin), normhisto);
742 residue_name = gmx_residuetype_get_name(rt, i);
746 sprintf(hisfile, "histo-phi%s", residue_name);
747 sprintf(title, "\\xf\\f{} Distribution for %s", residue_name);
750 sprintf(hisfile, "histo-psi%s", residue_name);
751 sprintf(title, "\\xy\\f{} Distribution for %s", residue_name);
754 sprintf(hisfile, "histo-omega%s", residue_name);
755 sprintf(title, "\\xw\\f{} Distribution for %s", residue_name);
758 sprintf(hisfile, "histo-chi%d%s", Dih-NONCHI+1, residue_name);
759 sprintf(title, "\\xc\\f{}\\s%d\\N Distribution for %s",
760 Dih-NONCHI+1, residue_name);
762 strcpy(hhisfile, hisfile);
763 strcat(hhisfile, ".xvg");
764 fp = xvgropen(hhisfile, title, "Degrees", "", oenv);
765 if (output_env_get_print_xvgr_codes(oenv))
767 fprintf(fp, "@ with g0\n");
769 xvgr_world(fp, -180, 0, 180, 0.1, oenv);
770 if (output_env_get_print_xvgr_codes(oenv))
772 fprintf(fp, "# this effort to set graph size fails unless you run with -autoscale none or -autoscale y flags\n");
773 fprintf(fp, "@ xaxis tick on\n");
774 fprintf(fp, "@ xaxis tick major 90\n");
775 fprintf(fp, "@ xaxis tick minor 30\n");
776 fprintf(fp, "@ xaxis ticklabel prec 0\n");
777 fprintf(fp, "@ yaxis tick off\n");
778 fprintf(fp, "@ yaxis ticklabel off\n");
779 fprintf(fp, "@ type xy\n");
783 for (k = 0; (k < 3); k++)
785 sprintf(sshisfile, "%s-%s.xvg", hisfile, sss[k]);
786 ssfp[k] = gmx_ffopen(sshisfile, "w");
789 for (j = 0; (j < nbin); j++)
791 angle = -180 + (360/nbin)*j;
794 fprintf(fp, "%5d %10g\n", angle, normhisto[j]);
798 fprintf(fp, "%5d %10d\n", angle, his_aa[Dih][i][j]);
802 for (k = 0; (k < 3); k++)
804 fprintf(ssfp[k], "%5d %10d\n", angle,
805 his_aa_ss[k][i][Dih][j]);
809 fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
813 for (k = 0; (k < 3); k++)
815 fprintf(ssfp[k], "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
816 gmx_ffclose(ssfp[k]);
826 /* Four dimensional array... Very cool */
827 for (i = 0; (i < 3); i++)
829 for (j = 0; (j <= rt_size); j++)
831 for (Dih = 0; (Dih < edMax); Dih++)
833 sfree(his_aa_ss[i][j][Dih]);
835 sfree(his_aa_ss[i][j]);
844 static FILE *rama_file(const char *fn, const char *title, const char *xaxis,
845 const char *yaxis, const output_env_t oenv)
849 fp = xvgropen(fn, title, xaxis, yaxis, oenv);
850 if (output_env_get_print_xvgr_codes(oenv))
852 fprintf(fp, "@ with g0\n");
854 xvgr_world(fp, -180, -180, 180, 180, oenv);
855 if (output_env_get_print_xvgr_codes(oenv))
857 fprintf(fp, "@ xaxis tick on\n");
858 fprintf(fp, "@ xaxis tick major 90\n");
859 fprintf(fp, "@ xaxis tick minor 30\n");
860 fprintf(fp, "@ xaxis ticklabel prec 0\n");
861 fprintf(fp, "@ yaxis tick on\n");
862 fprintf(fp, "@ yaxis tick major 90\n");
863 fprintf(fp, "@ yaxis tick minor 30\n");
864 fprintf(fp, "@ yaxis ticklabel prec 0\n");
865 fprintf(fp, "@ s0 type xy\n");
866 fprintf(fp, "@ s0 symbol 2\n");
867 fprintf(fp, "@ s0 symbol size 0.410000\n");
868 fprintf(fp, "@ s0 symbol fill 1\n");
869 fprintf(fp, "@ s0 symbol color 1\n");
870 fprintf(fp, "@ s0 symbol linewidth 1\n");
871 fprintf(fp, "@ s0 symbol linestyle 1\n");
872 fprintf(fp, "@ s0 symbol center false\n");
873 fprintf(fp, "@ s0 symbol char 0\n");
874 fprintf(fp, "@ s0 skip 0\n");
875 fprintf(fp, "@ s0 linestyle 0\n");
876 fprintf(fp, "@ s0 linewidth 1\n");
877 fprintf(fp, "@ type xy\n");
882 static void do_rama(int nf, int nlist, t_dlist dlist[], real **dih,
883 gmx_bool bViol, gmx_bool bRamOmega, const output_env_t oenv)
885 FILE *fp, *gp = NULL;
888 int i, j, k, Xi1, Xi2, Phi, Psi, Om = 0, nlevels;
890 real **mat = NULL, phi, psi, omega, axis[NMAT], lo, hi;
891 t_rgb rlo = { 1.0, 0.0, 0.0 };
892 t_rgb rmid = { 1.0, 1.0, 1.0 };
893 t_rgb rhi = { 0.0, 0.0, 1.0 };
895 for (i = 0; (i < nlist); i++)
897 if ((has_dihedral(edPhi, &(dlist[i]))) &&
898 (has_dihedral(edPsi, &(dlist[i]))))
900 sprintf(fn, "ramaPhiPsi%s.xvg", dlist[i].name);
901 fp = rama_file(fn, "Ramachandran Plot",
902 "\\8f\\4 (deg)", "\\8y\\4 (deg)", oenv);
903 bOm = bRamOmega && has_dihedral(edOmega, &(dlist[i]));
906 Om = dlist[i].j0[edOmega];
908 for (j = 0; (j < NMAT); j++)
911 axis[j] = -180+(360*j)/NMAT;
916 sprintf(fn, "violPhiPsi%s.xvg", dlist[i].name);
917 gp = gmx_ffopen(fn, "w");
919 Phi = dlist[i].j0[edPhi];
920 Psi = dlist[i].j0[edPsi];
921 for (j = 0; (j < nf); j++)
923 phi = RAD2DEG*dih[Phi][j];
924 psi = RAD2DEG*dih[Psi][j];
925 fprintf(fp, "%10g %10g\n", phi, psi);
928 fprintf(gp, "%d\n", !bAllowed(dih[Phi][j], RAD2DEG*dih[Psi][j]));
932 omega = RAD2DEG*dih[Om][j];
933 mat[(int)((phi*NMAT)/360)+NMAT/2][(int)((psi*NMAT)/360)+NMAT/2]
944 sprintf(fn, "ramomega%s.xpm", dlist[i].name);
945 fp = gmx_ffopen(fn, "w");
947 for (j = 0; (j < NMAT); j++)
949 for (k = 0; (k < NMAT); k++)
952 lo = min(mat[j][k], lo);
953 hi = max(mat[j][k], hi);
957 if (fabs(lo) > fabs(hi))
966 for (j = 0; (j < NMAT); j++)
968 for (k = 0; (k < NMAT); k++)
976 write_xpm3(fp, 0, "Omega/Ramachandran Plot", "Deg", "Phi", "Psi",
977 NMAT, NMAT, axis, axis, mat, lo, 180.0, hi, rlo, rmid, rhi, &nlevels);
979 for (j = 0; (j < NMAT); j++)
986 if ((has_dihedral(edChi1, &(dlist[i]))) &&
987 (has_dihedral(edChi2, &(dlist[i]))))
989 sprintf(fn, "ramaX1X2%s.xvg", dlist[i].name);
990 fp = rama_file(fn, "\\8c\\4\\s1\\N-\\8c\\4\\s2\\N Ramachandran Plot",
991 "\\8c\\4\\s1\\N (deg)", "\\8c\\4\\s2\\N (deg)", oenv);
992 Xi1 = dlist[i].j0[edChi1];
993 Xi2 = dlist[i].j0[edChi2];
994 for (j = 0; (j < nf); j++)
996 fprintf(fp, "%10g %10g\n", RAD2DEG*dih[Xi1][j], RAD2DEG*dih[Xi2][j]);
1002 fprintf(stderr, "No chi1 & chi2 angle for %s\n", dlist[i].name);
1008 static void print_transitions(const char *fn, int maxchi, int nlist,
1009 t_dlist dlist[], real dt,
1010 const output_env_t oenv)
1012 /* based on order_params below */
1017 /* must correspond with enum in pp2shift.h:38 */
1019 #define NLEG asize(leg)
1021 leg[0] = strdup("Phi");
1022 leg[1] = strdup("Psi");
1023 leg[2] = strdup("Omega");
1024 leg[3] = strdup("Chi1");
1025 leg[4] = strdup("Chi2");
1026 leg[5] = strdup("Chi3");
1027 leg[6] = strdup("Chi4");
1028 leg[7] = strdup("Chi5");
1029 leg[8] = strdup("Chi6");
1031 /* Print order parameters */
1032 fp = xvgropen(fn, "Dihedral Rotamer Transitions", "Residue", "Transitions/ns",
1034 xvgr_legend(fp, NONCHI+maxchi, (const char**)leg, oenv);
1036 for (Dih = 0; (Dih < edMax); Dih++)
1041 fprintf(fp, "%5s ", "#Res.");
1042 fprintf(fp, "%10s %10s %10s ", leg[edPhi], leg[edPsi], leg[edOmega]);
1043 for (Xi = 0; Xi < maxchi; Xi++)
1045 fprintf(fp, "%10s ", leg[NONCHI+Xi]);
1049 for (i = 0; (i < nlist); i++)
1051 fprintf(fp, "%5d ", dlist[i].resnr);
1052 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1054 fprintf(fp, "%10.3f ", dlist[i].ntr[Dih]/dt);
1056 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
1062 static void order_params(FILE *log,
1063 const char *fn, int maxchi, int nlist, t_dlist dlist[],
1064 const char *pdbfn, real bfac_init,
1065 t_atoms *atoms, rvec x[], int ePBC, matrix box,
1066 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, const output_env_t oenv)
1073 /* except for S2Min/Max, must correspond with enum in pp2shift.h:38 */
1074 const char *const_leg[2+edMax] = {
1075 "S2Min", "S2Max", "Phi", "Psi", "Omega",
1076 "Chi1", "Chi2", "Chi3", "Chi4", "Chi5",
1079 #define NLEG asize(leg)
1083 for (i = 0; i < NLEG; i++)
1085 leg[i] = strdup(const_leg[i]);
1088 /* Print order parameters */
1089 fp = xvgropen(fn, "Dihedral Order Parameters", "Residue", "S2", oenv);
1090 xvgr_legend(fp, 2+NONCHI+maxchi, const_leg, oenv);
1092 for (Dih = 0; (Dih < edMax); Dih++)
1097 fprintf(fp, "%5s ", "#Res.");
1098 fprintf(fp, "%10s %10s ", leg[0], leg[1]);
1099 fprintf(fp, "%10s %10s %10s ", leg[2+edPhi], leg[2+edPsi], leg[2+edOmega]);
1100 for (Xi = 0; Xi < maxchi; Xi++)
1102 fprintf(fp, "%10s ", leg[2+NONCHI+Xi]);
1106 for (i = 0; (i < nlist); i++)
1110 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1112 if (dlist[i].S2[Dih] != 0)
1114 if (dlist[i].S2[Dih] > S2Max)
1116 S2Max = dlist[i].S2[Dih];
1118 if (dlist[i].S2[Dih] < S2Min)
1120 S2Min = dlist[i].S2[Dih];
1123 if (dlist[i].S2[Dih] > 0.8)
1128 fprintf(fp, "%5d ", dlist[i].resnr);
1129 fprintf(fp, "%10.3f %10.3f ", S2Min, S2Max);
1130 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1132 fprintf(fp, "%10.3f ", dlist[i].S2[Dih]);
1135 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
1143 if (NULL == atoms->pdbinfo)
1145 snew(atoms->pdbinfo, atoms->nr);
1147 for (i = 0; (i < atoms->nr); i++)
1149 atoms->pdbinfo[i].bfac = bfac_init;
1152 for (i = 0; (i < nlist); i++)
1154 atoms->pdbinfo[dlist[i].atm.N].bfac = -dlist[i].S2[0]; /* Phi */
1155 atoms->pdbinfo[dlist[i].atm.H].bfac = -dlist[i].S2[0]; /* Phi */
1156 atoms->pdbinfo[dlist[i].atm.C].bfac = -dlist[i].S2[1]; /* Psi */
1157 atoms->pdbinfo[dlist[i].atm.O].bfac = -dlist[i].S2[1]; /* Psi */
1158 for (Xi = 0; (Xi < maxchi); Xi++) /* Chi's */
1160 if (dlist[i].atm.Cn[Xi+3] != -1)
1162 atoms->pdbinfo[dlist[i].atm.Cn[Xi+1]].bfac = -dlist[i].S2[NONCHI+Xi];
1167 fp = gmx_ffopen(pdbfn, "w");
1168 fprintf(fp, "REMARK generated by g_chi\n");
1169 fprintf(fp, "REMARK "
1170 "B-factor field contains negative of dihedral order parameters\n");
1171 write_pdbfile(fp, NULL, atoms, x, ePBC, box, ' ', 0, NULL, TRUE);
1172 x0 = y0 = z0 = 1000.0;
1173 for (i = 0; (i < atoms->nr); i++)
1175 x0 = min(x0, x[i][XX]);
1176 y0 = min(y0, x[i][YY]);
1177 z0 = min(z0, x[i][ZZ]);
1179 x0 *= 10.0; /* nm -> angstrom */
1180 y0 *= 10.0; /* nm -> angstrom */
1181 z0 *= 10.0; /* nm -> angstrom */
1182 for (i = 0; (i < 10); i++)
1184 gmx_fprintf_pdb_atomline(fp, epdbATOM, atoms->nr+1+i, "CA", ' ', "LEG", ' ', atoms->nres+1, ' ',
1185 x0, y0, z0+(1.2*i), 0.0, -0.1*i, "");
1190 fprintf(log, "Dihedrals with S2 > 0.8\n");
1191 fprintf(log, "Dihedral: ");
1194 fprintf(log, " Phi ");
1198 fprintf(log, " Psi ");
1202 for (Xi = 0; (Xi < maxchi); Xi++)
1204 fprintf(log, " %s ", leg[2+NONCHI+Xi]);
1207 fprintf(log, "\nNumber: ");
1210 fprintf(log, "%4d ", nh[0]);
1214 fprintf(log, "%4d ", nh[1]);
1218 for (Xi = 0; (Xi < maxchi); Xi++)
1220 fprintf(log, "%4d ", nh[NONCHI+Xi]);
1225 for (i = 0; i < NLEG; i++)
1232 int gmx_chi(int argc, char *argv[])
1234 const char *desc[] = {
1235 "[THISMODULE] computes [GRK]phi[grk], [GRK]psi[grk], [GRK]omega[grk],",
1236 "and [GRK]chi[grk] dihedrals for all your",
1237 "amino acid backbone and sidechains.",
1238 "It can compute dihedral angle as a function of time, and as",
1239 "histogram distributions.",
1240 "The distributions [TT](histo-(dihedral)(RESIDUE).xvg[tt]) are cumulative over all residues of each type.[PAR]",
1241 "If option [TT]-corr[tt] is given, the program will",
1242 "calculate dihedral autocorrelation functions. The function used",
1243 "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",
1244 "rather than angles themselves, resolves the problem of periodicity.",
1245 "(Van der Spoel & Berendsen (1997), Biophys. J. 72, 2032-2041).",
1246 "Separate files for each dihedral of each residue",
1247 "[TT](corr(dihedral)(RESIDUE)(nresnr).xvg[tt]) are output, as well as a",
1248 "file containing the information for all residues (argument of [TT]-corr[tt]).[PAR]",
1249 "With option [TT]-all[tt], the angles themselves as a function of time for",
1250 "each residue are printed to separate files [TT](dihedral)(RESIDUE)(nresnr).xvg[tt].",
1251 "These can be in radians or degrees.[PAR]",
1252 "A log file (argument [TT]-g[tt]) is also written. This contains [BR]",
1253 "(a) information about the number of residues of each type.[BR]",
1254 "(b) The NMR ^3J coupling constants from the Karplus equation.[BR]",
1255 "(c) a table for each residue of the number of transitions between ",
1256 "rotamers per nanosecond, and the order parameter S^2 of each dihedral.[BR]",
1257 "(d) a table for each residue of the rotamer occupancy.[PAR]",
1258 "All rotamers are taken as 3-fold, except for [GRK]omega[grk] and [GRK]chi[grk] dihedrals",
1259 "to planar groups (i.e. [GRK]chi[grk][SUB]2[sub] of aromatics, Asp and Asn; [GRK]chi[grk][SUB]3[sub] of Glu",
1260 "and Gln; and [GRK]chi[grk][SUB]4[sub] of Arg), which are 2-fold. \"rotamer 0\" means ",
1261 "that the dihedral was not in the core region of each rotamer. ",
1262 "The width of the core region can be set with [TT]-core_rotamer[tt][PAR]",
1264 "The S^2 order parameters are also output to an [TT].xvg[tt] file",
1265 "(argument [TT]-o[tt] ) and optionally as a [TT].pdb[tt] file with",
1266 "the S^2 values as B-factor (argument [TT]-p[tt]). ",
1267 "The total number of rotamer transitions per timestep",
1268 "(argument [TT]-ot[tt]), the number of transitions per rotamer",
1269 "(argument [TT]-rt[tt]), and the ^3J couplings (argument [TT]-jc[tt]), ",
1270 "can also be written to [TT].xvg[tt] files. Note that the analysis",
1271 "of rotamer transitions assumes that the supplied trajectory frames",
1272 "are equally spaced in time.[PAR]",
1274 "If [TT]-chi_prod[tt] is set (and [TT]-maxchi[tt] > 0), cumulative rotamers, e.g.",
1275 "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 ",
1276 "dihedrals and [TT]-maxchi[tt] >= 3)",
1277 "are calculated. As before, if any dihedral is not in the core region,",
1278 "the rotamer is taken to be 0. The occupancies of these cumulative ",
1279 "rotamers (starting with rotamer 0) are written to the file",
1280 "that is the argument of [TT]-cp[tt], and if the [TT]-all[tt] flag",
1281 "is given, the rotamers as functions of time",
1282 "are written to [TT]chiproduct(RESIDUE)(nresnr).xvg[tt] ",
1283 "and their occupancies to [TT]histo-chiproduct(RESIDUE)(nresnr).xvg[tt].[PAR]",
1285 "The option [TT]-r[tt] generates a contour plot of the average [GRK]omega[grk] angle",
1286 "as a function of the [GRK]phi[grk] and [GRK]psi[grk] angles, that is, in a Ramachandran plot",
1287 "the average [GRK]omega[grk] angle is plotted using color coding.",
1291 const char *bugs[] = {
1292 "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.",
1293 "[GRK]phi[grk] and [GRK]psi[grk] dihedrals are calculated in a "
1294 "non-standard way, using H-N-CA-C for [GRK]phi[grk] instead of "
1295 "C(-)-N-CA-C, and N-CA-C-O for [GRK]psi[grk] instead of N-CA-C-N(+). "
1296 "This causes (usually small) discrepancies with the output of other "
1297 "tools like [gmx-rama].",
1298 "[TT]-r0[tt] option does not work properly",
1299 "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"
1303 static int r0 = 1, ndeg = 1, maxchi = 2;
1304 static gmx_bool bAll = FALSE;
1305 static gmx_bool bPhi = FALSE, bPsi = FALSE, bOmega = FALSE;
1306 static real bfac_init = -1.0, bfac_max = 0;
1307 static const char *maxchistr[] = { NULL, "0", "1", "2", "3", "4", "5", "6", NULL };
1308 static gmx_bool bRama = FALSE, bShift = FALSE, bViol = FALSE, bRamOmega = FALSE;
1309 static gmx_bool bNormHisto = TRUE, bChiProduct = FALSE, bHChi = FALSE, bRAD = FALSE, bPBC = TRUE;
1310 static real core_frac = 0.5;
1312 { "-r0", FALSE, etINT, {&r0},
1313 "starting residue" },
1314 { "-phi", FALSE, etBOOL, {&bPhi},
1315 "Output for [GRK]phi[grk] dihedral angles" },
1316 { "-psi", FALSE, etBOOL, {&bPsi},
1317 "Output for [GRK]psi[grk] dihedral angles" },
1318 { "-omega", FALSE, etBOOL, {&bOmega},
1319 "Output for [GRK]omega[grk] dihedrals (peptide bonds)" },
1320 { "-rama", FALSE, etBOOL, {&bRama},
1321 "Generate [GRK]phi[grk]/[GRK]psi[grk] and [GRK]chi[grk][SUB]1[sub]/[GRK]chi[grk][SUB]2[sub] Ramachandran plots" },
1322 { "-viol", FALSE, etBOOL, {&bViol},
1323 "Write a file that gives 0 or 1 for violated Ramachandran angles" },
1324 { "-periodic", FALSE, etBOOL, {&bPBC},
1325 "Print dihedral angles modulo 360 degrees" },
1326 { "-all", FALSE, etBOOL, {&bAll},
1327 "Output separate files for every dihedral." },
1328 { "-rad", FALSE, etBOOL, {&bRAD},
1329 "in angle vs time files, use radians rather than degrees."},
1330 { "-shift", FALSE, etBOOL, {&bShift},
1331 "Compute chemical shifts from [GRK]phi[grk]/[GRK]psi[grk] angles" },
1332 { "-binwidth", FALSE, etINT, {&ndeg},
1333 "bin width for histograms (degrees)" },
1334 { "-core_rotamer", FALSE, etREAL, {&core_frac},
1335 "only the central [TT]-core_rotamer[tt]*(360/multiplicity) belongs to each rotamer (the rest is assigned to rotamer 0)" },
1336 { "-maxchi", FALSE, etENUM, {maxchistr},
1337 "calculate first ndih [GRK]chi[grk] dihedrals" },
1338 { "-normhisto", FALSE, etBOOL, {&bNormHisto},
1339 "Normalize histograms" },
1340 { "-ramomega", FALSE, etBOOL, {&bRamOmega},
1341 "compute average omega as a function of [GRK]phi[grk]/[GRK]psi[grk] and plot it in an [TT].xpm[tt] plot" },
1342 { "-bfact", FALSE, etREAL, {&bfac_init},
1343 "B-factor value for [TT].pdb[tt] file for atoms with no calculated dihedral order parameter"},
1344 { "-chi_prod", FALSE, etBOOL, {&bChiProduct},
1345 "compute a single cumulative rotamer for each residue"},
1346 { "-HChi", FALSE, etBOOL, {&bHChi},
1347 "Include dihedrals to sidechain hydrogens"},
1348 { "-bmax", FALSE, etREAL, {&bfac_max},
1349 "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." }
1353 int natoms, nlist, idum, nbin;
1358 char title[256], grpname[256];
1360 gmx_bool bChi, bCorr, bSSHisto;
1361 gmx_bool bDo_rt, bDo_oh, bDo_ot, bDo_jc;
1362 real dt = 0, traj_t_ns;
1364 gmx_residuetype_t *rt;
1366 atom_id isize, *index;
1367 int ndih, nactdih, nf;
1368 real **dih, *trans_frac, *aver_angle, *time;
1369 int i, j, **chi_lookup, *multiplicity;
1372 { efSTX, "-s", NULL, ffREAD },
1373 { efTRX, "-f", NULL, ffREAD },
1374 { efXVG, "-o", "order", ffWRITE },
1375 { efPDB, "-p", "order", ffOPTWR },
1376 { efDAT, "-ss", "ssdump", ffOPTRD },
1377 { efXVG, "-jc", "Jcoupling", ffWRITE },
1378 { efXVG, "-corr", "dihcorr", ffOPTWR },
1379 { efLOG, "-g", "chi", ffWRITE },
1380 /* add two more arguments copying from g_angle */
1381 { efXVG, "-ot", "dihtrans", ffOPTWR },
1382 { efXVG, "-oh", "trhisto", ffOPTWR },
1383 { efXVG, "-rt", "restrans", ffOPTWR },
1384 { efXVG, "-cp", "chiprodhisto", ffOPTWR }
1386 #define NFILE asize(fnm)
1391 ppa = add_acf_pargs(&npargs, pa);
1392 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
1393 NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs,
1399 /* Handle result from enumerated type */
1400 sscanf(maxchistr[0], "%d", &maxchi);
1401 bChi = (maxchi > 0);
1403 log = gmx_ffopen(ftp2fn(efLOG, NFILE, fnm), "w");
1412 /* set some options */
1413 bDo_rt = (opt2bSet("-rt", NFILE, fnm));
1414 bDo_oh = (opt2bSet("-oh", NFILE, fnm));
1415 bDo_ot = (opt2bSet("-ot", NFILE, fnm));
1416 bDo_jc = (opt2bSet("-jc", NFILE, fnm));
1417 bCorr = (opt2bSet("-corr", NFILE, fnm));
1420 fprintf(stderr, "Will calculate autocorrelation\n");
1423 if (core_frac > 1.0)
1425 fprintf(stderr, "core_rotamer fraction > 1.0 ; will use 1.0\n");
1428 if (core_frac < 0.0)
1430 fprintf(stderr, "core_rotamer fraction < 0.0 ; will use 0.0\n");
1434 if (maxchi > MAXCHI)
1437 "Will only calculate first %d Chi dihedrals in stead of %d.\n",
1441 bSSHisto = ftp2bSet(efDAT, NFILE, fnm);
1444 /* Find the chi angles using atoms struct and a list of amino acids */
1445 get_stx_coordnum(ftp2fn(efSTX, NFILE, fnm), &natoms);
1446 init_t_atoms(&atoms, natoms, TRUE);
1448 read_stx_conf(ftp2fn(efSTX, NFILE, fnm), title, &atoms, x, NULL, &ePBC, box);
1449 fprintf(log, "Title: %s\n", title);
1451 gmx_residuetype_init(&rt);
1452 dlist = mk_dlist(log, &atoms, &nlist, bPhi, bPsi, bChi, bHChi, maxchi, r0, rt);
1453 fprintf(stderr, "%d residues with dihedrals found\n", nlist);
1457 gmx_fatal(FARGS, "No dihedrals in your structure!\n");
1460 /* Make a linear index for reading all. */
1461 index = make_chi_ind(nlist, dlist, &ndih);
1463 fprintf(stderr, "%d dihedrals found\n", ndih);
1467 /* COMPUTE ALL DIHEDRALS! */
1468 read_ang_dih(ftp2fn(efTRX, NFILE, fnm), FALSE, TRUE, FALSE, bPBC, 1, &idum,
1469 &nf, &time, isize, index, &trans_frac, &aver_angle, dih, oenv);
1471 dt = (time[nf-1]-time[0])/(nf-1); /* might want this for corr or n. transit*/
1476 gmx_fatal(FARGS, "Need at least 2 frames for correlation");
1480 /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi
1481 * pass nactdih instead of ndih to low_ana_dih_trans and get_chi_product_traj
1482 * to prevent accessing off end of arrays when maxchi < 5 or 6. */
1483 nactdih = reset_em_all(nlist, dlist, nf, dih, maxchi);
1487 dump_em_all(nlist, dlist, nf, time, dih, maxchi, bPhi, bPsi, bChi, bOmega, bRAD, oenv);
1490 /* Histogramming & J coupling constants & calc of S2 order params */
1491 histogramming(log, nbin, rt, nf, maxchi, dih, nlist, dlist, index,
1492 bPhi, bPsi, bOmega, bChi,
1493 bNormHisto, bSSHisto, ftp2fn(efDAT, NFILE, fnm), bfac_max, &atoms,
1494 bDo_jc, opt2fn("-jc", NFILE, fnm), oenv);
1498 * added multiplicity */
1500 snew(multiplicity, ndih);
1501 mk_multiplicity_lookup(multiplicity, maxchi, nlist, dlist, ndih);
1503 strcpy(grpname, "All residues, ");
1506 strcat(grpname, "Phi ");
1510 strcat(grpname, "Psi ");
1514 strcat(grpname, "Omega ");
1518 strcat(grpname, "Chi 1-");
1519 sprintf(grpname + strlen(grpname), "%i", maxchi);
1523 low_ana_dih_trans(bDo_ot, opt2fn("-ot", NFILE, fnm),
1524 bDo_oh, opt2fn("-oh", NFILE, fnm), maxchi,
1525 dih, nlist, dlist, nf, nactdih, grpname, multiplicity,
1526 time, FALSE, core_frac, oenv);
1528 /* Order parameters */
1529 order_params(log, opt2fn("-o", NFILE, fnm), maxchi, nlist, dlist,
1530 ftp2fn_null(efPDB, NFILE, fnm), bfac_init,
1531 &atoms, x, ePBC, box, bPhi, bPsi, bChi, oenv);
1533 /* Print ramachandran maps! */
1536 do_rama(nf, nlist, dlist, dih, bViol, bRamOmega, oenv);
1541 do_pp2shifts(log, nf, nlist, dlist, dih);
1544 /* rprint S^2, transitions, and rotamer occupancies to log */
1545 traj_t_ns = 0.001 * (time[nf-1]-time[0]);
1546 pr_dlist(log, nlist, dlist, traj_t_ns, edPrintST, bPhi, bPsi, bChi, bOmega, maxchi);
1547 pr_dlist(log, nlist, dlist, traj_t_ns, edPrintRO, bPhi, bPsi, bChi, bOmega, maxchi);
1549 /* transitions to xvg */
1552 print_transitions(opt2fn("-rt", NFILE, fnm), maxchi, nlist, dlist, traj_t_ns, oenv);
1555 /* chi_product trajectories (ie one "rotamer number" for each residue) */
1556 if (bChiProduct && bChi)
1558 snew(chi_lookup, nlist);
1559 for (i = 0; i < nlist; i++)
1561 snew(chi_lookup[i], maxchi);
1563 mk_chi_lookup(chi_lookup, maxchi, nlist, dlist);
1565 get_chi_product_traj(dih, nf, nactdih,
1566 maxchi, dlist, time, chi_lookup, multiplicity,
1567 FALSE, bNormHisto, core_frac, bAll,
1568 opt2fn("-cp", NFILE, fnm), oenv);
1570 for (i = 0; i < nlist; i++)
1572 sfree(chi_lookup[i]);
1576 /* Correlation comes last because it messes up the angles */
1579 do_dihcorr(opt2fn("-corr", NFILE, fnm), nf, ndih, dih, dt, nlist, dlist, time,
1580 maxchi, bPhi, bPsi, bChi, bOmega, oenv);
1584 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1585 do_view(oenv, opt2fn("-jc", NFILE, fnm), "-nxy");
1588 do_view(oenv, opt2fn("-corr", NFILE, fnm), "-nxy");
1591 gmx_residuetype_destroy(rt);