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