?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Change gmx_bool from int to bool
[gromacs.git] / src / gromacs / gmxana / gmx_chi.cpp
blob131fb59fa379622753778f4d000f5e13104f651f
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,2015,2016,2017,2018, 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 #include "gmxpre.h"
38
39 #include <cmath>
40 #include <cstdio>
41 #include <cstring>
42
43 #include <algorithm>
44
45 #include "gromacs/commandline/pargs.h"
46 #include "gromacs/commandline/viewit.h"
47 #include "gromacs/correlationfunctions/autocorr.h"
48 #include "gromacs/fileio/confio.h"
49 #include "gromacs/fileio/matio.h"
50 #include "gromacs/fileio/pdbio.h"
51 #include "gromacs/fileio/xvgr.h"
52 #include "gromacs/gmxana/gmx_ana.h"
53 #include "gromacs/gmxana/gstat.h"
54 #include "gromacs/math/units.h"
55 #include "gromacs/math/utilities.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/topology/residuetypes.h"
58 #include "gromacs/topology/topology.h"
59 #include "gromacs/utility/arraysize.h"
60 #include "gromacs/utility/cstringutil.h"
61 #include "gromacs/utility/fatalerror.h"
62 #include "gromacs/utility/futil.h"
63 #include "gromacs/utility/smalloc.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) (((static_cast<int> (360+(ppp)*RAD2DEG)) % 360)/6)
134 x = INDEX(phi);
135 y = INDEX(psi);
136 #undef INDEX
137
138 return map[x][y] == '1';
139 }
140
141 static int *make_chi_ind(int nl, t_dlist dl[], int *ndih)
142 {
143 int *id;
144 int i, Xi, n;
145
146 /* There are nl residues with max edMax dihedrals with 4 atoms each */
147 snew(id, nl*edMax*4);
148
149 n = 0;
150 for (i = 0; (i < nl); i++)
151 {
152 /* Phi, fake the first one */
153 dl[i].j0[edPhi] = n/4;
154 if (dl[i].atm.minC >= 0)
155 {
156 id[n++] = dl[i].atm.minC;
157 }
158 else
159 {
160 id[n++] = dl[i].atm.H;
161 }
162 id[n++] = dl[i].atm.N;
163 id[n++] = dl[i].atm.Cn[1];
164 id[n++] = dl[i].atm.C;
165 }
166 for (i = 0; (i < nl); i++)
167 {
168 /* Psi, fake the last one */
169 dl[i].j0[edPsi] = n/4;
170 id[n++] = dl[i].atm.N;
171 id[n++] = dl[i].atm.Cn[1];
172 id[n++] = dl[i].atm.C;
173 if (i < (nl-1) )
174 {
175 id[n++] = dl[i+1].atm.N;
176 }
177 else
178 {
179 id[n++] = dl[i].atm.O;
180 }
181 }
182 for (i = 0; (i < nl); i++)
183 {
184 /* Omega */
185 if (has_dihedral(edOmega, &(dl[i])))
186 {
187 dl[i].j0[edOmega] = n/4;
188 id[n++] = dl[i].atm.minCalpha;
189 id[n++] = dl[i].atm.minC;
190 id[n++] = dl[i].atm.N;
191 id[n++] = dl[i].atm.Cn[1];
192 }
193 }
194 for (Xi = 0; (Xi < MAXCHI); Xi++)
195 {
196 /* Chi# */
197 for (i = 0; (i < nl); i++)
198 {
199 if (dl[i].atm.Cn[Xi+3] != -1)
200 {
201 dl[i].j0[edChi1+Xi] = n/4;
202 id[n++] = dl[i].atm.Cn[Xi];
203 id[n++] = dl[i].atm.Cn[Xi+1];
204 id[n++] = dl[i].atm.Cn[Xi+2];
205 id[n++] = dl[i].atm.Cn[Xi+3];
206 }
207 }
208 }
209 *ndih = n/4;
210
211 return id;
212 }
213
214 static void do_dihcorr(const char *fn, int nf, int ndih, real **dih, real dt,
215 int nlist, t_dlist dlist[], real time[], int maxchi,
216 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bOmega,
217 const gmx_output_env_t *oenv)
218 {
219 char name1[256], name2[256];
220 int i, j, Xi;
221
222 do_autocorr(fn, oenv, "Dihedral Autocorrelation Function",
223 nf, ndih, dih, dt, eacCos, FALSE);
224 /* Dump em all */
225 j = 0;
226 for (i = 0; (i < nlist); i++)
227 {
228 if (bPhi)
229 {
230 print_one(oenv, "corrphi", dlist[i].name, "Phi ACF for", "C(t)", nf/2, time,
231 dih[j]);
232 }
233 j++;
234 }
235 for (i = 0; (i < nlist); i++)
236 {
237 if (bPsi)
238 {
239 print_one(oenv, "corrpsi", dlist[i].name, "Psi ACF for", "C(t)", nf/2, time,
240 dih[j]);
241 }
242 j++;
243 }
244 for (i = 0; (i < nlist); i++)
245 {
246 if (has_dihedral(edOmega, &dlist[i]))
247 {
248 if (bOmega)
249 {
250 print_one(oenv, "corromega", dlist[i].name, "Omega ACF for", "C(t)",
251 nf/2, time, dih[j]);
252 }
253 j++;
254 }
255 }
256 for (Xi = 0; (Xi < maxchi); Xi++)
257 {
258 sprintf(name1, "corrchi%d", Xi+1);
259 sprintf(name2, "Chi%d ACF for", Xi+1);
260 for (i = 0; (i < nlist); i++)
261 {
262 if (dlist[i].atm.Cn[Xi+3] != -1)
263 {
264 if (bChi)
265 {
266 print_one(oenv, name1, dlist[i].name, name2, "C(t)", nf/2, time, dih[j]);
267 }
268 j++;
269 }
270 }
271 }
272 fprintf(stderr, "\n");
273 }
274
275 static void copy_dih_data(const real in[], real out[], int nf, gmx_bool bLEAVE)
276 {
277 /* if bLEAVE, do nothing to data in copying to out
278 * otherwise multiply by 180/pi to convert rad to deg */
279 int i;
280 real mult;
281 if (bLEAVE)
282 {
283 mult = 1;
284 }
285 else
286 {
287 mult = (180.0/M_PI);
288 }
289 for (i = 0; (i < nf); i++)
290 {
291 out[i] = in[i]*mult;
292 }
293 }
294
295 static void dump_em_all(int nlist, t_dlist dlist[], int nf, real time[],
296 real **dih, int maxchi,
297 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bOmega, gmx_bool bRAD,
298 const gmx_output_env_t *oenv)
299 {
300 char name[256], titlestr[256], ystr[256];
301 real *data;
302 int i, j, Xi;
303
304 snew(data, nf);
305 if (bRAD)
306 {
307 std::strcpy(ystr, "Angle (rad)");
308 }
309 else
310 {
311 std::strcpy(ystr, "Angle (degrees)");
312 }
313
314 /* Dump em all */
315 j = 0;
316 for (i = 0; (i < nlist); i++)
317 {
318 /* grs debug printf("OK i %d j %d\n", i, j) ; */
319 if (bPhi)
320 {
321 copy_dih_data(dih[j], data, nf, bRAD);
322 print_one(oenv, "phi", dlist[i].name, "\\xf\\f{}", ystr, nf, time, data);
323 }
324 j++;
325 }
326 for (i = 0; (i < nlist); i++)
327 {
328 if (bPsi)
329 {
330 copy_dih_data(dih[j], data, nf, bRAD);
331 print_one(oenv, "psi", dlist[i].name, "\\xy\\f{}", ystr, nf, time, data);
332 }
333 j++;
334 }
335 for (i = 0; (i < nlist); i++)
336 {
337 if (has_dihedral(edOmega, &(dlist[i])))
338 {
339 if (bOmega)
340 {
341 copy_dih_data(dih[j], data, nf, bRAD);
342 print_one(oenv, "omega", dlist[i].name, "\\xw\\f{}", ystr, nf, time, data);
343 }
344 j++;
345 }
346 }
347
348 for (Xi = 0; (Xi < maxchi); Xi++)
349 {
350 for (i = 0; (i < nlist); i++)
351 {
352 if (dlist[i].atm.Cn[Xi+3] != -1)
353 {
354 if (bChi)
355 {
356 sprintf(name, "chi%d", Xi+1);
357 sprintf(titlestr, "\\xc\\f{}\\s%d\\N", Xi+1);
358 copy_dih_data(dih[j], data, nf, bRAD);
359 print_one(oenv, name, dlist[i].name, titlestr, ystr, nf, time, data);
360 }
361 j++;
362 }
363 }
364 }
365 fprintf(stderr, "\n");
366 }
367
368 static void reset_one(real dih[], int nf, real phase)
369 {
370 int j;
371
372 for (j = 0; (j < nf); j++)
373 {
374 dih[j] += phase;
375 while (dih[j] < -M_PI)
376 {
377 dih[j] += 2*M_PI;
378 }
379 while (dih[j] >= M_PI)
380 {
381 dih[j] -= 2*M_PI;
382 }
383 }
384 }
385
386 static int reset_em_all(int nlist, t_dlist dlist[], int nf,
387 real **dih, int maxchi)
388 {
389 int i, j, Xi;
390
391 /* Reset em all */
392 j = 0;
393 /* Phi */
394 for (i = 0; (i < nlist); i++)
395 {
396 if (dlist[i].atm.minC == -1)
397 {
398 reset_one(dih[j++], nf, M_PI);
399 }
400 else
401 {
402 reset_one(dih[j++], nf, 0);
403 }
404 }
405 /* Psi */
406 for (i = 0; (i < nlist-1); i++)
407 {
408 reset_one(dih[j++], nf, 0);
409 }
410 /* last Psi is faked from O */
411 reset_one(dih[j++], nf, M_PI);
412
413 /* Omega */
414 for (i = 0; (i < nlist); i++)
415 {
416 if (has_dihedral(edOmega, &dlist[i]))
417 {
418 reset_one(dih[j++], nf, 0);
419 }
420 }
421 /* Chi 1 thru maxchi */
422 for (Xi = 0; (Xi < maxchi); Xi++)
423 {
424 for (i = 0; (i < nlist); i++)
425 {
426 if (dlist[i].atm.Cn[Xi+3] != -1)
427 {
428 reset_one(dih[j], nf, 0);
429 j++;
430 }
431 }
432 }
433 fprintf(stderr, "j after resetting (nr. active dihedrals) = %d\n", j);
434 return j;
435 }
436
437 static void histogramming(FILE *log, int nbin, gmx_residuetype_t *rt,
438 int nf, int maxchi, real **dih,
439 int nlist, t_dlist dlist[],
440 const int index[],
441 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bOmega, gmx_bool bChi,
442 gmx_bool bNormalize, gmx_bool bSSHisto, const char *ssdump,
443 real bfac_max, const t_atoms *atoms,
444 gmx_bool bDo_jc, const char *fn,
445 const gmx_output_env_t *oenv)
446 {
447 /* also gets 3J couplings and order parameters S2 */
448 // Avoid warnings about narrowing conversions from double to real
449 #ifdef _MSC_VER
450 #pragma warning(disable: 4838)
451 #endif
452 t_karplus kkkphi[] = {
453 { "J_NHa1", 6.51, -1.76, 1.6, -M_PI/3, 0.0, 0.0 },
454 { "J_NHa2", 6.51, -1.76, 1.6, M_PI/3, 0.0, 0.0 },
455 { "J_HaC'", 4.0, 1.1, 0.1, 0.0, 0.0, 0.0 },
456 { "J_NHCb", 4.7, -1.5, -0.2, M_PI/3, 0.0, 0.0 },
457 { "J_Ci-1Hai", 4.5, -1.3, -1.2, 2*M_PI/3, 0.0, 0.0 }
458 };
459 t_karplus kkkpsi[] = {
460 { "J_HaN", -0.88, -0.61, -0.27, M_PI/3, 0.0, 0.0 }
461 };
462 t_karplus kkkchi1[] = {
463 { "JHaHb2", 9.5, -1.6, 1.8, -M_PI/3, 0, 0.0 },
464 { "JHaHb3", 9.5, -1.6, 1.8, 0, 0, 0.0 }
465 };
466 #ifdef _MSC_VER
467 #pragma warning(default: 4838)
468 #endif
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] = {nullptr, nullptr, nullptr};
475 const char *sss[3] = { "sheet", "helix", "coil" };
476 real S2;
477 real *normhisto;
478 real **Jc, **Jcsig;
479 int ****his_aa_ss = nullptr;
480 int ***his_aa, *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 = nullptr;
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) || bBfac))
565 {
566 hindex = static_cast<int>(((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, nullptr, &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] = gmx_strdup(kkkphi[i].name);
687 }
688 for (i = 0; (i < NKKKPSI); i++)
689 {
690 leg[i+NKKKPHI] = gmx_strdup(kkkpsi[i].name);
691 }
692 for (i = 0; (i < NKKKCHI); i++)
693 {
694 leg[i+NKKKPHI+NKKKPSI] = gmx_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 xvgrclose(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 std::strcpy(hhisfile, hisfile);
765 std::strcat(hhisfile, ".xvg");
766 fp = xvgropen(hhisfile, title, "Degrees", "", oenv);
767 if (output_env_get_print_xvgr_codes(oenv))
768 {
769 fprintf(fp, "@ with g0\n");
770 }
771 xvgr_world(fp, -180, 0, 180, 0.1, oenv);
772 if (output_env_get_print_xvgr_codes(oenv))
773 {
774 fprintf(fp, "# this effort to set graph size fails unless you run with -autoscale none or -autoscale y flags\n");
775 fprintf(fp, "@ xaxis tick on\n");
776 fprintf(fp, "@ xaxis tick major 90\n");
777 fprintf(fp, "@ xaxis tick minor 30\n");
778 fprintf(fp, "@ xaxis ticklabel prec 0\n");
779 fprintf(fp, "@ yaxis tick off\n");
780 fprintf(fp, "@ yaxis ticklabel off\n");
781 fprintf(fp, "@ type xy\n");
782 }
783 if (bSSHisto)
784 {
785 for (k = 0; (k < 3); k++)
786 {
787 sprintf(sshisfile, "%s-%s.xvg", hisfile, sss[k]);
788 ssfp[k] = gmx_ffopen(sshisfile, "w");
789 }
790 }
791 for (j = 0; (j < nbin); j++)
792 {
793 angle = -180 + (360/nbin)*j;
794 if (bNormalize)
795 {
796 fprintf(fp, "%5d %10g\n", angle, normhisto[j]);
797 }
798 else
799 {
800 fprintf(fp, "%5d %10d\n", angle, his_aa[Dih][i][j]);
801 }
802 if (bSSHisto)
803 {
804 for (k = 0; (k < 3); k++)
805 {
806 fprintf(ssfp[k], "%5d %10d\n", angle,
807 his_aa_ss[k][i][Dih][j]);
808 }
809 }
810 }
811 fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
812 xvgrclose(fp);
813 if (bSSHisto)
814 {
815 for (k = 0; (k < 3); k++)
816 {
817 fprintf(ssfp[k], "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
818 gmx_ffclose(ssfp[k]);
819 }
820 }
821 }
822 }
823 }
824 sfree(normhisto);
825
826 if (bSSHisto)
827 {
828 /* Four dimensional array... Very cool */
829 for (i = 0; (i < 3); i++)
830 {
831 for (j = 0; (j <= rt_size); j++)
832 {
833 for (Dih = 0; (Dih < edMax); Dih++)
834 {
835 sfree(his_aa_ss[i][j][Dih]);
836 }
837 sfree(his_aa_ss[i][j]);
838 }
839 sfree(his_aa_ss[i]);
840 }
841 sfree(his_aa_ss);
842 sfree(ss_str);
843 }
844 }
845
846 static FILE *rama_file(const char *fn, const char *title, const char *xaxis,
847 const char *yaxis, const gmx_output_env_t *oenv)
848 {
849 FILE *fp;
850
851 fp = xvgropen(fn, title, xaxis, yaxis, oenv);
852 if (output_env_get_print_xvgr_codes(oenv))
853 {
854 fprintf(fp, "@ with g0\n");
855 }
856 xvgr_world(fp, -180, -180, 180, 180, oenv);
857 if (output_env_get_print_xvgr_codes(oenv))
858 {
859 fprintf(fp, "@ xaxis tick on\n");
860 fprintf(fp, "@ xaxis tick major 90\n");
861 fprintf(fp, "@ xaxis tick minor 30\n");
862 fprintf(fp, "@ xaxis ticklabel prec 0\n");
863 fprintf(fp, "@ yaxis tick on\n");
864 fprintf(fp, "@ yaxis tick major 90\n");
865 fprintf(fp, "@ yaxis tick minor 30\n");
866 fprintf(fp, "@ yaxis ticklabel prec 0\n");
867 fprintf(fp, "@ s0 type xy\n");
868 fprintf(fp, "@ s0 symbol 2\n");
869 fprintf(fp, "@ s0 symbol size 0.410000\n");
870 fprintf(fp, "@ s0 symbol fill 1\n");
871 fprintf(fp, "@ s0 symbol color 1\n");
872 fprintf(fp, "@ s0 symbol linewidth 1\n");
873 fprintf(fp, "@ s0 symbol linestyle 1\n");
874 fprintf(fp, "@ s0 symbol center false\n");
875 fprintf(fp, "@ s0 symbol char 0\n");
876 fprintf(fp, "@ s0 skip 0\n");
877 fprintf(fp, "@ s0 linestyle 0\n");
878 fprintf(fp, "@ s0 linewidth 1\n");
879 fprintf(fp, "@ type xy\n");
880 }
881 return fp;
882 }
883
884 static void do_rama(int nf, int nlist, t_dlist dlist[], real **dih,
885 gmx_bool bViol, gmx_bool bRamOmega, const gmx_output_env_t *oenv)
886 {
887 FILE *fp, *gp = nullptr;
888 gmx_bool bOm;
889 char fn[256];
890 int i, j, k, Xi1, Xi2, Phi, Psi, Om = 0, nlevels;
891 constexpr int NMAT = 120;
892 real **mat = nullptr, phi, psi, omega, axis[NMAT], lo, hi;
893 t_rgb rlo = { 1.0, 0.0, 0.0 };
894 t_rgb rmid = { 1.0, 1.0, 1.0 };
895 t_rgb rhi = { 0.0, 0.0, 1.0 };
896
897 for (i = 0; (i < nlist); i++)
898 {
899 if ((has_dihedral(edPhi, &(dlist[i]))) &&
900 (has_dihedral(edPsi, &(dlist[i]))))
901 {
902 sprintf(fn, "ramaPhiPsi%s.xvg", dlist[i].name);
903 fp = rama_file(fn, "Ramachandran Plot",
904 "\\8f\\4 (deg)", "\\8y\\4 (deg)", oenv);
905 bOm = bRamOmega && has_dihedral(edOmega, &(dlist[i]));
906 if (bOm)
907 {
908 Om = dlist[i].j0[edOmega];
909 snew(mat, NMAT);
910 for (j = 0; (j < NMAT); j++)
911 {
912 snew(mat[j], NMAT);
913 axis[j] = -180+gmx::exactDiv(360*j, NMAT);
914 }
915 }
916 if (bViol)
917 {
918 sprintf(fn, "violPhiPsi%s.xvg", dlist[i].name);
919 gp = gmx_ffopen(fn, "w");
920 }
921 Phi = dlist[i].j0[edPhi];
922 Psi = dlist[i].j0[edPsi];
923 for (j = 0; (j < nf); j++)
924 {
925 phi = RAD2DEG*dih[Phi][j];
926 psi = RAD2DEG*dih[Psi][j];
927 fprintf(fp, "%10g %10g\n", phi, psi);
928 if (bViol)
929 {
930 fprintf(gp, "%d\n", int{(bAllowed(dih[Phi][j], RAD2DEG*dih[Psi][j]) == FALSE)} );
931 }
932 if (bOm)
933 {
934 omega = RAD2DEG*dih[Om][j];
935 mat[static_cast<int>(((phi*NMAT)/360)+gmx::exactDiv(NMAT, 2))][static_cast<int>(((psi*NMAT)/360)+gmx::exactDiv(NMAT, 2))]
936 += omega;
937 }
938 }
939 if (bViol)
940 {
941 gmx_ffclose(gp);
942 }
943 xvgrclose(fp);
944 if (bOm)
945 {
946 sprintf(fn, "ramomega%s.xpm", dlist[i].name);
947 fp = gmx_ffopen(fn, "w");
948 lo = hi = 0;
949 for (j = 0; (j < NMAT); j++)
950 {
951 for (k = 0; (k < NMAT); k++)
952 {
953 mat[j][k] /= nf;
954 lo = std::min(mat[j][k], lo);
955 hi = std::max(mat[j][k], hi);
956 }
957 }
958 /* Symmetrise */
959 if (std::abs(lo) > std::abs(hi))
960 {
961 hi = -lo;
962 }
963 else
964 {
965 lo = -hi;
966 }
967 /* Add 180 */
968 for (j = 0; (j < NMAT); j++)
969 {
970 for (k = 0; (k < NMAT); k++)
971 {
972 mat[j][k] += 180;
973 }
974 }
975 lo += 180;
976 hi += 180;
977 nlevels = 20;
978 write_xpm3(fp, 0, "Omega/Ramachandran Plot", "Deg", "Phi", "Psi",
979 NMAT, NMAT, axis, axis, mat, lo, 180.0, hi, rlo, rmid, rhi, &nlevels);
980 gmx_ffclose(fp);
981 for (j = 0; (j < NMAT); j++)
982 {
983 sfree(mat[j]);
984 }
985 sfree(mat);
986 }
987 }
988 if ((has_dihedral(edChi1, &(dlist[i]))) &&
989 (has_dihedral(edChi2, &(dlist[i]))))
990 {
991 sprintf(fn, "ramaX1X2%s.xvg", dlist[i].name);
992 fp = rama_file(fn, "\\8c\\4\\s1\\N-\\8c\\4\\s2\\N Ramachandran Plot",
993 "\\8c\\4\\s1\\N (deg)", "\\8c\\4\\s2\\N (deg)", oenv);
994 Xi1 = dlist[i].j0[edChi1];
995 Xi2 = dlist[i].j0[edChi2];
996 for (j = 0; (j < nf); j++)
997 {
998 fprintf(fp, "%10g %10g\n", RAD2DEG*dih[Xi1][j], RAD2DEG*dih[Xi2][j]);
999 }
1000 xvgrclose(fp);
1001 }
1002 else
1003 {
1004 fprintf(stderr, "No chi1 & chi2 angle for %s\n", dlist[i].name);
1005 }
1006 }
1007 }
1008
1009
1010 static void print_transitions(const char *fn, int maxchi, int nlist,
1011 t_dlist dlist[], real dt,
1012 const gmx_output_env_t *oenv)
1013 {
1014 /* based on order_params below */
1015 FILE *fp;
1016 int i, Dih, Xi;
1017
1018 /* must correspond with enum in pp2shift.h:38 */
1019 char *leg[edMax];
1020 #define NLEG asize(leg)
1021
1022 leg[0] = gmx_strdup("Phi");
1023 leg[1] = gmx_strdup("Psi");
1024 leg[2] = gmx_strdup("Omega");
1025 leg[3] = gmx_strdup("Chi1");
1026 leg[4] = gmx_strdup("Chi2");
1027 leg[5] = gmx_strdup("Chi3");
1028 leg[6] = gmx_strdup("Chi4");
1029 leg[7] = gmx_strdup("Chi5");
1030 leg[8] = gmx_strdup("Chi6");
1031
1032 /* Print order parameters */
1033 fp = xvgropen(fn, "Dihedral Rotamer Transitions", "Residue", "Transitions/ns",
1034 oenv);
1035 xvgr_legend(fp, NONCHI+maxchi, (const char**)leg, oenv);
1036
1037 fprintf(fp, "%5s ", "#Res.");
1038 fprintf(fp, "%10s %10s %10s ", leg[edPhi], leg[edPsi], leg[edOmega]);
1039 for (Xi = 0; Xi < maxchi; Xi++)
1040 {
1041 fprintf(fp, "%10s ", leg[NONCHI+Xi]);
1042 }
1043 fprintf(fp, "\n");
1044
1045 for (i = 0; (i < nlist); i++)
1046 {
1047 fprintf(fp, "%5d ", dlist[i].resnr);
1048 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1049 {
1050 fprintf(fp, "%10.3f ", dlist[i].ntr[Dih]/dt);
1051 }
1052 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
1053 fprintf(fp, "\n");
1054 }
1055 xvgrclose(fp);
1056 }
1057
1058 static void order_params(FILE *log,
1059 const char *fn, int maxchi, int nlist, t_dlist dlist[],
1060 const char *pdbfn, real bfac_init,
1061 t_atoms *atoms, const rvec x[], int ePBC, matrix box,
1062 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, const gmx_output_env_t *oenv)
1063 {
1064 FILE *fp;
1065 int nh[edMax];
1066 int i, Dih, Xi;
1067 real S2Max, S2Min;
1068
1069 /* except for S2Min/Max, must correspond with enum in pp2shift.h:38 */
1070 const char *const_leg[2+edMax] = {
1071 "S2Min", "S2Max", "Phi", "Psi", "Omega",
1072 "Chi1", "Chi2", "Chi3", "Chi4", "Chi5",
1073 "Chi6"
1074 };
1075 #define NLEG asize(leg)
1076
1077 char *leg[2+edMax];
1078
1079 for (i = 0; i < NLEG; i++)
1080 {
1081 leg[i] = gmx_strdup(const_leg[i]);
1082 }
1083
1084 /* Print order parameters */
1085 fp = xvgropen(fn, "Dihedral Order Parameters", "Residue", "S2", oenv);
1086 xvgr_legend(fp, 2+NONCHI+maxchi, const_leg, oenv);
1087
1088 for (Dih = 0; (Dih < edMax); Dih++)
1089 {
1090 nh[Dih] = 0;
1091 }
1092
1093 fprintf(fp, "%5s ", "#Res.");
1094 fprintf(fp, "%10s %10s ", leg[0], leg[1]);
1095 fprintf(fp, "%10s %10s %10s ", leg[2+edPhi], leg[2+edPsi], leg[2+edOmega]);
1096 for (Xi = 0; Xi < maxchi; Xi++)
1097 {
1098 fprintf(fp, "%10s ", leg[2+NONCHI+Xi]);
1099 }
1100 fprintf(fp, "\n");
1101
1102 for (i = 0; (i < nlist); i++)
1103 {
1104 S2Max = -10;
1105 S2Min = 10;
1106 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1107 {
1108 if (dlist[i].S2[Dih] != 0)
1109 {
1110 if (dlist[i].S2[Dih] > S2Max)
1111 {
1112 S2Max = dlist[i].S2[Dih];
1113 }
1114 if (dlist[i].S2[Dih] < S2Min)
1115 {
1116 S2Min = dlist[i].S2[Dih];
1117 }
1118 }
1119 if (dlist[i].S2[Dih] > 0.8)
1120 {
1121 nh[Dih]++;
1122 }
1123 }
1124 fprintf(fp, "%5d ", dlist[i].resnr);
1125 fprintf(fp, "%10.3f %10.3f ", S2Min, S2Max);
1126 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1127 {
1128 fprintf(fp, "%10.3f ", dlist[i].S2[Dih]);
1129 }
1130 fprintf(fp, "\n");
1131 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
1132 }
1133 xvgrclose(fp);
1134
1135 if (nullptr != pdbfn)
1136 {
1137 real x0, y0, z0;
1138
1139 atoms->havePdbInfo = TRUE;
1140
1141 if (nullptr == atoms->pdbinfo)
1142 {
1143 snew(atoms->pdbinfo, atoms->nr);
1144 }
1145 for (i = 0; (i < atoms->nr); i++)
1146 {
1147 atoms->pdbinfo[i].bfac = bfac_init;
1148 }
1149
1150 for (i = 0; (i < nlist); i++)
1151 {
1152 atoms->pdbinfo[dlist[i].atm.N].bfac = -dlist[i].S2[0]; /* Phi */
1153 atoms->pdbinfo[dlist[i].atm.H].bfac = -dlist[i].S2[0]; /* Phi */
1154 atoms->pdbinfo[dlist[i].atm.C].bfac = -dlist[i].S2[1]; /* Psi */
1155 atoms->pdbinfo[dlist[i].atm.O].bfac = -dlist[i].S2[1]; /* Psi */
1156 for (Xi = 0; (Xi < maxchi); Xi++) /* Chi's */
1157 {
1158 if (dlist[i].atm.Cn[Xi+3] != -1)
1159 {
1160 atoms->pdbinfo[dlist[i].atm.Cn[Xi+1]].bfac = -dlist[i].S2[NONCHI+Xi];
1161 }
1162 }
1163 }
1164
1165 fp = gmx_ffopen(pdbfn, "w");
1166 fprintf(fp, "REMARK generated by g_chi\n");
1167 fprintf(fp, "REMARK "
1168 "B-factor field contains negative of dihedral order parameters\n");
1169 write_pdbfile(fp, nullptr, atoms, x, ePBC, box, ' ', 0, nullptr, TRUE);
1170 x0 = y0 = z0 = 1000.0;
1171 for (i = 0; (i < atoms->nr); i++)
1172 {
1173 x0 = std::min(x0, x[i][XX]);
1174 y0 = std::min(y0, x[i][YY]);
1175 z0 = std::min(z0, x[i][ZZ]);
1176 }
1177 x0 *= 10.0; /* nm -> angstrom */
1178 y0 *= 10.0; /* nm -> angstrom */
1179 z0 *= 10.0; /* nm -> angstrom */
1180 for (i = 0; (i < 10); i++)
1181 {
1182 gmx_fprintf_pdb_atomline(fp, epdbATOM, atoms->nr+1+i, "CA", ' ', "LEG", ' ', atoms->nres+1, ' ',
1183 x0, y0, z0+(1.2*i), 0.0, -0.1*i, "");
1184 }
1185 gmx_ffclose(fp);
1186 }
1187
1188 fprintf(log, "Dihedrals with S2 > 0.8\n");
1189 fprintf(log, "Dihedral: ");
1190 if (bPhi)
1191 {
1192 fprintf(log, " Phi ");
1193 }
1194 if (bPsi)
1195 {
1196 fprintf(log, " Psi ");
1197 }
1198 if (bChi)
1199 {
1200 for (Xi = 0; (Xi < maxchi); Xi++)
1201 {
1202 fprintf(log, " %s ", leg[2+NONCHI+Xi]);
1203 }
1204 }
1205 fprintf(log, "\nNumber: ");
1206 if (bPhi)
1207 {
1208 fprintf(log, "%4d ", nh[0]);
1209 }
1210 if (bPsi)
1211 {
1212 fprintf(log, "%4d ", nh[1]);
1213 }
1214 if (bChi)
1215 {
1216 for (Xi = 0; (Xi < maxchi); Xi++)
1217 {
1218 fprintf(log, "%4d ", nh[NONCHI+Xi]);
1219 }
1220 }
1221 fprintf(log, "\n");
1222
1223 for (i = 0; i < NLEG; i++)
1224 {
1225 sfree(leg[i]);
1226 }
1227
1228 }
1229
1230 int gmx_chi(int argc, char *argv[])
1231 {
1232 const char *desc[] = {
1233 "[THISMODULE] computes [GRK]phi[grk], [GRK]psi[grk], [GRK]omega[grk],",
1234 "and [GRK]chi[grk] dihedrals for all your",
1235 "amino acid backbone and sidechains.",
1236 "It can compute dihedral angle as a function of time, and as",
1237 "histogram distributions.",
1238 "The distributions [TT](histo-(dihedral)(RESIDUE).xvg[tt]) are cumulative over all residues of each type.[PAR]",
1239 "If option [TT]-corr[tt] is given, the program will",
1240 "calculate dihedral autocorrelation functions. The function used",
1241 "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",
1242 "rather than angles themselves, resolves the problem of periodicity.",
1243 "(Van der Spoel & Berendsen (1997), Biophys. J. 72, 2032-2041).",
1244 "Separate files for each dihedral of each residue",
1245 "[TT](corr(dihedral)(RESIDUE)(nresnr).xvg[tt]) are output, as well as a",
1246 "file containing the information for all residues (argument of [TT]-corr[tt]).[PAR]",
1247 "With option [TT]-all[tt], the angles themselves as a function of time for",
1248 "each residue are printed to separate files [TT](dihedral)(RESIDUE)(nresnr).xvg[tt].",
1249 "These can be in radians or degrees.[PAR]",
1250 "A log file (argument [TT]-g[tt]) is also written. This contains",
1251 "",
1252 " * information about the number of residues of each type.",
1253 " * The NMR ^3J coupling constants from the Karplus equation.",
1254 " * a table for each residue of the number of transitions between ",
1255 " rotamers per nanosecond, and the order parameter S^2 of each dihedral.",
1256 " * a table for each residue of the rotamer occupancy.",
1257 "",
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]",
1263
1264 "The S^2 order parameters are also output to an [REF].xvg[ref] file",
1265 "(argument [TT]-o[tt] ) and optionally as a [REF].pdb[ref] 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 [REF].xvg[ref] files. Note that the analysis",
1271 "of rotamer transitions assumes that the supplied trajectory frames",
1272 "are equally spaced in time.[PAR]",
1273
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]",
1284
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.",
1288
1289 };
1290
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"
1300 };
1301
1302 /* defaults */
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[] = { nullptr, "0", "1", "2", "3", "4", "5", "6", nullptr };
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;
1311 t_pargs pa[] = {
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 [REF].xpm[ref] plot" },
1342 { "-bfact", FALSE, etREAL, {&bfac_init},
1343 "B-factor value for [REF].pdb[ref] 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." }
1350 };
1351
1352 FILE *log;
1353 int nlist, idum, nbin;
1354 rvec *x;
1355 int ePBC;
1356 matrix box;
1357 char grpname[256];
1358 t_dlist *dlist;
1359 gmx_bool bChi, bCorr, bSSHisto;
1360 gmx_bool bDo_rt, bDo_oh, bDo_ot, bDo_jc;
1361 real dt = 0, traj_t_ns;
1362 gmx_output_env_t *oenv;
1363 gmx_residuetype_t *rt;
1364
1365 int isize, *index;
1366 int ndih, nactdih, nf;
1367 real **dih, *trans_frac, *aver_angle, *time;
1368 int i, **chi_lookup, *multiplicity;
1369
1370 t_filenm fnm[] = {
1371 { efSTX, "-s", nullptr, ffREAD },
1372 { efTRX, "-f", nullptr, ffREAD },
1373 { efXVG, "-o", "order", ffWRITE },
1374 { efPDB, "-p", "order", ffOPTWR },
1375 { efDAT, "-ss", "ssdump", ffOPTRD },
1376 { efXVG, "-jc", "Jcoupling", ffWRITE },
1377 { efXVG, "-corr", "dihcorr", ffOPTWR },
1378 { efLOG, "-g", "chi", ffWRITE },
1379 /* add two more arguments copying from g_angle */
1380 { efXVG, "-ot", "dihtrans", ffOPTWR },
1381 { efXVG, "-oh", "trhisto", ffOPTWR },
1382 { efXVG, "-rt", "restrans", ffOPTWR },
1383 { efXVG, "-cp", "chiprodhisto", ffOPTWR }
1384 };
1385 #define NFILE asize(fnm)
1386 int npargs;
1387 t_pargs *ppa;
1388
1389 npargs = asize(pa);
1390 ppa = add_acf_pargs(&npargs, pa);
1391 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
1392 NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs,
1393 &oenv))
1394 {
1395 sfree(ppa);
1396 return 0;
1397 }
1398
1399 /* Handle result from enumerated type */
1400 sscanf(maxchistr[0], "%d", &maxchi);
1401 bChi = (maxchi > 0);
1402
1403 log = gmx_ffopen(ftp2fn(efLOG, NFILE, fnm), "w");
1404
1405 if (bRamOmega)
1406 {
1407 bOmega = TRUE;
1408 bPhi = TRUE;
1409 bPsi = TRUE;
1410 }
1411
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));
1418 if (bCorr)
1419 {
1420 fprintf(stderr, "Will calculate autocorrelation\n");
1421 }
1422
1423 if (core_frac > 1.0)
1424 {
1425 fprintf(stderr, "core_rotamer fraction > 1.0 ; will use 1.0\n");
1426 core_frac = 1.0;
1427 }
1428 if (core_frac < 0.0)
1429 {
1430 fprintf(stderr, "core_rotamer fraction < 0.0 ; will use 0.0\n");
1431 core_frac = 0.0;
1432 }
1433
1434 if (maxchi > MAXCHI)
1435 {
1436 fprintf(stderr,
1437 "Will only calculate first %d Chi dihedrals instead of %d.\n",
1438 MAXCHI, maxchi);
1439 maxchi = MAXCHI;
1440 }
1441 bSSHisto = ftp2bSet(efDAT, NFILE, fnm);
1442 nbin = 360/ndeg;
1443
1444 /* Find the chi angles using atoms struct and a list of amino acids */
1445 t_topology *top;
1446 snew(top, 1);
1447 read_tps_conf(ftp2fn(efSTX, NFILE, fnm), top, &ePBC, &x, nullptr, box, FALSE);
1448 t_atoms &atoms = top->atoms;
1449 if (atoms.pdbinfo == nullptr)
1450 {
1451 snew(atoms.pdbinfo, atoms.nr);
1452 }
1453 fprintf(log, "Title: %s\n", *top->name);
1454
1455 gmx_residuetype_init(&rt);
1456 dlist = mk_dlist(log, &atoms, &nlist, bPhi, bPsi, bChi, bHChi, maxchi, r0, rt);
1457 fprintf(stderr, "%d residues with dihedrals found\n", nlist);
1458
1459 if (nlist == 0)
1460 {
1461 gmx_fatal(FARGS, "No dihedrals in your structure!\n");
1462 }
1463
1464 /* Make a linear index for reading all. */
1465 index = make_chi_ind(nlist, dlist, &ndih);
1466 isize = 4*ndih;
1467 fprintf(stderr, "%d dihedrals found\n", ndih);
1468
1469 snew(dih, ndih);
1470
1471 /* COMPUTE ALL DIHEDRALS! */
1472 read_ang_dih(ftp2fn(efTRX, NFILE, fnm), FALSE, TRUE, FALSE, bPBC, 1, &idum,
1473 &nf, &time, isize, index, &trans_frac, &aver_angle, dih, oenv);
1474
1475 dt = (time[nf-1]-time[0])/(nf-1); /* might want this for corr or n. transit*/
1476 if (bCorr)
1477 {
1478 if (nf < 2)
1479 {
1480 gmx_fatal(FARGS, "Need at least 2 frames for correlation");
1481 }
1482 }
1483
1484 /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi
1485 * pass nactdih instead of ndih to low_ana_dih_trans and get_chi_product_traj
1486 * to prevent accessing off end of arrays when maxchi < 5 or 6. */
1487 nactdih = reset_em_all(nlist, dlist, nf, dih, maxchi);
1488
1489 if (bAll)
1490 {
1491 dump_em_all(nlist, dlist, nf, time, dih, maxchi, bPhi, bPsi, bChi, bOmega, bRAD, oenv);
1492 }
1493
1494 /* Histogramming & J coupling constants & calc of S2 order params */
1495 histogramming(log, nbin, rt, nf, maxchi, dih, nlist, dlist, index,
1496 bPhi, bPsi, bOmega, bChi,
1497 bNormHisto, bSSHisto, ftp2fn(efDAT, NFILE, fnm), bfac_max, &atoms,
1498 bDo_jc, opt2fn("-jc", NFILE, fnm), oenv);
1499
1500 /* transitions
1501 *
1502 * added multiplicity */
1503
1504 snew(multiplicity, ndih);
1505 mk_multiplicity_lookup(multiplicity, maxchi, nlist, dlist, ndih);
1506
1507 std::strcpy(grpname, "All residues, ");
1508 if (bPhi)
1509 {
1510 std::strcat(grpname, "Phi ");
1511 }
1512 if (bPsi)
1513 {
1514 std::strcat(grpname, "Psi ");
1515 }
1516 if (bOmega)
1517 {
1518 std::strcat(grpname, "Omega ");
1519 }
1520 if (bChi)
1521 {
1522 std::strcat(grpname, "Chi 1-");
1523 sprintf(grpname + std::strlen(grpname), "%i", maxchi);
1524 }
1525
1526
1527 low_ana_dih_trans(bDo_ot, opt2fn("-ot", NFILE, fnm),
1528 bDo_oh, opt2fn("-oh", NFILE, fnm), maxchi,
1529 dih, nlist, dlist, nf, nactdih, grpname, multiplicity,
1530 time, FALSE, core_frac, oenv);
1531
1532 /* Order parameters */
1533 order_params(log, opt2fn("-o", NFILE, fnm), maxchi, nlist, dlist,
1534 ftp2fn_null(efPDB, NFILE, fnm), bfac_init,
1535 &atoms, x, ePBC, box, bPhi, bPsi, bChi, oenv);
1536
1537 /* Print ramachandran maps! */
1538 if (bRama)
1539 {
1540 do_rama(nf, nlist, dlist, dih, bViol, bRamOmega, oenv);
1541 }
1542
1543 if (bShift)
1544 {
1545 do_pp2shifts(log, nf, nlist, dlist, dih);
1546 }
1547
1548 /* rprint S^2, transitions, and rotamer occupancies to log */
1549 traj_t_ns = 0.001 * (time[nf-1]-time[0]);
1550 pr_dlist(log, nlist, dlist, traj_t_ns, edPrintST, bPhi, bPsi, bChi, bOmega, maxchi);
1551 pr_dlist(log, nlist, dlist, traj_t_ns, edPrintRO, bPhi, bPsi, bChi, bOmega, maxchi);
1552 gmx_ffclose(log);
1553 /* transitions to xvg */
1554 if (bDo_rt)
1555 {
1556 print_transitions(opt2fn("-rt", NFILE, fnm), maxchi, nlist, dlist, traj_t_ns, oenv);
1557 }
1558
1559 /* chi_product trajectories (ie one "rotamer number" for each residue) */
1560 if (bChiProduct && bChi)
1561 {
1562 snew(chi_lookup, nlist);
1563 for (i = 0; i < nlist; i++)
1564 {
1565 snew(chi_lookup[i], maxchi);
1566 }
1567 mk_chi_lookup(chi_lookup, maxchi, nlist, dlist);
1568
1569 get_chi_product_traj(dih, nf, nactdih,
1570 maxchi, dlist, time, chi_lookup, multiplicity,
1571 FALSE, bNormHisto, core_frac, bAll,
1572 opt2fn("-cp", NFILE, fnm), oenv);
1573
1574 for (i = 0; i < nlist; i++)
1575 {
1576 sfree(chi_lookup[i]);
1577 }
1578 }
1579
1580 /* Correlation comes last because it messes up the angles */
1581 if (bCorr)
1582 {
1583 do_dihcorr(opt2fn("-corr", NFILE, fnm), nf, ndih, dih, dt, nlist, dlist, time,
1584 maxchi, bPhi, bPsi, bChi, bOmega, oenv);
1585 }
1586
1587
1588 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1589 do_view(oenv, opt2fn("-jc", NFILE, fnm), "-nxy");
1590 if (bCorr)
1591 {
1592 do_view(oenv, opt2fn("-corr", NFILE, fnm), "-nxy");
1593 }
1594
1595 gmx_residuetype_destroy(rt);
1596
1597 return 0;
1598 }
The ultimate molecular dynamics simulation package