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 "gmx_fatal.h"
46 #include "gromacs/fileio/futil.h"
49 #include "gromacs/math/utilities.h"
52 #include "gromacs/utility/smalloc.h"
53 #include "gromacs/commandline/pargs.h"
54 #include "gromacs/fileio/tpxio.h"
61 #include "gromacs/fileio/matio.h"
64 static gmx_bool bAllowed(real phi, real psi)
66 static const char *map[] = {
67 "1100000000000000001111111000000000001111111111111111111111111",
68 "1100000000000000001111110000000000011111111111111111111111111",
69 "1100000000000000001111110000000000011111111111111111111111111",
70 "1100000000000000001111100000000000111111111111111111111111111",
71 "1100000000000000001111100000000000111111111111111111111111111",
72 "1100000000000000001111100000000001111111111111111111111111111",
73 "1100000000000000001111100000000001111111111111111111111111111",
74 "1100000000000000001111100000000011111111111111111111111111111",
75 "1110000000000000001111110000000111111111111111111111111111111",
76 "1110000000000000001111110000001111111111111111111111111111111",
77 "1110000000000000001111111000011111111111111111111111111111111",
78 "1110000000000000001111111100111111111111111111111111111111111",
79 "1110000000000000001111111111111111111111111111111111111111111",
80 "1110000000000000001111111111111111111111111111111111111111111",
81 "1110000000000000001111111111111111111111111111111111111111111",
82 "1110000000000000001111111111111111111111111111111111111111111",
83 "1110000000000000001111111111111110011111111111111111111111111",
84 "1110000000000000001111111111111100000111111111111111111111111",
85 "1110000000000000001111111111111000000000001111111111111111111",
86 "1100000000000000001111111111110000000000000011111111111111111",
87 "1100000000000000001111111111100000000000000011111111111111111",
88 "1000000000000000001111111111000000000000000001111111111111110",
89 "0000000000000000001111111110000000000000000000111111111111100",
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 "0000000000000000000000000000000000000000000000000000000000000",
105 "0000000000000000000000000000000000111111111111000000000000000",
106 "1100000000000000000000000000000001111111111111100000000000111",
107 "1100000000000000000000000000000001111111111111110000000000111",
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",
127 "0000000000000000000000000000000000000000000000000000000000000"
129 #define NPP asize(map)
132 #define INDEX(ppp) ((((int) (360+ppp*RAD2DEG)) % 360)/6)
136 return (gmx_bool) map[x][y];
139 atom_id *make_chi_ind(int nl, t_dlist dl[], int *ndih)
144 /* There are nl residues with max edMax dihedrals with 4 atoms each */
145 snew(id, nl*edMax*4);
148 for (i = 0; (i < nl); i++)
150 /* Phi, fake the first one */
151 dl[i].j0[edPhi] = n/4;
152 if (dl[i].atm.minC >= 0)
154 id[n++] = dl[i].atm.minC;
158 id[n++] = dl[i].atm.H;
160 id[n++] = dl[i].atm.N;
161 id[n++] = dl[i].atm.Cn[1];
162 id[n++] = dl[i].atm.C;
164 for (i = 0; (i < nl); i++)
166 /* Psi, fake the last one */
167 dl[i].j0[edPsi] = n/4;
168 id[n++] = dl[i].atm.N;
169 id[n++] = dl[i].atm.Cn[1];
170 id[n++] = dl[i].atm.C;
173 id[n++] = dl[i+1].atm.N;
177 id[n++] = dl[i].atm.O;
180 for (i = 1; (i < nl); i++)
183 if (has_dihedral(edOmega, &(dl[i])))
185 dl[i].j0[edOmega] = n/4;
186 id[n++] = dl[i].atm.minCalpha;
187 id[n++] = dl[i].atm.minC;
188 id[n++] = dl[i].atm.N;
189 id[n++] = dl[i].atm.Cn[1];
192 for (Xi = 0; (Xi < MAXCHI); Xi++)
195 for (i = 0; (i < nl); i++)
197 if (dl[i].atm.Cn[Xi+3] != -1)
199 dl[i].j0[edChi1+Xi] = n/4;
200 id[n++] = dl[i].atm.Cn[Xi];
201 id[n++] = dl[i].atm.Cn[Xi+1];
202 id[n++] = dl[i].atm.Cn[Xi+2];
203 id[n++] = dl[i].atm.Cn[Xi+3];
212 int bin(real chi, int mult)
216 return (int) (chi*mult/360.0);
220 static void do_dihcorr(const char *fn, int nf, int ndih, real **dih, real dt,
221 int nlist, t_dlist dlist[], real time[], int maxchi,
222 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bOmega,
223 const output_env_t oenv)
225 char name1[256], name2[256];
228 do_autocorr(fn, oenv, "Dihedral Autocorrelation Function",
229 nf, ndih, dih, dt, eacCos, FALSE);
232 for (i = 0; (i < nlist); i++)
236 print_one(oenv, "corrphi", dlist[i].name, "Phi ACF for", "C(t)", nf/2, time,
241 for (i = 0; (i < nlist); i++)
245 print_one(oenv, "corrpsi", dlist[i].name, "Psi ACF for", "C(t)", nf/2, time,
250 for (i = 0; (i < nlist); i++)
252 if (has_dihedral(edOmega, &dlist[i]))
256 print_one(oenv, "corromega", dlist[i].name, "Omega ACF for", "C(t)",
262 for (Xi = 0; (Xi < maxchi); Xi++)
264 sprintf(name1, "corrchi%d", Xi+1);
265 sprintf(name2, "Chi%d ACF for", Xi+1);
266 for (i = 0; (i < nlist); i++)
268 if (dlist[i].atm.Cn[Xi+3] != -1)
272 print_one(oenv, name1, dlist[i].name, name2, "C(t)", nf/2, time, dih[j]);
278 fprintf(stderr, "\n");
281 static void copy_dih_data(real in[], real out[], int nf, gmx_bool bLEAVE)
283 /* if bLEAVE, do nothing to data in copying to out
284 * otherwise multiply by 180/pi to convert rad to deg */
295 for (i = 0; (i < nf); i++)
301 static void dump_em_all(int nlist, t_dlist dlist[], int nf, real time[],
302 real **dih, int maxchi,
303 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bOmega, gmx_bool bRAD,
304 const output_env_t oenv)
306 char name[256], titlestr[256], ystr[256];
313 strcpy(ystr, "Angle (rad)");
317 strcpy(ystr, "Angle (degrees)");
322 for (i = 0; (i < nlist); i++)
324 /* grs debug printf("OK i %d j %d\n", i, j) ; */
327 copy_dih_data(dih[j], data, nf, bRAD);
328 print_one(oenv, "phi", dlist[i].name, "\\xf\\f{}", ystr, nf, time, data);
332 for (i = 0; (i < nlist); i++)
336 copy_dih_data(dih[j], data, nf, bRAD);
337 print_one(oenv, "psi", dlist[i].name, "\\xy\\f{}", ystr, nf, time, data);
341 for (i = 0; (i < nlist); i++)
343 if (has_dihedral(edOmega, &(dlist[i])))
347 copy_dih_data(dih[j], data, nf, bRAD);
348 print_one(oenv, "omega", dlist[i].name, "\\xw\\f{}", ystr, nf, time, data);
354 for (Xi = 0; (Xi < maxchi); Xi++)
356 for (i = 0; (i < nlist); i++)
358 if (dlist[i].atm.Cn[Xi+3] != -1)
362 sprintf(name, "chi%d", Xi+1);
363 sprintf(titlestr, "\\xc\\f{}\\s%d\\N", Xi+1);
364 copy_dih_data(dih[j], data, nf, bRAD);
365 print_one(oenv, name, dlist[i].name, titlestr, ystr, nf, time, data);
371 fprintf(stderr, "\n");
374 static void reset_one(real dih[], int nf, real phase)
378 for (j = 0; (j < nf); j++)
381 while (dih[j] < -M_PI)
385 while (dih[j] >= M_PI)
392 static int reset_em_all(int nlist, t_dlist dlist[], int nf,
393 real **dih, int maxchi)
400 for (i = 0; (i < nlist); i++)
402 if (dlist[i].atm.minC == -1)
404 reset_one(dih[j++], nf, M_PI);
408 reset_one(dih[j++], nf, 0);
412 for (i = 0; (i < nlist-1); i++)
414 reset_one(dih[j++], nf, 0);
416 /* last Psi is faked from O */
417 reset_one(dih[j++], nf, M_PI);
420 for (i = 0; (i < nlist); i++)
422 if (has_dihedral(edOmega, &dlist[i]))
424 reset_one(dih[j++], nf, 0);
427 /* Chi 1 thru maxchi */
428 for (Xi = 0; (Xi < maxchi); Xi++)
430 for (i = 0; (i < nlist); i++)
432 if (dlist[i].atm.Cn[Xi+3] != -1)
434 reset_one(dih[j], nf, 0);
439 fprintf(stderr, "j after resetting (nr. active dihedrals) = %d\n", j);
443 static void histogramming(FILE *log, int nbin, gmx_residuetype_t rt,
444 int nf, int maxchi, real **dih,
445 int nlist, t_dlist dlist[],
447 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bOmega, gmx_bool bChi,
448 gmx_bool bNormalize, gmx_bool bSSHisto, const char *ssdump,
449 real bfac_max, t_atoms *atoms,
450 gmx_bool bDo_jc, const char *fn,
451 const output_env_t oenv)
453 /* also gets 3J couplings and order parameters S2 */
454 t_karplus kkkphi[] = {
455 { "J_NHa1", 6.51, -1.76, 1.6, -M_PI/3, 0.0, 0.0 },
456 { "J_NHa2", 6.51, -1.76, 1.6, M_PI/3, 0.0, 0.0 },
457 { "J_HaC'", 4.0, 1.1, 0.1, 0.0, 0.0, 0.0 },
458 { "J_NHCb", 4.7, -1.5, -0.2, M_PI/3, 0.0, 0.0 },
459 { "J_Ci-1Hai", 4.5, -1.3, -1.2, 2*M_PI/3, 0.0, 0.0 }
461 t_karplus kkkpsi[] = {
462 { "J_HaN", -0.88, -0.61, -0.27, M_PI/3, 0.0, 0.0 }
464 t_karplus kkkchi1[] = {
465 { "JHaHb2", 9.5, -1.6, 1.8, -M_PI/3, 0, 0.0 },
466 { "JHaHb3", 9.5, -1.6, 1.8, 0, 0, 0.0 }
468 #define NKKKPHI asize(kkkphi)
469 #define NKKKPSI asize(kkkpsi)
470 #define NKKKCHI asize(kkkchi1)
471 #define NJC (NKKKPHI+NKKKPSI+NKKKCHI)
473 FILE *fp, *ssfp[3] = {NULL, NULL, NULL};
474 const char *sss[3] = { "sheet", "helix", "coil" };
478 int ****his_aa_ss = NULL;
479 int ***his_aa, **his_aa1, *histmp;
480 int i, j, k, m, n, nn, Dih, nres, hindex, angle;
481 gmx_bool bBfac, bOccup;
482 char hisfile[256], hhisfile[256], sshisfile[256], title[256], *ss_str = NULL;
484 const char *residue_name;
487 rt_size = gmx_residuetype_get_size(rt);
490 fp = gmx_ffopen(ssdump, "r");
491 if (1 != fscanf(fp, "%d", &nres))
493 gmx_fatal(FARGS, "Error reading from file %s", ssdump);
496 snew(ss_str, nres+1);
497 if (1 != fscanf(fp, "%s", ss_str))
499 gmx_fatal(FARGS, "Error reading from file %s", ssdump);
503 /* Four dimensional array... Very cool */
505 for (i = 0; (i < 3); i++)
507 snew(his_aa_ss[i], rt_size+1);
508 for (j = 0; (j <= rt_size); j++)
510 snew(his_aa_ss[i][j], edMax);
511 for (Dih = 0; (Dih < edMax); Dih++)
513 snew(his_aa_ss[i][j][Dih], nbin+1);
519 for (Dih = 0; (Dih < edMax); Dih++)
521 snew(his_aa[Dih], rt_size+1);
522 for (i = 0; (i <= rt_size); i++)
524 snew(his_aa[Dih][i], nbin+1);
531 for (i = 0; (i < nlist); i++)
539 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
541 for (i = 0; (i < nlist); i++)
543 if (((Dih < edOmega) ) ||
544 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
545 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)))
547 make_histo(log, nf, dih[j], nbin, histmp, -M_PI, M_PI);
551 /* Assume there is only one structure, the first.
552 * Compute index in histogram.
554 /* Check the atoms to see whether their B-factors are low enough
555 * Check atoms to see their occupancy is 1.
557 bBfac = bOccup = TRUE;
558 for (nn = 0; (nn < 4); nn++, n++)
560 bBfac = bBfac && (atoms->pdbinfo[index[n]].bfac <= bfac_max);
561 bOccup = bOccup && (atoms->pdbinfo[index[n]].occup == 1);
563 if (bOccup && ((bfac_max <= 0) || ((bfac_max > 0) && bBfac)))
565 hindex = ((dih[j][0]+M_PI)*nbin)/(2*M_PI);
566 range_check(hindex, 0, nbin);
568 /* Assign dihedral to either of the structure determined
571 switch (ss_str[dlist[i].resnr])
574 his_aa_ss[0][dlist[i].index][Dih][hindex]++;
577 his_aa_ss[1][dlist[i].index][Dih][hindex]++;
580 his_aa_ss[2][dlist[i].index][Dih][hindex]++;
586 fprintf(debug, "Res. %d has imcomplete occupancy or bfacs > %g\n",
587 dlist[i].resnr, bfac_max);
598 calc_distribution_props(nbin, histmp, -M_PI, NKKKPHI, kkkphi, &S2);
600 for (m = 0; (m < NKKKPHI); m++)
602 Jc[i][m] = kkkphi[m].Jc;
603 Jcsig[i][m] = kkkphi[m].Jcsig;
607 calc_distribution_props(nbin, histmp, -M_PI, NKKKPSI, kkkpsi, &S2);
609 for (m = 0; (m < NKKKPSI); m++)
611 Jc[i][NKKKPHI+m] = kkkpsi[m].Jc;
612 Jcsig[i][NKKKPHI+m] = kkkpsi[m].Jcsig;
616 calc_distribution_props(nbin, histmp, -M_PI, NKKKCHI, kkkchi1, &S2);
617 for (m = 0; (m < NKKKCHI); m++)
619 Jc[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jc;
620 Jcsig[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jcsig;
623 default: /* covers edOmega and higher Chis than Chi1 */
624 calc_distribution_props(nbin, histmp, -M_PI, 0, NULL, &S2);
627 dlist[i].S2[Dih] = S2;
629 /* Sum distribution per amino acid type as well */
630 for (k = 0; (k < nbin); k++)
632 his_aa[Dih][dlist[i].index][k] += histmp[k];
637 else /* dihed not defined */
639 dlist[i].S2[Dih] = 0.0;
645 /* Print out Jcouplings */
646 fprintf(log, "\n *** J-Couplings from simulation (plus std. dev.) ***\n\n");
647 fprintf(log, "Residue ");
648 for (i = 0; (i < NKKKPHI); i++)
650 fprintf(log, "%7s SD", kkkphi[i].name);
652 for (i = 0; (i < NKKKPSI); i++)
654 fprintf(log, "%7s SD", kkkpsi[i].name);
656 for (i = 0; (i < NKKKCHI); i++)
658 fprintf(log, "%7s SD", kkkchi1[i].name);
661 for (i = 0; (i < NJC+1); i++)
663 fprintf(log, "------------");
666 for (i = 0; (i < nlist); i++)
668 fprintf(log, "%-10s", dlist[i].name);
669 for (j = 0; (j < NJC); j++)
671 fprintf(log, " %5.2f %4.2f", Jc[i][j], Jcsig[i][j]);
677 /* and to -jc file... */
680 fp = xvgropen(fn, "\\S3\\NJ-Couplings from Karplus Equation", "Residue",
683 for (i = 0; (i < NKKKPHI); i++)
685 leg[i] = strdup(kkkphi[i].name);
687 for (i = 0; (i < NKKKPSI); i++)
689 leg[i+NKKKPHI] = strdup(kkkpsi[i].name);
691 for (i = 0; (i < NKKKCHI); i++)
693 leg[i+NKKKPHI+NKKKPSI] = strdup(kkkchi1[i].name);
695 xvgr_legend(fp, NJC, (const char**)leg, oenv);
696 fprintf(fp, "%5s ", "#Res.");
697 for (i = 0; (i < NJC); i++)
699 fprintf(fp, "%10s ", leg[i]);
702 for (i = 0; (i < nlist); i++)
704 fprintf(fp, "%5d ", dlist[i].resnr);
705 for (j = 0; (j < NJC); j++)
707 fprintf(fp, " %8.3f", Jc[i][j]);
712 for (i = 0; (i < NJC); i++)
717 /* finished -jc stuff */
719 snew(normhisto, nbin);
720 for (i = 0; (i < rt_size); i++)
722 for (Dih = 0; (Dih < edMax); Dih++)
724 /* First check whether something is in there */
725 for (j = 0; (j < nbin); j++)
727 if (his_aa[Dih][i][j] != 0)
733 ((bPhi && (Dih == edPhi)) ||
734 (bPsi && (Dih == edPsi)) ||
735 (bOmega && (Dih == edOmega)) ||
736 (bChi && (Dih >= edChi1))))
740 normalize_histo(nbin, his_aa[Dih][i], (360.0/nbin), normhisto);
743 residue_name = gmx_residuetype_get_name(rt, i);
747 sprintf(hisfile, "histo-phi%s", residue_name);
748 sprintf(title, "\\xf\\f{} Distribution for %s", residue_name);
751 sprintf(hisfile, "histo-psi%s", residue_name);
752 sprintf(title, "\\xy\\f{} Distribution for %s", residue_name);
755 sprintf(hisfile, "histo-omega%s", residue_name);
756 sprintf(title, "\\xw\\f{} Distribution for %s", residue_name);
759 sprintf(hisfile, "histo-chi%d%s", Dih-NONCHI+1, residue_name);
760 sprintf(title, "\\xc\\f{}\\s%d\\N Distribution for %s",
761 Dih-NONCHI+1, residue_name);
763 strcpy(hhisfile, hisfile);
764 strcat(hhisfile, ".xvg");
765 fp = xvgropen(hhisfile, title, "Degrees", "", oenv);
766 if (output_env_get_print_xvgr_codes(oenv))
768 fprintf(fp, "@ with g0\n");
770 xvgr_world(fp, -180, 0, 180, 0.1, oenv);
771 if (output_env_get_print_xvgr_codes(oenv))
773 fprintf(fp, "# this effort to set graph size fails unless you run with -autoscale none or -autoscale y flags\n");
774 fprintf(fp, "@ xaxis tick on\n");
775 fprintf(fp, "@ xaxis tick major 90\n");
776 fprintf(fp, "@ xaxis tick minor 30\n");
777 fprintf(fp, "@ xaxis ticklabel prec 0\n");
778 fprintf(fp, "@ yaxis tick off\n");
779 fprintf(fp, "@ yaxis ticklabel off\n");
780 fprintf(fp, "@ type xy\n");
784 for (k = 0; (k < 3); k++)
786 sprintf(sshisfile, "%s-%s.xvg", hisfile, sss[k]);
787 ssfp[k] = gmx_ffopen(sshisfile, "w");
790 for (j = 0; (j < nbin); j++)
792 angle = -180 + (360/nbin)*j;
795 fprintf(fp, "%5d %10g\n", angle, normhisto[j]);
799 fprintf(fp, "%5d %10d\n", angle, his_aa[Dih][i][j]);
803 for (k = 0; (k < 3); k++)
805 fprintf(ssfp[k], "%5d %10d\n", angle,
806 his_aa_ss[k][i][Dih][j]);
810 fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
814 for (k = 0; (k < 3); k++)
816 fprintf(ssfp[k], "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
817 gmx_ffclose(ssfp[k]);
827 /* Four dimensional array... Very cool */
828 for (i = 0; (i < 3); i++)
830 for (j = 0; (j <= rt_size); j++)
832 for (Dih = 0; (Dih < edMax); Dih++)
834 sfree(his_aa_ss[i][j][Dih]);
836 sfree(his_aa_ss[i][j]);
845 static FILE *rama_file(const char *fn, const char *title, const char *xaxis,
846 const char *yaxis, const output_env_t oenv)
850 fp = xvgropen(fn, title, xaxis, yaxis, oenv);
851 if (output_env_get_print_xvgr_codes(oenv))
853 fprintf(fp, "@ with g0\n");
855 xvgr_world(fp, -180, -180, 180, 180, oenv);
856 if (output_env_get_print_xvgr_codes(oenv))
858 fprintf(fp, "@ xaxis tick on\n");
859 fprintf(fp, "@ xaxis tick major 90\n");
860 fprintf(fp, "@ xaxis tick minor 30\n");
861 fprintf(fp, "@ xaxis ticklabel prec 0\n");
862 fprintf(fp, "@ yaxis tick on\n");
863 fprintf(fp, "@ yaxis tick major 90\n");
864 fprintf(fp, "@ yaxis tick minor 30\n");
865 fprintf(fp, "@ yaxis ticklabel prec 0\n");
866 fprintf(fp, "@ s0 type xy\n");
867 fprintf(fp, "@ s0 symbol 2\n");
868 fprintf(fp, "@ s0 symbol size 0.410000\n");
869 fprintf(fp, "@ s0 symbol fill 1\n");
870 fprintf(fp, "@ s0 symbol color 1\n");
871 fprintf(fp, "@ s0 symbol linewidth 1\n");
872 fprintf(fp, "@ s0 symbol linestyle 1\n");
873 fprintf(fp, "@ s0 symbol center false\n");
874 fprintf(fp, "@ s0 symbol char 0\n");
875 fprintf(fp, "@ s0 skip 0\n");
876 fprintf(fp, "@ s0 linestyle 0\n");
877 fprintf(fp, "@ s0 linewidth 1\n");
878 fprintf(fp, "@ type xy\n");
883 static void do_rama(int nf, int nlist, t_dlist dlist[], real **dih,
884 gmx_bool bViol, gmx_bool bRamOmega, const output_env_t oenv)
886 FILE *fp, *gp = NULL;
889 int i, j, k, Xi1, Xi2, Phi, Psi, Om = 0, nlevels;
891 real **mat = NULL, phi, psi, omega, axis[NMAT], lo, hi;
892 t_rgb rlo = { 1.0, 0.0, 0.0 };
893 t_rgb rmid = { 1.0, 1.0, 1.0 };
894 t_rgb rhi = { 0.0, 0.0, 1.0 };
896 for (i = 0; (i < nlist); i++)
898 if ((has_dihedral(edPhi, &(dlist[i]))) &&
899 (has_dihedral(edPsi, &(dlist[i]))))
901 sprintf(fn, "ramaPhiPsi%s.xvg", dlist[i].name);
902 fp = rama_file(fn, "Ramachandran Plot",
903 "\\8f\\4 (deg)", "\\8y\\4 (deg)", oenv);
904 bOm = bRamOmega && has_dihedral(edOmega, &(dlist[i]));
907 Om = dlist[i].j0[edOmega];
909 for (j = 0; (j < NMAT); j++)
912 axis[j] = -180+(360*j)/NMAT;
917 sprintf(fn, "violPhiPsi%s.xvg", dlist[i].name);
918 gp = gmx_ffopen(fn, "w");
920 Phi = dlist[i].j0[edPhi];
921 Psi = dlist[i].j0[edPsi];
922 for (j = 0; (j < nf); j++)
924 phi = RAD2DEG*dih[Phi][j];
925 psi = RAD2DEG*dih[Psi][j];
926 fprintf(fp, "%10g %10g\n", phi, psi);
929 fprintf(gp, "%d\n", !bAllowed(dih[Phi][j], RAD2DEG*dih[Psi][j]));
933 omega = RAD2DEG*dih[Om][j];
934 mat[(int)((phi*NMAT)/360)+NMAT/2][(int)((psi*NMAT)/360)+NMAT/2]
945 sprintf(fn, "ramomega%s.xpm", dlist[i].name);
946 fp = gmx_ffopen(fn, "w");
948 for (j = 0; (j < NMAT); j++)
950 for (k = 0; (k < NMAT); k++)
953 lo = min(mat[j][k], lo);
954 hi = max(mat[j][k], hi);
958 if (fabs(lo) > fabs(hi))
967 for (j = 0; (j < NMAT); j++)
969 for (k = 0; (k < NMAT); k++)
977 write_xpm3(fp, 0, "Omega/Ramachandran Plot", "Deg", "Phi", "Psi",
978 NMAT, NMAT, axis, axis, mat, lo, 180.0, hi, rlo, rmid, rhi, &nlevels);
980 for (j = 0; (j < NMAT); j++)
987 if ((has_dihedral(edChi1, &(dlist[i]))) &&
988 (has_dihedral(edChi2, &(dlist[i]))))
990 sprintf(fn, "ramaX1X2%s.xvg", dlist[i].name);
991 fp = rama_file(fn, "\\8c\\4\\s1\\N-\\8c\\4\\s2\\N Ramachandran Plot",
992 "\\8c\\4\\s1\\N (deg)", "\\8c\\4\\s2\\N (deg)", oenv);
993 Xi1 = dlist[i].j0[edChi1];
994 Xi2 = dlist[i].j0[edChi2];
995 for (j = 0; (j < nf); j++)
997 fprintf(fp, "%10g %10g\n", RAD2DEG*dih[Xi1][j], RAD2DEG*dih[Xi2][j]);
1003 fprintf(stderr, "No chi1 & chi2 angle for %s\n", dlist[i].name);
1009 static void print_transitions(const char *fn, int maxchi, int nlist,
1010 t_dlist dlist[], real dt,
1011 const output_env_t oenv)
1013 /* based on order_params below */
1018 /* must correspond with enum in pp2shift.h:38 */
1020 #define NLEG asize(leg)
1022 leg[0] = strdup("Phi");
1023 leg[1] = strdup("Psi");
1024 leg[2] = strdup("Omega");
1025 leg[3] = strdup("Chi1");
1026 leg[4] = strdup("Chi2");
1027 leg[5] = strdup("Chi3");
1028 leg[6] = strdup("Chi4");
1029 leg[7] = strdup("Chi5");
1030 leg[8] = strdup("Chi6");
1032 /* Print order parameters */
1033 fp = xvgropen(fn, "Dihedral Rotamer Transitions", "Residue", "Transitions/ns",
1035 xvgr_legend(fp, NONCHI+maxchi, (const char**)leg, oenv);
1037 for (Dih = 0; (Dih < edMax); Dih++)
1042 fprintf(fp, "%5s ", "#Res.");
1043 fprintf(fp, "%10s %10s %10s ", leg[edPhi], leg[edPsi], leg[edOmega]);
1044 for (Xi = 0; Xi < maxchi; Xi++)
1046 fprintf(fp, "%10s ", leg[NONCHI+Xi]);
1050 for (i = 0; (i < nlist); i++)
1052 fprintf(fp, "%5d ", dlist[i].resnr);
1053 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1055 fprintf(fp, "%10.3f ", dlist[i].ntr[Dih]/dt);
1057 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
1063 static void order_params(FILE *log,
1064 const char *fn, int maxchi, int nlist, t_dlist dlist[],
1065 const char *pdbfn, real bfac_init,
1066 t_atoms *atoms, rvec x[], int ePBC, matrix box,
1067 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, const output_env_t oenv)
1074 /* except for S2Min/Max, must correspond with enum in pp2shift.h:38 */
1075 const char *const_leg[2+edMax] = {
1076 "S2Min", "S2Max", "Phi", "Psi", "Omega",
1077 "Chi1", "Chi2", "Chi3", "Chi4", "Chi5",
1080 #define NLEG asize(leg)
1084 for (i = 0; i < NLEG; i++)
1086 leg[i] = strdup(const_leg[i]);
1089 /* Print order parameters */
1090 fp = xvgropen(fn, "Dihedral Order Parameters", "Residue", "S2", oenv);
1091 xvgr_legend(fp, 2+NONCHI+maxchi, const_leg, oenv);
1093 for (Dih = 0; (Dih < edMax); Dih++)
1098 fprintf(fp, "%5s ", "#Res.");
1099 fprintf(fp, "%10s %10s ", leg[0], leg[1]);
1100 fprintf(fp, "%10s %10s %10s ", leg[2+edPhi], leg[2+edPsi], leg[2+edOmega]);
1101 for (Xi = 0; Xi < maxchi; Xi++)
1103 fprintf(fp, "%10s ", leg[2+NONCHI+Xi]);
1107 for (i = 0; (i < nlist); i++)
1111 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1113 if (dlist[i].S2[Dih] != 0)
1115 if (dlist[i].S2[Dih] > S2Max)
1117 S2Max = dlist[i].S2[Dih];
1119 if (dlist[i].S2[Dih] < S2Min)
1121 S2Min = dlist[i].S2[Dih];
1124 if (dlist[i].S2[Dih] > 0.8)
1129 fprintf(fp, "%5d ", dlist[i].resnr);
1130 fprintf(fp, "%10.3f %10.3f ", S2Min, S2Max);
1131 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1133 fprintf(fp, "%10.3f ", dlist[i].S2[Dih]);
1136 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
1144 if (NULL == atoms->pdbinfo)
1146 snew(atoms->pdbinfo, atoms->nr);
1148 for (i = 0; (i < atoms->nr); i++)
1150 atoms->pdbinfo[i].bfac = bfac_init;
1153 for (i = 0; (i < nlist); i++)
1155 atoms->pdbinfo[dlist[i].atm.N].bfac = -dlist[i].S2[0]; /* Phi */
1156 atoms->pdbinfo[dlist[i].atm.H].bfac = -dlist[i].S2[0]; /* Phi */
1157 atoms->pdbinfo[dlist[i].atm.C].bfac = -dlist[i].S2[1]; /* Psi */
1158 atoms->pdbinfo[dlist[i].atm.O].bfac = -dlist[i].S2[1]; /* Psi */
1159 for (Xi = 0; (Xi < maxchi); Xi++) /* Chi's */
1161 if (dlist[i].atm.Cn[Xi+3] != -1)
1163 atoms->pdbinfo[dlist[i].atm.Cn[Xi+1]].bfac = -dlist[i].S2[NONCHI+Xi];
1168 fp = gmx_ffopen(pdbfn, "w");
1169 fprintf(fp, "REMARK generated by g_chi\n");
1170 fprintf(fp, "REMARK "
1171 "B-factor field contains negative of dihedral order parameters\n");
1172 write_pdbfile(fp, NULL, atoms, x, ePBC, box, ' ', 0, NULL, TRUE);
1173 x0 = y0 = z0 = 1000.0;
1174 for (i = 0; (i < atoms->nr); i++)
1176 x0 = min(x0, x[i][XX]);
1177 y0 = min(y0, x[i][YY]);
1178 z0 = min(z0, x[i][ZZ]);
1180 x0 *= 10.0; /* nm -> angstrom */
1181 y0 *= 10.0; /* nm -> angstrom */
1182 z0 *= 10.0; /* nm -> angstrom */
1183 for (i = 0; (i < 10); i++)
1185 gmx_fprintf_pdb_atomline(fp, epdbATOM, atoms->nr+1+i, "CA", ' ', "LEG", ' ', atoms->nres+1, ' ',
1186 x0, y0, z0+(1.2*i), 0.0, -0.1*i, "");
1191 fprintf(log, "Dihedrals with S2 > 0.8\n");
1192 fprintf(log, "Dihedral: ");
1195 fprintf(log, " Phi ");
1199 fprintf(log, " Psi ");
1203 for (Xi = 0; (Xi < maxchi); Xi++)
1205 fprintf(log, " %s ", leg[2+NONCHI+Xi]);
1208 fprintf(log, "\nNumber: ");
1211 fprintf(log, "%4d ", nh[0]);
1215 fprintf(log, "%4d ", nh[1]);
1219 for (Xi = 0; (Xi < maxchi); Xi++)
1221 fprintf(log, "%4d ", nh[NONCHI+Xi]);
1226 for (i = 0; i < NLEG; i++)
1233 int gmx_chi(int argc, char *argv[])
1235 const char *desc[] = {
1236 "[THISMODULE] computes [GRK]phi[grk], [GRK]psi[grk], [GRK]omega[grk],",
1237 "and [GRK]chi[grk] dihedrals for all your",
1238 "amino acid backbone and sidechains.",
1239 "It can compute dihedral angle as a function of time, and as",
1240 "histogram distributions.",
1241 "The distributions [TT](histo-(dihedral)(RESIDUE).xvg[tt]) are cumulative over all residues of each type.[PAR]",
1242 "If option [TT]-corr[tt] is given, the program will",
1243 "calculate dihedral autocorrelation functions. The function used",
1244 "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",
1245 "rather than angles themselves, resolves the problem of periodicity.",
1246 "(Van der Spoel & Berendsen (1997), Biophys. J. 72, 2032-2041).",
1247 "Separate files for each dihedral of each residue",
1248 "[TT](corr(dihedral)(RESIDUE)(nresnr).xvg[tt]) are output, as well as a",
1249 "file containing the information for all residues (argument of [TT]-corr[tt]).[PAR]",
1250 "With option [TT]-all[tt], the angles themselves as a function of time for",
1251 "each residue are printed to separate files [TT](dihedral)(RESIDUE)(nresnr).xvg[tt].",
1252 "These can be in radians or degrees.[PAR]",
1253 "A log file (argument [TT]-g[tt]) is also written. This contains [BR]",
1254 "(a) information about the number of residues of each type.[BR]",
1255 "(b) The NMR ^3J coupling constants from the Karplus equation.[BR]",
1256 "(c) a table for each residue of the number of transitions between ",
1257 "rotamers per nanosecond, and the order parameter S^2 of each dihedral.[BR]",
1258 "(d) a table for each residue of the rotamer occupancy.[PAR]",
1259 "All rotamers are taken as 3-fold, except for [GRK]omega[grk] and [GRK]chi[grk] dihedrals",
1260 "to planar groups (i.e. [GRK]chi[grk][SUB]2[sub] of aromatics, Asp and Asn; [GRK]chi[grk][SUB]3[sub] of Glu",
1261 "and Gln; and [GRK]chi[grk][SUB]4[sub] of Arg), which are 2-fold. \"rotamer 0\" means ",
1262 "that the dihedral was not in the core region of each rotamer. ",
1263 "The width of the core region can be set with [TT]-core_rotamer[tt][PAR]",
1265 "The S^2 order parameters are also output to an [TT].xvg[tt] file",
1266 "(argument [TT]-o[tt] ) and optionally as a [TT].pdb[tt] file with",
1267 "the S^2 values as B-factor (argument [TT]-p[tt]). ",
1268 "The total number of rotamer transitions per timestep",
1269 "(argument [TT]-ot[tt]), the number of transitions per rotamer",
1270 "(argument [TT]-rt[tt]), and the ^3J couplings (argument [TT]-jc[tt]), ",
1271 "can also be written to [TT].xvg[tt] files. Note that the analysis",
1272 "of rotamer transitions assumes that the supplied trajectory frames",
1273 "are equally spaced in time.[PAR]",
1275 "If [TT]-chi_prod[tt] is set (and [TT]-maxchi[tt] > 0), cumulative rotamers, e.g.",
1276 "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 ",
1277 "dihedrals and [TT]-maxchi[tt] >= 3)",
1278 "are calculated. As before, if any dihedral is not in the core region,",
1279 "the rotamer is taken to be 0. The occupancies of these cumulative ",
1280 "rotamers (starting with rotamer 0) are written to the file",
1281 "that is the argument of [TT]-cp[tt], and if the [TT]-all[tt] flag",
1282 "is given, the rotamers as functions of time",
1283 "are written to [TT]chiproduct(RESIDUE)(nresnr).xvg[tt] ",
1284 "and their occupancies to [TT]histo-chiproduct(RESIDUE)(nresnr).xvg[tt].[PAR]",
1286 "The option [TT]-r[tt] generates a contour plot of the average [GRK]omega[grk] angle",
1287 "as a function of the [GRK]phi[grk] and [GRK]psi[grk] angles, that is, in a Ramachandran plot",
1288 "the average [GRK]omega[grk] angle is plotted using color coding.",
1292 const char *bugs[] = {
1293 "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.",
1294 "[GRK]phi[grk] and [GRK]psi[grk] dihedrals are calculated in a "
1295 "non-standard way, using H-N-CA-C for [GRK]phi[grk] instead of "
1296 "C(-)-N-CA-C, and N-CA-C-O for [GRK]psi[grk] instead of N-CA-C-N(+). "
1297 "This causes (usually small) discrepancies with the output of other "
1298 "tools like [gmx-rama].",
1299 "[TT]-r0[tt] option does not work properly",
1300 "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"
1304 static int r0 = 1, ndeg = 1, maxchi = 2;
1305 static gmx_bool bAll = FALSE;
1306 static gmx_bool bPhi = FALSE, bPsi = FALSE, bOmega = FALSE;
1307 static real bfac_init = -1.0, bfac_max = 0;
1308 static const char *maxchistr[] = { NULL, "0", "1", "2", "3", "4", "5", "6", NULL };
1309 static gmx_bool bRama = FALSE, bShift = FALSE, bViol = FALSE, bRamOmega = FALSE;
1310 static gmx_bool bNormHisto = TRUE, bChiProduct = FALSE, bHChi = FALSE, bRAD = FALSE, bPBC = TRUE;
1311 static real core_frac = 0.5;
1313 { "-r0", FALSE, etINT, {&r0},
1314 "starting residue" },
1315 { "-phi", FALSE, etBOOL, {&bPhi},
1316 "Output for [GRK]phi[grk] dihedral angles" },
1317 { "-psi", FALSE, etBOOL, {&bPsi},
1318 "Output for [GRK]psi[grk] dihedral angles" },
1319 { "-omega", FALSE, etBOOL, {&bOmega},
1320 "Output for [GRK]omega[grk] dihedrals (peptide bonds)" },
1321 { "-rama", FALSE, etBOOL, {&bRama},
1322 "Generate [GRK]phi[grk]/[GRK]psi[grk] and [GRK]chi[grk][SUB]1[sub]/[GRK]chi[grk][SUB]2[sub] Ramachandran plots" },
1323 { "-viol", FALSE, etBOOL, {&bViol},
1324 "Write a file that gives 0 or 1 for violated Ramachandran angles" },
1325 { "-periodic", FALSE, etBOOL, {&bPBC},
1326 "Print dihedral angles modulo 360 degrees" },
1327 { "-all", FALSE, etBOOL, {&bAll},
1328 "Output separate files for every dihedral." },
1329 { "-rad", FALSE, etBOOL, {&bRAD},
1330 "in angle vs time files, use radians rather than degrees."},
1331 { "-shift", FALSE, etBOOL, {&bShift},
1332 "Compute chemical shifts from [GRK]phi[grk]/[GRK]psi[grk] angles" },
1333 { "-binwidth", FALSE, etINT, {&ndeg},
1334 "bin width for histograms (degrees)" },
1335 { "-core_rotamer", FALSE, etREAL, {&core_frac},
1336 "only the central [TT]-core_rotamer[tt]*(360/multiplicity) belongs to each rotamer (the rest is assigned to rotamer 0)" },
1337 { "-maxchi", FALSE, etENUM, {maxchistr},
1338 "calculate first ndih [GRK]chi[grk] dihedrals" },
1339 { "-normhisto", FALSE, etBOOL, {&bNormHisto},
1340 "Normalize histograms" },
1341 { "-ramomega", FALSE, etBOOL, {&bRamOmega},
1342 "compute average omega as a function of [GRK]phi[grk]/[GRK]psi[grk] and plot it in an [TT].xpm[tt] plot" },
1343 { "-bfact", FALSE, etREAL, {&bfac_init},
1344 "B-factor value for [TT].pdb[tt] file for atoms with no calculated dihedral order parameter"},
1345 { "-chi_prod", FALSE, etBOOL, {&bChiProduct},
1346 "compute a single cumulative rotamer for each residue"},
1347 { "-HChi", FALSE, etBOOL, {&bHChi},
1348 "Include dihedrals to sidechain hydrogens"},
1349 { "-bmax", FALSE, etREAL, {&bfac_max},
1350 "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." }
1354 int natoms, nlist, idum, nbin;
1359 char title[256], grpname[256];
1361 gmx_bool bChi, bCorr, bSSHisto;
1362 gmx_bool bDo_rt, bDo_oh, bDo_ot, bDo_jc;
1363 real dt = 0, traj_t_ns;
1365 gmx_residuetype_t rt;
1367 atom_id isize, *index;
1368 int ndih, nactdih, nf;
1369 real **dih, *trans_frac, *aver_angle, *time;
1370 int i, j, **chi_lookup, *multiplicity;
1373 { efSTX, "-s", NULL, ffREAD },
1374 { efTRX, "-f", NULL, ffREAD },
1375 { efXVG, "-o", "order", ffWRITE },
1376 { efPDB, "-p", "order", ffOPTWR },
1377 { efDAT, "-ss", "ssdump", ffOPTRD },
1378 { efXVG, "-jc", "Jcoupling", ffWRITE },
1379 { efXVG, "-corr", "dihcorr", ffOPTWR },
1380 { efLOG, "-g", "chi", ffWRITE },
1381 /* add two more arguments copying from g_angle */
1382 { efXVG, "-ot", "dihtrans", ffOPTWR },
1383 { efXVG, "-oh", "trhisto", ffOPTWR },
1384 { efXVG, "-rt", "restrans", ffOPTWR },
1385 { efXVG, "-cp", "chiprodhisto", ffOPTWR }
1387 #define NFILE asize(fnm)
1392 ppa = add_acf_pargs(&npargs, pa);
1393 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
1394 NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs,
1400 /* Handle result from enumerated type */
1401 sscanf(maxchistr[0], "%d", &maxchi);
1402 bChi = (maxchi > 0);
1404 log = gmx_ffopen(ftp2fn(efLOG, NFILE, fnm), "w");
1413 /* set some options */
1414 bDo_rt = (opt2bSet("-rt", NFILE, fnm));
1415 bDo_oh = (opt2bSet("-oh", NFILE, fnm));
1416 bDo_ot = (opt2bSet("-ot", NFILE, fnm));
1417 bDo_jc = (opt2bSet("-jc", NFILE, fnm));
1418 bCorr = (opt2bSet("-corr", NFILE, fnm));
1421 fprintf(stderr, "Will calculate autocorrelation\n");
1424 if (core_frac > 1.0)
1426 fprintf(stderr, "core_rotamer fraction > 1.0 ; will use 1.0\n");
1429 if (core_frac < 0.0)
1431 fprintf(stderr, "core_rotamer fraction < 0.0 ; will use 0.0\n");
1435 if (maxchi > MAXCHI)
1438 "Will only calculate first %d Chi dihedrals in stead of %d.\n",
1442 bSSHisto = ftp2bSet(efDAT, NFILE, fnm);
1445 /* Find the chi angles using atoms struct and a list of amino acids */
1446 get_stx_coordnum(ftp2fn(efSTX, NFILE, fnm), &natoms);
1447 init_t_atoms(&atoms, natoms, TRUE);
1449 read_stx_conf(ftp2fn(efSTX, NFILE, fnm), title, &atoms, x, NULL, &ePBC, box);
1450 fprintf(log, "Title: %s\n", title);
1452 gmx_residuetype_init(&rt);
1453 dlist = mk_dlist(log, &atoms, &nlist, bPhi, bPsi, bChi, bHChi, maxchi, r0, rt);
1454 fprintf(stderr, "%d residues with dihedrals found\n", nlist);
1458 gmx_fatal(FARGS, "No dihedrals in your structure!\n");
1461 /* Make a linear index for reading all. */
1462 index = make_chi_ind(nlist, dlist, &ndih);
1464 fprintf(stderr, "%d dihedrals found\n", ndih);
1468 /* COMPUTE ALL DIHEDRALS! */
1469 read_ang_dih(ftp2fn(efTRX, NFILE, fnm), FALSE, TRUE, FALSE, bPBC, 1, &idum,
1470 &nf, &time, isize, index, &trans_frac, &aver_angle, dih, oenv);
1472 dt = (time[nf-1]-time[0])/(nf-1); /* might want this for corr or n. transit*/
1477 gmx_fatal(FARGS, "Need at least 2 frames for correlation");
1481 /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi
1482 * pass nactdih instead of ndih to low_ana_dih_trans and get_chi_product_traj
1483 * to prevent accessing off end of arrays when maxchi < 5 or 6. */
1484 nactdih = reset_em_all(nlist, dlist, nf, dih, maxchi);
1488 dump_em_all(nlist, dlist, nf, time, dih, maxchi, bPhi, bPsi, bChi, bOmega, bRAD, oenv);
1491 /* Histogramming & J coupling constants & calc of S2 order params */
1492 histogramming(log, nbin, rt, nf, maxchi, dih, nlist, dlist, index,
1493 bPhi, bPsi, bOmega, bChi,
1494 bNormHisto, bSSHisto, ftp2fn(efDAT, NFILE, fnm), bfac_max, &atoms,
1495 bDo_jc, opt2fn("-jc", NFILE, fnm), oenv);
1499 * added multiplicity */
1501 snew(multiplicity, ndih);
1502 mk_multiplicity_lookup(multiplicity, maxchi, nlist, dlist, ndih);
1504 strcpy(grpname, "All residues, ");
1507 strcat(grpname, "Phi ");
1511 strcat(grpname, "Psi ");
1515 strcat(grpname, "Omega ");
1519 strcat(grpname, "Chi 1-");
1520 sprintf(grpname + strlen(grpname), "%i", maxchi);
1524 low_ana_dih_trans(bDo_ot, opt2fn("-ot", NFILE, fnm),
1525 bDo_oh, opt2fn("-oh", NFILE, fnm), maxchi,
1526 dih, nlist, dlist, nf, nactdih, grpname, multiplicity,
1527 time, FALSE, core_frac, oenv);
1529 /* Order parameters */
1530 order_params(log, opt2fn("-o", NFILE, fnm), maxchi, nlist, dlist,
1531 ftp2fn_null(efPDB, NFILE, fnm), bfac_init,
1532 &atoms, x, ePBC, box, bPhi, bPsi, bChi, oenv);
1534 /* Print ramachandran maps! */
1537 do_rama(nf, nlist, dlist, dih, bViol, bRamOmega, oenv);
1542 do_pp2shifts(log, nf, nlist, dlist, dih);
1545 /* rprint S^2, transitions, and rotamer occupancies to log */
1546 traj_t_ns = 0.001 * (time[nf-1]-time[0]);
1547 pr_dlist(log, nlist, dlist, traj_t_ns, edPrintST, bPhi, bPsi, bChi, bOmega, maxchi);
1548 pr_dlist(log, nlist, dlist, traj_t_ns, edPrintRO, bPhi, bPsi, bChi, bOmega, maxchi);
1550 /* transitions to xvg */
1553 print_transitions(opt2fn("-rt", NFILE, fnm), maxchi, nlist, dlist, traj_t_ns, oenv);
1556 /* chi_product trajectories (ie one "rotamer number" for each residue) */
1557 if (bChiProduct && bChi)
1559 snew(chi_lookup, nlist);
1560 for (i = 0; i < nlist; i++)
1562 snew(chi_lookup[i], maxchi);
1564 mk_chi_lookup(chi_lookup, maxchi, nlist, dlist);
1566 get_chi_product_traj(dih, nf, nactdih,
1567 maxchi, dlist, time, chi_lookup, multiplicity,
1568 FALSE, bNormHisto, core_frac, bAll,
1569 opt2fn("-cp", NFILE, fnm), oenv);
1571 for (i = 0; i < nlist; i++)
1573 sfree(chi_lookup[i]);
1577 /* Correlation comes last because it messes up the angles */
1580 do_dihcorr(opt2fn("-corr", NFILE, fnm), nf, ndih, dih, dt, nlist, dlist, time,
1581 maxchi, bPhi, bPsi, bChi, bOmega, oenv);
1585 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1586 do_view(oenv, opt2fn("-jc", NFILE, fnm), "-nxy");
1589 do_view(oenv, opt2fn("-corr", NFILE, fnm), "-nxy");
1592 gmx_residuetype_destroy(rt);