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.
41 #include "gromacs/legacyheaders/typedefs.h"
42 #include "gromacs/legacyheaders/macros.h"
43 #include "gromacs/math/vec.h"
44 #include "gromacs/topology/index.h"
48 #include "gromacs/commandline/pargs.h"
49 #include "gromacs/fileio/trxio.h"
50 #include "gromacs/fileio/xvgr.h"
51 #include "gromacs/math/do_fit.h"
52 #include "gromacs/pbcutil/pbc.h"
53 #include "gromacs/pbcutil/rmpbc.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 [IT]y[it], and the",
73 "[IT]z[it]-axis from their cross product. We use the Euler Y-Z-X rotation, meaning",
74 "we first tilt the helix axis (1) around and (2) orthogonal to the residues",
75 "vector, and finally apply the (3) rotation around it. For debugging or other",
76 "purposes, we also write out the actual Euler rotation angles as [TT]theta[1-3].xvg[tt]"
79 t_topology *top = NULL;
85 real theta1, theta2, theta3;
87 int d, i, j, teller = 0;
94 rvec v1, v2, p1, p2, vtmp, vproj;
102 rvec *residuehelixaxis;
105 rvec *sidechainvector;
109 rvec *residuehelixaxis_t0;
110 rvec *residuevector_t0;
112 rvec *residuehelixaxis_tlast;
113 rvec *residuevector_tlast;
115 rvec refaxes[3], newaxes[3];
117 rvec rot_refaxes[3], rot_newaxes[3];
121 real *twist, *residuetwist;
122 real *radius, *residueradius;
123 real *rise, *residuerise;
124 real *residuebending;
131 FILE *fpaxis, *fpcenter, *fptilt, *fprotation;
132 FILE *fpradius, *fprise, *fptwist;
133 FILE *fptheta1, *fptheta2, *fptheta3;
138 gmx_rmpbc_t gpbc = NULL;
140 static gmx_bool bSC = FALSE;
141 static gmx_bool bIncremental = FALSE;
143 static t_pargs pa[] = {
144 { "-sidechain", FALSE, etBOOL, {&bSC},
145 "Calculate sidechain directions relative to helix axis too." },
146 { "-incremental", FALSE, etBOOL, {&bIncremental},
147 "Calculate incremental rather than total rotation/tilt." },
149 #define NPA asize(pa)
152 { efTPR, NULL, NULL, ffREAD },
153 { efTRX, "-f", NULL, ffREAD },
154 { efNDX, NULL, NULL, ffOPTRD },
155 { efDAT, "-oaxis", "helixaxis", ffWRITE },
156 { efDAT, "-ocenter", "center", ffWRITE },
157 { efXVG, "-orise", "rise", ffWRITE },
158 { efXVG, "-oradius", "radius", ffWRITE },
159 { efXVG, "-otwist", "twist", ffWRITE },
160 { efXVG, "-obending", "bending", ffWRITE },
161 { efXVG, "-otilt", "tilt", ffWRITE },
162 { efXVG, "-orot", "rotation", ffWRITE }
164 #define NFILE asize(fnm)
166 if (!parse_common_args(&argc, argv, PCA_CAN_TIME,
167 NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
172 top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC);
174 for (i = 0; i < 3; i++)
179 /* read index files */
180 printf("Select a group of Calpha atoms corresponding to a single continuous helix:\n");
181 get_index(&(top->atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &iCA, &ind_CA, &gn_CA);
183 snew(x_SC, iCA); /* sic! */
190 snew(helixaxis, iCA-3);
192 snew(residuetwist, iCA);
194 snew(residueradius, iCA);
196 snew(residuerise, iCA);
197 snew(residueorigin, iCA);
198 snew(residuehelixaxis, iCA);
199 snew(residuevector, iCA);
200 snew(sidechainvector, iCA);
201 snew(residuebending, iCA);
202 snew(residuehelixaxis_t0, iCA);
203 snew(residuevector_t0, iCA);
205 snew(residuehelixaxis_tlast, iCA);
206 snew(residuevector_tlast, iCA);
207 snew(axis3_tlast, iCA);
212 printf("Select a group of atoms defining the sidechain direction (1/residue):\n");
213 get_index(&(top->atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &iSC, &ind_SC, &gn_SC);
216 gmx_fatal(FARGS, "Number of sidechain atoms (%d) != number of CA atoms (%d)", iSC, iCA);
221 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
223 fpaxis = gmx_ffopen(opt2fn("-oaxis", NFILE, fnm), "w");
224 fpcenter = gmx_ffopen(opt2fn("-ocenter", NFILE, fnm), "w");
225 fprise = gmx_ffopen(opt2fn("-orise", NFILE, fnm), "w");
226 fpradius = gmx_ffopen(opt2fn("-oradius", NFILE, fnm), "w");
227 fptwist = gmx_ffopen(opt2fn("-otwist", NFILE, fnm), "w");
228 fpbending = gmx_ffopen(opt2fn("-obending", NFILE, fnm), "w");
230 fptheta1 = gmx_ffopen("theta1.xvg", "w");
231 fptheta2 = gmx_ffopen("theta2.xvg", "w");
232 fptheta3 = gmx_ffopen("theta3.xvg", "w");
236 fptilt = xvgropen(opt2fn("-otilt", NFILE, fnm),
237 "Incremental local helix tilt", "Time(ps)", "Tilt (degrees)",
239 fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
240 "Incremental local helix rotation", "Time(ps)",
241 "Rotation (degrees)", oenv);
245 fptilt = xvgropen(opt2fn("-otilt", NFILE, fnm),
246 "Cumulative local helix tilt", "Time(ps)", "Tilt (degrees)", oenv);
247 fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
248 "Cumulative local helix rotation", "Time(ps)",
249 "Rotation (degrees)", oenv);
252 clear_rvecs(3, unitaxes);
257 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
261 /* initialisation for correct distance calculations */
262 set_pbc(&pbc, ePBC, box);
263 /* make molecules whole again */
264 gmx_rmpbc(gpbc, natoms, box, x);
266 /* copy coords to our smaller arrays */
267 for (i = 0; i < iCA; i++)
269 copy_rvec(x[ind_CA[i]], x_CA[i]);
272 copy_rvec(x[ind_SC[i]], x_SC[i]);
276 for (i = 0; i < iCA-3; i++)
278 rvec_sub(x_CA[i+1], x_CA[i], r12[i]);
279 rvec_sub(x_CA[i+2], x_CA[i+1], r23[i]);
280 rvec_sub(x_CA[i+3], x_CA[i+2], r34[i]);
281 rvec_sub(r12[i], r23[i], diff13[i]);
282 rvec_sub(r23[i], r34[i], diff24[i]);
283 /* calculate helix axis */
284 cprod(diff13[i], diff24[i], helixaxis[i]);
285 svmul(1.0/norm(helixaxis[i]), helixaxis[i], helixaxis[i]);
287 tmp = cos_angle(diff13[i], diff24[i]);
288 twist[i] = 180.0/M_PI * acos( tmp );
289 radius[i] = sqrt( norm(diff13[i])*norm(diff24[i]) ) / (2.0* (1.0-tmp) );
290 rise[i] = fabs(iprod(r23[i], helixaxis[i]));
292 svmul(radius[i]/norm(diff13[i]), diff13[i], v1);
293 svmul(radius[i]/norm(diff24[i]), diff24[i], v2);
295 rvec_sub(x_CA[i+1], v1, residueorigin[i+1]);
296 rvec_sub(x_CA[i+2], v2, residueorigin[i+2]);
298 residueradius[0] = residuetwist[0] = residuerise[0] = 0;
300 residueradius[1] = radius[0];
301 residuetwist[1] = twist[0];
302 residuerise[1] = rise[0];
304 residuebending[0] = residuebending[1] = 0;
305 for (i = 2; i < iCA-2; i++)
307 residueradius[i] = 0.5*(radius[i-2]+radius[i-1]);
308 residuetwist[i] = 0.5*(twist[i-2]+twist[i-1]);
309 residuerise[i] = 0.5*(rise[i-2]+rise[i-1]);
310 residuebending[i] = 180.0/M_PI*acos( cos_angle(helixaxis[i-2], helixaxis[i-1]) );
312 residueradius[iCA-2] = radius[iCA-4];
313 residuetwist[iCA-2] = twist[iCA-4];
314 residuerise[iCA-2] = rise[iCA-4];
315 residueradius[iCA-1] = residuetwist[iCA-1] = residuerise[iCA-1] = 0;
316 residuebending[iCA-2] = residuebending[iCA-1] = 0;
318 clear_rvec(residueorigin[0]);
319 clear_rvec(residueorigin[iCA-1]);
321 /* average helix axes to define them on the residues.
322 * Just extrapolate second first/list atom.
324 copy_rvec(helixaxis[0], residuehelixaxis[0]);
325 copy_rvec(helixaxis[0], residuehelixaxis[1]);
327 for (i = 2; i < iCA-2; i++)
329 rvec_add(helixaxis[i-2], helixaxis[i-1], residuehelixaxis[i]);
330 svmul(0.5, residuehelixaxis[i], residuehelixaxis[i]);
332 copy_rvec(helixaxis[iCA-4], residuehelixaxis[iCA-2]);
333 copy_rvec(helixaxis[iCA-4], residuehelixaxis[iCA-1]);
335 /* Normalize the axis */
336 for (i = 0; i < iCA; i++)
338 svmul(1.0/norm(residuehelixaxis[i]), residuehelixaxis[i], residuehelixaxis[i]);
341 /* calculate vector from origin to residue CA */
342 fprintf(fpaxis, "%15.12g ", t);
343 fprintf(fpcenter, "%15.12g ", t);
344 fprintf(fprise, "%15.12g ", t);
345 fprintf(fpradius, "%15.12g ", t);
346 fprintf(fptwist, "%15.12g ", t);
347 fprintf(fpbending, "%15.12g ", t);
349 for (i = 0; i < iCA; i++)
351 if (i == 0 || i == iCA-1)
353 fprintf(fpaxis, "%15.12g %15.12g %15.12g ", 0.0, 0.0, 0.0);
354 fprintf(fpcenter, "%15.12g %15.12g %15.12g ", 0.0, 0.0, 0.0);
355 fprintf(fprise, "%15.12g ", 0.0);
356 fprintf(fpradius, "%15.12g ", 0.0);
357 fprintf(fptwist, "%15.12g ", 0.0);
358 fprintf(fpbending, "%15.12g ", 0.0);
362 rvec_sub( bSC ? x_SC[i] : x_CA[i], residueorigin[i], residuevector[i]);
363 svmul(1.0/norm(residuevector[i]), residuevector[i], residuevector[i]);
364 cprod(residuehelixaxis[i], residuevector[i], axis3[i]);
365 fprintf(fpaxis, "%15.12g %15.12g %15.12g ", residuehelixaxis[i][0], residuehelixaxis[i][1], residuehelixaxis[i][2]);
366 fprintf(fpcenter, "%15.12g %15.12g %15.12g ", residueorigin[i][0], residueorigin[i][1], residueorigin[i][2]);
368 fprintf(fprise, "%15.12g ", residuerise[i]);
369 fprintf(fpradius, "%15.12g ", residueradius[i]);
370 fprintf(fptwist, "%15.12g ", residuetwist[i]);
371 fprintf(fpbending, "%15.12g ", residuebending[i]);
373 /* angle with local vector? */
375 printf("res[%2d]: axis: %g %g %g origin: %g %g %g vector: %g %g %g angle: %g\n",i,
376 residuehelixaxis[i][0],
377 residuehelixaxis[i][1],
378 residuehelixaxis[i][2],
385 180.0/M_PI*acos( cos_angle(residuevector[i],residuehelixaxis[i]) ));
387 /* fprintf(fp,"%15.12g %15.12g %15.12g %15.12g %15.12g %15.12g\n",
388 residuehelixaxis[i][0],
389 residuehelixaxis[i][1],
390 residuehelixaxis[i][2],
393 residuevector[i][2]);
397 fprintf(fprise, "\n");
398 fprintf(fpradius, "\n");
399 fprintf(fpaxis, "\n");
400 fprintf(fpcenter, "\n");
401 fprintf(fptwist, "\n");
402 fprintf(fpbending, "\n");
406 for (i = 0; i < iCA; i++)
408 copy_rvec(residuehelixaxis[i], residuehelixaxis_t0[i]);
409 copy_rvec(residuevector[i], residuevector_t0[i]);
410 copy_rvec(axis3[i], axis3_t0[i]);
415 fprintf(fptilt, "%15.12g ", t);
416 fprintf(fprotation, "%15.12g ", t);
417 fprintf(fptheta1, "%15.12g ", t);
418 fprintf(fptheta2, "%15.12g ", t);
419 fprintf(fptheta3, "%15.12g ", t);
421 for (i = 0; i < iCA; i++)
423 if (i == 0 || i == iCA-1)
431 /* Total rotation & tilt */
432 copy_rvec(residuehelixaxis_t0[i], refaxes[0]);
433 copy_rvec(residuevector_t0[i], refaxes[1]);
434 copy_rvec(axis3_t0[i], refaxes[2]);
438 /* Rotation/tilt since last step */
439 copy_rvec(residuehelixaxis_tlast[i], refaxes[0]);
440 copy_rvec(residuevector_tlast[i], refaxes[1]);
441 copy_rvec(axis3_tlast[i], refaxes[2]);
443 copy_rvec(residuehelixaxis[i], newaxes[0]);
444 copy_rvec(residuevector[i], newaxes[1]);
445 copy_rvec(axis3[i], newaxes[2]);
448 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",
450 refaxes[0][0],refaxes[0][1],refaxes[0][2],
451 refaxes[1][0],refaxes[1][1],refaxes[1][2],
452 refaxes[2][0],refaxes[2][1],refaxes[2][2],
453 newaxes[0][0],newaxes[0][1],newaxes[0][2],
454 newaxes[1][0],newaxes[1][1],newaxes[1][2],
455 newaxes[2][0],newaxes[2][1],newaxes[2][2]);
458 /* rotate reference frame onto unit axes */
459 calc_fit_R(3, 3, weight, unitaxes, refaxes, A);
460 for (j = 0; j < 3; j++)
462 mvmul(A, refaxes[j], rot_refaxes[j]);
463 mvmul(A, newaxes[j], rot_newaxes[j]);
466 /* Determine local rotation matrix A */
467 calc_fit_R(3, 3, weight, rot_newaxes, rot_refaxes, A);
468 /* Calculate euler angles, from rotation order y-z-x, where
469 * x is helixaxis, y residuevector, and z axis3.
471 * A contains rotation column vectors.
475 printf("frame %d, i=%d, A: %g %g %g , %g %g %g , %g %g %g\n",
476 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]);
479 theta1 = 180.0/M_PI*atan2(A[0][2], A[0][0]);
480 theta2 = 180.0/M_PI*asin(-A[0][1]);
481 theta3 = 180.0/M_PI*atan2(A[2][1], A[1][1]);
483 tilt = sqrt(theta1*theta1+theta2*theta2);
485 fprintf(fptheta1, "%15.12g ", theta1);
486 fprintf(fptheta2, "%15.12g ", theta2);
487 fprintf(fptheta3, "%15.12g ", theta3);
490 fprintf(fptilt, "%15.12g ", tilt);
491 fprintf(fprotation, "%15.12g ", rotation);
493 fprintf(fptilt, "\n");
494 fprintf(fprotation, "\n");
495 fprintf(fptheta1, "\n");
496 fprintf(fptheta2, "\n");
497 fprintf(fptheta3, "\n");
500 for (i = 0; i < iCA; i++)
502 copy_rvec(residuehelixaxis[i], residuehelixaxis_tlast[i]);
503 copy_rvec(residuevector[i], residuevector_tlast[i]);
504 copy_rvec(axis3[i], axis3_tlast[i]);
509 while (read_next_x(oenv, status, &t, x, box));
511 gmx_rmpbc_done(gpbc);
514 gmx_ffclose(fpcenter);
516 gmx_ffclose(fprotation);
518 gmx_ffclose(fpradius);
519 gmx_ffclose(fptwist);
520 gmx_ffclose(fpbending);
521 gmx_ffclose(fptheta1);
522 gmx_ffclose(fptheta2);
523 gmx_ffclose(fptheta3);