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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
53 #include "gmx_fatal.h"
61 int gmx_helixorient(int argc,char *argv[])
63 const char *desc[] = {
64 "[TT]g_helixorient[tt] 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]"
88 real theta1,theta2,theta3;
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)
170 CopyRight(stderr,argv[0]);
172 parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,
173 NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
175 top=read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);
180 /* read index files */
181 printf("Select a group of Calpha atoms corresponding to a single continuous helix:\n");
182 get_index(&(top->atoms),ftp2fn_null(efNDX,NFILE,fnm),1,&iCA,&ind_CA,&gn_CA);
184 snew(x_SC,iCA); /* sic! */
191 snew(helixaxis,iCA-3);
193 snew(residuetwist,iCA);
195 snew(residueradius,iCA);
197 snew(residuerise,iCA);
198 snew(residueorigin,iCA);
199 snew(residuehelixaxis,iCA);
200 snew(residuevector,iCA);
201 snew(sidechainvector,iCA);
202 snew(residuebending,iCA);
203 snew(residuehelixaxis_t0,iCA);
204 snew(residuevector_t0,iCA);
206 snew(residuehelixaxis_tlast,iCA);
207 snew(residuevector_tlast,iCA);
208 snew(axis3_tlast,iCA);
213 printf("Select a group of atoms defining the sidechain direction (1/residue):\n");
214 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);
220 natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
222 fpaxis=ffopen(opt2fn("-oaxis",NFILE,fnm),"w");
223 fpcenter=ffopen(opt2fn("-ocenter",NFILE,fnm),"w");
224 fprise=ffopen(opt2fn("-orise",NFILE,fnm),"w");
225 fpradius=ffopen(opt2fn("-oradius",NFILE,fnm),"w");
226 fptwist=ffopen(opt2fn("-otwist",NFILE,fnm),"w");
227 fpbending=ffopen(opt2fn("-obending",NFILE,fnm),"w");
229 fptheta1=ffopen("theta1.xvg","w");
230 fptheta2=ffopen("theta2.xvg","w");
231 fptheta3=ffopen("theta3.xvg","w");
235 fptilt=xvgropen(opt2fn("-otilt",NFILE,fnm),
236 "Incremental local helix tilt","Time(ps)","Tilt (degrees)",
238 fprotation=xvgropen(opt2fn("-orot",NFILE,fnm),
239 "Incremental local helix rotation","Time(ps)",
240 "Rotation (degrees)",oenv);
244 fptilt=xvgropen(opt2fn("-otilt",NFILE,fnm),
245 "Cumulative local helix tilt","Time(ps)","Tilt (degrees)",oenv);
246 fprotation=xvgropen(opt2fn("-orot",NFILE,fnm),
247 "Cumulative local helix rotation","Time(ps)",
248 "Rotation (degrees)",oenv);
251 clear_rvecs(3,unitaxes);
256 gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
260 /* initialisation for correct distance calculations */
261 set_pbc(&pbc,ePBC,box);
262 /* make molecules whole again */
263 gmx_rmpbc(gpbc,natoms,box,x);
265 /* copy coords to our smaller arrays */
268 copy_rvec(x[ind_CA[i]],x_CA[i]);
271 copy_rvec(x[ind_SC[i]],x_SC[i]);
277 rvec_sub(x_CA[i+1],x_CA[i],r12[i]);
278 rvec_sub(x_CA[i+2],x_CA[i+1],r23[i]);
279 rvec_sub(x_CA[i+3],x_CA[i+2],r34[i]);
280 rvec_sub(r12[i],r23[i],diff13[i]);
281 rvec_sub(r23[i],r34[i],diff24[i]);
282 /* calculate helix axis */
283 cprod(diff13[i],diff24[i],helixaxis[i]);
284 svmul(1.0/norm(helixaxis[i]),helixaxis[i],helixaxis[i]);
286 tmp = cos_angle(diff13[i],diff24[i]);
287 twist[i] = 180.0/M_PI * acos( tmp );
288 radius[i] = sqrt( norm(diff13[i])*norm(diff24[i]) ) / (2.0* (1.0-tmp) );
289 rise[i]=fabs(iprod(r23[i],helixaxis[i]));
291 svmul(radius[i]/norm(diff13[i]),diff13[i],v1);
292 svmul(radius[i]/norm(diff24[i]),diff24[i],v2);
294 rvec_sub(x_CA[i+1],v1,residueorigin[i+1]);
295 rvec_sub(x_CA[i+2],v2,residueorigin[i+2]);
297 residueradius[0]=residuetwist[0]=residuerise[0]=0;
299 residueradius[1]=radius[0];
300 residuetwist[1]=twist[0];
301 residuerise[1]=rise[0];
303 residuebending[0]=residuebending[1]=0;
306 residueradius[i]=0.5*(radius[i-2]+radius[i-1]);
307 residuetwist[i]=0.5*(twist[i-2]+twist[i-1]);
308 residuerise[i]=0.5*(rise[i-2]+rise[i-1]);
309 residuebending[i] = 180.0/M_PI*acos( cos_angle(helixaxis[i-2],helixaxis[i-1]) );
311 residueradius[iCA-2]=radius[iCA-4];
312 residuetwist[iCA-2]=twist[iCA-4];
313 residuerise[iCA-2]=rise[iCA-4];
314 residueradius[iCA-1]=residuetwist[iCA-1]=residuerise[iCA-1]=0;
315 residuebending[iCA-2]=residuebending[iCA-1]=0;
317 clear_rvec(residueorigin[0]);
318 clear_rvec(residueorigin[iCA-1]);
320 /* average helix axes to define them on the residues.
321 * Just extrapolate second first/list atom.
323 copy_rvec(helixaxis[0],residuehelixaxis[0]);
324 copy_rvec(helixaxis[0],residuehelixaxis[1]);
328 rvec_add(helixaxis[i-2],helixaxis[i-1],residuehelixaxis[i]);
329 svmul(0.5,residuehelixaxis[i],residuehelixaxis[i]);
331 copy_rvec(helixaxis[iCA-4],residuehelixaxis[iCA-2]);
332 copy_rvec(helixaxis[iCA-4],residuehelixaxis[iCA-1]);
334 /* Normalize the axis */
337 svmul(1.0/norm(residuehelixaxis[i]),residuehelixaxis[i],residuehelixaxis[i]);
340 /* calculate vector from origin to residue CA */
341 fprintf(fpaxis,"%15.12g ",t);
342 fprintf(fpcenter,"%15.12g ",t);
343 fprintf(fprise,"%15.12g ",t);
344 fprintf(fpradius,"%15.12g ",t);
345 fprintf(fptwist,"%15.12g ",t);
346 fprintf(fpbending,"%15.12g ",t);
352 fprintf(fpaxis,"%15.12g %15.12g %15.12g ",0.0,0.0,0.0);
353 fprintf(fpcenter,"%15.12g %15.12g %15.12g ",0.0,0.0,0.0);
354 fprintf(fprise,"%15.12g ",0.0);
355 fprintf(fpradius,"%15.12g ",0.0);
356 fprintf(fptwist,"%15.12g ",0.0);
357 fprintf(fpbending,"%15.12g ",0.0);
361 rvec_sub( bSC ? x_SC[i] : x_CA[i] ,residueorigin[i], residuevector[i]);
362 svmul(1.0/norm(residuevector[i]),residuevector[i],residuevector[i]);
363 cprod(residuehelixaxis[i],residuevector[i],axis3[i]);
364 fprintf(fpaxis,"%15.12g %15.12g %15.12g ",residuehelixaxis[i][0],residuehelixaxis[i][1],residuehelixaxis[i][2]);
365 fprintf(fpcenter,"%15.12g %15.12g %15.12g ",residueorigin[i][0],residueorigin[i][1],residueorigin[i][2]);
367 fprintf(fprise,"%15.12g ",residuerise[i]);
368 fprintf(fpradius,"%15.12g ",residueradius[i]);
369 fprintf(fptwist,"%15.12g ",residuetwist[i]);
370 fprintf(fpbending,"%15.12g ",residuebending[i]);
372 /* angle with local vector? */
374 printf("res[%2d]: axis: %g %g %g origin: %g %g %g vector: %g %g %g angle: %g\n",i,
375 residuehelixaxis[i][0],
376 residuehelixaxis[i][1],
377 residuehelixaxis[i][2],
384 180.0/M_PI*acos( cos_angle(residuevector[i],residuehelixaxis[i]) ));
386 /* fprintf(fp,"%15.12g %15.12g %15.12g %15.12g %15.12g %15.12g\n",
387 residuehelixaxis[i][0],
388 residuehelixaxis[i][1],
389 residuehelixaxis[i][2],
392 residuevector[i][2]);
396 fprintf(fprise,"\n");
397 fprintf(fpradius,"\n");
398 fprintf(fpaxis,"\n");
399 fprintf(fpcenter,"\n");
400 fprintf(fptwist,"\n");
401 fprintf(fpbending,"\n");
407 copy_rvec(residuehelixaxis[i],residuehelixaxis_t0[i]);
408 copy_rvec(residuevector[i],residuevector_t0[i]);
409 copy_rvec(axis3[i],axis3_t0[i]);
414 fprintf(fptilt,"%15.12g ",t);
415 fprintf(fprotation,"%15.12g ",t);
416 fprintf(fptheta1,"%15.12g ",t);
417 fprintf(fptheta2,"%15.12g ",t);
418 fprintf(fptheta3,"%15.12g ",t);
430 /* Total rotation & tilt */
431 copy_rvec(residuehelixaxis_t0[i],refaxes[0]);
432 copy_rvec(residuevector_t0[i],refaxes[1]);
433 copy_rvec(axis3_t0[i],refaxes[2]);
437 /* Rotation/tilt since last step */
438 copy_rvec(residuehelixaxis_tlast[i],refaxes[0]);
439 copy_rvec(residuevector_tlast[i],refaxes[1]);
440 copy_rvec(axis3_tlast[i],refaxes[2]);
442 copy_rvec(residuehelixaxis[i],newaxes[0]);
443 copy_rvec(residuevector[i],newaxes[1]);
444 copy_rvec(axis3[i],newaxes[2]);
447 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",
449 refaxes[0][0],refaxes[0][1],refaxes[0][2],
450 refaxes[1][0],refaxes[1][1],refaxes[1][2],
451 refaxes[2][0],refaxes[2][1],refaxes[2][2],
452 newaxes[0][0],newaxes[0][1],newaxes[0][2],
453 newaxes[1][0],newaxes[1][1],newaxes[1][2],
454 newaxes[2][0],newaxes[2][1],newaxes[2][2]);
457 /* rotate reference frame onto unit axes */
458 calc_fit_R(3,3,weight,unitaxes,refaxes,A);
461 mvmul(A,refaxes[j],rot_refaxes[j]);
462 mvmul(A,newaxes[j],rot_newaxes[j]);
465 /* Determine local rotation matrix A */
466 calc_fit_R(3,3,weight,rot_newaxes,rot_refaxes,A);
467 /* Calculate euler angles, from rotation order y-z-x, where
468 * x is helixaxis, y residuevector, and z axis3.
470 * A contains rotation column vectors.
474 printf("frame %d, i=%d, A: %g %g %g , %g %g %g , %g %g %g\n",
475 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]);
478 theta1 = 180.0/M_PI*atan2(A[0][2],A[0][0]);
479 theta2 = 180.0/M_PI*asin(-A[0][1]);
480 theta3 = 180.0/M_PI*atan2(A[2][1],A[1][1]);
482 tilt = sqrt(theta1*theta1+theta2*theta2);
484 fprintf(fptheta1,"%15.12g ",theta1);
485 fprintf(fptheta2,"%15.12g ",theta2);
486 fprintf(fptheta3,"%15.12g ",theta3);
489 fprintf(fptilt,"%15.12g ",tilt);
490 fprintf(fprotation,"%15.12g ",rotation);
492 fprintf(fptilt,"\n");
493 fprintf(fprotation,"\n");
494 fprintf(fptheta1,"\n");
495 fprintf(fptheta2,"\n");
496 fprintf(fptheta3,"\n");
501 copy_rvec(residuehelixaxis[i],residuehelixaxis_tlast[i]);
502 copy_rvec(residuevector[i],residuevector_tlast[i]);
503 copy_rvec(axis3[i],axis3_tlast[i]);
507 } while (read_next_x(oenv,status,&t,natoms,x,box));
509 gmx_rmpbc_done(gpbc);