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, 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.
45 #include "gromacs/math/vec.h"
52 #include "gromacs/commandline/pargs.h"
53 #include "gromacs/fileio/trxio.h"
54 #include "gromacs/fileio/xvgr.h"
55 #include "gromacs/math/do_fit.h"
56 #include "gromacs/pbcutil/rmpbc.h"
57 #include "gromacs/utility/fatalerror.h"
58 #include "gromacs/utility/futil.h"
59 #include "gromacs/utility/smalloc.h"
61 int gmx_helixorient(int argc, char *argv[])
63 const char *desc[] = {
64 "[THISMODULE] calculates the coordinates and direction of the average",
65 "axis inside an alpha helix, and the direction/vectors of both the",
66 "C[GRK]alpha[grk] and (optionally) a sidechain atom relative to the axis.[PAR]",
67 "As input, you need to specify an index group with C[GRK]alpha[grk] atoms",
68 "corresponding to an [GRK]alpha[grk]-helix of continuous residues. Sidechain",
69 "directions require a second index group of the same size, containing",
70 "the heavy atom in each residue that should represent the sidechain.[PAR]",
71 "[BB]Note[bb] that this program does not do any fitting of structures.[PAR]",
72 "We need four C[GRK]alpha[grk] coordinates to define the local direction of the helix",
74 "The tilt/rotation is calculated from Euler rotations, where we define",
75 "the helix axis as the local [IT]x[it]-axis, the residues/C[GRK]alpha[grk] vector as [IT]y[it], and the",
76 "[IT]z[it]-axis from their cross product. We use the Euler Y-Z-X rotation, meaning",
77 "we first tilt the helix axis (1) around and (2) orthogonal to the residues",
78 "vector, and finally apply the (3) rotation around it. For debugging or other",
79 "purposes, we also write out the actual Euler rotation angles as [TT]theta[1-3].xvg[tt]"
82 t_topology *top = NULL;
88 real theta1, theta2, theta3;
90 int d, i, j, teller = 0;
97 rvec v1, v2, p1, p2, vtmp, vproj;
105 rvec *residuehelixaxis;
108 rvec *sidechainvector;
112 rvec *residuehelixaxis_t0;
113 rvec *residuevector_t0;
115 rvec *residuehelixaxis_tlast;
116 rvec *residuevector_tlast;
118 rvec refaxes[3], newaxes[3];
120 rvec rot_refaxes[3], rot_newaxes[3];
124 real *twist, *residuetwist;
125 real *radius, *residueradius;
126 real *rise, *residuerise;
127 real *residuebending;
134 FILE *fpaxis, *fpcenter, *fptilt, *fprotation;
135 FILE *fpradius, *fprise, *fptwist;
136 FILE *fptheta1, *fptheta2, *fptheta3;
141 gmx_rmpbc_t gpbc = NULL;
143 static gmx_bool bSC = FALSE;
144 static gmx_bool bIncremental = FALSE;
146 static t_pargs pa[] = {
147 { "-sidechain", FALSE, etBOOL, {&bSC},
148 "Calculate sidechain directions relative to helix axis too." },
149 { "-incremental", FALSE, etBOOL, {&bIncremental},
150 "Calculate incremental rather than total rotation/tilt." },
152 #define NPA asize(pa)
155 { efTPX, NULL, NULL, ffREAD },
156 { efTRX, "-f", NULL, ffREAD },
157 { efNDX, NULL, NULL, ffOPTRD },
158 { efDAT, "-oaxis", "helixaxis", ffWRITE },
159 { efDAT, "-ocenter", "center", ffWRITE },
160 { efXVG, "-orise", "rise", ffWRITE },
161 { efXVG, "-oradius", "radius", ffWRITE },
162 { efXVG, "-otwist", "twist", ffWRITE },
163 { efXVG, "-obending", "bending", ffWRITE },
164 { efXVG, "-otilt", "tilt", ffWRITE },
165 { efXVG, "-orot", "rotation", ffWRITE }
167 #define NFILE asize(fnm)
169 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_BE_NICE,
170 NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
175 top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC);
177 for (i = 0; i < 3; i++)
182 /* read index files */
183 printf("Select a group of Calpha atoms corresponding to a single continuous helix:\n");
184 get_index(&(top->atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &iCA, &ind_CA, &gn_CA);
186 snew(x_SC, iCA); /* sic! */
193 snew(helixaxis, iCA-3);
195 snew(residuetwist, iCA);
197 snew(residueradius, iCA);
199 snew(residuerise, iCA);
200 snew(residueorigin, iCA);
201 snew(residuehelixaxis, iCA);
202 snew(residuevector, iCA);
203 snew(sidechainvector, iCA);
204 snew(residuebending, iCA);
205 snew(residuehelixaxis_t0, iCA);
206 snew(residuevector_t0, iCA);
208 snew(residuehelixaxis_tlast, iCA);
209 snew(residuevector_tlast, iCA);
210 snew(axis3_tlast, iCA);
215 printf("Select a group of atoms defining the sidechain direction (1/residue):\n");
216 get_index(&(top->atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &iSC, &ind_SC, &gn_SC);
219 gmx_fatal(FARGS, "Number of sidechain atoms (%d) != number of CA atoms (%d)", iSC, iCA);
224 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
226 fpaxis = gmx_ffopen(opt2fn("-oaxis", NFILE, fnm), "w");
227 fpcenter = gmx_ffopen(opt2fn("-ocenter", NFILE, fnm), "w");
228 fprise = gmx_ffopen(opt2fn("-orise", NFILE, fnm), "w");
229 fpradius = gmx_ffopen(opt2fn("-oradius", NFILE, fnm), "w");
230 fptwist = gmx_ffopen(opt2fn("-otwist", NFILE, fnm), "w");
231 fpbending = gmx_ffopen(opt2fn("-obending", NFILE, fnm), "w");
233 fptheta1 = gmx_ffopen("theta1.xvg", "w");
234 fptheta2 = gmx_ffopen("theta2.xvg", "w");
235 fptheta3 = gmx_ffopen("theta3.xvg", "w");
239 fptilt = xvgropen(opt2fn("-otilt", NFILE, fnm),
240 "Incremental local helix tilt", "Time(ps)", "Tilt (degrees)",
242 fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
243 "Incremental local helix rotation", "Time(ps)",
244 "Rotation (degrees)", oenv);
248 fptilt = xvgropen(opt2fn("-otilt", NFILE, fnm),
249 "Cumulative local helix tilt", "Time(ps)", "Tilt (degrees)", oenv);
250 fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
251 "Cumulative local helix rotation", "Time(ps)",
252 "Rotation (degrees)", oenv);
255 clear_rvecs(3, unitaxes);
260 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
264 /* initialisation for correct distance calculations */
265 set_pbc(&pbc, ePBC, box);
266 /* make molecules whole again */
267 gmx_rmpbc(gpbc, natoms, box, x);
269 /* copy coords to our smaller arrays */
270 for (i = 0; i < iCA; i++)
272 copy_rvec(x[ind_CA[i]], x_CA[i]);
275 copy_rvec(x[ind_SC[i]], x_SC[i]);
279 for (i = 0; i < iCA-3; i++)
281 rvec_sub(x_CA[i+1], x_CA[i], r12[i]);
282 rvec_sub(x_CA[i+2], x_CA[i+1], r23[i]);
283 rvec_sub(x_CA[i+3], x_CA[i+2], r34[i]);
284 rvec_sub(r12[i], r23[i], diff13[i]);
285 rvec_sub(r23[i], r34[i], diff24[i]);
286 /* calculate helix axis */
287 cprod(diff13[i], diff24[i], helixaxis[i]);
288 svmul(1.0/norm(helixaxis[i]), helixaxis[i], helixaxis[i]);
290 tmp = cos_angle(diff13[i], diff24[i]);
291 twist[i] = 180.0/M_PI * acos( tmp );
292 radius[i] = sqrt( norm(diff13[i])*norm(diff24[i]) ) / (2.0* (1.0-tmp) );
293 rise[i] = fabs(iprod(r23[i], helixaxis[i]));
295 svmul(radius[i]/norm(diff13[i]), diff13[i], v1);
296 svmul(radius[i]/norm(diff24[i]), diff24[i], v2);
298 rvec_sub(x_CA[i+1], v1, residueorigin[i+1]);
299 rvec_sub(x_CA[i+2], v2, residueorigin[i+2]);
301 residueradius[0] = residuetwist[0] = residuerise[0] = 0;
303 residueradius[1] = radius[0];
304 residuetwist[1] = twist[0];
305 residuerise[1] = rise[0];
307 residuebending[0] = residuebending[1] = 0;
308 for (i = 2; i < iCA-2; i++)
310 residueradius[i] = 0.5*(radius[i-2]+radius[i-1]);
311 residuetwist[i] = 0.5*(twist[i-2]+twist[i-1]);
312 residuerise[i] = 0.5*(rise[i-2]+rise[i-1]);
313 residuebending[i] = 180.0/M_PI*acos( cos_angle(helixaxis[i-2], helixaxis[i-1]) );
315 residueradius[iCA-2] = radius[iCA-4];
316 residuetwist[iCA-2] = twist[iCA-4];
317 residuerise[iCA-2] = rise[iCA-4];
318 residueradius[iCA-1] = residuetwist[iCA-1] = residuerise[iCA-1] = 0;
319 residuebending[iCA-2] = residuebending[iCA-1] = 0;
321 clear_rvec(residueorigin[0]);
322 clear_rvec(residueorigin[iCA-1]);
324 /* average helix axes to define them on the residues.
325 * Just extrapolate second first/list atom.
327 copy_rvec(helixaxis[0], residuehelixaxis[0]);
328 copy_rvec(helixaxis[0], residuehelixaxis[1]);
330 for (i = 2; i < iCA-2; i++)
332 rvec_add(helixaxis[i-2], helixaxis[i-1], residuehelixaxis[i]);
333 svmul(0.5, residuehelixaxis[i], residuehelixaxis[i]);
335 copy_rvec(helixaxis[iCA-4], residuehelixaxis[iCA-2]);
336 copy_rvec(helixaxis[iCA-4], residuehelixaxis[iCA-1]);
338 /* Normalize the axis */
339 for (i = 0; i < iCA; i++)
341 svmul(1.0/norm(residuehelixaxis[i]), residuehelixaxis[i], residuehelixaxis[i]);
344 /* calculate vector from origin to residue CA */
345 fprintf(fpaxis, "%15.12g ", t);
346 fprintf(fpcenter, "%15.12g ", t);
347 fprintf(fprise, "%15.12g ", t);
348 fprintf(fpradius, "%15.12g ", t);
349 fprintf(fptwist, "%15.12g ", t);
350 fprintf(fpbending, "%15.12g ", t);
352 for (i = 0; i < iCA; i++)
354 if (i == 0 || i == iCA-1)
356 fprintf(fpaxis, "%15.12g %15.12g %15.12g ", 0.0, 0.0, 0.0);
357 fprintf(fpcenter, "%15.12g %15.12g %15.12g ", 0.0, 0.0, 0.0);
358 fprintf(fprise, "%15.12g ", 0.0);
359 fprintf(fpradius, "%15.12g ", 0.0);
360 fprintf(fptwist, "%15.12g ", 0.0);
361 fprintf(fpbending, "%15.12g ", 0.0);
365 rvec_sub( bSC ? x_SC[i] : x_CA[i], residueorigin[i], residuevector[i]);
366 svmul(1.0/norm(residuevector[i]), residuevector[i], residuevector[i]);
367 cprod(residuehelixaxis[i], residuevector[i], axis3[i]);
368 fprintf(fpaxis, "%15.12g %15.12g %15.12g ", residuehelixaxis[i][0], residuehelixaxis[i][1], residuehelixaxis[i][2]);
369 fprintf(fpcenter, "%15.12g %15.12g %15.12g ", residueorigin[i][0], residueorigin[i][1], residueorigin[i][2]);
371 fprintf(fprise, "%15.12g ", residuerise[i]);
372 fprintf(fpradius, "%15.12g ", residueradius[i]);
373 fprintf(fptwist, "%15.12g ", residuetwist[i]);
374 fprintf(fpbending, "%15.12g ", residuebending[i]);
376 /* angle with local vector? */
378 printf("res[%2d]: axis: %g %g %g origin: %g %g %g vector: %g %g %g angle: %g\n",i,
379 residuehelixaxis[i][0],
380 residuehelixaxis[i][1],
381 residuehelixaxis[i][2],
388 180.0/M_PI*acos( cos_angle(residuevector[i],residuehelixaxis[i]) ));
390 /* fprintf(fp,"%15.12g %15.12g %15.12g %15.12g %15.12g %15.12g\n",
391 residuehelixaxis[i][0],
392 residuehelixaxis[i][1],
393 residuehelixaxis[i][2],
396 residuevector[i][2]);
400 fprintf(fprise, "\n");
401 fprintf(fpradius, "\n");
402 fprintf(fpaxis, "\n");
403 fprintf(fpcenter, "\n");
404 fprintf(fptwist, "\n");
405 fprintf(fpbending, "\n");
409 for (i = 0; i < iCA; i++)
411 copy_rvec(residuehelixaxis[i], residuehelixaxis_t0[i]);
412 copy_rvec(residuevector[i], residuevector_t0[i]);
413 copy_rvec(axis3[i], axis3_t0[i]);
418 fprintf(fptilt, "%15.12g ", t);
419 fprintf(fprotation, "%15.12g ", t);
420 fprintf(fptheta1, "%15.12g ", t);
421 fprintf(fptheta2, "%15.12g ", t);
422 fprintf(fptheta3, "%15.12g ", t);
424 for (i = 0; i < iCA; i++)
426 if (i == 0 || i == iCA-1)
434 /* Total rotation & tilt */
435 copy_rvec(residuehelixaxis_t0[i], refaxes[0]);
436 copy_rvec(residuevector_t0[i], refaxes[1]);
437 copy_rvec(axis3_t0[i], refaxes[2]);
441 /* Rotation/tilt since last step */
442 copy_rvec(residuehelixaxis_tlast[i], refaxes[0]);
443 copy_rvec(residuevector_tlast[i], refaxes[1]);
444 copy_rvec(axis3_tlast[i], refaxes[2]);
446 copy_rvec(residuehelixaxis[i], newaxes[0]);
447 copy_rvec(residuevector[i], newaxes[1]);
448 copy_rvec(axis3[i], newaxes[2]);
451 printf("frame %d, i=%d:\n old: %g %g %g , %g %g %g , %g %g %g\n new: %g %g %g , %g %g %g , %g %g %g\n",
453 refaxes[0][0],refaxes[0][1],refaxes[0][2],
454 refaxes[1][0],refaxes[1][1],refaxes[1][2],
455 refaxes[2][0],refaxes[2][1],refaxes[2][2],
456 newaxes[0][0],newaxes[0][1],newaxes[0][2],
457 newaxes[1][0],newaxes[1][1],newaxes[1][2],
458 newaxes[2][0],newaxes[2][1],newaxes[2][2]);
461 /* rotate reference frame onto unit axes */
462 calc_fit_R(3, 3, weight, unitaxes, refaxes, A);
463 for (j = 0; j < 3; j++)
465 mvmul(A, refaxes[j], rot_refaxes[j]);
466 mvmul(A, newaxes[j], rot_newaxes[j]);
469 /* Determine local rotation matrix A */
470 calc_fit_R(3, 3, weight, rot_newaxes, rot_refaxes, A);
471 /* Calculate euler angles, from rotation order y-z-x, where
472 * x is helixaxis, y residuevector, and z axis3.
474 * A contains rotation column vectors.
478 printf("frame %d, i=%d, A: %g %g %g , %g %g %g , %g %g %g\n",
479 teller,i,A[0][0],A[0][1],A[0][2],A[1][0],A[1][1],A[1][2],A[2][0],A[2][1],A[2][2]);
482 theta1 = 180.0/M_PI*atan2(A[0][2], A[0][0]);
483 theta2 = 180.0/M_PI*asin(-A[0][1]);
484 theta3 = 180.0/M_PI*atan2(A[2][1], A[1][1]);
486 tilt = sqrt(theta1*theta1+theta2*theta2);
488 fprintf(fptheta1, "%15.12g ", theta1);
489 fprintf(fptheta2, "%15.12g ", theta2);
490 fprintf(fptheta3, "%15.12g ", theta3);
493 fprintf(fptilt, "%15.12g ", tilt);
494 fprintf(fprotation, "%15.12g ", rotation);
496 fprintf(fptilt, "\n");
497 fprintf(fprotation, "\n");
498 fprintf(fptheta1, "\n");
499 fprintf(fptheta2, "\n");
500 fprintf(fptheta3, "\n");
503 for (i = 0; i < iCA; i++)
505 copy_rvec(residuehelixaxis[i], residuehelixaxis_tlast[i]);
506 copy_rvec(residuevector[i], residuevector_tlast[i]);
507 copy_rvec(axis3[i], axis3_tlast[i]);
512 while (read_next_x(oenv, status, &t, x, box));
514 gmx_rmpbc_done(gpbc);
517 gmx_ffclose(fpcenter);
519 gmx_ffclose(fprotation);
521 gmx_ffclose(fpradius);
522 gmx_ffclose(fptwist);
523 gmx_ffclose(fpbending);
524 gmx_ffclose(fptheta1);
525 gmx_ffclose(fptheta2);
526 gmx_ffclose(fptheta3);