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]"
85 real theta1
,theta2
,theta3
;
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)
167 CopyRight(stderr
,argv
[0]);
169 parse_common_args(&argc
,argv
,PCA_CAN_TIME
| PCA_BE_NICE
,
170 NFILE
,fnm
,NPA
,pa
,asize(desc
),desc
,0,NULL
,&oenv
);
172 top
=read_top(ftp2fn(efTPX
,NFILE
,fnm
),&ePBC
);
177 /* read index files */
178 printf("Select a group of Calpha atoms corresponding to a single continuous helix:\n");
179 get_index(&(top
->atoms
),ftp2fn_null(efNDX
,NFILE
,fnm
),1,&iCA
,&ind_CA
,&gn_CA
);
181 snew(x_SC
,iCA
); /* sic! */
188 snew(helixaxis
,iCA
-3);
190 snew(residuetwist
,iCA
);
192 snew(residueradius
,iCA
);
194 snew(residuerise
,iCA
);
195 snew(residueorigin
,iCA
);
196 snew(residuehelixaxis
,iCA
);
197 snew(residuevector
,iCA
);
198 snew(sidechainvector
,iCA
);
199 snew(residuebending
,iCA
);
200 snew(residuehelixaxis_t0
,iCA
);
201 snew(residuevector_t0
,iCA
);
203 snew(residuehelixaxis_tlast
,iCA
);
204 snew(residuevector_tlast
,iCA
);
205 snew(axis3_tlast
,iCA
);
210 printf("Select a group of atoms defining the sidechain direction (1/residue):\n");
211 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
);
217 natoms
=read_first_x(oenv
,&status
,ftp2fn(efTRX
,NFILE
,fnm
),&t
,&x
,box
);
219 fpaxis
=ffopen(opt2fn("-oaxis",NFILE
,fnm
),"w");
220 fpcenter
=ffopen(opt2fn("-ocenter",NFILE
,fnm
),"w");
221 fprise
=ffopen(opt2fn("-orise",NFILE
,fnm
),"w");
222 fpradius
=ffopen(opt2fn("-oradius",NFILE
,fnm
),"w");
223 fptwist
=ffopen(opt2fn("-otwist",NFILE
,fnm
),"w");
224 fpbending
=ffopen(opt2fn("-obending",NFILE
,fnm
),"w");
226 fptheta1
=ffopen("theta1.xvg","w");
227 fptheta2
=ffopen("theta2.xvg","w");
228 fptheta3
=ffopen("theta3.xvg","w");
232 fptilt
=xvgropen(opt2fn("-otilt",NFILE
,fnm
),
233 "Incremental local helix tilt","Time(ps)","Tilt (degrees)",
235 fprotation
=xvgropen(opt2fn("-orot",NFILE
,fnm
),
236 "Incremental local helix rotation","Time(ps)",
237 "Rotation (degrees)",oenv
);
241 fptilt
=xvgropen(opt2fn("-otilt",NFILE
,fnm
),
242 "Cumulative local helix tilt","Time(ps)","Tilt (degrees)",oenv
);
243 fprotation
=xvgropen(opt2fn("-orot",NFILE
,fnm
),
244 "Cumulative local helix rotation","Time(ps)",
245 "Rotation (degrees)",oenv
);
248 clear_rvecs(3,unitaxes
);
253 gpbc
= gmx_rmpbc_init(&top
->idef
,ePBC
,natoms
,box
);
257 /* initialisation for correct distance calculations */
258 set_pbc(&pbc
,ePBC
,box
);
259 /* make molecules whole again */
260 gmx_rmpbc(gpbc
,natoms
,box
,x
);
262 /* copy coords to our smaller arrays */
265 copy_rvec(x
[ind_CA
[i
]],x_CA
[i
]);
268 copy_rvec(x
[ind_SC
[i
]],x_SC
[i
]);
274 rvec_sub(x_CA
[i
+1],x_CA
[i
],r12
[i
]);
275 rvec_sub(x_CA
[i
+2],x_CA
[i
+1],r23
[i
]);
276 rvec_sub(x_CA
[i
+3],x_CA
[i
+2],r34
[i
]);
277 rvec_sub(r12
[i
],r23
[i
],diff13
[i
]);
278 rvec_sub(r23
[i
],r34
[i
],diff24
[i
]);
279 /* calculate helix axis */
280 cprod(diff13
[i
],diff24
[i
],helixaxis
[i
]);
281 svmul(1.0/norm(helixaxis
[i
]),helixaxis
[i
],helixaxis
[i
]);
283 tmp
= cos_angle(diff13
[i
],diff24
[i
]);
284 twist
[i
] = 180.0/M_PI
* acos( tmp
);
285 radius
[i
] = sqrt( norm(diff13
[i
])*norm(diff24
[i
]) ) / (2.0* (1.0-tmp
) );
286 rise
[i
]=fabs(iprod(r23
[i
],helixaxis
[i
]));
288 svmul(radius
[i
]/norm(diff13
[i
]),diff13
[i
],v1
);
289 svmul(radius
[i
]/norm(diff24
[i
]),diff24
[i
],v2
);
291 rvec_sub(x_CA
[i
+1],v1
,residueorigin
[i
+1]);
292 rvec_sub(x_CA
[i
+2],v2
,residueorigin
[i
+2]);
294 residueradius
[0]=residuetwist
[0]=residuerise
[0]=0;
296 residueradius
[1]=radius
[0];
297 residuetwist
[1]=twist
[0];
298 residuerise
[1]=rise
[0];
300 residuebending
[0]=residuebending
[1]=0;
303 residueradius
[i
]=0.5*(radius
[i
-2]+radius
[i
-1]);
304 residuetwist
[i
]=0.5*(twist
[i
-2]+twist
[i
-1]);
305 residuerise
[i
]=0.5*(rise
[i
-2]+rise
[i
-1]);
306 residuebending
[i
] = 180.0/M_PI
*acos( cos_angle(helixaxis
[i
-2],helixaxis
[i
-1]) );
308 residueradius
[iCA
-2]=radius
[iCA
-4];
309 residuetwist
[iCA
-2]=twist
[iCA
-4];
310 residuerise
[iCA
-2]=rise
[iCA
-4];
311 residueradius
[iCA
-1]=residuetwist
[iCA
-1]=residuerise
[iCA
-1]=0;
312 residuebending
[iCA
-2]=residuebending
[iCA
-1]=0;
314 clear_rvec(residueorigin
[0]);
315 clear_rvec(residueorigin
[iCA
-1]);
317 /* average helix axes to define them on the residues.
318 * Just extrapolate second first/list atom.
320 copy_rvec(helixaxis
[0],residuehelixaxis
[0]);
321 copy_rvec(helixaxis
[0],residuehelixaxis
[1]);
325 rvec_add(helixaxis
[i
-2],helixaxis
[i
-1],residuehelixaxis
[i
]);
326 svmul(0.5,residuehelixaxis
[i
],residuehelixaxis
[i
]);
328 copy_rvec(helixaxis
[iCA
-4],residuehelixaxis
[iCA
-2]);
329 copy_rvec(helixaxis
[iCA
-4],residuehelixaxis
[iCA
-1]);
331 /* Normalize the axis */
334 svmul(1.0/norm(residuehelixaxis
[i
]),residuehelixaxis
[i
],residuehelixaxis
[i
]);
337 /* calculate vector from origin to residue CA */
338 fprintf(fpaxis
,"%15.12g ",t
);
339 fprintf(fpcenter
,"%15.12g ",t
);
340 fprintf(fprise
,"%15.12g ",t
);
341 fprintf(fpradius
,"%15.12g ",t
);
342 fprintf(fptwist
,"%15.12g ",t
);
343 fprintf(fpbending
,"%15.12g ",t
);
349 fprintf(fpaxis
,"%15.12g %15.12g %15.12g ",0.0,0.0,0.0);
350 fprintf(fpcenter
,"%15.12g %15.12g %15.12g ",0.0,0.0,0.0);
351 fprintf(fprise
,"%15.12g ",0.0);
352 fprintf(fpradius
,"%15.12g ",0.0);
353 fprintf(fptwist
,"%15.12g ",0.0);
354 fprintf(fpbending
,"%15.12g ",0.0);
358 rvec_sub( bSC
? x_SC
[i
] : x_CA
[i
] ,residueorigin
[i
], residuevector
[i
]);
359 svmul(1.0/norm(residuevector
[i
]),residuevector
[i
],residuevector
[i
]);
360 cprod(residuehelixaxis
[i
],residuevector
[i
],axis3
[i
]);
361 fprintf(fpaxis
,"%15.12g %15.12g %15.12g ",residuehelixaxis
[i
][0],residuehelixaxis
[i
][1],residuehelixaxis
[i
][2]);
362 fprintf(fpcenter
,"%15.12g %15.12g %15.12g ",residueorigin
[i
][0],residueorigin
[i
][1],residueorigin
[i
][2]);
364 fprintf(fprise
,"%15.12g ",residuerise
[i
]);
365 fprintf(fpradius
,"%15.12g ",residueradius
[i
]);
366 fprintf(fptwist
,"%15.12g ",residuetwist
[i
]);
367 fprintf(fpbending
,"%15.12g ",residuebending
[i
]);
369 /* angle with local vector? */
371 printf("res[%2d]: axis: %g %g %g origin: %g %g %g vector: %g %g %g angle: %g\n",i,
372 residuehelixaxis[i][0],
373 residuehelixaxis[i][1],
374 residuehelixaxis[i][2],
381 180.0/M_PI*acos( cos_angle(residuevector[i],residuehelixaxis[i]) ));
383 /* fprintf(fp,"%15.12g %15.12g %15.12g %15.12g %15.12g %15.12g\n",
384 residuehelixaxis[i][0],
385 residuehelixaxis[i][1],
386 residuehelixaxis[i][2],
389 residuevector[i][2]);
393 fprintf(fprise
,"\n");
394 fprintf(fpradius
,"\n");
395 fprintf(fpaxis
,"\n");
396 fprintf(fpcenter
,"\n");
397 fprintf(fptwist
,"\n");
398 fprintf(fpbending
,"\n");
404 copy_rvec(residuehelixaxis
[i
],residuehelixaxis_t0
[i
]);
405 copy_rvec(residuevector
[i
],residuevector_t0
[i
]);
406 copy_rvec(axis3
[i
],axis3_t0
[i
]);
411 fprintf(fptilt
,"%15.12g ",t
);
412 fprintf(fprotation
,"%15.12g ",t
);
413 fprintf(fptheta1
,"%15.12g ",t
);
414 fprintf(fptheta2
,"%15.12g ",t
);
415 fprintf(fptheta3
,"%15.12g ",t
);
427 /* Total rotation & tilt */
428 copy_rvec(residuehelixaxis_t0
[i
],refaxes
[0]);
429 copy_rvec(residuevector_t0
[i
],refaxes
[1]);
430 copy_rvec(axis3_t0
[i
],refaxes
[2]);
434 /* Rotation/tilt since last step */
435 copy_rvec(residuehelixaxis_tlast
[i
],refaxes
[0]);
436 copy_rvec(residuevector_tlast
[i
],refaxes
[1]);
437 copy_rvec(axis3_tlast
[i
],refaxes
[2]);
439 copy_rvec(residuehelixaxis
[i
],newaxes
[0]);
440 copy_rvec(residuevector
[i
],newaxes
[1]);
441 copy_rvec(axis3
[i
],newaxes
[2]);
444 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",
446 refaxes[0][0],refaxes[0][1],refaxes[0][2],
447 refaxes[1][0],refaxes[1][1],refaxes[1][2],
448 refaxes[2][0],refaxes[2][1],refaxes[2][2],
449 newaxes[0][0],newaxes[0][1],newaxes[0][2],
450 newaxes[1][0],newaxes[1][1],newaxes[1][2],
451 newaxes[2][0],newaxes[2][1],newaxes[2][2]);
454 /* rotate reference frame onto unit axes */
455 calc_fit_R(3,3,weight
,unitaxes
,refaxes
,A
);
458 mvmul(A
,refaxes
[j
],rot_refaxes
[j
]);
459 mvmul(A
,newaxes
[j
],rot_newaxes
[j
]);
462 /* Determine local rotation matrix A */
463 calc_fit_R(3,3,weight
,rot_newaxes
,rot_refaxes
,A
);
464 /* Calculate euler angles, from rotation order y-z-x, where
465 * x is helixaxis, y residuevector, and z axis3.
467 * A contains rotation column vectors.
471 printf("frame %d, i=%d, A: %g %g %g , %g %g %g , %g %g %g\n",
472 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]);
475 theta1
= 180.0/M_PI
*atan2(A
[0][2],A
[0][0]);
476 theta2
= 180.0/M_PI
*asin(-A
[0][1]);
477 theta3
= 180.0/M_PI
*atan2(A
[2][1],A
[1][1]);
479 tilt
= sqrt(theta1
*theta1
+theta2
*theta2
);
481 fprintf(fptheta1
,"%15.12g ",theta1
);
482 fprintf(fptheta2
,"%15.12g ",theta2
);
483 fprintf(fptheta3
,"%15.12g ",theta3
);
486 fprintf(fptilt
,"%15.12g ",tilt
);
487 fprintf(fprotation
,"%15.12g ",rotation
);
489 fprintf(fptilt
,"\n");
490 fprintf(fprotation
,"\n");
491 fprintf(fptheta1
,"\n");
492 fprintf(fptheta2
,"\n");
493 fprintf(fptheta3
,"\n");
498 copy_rvec(residuehelixaxis
[i
],residuehelixaxis_tlast
[i
]);
499 copy_rvec(residuevector
[i
],residuevector_tlast
[i
]);
500 copy_rvec(axis3
[i
],axis3_tlast
[i
]);
504 } while (read_next_x(oenv
,status
,&t
,natoms
,x
,box
));
506 gmx_rmpbc_done(gpbc
);