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.
43 #include "gromacs/legacyheaders/typedefs.h"
44 #include "gromacs/legacyheaders/macros.h"
45 #include "gromacs/math/vec.h"
46 #include "gromacs/topology/index.h"
50 #include "gromacs/commandline/pargs.h"
51 #include "gromacs/fileio/trxio.h"
52 #include "gromacs/fileio/xvgr.h"
53 #include "gromacs/math/do_fit.h"
54 #include "gromacs/pbcutil/pbc.h"
55 #include "gromacs/pbcutil/rmpbc.h"
56 #include "gromacs/utility/fatalerror.h"
57 #include "gromacs/utility/futil.h"
58 #include "gromacs/utility/smalloc.h"
60 int gmx_helixorient(int argc
, char *argv
[])
62 const char *desc
[] = {
63 "[THISMODULE] calculates the coordinates and direction of the average",
64 "axis inside an alpha helix, and the direction/vectors of both the",
65 "C[GRK]alpha[grk] and (optionally) a sidechain atom relative to the axis.[PAR]",
66 "As input, you need to specify an index group with C[GRK]alpha[grk] atoms",
67 "corresponding to an [GRK]alpha[grk]-helix of continuous residues. Sidechain",
68 "directions require a second index group of the same size, containing",
69 "the heavy atom in each residue that should represent the sidechain.[PAR]",
70 "[BB]Note[bb] that this program does not do any fitting of structures.[PAR]",
71 "We need four C[GRK]alpha[grk] coordinates to define the local direction of the helix",
73 "The tilt/rotation is calculated from Euler rotations, where we define",
74 "the helix axis as the local [IT]x[it]-axis, the residues/C[GRK]alpha[grk] vector as [IT]y[it], and the",
75 "[IT]z[it]-axis from their cross product. We use the Euler Y-Z-X rotation, meaning",
76 "we first tilt the helix axis (1) around and (2) orthogonal to the residues",
77 "vector, and finally apply the (3) rotation around it. For debugging or other",
78 "purposes, we also write out the actual Euler rotation angles as [TT]theta[1-3].xvg[tt]"
81 t_topology
*top
= NULL
;
87 real theta1
, theta2
, theta3
;
89 int d
, i
, j
, teller
= 0;
96 rvec v1
, v2
, p1
, p2
, vtmp
, vproj
;
104 rvec
*residuehelixaxis
;
107 rvec
*sidechainvector
;
111 rvec
*residuehelixaxis_t0
;
112 rvec
*residuevector_t0
;
114 rvec
*residuehelixaxis_tlast
;
115 rvec
*residuevector_tlast
;
117 rvec refaxes
[3], newaxes
[3];
119 rvec rot_refaxes
[3], rot_newaxes
[3];
123 real
*twist
, *residuetwist
;
124 real
*radius
, *residueradius
;
125 real
*rise
, *residuerise
;
126 real
*residuebending
;
133 FILE *fpaxis
, *fpcenter
, *fptilt
, *fprotation
;
134 FILE *fpradius
, *fprise
, *fptwist
;
135 FILE *fptheta1
, *fptheta2
, *fptheta3
;
140 gmx_rmpbc_t gpbc
= NULL
;
142 static gmx_bool bSC
= FALSE
;
143 static gmx_bool bIncremental
= FALSE
;
145 static t_pargs pa
[] = {
146 { "-sidechain", FALSE
, etBOOL
, {&bSC
},
147 "Calculate sidechain directions relative to helix axis too." },
148 { "-incremental", FALSE
, etBOOL
, {&bIncremental
},
149 "Calculate incremental rather than total rotation/tilt." },
151 #define NPA asize(pa)
154 { efTPX
, NULL
, NULL
, ffREAD
},
155 { efTRX
, "-f", NULL
, ffREAD
},
156 { efNDX
, NULL
, NULL
, ffOPTRD
},
157 { efDAT
, "-oaxis", "helixaxis", ffWRITE
},
158 { efDAT
, "-ocenter", "center", ffWRITE
},
159 { efXVG
, "-orise", "rise", ffWRITE
},
160 { efXVG
, "-oradius", "radius", ffWRITE
},
161 { efXVG
, "-otwist", "twist", ffWRITE
},
162 { efXVG
, "-obending", "bending", ffWRITE
},
163 { efXVG
, "-otilt", "tilt", ffWRITE
},
164 { efXVG
, "-orot", "rotation", ffWRITE
}
166 #define NFILE asize(fnm)
168 if (!parse_common_args(&argc
, argv
, PCA_CAN_TIME
| PCA_BE_NICE
,
169 NFILE
, fnm
, NPA
, pa
, asize(desc
), desc
, 0, NULL
, &oenv
))
174 top
= read_top(ftp2fn(efTPX
, NFILE
, fnm
), &ePBC
);
176 for (i
= 0; i
< 3; i
++)
181 /* read index files */
182 printf("Select a group of Calpha atoms corresponding to a single continuous helix:\n");
183 get_index(&(top
->atoms
), ftp2fn_null(efNDX
, NFILE
, fnm
), 1, &iCA
, &ind_CA
, &gn_CA
);
185 snew(x_SC
, iCA
); /* sic! */
192 snew(helixaxis
, iCA
-3);
194 snew(residuetwist
, iCA
);
196 snew(residueradius
, iCA
);
198 snew(residuerise
, iCA
);
199 snew(residueorigin
, iCA
);
200 snew(residuehelixaxis
, iCA
);
201 snew(residuevector
, iCA
);
202 snew(sidechainvector
, iCA
);
203 snew(residuebending
, iCA
);
204 snew(residuehelixaxis_t0
, iCA
);
205 snew(residuevector_t0
, iCA
);
207 snew(residuehelixaxis_tlast
, iCA
);
208 snew(residuevector_tlast
, iCA
);
209 snew(axis3_tlast
, iCA
);
214 printf("Select a group of atoms defining the sidechain direction (1/residue):\n");
215 get_index(&(top
->atoms
), ftp2fn_null(efNDX
, NFILE
, fnm
), 1, &iSC
, &ind_SC
, &gn_SC
);
218 gmx_fatal(FARGS
, "Number of sidechain atoms (%d) != number of CA atoms (%d)", iSC
, iCA
);
223 natoms
= read_first_x(oenv
, &status
, ftp2fn(efTRX
, NFILE
, fnm
), &t
, &x
, box
);
225 fpaxis
= gmx_ffopen(opt2fn("-oaxis", NFILE
, fnm
), "w");
226 fpcenter
= gmx_ffopen(opt2fn("-ocenter", NFILE
, fnm
), "w");
227 fprise
= gmx_ffopen(opt2fn("-orise", NFILE
, fnm
), "w");
228 fpradius
= gmx_ffopen(opt2fn("-oradius", NFILE
, fnm
), "w");
229 fptwist
= gmx_ffopen(opt2fn("-otwist", NFILE
, fnm
), "w");
230 fpbending
= gmx_ffopen(opt2fn("-obending", NFILE
, fnm
), "w");
232 fptheta1
= gmx_ffopen("theta1.xvg", "w");
233 fptheta2
= gmx_ffopen("theta2.xvg", "w");
234 fptheta3
= gmx_ffopen("theta3.xvg", "w");
238 fptilt
= xvgropen(opt2fn("-otilt", NFILE
, fnm
),
239 "Incremental local helix tilt", "Time(ps)", "Tilt (degrees)",
241 fprotation
= xvgropen(opt2fn("-orot", NFILE
, fnm
),
242 "Incremental local helix rotation", "Time(ps)",
243 "Rotation (degrees)", oenv
);
247 fptilt
= xvgropen(opt2fn("-otilt", NFILE
, fnm
),
248 "Cumulative local helix tilt", "Time(ps)", "Tilt (degrees)", oenv
);
249 fprotation
= xvgropen(opt2fn("-orot", NFILE
, fnm
),
250 "Cumulative local helix rotation", "Time(ps)",
251 "Rotation (degrees)", oenv
);
254 clear_rvecs(3, unitaxes
);
259 gpbc
= gmx_rmpbc_init(&top
->idef
, ePBC
, natoms
);
263 /* initialisation for correct distance calculations */
264 set_pbc(&pbc
, ePBC
, box
);
265 /* make molecules whole again */
266 gmx_rmpbc(gpbc
, natoms
, box
, x
);
268 /* copy coords to our smaller arrays */
269 for (i
= 0; i
< iCA
; i
++)
271 copy_rvec(x
[ind_CA
[i
]], x_CA
[i
]);
274 copy_rvec(x
[ind_SC
[i
]], x_SC
[i
]);
278 for (i
= 0; i
< iCA
-3; i
++)
280 rvec_sub(x_CA
[i
+1], x_CA
[i
], r12
[i
]);
281 rvec_sub(x_CA
[i
+2], x_CA
[i
+1], r23
[i
]);
282 rvec_sub(x_CA
[i
+3], x_CA
[i
+2], r34
[i
]);
283 rvec_sub(r12
[i
], r23
[i
], diff13
[i
]);
284 rvec_sub(r23
[i
], r34
[i
], diff24
[i
]);
285 /* calculate helix axis */
286 cprod(diff13
[i
], diff24
[i
], helixaxis
[i
]);
287 svmul(1.0/norm(helixaxis
[i
]), helixaxis
[i
], helixaxis
[i
]);
289 tmp
= cos_angle(diff13
[i
], diff24
[i
]);
290 twist
[i
] = 180.0/M_PI
* acos( tmp
);
291 radius
[i
] = sqrt( norm(diff13
[i
])*norm(diff24
[i
]) ) / (2.0* (1.0-tmp
) );
292 rise
[i
] = fabs(iprod(r23
[i
], helixaxis
[i
]));
294 svmul(radius
[i
]/norm(diff13
[i
]), diff13
[i
], v1
);
295 svmul(radius
[i
]/norm(diff24
[i
]), diff24
[i
], v2
);
297 rvec_sub(x_CA
[i
+1], v1
, residueorigin
[i
+1]);
298 rvec_sub(x_CA
[i
+2], v2
, residueorigin
[i
+2]);
300 residueradius
[0] = residuetwist
[0] = residuerise
[0] = 0;
302 residueradius
[1] = radius
[0];
303 residuetwist
[1] = twist
[0];
304 residuerise
[1] = rise
[0];
306 residuebending
[0] = residuebending
[1] = 0;
307 for (i
= 2; i
< iCA
-2; i
++)
309 residueradius
[i
] = 0.5*(radius
[i
-2]+radius
[i
-1]);
310 residuetwist
[i
] = 0.5*(twist
[i
-2]+twist
[i
-1]);
311 residuerise
[i
] = 0.5*(rise
[i
-2]+rise
[i
-1]);
312 residuebending
[i
] = 180.0/M_PI
*acos( cos_angle(helixaxis
[i
-2], helixaxis
[i
-1]) );
314 residueradius
[iCA
-2] = radius
[iCA
-4];
315 residuetwist
[iCA
-2] = twist
[iCA
-4];
316 residuerise
[iCA
-2] = rise
[iCA
-4];
317 residueradius
[iCA
-1] = residuetwist
[iCA
-1] = residuerise
[iCA
-1] = 0;
318 residuebending
[iCA
-2] = residuebending
[iCA
-1] = 0;
320 clear_rvec(residueorigin
[0]);
321 clear_rvec(residueorigin
[iCA
-1]);
323 /* average helix axes to define them on the residues.
324 * Just extrapolate second first/list atom.
326 copy_rvec(helixaxis
[0], residuehelixaxis
[0]);
327 copy_rvec(helixaxis
[0], residuehelixaxis
[1]);
329 for (i
= 2; i
< iCA
-2; i
++)
331 rvec_add(helixaxis
[i
-2], helixaxis
[i
-1], residuehelixaxis
[i
]);
332 svmul(0.5, residuehelixaxis
[i
], residuehelixaxis
[i
]);
334 copy_rvec(helixaxis
[iCA
-4], residuehelixaxis
[iCA
-2]);
335 copy_rvec(helixaxis
[iCA
-4], residuehelixaxis
[iCA
-1]);
337 /* Normalize the axis */
338 for (i
= 0; i
< iCA
; i
++)
340 svmul(1.0/norm(residuehelixaxis
[i
]), residuehelixaxis
[i
], residuehelixaxis
[i
]);
343 /* calculate vector from origin to residue CA */
344 fprintf(fpaxis
, "%15.12g ", t
);
345 fprintf(fpcenter
, "%15.12g ", t
);
346 fprintf(fprise
, "%15.12g ", t
);
347 fprintf(fpradius
, "%15.12g ", t
);
348 fprintf(fptwist
, "%15.12g ", t
);
349 fprintf(fpbending
, "%15.12g ", t
);
351 for (i
= 0; i
< iCA
; i
++)
353 if (i
== 0 || i
== iCA
-1)
355 fprintf(fpaxis
, "%15.12g %15.12g %15.12g ", 0.0, 0.0, 0.0);
356 fprintf(fpcenter
, "%15.12g %15.12g %15.12g ", 0.0, 0.0, 0.0);
357 fprintf(fprise
, "%15.12g ", 0.0);
358 fprintf(fpradius
, "%15.12g ", 0.0);
359 fprintf(fptwist
, "%15.12g ", 0.0);
360 fprintf(fpbending
, "%15.12g ", 0.0);
364 rvec_sub( bSC
? x_SC
[i
] : x_CA
[i
], residueorigin
[i
], residuevector
[i
]);
365 svmul(1.0/norm(residuevector
[i
]), residuevector
[i
], residuevector
[i
]);
366 cprod(residuehelixaxis
[i
], residuevector
[i
], axis3
[i
]);
367 fprintf(fpaxis
, "%15.12g %15.12g %15.12g ", residuehelixaxis
[i
][0], residuehelixaxis
[i
][1], residuehelixaxis
[i
][2]);
368 fprintf(fpcenter
, "%15.12g %15.12g %15.12g ", residueorigin
[i
][0], residueorigin
[i
][1], residueorigin
[i
][2]);
370 fprintf(fprise
, "%15.12g ", residuerise
[i
]);
371 fprintf(fpradius
, "%15.12g ", residueradius
[i
]);
372 fprintf(fptwist
, "%15.12g ", residuetwist
[i
]);
373 fprintf(fpbending
, "%15.12g ", residuebending
[i
]);
375 /* angle with local vector? */
377 printf("res[%2d]: axis: %g %g %g origin: %g %g %g vector: %g %g %g angle: %g\n",i,
378 residuehelixaxis[i][0],
379 residuehelixaxis[i][1],
380 residuehelixaxis[i][2],
387 180.0/M_PI*acos( cos_angle(residuevector[i],residuehelixaxis[i]) ));
389 /* fprintf(fp,"%15.12g %15.12g %15.12g %15.12g %15.12g %15.12g\n",
390 residuehelixaxis[i][0],
391 residuehelixaxis[i][1],
392 residuehelixaxis[i][2],
395 residuevector[i][2]);
399 fprintf(fprise
, "\n");
400 fprintf(fpradius
, "\n");
401 fprintf(fpaxis
, "\n");
402 fprintf(fpcenter
, "\n");
403 fprintf(fptwist
, "\n");
404 fprintf(fpbending
, "\n");
408 for (i
= 0; i
< iCA
; i
++)
410 copy_rvec(residuehelixaxis
[i
], residuehelixaxis_t0
[i
]);
411 copy_rvec(residuevector
[i
], residuevector_t0
[i
]);
412 copy_rvec(axis3
[i
], axis3_t0
[i
]);
417 fprintf(fptilt
, "%15.12g ", t
);
418 fprintf(fprotation
, "%15.12g ", t
);
419 fprintf(fptheta1
, "%15.12g ", t
);
420 fprintf(fptheta2
, "%15.12g ", t
);
421 fprintf(fptheta3
, "%15.12g ", t
);
423 for (i
= 0; i
< iCA
; i
++)
425 if (i
== 0 || i
== iCA
-1)
433 /* Total rotation & tilt */
434 copy_rvec(residuehelixaxis_t0
[i
], refaxes
[0]);
435 copy_rvec(residuevector_t0
[i
], refaxes
[1]);
436 copy_rvec(axis3_t0
[i
], refaxes
[2]);
440 /* Rotation/tilt since last step */
441 copy_rvec(residuehelixaxis_tlast
[i
], refaxes
[0]);
442 copy_rvec(residuevector_tlast
[i
], refaxes
[1]);
443 copy_rvec(axis3_tlast
[i
], refaxes
[2]);
445 copy_rvec(residuehelixaxis
[i
], newaxes
[0]);
446 copy_rvec(residuevector
[i
], newaxes
[1]);
447 copy_rvec(axis3
[i
], newaxes
[2]);
450 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",
452 refaxes[0][0],refaxes[0][1],refaxes[0][2],
453 refaxes[1][0],refaxes[1][1],refaxes[1][2],
454 refaxes[2][0],refaxes[2][1],refaxes[2][2],
455 newaxes[0][0],newaxes[0][1],newaxes[0][2],
456 newaxes[1][0],newaxes[1][1],newaxes[1][2],
457 newaxes[2][0],newaxes[2][1],newaxes[2][2]);
460 /* rotate reference frame onto unit axes */
461 calc_fit_R(3, 3, weight
, unitaxes
, refaxes
, A
);
462 for (j
= 0; j
< 3; j
++)
464 mvmul(A
, refaxes
[j
], rot_refaxes
[j
]);
465 mvmul(A
, newaxes
[j
], rot_newaxes
[j
]);
468 /* Determine local rotation matrix A */
469 calc_fit_R(3, 3, weight
, rot_newaxes
, rot_refaxes
, A
);
470 /* Calculate euler angles, from rotation order y-z-x, where
471 * x is helixaxis, y residuevector, and z axis3.
473 * A contains rotation column vectors.
477 printf("frame %d, i=%d, A: %g %g %g , %g %g %g , %g %g %g\n",
478 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]);
481 theta1
= 180.0/M_PI
*atan2(A
[0][2], A
[0][0]);
482 theta2
= 180.0/M_PI
*asin(-A
[0][1]);
483 theta3
= 180.0/M_PI
*atan2(A
[2][1], A
[1][1]);
485 tilt
= sqrt(theta1
*theta1
+theta2
*theta2
);
487 fprintf(fptheta1
, "%15.12g ", theta1
);
488 fprintf(fptheta2
, "%15.12g ", theta2
);
489 fprintf(fptheta3
, "%15.12g ", theta3
);
492 fprintf(fptilt
, "%15.12g ", tilt
);
493 fprintf(fprotation
, "%15.12g ", rotation
);
495 fprintf(fptilt
, "\n");
496 fprintf(fprotation
, "\n");
497 fprintf(fptheta1
, "\n");
498 fprintf(fptheta2
, "\n");
499 fprintf(fptheta3
, "\n");
502 for (i
= 0; i
< iCA
; i
++)
504 copy_rvec(residuehelixaxis
[i
], residuehelixaxis_tlast
[i
]);
505 copy_rvec(residuevector
[i
], residuevector_tlast
[i
]);
506 copy_rvec(axis3
[i
], axis3_tlast
[i
]);
511 while (read_next_x(oenv
, status
, &t
, x
, box
));
513 gmx_rmpbc_done(gpbc
);
516 gmx_ffclose(fpcenter
);
518 gmx_ffclose(fprotation
);
520 gmx_ffclose(fpradius
);
521 gmx_ffclose(fptwist
);
522 gmx_ffclose(fpbending
);
523 gmx_ffclose(fptheta1
);
524 gmx_ffclose(fptheta2
);
525 gmx_ffclose(fptheta3
);