3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
50 #include "gmx_fatal.h"
58 int gmx_helixorient(int argc, char *argv[])
60 const char *desc[] = {
61 "[TT]g_helixorient[tt] 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 { efTPX, 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 parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_BE_NICE,
167 NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv);
169 top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC);
171 for (i = 0; i < 3; i++)
176 /* read index files */
177 printf("Select a group of Calpha atoms corresponding to a single continuous helix:\n");
178 get_index(&(top->atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &iCA, &ind_CA, &gn_CA);
180 snew(x_SC, iCA); /* sic! */
187 snew(helixaxis, iCA-3);
189 snew(residuetwist, iCA);
191 snew(residueradius, iCA);
193 snew(residuerise, iCA);
194 snew(residueorigin, iCA);
195 snew(residuehelixaxis, iCA);
196 snew(residuevector, iCA);
197 snew(sidechainvector, iCA);
198 snew(residuebending, iCA);
199 snew(residuehelixaxis_t0, iCA);
200 snew(residuevector_t0, iCA);
202 snew(residuehelixaxis_tlast, iCA);
203 snew(residuevector_tlast, iCA);
204 snew(axis3_tlast, iCA);
209 printf("Select a group of atoms defining the sidechain direction (1/residue):\n");
210 get_index(&(top->atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &iSC, &ind_SC, &gn_SC);
213 gmx_fatal(FARGS, "Number of sidechain atoms (%d) != number of CA atoms (%d)", iSC, iCA);
218 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
220 fpaxis = ffopen(opt2fn("-oaxis", NFILE, fnm), "w");
221 fpcenter = ffopen(opt2fn("-ocenter", NFILE, fnm), "w");
222 fprise = ffopen(opt2fn("-orise", NFILE, fnm), "w");
223 fpradius = ffopen(opt2fn("-oradius", NFILE, fnm), "w");
224 fptwist = ffopen(opt2fn("-otwist", NFILE, fnm), "w");
225 fpbending = ffopen(opt2fn("-obending", NFILE, fnm), "w");
227 fptheta1 = ffopen("theta1.xvg", "w");
228 fptheta2 = ffopen("theta2.xvg", "w");
229 fptheta3 = ffopen("theta3.xvg", "w");
233 fptilt = xvgropen(opt2fn("-otilt", NFILE, fnm),
234 "Incremental local helix tilt", "Time(ps)", "Tilt (degrees)",
236 fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
237 "Incremental local helix rotation", "Time(ps)",
238 "Rotation (degrees)", oenv);
242 fptilt = xvgropen(opt2fn("-otilt", NFILE, fnm),
243 "Cumulative local helix tilt", "Time(ps)", "Tilt (degrees)", oenv);
244 fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
245 "Cumulative local helix rotation", "Time(ps)",
246 "Rotation (degrees)", oenv);
249 clear_rvecs(3, unitaxes);
254 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms, box);
258 /* initialisation for correct distance calculations */
259 set_pbc(&pbc, ePBC, box);
260 /* make molecules whole again */
261 gmx_rmpbc(gpbc, natoms, box, x);
263 /* copy coords to our smaller arrays */
264 for (i = 0; i < iCA; i++)
266 copy_rvec(x[ind_CA[i]], x_CA[i]);
269 copy_rvec(x[ind_SC[i]], x_SC[i]);
273 for (i = 0; i < iCA-3; i++)
275 rvec_sub(x_CA[i+1], x_CA[i], r12[i]);
276 rvec_sub(x_CA[i+2], x_CA[i+1], r23[i]);
277 rvec_sub(x_CA[i+3], x_CA[i+2], r34[i]);
278 rvec_sub(r12[i], r23[i], diff13[i]);
279 rvec_sub(r23[i], r34[i], diff24[i]);
280 /* calculate helix axis */
281 cprod(diff13[i], diff24[i], helixaxis[i]);
282 svmul(1.0/norm(helixaxis[i]), helixaxis[i], helixaxis[i]);
284 tmp = cos_angle(diff13[i], diff24[i]);
285 twist[i] = 180.0/M_PI * acos( tmp );
286 radius[i] = sqrt( norm(diff13[i])*norm(diff24[i]) ) / (2.0* (1.0-tmp) );
287 rise[i] = fabs(iprod(r23[i], helixaxis[i]));
289 svmul(radius[i]/norm(diff13[i]), diff13[i], v1);
290 svmul(radius[i]/norm(diff24[i]), diff24[i], v2);
292 rvec_sub(x_CA[i+1], v1, residueorigin[i+1]);
293 rvec_sub(x_CA[i+2], v2, residueorigin[i+2]);
295 residueradius[0] = residuetwist[0] = residuerise[0] = 0;
297 residueradius[1] = radius[0];
298 residuetwist[1] = twist[0];
299 residuerise[1] = rise[0];
301 residuebending[0] = residuebending[1] = 0;
302 for (i = 2; i < iCA-2; i++)
304 residueradius[i] = 0.5*(radius[i-2]+radius[i-1]);
305 residuetwist[i] = 0.5*(twist[i-2]+twist[i-1]);
306 residuerise[i] = 0.5*(rise[i-2]+rise[i-1]);
307 residuebending[i] = 180.0/M_PI*acos( cos_angle(helixaxis[i-2], helixaxis[i-1]) );
309 residueradius[iCA-2] = radius[iCA-4];
310 residuetwist[iCA-2] = twist[iCA-4];
311 residuerise[iCA-2] = rise[iCA-4];
312 residueradius[iCA-1] = residuetwist[iCA-1] = residuerise[iCA-1] = 0;
313 residuebending[iCA-2] = residuebending[iCA-1] = 0;
315 clear_rvec(residueorigin[0]);
316 clear_rvec(residueorigin[iCA-1]);
318 /* average helix axes to define them on the residues.
319 * Just extrapolate second first/list atom.
321 copy_rvec(helixaxis[0], residuehelixaxis[0]);
322 copy_rvec(helixaxis[0], residuehelixaxis[1]);
324 for (i = 2; i < iCA-2; i++)
326 rvec_add(helixaxis[i-2], helixaxis[i-1], residuehelixaxis[i]);
327 svmul(0.5, residuehelixaxis[i], residuehelixaxis[i]);
329 copy_rvec(helixaxis[iCA-4], residuehelixaxis[iCA-2]);
330 copy_rvec(helixaxis[iCA-4], residuehelixaxis[iCA-1]);
332 /* Normalize the axis */
333 for (i = 0; i < iCA; i++)
335 svmul(1.0/norm(residuehelixaxis[i]), residuehelixaxis[i], residuehelixaxis[i]);
338 /* calculate vector from origin to residue CA */
339 fprintf(fpaxis, "%15.12g ", t);
340 fprintf(fpcenter, "%15.12g ", t);
341 fprintf(fprise, "%15.12g ", t);
342 fprintf(fpradius, "%15.12g ", t);
343 fprintf(fptwist, "%15.12g ", t);
344 fprintf(fpbending, "%15.12g ", t);
346 for (i = 0; i < iCA; i++)
348 if (i == 0 || i == iCA-1)
350 fprintf(fpaxis, "%15.12g %15.12g %15.12g ", 0.0, 0.0, 0.0);
351 fprintf(fpcenter, "%15.12g %15.12g %15.12g ", 0.0, 0.0, 0.0);
352 fprintf(fprise, "%15.12g ", 0.0);
353 fprintf(fpradius, "%15.12g ", 0.0);
354 fprintf(fptwist, "%15.12g ", 0.0);
355 fprintf(fpbending, "%15.12g ", 0.0);
359 rvec_sub( bSC ? x_SC[i] : x_CA[i], residueorigin[i], residuevector[i]);
360 svmul(1.0/norm(residuevector[i]), residuevector[i], residuevector[i]);
361 cprod(residuehelixaxis[i], residuevector[i], axis3[i]);
362 fprintf(fpaxis, "%15.12g %15.12g %15.12g ", residuehelixaxis[i][0], residuehelixaxis[i][1], residuehelixaxis[i][2]);
363 fprintf(fpcenter, "%15.12g %15.12g %15.12g ", residueorigin[i][0], residueorigin[i][1], residueorigin[i][2]);
365 fprintf(fprise, "%15.12g ", residuerise[i]);
366 fprintf(fpradius, "%15.12g ", residueradius[i]);
367 fprintf(fptwist, "%15.12g ", residuetwist[i]);
368 fprintf(fpbending, "%15.12g ", residuebending[i]);
370 /* angle with local vector? */
372 printf("res[%2d]: axis: %g %g %g origin: %g %g %g vector: %g %g %g angle: %g\n",i,
373 residuehelixaxis[i][0],
374 residuehelixaxis[i][1],
375 residuehelixaxis[i][2],
382 180.0/M_PI*acos( cos_angle(residuevector[i],residuehelixaxis[i]) ));
384 /* fprintf(fp,"%15.12g %15.12g %15.12g %15.12g %15.12g %15.12g\n",
385 residuehelixaxis[i][0],
386 residuehelixaxis[i][1],
387 residuehelixaxis[i][2],
390 residuevector[i][2]);
394 fprintf(fprise, "\n");
395 fprintf(fpradius, "\n");
396 fprintf(fpaxis, "\n");
397 fprintf(fpcenter, "\n");
398 fprintf(fptwist, "\n");
399 fprintf(fpbending, "\n");
403 for (i = 0; i < iCA; i++)
405 copy_rvec(residuehelixaxis[i], residuehelixaxis_t0[i]);
406 copy_rvec(residuevector[i], residuevector_t0[i]);
407 copy_rvec(axis3[i], axis3_t0[i]);
412 fprintf(fptilt, "%15.12g ", t);
413 fprintf(fprotation, "%15.12g ", t);
414 fprintf(fptheta1, "%15.12g ", t);
415 fprintf(fptheta2, "%15.12g ", t);
416 fprintf(fptheta3, "%15.12g ", t);
418 for (i = 0; i < iCA; i++)
420 if (i == 0 || i == iCA-1)
428 /* Total rotation & tilt */
429 copy_rvec(residuehelixaxis_t0[i], refaxes[0]);
430 copy_rvec(residuevector_t0[i], refaxes[1]);
431 copy_rvec(axis3_t0[i], refaxes[2]);
435 /* Rotation/tilt since last step */
436 copy_rvec(residuehelixaxis_tlast[i], refaxes[0]);
437 copy_rvec(residuevector_tlast[i], refaxes[1]);
438 copy_rvec(axis3_tlast[i], refaxes[2]);
440 copy_rvec(residuehelixaxis[i], newaxes[0]);
441 copy_rvec(residuevector[i], newaxes[1]);
442 copy_rvec(axis3[i], newaxes[2]);
445 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",
447 refaxes[0][0],refaxes[0][1],refaxes[0][2],
448 refaxes[1][0],refaxes[1][1],refaxes[1][2],
449 refaxes[2][0],refaxes[2][1],refaxes[2][2],
450 newaxes[0][0],newaxes[0][1],newaxes[0][2],
451 newaxes[1][0],newaxes[1][1],newaxes[1][2],
452 newaxes[2][0],newaxes[2][1],newaxes[2][2]);
455 /* rotate reference frame onto unit axes */
456 calc_fit_R(3, 3, weight, unitaxes, refaxes, A);
457 for (j = 0; j < 3; j++)
459 mvmul(A, refaxes[j], rot_refaxes[j]);
460 mvmul(A, newaxes[j], rot_newaxes[j]);
463 /* Determine local rotation matrix A */
464 calc_fit_R(3, 3, weight, rot_newaxes, rot_refaxes, A);
465 /* Calculate euler angles, from rotation order y-z-x, where
466 * x is helixaxis, y residuevector, and z axis3.
468 * A contains rotation column vectors.
472 printf("frame %d, i=%d, A: %g %g %g , %g %g %g , %g %g %g\n",
473 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]);
476 theta1 = 180.0/M_PI*atan2(A[0][2], A[0][0]);
477 theta2 = 180.0/M_PI*asin(-A[0][1]);
478 theta3 = 180.0/M_PI*atan2(A[2][1], A[1][1]);
480 tilt = sqrt(theta1*theta1+theta2*theta2);
482 fprintf(fptheta1, "%15.12g ", theta1);
483 fprintf(fptheta2, "%15.12g ", theta2);
484 fprintf(fptheta3, "%15.12g ", theta3);
487 fprintf(fptilt, "%15.12g ", tilt);
488 fprintf(fprotation, "%15.12g ", rotation);
490 fprintf(fptilt, "\n");
491 fprintf(fprotation, "\n");
492 fprintf(fptheta1, "\n");
493 fprintf(fptheta2, "\n");
494 fprintf(fptheta3, "\n");
497 for (i = 0; i < iCA; i++)
499 copy_rvec(residuehelixaxis[i], residuehelixaxis_tlast[i]);
500 copy_rvec(residuevector[i], residuevector_tlast[i]);
501 copy_rvec(axis3[i], axis3_tlast[i]);
506 while (read_next_x(oenv, status, &t, natoms, x, box));
508 gmx_rmpbc_done(gpbc);