?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Move residuetype handling to its own file.
[gromacs.git] / src / gromacs / gmxana / gmx_chi.c
blobb0f46c0a0906d0535826a45313430503c7c9a94b
1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
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.
10 *
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.
15 *
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.
20 *
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.
25 *
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.
33 *
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.
36 */
37 #ifdef HAVE_CONFIG_H
38 #include <config.h>
39 #endif
40
41 #include <math.h>
42 #include <stdio.h>
43 #include <string.h>
44
45 #include "gromacs/fileio/confio.h"
46 #include "gromacs/fileio/pdbio.h"
47 #include "gromacs/utility/fatalerror.h"
48 #include "gromacs/utility/futil.h"
49 #include "gstat.h"
50 #include "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 "txtdump.h"
58 #include "typedefs.h"
59 #include "gromacs/math/vec.h"
60 #include "gromacs/fileio/xvgr.h"
61 #include "viewit.h"
62 #include "gromacs/fileio/matio.h"
63 #include "gmx_ana.h"
64
65 static gmx_bool bAllowed(real phi, real psi)
66 {
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"
129 };
130 #define NPP asize(map)
131 int x, y;
132
133 #define INDEX(ppp) ((((int) (360+ppp*RAD2DEG)) % 360)/6)
134 x = INDEX(phi);
135 y = INDEX(psi);
136 #undef INDEX
137 return (gmx_bool) map[x][y];
138 }
139
140 atom_id *make_chi_ind(int nl, t_dlist dl[], int *ndih)
141 {
142 atom_id *id;
143 int i, Xi, n;
144
145 /* There are nl residues with max edMax dihedrals with 4 atoms each */
146 snew(id, nl*edMax*4);
147
148 n = 0;
149 for (i = 0; (i < nl); i++)
150 {
151 /* Phi, fake the first one */
152 dl[i].j0[edPhi] = n/4;
153 if (dl[i].atm.minC >= 0)
154 {
155 id[n++] = dl[i].atm.minC;
156 }
157 else
158 {
159 id[n++] = dl[i].atm.H;
160 }
161 id[n++] = dl[i].atm.N;
162 id[n++] = dl[i].atm.Cn[1];
163 id[n++] = dl[i].atm.C;
164 }
165 for (i = 0; (i < nl); i++)
166 {
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;
172 if (i < (nl-1) )
173 {
174 id[n++] = dl[i+1].atm.N;
175 }
176 else
177 {
178 id[n++] = dl[i].atm.O;
179 }
180 }
181 for (i = 1; (i < nl); i++)
182 {
183 /* Omega */
184 if (has_dihedral(edOmega, &(dl[i])))
185 {
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];
191 }
192 }
193 for (Xi = 0; (Xi < MAXCHI); Xi++)
194 {
195 /* Chi# */
196 for (i = 0; (i < nl); i++)
197 {
198 if (dl[i].atm.Cn[Xi+3] != -1)
199 {
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];
205 }
206 }
207 }
208 *ndih = n/4;
209
210 return id;
211 }
212
213 int bin(real chi, int mult)
214 {
215 mult = 3;
216
217 return (int) (chi*mult/360.0);
218 }
219
220
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)
225 {
226 char name1[256], name2[256];
227 int i, j, Xi;
228
229 do_autocorr(fn, oenv, "Dihedral Autocorrelation Function",
230 nf, ndih, dih, dt, eacCos, FALSE);
231 /* Dump em all */
232 j = 0;
233 for (i = 0; (i < nlist); i++)
234 {
235 if (bPhi)
236 {
237 print_one(oenv, "corrphi", dlist[i].name, "Phi ACF for", "C(t)", nf/2, time,
238 dih[j]);
239 }
240 j++;
241 }
242 for (i = 0; (i < nlist); i++)
243 {
244 if (bPsi)
245 {
246 print_one(oenv, "corrpsi", dlist[i].name, "Psi ACF for", "C(t)", nf/2, time,
247 dih[j]);
248 }
249 j++;
250 }
251 for (i = 0; (i < nlist); i++)
252 {
253 if (has_dihedral(edOmega, &dlist[i]))
254 {
255 if (bOmega)
256 {
257 print_one(oenv, "corromega", dlist[i].name, "Omega ACF for", "C(t)",
258 nf/2, time, dih[j]);
259 }
260 j++;
261 }
262 }
263 for (Xi = 0; (Xi < maxchi); Xi++)
264 {
265 sprintf(name1, "corrchi%d", Xi+1);
266 sprintf(name2, "Chi%d ACF for", Xi+1);
267 for (i = 0; (i < nlist); i++)
268 {
269 if (dlist[i].atm.Cn[Xi+3] != -1)
270 {
271 if (bChi)
272 {
273 print_one(oenv, name1, dlist[i].name, name2, "C(t)", nf/2, time, dih[j]);
274 }
275 j++;
276 }
277 }
278 }
279 fprintf(stderr, "\n");
280 }
281
282 static void copy_dih_data(real in[], real out[], int nf, gmx_bool bLEAVE)
283 {
284 /* if bLEAVE, do nothing to data in copying to out
285 * otherwise multiply by 180/pi to convert rad to deg */
286 int i;
287 real mult;
288 if (bLEAVE)
289 {
290 mult = 1;
291 }
292 else
293 {
294 mult = (180.0/M_PI);
295 }
296 for (i = 0; (i < nf); i++)
297 {
298 out[i] = in[i]*mult;
299 }
300 }
301
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)
306 {
307 char name[256], titlestr[256], ystr[256];
308 real *data;
309 int i, j, Xi;
310
311 snew(data, nf);
312 if (bRAD)
313 {
314 strcpy(ystr, "Angle (rad)");
315 }
316 else
317 {
318 strcpy(ystr, "Angle (degrees)");
319 }
320
321 /* Dump em all */
322 j = 0;
323 for (i = 0; (i < nlist); i++)
324 {
325 /* grs debug printf("OK i %d j %d\n", i, j) ; */
326 if (bPhi)
327 {
328 copy_dih_data(dih[j], data, nf, bRAD);
329 print_one(oenv, "phi", dlist[i].name, "\\xf\\f{}", ystr, nf, time, data);
330 }
331 j++;
332 }
333 for (i = 0; (i < nlist); i++)
334 {
335 if (bPsi)
336 {
337 copy_dih_data(dih[j], data, nf, bRAD);
338 print_one(oenv, "psi", dlist[i].name, "\\xy\\f{}", ystr, nf, time, data);
339 }
340 j++;
341 }
342 for (i = 0; (i < nlist); i++)
343 {
344 if (has_dihedral(edOmega, &(dlist[i])))
345 {
346 if (bOmega)
347 {
348 copy_dih_data(dih[j], data, nf, bRAD);
349 print_one(oenv, "omega", dlist[i].name, "\\xw\\f{}", ystr, nf, time, data);
350 }
351 j++;
352 }
353 }
354
355 for (Xi = 0; (Xi < maxchi); Xi++)
356 {
357 for (i = 0; (i < nlist); i++)
358 {
359 if (dlist[i].atm.Cn[Xi+3] != -1)
360 {
361 if (bChi)
362 {
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);
367 }
368 j++;
369 }
370 }
371 }
372 fprintf(stderr, "\n");
373 }
374
375 static void reset_one(real dih[], int nf, real phase)
376 {
377 int j;
378
379 for (j = 0; (j < nf); j++)
380 {
381 dih[j] += phase;
382 while (dih[j] < -M_PI)
383 {
384 dih[j] += 2*M_PI;
385 }
386 while (dih[j] >= M_PI)
387 {
388 dih[j] -= 2*M_PI;
389 }
390 }
391 }
392
393 static int reset_em_all(int nlist, t_dlist dlist[], int nf,
394 real **dih, int maxchi)
395 {
396 int i, j, Xi;
397
398 /* Reset em all */
399 j = 0;
400 /* Phi */
401 for (i = 0; (i < nlist); i++)
402 {
403 if (dlist[i].atm.minC == -1)
404 {
405 reset_one(dih[j++], nf, M_PI);
406 }
407 else
408 {
409 reset_one(dih[j++], nf, 0);
410 }
411 }
412 /* Psi */
413 for (i = 0; (i < nlist-1); i++)
414 {
415 reset_one(dih[j++], nf, 0);
416 }
417 /* last Psi is faked from O */
418 reset_one(dih[j++], nf, M_PI);
419
420 /* Omega */
421 for (i = 0; (i < nlist); i++)
422 {
423 if (has_dihedral(edOmega, &dlist[i]))
424 {
425 reset_one(dih[j++], nf, 0);
426 }
427 }
428 /* Chi 1 thru maxchi */
429 for (Xi = 0; (Xi < maxchi); Xi++)
430 {
431 for (i = 0; (i < nlist); i++)
432 {
433 if (dlist[i].atm.Cn[Xi+3] != -1)
434 {
435 reset_one(dih[j], nf, 0);
436 j++;
437 }
438 }
439 }
440 fprintf(stderr, "j after resetting (nr. active dihedrals) = %d\n", j);
441 return j;
442 }
443
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[],
447 atom_id index[],
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)
453 {
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 }
461 };
462 t_karplus kkkpsi[] = {
463 { "J_HaN", -0.88, -0.61, -0.27, M_PI/3, 0.0, 0.0 }
464 };
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 }
468 };
469 #define NKKKPHI asize(kkkphi)
470 #define NKKKPSI asize(kkkpsi)
471 #define NKKKCHI asize(kkkchi1)
472 #define NJC (NKKKPHI+NKKKPSI+NKKKCHI)
473
474 FILE *fp, *ssfp[3] = {NULL, NULL, NULL};
475 const char *sss[3] = { "sheet", "helix", "coil" };
476 real S2;
477 real *normhisto;
478 real **Jc, **Jcsig;
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;
484 char **leg;
485 const char *residue_name;
486 int rt_size;
487
488 rt_size = gmx_residuetype_get_size(rt);
489 if (bSSHisto)
490 {
491 fp = gmx_ffopen(ssdump, "r");
492 if (1 != fscanf(fp, "%d", &nres))
493 {
494 gmx_fatal(FARGS, "Error reading from file %s", ssdump);
495 }
496
497 snew(ss_str, nres+1);
498 if (1 != fscanf(fp, "%s", ss_str))
499 {
500 gmx_fatal(FARGS, "Error reading from file %s", ssdump);
501 }
502
503 gmx_ffclose(fp);
504 /* Four dimensional array... Very cool */
505 snew(his_aa_ss, 3);
506 for (i = 0; (i < 3); i++)
507 {
508 snew(his_aa_ss[i], rt_size+1);
509 for (j = 0; (j <= rt_size); j++)
510 {
511 snew(his_aa_ss[i][j], edMax);
512 for (Dih = 0; (Dih < edMax); Dih++)
513 {
514 snew(his_aa_ss[i][j][Dih], nbin+1);
515 }
516 }
517 }
518 }
519 snew(his_aa, edMax);
520 for (Dih = 0; (Dih < edMax); Dih++)
521 {
522 snew(his_aa[Dih], rt_size+1);
523 for (i = 0; (i <= rt_size); i++)
524 {
525 snew(his_aa[Dih][i], nbin+1);
526 }
527 }
528 snew(histmp, nbin);
529
530 snew(Jc, nlist);
531 snew(Jcsig, nlist);
532 for (i = 0; (i < nlist); i++)
533 {
534 snew(Jc[i], NJC);
535 snew(Jcsig[i], NJC);
536 }
537
538 j = 0;
539 n = 0;
540 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
541 {
542 for (i = 0; (i < nlist); i++)
543 {
544 if (((Dih < edOmega) ) ||
545 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
546 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)))
547 {
548 make_histo(log, nf, dih[j], nbin, histmp, -M_PI, M_PI);
549
550 if (bSSHisto)
551 {
552 /* Assume there is only one structure, the first.
553 * Compute index in histogram.
554 */
555 /* Check the atoms to see whether their B-factors are low enough
556 * Check atoms to see their occupancy is 1.
557 */
558 bBfac = bOccup = TRUE;
559 for (nn = 0; (nn < 4); nn++, n++)
560 {
561 bBfac = bBfac && (atoms->pdbinfo[index[n]].bfac <= bfac_max);
562 bOccup = bOccup && (atoms->pdbinfo[index[n]].occup == 1);
563 }
564 if (bOccup && ((bfac_max <= 0) || ((bfac_max > 0) && bBfac)))
565 {
566 hindex = ((dih[j][0]+M_PI)*nbin)/(2*M_PI);
567 range_check(hindex, 0, nbin);
568
569 /* Assign dihedral to either of the structure determined
570 * histograms
571 */
572 switch (ss_str[dlist[i].resnr])
573 {
574 case 'E':
575 his_aa_ss[0][dlist[i].index][Dih][hindex]++;
576 break;
577 case 'H':
578 his_aa_ss[1][dlist[i].index][Dih][hindex]++;
579 break;
580 default:
581 his_aa_ss[2][dlist[i].index][Dih][hindex]++;
582 break;
583 }
584 }
585 else if (debug)
586 {
587 fprintf(debug, "Res. %d has imcomplete occupancy or bfacs > %g\n",
588 dlist[i].resnr, bfac_max);
589 }
590 }
591 else
592 {
593 n += 4;
594 }
595
596 switch (Dih)
597 {
598 case edPhi:
599 calc_distribution_props(nbin, histmp, -M_PI, NKKKPHI, kkkphi, &S2);
600
601 for (m = 0; (m < NKKKPHI); m++)
602 {
603 Jc[i][m] = kkkphi[m].Jc;
604 Jcsig[i][m] = kkkphi[m].Jcsig;
605 }
606 break;
607 case edPsi:
608 calc_distribution_props(nbin, histmp, -M_PI, NKKKPSI, kkkpsi, &S2);
609
610 for (m = 0; (m < NKKKPSI); m++)
611 {
612 Jc[i][NKKKPHI+m] = kkkpsi[m].Jc;
613 Jcsig[i][NKKKPHI+m] = kkkpsi[m].Jcsig;
614 }
615 break;
616 case edChi1:
617 calc_distribution_props(nbin, histmp, -M_PI, NKKKCHI, kkkchi1, &S2);
618 for (m = 0; (m < NKKKCHI); m++)
619 {
620 Jc[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jc;
621 Jcsig[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jcsig;
622 }
623 break;
624 default: /* covers edOmega and higher Chis than Chi1 */
625 calc_distribution_props(nbin, histmp, -M_PI, 0, NULL, &S2);
626 break;
627 }
628 dlist[i].S2[Dih] = S2;
629
630 /* Sum distribution per amino acid type as well */
631 for (k = 0; (k < nbin); k++)
632 {
633 his_aa[Dih][dlist[i].index][k] += histmp[k];
634 histmp[k] = 0;
635 }
636 j++;
637 }
638 else /* dihed not defined */
639 {
640 dlist[i].S2[Dih] = 0.0;
641 }
642 }
643 }
644 sfree(histmp);
645
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++)
650 {
651 fprintf(log, "%7s SD", kkkphi[i].name);
652 }
653 for (i = 0; (i < NKKKPSI); i++)
654 {
655 fprintf(log, "%7s SD", kkkpsi[i].name);
656 }
657 for (i = 0; (i < NKKKCHI); i++)
658 {
659 fprintf(log, "%7s SD", kkkchi1[i].name);
660 }
661 fprintf(log, "\n");
662 for (i = 0; (i < NJC+1); i++)
663 {
664 fprintf(log, "------------");
665 }
666 fprintf(log, "\n");
667 for (i = 0; (i < nlist); i++)
668 {
669 fprintf(log, "%-10s", dlist[i].name);
670 for (j = 0; (j < NJC); j++)
671 {
672 fprintf(log, " %5.2f %4.2f", Jc[i][j], Jcsig[i][j]);
673 }
674 fprintf(log, "\n");
675 }
676 fprintf(log, "\n");
677
678 /* and to -jc file... */
679 if (bDo_jc)
680 {
681 fp = xvgropen(fn, "\\S3\\NJ-Couplings from Karplus Equation", "Residue",
682 "Coupling", oenv);
683 snew(leg, NJC);
684 for (i = 0; (i < NKKKPHI); i++)
685 {
686 leg[i] = strdup(kkkphi[i].name);
687 }
688 for (i = 0; (i < NKKKPSI); i++)
689 {
690 leg[i+NKKKPHI] = strdup(kkkpsi[i].name);
691 }
692 for (i = 0; (i < NKKKCHI); i++)
693 {
694 leg[i+NKKKPHI+NKKKPSI] = strdup(kkkchi1[i].name);
695 }
696 xvgr_legend(fp, NJC, (const char**)leg, oenv);
697 fprintf(fp, "%5s ", "#Res.");
698 for (i = 0; (i < NJC); i++)
699 {
700 fprintf(fp, "%10s ", leg[i]);
701 }
702 fprintf(fp, "\n");
703 for (i = 0; (i < nlist); i++)
704 {
705 fprintf(fp, "%5d ", dlist[i].resnr);
706 for (j = 0; (j < NJC); j++)
707 {
708 fprintf(fp, " %8.3f", Jc[i][j]);
709 }
710 fprintf(fp, "\n");
711 }
712 gmx_ffclose(fp);
713 for (i = 0; (i < NJC); i++)
714 {
715 sfree(leg[i]);
716 }
717 }
718 /* finished -jc stuff */
719
720 snew(normhisto, nbin);
721 for (i = 0; (i < rt_size); i++)
722 {
723 for (Dih = 0; (Dih < edMax); Dih++)
724 {
725 /* First check whether something is in there */
726 for (j = 0; (j < nbin); j++)
727 {
728 if (his_aa[Dih][i][j] != 0)
729 {
730 break;
731 }
732 }
733 if ((j < nbin) &&
734 ((bPhi && (Dih == edPhi)) ||
735 (bPsi && (Dih == edPsi)) ||
736 (bOmega && (Dih == edOmega)) ||
737 (bChi && (Dih >= edChi1))))
738 {
739 if (bNormalize)
740 {
741 normalize_histo(nbin, his_aa[Dih][i], (360.0/nbin), normhisto);
742 }
743
744 residue_name = gmx_residuetype_get_name(rt, i);
745 switch (Dih)
746 {
747 case edPhi:
748 sprintf(hisfile, "histo-phi%s", residue_name);
749 sprintf(title, "\\xf\\f{} Distribution for %s", residue_name);
750 break;
751 case edPsi:
752 sprintf(hisfile, "histo-psi%s", residue_name);
753 sprintf(title, "\\xy\\f{} Distribution for %s", residue_name);
754 break;
755 case edOmega:
756 sprintf(hisfile, "histo-omega%s", residue_name);
757 sprintf(title, "\\xw\\f{} Distribution for %s", residue_name);
758 break;
759 default:
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);
763 }
764 strcpy(hhisfile, hisfile);
765 strcat(hhisfile, ".xvg");
766 fp = xvgropen(hhisfile, title, "Degrees", "", oenv);
767 fprintf(fp, "@ with g0\n");
768 xvgr_world(fp, -180, 0, 180, 0.1, oenv);
769 fprintf(fp, "# this effort to set graph size fails unless you run with -autoscale none or -autoscale y flags\n");
770 fprintf(fp, "@ xaxis tick on\n");
771 fprintf(fp, "@ xaxis tick major 90\n");
772 fprintf(fp, "@ xaxis tick minor 30\n");
773 fprintf(fp, "@ xaxis ticklabel prec 0\n");
774 fprintf(fp, "@ yaxis tick off\n");
775 fprintf(fp, "@ yaxis ticklabel off\n");
776 fprintf(fp, "@ type xy\n");
777 if (bSSHisto)
778 {
779 for (k = 0; (k < 3); k++)
780 {
781 sprintf(sshisfile, "%s-%s.xvg", hisfile, sss[k]);
782 ssfp[k] = gmx_ffopen(sshisfile, "w");
783 }
784 }
785 for (j = 0; (j < nbin); j++)
786 {
787 angle = -180 + (360/nbin)*j;
788 if (bNormalize)
789 {
790 fprintf(fp, "%5d %10g\n", angle, normhisto[j]);
791 }
792 else
793 {
794 fprintf(fp, "%5d %10d\n", angle, his_aa[Dih][i][j]);
795 }
796 if (bSSHisto)
797 {
798 for (k = 0; (k < 3); k++)
799 {
800 fprintf(ssfp[k], "%5d %10d\n", angle,
801 his_aa_ss[k][i][Dih][j]);
802 }
803 }
804 }
805 fprintf(fp, "&\n");
806 gmx_ffclose(fp);
807 if (bSSHisto)
808 {
809 for (k = 0; (k < 3); k++)
810 {
811 fprintf(ssfp[k], "&\n");
812 gmx_ffclose(ssfp[k]);
813 }
814 }
815 }
816 }
817 }
818 sfree(normhisto);
819
820 if (bSSHisto)
821 {
822 /* Four dimensional array... Very cool */
823 for (i = 0; (i < 3); i++)
824 {
825 for (j = 0; (j <= rt_size); j++)
826 {
827 for (Dih = 0; (Dih < edMax); Dih++)
828 {
829 sfree(his_aa_ss[i][j][Dih]);
830 }
831 sfree(his_aa_ss[i][j]);
832 }
833 sfree(his_aa_ss[i]);
834 }
835 sfree(his_aa_ss);
836 sfree(ss_str);
837 }
838 }
839
840 static FILE *rama_file(const char *fn, const char *title, const char *xaxis,
841 const char *yaxis, const output_env_t oenv)
842 {
843 FILE *fp;
844
845 fp = xvgropen(fn, title, xaxis, yaxis, oenv);
846 fprintf(fp, "@ with g0\n");
847 xvgr_world(fp, -180, -180, 180, 180, oenv);
848 fprintf(fp, "@ xaxis tick on\n");
849 fprintf(fp, "@ xaxis tick major 90\n");
850 fprintf(fp, "@ xaxis tick minor 30\n");
851 fprintf(fp, "@ xaxis ticklabel prec 0\n");
852 fprintf(fp, "@ yaxis tick on\n");
853 fprintf(fp, "@ yaxis tick major 90\n");
854 fprintf(fp, "@ yaxis tick minor 30\n");
855 fprintf(fp, "@ yaxis ticklabel prec 0\n");
856 fprintf(fp, "@ s0 type xy\n");
857 fprintf(fp, "@ s0 symbol 2\n");
858 fprintf(fp, "@ s0 symbol size 0.410000\n");
859 fprintf(fp, "@ s0 symbol fill 1\n");
860 fprintf(fp, "@ s0 symbol color 1\n");
861 fprintf(fp, "@ s0 symbol linewidth 1\n");
862 fprintf(fp, "@ s0 symbol linestyle 1\n");
863 fprintf(fp, "@ s0 symbol center false\n");
864 fprintf(fp, "@ s0 symbol char 0\n");
865 fprintf(fp, "@ s0 skip 0\n");
866 fprintf(fp, "@ s0 linestyle 0\n");
867 fprintf(fp, "@ s0 linewidth 1\n");
868 fprintf(fp, "@ type xy\n");
869
870 return fp;
871 }
872
873 static void do_rama(int nf, int nlist, t_dlist dlist[], real **dih,
874 gmx_bool bViol, gmx_bool bRamOmega, const output_env_t oenv)
875 {
876 FILE *fp, *gp = NULL;
877 gmx_bool bOm;
878 char fn[256];
879 int i, j, k, Xi1, Xi2, Phi, Psi, Om = 0, nlevels;
880 #define NMAT 120
881 real **mat = NULL, phi, psi, omega, axis[NMAT], lo, hi;
882 t_rgb rlo = { 1.0, 0.0, 0.0 };
883 t_rgb rmid = { 1.0, 1.0, 1.0 };
884 t_rgb rhi = { 0.0, 0.0, 1.0 };
885
886 for (i = 0; (i < nlist); i++)
887 {
888 if ((has_dihedral(edPhi, &(dlist[i]))) &&
889 (has_dihedral(edPsi, &(dlist[i]))))
890 {
891 sprintf(fn, "ramaPhiPsi%s.xvg", dlist[i].name);
892 fp = rama_file(fn, "Ramachandran Plot",
893 "\\8f\\4 (deg)", "\\8y\\4 (deg)", oenv);
894 bOm = bRamOmega && has_dihedral(edOmega, &(dlist[i]));
895 if (bOm)
896 {
897 Om = dlist[i].j0[edOmega];
898 snew(mat, NMAT);
899 for (j = 0; (j < NMAT); j++)
900 {
901 snew(mat[j], NMAT);
902 axis[j] = -180+(360*j)/NMAT;
903 }
904 }
905 if (bViol)
906 {
907 sprintf(fn, "violPhiPsi%s.xvg", dlist[i].name);
908 gp = gmx_ffopen(fn, "w");
909 }
910 Phi = dlist[i].j0[edPhi];
911 Psi = dlist[i].j0[edPsi];
912 for (j = 0; (j < nf); j++)
913 {
914 phi = RAD2DEG*dih[Phi][j];
915 psi = RAD2DEG*dih[Psi][j];
916 fprintf(fp, "%10g %10g\n", phi, psi);
917 if (bViol)
918 {
919 fprintf(gp, "%d\n", !bAllowed(dih[Phi][j], RAD2DEG*dih[Psi][j]));
920 }
921 if (bOm)
922 {
923 omega = RAD2DEG*dih[Om][j];
924 mat[(int)((phi*NMAT)/360)+NMAT/2][(int)((psi*NMAT)/360)+NMAT/2]
925 += omega;
926 }
927 }
928 if (bViol)
929 {
930 gmx_ffclose(gp);
931 }
932 gmx_ffclose(fp);
933 if (bOm)
934 {
935 sprintf(fn, "ramomega%s.xpm", dlist[i].name);
936 fp = gmx_ffopen(fn, "w");
937 lo = hi = 0;
938 for (j = 0; (j < NMAT); j++)
939 {
940 for (k = 0; (k < NMAT); k++)
941 {
942 mat[j][k] /= nf;
943 lo = min(mat[j][k], lo);
944 hi = max(mat[j][k], hi);
945 }
946 }
947 /* Symmetrise */
948 if (fabs(lo) > fabs(hi))
949 {
950 hi = -lo;
951 }
952 else
953 {
954 lo = -hi;
955 }
956 /* Add 180 */
957 for (j = 0; (j < NMAT); j++)
958 {
959 for (k = 0; (k < NMAT); k++)
960 {
961 mat[j][k] += 180;
962 }
963 }
964 lo += 180;
965 hi += 180;
966 nlevels = 20;
967 write_xpm3(fp, 0, "Omega/Ramachandran Plot", "Deg", "Phi", "Psi",
968 NMAT, NMAT, axis, axis, mat, lo, 180.0, hi, rlo, rmid, rhi, &nlevels);
969 gmx_ffclose(fp);
970 for (j = 0; (j < NMAT); j++)
971 {
972 sfree(mat[j]);
973 }
974 sfree(mat);
975 }
976 }
977 if ((has_dihedral(edChi1, &(dlist[i]))) &&
978 (has_dihedral(edChi2, &(dlist[i]))))
979 {
980 sprintf(fn, "ramaX1X2%s.xvg", dlist[i].name);
981 fp = rama_file(fn, "\\8c\\4\\s1\\N-\\8c\\4\\s2\\N Ramachandran Plot",
982 "\\8c\\4\\s1\\N (deg)", "\\8c\\4\\s2\\N (deg)", oenv);
983 Xi1 = dlist[i].j0[edChi1];
984 Xi2 = dlist[i].j0[edChi2];
985 for (j = 0; (j < nf); j++)
986 {
987 fprintf(fp, "%10g %10g\n", RAD2DEG*dih[Xi1][j], RAD2DEG*dih[Xi2][j]);
988 }
989 gmx_ffclose(fp);
990 }
991 else
992 {
993 fprintf(stderr, "No chi1 & chi2 angle for %s\n", dlist[i].name);
994 }
995 }
996 }
997
998
999 static void print_transitions(const char *fn, int maxchi, int nlist,
1000 t_dlist dlist[], real dt,
1001 const output_env_t oenv)
1002 {
1003 /* based on order_params below */
1004 FILE *fp;
1005 int nh[edMax];
1006 int i, Dih, Xi;
1007
1008 /* must correspond with enum in pp2shift.h:38 */
1009 char *leg[edMax];
1010 #define NLEG asize(leg)
1011
1012 leg[0] = strdup("Phi");
1013 leg[1] = strdup("Psi");
1014 leg[2] = strdup("Omega");
1015 leg[3] = strdup("Chi1");
1016 leg[4] = strdup("Chi2");
1017 leg[5] = strdup("Chi3");
1018 leg[6] = strdup("Chi4");
1019 leg[7] = strdup("Chi5");
1020 leg[8] = strdup("Chi6");
1021
1022 /* Print order parameters */
1023 fp = xvgropen(fn, "Dihedral Rotamer Transitions", "Residue", "Transitions/ns",
1024 oenv);
1025 xvgr_legend(fp, NONCHI+maxchi, (const char**)leg, oenv);
1026
1027 for (Dih = 0; (Dih < edMax); Dih++)
1028 {
1029 nh[Dih] = 0;
1030 }
1031
1032 fprintf(fp, "%5s ", "#Res.");
1033 fprintf(fp, "%10s %10s %10s ", leg[edPhi], leg[edPsi], leg[edOmega]);
1034 for (Xi = 0; Xi < maxchi; Xi++)
1035 {
1036 fprintf(fp, "%10s ", leg[NONCHI+Xi]);
1037 }
1038 fprintf(fp, "\n");
1039
1040 for (i = 0; (i < nlist); i++)
1041 {
1042 fprintf(fp, "%5d ", dlist[i].resnr);
1043 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1044 {
1045 fprintf(fp, "%10.3f ", dlist[i].ntr[Dih]/dt);
1046 }
1047 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
1048 fprintf(fp, "\n");
1049 }
1050 gmx_ffclose(fp);
1051 }
1052
1053 static void order_params(FILE *log,
1054 const char *fn, int maxchi, int nlist, t_dlist dlist[],
1055 const char *pdbfn, real bfac_init,
1056 t_atoms *atoms, rvec x[], int ePBC, matrix box,
1057 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, const output_env_t oenv)
1058 {
1059 FILE *fp;
1060 int nh[edMax];
1061 char buf[STRLEN];
1062 int i, Dih, Xi;
1063 real S2Max, S2Min;
1064
1065 /* except for S2Min/Max, must correspond with enum in pp2shift.h:38 */
1066 const char *const_leg[2+edMax] = {
1067 "S2Min", "S2Max", "Phi", "Psi", "Omega",
1068 "Chi1", "Chi2", "Chi3", "Chi4", "Chi5",
1069 "Chi6"
1070 };
1071 #define NLEG asize(leg)
1072
1073 char *leg[2+edMax];
1074
1075 for (i = 0; i < NLEG; i++)
1076 {
1077 leg[i] = strdup(const_leg[i]);
1078 }
1079
1080 /* Print order parameters */
1081 fp = xvgropen(fn, "Dihedral Order Parameters", "Residue", "S2", oenv);
1082 xvgr_legend(fp, 2+NONCHI+maxchi, const_leg, oenv);
1083
1084 for (Dih = 0; (Dih < edMax); Dih++)
1085 {
1086 nh[Dih] = 0;
1087 }
1088
1089 fprintf(fp, "%5s ", "#Res.");
1090 fprintf(fp, "%10s %10s ", leg[0], leg[1]);
1091 fprintf(fp, "%10s %10s %10s ", leg[2+edPhi], leg[2+edPsi], leg[2+edOmega]);
1092 for (Xi = 0; Xi < maxchi; Xi++)
1093 {
1094 fprintf(fp, "%10s ", leg[2+NONCHI+Xi]);
1095 }
1096 fprintf(fp, "\n");
1097
1098 for (i = 0; (i < nlist); i++)
1099 {
1100 S2Max = -10;
1101 S2Min = 10;
1102 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1103 {
1104 if (dlist[i].S2[Dih] != 0)
1105 {
1106 if (dlist[i].S2[Dih] > S2Max)
1107 {
1108 S2Max = dlist[i].S2[Dih];
1109 }
1110 if (dlist[i].S2[Dih] < S2Min)
1111 {
1112 S2Min = dlist[i].S2[Dih];
1113 }
1114 }
1115 if (dlist[i].S2[Dih] > 0.8)
1116 {
1117 nh[Dih]++;
1118 }
1119 }
1120 fprintf(fp, "%5d ", dlist[i].resnr);
1121 fprintf(fp, "%10.3f %10.3f ", S2Min, S2Max);
1122 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1123 {
1124 fprintf(fp, "%10.3f ", dlist[i].S2[Dih]);
1125 }
1126 fprintf(fp, "\n");
1127 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
1128 }
1129 gmx_ffclose(fp);
1130
1131 if (NULL != pdbfn)
1132 {
1133 real x0, y0, z0;
1134
1135 if (NULL == atoms->pdbinfo)
1136 {
1137 snew(atoms->pdbinfo, atoms->nr);
1138 }
1139 for (i = 0; (i < atoms->nr); i++)
1140 {
1141 atoms->pdbinfo[i].bfac = bfac_init;
1142 }
1143
1144 for (i = 0; (i < nlist); i++)
1145 {
1146 atoms->pdbinfo[dlist[i].atm.N].bfac = -dlist[i].S2[0]; /* Phi */
1147 atoms->pdbinfo[dlist[i].atm.H].bfac = -dlist[i].S2[0]; /* Phi */
1148 atoms->pdbinfo[dlist[i].atm.C].bfac = -dlist[i].S2[1]; /* Psi */
1149 atoms->pdbinfo[dlist[i].atm.O].bfac = -dlist[i].S2[1]; /* Psi */
1150 for (Xi = 0; (Xi < maxchi); Xi++) /* Chi's */
1151 {
1152 if (dlist[i].atm.Cn[Xi+3] != -1)
1153 {
1154 atoms->pdbinfo[dlist[i].atm.Cn[Xi+1]].bfac = -dlist[i].S2[NONCHI+Xi];
1155 }
1156 }
1157 }
1158
1159 fp = gmx_ffopen(pdbfn, "w");
1160 fprintf(fp, "REMARK generated by g_chi\n");
1161 fprintf(fp, "REMARK "
1162 "B-factor field contains negative of dihedral order parameters\n");
1163 write_pdbfile(fp, NULL, atoms, x, ePBC, box, ' ', 0, NULL, TRUE);
1164 x0 = y0 = z0 = 1000.0;
1165 for (i = 0; (i < atoms->nr); i++)
1166 {
1167 x0 = min(x0, x[i][XX]);
1168 y0 = min(y0, x[i][YY]);
1169 z0 = min(z0, x[i][ZZ]);
1170 }
1171 x0 *= 10.0; /* nm -> angstrom */
1172 y0 *= 10.0; /* nm -> angstrom */
1173 z0 *= 10.0; /* nm -> angstrom */
1174 sprintf(buf, "%s%%6.f%%6.2f\n", get_pdbformat());
1175 for (i = 0; (i < 10); i++)
1176 {
1177 fprintf(fp, buf, "ATOM ", atoms->nr+1+i, "CA", "LEG", ' ',
1178 atoms->nres+1, ' ', x0, y0, z0+(1.2*i), 0.0, -0.1*i);
1179 }
1180 gmx_ffclose(fp);
1181 }
1182
1183 fprintf(log, "Dihedrals with S2 > 0.8\n");
1184 fprintf(log, "Dihedral: ");
1185 if (bPhi)
1186 {
1187 fprintf(log, " Phi ");
1188 }
1189 if (bPsi)
1190 {
1191 fprintf(log, " Psi ");
1192 }
1193 if (bChi)
1194 {
1195 for (Xi = 0; (Xi < maxchi); Xi++)
1196 {
1197 fprintf(log, " %s ", leg[2+NONCHI+Xi]);
1198 }
1199 }
1200 fprintf(log, "\nNumber: ");
1201 if (bPhi)
1202 {
1203 fprintf(log, "%4d ", nh[0]);
1204 }
1205 if (bPsi)
1206 {
1207 fprintf(log, "%4d ", nh[1]);
1208 }
1209 if (bChi)
1210 {
1211 for (Xi = 0; (Xi < maxchi); Xi++)
1212 {
1213 fprintf(log, "%4d ", nh[NONCHI+Xi]);
1214 }
1215 }
1216 fprintf(log, "\n");
1217
1218 for (i = 0; i < NLEG; i++)
1219 {
1220 sfree(leg[i]);
1221 }
1222
1223 }
1224
1225 int gmx_chi(int argc, char *argv[])
1226 {
1227 const char *desc[] = {
1228 "[THISMODULE] computes [GRK]phi[grk], [GRK]psi[grk], [GRK]omega[grk],",
1229 "and [GRK]chi[grk] dihedrals for all your",
1230 "amino acid backbone and sidechains.",
1231 "It can compute dihedral angle as a function of time, and as",
1232 "histogram distributions.",
1233 "The distributions [TT](histo-(dihedral)(RESIDUE).xvg[tt]) are cumulative over all residues of each type.[PAR]",
1234 "If option [TT]-corr[tt] is given, the program will",
1235 "calculate dihedral autocorrelation functions. The function used",
1236 "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",
1237 "rather than angles themselves, resolves the problem of periodicity.",
1238 "(Van der Spoel & Berendsen (1997), Biophys. J. 72, 2032-2041).",
1239 "Separate files for each dihedral of each residue",
1240 "[TT](corr(dihedral)(RESIDUE)(nresnr).xvg[tt]) are output, as well as a",
1241 "file containing the information for all residues (argument of [TT]-corr[tt]).[PAR]",
1242 "With option [TT]-all[tt], the angles themselves as a function of time for",
1243 "each residue are printed to separate files [TT](dihedral)(RESIDUE)(nresnr).xvg[tt].",
1244 "These can be in radians or degrees.[PAR]",
1245 "A log file (argument [TT]-g[tt]) is also written. This contains [BR]",
1246 "(a) information about the number of residues of each type.[BR]",
1247 "(b) The NMR ^3J coupling constants from the Karplus equation.[BR]",
1248 "(c) a table for each residue of the number of transitions between ",
1249 "rotamers per nanosecond, and the order parameter S^2 of each dihedral.[BR]",
1250 "(d) a table for each residue of the rotamer occupancy.[PAR]",
1251 "All rotamers are taken as 3-fold, except for [GRK]omega[grk] and [GRK]chi[grk] dihedrals",
1252 "to planar groups (i.e. [GRK]chi[grk][SUB]2[sub] of aromatics, Asp and Asn; [GRK]chi[grk][SUB]3[sub] of Glu",
1253 "and Gln; and [GRK]chi[grk][SUB]4[sub] of Arg), which are 2-fold. \"rotamer 0\" means ",
1254 "that the dihedral was not in the core region of each rotamer. ",
1255 "The width of the core region can be set with [TT]-core_rotamer[tt][PAR]",
1256
1257 "The S^2 order parameters are also output to an [TT].xvg[tt] file",
1258 "(argument [TT]-o[tt] ) and optionally as a [TT].pdb[tt] file with",
1259 "the S^2 values as B-factor (argument [TT]-p[tt]). ",
1260 "The total number of rotamer transitions per timestep",
1261 "(argument [TT]-ot[tt]), the number of transitions per rotamer",
1262 "(argument [TT]-rt[tt]), and the ^3J couplings (argument [TT]-jc[tt]), ",
1263 "can also be written to [TT].xvg[tt] files. Note that the analysis",
1264 "of rotamer transitions assumes that the supplied trajectory frames",
1265 "are equally spaced in time.[PAR]",
1266
1267 "If [TT]-chi_prod[tt] is set (and [TT]-maxchi[tt] > 0), cumulative rotamers, e.g.",
1268 "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 ",
1269 "dihedrals and [TT]-maxchi[tt] >= 3)",
1270 "are calculated. As before, if any dihedral is not in the core region,",
1271 "the rotamer is taken to be 0. The occupancies of these cumulative ",
1272 "rotamers (starting with rotamer 0) are written to the file",
1273 "that is the argument of [TT]-cp[tt], and if the [TT]-all[tt] flag",
1274 "is given, the rotamers as functions of time",
1275 "are written to [TT]chiproduct(RESIDUE)(nresnr).xvg[tt] ",
1276 "and their occupancies to [TT]histo-chiproduct(RESIDUE)(nresnr).xvg[tt].[PAR]",
1277
1278 "The option [TT]-r[tt] generates a contour plot of the average [GRK]omega[grk] angle",
1279 "as a function of the [GRK]phi[grk] and [GRK]psi[grk] angles, that is, in a Ramachandran plot",
1280 "the average [GRK]omega[grk] angle is plotted using color coding.",
1281
1282 };
1283
1284 const char *bugs[] = {
1285 "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.",
1286 "[GRK]phi[grk] and [GRK]psi[grk] dihedrals are calculated in a "
1287 "non-standard way, using H-N-CA-C for [GRK]phi[grk] instead of "
1288 "C(-)-N-CA-C, and N-CA-C-O for [GRK]psi[grk] instead of N-CA-C-N(+). "
1289 "This causes (usually small) discrepancies with the output of other "
1290 "tools like [gmx-rama].",
1291 "[TT]-r0[tt] option does not work properly",
1292 "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"
1293 };
1294
1295 /* defaults */
1296 static int r0 = 1, ndeg = 1, maxchi = 2;
1297 static gmx_bool bAll = FALSE;
1298 static gmx_bool bPhi = FALSE, bPsi = FALSE, bOmega = FALSE;
1299 static real bfac_init = -1.0, bfac_max = 0;
1300 static const char *maxchistr[] = { NULL, "0", "1", "2", "3", "4", "5", "6", NULL };
1301 static gmx_bool bRama = FALSE, bShift = FALSE, bViol = FALSE, bRamOmega = FALSE;
1302 static gmx_bool bNormHisto = TRUE, bChiProduct = FALSE, bHChi = FALSE, bRAD = FALSE, bPBC = TRUE;
1303 static real core_frac = 0.5;
1304 t_pargs pa[] = {
1305 { "-r0", FALSE, etINT, {&r0},
1306 "starting residue" },
1307 { "-phi", FALSE, etBOOL, {&bPhi},
1308 "Output for [GRK]phi[grk] dihedral angles" },
1309 { "-psi", FALSE, etBOOL, {&bPsi},
1310 "Output for [GRK]psi[grk] dihedral angles" },
1311 { "-omega", FALSE, etBOOL, {&bOmega},
1312 "Output for [GRK]omega[grk] dihedrals (peptide bonds)" },
1313 { "-rama", FALSE, etBOOL, {&bRama},
1314 "Generate [GRK]phi[grk]/[GRK]psi[grk] and [GRK]chi[grk][SUB]1[sub]/[GRK]chi[grk][SUB]2[sub] Ramachandran plots" },
1315 { "-viol", FALSE, etBOOL, {&bViol},
1316 "Write a file that gives 0 or 1 for violated Ramachandran angles" },
1317 { "-periodic", FALSE, etBOOL, {&bPBC},
1318 "Print dihedral angles modulo 360 degrees" },
1319 { "-all", FALSE, etBOOL, {&bAll},
1320 "Output separate files for every dihedral." },
1321 { "-rad", FALSE, etBOOL, {&bRAD},
1322 "in angle vs time files, use radians rather than degrees."},
1323 { "-shift", FALSE, etBOOL, {&bShift},
1324 "Compute chemical shifts from [GRK]phi[grk]/[GRK]psi[grk] angles" },
1325 { "-binwidth", FALSE, etINT, {&ndeg},
1326 "bin width for histograms (degrees)" },
1327 { "-core_rotamer", FALSE, etREAL, {&core_frac},
1328 "only the central [TT]-core_rotamer[tt]*(360/multiplicity) belongs to each rotamer (the rest is assigned to rotamer 0)" },
1329 { "-maxchi", FALSE, etENUM, {maxchistr},
1330 "calculate first ndih [GRK]chi[grk] dihedrals" },
1331 { "-normhisto", FALSE, etBOOL, {&bNormHisto},
1332 "Normalize histograms" },
1333 { "-ramomega", FALSE, etBOOL, {&bRamOmega},
1334 "compute average omega as a function of [GRK]phi[grk]/[GRK]psi[grk] and plot it in an [TT].xpm[tt] plot" },
1335 { "-bfact", FALSE, etREAL, {&bfac_init},
1336 "B-factor value for [TT].pdb[tt] file for atoms with no calculated dihedral order parameter"},
1337 { "-chi_prod", FALSE, etBOOL, {&bChiProduct},
1338 "compute a single cumulative rotamer for each residue"},
1339 { "-HChi", FALSE, etBOOL, {&bHChi},
1340 "Include dihedrals to sidechain hydrogens"},
1341 { "-bmax", FALSE, etREAL, {&bfac_max},
1342 "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." }
1343 };
1344
1345 FILE *log;
1346 int natoms, nlist, idum, nbin;
1347 t_atoms atoms;
1348 rvec *x;
1349 int ePBC;
1350 matrix box;
1351 char title[256], grpname[256];
1352 t_dlist *dlist;
1353 gmx_bool bChi, bCorr, bSSHisto;
1354 gmx_bool bDo_rt, bDo_oh, bDo_ot, bDo_jc;
1355 real dt = 0, traj_t_ns;
1356 output_env_t oenv;
1357 gmx_residuetype_t rt;
1358
1359 atom_id isize, *index;
1360 int ndih, nactdih, nf;
1361 real **dih, *trans_frac, *aver_angle, *time;
1362 int i, j, **chi_lookup, *multiplicity;
1363
1364 t_filenm fnm[] = {
1365 { efSTX, "-s", NULL, ffREAD },
1366 { efTRX, "-f", NULL, ffREAD },
1367 { efXVG, "-o", "order", ffWRITE },
1368 { efPDB, "-p", "order", ffOPTWR },
1369 { efDAT, "-ss", "ssdump", ffOPTRD },
1370 { efXVG, "-jc", "Jcoupling", ffWRITE },
1371 { efXVG, "-corr", "dihcorr", ffOPTWR },
1372 { efLOG, "-g", "chi", ffWRITE },
1373 /* add two more arguments copying from g_angle */
1374 { efXVG, "-ot", "dihtrans", ffOPTWR },
1375 { efXVG, "-oh", "trhisto", ffOPTWR },
1376 { efXVG, "-rt", "restrans", ffOPTWR },
1377 { efXVG, "-cp", "chiprodhisto", ffOPTWR }
1378 };
1379 #define NFILE asize(fnm)
1380 int npargs;
1381 t_pargs *ppa;
1382
1383 npargs = asize(pa);
1384 ppa = add_acf_pargs(&npargs, pa);
1385 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
1386 NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs,
1387 &oenv))
1388 {
1389 return 0;
1390 }
1391
1392 /* Handle result from enumerated type */
1393 sscanf(maxchistr[0], "%d", &maxchi);
1394 bChi = (maxchi > 0);
1395
1396 log = gmx_ffopen(ftp2fn(efLOG, NFILE, fnm), "w");
1397
1398 if (bRamOmega)
1399 {
1400 bOmega = TRUE;
1401 bPhi = TRUE;
1402 bPsi = TRUE;
1403 }
1404
1405 /* set some options */
1406 bDo_rt = (opt2bSet("-rt", NFILE, fnm));
1407 bDo_oh = (opt2bSet("-oh", NFILE, fnm));
1408 bDo_ot = (opt2bSet("-ot", NFILE, fnm));
1409 bDo_jc = (opt2bSet("-jc", NFILE, fnm));
1410 bCorr = (opt2bSet("-corr", NFILE, fnm));
1411 if (bCorr)
1412 {
1413 fprintf(stderr, "Will calculate autocorrelation\n");
1414 }
1415
1416 if (core_frac > 1.0)
1417 {
1418 fprintf(stderr, "core_rotamer fraction > 1.0 ; will use 1.0\n");
1419 core_frac = 1.0;
1420 }
1421 if (core_frac < 0.0)
1422 {
1423 fprintf(stderr, "core_rotamer fraction < 0.0 ; will use 0.0\n");
1424 core_frac = 0.0;
1425 }
1426
1427 if (maxchi > MAXCHI)
1428 {
1429 fprintf(stderr,
1430 "Will only calculate first %d Chi dihedrals in stead of %d.\n",
1431 MAXCHI, maxchi);
1432 maxchi = MAXCHI;
1433 }
1434 bSSHisto = ftp2bSet(efDAT, NFILE, fnm);
1435 nbin = 360/ndeg;
1436
1437 /* Find the chi angles using atoms struct and a list of amino acids */
1438 get_stx_coordnum(ftp2fn(efSTX, NFILE, fnm), &natoms);
1439 init_t_atoms(&atoms, natoms, TRUE);
1440 snew(x, natoms);
1441 read_stx_conf(ftp2fn(efSTX, NFILE, fnm), title, &atoms, x, NULL, &ePBC, box);
1442 fprintf(log, "Title: %s\n", title);
1443
1444 gmx_residuetype_init(&rt);
1445 dlist = mk_dlist(log, &atoms, &nlist, bPhi, bPsi, bChi, bHChi, maxchi, r0, rt);
1446 fprintf(stderr, "%d residues with dihedrals found\n", nlist);
1447
1448 if (nlist == 0)
1449 {
1450 gmx_fatal(FARGS, "No dihedrals in your structure!\n");
1451 }
1452
1453 /* Make a linear index for reading all. */
1454 index = make_chi_ind(nlist, dlist, &ndih);
1455 isize = 4*ndih;
1456 fprintf(stderr, "%d dihedrals found\n", ndih);
1457
1458 snew(dih, ndih);
1459
1460 /* COMPUTE ALL DIHEDRALS! */
1461 read_ang_dih(ftp2fn(efTRX, NFILE, fnm), FALSE, TRUE, FALSE, bPBC, 1, &idum,
1462 &nf, &time, isize, index, &trans_frac, &aver_angle, dih, oenv);
1463
1464 dt = (time[nf-1]-time[0])/(nf-1); /* might want this for corr or n. transit*/
1465 if (bCorr)
1466 {
1467 if (nf < 2)
1468 {
1469 gmx_fatal(FARGS, "Need at least 2 frames for correlation");
1470 }
1471 }
1472
1473 /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi
1474 * pass nactdih instead of ndih to low_ana_dih_trans and get_chi_product_traj
1475 * to prevent accessing off end of arrays when maxchi < 5 or 6. */
1476 nactdih = reset_em_all(nlist, dlist, nf, dih, maxchi);
1477
1478 if (bAll)
1479 {
1480 dump_em_all(nlist, dlist, nf, time, dih, maxchi, bPhi, bPsi, bChi, bOmega, bRAD, oenv);
1481 }
1482
1483 /* Histogramming & J coupling constants & calc of S2 order params */
1484 histogramming(log, nbin, rt, nf, maxchi, dih, nlist, dlist, index,
1485 bPhi, bPsi, bOmega, bChi,
1486 bNormHisto, bSSHisto, ftp2fn(efDAT, NFILE, fnm), bfac_max, &atoms,
1487 bDo_jc, opt2fn("-jc", NFILE, fnm), oenv);
1488
1489 /* transitions
1490 *
1491 * added multiplicity */
1492
1493 snew(multiplicity, ndih);
1494 mk_multiplicity_lookup(multiplicity, maxchi, nlist, dlist, ndih);
1495
1496 strcpy(grpname, "All residues, ");
1497 if (bPhi)
1498 {
1499 strcat(grpname, "Phi ");
1500 }
1501 if (bPsi)
1502 {
1503 strcat(grpname, "Psi ");
1504 }
1505 if (bOmega)
1506 {
1507 strcat(grpname, "Omega ");
1508 }
1509 if (bChi)
1510 {
1511 strcat(grpname, "Chi 1-");
1512 sprintf(grpname + strlen(grpname), "%i", maxchi);
1513 }
1514
1515
1516 low_ana_dih_trans(bDo_ot, opt2fn("-ot", NFILE, fnm),
1517 bDo_oh, opt2fn("-oh", NFILE, fnm), maxchi,
1518 dih, nlist, dlist, nf, nactdih, grpname, multiplicity,
1519 time, FALSE, core_frac, oenv);
1520
1521 /* Order parameters */
1522 order_params(log, opt2fn("-o", NFILE, fnm), maxchi, nlist, dlist,
1523 ftp2fn_null(efPDB, NFILE, fnm), bfac_init,
1524 &atoms, x, ePBC, box, bPhi, bPsi, bChi, oenv);
1525
1526 /* Print ramachandran maps! */
1527 if (bRama)
1528 {
1529 do_rama(nf, nlist, dlist, dih, bViol, bRamOmega, oenv);
1530 }
1531
1532 if (bShift)
1533 {
1534 do_pp2shifts(log, nf, nlist, dlist, dih);
1535 }
1536
1537 /* rprint S^2, transitions, and rotamer occupancies to log */
1538 traj_t_ns = 0.001 * (time[nf-1]-time[0]);
1539 pr_dlist(log, nlist, dlist, traj_t_ns, edPrintST, bPhi, bPsi, bChi, bOmega, maxchi);
1540 pr_dlist(log, nlist, dlist, traj_t_ns, edPrintRO, bPhi, bPsi, bChi, bOmega, maxchi);
1541 gmx_ffclose(log);
1542 /* transitions to xvg */
1543 if (bDo_rt)
1544 {
1545 print_transitions(opt2fn("-rt", NFILE, fnm), maxchi, nlist, dlist, traj_t_ns, oenv);
1546 }
1547
1548 /* chi_product trajectories (ie one "rotamer number" for each residue) */
1549 if (bChiProduct && bChi)
1550 {
1551 snew(chi_lookup, nlist);
1552 for (i = 0; i < nlist; i++)
1553 {
1554 snew(chi_lookup[i], maxchi);
1555 }
1556 mk_chi_lookup(chi_lookup, maxchi, nlist, dlist);
1557
1558 get_chi_product_traj(dih, nf, nactdih,
1559 maxchi, dlist, time, chi_lookup, multiplicity,
1560 FALSE, bNormHisto, core_frac, bAll,
1561 opt2fn("-cp", NFILE, fnm), oenv);
1562
1563 for (i = 0; i < nlist; i++)
1564 {
1565 sfree(chi_lookup[i]);
1566 }
1567 }
1568
1569 /* Correlation comes last because it messes up the angles */
1570 if (bCorr)
1571 {
1572 do_dihcorr(opt2fn("-corr", NFILE, fnm), nf, ndih, dih, dt, nlist, dlist, time,
1573 maxchi, bPhi, bPsi, bChi, bOmega, oenv);
1574 }
1575
1576
1577 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1578 do_view(oenv, opt2fn("-jc", NFILE, fnm), "-nxy");
1579 if (bCorr)
1580 {
1581 do_view(oenv, opt2fn("-corr", NFILE, fnm), "-nxy");
1582 }
1583
1584 gmx_residuetype_destroy(rt);
1585
1586 return 0;
1587 }
The ultimate molecular dynamics simulation package