2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2017,2019,2020, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
41 #include "gromacs/commandline/pargs.h"
42 #include "gromacs/fileio/trxio.h"
43 #include "gromacs/fileio/xvgr.h"
44 #include "gromacs/gmxana/gmx_ana.h"
45 #include "gromacs/gmxana/gstat.h"
46 #include "gromacs/math/do_fit.h"
47 #include "gromacs/math/utilities.h"
48 #include "gromacs/math/vec.h"
49 #include "gromacs/pbcutil/pbc.h"
50 #include "gromacs/pbcutil/rmpbc.h"
51 #include "gromacs/topology/index.h"
52 #include "gromacs/topology/topology.h"
53 #include "gromacs/utility/arraysize.h"
54 #include "gromacs/utility/fatalerror.h"
55 #include "gromacs/utility/futil.h"
56 #include "gromacs/utility/smalloc.h"
58 int gmx_helixorient(int argc, char* argv[])
60 const char* desc[] = {
61 "[THISMODULE] calculates the coordinates and direction of the average",
62 "axis inside an alpha helix, and the direction/vectors of both the",
63 "C[GRK]alpha[grk] and (optionally) a sidechain atom relative to the axis.[PAR]",
64 "As input, you need to specify an index group with C[GRK]alpha[grk] atoms",
65 "corresponding to an [GRK]alpha[grk]-helix of continuous residues. Sidechain",
66 "directions require a second index group of the same size, containing",
67 "the heavy atom in each residue that should represent the sidechain.[PAR]",
68 "[BB]Note[bb] that this program does not do any fitting of structures.[PAR]",
69 "We need four C[GRK]alpha[grk] coordinates to define the local direction of the helix",
71 "The tilt/rotation is calculated from Euler rotations, where we define",
72 "the helix axis as the local [IT]x[it]-axis, the residues/C[GRK]alpha[grk] vector as ",
74 "[IT]z[it]-axis from their cross product. We use the Euler Y-Z-X rotation, meaning",
75 "we first tilt the helix axis (1) around and (2) orthogonal to the residues",
76 "vector, and finally apply the (3) rotation around it. For debugging or other",
77 "purposes, we also write out the actual Euler rotation angles as [TT]theta[1-3].xvg[tt]"
80 t_topology* top = nullptr;
86 real theta1, theta2, theta3;
102 rvec* residuehelixaxis;
105 rvec* sidechainvector;
107 rvec* residuehelixaxis_t0;
108 rvec* residuevector_t0;
110 rvec* residuehelixaxis_tlast;
111 rvec* residuevector_tlast;
113 rvec refaxes[3], newaxes[3];
115 rvec rot_refaxes[3], rot_newaxes[3];
119 real *twist, *residuetwist;
120 real *radius, *residueradius;
121 real *rise, *residuerise;
122 real* residuebending;
129 FILE * fpaxis, *fpcenter, *fptilt, *fprotation;
130 FILE * fpradius, *fprise, *fptwist;
131 FILE * fptheta1, *fptheta2, *fptheta3;
135 gmx_output_env_t* oenv;
136 gmx_rmpbc_t gpbc = nullptr;
138 static gmx_bool bSC = FALSE;
139 static gmx_bool bIncremental = FALSE;
141 static t_pargs pa[] = {
146 "Calculate sidechain directions relative to helix axis too." },
151 "Calculate incremental rather than total rotation/tilt." },
153 #define NPA asize(pa)
156 { efTPR, nullptr, nullptr, ffREAD }, { efTRX, "-f", nullptr, ffREAD },
157 { efNDX, nullptr, nullptr, ffOPTRD }, { efDAT, "-oaxis", "helixaxis", ffWRITE },
158 { efDAT, "-ocenter", "center", ffWRITE }, { efXVG, "-orise", "rise", ffWRITE },
159 { efXVG, "-oradius", "radius", ffWRITE }, { efXVG, "-otwist", "twist", ffWRITE },
160 { efXVG, "-obending", "bending", ffWRITE }, { efXVG, "-otilt", "tilt", ffWRITE },
161 { efXVG, "-orot", "rotation", ffWRITE }
163 #define NFILE asize(fnm)
165 if (!parse_common_args(
166 &argc, argv, PCA_CAN_TIME, NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
171 top = read_top(ftp2fn(efTPR, NFILE, fnm), &pbcType);
173 for (i = 0; i < 3; i++)
178 /* read index files */
179 printf("Select a group of Calpha atoms corresponding to a single continuous helix:\n");
180 get_index(&(top->atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &iCA, &ind_CA, &gn_CA);
182 snew(x_SC, iCA); /* sic! */
187 snew(diff13, iCA - 3);
188 snew(diff24, iCA - 3);
189 snew(helixaxis, iCA - 3);
191 snew(residuetwist, iCA);
193 snew(residueradius, iCA);
195 snew(residuerise, iCA);
196 snew(residueorigin, iCA);
197 snew(residuehelixaxis, iCA);
198 snew(residuevector, iCA);
199 snew(sidechainvector, iCA);
200 snew(residuebending, iCA);
201 snew(residuehelixaxis_t0, iCA);
202 snew(residuevector_t0, iCA);
204 snew(residuehelixaxis_tlast, iCA);
205 snew(residuevector_tlast, iCA);
206 snew(axis3_tlast, iCA);
211 printf("Select a group of atoms defining the sidechain direction (1/residue):\n");
212 get_index(&(top->atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &iSC, &ind_SC, &gn_SC);
215 gmx_fatal(FARGS, "Number of sidechain atoms (%d) != number of CA atoms (%d)", iSC, iCA);
219 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
221 fpaxis = gmx_ffopen(opt2fn("-oaxis", NFILE, fnm), "w");
222 fpcenter = gmx_ffopen(opt2fn("-ocenter", NFILE, fnm), "w");
223 fprise = gmx_ffopen(opt2fn("-orise", NFILE, fnm), "w");
224 fpradius = gmx_ffopen(opt2fn("-oradius", NFILE, fnm), "w");
225 fptwist = gmx_ffopen(opt2fn("-otwist", NFILE, fnm), "w");
226 fpbending = gmx_ffopen(opt2fn("-obending", NFILE, fnm), "w");
228 fptheta1 = gmx_ffopen("theta1.xvg", "w");
229 fptheta2 = gmx_ffopen("theta2.xvg", "w");
230 fptheta3 = gmx_ffopen("theta3.xvg", "w");
234 fptilt = xvgropen(opt2fn("-otilt", NFILE, fnm),
235 "Incremental local helix tilt",
239 fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
240 "Incremental local helix rotation",
242 "Rotation (degrees)",
247 fptilt = xvgropen(opt2fn("-otilt", NFILE, fnm),
248 "Cumulative local helix tilt",
252 fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
253 "Cumulative local helix rotation",
255 "Rotation (degrees)",
259 clear_rvecs(3, unitaxes);
264 gpbc = gmx_rmpbc_init(&top->idef, pbcType, natoms);
268 /* initialisation for correct distance calculations */
269 set_pbc(&pbc, pbcType, box);
270 /* make molecules whole again */
271 gmx_rmpbc(gpbc, natoms, box, x);
273 /* copy coords to our smaller arrays */
274 for (i = 0; i < iCA; i++)
276 copy_rvec(x[ind_CA[i]], x_CA[i]);
279 copy_rvec(x[ind_SC[i]], x_SC[i]);
283 for (i = 0; i < iCA - 3; i++)
285 rvec_sub(x_CA[i + 1], x_CA[i], r12[i]);
286 rvec_sub(x_CA[i + 2], x_CA[i + 1], r23[i]);
287 rvec_sub(x_CA[i + 3], x_CA[i + 2], r34[i]);
288 rvec_sub(r12[i], r23[i], diff13[i]);
289 rvec_sub(r23[i], r34[i], diff24[i]);
290 /* calculate helix axis */
291 cprod(diff13[i], diff24[i], helixaxis[i]);
292 svmul(1.0 / norm(helixaxis[i]), helixaxis[i], helixaxis[i]);
294 tmp = cos_angle(diff13[i], diff24[i]);
295 twist[i] = 180.0 / M_PI * std::acos(tmp);
296 radius[i] = std::sqrt(norm(diff13[i]) * norm(diff24[i])) / (2.0 * (1.0 - tmp));
297 rise[i] = std::abs(iprod(r23[i], helixaxis[i]));
299 svmul(radius[i] / norm(diff13[i]), diff13[i], v1);
300 svmul(radius[i] / norm(diff24[i]), diff24[i], v2);
302 rvec_sub(x_CA[i + 1], v1, residueorigin[i + 1]);
303 rvec_sub(x_CA[i + 2], v2, residueorigin[i + 2]);
305 residueradius[0] = residuetwist[0] = residuerise[0] = 0;
307 residueradius[1] = radius[0];
308 residuetwist[1] = twist[0];
309 residuerise[1] = rise[0];
311 residuebending[0] = residuebending[1] = 0;
312 for (i = 2; i < iCA - 2; i++)
314 residueradius[i] = 0.5 * (radius[i - 2] + radius[i - 1]);
315 residuetwist[i] = 0.5 * (twist[i - 2] + twist[i - 1]);
316 residuerise[i] = 0.5 * (rise[i - 2] + rise[i - 1]);
317 residuebending[i] = 180.0 / M_PI * std::acos(cos_angle(helixaxis[i - 2], helixaxis[i - 1]));
319 residueradius[iCA - 2] = radius[iCA - 4];
320 residuetwist[iCA - 2] = twist[iCA - 4];
321 residuerise[iCA - 2] = rise[iCA - 4];
322 residueradius[iCA - 1] = residuetwist[iCA - 1] = residuerise[iCA - 1] = 0;
323 residuebending[iCA - 2] = residuebending[iCA - 1] = 0;
325 clear_rvec(residueorigin[0]);
326 clear_rvec(residueorigin[iCA - 1]);
328 /* average helix axes to define them on the residues.
329 * Just extrapolate second first/list atom.
331 copy_rvec(helixaxis[0], residuehelixaxis[0]);
332 copy_rvec(helixaxis[0], residuehelixaxis[1]);
334 for (i = 2; i < iCA - 2; i++)
336 rvec_add(helixaxis[i - 2], helixaxis[i - 1], residuehelixaxis[i]);
337 svmul(0.5, residuehelixaxis[i], residuehelixaxis[i]);
339 copy_rvec(helixaxis[iCA - 4], residuehelixaxis[iCA - 2]);
340 copy_rvec(helixaxis[iCA - 4], residuehelixaxis[iCA - 1]);
342 /* Normalize the axis */
343 for (i = 0; i < iCA; i++)
345 svmul(1.0 / norm(residuehelixaxis[i]), residuehelixaxis[i], residuehelixaxis[i]);
348 /* calculate vector from origin to residue CA */
349 fprintf(fpaxis, "%15.12g ", t);
350 fprintf(fpcenter, "%15.12g ", t);
351 fprintf(fprise, "%15.12g ", t);
352 fprintf(fpradius, "%15.12g ", t);
353 fprintf(fptwist, "%15.12g ", t);
354 fprintf(fpbending, "%15.12g ", t);
356 for (i = 0; i < iCA; i++)
358 if (i == 0 || i == iCA - 1)
360 fprintf(fpaxis, "%15.12g %15.12g %15.12g ", 0.0, 0.0, 0.0);
361 fprintf(fpcenter, "%15.12g %15.12g %15.12g ", 0.0, 0.0, 0.0);
362 fprintf(fprise, "%15.12g ", 0.0);
363 fprintf(fpradius, "%15.12g ", 0.0);
364 fprintf(fptwist, "%15.12g ", 0.0);
365 fprintf(fpbending, "%15.12g ", 0.0);
369 rvec_sub(bSC ? x_SC[i] : x_CA[i], residueorigin[i], residuevector[i]);
370 svmul(1.0 / norm(residuevector[i]), residuevector[i], residuevector[i]);
371 cprod(residuehelixaxis[i], residuevector[i], axis3[i]);
373 "%15.12g %15.12g %15.12g ",
374 residuehelixaxis[i][0],
375 residuehelixaxis[i][1],
376 residuehelixaxis[i][2]);
378 "%15.12g %15.12g %15.12g ",
381 residueorigin[i][2]);
383 fprintf(fprise, "%15.12g ", residuerise[i]);
384 fprintf(fpradius, "%15.12g ", residueradius[i]);
385 fprintf(fptwist, "%15.12g ", residuetwist[i]);
386 fprintf(fpbending, "%15.12g ", residuebending[i]);
389 fprintf(fprise, "\n");
390 fprintf(fpradius, "\n");
391 fprintf(fpaxis, "\n");
392 fprintf(fpcenter, "\n");
393 fprintf(fptwist, "\n");
394 fprintf(fpbending, "\n");
398 for (i = 0; i < iCA; i++)
400 copy_rvec(residuehelixaxis[i], residuehelixaxis_t0[i]);
401 copy_rvec(residuevector[i], residuevector_t0[i]);
402 copy_rvec(axis3[i], axis3_t0[i]);
407 fprintf(fptilt, "%15.12g ", t);
408 fprintf(fprotation, "%15.12g ", t);
409 fprintf(fptheta1, "%15.12g ", t);
410 fprintf(fptheta2, "%15.12g ", t);
411 fprintf(fptheta3, "%15.12g ", t);
413 for (i = 0; i < iCA; i++)
415 if (i == 0 || i == iCA - 1)
423 /* Total rotation & tilt */
424 copy_rvec(residuehelixaxis_t0[i], refaxes[0]);
425 copy_rvec(residuevector_t0[i], refaxes[1]);
426 copy_rvec(axis3_t0[i], refaxes[2]);
430 /* Rotation/tilt since last step */
431 copy_rvec(residuehelixaxis_tlast[i], refaxes[0]);
432 copy_rvec(residuevector_tlast[i], refaxes[1]);
433 copy_rvec(axis3_tlast[i], refaxes[2]);
435 copy_rvec(residuehelixaxis[i], newaxes[0]);
436 copy_rvec(residuevector[i], newaxes[1]);
437 copy_rvec(axis3[i], newaxes[2]);
439 /* rotate reference frame onto unit axes */
440 calc_fit_R(3, 3, weight, unitaxes, refaxes, A);
441 for (j = 0; j < 3; j++)
443 mvmul(A, refaxes[j], rot_refaxes[j]);
444 mvmul(A, newaxes[j], rot_newaxes[j]);
447 /* Determine local rotation matrix A */
448 calc_fit_R(3, 3, weight, rot_newaxes, rot_refaxes, A);
449 /* Calculate euler angles, from rotation order y-z-x, where
450 * x is helixaxis, y residuevector, and z axis3.
452 * A contains rotation column vectors.
455 theta1 = 180.0 / M_PI * std::atan2(A[0][2], A[0][0]);
456 theta2 = 180.0 / M_PI * std::asin(-A[0][1]);
457 theta3 = 180.0 / M_PI * std::atan2(A[2][1], A[1][1]);
459 tilt = std::sqrt(theta1 * theta1 + theta2 * theta2);
461 fprintf(fptheta1, "%15.12g ", theta1);
462 fprintf(fptheta2, "%15.12g ", theta2);
463 fprintf(fptheta3, "%15.12g ", theta3);
465 fprintf(fptilt, "%15.12g ", tilt);
466 fprintf(fprotation, "%15.12g ", rotation);
468 fprintf(fptilt, "\n");
469 fprintf(fprotation, "\n");
470 fprintf(fptheta1, "\n");
471 fprintf(fptheta2, "\n");
472 fprintf(fptheta3, "\n");
475 for (i = 0; i < iCA; i++)
477 copy_rvec(residuehelixaxis[i], residuehelixaxis_tlast[i]);
478 copy_rvec(residuevector[i], residuevector_tlast[i]);
479 copy_rvec(axis3[i], axis3_tlast[i]);
483 } while (read_next_x(oenv, status, &t, x, box));
485 gmx_rmpbc_done(gpbc);
488 gmx_ffclose(fpcenter);
490 xvgrclose(fprotation);
492 gmx_ffclose(fpradius);
493 gmx_ffclose(fptwist);
494 gmx_ffclose(fpbending);
495 gmx_ffclose(fptheta1);
496 gmx_ffclose(fptheta2);
497 gmx_ffclose(fptheta3);