?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Add replacements for pbc enumerations
[gromacs.git] / src / gromacs / gmxana / gmx_chi.cpp
blobbda1e768f9f927259e4717d3564959e363b8cafa
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,2019, 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 #include "gromacs/utility/stringutil.h"
65
66 static gmx_bool bAllowed(real phi, real psi)
67 {
68 static const char *map[] = {
69 "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 };
131 #define NPP asize(map)
132 int x, y;
133
134 #define INDEX(ppp) (((static_cast<int> (360+(ppp)*RAD2DEG)) % 360)/6)
135 x = INDEX(phi);
136 y = INDEX(psi);
137 #undef INDEX
138
139 return map[x][y] == '1';
140 }
141
142 static int *make_chi_ind(int nl, t_dlist dl[], int *ndih)
143 {
144 int *id;
145 int i, Xi, n;
146
147 /* There are nl residues with max edMax dihedrals with 4 atoms each */
148 snew(id, nl*edMax*4);
149
150 n = 0;
151 for (i = 0; (i < nl); i++)
152 {
153 /* Phi, fake the first one */
154 dl[i].j0[edPhi] = n/4;
155 if (dl[i].atm.minC >= 0)
156 {
157 id[n++] = dl[i].atm.minC;
158 }
159 else
160 {
161 id[n++] = dl[i].atm.H;
162 }
163 id[n++] = dl[i].atm.N;
164 id[n++] = dl[i].atm.Cn[1];
165 id[n++] = dl[i].atm.C;
166 }
167 for (i = 0; (i < nl); i++)
168 {
169 /* Psi, fake the last one */
170 dl[i].j0[edPsi] = n/4;
171 id[n++] = dl[i].atm.N;
172 id[n++] = dl[i].atm.Cn[1];
173 id[n++] = dl[i].atm.C;
174 if (i < (nl-1) )
175 {
176 id[n++] = dl[i+1].atm.N;
177 }
178 else
179 {
180 id[n++] = dl[i].atm.O;
181 }
182 }
183 for (i = 0; (i < nl); i++)
184 {
185 /* Omega */
186 if (has_dihedral(edOmega, &(dl[i])))
187 {
188 dl[i].j0[edOmega] = n/4;
189 id[n++] = dl[i].atm.minCalpha;
190 id[n++] = dl[i].atm.minC;
191 id[n++] = dl[i].atm.N;
192 id[n++] = dl[i].atm.Cn[1];
193 }
194 }
195 for (Xi = 0; (Xi < MAXCHI); Xi++)
196 {
197 /* Chi# */
198 for (i = 0; (i < nl); i++)
199 {
200 if (dl[i].atm.Cn[Xi+3] != -1)
201 {
202 dl[i].j0[edChi1+Xi] = n/4;
203 id[n++] = dl[i].atm.Cn[Xi];
204 id[n++] = dl[i].atm.Cn[Xi+1];
205 id[n++] = dl[i].atm.Cn[Xi+2];
206 id[n++] = dl[i].atm.Cn[Xi+3];
207 }
208 }
209 }
210 *ndih = n/4;
211
212 return id;
213 }
214
215 static void do_dihcorr(const char *fn, int nf, int ndih, real **dih, real dt,
216 int nlist, t_dlist dlist[], real time[], int maxchi,
217 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bOmega,
218 const gmx_output_env_t *oenv)
219 {
220 char name1[256], name2[256];
221 int i, j, Xi;
222
223 do_autocorr(fn, oenv, "Dihedral Autocorrelation Function",
224 nf, ndih, dih, dt, eacCos, FALSE);
225 /* Dump em all */
226 j = 0;
227 for (i = 0; (i < nlist); i++)
228 {
229 if (bPhi)
230 {
231 print_one(oenv, "corrphi", dlist[i].name, "Phi ACF for", "C(t)", nf/2, time,
232 dih[j]);
233 }
234 j++;
235 }
236 for (i = 0; (i < nlist); i++)
237 {
238 if (bPsi)
239 {
240 print_one(oenv, "corrpsi", dlist[i].name, "Psi ACF for", "C(t)", nf/2, time,
241 dih[j]);
242 }
243 j++;
244 }
245 for (i = 0; (i < nlist); i++)
246 {
247 if (has_dihedral(edOmega, &dlist[i]))
248 {
249 if (bOmega)
250 {
251 print_one(oenv, "corromega", dlist[i].name, "Omega ACF for", "C(t)",
252 nf/2, time, dih[j]);
253 }
254 j++;
255 }
256 }
257 for (Xi = 0; (Xi < maxchi); Xi++)
258 {
259 sprintf(name1, "corrchi%d", Xi+1);
260 sprintf(name2, "Chi%d ACF for", Xi+1);
261 for (i = 0; (i < nlist); i++)
262 {
263 if (dlist[i].atm.Cn[Xi+3] != -1)
264 {
265 if (bChi)
266 {
267 print_one(oenv, name1, dlist[i].name, name2, "C(t)", nf/2, time, dih[j]);
268 }
269 j++;
270 }
271 }
272 }
273 fprintf(stderr, "\n");
274 }
275
276 static void copy_dih_data(const real in[], real out[], int nf, gmx_bool bLEAVE)
277 {
278 /* if bLEAVE, do nothing to data in copying to out
279 * otherwise multiply by 180/pi to convert rad to deg */
280 int i;
281 real mult;
282 if (bLEAVE)
283 {
284 mult = 1;
285 }
286 else
287 {
288 mult = (180.0/M_PI);
289 }
290 for (i = 0; (i < nf); i++)
291 {
292 out[i] = in[i]*mult;
293 }
294 }
295
296 static void dump_em_all(int nlist, t_dlist dlist[], int nf, real time[],
297 real **dih, int maxchi,
298 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bOmega, gmx_bool bRAD,
299 const gmx_output_env_t *oenv)
300 {
301 char name[256], titlestr[256], ystr[256];
302 real *data;
303 int i, j, Xi;
304
305 snew(data, nf);
306 if (bRAD)
307 {
308 std::strcpy(ystr, "Angle (rad)");
309 }
310 else
311 {
312 std::strcpy(ystr, "Angle (degrees)");
313 }
314
315 /* Dump em all */
316 j = 0;
317 for (i = 0; (i < nlist); i++)
318 {
319 /* grs debug printf("OK i %d j %d\n", i, j) ; */
320 if (bPhi)
321 {
322 copy_dih_data(dih[j], data, nf, bRAD);
323 print_one(oenv, "phi", dlist[i].name, "\\xf\\f{}", ystr, nf, time, data);
324 }
325 j++;
326 }
327 for (i = 0; (i < nlist); i++)
328 {
329 if (bPsi)
330 {
331 copy_dih_data(dih[j], data, nf, bRAD);
332 print_one(oenv, "psi", dlist[i].name, "\\xy\\f{}", ystr, nf, time, data);
333 }
334 j++;
335 }
336 for (i = 0; (i < nlist); i++)
337 {
338 if (has_dihedral(edOmega, &(dlist[i])))
339 {
340 if (bOmega)
341 {
342 copy_dih_data(dih[j], data, nf, bRAD);
343 print_one(oenv, "omega", dlist[i].name, "\\xw\\f{}", ystr, nf, time, data);
344 }
345 j++;
346 }
347 }
348
349 for (Xi = 0; (Xi < maxchi); Xi++)
350 {
351 for (i = 0; (i < nlist); i++)
352 {
353 if (dlist[i].atm.Cn[Xi+3] != -1)
354 {
355 if (bChi)
356 {
357 sprintf(name, "chi%d", Xi+1);
358 sprintf(titlestr, "\\xc\\f{}\\s%d\\N", Xi+1);
359 copy_dih_data(dih[j], data, nf, bRAD);
360 print_one(oenv, name, dlist[i].name, titlestr, ystr, nf, time, data);
361 }
362 j++;
363 }
364 }
365 }
366 fprintf(stderr, "\n");
367 }
368
369 static void reset_one(real dih[], int nf, real phase)
370 {
371 int j;
372
373 for (j = 0; (j < nf); j++)
374 {
375 dih[j] += phase;
376 while (dih[j] < -M_PI)
377 {
378 dih[j] += 2*M_PI;
379 }
380 while (dih[j] >= M_PI)
381 {
382 dih[j] -= 2*M_PI;
383 }
384 }
385 }
386
387 static int reset_em_all(int nlist, t_dlist dlist[], int nf,
388 real **dih, int maxchi)
389 {
390 int i, j, Xi;
391
392 /* Reset em all */
393 j = 0;
394 /* Phi */
395 for (i = 0; (i < nlist); i++)
396 {
397 if (dlist[i].atm.minC == -1)
398 {
399 reset_one(dih[j++], nf, M_PI);
400 }
401 else
402 {
403 reset_one(dih[j++], nf, 0);
404 }
405 }
406 /* Psi */
407 for (i = 0; (i < nlist-1); i++)
408 {
409 reset_one(dih[j++], nf, 0);
410 }
411 /* last Psi is faked from O */
412 reset_one(dih[j++], nf, M_PI);
413
414 /* Omega */
415 for (i = 0; (i < nlist); i++)
416 {
417 if (has_dihedral(edOmega, &dlist[i]))
418 {
419 reset_one(dih[j++], nf, 0);
420 }
421 }
422 /* Chi 1 thru maxchi */
423 for (Xi = 0; (Xi < maxchi); Xi++)
424 {
425 for (i = 0; (i < nlist); i++)
426 {
427 if (dlist[i].atm.Cn[Xi+3] != -1)
428 {
429 reset_one(dih[j], nf, 0);
430 j++;
431 }
432 }
433 }
434 fprintf(stderr, "j after resetting (nr. active dihedrals) = %d\n", j);
435 return j;
436 }
437
438 static void histogramming(FILE *log, int nbin, ResidueType *rt,
439 int nf, int maxchi, real **dih,
440 int nlist, t_dlist dlist[],
441 const int index[],
442 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bOmega, gmx_bool bChi,
443 gmx_bool bNormalize, gmx_bool bSSHisto, const char *ssdump,
444 real bfac_max, const t_atoms *atoms,
445 gmx_bool bDo_jc, const char *fn,
446 const gmx_output_env_t *oenv)
447 {
448 /* also gets 3J couplings and order parameters S2 */
449 // Avoid warnings about narrowing conversions from double to real
450 #ifdef _MSC_VER
451 #pragma warning(disable: 4838)
452 #endif
453 t_karplus kkkphi[] = {
454 { "J_NHa1", 6.51, -1.76, 1.6, -M_PI/3, 0.0, 0.0 },
455 { "J_NHa2", 6.51, -1.76, 1.6, M_PI/3, 0.0, 0.0 },
456 { "J_HaC'", 4.0, 1.1, 0.1, 0.0, 0.0, 0.0 },
457 { "J_NHCb", 4.7, -1.5, -0.2, M_PI/3, 0.0, 0.0 },
458 { "J_Ci-1Hai", 4.5, -1.3, -1.2, 2*M_PI/3, 0.0, 0.0 }
459 };
460 t_karplus kkkpsi[] = {
461 { "J_HaN", -0.88, -0.61, -0.27, M_PI/3, 0.0, 0.0 }
462 };
463 t_karplus kkkchi1[] = {
464 { "JHaHb2", 9.5, -1.6, 1.8, -M_PI/3, 0, 0.0 },
465 { "JHaHb3", 9.5, -1.6, 1.8, 0, 0, 0.0 }
466 };
467 #ifdef _MSC_VER
468 #pragma warning(default: 4838)
469 #endif
470 #define NKKKPHI asize(kkkphi)
471 #define NKKKPSI asize(kkkpsi)
472 #define NKKKCHI asize(kkkchi1)
473 #define NJC (NKKKPHI+NKKKPSI+NKKKCHI)
474
475 FILE *fp, *ssfp[3] = {nullptr, nullptr, nullptr};
476 const char *sss[3] = { "sheet", "helix", "coil" };
477 real S2;
478 real *normhisto;
479 real **Jc, **Jcsig;
480 int ****his_aa_ss = nullptr;
481 int ***his_aa, *histmp;
482 int i, j, k, m, n, nn, Dih, nres, hindex, angle;
483 gmx_bool bBfac, bOccup;
484 char hisfile[256], hhisfile[256], title[256], *ss_str = nullptr;
485 char **leg;
486 const char *residue_name;
487 int rt_size;
488
489 rt_size = rt->numberOfEntries();
490 if (bSSHisto)
491 {
492 fp = gmx_ffopen(ssdump, "r");
493 if (1 != fscanf(fp, "%d", &nres))
494 {
495 gmx_fatal(FARGS, "Error reading from file %s", ssdump);
496 }
497
498 snew(ss_str, nres+1);
499 if (1 != fscanf(fp, "%s", ss_str))
500 {
501 gmx_fatal(FARGS, "Error reading from file %s", ssdump);
502 }
503
504 gmx_ffclose(fp);
505 /* Four dimensional array... Very cool */
506 snew(his_aa_ss, 3);
507 for (i = 0; (i < 3); i++)
508 {
509 snew(his_aa_ss[i], rt_size+1);
510 for (j = 0; (j <= rt_size); j++)
511 {
512 snew(his_aa_ss[i][j], edMax);
513 for (Dih = 0; (Dih < edMax); Dih++)
514 {
515 snew(his_aa_ss[i][j][Dih], nbin+1);
516 }
517 }
518 }
519 }
520 snew(his_aa, edMax);
521 for (Dih = 0; (Dih < edMax); Dih++)
522 {
523 snew(his_aa[Dih], rt_size+1);
524 for (i = 0; (i <= rt_size); i++)
525 {
526 snew(his_aa[Dih][i], nbin+1);
527 }
528 }
529 snew(histmp, nbin);
530
531 snew(Jc, nlist);
532 snew(Jcsig, nlist);
533 for (i = 0; (i < nlist); i++)
534 {
535 snew(Jc[i], NJC);
536 snew(Jcsig[i], NJC);
537 }
538
539 j = 0;
540 n = 0;
541 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
542 {
543 for (i = 0; (i < nlist); i++)
544 {
545 if (((Dih < edOmega) ) ||
546 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
547 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)))
548 {
549 make_histo(log, nf, dih[j], nbin, histmp, -M_PI, M_PI);
550
551 if (bSSHisto)
552 {
553 /* Assume there is only one structure, the first.
554 * Compute index in histogram.
555 */
556 /* Check the atoms to see whether their B-factors are low enough
557 * Check atoms to see their occupancy is 1.
558 */
559 bBfac = bOccup = TRUE;
560 for (nn = 0; (nn < 4); nn++, n++)
561 {
562 bBfac = bBfac && (atoms->pdbinfo[index[n]].bfac <= bfac_max);
563 bOccup = bOccup && (atoms->pdbinfo[index[n]].occup == 1);
564 }
565 if (bOccup && ((bfac_max <= 0) || bBfac))
566 {
567 hindex = static_cast<int>(((dih[j][0]+M_PI)*nbin)/(2*M_PI));
568 range_check(hindex, 0, nbin);
569
570 /* Assign dihedral to either of the structure determined
571 * histograms
572 */
573 switch (ss_str[dlist[i].resnr])
574 {
575 case 'E':
576 his_aa_ss[0][dlist[i].index][Dih][hindex]++;
577 break;
578 case 'H':
579 his_aa_ss[1][dlist[i].index][Dih][hindex]++;
580 break;
581 default:
582 his_aa_ss[2][dlist[i].index][Dih][hindex]++;
583 break;
584 }
585 }
586 else if (debug)
587 {
588 fprintf(debug, "Res. %d has imcomplete occupancy or bfacs > %g\n",
589 dlist[i].resnr, bfac_max);
590 }
591 }
592 else
593 {
594 n += 4;
595 }
596
597 switch (Dih)
598 {
599 case edPhi:
600 calc_distribution_props(nbin, histmp, -M_PI, NKKKPHI, kkkphi, &S2);
601
602 for (m = 0; (m < NKKKPHI); m++)
603 {
604 Jc[i][m] = kkkphi[m].Jc;
605 Jcsig[i][m] = kkkphi[m].Jcsig;
606 }
607 break;
608 case edPsi:
609 calc_distribution_props(nbin, histmp, -M_PI, NKKKPSI, kkkpsi, &S2);
610
611 for (m = 0; (m < NKKKPSI); m++)
612 {
613 Jc[i][NKKKPHI+m] = kkkpsi[m].Jc;
614 Jcsig[i][NKKKPHI+m] = kkkpsi[m].Jcsig;
615 }
616 break;
617 case edChi1:
618 calc_distribution_props(nbin, histmp, -M_PI, NKKKCHI, kkkchi1, &S2);
619 for (m = 0; (m < NKKKCHI); m++)
620 {
621 Jc[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jc;
622 Jcsig[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jcsig;
623 }
624 break;
625 default: /* covers edOmega and higher Chis than Chi1 */
626 calc_distribution_props(nbin, histmp, -M_PI, 0, nullptr, &S2);
627 break;
628 }
629 dlist[i].S2[Dih] = S2;
630
631 /* Sum distribution per amino acid type as well */
632 for (k = 0; (k < nbin); k++)
633 {
634 his_aa[Dih][dlist[i].index][k] += histmp[k];
635 histmp[k] = 0;
636 }
637 j++;
638 }
639 else /* dihed not defined */
640 {
641 dlist[i].S2[Dih] = 0.0;
642 }
643 }
644 }
645 sfree(histmp);
646
647 /* Print out Jcouplings */
648 fprintf(log, "\n *** J-Couplings from simulation (plus std. dev.) ***\n\n");
649 fprintf(log, "Residue ");
650 for (i = 0; (i < NKKKPHI); i++)
651 {
652 fprintf(log, "%7s SD", kkkphi[i].name);
653 }
654 for (i = 0; (i < NKKKPSI); i++)
655 {
656 fprintf(log, "%7s SD", kkkpsi[i].name);
657 }
658 for (i = 0; (i < NKKKCHI); i++)
659 {
660 fprintf(log, "%7s SD", kkkchi1[i].name);
661 }
662 fprintf(log, "\n");
663 for (i = 0; (i < NJC+1); i++)
664 {
665 fprintf(log, "------------");
666 }
667 fprintf(log, "\n");
668 for (i = 0; (i < nlist); i++)
669 {
670 fprintf(log, "%-10s", dlist[i].name);
671 for (j = 0; (j < NJC); j++)
672 {
673 fprintf(log, " %5.2f %4.2f", Jc[i][j], Jcsig[i][j]);
674 }
675 fprintf(log, "\n");
676 }
677 fprintf(log, "\n");
678
679 /* and to -jc file... */
680 if (bDo_jc)
681 {
682 fp = xvgropen(fn, "\\S3\\NJ-Couplings from Karplus Equation", "Residue",
683 "Coupling", oenv);
684 snew(leg, NJC);
685 for (i = 0; (i < NKKKPHI); i++)
686 {
687 leg[i] = gmx_strdup(kkkphi[i].name);
688 }
689 for (i = 0; (i < NKKKPSI); i++)
690 {
691 leg[i+NKKKPHI] = gmx_strdup(kkkpsi[i].name);
692 }
693 for (i = 0; (i < NKKKCHI); i++)
694 {
695 leg[i+NKKKPHI+NKKKPSI] = gmx_strdup(kkkchi1[i].name);
696 }
697 xvgr_legend(fp, NJC, leg, oenv);
698 fprintf(fp, "%5s ", "#Res.");
699 for (i = 0; (i < NJC); i++)
700 {
701 fprintf(fp, "%10s ", leg[i]);
702 }
703 fprintf(fp, "\n");
704 for (i = 0; (i < nlist); i++)
705 {
706 fprintf(fp, "%5d ", dlist[i].resnr);
707 for (j = 0; (j < NJC); j++)
708 {
709 fprintf(fp, " %8.3f", Jc[i][j]);
710 }
711 fprintf(fp, "\n");
712 }
713 xvgrclose(fp);
714 for (i = 0; (i < NJC); i++)
715 {
716 sfree(leg[i]);
717 }
718 }
719 /* finished -jc stuff */
720
721 snew(normhisto, nbin);
722 for (i = 0; (i < rt_size); i++)
723 {
724 for (Dih = 0; (Dih < edMax); Dih++)
725 {
726 /* First check whether something is in there */
727 for (j = 0; (j < nbin); j++)
728 {
729 if (his_aa[Dih][i][j] != 0)
730 {
731 break;
732 }
733 }
734 if ((j < nbin) &&
735 ((bPhi && (Dih == edPhi)) ||
736 (bPsi && (Dih == edPsi)) ||
737 (bOmega && (Dih == edOmega)) ||
738 (bChi && (Dih >= edChi1))))
739 {
740 if (bNormalize)
741 {
742 normalize_histo(nbin, his_aa[Dih][i], (360.0/nbin), normhisto);
743 }
744
745 residue_name = rt->nameFromResidueIndex(i).c_str();
746 switch (Dih)
747 {
748 case edPhi:
749 sprintf(hisfile, "histo-phi%s", residue_name);
750 sprintf(title, "\\xf\\f{} Distribution for %s", residue_name);
751 break;
752 case edPsi:
753 sprintf(hisfile, "histo-psi%s", residue_name);
754 sprintf(title, "\\xy\\f{} Distribution for %s", residue_name);
755 break;
756 case edOmega:
757 sprintf(hisfile, "histo-omega%s", residue_name);
758 sprintf(title, "\\xw\\f{} Distribution for %s", residue_name);
759 break;
760 default:
761 sprintf(hisfile, "histo-chi%d%s", Dih-NONCHI+1, residue_name);
762 sprintf(title, "\\xc\\f{}\\s%d\\N Distribution for %s",
763 Dih-NONCHI+1, residue_name);
764 }
765 std::strcpy(hhisfile, hisfile);
766 std::strcat(hhisfile, ".xvg");
767 fp = xvgropen(hhisfile, title, "Degrees", "", oenv);
768 if (output_env_get_print_xvgr_codes(oenv))
769 {
770 fprintf(fp, "@ with g0\n");
771 }
772 xvgr_world(fp, -180, 0, 180, 0.1, oenv);
773 if (output_env_get_print_xvgr_codes(oenv))
774 {
775 fprintf(fp, "# this effort to set graph size fails unless you run with -autoscale none or -autoscale y flags\n");
776 fprintf(fp, "@ xaxis tick on\n");
777 fprintf(fp, "@ xaxis tick major 90\n");
778 fprintf(fp, "@ xaxis tick minor 30\n");
779 fprintf(fp, "@ xaxis ticklabel prec 0\n");
780 fprintf(fp, "@ yaxis tick off\n");
781 fprintf(fp, "@ yaxis ticklabel off\n");
782 fprintf(fp, "@ type xy\n");
783 }
784 if (bSSHisto)
785 {
786 for (k = 0; (k < 3); k++)
787 {
788 std::string sshisfile = gmx::formatString("%s-%s.xvg", hisfile, sss[k]);
789 ssfp[k] = gmx_ffopen(sshisfile, "w");
790 }
791 }
792 for (j = 0; (j < nbin); j++)
793 {
794 angle = -180 + (360/nbin)*j;
795 if (bNormalize)
796 {
797 fprintf(fp, "%5d %10g\n", angle, normhisto[j]);
798 }
799 else
800 {
801 fprintf(fp, "%5d %10d\n", angle, his_aa[Dih][i][j]);
802 }
803 if (bSSHisto)
804 {
805 for (k = 0; (k < 3); k++)
806 {
807 fprintf(ssfp[k], "%5d %10d\n", angle,
808 his_aa_ss[k][i][Dih][j]);
809 }
810 }
811 }
812 fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
813 xvgrclose(fp);
814 if (bSSHisto)
815 {
816 for (k = 0; (k < 3); k++)
817 {
818 fprintf(ssfp[k], "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
819 gmx_ffclose(ssfp[k]);
820 }
821 }
822 }
823 }
824 }
825 sfree(normhisto);
826
827 if (bSSHisto)
828 {
829 /* Four dimensional array... Very cool */
830 for (i = 0; (i < 3); i++)
831 {
832 for (j = 0; (j <= rt_size); j++)
833 {
834 for (Dih = 0; (Dih < edMax); Dih++)
835 {
836 sfree(his_aa_ss[i][j][Dih]);
837 }
838 sfree(his_aa_ss[i][j]);
839 }
840 sfree(his_aa_ss[i]);
841 }
842 sfree(his_aa_ss);
843 sfree(ss_str);
844 }
845 }
846
847 static FILE *rama_file(const char *fn, const char *title, const char *xaxis,
848 const char *yaxis, const gmx_output_env_t *oenv)
849 {
850 FILE *fp;
851
852 fp = xvgropen(fn, title, xaxis, yaxis, oenv);
853 if (output_env_get_print_xvgr_codes(oenv))
854 {
855 fprintf(fp, "@ with g0\n");
856 }
857 xvgr_world(fp, -180, -180, 180, 180, oenv);
858 if (output_env_get_print_xvgr_codes(oenv))
859 {
860 fprintf(fp, "@ xaxis tick on\n");
861 fprintf(fp, "@ xaxis tick major 90\n");
862 fprintf(fp, "@ xaxis tick minor 30\n");
863 fprintf(fp, "@ xaxis ticklabel prec 0\n");
864 fprintf(fp, "@ yaxis tick on\n");
865 fprintf(fp, "@ yaxis tick major 90\n");
866 fprintf(fp, "@ yaxis tick minor 30\n");
867 fprintf(fp, "@ yaxis ticklabel prec 0\n");
868 fprintf(fp, "@ s0 type xy\n");
869 fprintf(fp, "@ s0 symbol 2\n");
870 fprintf(fp, "@ s0 symbol size 0.410000\n");
871 fprintf(fp, "@ s0 symbol fill 1\n");
872 fprintf(fp, "@ s0 symbol color 1\n");
873 fprintf(fp, "@ s0 symbol linewidth 1\n");
874 fprintf(fp, "@ s0 symbol linestyle 1\n");
875 fprintf(fp, "@ s0 symbol center false\n");
876 fprintf(fp, "@ s0 symbol char 0\n");
877 fprintf(fp, "@ s0 skip 0\n");
878 fprintf(fp, "@ s0 linestyle 0\n");
879 fprintf(fp, "@ s0 linewidth 1\n");
880 fprintf(fp, "@ type xy\n");
881 }
882 return fp;
883 }
884
885 static void do_rama(int nf, int nlist, t_dlist dlist[], real **dih,
886 gmx_bool bViol, gmx_bool bRamOmega, const gmx_output_env_t *oenv)
887 {
888 FILE *fp, *gp = nullptr;
889 gmx_bool bOm;
890 char fn[256];
891 int i, j, k, Xi1, Xi2, Phi, Psi, Om = 0, nlevels;
892 constexpr int NMAT = 120;
893 real **mat = nullptr, phi, psi, omega, axis[NMAT], lo, hi;
894 t_rgb rlo = { 1.0, 0.0, 0.0 };
895 t_rgb rmid = { 1.0, 1.0, 1.0 };
896 t_rgb rhi = { 0.0, 0.0, 1.0 };
897
898 for (i = 0; (i < nlist); i++)
899 {
900 if ((has_dihedral(edPhi, &(dlist[i]))) &&
901 (has_dihedral(edPsi, &(dlist[i]))))
902 {
903 sprintf(fn, "ramaPhiPsi%s.xvg", dlist[i].name);
904 fp = rama_file(fn, "Ramachandran Plot",
905 "\\8f\\4 (deg)", "\\8y\\4 (deg)", oenv);
906 bOm = bRamOmega && has_dihedral(edOmega, &(dlist[i]));
907 if (bOm)
908 {
909 Om = dlist[i].j0[edOmega];
910 snew(mat, NMAT);
911 for (j = 0; (j < NMAT); j++)
912 {
913 snew(mat[j], NMAT);
914 axis[j] = -180+gmx::exactDiv(360*j, NMAT);
915 }
916 }
917 if (bViol)
918 {
919 sprintf(fn, "violPhiPsi%s.xvg", dlist[i].name);
920 gp = gmx_ffopen(fn, "w");
921 }
922 Phi = dlist[i].j0[edPhi];
923 Psi = dlist[i].j0[edPsi];
924 for (j = 0; (j < nf); j++)
925 {
926 phi = RAD2DEG*dih[Phi][j];
927 psi = RAD2DEG*dih[Psi][j];
928 fprintf(fp, "%10g %10g\n", phi, psi);
929 if (bViol)
930 {
931 fprintf(gp, "%d\n", static_cast<int>(!bAllowed(dih[Phi][j], RAD2DEG*dih[Psi][j])) );
932 }
933 if (bOm)
934 {
935 omega = RAD2DEG*dih[Om][j];
936 mat[static_cast<int>(((phi*NMAT)/360)+gmx::exactDiv(NMAT, 2))][static_cast<int>(((psi*NMAT)/360)+gmx::exactDiv(NMAT, 2))]
937 += omega;
938 }
939 }
940 if (bViol)
941 {
942 gmx_ffclose(gp);
943 }
944 xvgrclose(fp);
945 if (bOm)
946 {
947 sprintf(fn, "ramomega%s.xpm", dlist[i].name);
948 fp = gmx_ffopen(fn, "w");
949 lo = hi = 0;
950 for (j = 0; (j < NMAT); j++)
951 {
952 for (k = 0; (k < NMAT); k++)
953 {
954 mat[j][k] /= nf;
955 lo = std::min(mat[j][k], lo);
956 hi = std::max(mat[j][k], hi);
957 }
958 }
959 /* Symmetrise */
960 if (std::abs(lo) > std::abs(hi))
961 {
962 hi = -lo;
963 }
964 else
965 {
966 lo = -hi;
967 }
968 /* Add 180 */
969 for (j = 0; (j < NMAT); j++)
970 {
971 for (k = 0; (k < NMAT); k++)
972 {
973 mat[j][k] += 180;
974 }
975 }
976 lo += 180;
977 hi += 180;
978 nlevels = 20;
979 write_xpm3(fp, 0, "Omega/Ramachandran Plot", "Deg", "Phi", "Psi",
980 NMAT, NMAT, axis, axis, mat, lo, 180.0, hi, rlo, rmid, rhi, &nlevels);
981 gmx_ffclose(fp);
982 for (j = 0; (j < NMAT); j++)
983 {
984 sfree(mat[j]);
985 }
986 sfree(mat);
987 }
988 }
989 if ((has_dihedral(edChi1, &(dlist[i]))) &&
990 (has_dihedral(edChi2, &(dlist[i]))))
991 {
992 sprintf(fn, "ramaX1X2%s.xvg", dlist[i].name);
993 fp = rama_file(fn, "\\8c\\4\\s1\\N-\\8c\\4\\s2\\N Ramachandran Plot",
994 "\\8c\\4\\s1\\N (deg)", "\\8c\\4\\s2\\N (deg)", oenv);
995 Xi1 = dlist[i].j0[edChi1];
996 Xi2 = dlist[i].j0[edChi2];
997 for (j = 0; (j < nf); j++)
998 {
999 fprintf(fp, "%10g %10g\n", RAD2DEG*dih[Xi1][j], RAD2DEG*dih[Xi2][j]);
1000 }
1001 xvgrclose(fp);
1002 }
1003 else
1004 {
1005 fprintf(stderr, "No chi1 & chi2 angle for %s\n", dlist[i].name);
1006 }
1007 }
1008 }
1009
1010
1011 static void print_transitions(const char *fn, int maxchi, int nlist,
1012 t_dlist dlist[], real dt,
1013 const gmx_output_env_t *oenv)
1014 {
1015 /* based on order_params below */
1016 FILE *fp;
1017 int i, Dih, Xi;
1018
1019 /* must correspond with enum in pp2shift.h:38 */
1020 char *leg[edMax];
1021 #define NLEG asize(leg)
1022
1023 leg[0] = gmx_strdup("Phi");
1024 leg[1] = gmx_strdup("Psi");
1025 leg[2] = gmx_strdup("Omega");
1026 leg[3] = gmx_strdup("Chi1");
1027 leg[4] = gmx_strdup("Chi2");
1028 leg[5] = gmx_strdup("Chi3");
1029 leg[6] = gmx_strdup("Chi4");
1030 leg[7] = gmx_strdup("Chi5");
1031 leg[8] = gmx_strdup("Chi6");
1032
1033 /* Print order parameters */
1034 fp = xvgropen(fn, "Dihedral Rotamer Transitions", "Residue", "Transitions/ns",
1035 oenv);
1036 xvgr_legend(fp, NONCHI+maxchi, leg, oenv);
1037
1038 fprintf(fp, "%5s ", "#Res.");
1039 fprintf(fp, "%10s %10s %10s ", leg[edPhi], leg[edPsi], leg[edOmega]);
1040 for (Xi = 0; Xi < maxchi; Xi++)
1041 {
1042 fprintf(fp, "%10s ", leg[NONCHI+Xi]);
1043 }
1044 fprintf(fp, "\n");
1045
1046 for (i = 0; (i < nlist); i++)
1047 {
1048 fprintf(fp, "%5d ", dlist[i].resnr);
1049 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1050 {
1051 fprintf(fp, "%10.3f ", dlist[i].ntr[Dih]/dt);
1052 }
1053 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
1054 fprintf(fp, "\n");
1055 }
1056 xvgrclose(fp);
1057 }
1058
1059 static void order_params(FILE *log,
1060 const char *fn, int maxchi, int nlist, t_dlist dlist[],
1061 const char *pdbfn, real bfac_init,
1062 t_atoms *atoms, const rvec x[], int ePBC, matrix box,
1063 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, const gmx_output_env_t *oenv)
1064 {
1065 FILE *fp;
1066 int nh[edMax];
1067 int i, Dih, Xi;
1068 real S2Max, S2Min;
1069
1070 /* except for S2Min/Max, must correspond with enum in pp2shift.h:38 */
1071 const char *const_leg[2+edMax] = {
1072 "S2Min", "S2Max", "Phi", "Psi", "Omega",
1073 "Chi1", "Chi2", "Chi3", "Chi4", "Chi5",
1074 "Chi6"
1075 };
1076 #define NLEG asize(leg)
1077
1078 char *leg[2+edMax];
1079
1080 for (i = 0; i < NLEG; i++)
1081 {
1082 leg[i] = gmx_strdup(const_leg[i]);
1083 }
1084
1085 /* Print order parameters */
1086 fp = xvgropen(fn, "Dihedral Order Parameters", "Residue", "S2", oenv);
1087 xvgr_legend(fp, 2+NONCHI+maxchi, const_leg, oenv);
1088
1089 for (Dih = 0; (Dih < edMax); Dih++)
1090 {
1091 nh[Dih] = 0;
1092 }
1093
1094 fprintf(fp, "%5s ", "#Res.");
1095 fprintf(fp, "%10s %10s ", leg[0], leg[1]);
1096 fprintf(fp, "%10s %10s %10s ", leg[2+edPhi], leg[2+edPsi], leg[2+edOmega]);
1097 for (Xi = 0; Xi < maxchi; Xi++)
1098 {
1099 fprintf(fp, "%10s ", leg[2+NONCHI+Xi]);
1100 }
1101 fprintf(fp, "\n");
1102
1103 for (i = 0; (i < nlist); i++)
1104 {
1105 S2Max = -10;
1106 S2Min = 10;
1107 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1108 {
1109 if (dlist[i].S2[Dih] != 0)
1110 {
1111 if (dlist[i].S2[Dih] > S2Max)
1112 {
1113 S2Max = dlist[i].S2[Dih];
1114 }
1115 if (dlist[i].S2[Dih] < S2Min)
1116 {
1117 S2Min = dlist[i].S2[Dih];
1118 }
1119 }
1120 if (dlist[i].S2[Dih] > 0.8)
1121 {
1122 nh[Dih]++;
1123 }
1124 }
1125 fprintf(fp, "%5d ", dlist[i].resnr);
1126 fprintf(fp, "%10.3f %10.3f ", S2Min, S2Max);
1127 for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
1128 {
1129 fprintf(fp, "%10.3f ", dlist[i].S2[Dih]);
1130 }
1131 fprintf(fp, "\n");
1132 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
1133 }
1134 xvgrclose(fp);
1135
1136 if (nullptr != pdbfn)
1137 {
1138 real x0, y0, z0;
1139
1140 atoms->havePdbInfo = TRUE;
1141
1142 if (nullptr == atoms->pdbinfo)
1143 {
1144 snew(atoms->pdbinfo, atoms->nr);
1145 }
1146 for (i = 0; (i < atoms->nr); i++)
1147 {
1148 atoms->pdbinfo[i].bfac = bfac_init;
1149 }
1150
1151 for (i = 0; (i < nlist); i++)
1152 {
1153 atoms->pdbinfo[dlist[i].atm.N].bfac = -dlist[i].S2[0]; /* Phi */
1154 atoms->pdbinfo[dlist[i].atm.H].bfac = -dlist[i].S2[0]; /* Phi */
1155 atoms->pdbinfo[dlist[i].atm.C].bfac = -dlist[i].S2[1]; /* Psi */
1156 atoms->pdbinfo[dlist[i].atm.O].bfac = -dlist[i].S2[1]; /* Psi */
1157 for (Xi = 0; (Xi < maxchi); Xi++) /* Chi's */
1158 {
1159 if (dlist[i].atm.Cn[Xi+3] != -1)
1160 {
1161 atoms->pdbinfo[dlist[i].atm.Cn[Xi+1]].bfac = -dlist[i].S2[NONCHI+Xi];
1162 }
1163 }
1164 }
1165
1166 fp = gmx_ffopen(pdbfn, "w");
1167 fprintf(fp, "REMARK generated by g_chi\n");
1168 fprintf(fp, "REMARK "
1169 "B-factor field contains negative of dihedral order parameters\n");
1170 write_pdbfile(fp, nullptr, atoms, x, ePBC, box, ' ', 0, nullptr);
1171 x0 = y0 = z0 = 1000.0;
1172 for (i = 0; (i < atoms->nr); i++)
1173 {
1174 x0 = std::min(x0, x[i][XX]);
1175 y0 = std::min(y0, x[i][YY]);
1176 z0 = std::min(z0, x[i][ZZ]);
1177 }
1178 x0 *= 10.0; /* nm -> angstrom */
1179 y0 *= 10.0; /* nm -> angstrom */
1180 z0 *= 10.0; /* nm -> angstrom */
1181 for (i = 0; (i < 10); i++)
1182 {
1183 gmx_fprintf_pdb_atomline(fp, epdbATOM, atoms->nr+1+i, "CA", ' ', "LEG", ' ', atoms->nres+1, ' ',
1184 x0, y0, z0+(1.2*i), 0.0, -0.1*i, "");
1185 }
1186 gmx_ffclose(fp);
1187 }
1188
1189 fprintf(log, "Dihedrals with S2 > 0.8\n");
1190 fprintf(log, "Dihedral: ");
1191 if (bPhi)
1192 {
1193 fprintf(log, " Phi ");
1194 }
1195 if (bPsi)
1196 {
1197 fprintf(log, " Psi ");
1198 }
1199 if (bChi)
1200 {
1201 for (Xi = 0; (Xi < maxchi); Xi++)
1202 {
1203 fprintf(log, " %s ", leg[2+NONCHI+Xi]);
1204 }
1205 }
1206 fprintf(log, "\nNumber: ");
1207 if (bPhi)
1208 {
1209 fprintf(log, "%4d ", nh[0]);
1210 }
1211 if (bPsi)
1212 {
1213 fprintf(log, "%4d ", nh[1]);
1214 }
1215 if (bChi)
1216 {
1217 for (Xi = 0; (Xi < maxchi); Xi++)
1218 {
1219 fprintf(log, "%4d ", nh[NONCHI+Xi]);
1220 }
1221 }
1222 fprintf(log, "\n");
1223
1224 for (i = 0; i < NLEG; i++)
1225 {
1226 sfree(leg[i]);
1227 }
1228
1229 }
1230
1231 int gmx_chi(int argc, char *argv[])
1232 {
1233 const char *desc[] = {
1234 "[THISMODULE] computes [GRK]phi[grk], [GRK]psi[grk], [GRK]omega[grk],",
1235 "and [GRK]chi[grk] dihedrals for all your",
1236 "amino acid backbone and sidechains.",
1237 "It can compute dihedral angle as a function of time, and as",
1238 "histogram distributions.",
1239 "The distributions [TT](histo-(dihedral)(RESIDUE).xvg[tt]) are cumulative over all residues of each type.[PAR]",
1240 "If option [TT]-corr[tt] is given, the program will",
1241 "calculate dihedral autocorrelation functions. The function used",
1242 "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",
1243 "rather than angles themselves, resolves the problem of periodicity.",
1244 "(Van der Spoel & Berendsen (1997), Biophys. J. 72, 2032-2041).",
1245 "Separate files for each dihedral of each residue",
1246 "[TT](corr(dihedral)(RESIDUE)(nresnr).xvg[tt]) are output, as well as a",
1247 "file containing the information for all residues (argument of [TT]-corr[tt]).[PAR]",
1248 "With option [TT]-all[tt], the angles themselves as a function of time for",
1249 "each residue are printed to separate files [TT](dihedral)(RESIDUE)(nresnr).xvg[tt].",
1250 "These can be in radians or degrees.[PAR]",
1251 "A log file (argument [TT]-g[tt]) is also written. This contains",
1252 "",
1253 " * information about the number of residues of each type.",
1254 " * The NMR ^3J coupling constants from the Karplus equation.",
1255 " * a table for each residue of the number of transitions between ",
1256 " rotamers per nanosecond, and the order parameter S^2 of each dihedral.",
1257 " * a table for each residue of the rotamer occupancy.",
1258 "",
1259 "All rotamers are taken as 3-fold, except for [GRK]omega[grk] and [GRK]chi[grk] dihedrals",
1260 "to planar groups (i.e. [GRK]chi[grk][SUB]2[sub] of aromatics, Asp and Asn; [GRK]chi[grk][SUB]3[sub] of Glu",
1261 "and Gln; and [GRK]chi[grk][SUB]4[sub] of Arg), which are 2-fold. \"rotamer 0\" means ",
1262 "that the dihedral was not in the core region of each rotamer. ",
1263 "The width of the core region can be set with [TT]-core_rotamer[tt][PAR]",
1264
1265 "The S^2 order parameters are also output to an [REF].xvg[ref] file",
1266 "(argument [TT]-o[tt] ) and optionally as a [REF].pdb[ref] file with",
1267 "the S^2 values as B-factor (argument [TT]-p[tt]). ",
1268 "The total number of rotamer transitions per timestep",
1269 "(argument [TT]-ot[tt]), the number of transitions per rotamer",
1270 "(argument [TT]-rt[tt]), and the ^3J couplings (argument [TT]-jc[tt]), ",
1271 "can also be written to [REF].xvg[ref] files. Note that the analysis",
1272 "of rotamer transitions assumes that the supplied trajectory frames",
1273 "are equally spaced in time.[PAR]",
1274
1275 "If [TT]-chi_prod[tt] is set (and [TT]-maxchi[tt] > 0), cumulative rotamers, e.g.",
1276 "1+9([GRK]chi[grk][SUB]1[sub]-1)+3([GRK]chi[grk][SUB]2[sub]-1)+([GRK]chi[grk][SUB]3[sub]-1) (if the residue has three 3-fold ",
1277 "dihedrals and [TT]-maxchi[tt] >= 3)",
1278 "are calculated. As before, if any dihedral is not in the core region,",
1279 "the rotamer is taken to be 0. The occupancies of these cumulative ",
1280 "rotamers (starting with rotamer 0) are written to the file",
1281 "that is the argument of [TT]-cp[tt], and if the [TT]-all[tt] flag",
1282 "is given, the rotamers as functions of time",
1283 "are written to [TT]chiproduct(RESIDUE)(nresnr).xvg[tt] ",
1284 "and their occupancies to [TT]histo-chiproduct(RESIDUE)(nresnr).xvg[tt].[PAR]",
1285
1286 "The option [TT]-r[tt] generates a contour plot of the average [GRK]omega[grk] angle",
1287 "as a function of the [GRK]phi[grk] and [GRK]psi[grk] angles, that is, in a Ramachandran plot",
1288 "the average [GRK]omega[grk] angle is plotted using color coding.",
1289
1290 };
1291
1292 const char *bugs[] = {
1293 "Produces MANY output files (up to about 4 times the number of residues in the protein, twice that if autocorrelation functions are calculated). Typically several hundred files are output.",
1294 "[GRK]phi[grk] and [GRK]psi[grk] dihedrals are calculated in a "
1295 "non-standard way, using H-N-CA-C for [GRK]phi[grk] instead of "
1296 "C(-)-N-CA-C, and N-CA-C-O for [GRK]psi[grk] instead of N-CA-C-N(+). "
1297 "This causes (usually small) discrepancies with the output of other "
1298 "tools like [gmx-rama].",
1299 "[TT]-r0[tt] option does not work properly",
1300 "Rotamers with multiplicity 2 are printed in [TT]chi.log[tt] as if they had multiplicity 3, with the 3rd (g(+)) always having probability 0"
1301 };
1302
1303 /* defaults */
1304 static int r0 = 1, ndeg = 1, maxchi = 2;
1305 static gmx_bool bAll = FALSE;
1306 static gmx_bool bPhi = FALSE, bPsi = FALSE, bOmega = FALSE;
1307 static real bfac_init = -1.0, bfac_max = 0;
1308 static const char *maxchistr[] = { nullptr, "0", "1", "2", "3", "4", "5", "6", nullptr };
1309 static gmx_bool bRama = FALSE, bShift = FALSE, bViol = FALSE, bRamOmega = FALSE;
1310 static gmx_bool bNormHisto = TRUE, bChiProduct = FALSE, bHChi = FALSE, bRAD = FALSE, bPBC = TRUE;
1311 static real core_frac = 0.5;
1312 t_pargs pa[] = {
1313 { "-r0", FALSE, etINT, {&r0},
1314 "starting residue" },
1315 { "-phi", FALSE, etBOOL, {&bPhi},
1316 "Output for [GRK]phi[grk] dihedral angles" },
1317 { "-psi", FALSE, etBOOL, {&bPsi},
1318 "Output for [GRK]psi[grk] dihedral angles" },
1319 { "-omega", FALSE, etBOOL, {&bOmega},
1320 "Output for [GRK]omega[grk] dihedrals (peptide bonds)" },
1321 { "-rama", FALSE, etBOOL, {&bRama},
1322 "Generate [GRK]phi[grk]/[GRK]psi[grk] and [GRK]chi[grk][SUB]1[sub]/[GRK]chi[grk][SUB]2[sub] Ramachandran plots" },
1323 { "-viol", FALSE, etBOOL, {&bViol},
1324 "Write a file that gives 0 or 1 for violated Ramachandran angles" },
1325 { "-periodic", FALSE, etBOOL, {&bPBC},
1326 "Print dihedral angles modulo 360 degrees" },
1327 { "-all", FALSE, etBOOL, {&bAll},
1328 "Output separate files for every dihedral." },
1329 { "-rad", FALSE, etBOOL, {&bRAD},
1330 "in angle vs time files, use radians rather than degrees."},
1331 { "-shift", FALSE, etBOOL, {&bShift},
1332 "Compute chemical shifts from [GRK]phi[grk]/[GRK]psi[grk] angles" },
1333 { "-binwidth", FALSE, etINT, {&ndeg},
1334 "bin width for histograms (degrees)" },
1335 { "-core_rotamer", FALSE, etREAL, {&core_frac},
1336 "only the central [TT]-core_rotamer[tt]\\*(360/multiplicity) belongs to each rotamer (the rest is assigned to rotamer 0)" },
1337 { "-maxchi", FALSE, etENUM, {maxchistr},
1338 "calculate first ndih [GRK]chi[grk] dihedrals" },
1339 { "-normhisto", FALSE, etBOOL, {&bNormHisto},
1340 "Normalize histograms" },
1341 { "-ramomega", FALSE, etBOOL, {&bRamOmega},
1342 "compute average omega as a function of [GRK]phi[grk]/[GRK]psi[grk] and plot it in an [REF].xpm[ref] plot" },
1343 { "-bfact", FALSE, etREAL, {&bfac_init},
1344 "B-factor value for [REF].pdb[ref] file for atoms with no calculated dihedral order parameter"},
1345 { "-chi_prod", FALSE, etBOOL, {&bChiProduct},
1346 "compute a single cumulative rotamer for each residue"},
1347 { "-HChi", FALSE, etBOOL, {&bHChi},
1348 "Include dihedrals to sidechain hydrogens"},
1349 { "-bmax", FALSE, etREAL, {&bfac_max},
1350 "Maximum B-factor on any of the atoms that make up a dihedral, for the dihedral angle to be considere in the statistics. Applies to database work where a number of X-Ray structures is analyzed. [TT]-bmax[tt] <= 0 means no limit." }
1351 };
1352
1353 FILE *log;
1354 int nlist, idum, nbin;
1355 rvec *x;
1356 int ePBC;
1357 matrix box;
1358 char grpname[256];
1359 t_dlist *dlist;
1360 gmx_bool bChi, bCorr, bSSHisto;
1361 gmx_bool bDo_rt, bDo_oh, bDo_ot, bDo_jc;
1362 real dt = 0, traj_t_ns;
1363 gmx_output_env_t *oenv;
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 ResidueType 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 return 0;
1596 }
The ultimate molecular dynamics simulation package