src/tools/gmx_helixorient.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  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.
  14
  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.
  19  *
  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.
  26  *
  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.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38 #include <typedefs.h>
  39
  40 #include "smalloc.h"
  41 #include "macros.h"
  42 #include "math.h"
  43 #include "xvgr.h"
  44 #include "copyrite.h"
  45 #include "statutil.h"
  46 #include "string2.h"
  47 #include "vec.h"
  48 #include "index.h"
  49 #include "pbc.h"
  50 #include "gmx_fatal.h"
  51 #include "futil.h"
  52 #include "gstat.h"
  53 #include "pbc.h"
  54 #include "do_fit.h"
  55 #include "gmx_ana.h"
  56
  57
  58 int gmx_helixorient(int argc,char *argv[])
  59 {
  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",
  70     "axis.[PAR]",
  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]"
  77     };
  78
  79     t_topology *top=NULL;
  80     real t;
  81     rvec *x=NULL,dx;
  82     matrix box;
  83     t_trxstatus *status;
  84     int natoms;
  85     real theta1,theta2,theta3;
  86
  87     int d,i,j,teller=0;
  88     int iCA,iSC;
  89     atom_id *ind_CA;
  90     atom_id *ind_SC;
  91     char *gn_CA;
  92     char *gn_SC;
  93     rvec averageaxis;
  94     rvec v1,v2,p1,p2,vtmp,vproj;
  95     rvec *x_CA,*x_SC;
  96     rvec *r12;
  97     rvec *r23;
  98     rvec *r34;
  99     rvec *diff13;
 100     rvec *diff24;
 101     rvec *helixaxis;
 102     rvec *residuehelixaxis;
 103     rvec *residueorigin;
 104     rvec *residuevector;
 105     rvec *sidechainvector;
 106
 107     rvec axes_t0[3];
 108     rvec axes[3];
 109     rvec *residuehelixaxis_t0;
 110     rvec *residuevector_t0;
 111     rvec *axis3_t0;
 112     rvec *residuehelixaxis_tlast;
 113     rvec *residuevector_tlast;
 114     rvec *axis3_tlast;
 115     rvec refaxes[3],newaxes[3];
 116     rvec unitaxes[3];
 117     rvec rot_refaxes[3],rot_newaxes[3];
 118
 119     real tilt,rotation;
 120     rvec *axis3;
 121     real *twist,*residuetwist;
 122     real *radius,*residueradius;
 123     real *rise,*residuerise;
 124     real *residuebending;
 125
 126     real tmp,rotangle;
 127     real weight[3];
 128     t_pbc   pbc;
 129     matrix A;
 130
 131     FILE *fpaxis,*fpcenter,*fptilt,*fprotation;
 132     FILE *fpradius,*fprise,*fptwist;
 133     FILE *fptheta1,*fptheta2,*fptheta3;
 134     FILE *fpbending;
 135     int ePBC;
 136
 137     output_env_t oenv;
 138     gmx_rmpbc_t  gpbc=NULL;
 139
 140     static  gmx_bool bSC=FALSE;
 141     static gmx_bool bIncremental = FALSE;
 142
 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." },
 148     };
 149 #define NPA asize(pa)
 150
 151   t_filenm fnm[] = {
 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 }
 163   };
 164 #define NFILE asize(fnm)
 165
 166
 167   CopyRight(stderr,argv[0]);
 168
 169   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,
 170                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
 171
 172   top=read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);
 173
 174     for(i=0;i<3;i++)
 175         weight[i]=1.0;
 176
 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);
 180     snew(x_CA,iCA);
 181     snew(x_SC,iCA); /* sic! */
 182
 183     snew(r12,iCA-3);
 184     snew(r23,iCA-3);
 185     snew(r34,iCA-3);
 186     snew(diff13,iCA-3);
 187     snew(diff24,iCA-3);
 188     snew(helixaxis,iCA-3);
 189     snew(twist,iCA);
 190     snew(residuetwist,iCA);
 191     snew(radius,iCA);
 192     snew(residueradius,iCA);
 193     snew(rise,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);
 202     snew(axis3_t0,iCA);
 203     snew(residuehelixaxis_tlast,iCA);
 204     snew(residuevector_tlast,iCA);
 205     snew(axis3_tlast,iCA);
 206     snew(axis3,iCA);
 207
 208     if(bSC)
 209     {
 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);
 212         if(iSC!=iCA)
 213             gmx_fatal(FARGS,"Number of sidechain atoms (%d) != number of CA atoms (%d)",iSC,iCA);
 214
 215     }
 216
 217     natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 218
 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");
 225
 226     fptheta1=ffopen("theta1.xvg","w");
 227     fptheta2=ffopen("theta2.xvg","w");
 228     fptheta3=ffopen("theta3.xvg","w");
 229
 230     if(bIncremental)
 231     {
 232         fptilt=xvgropen(opt2fn("-otilt",NFILE,fnm),
 233                 "Incremental local helix tilt","Time(ps)","Tilt (degrees)",
 234                 oenv);
 235         fprotation=xvgropen(opt2fn("-orot",NFILE,fnm),
 236                 "Incremental local helix rotation","Time(ps)",
 237                 "Rotation (degrees)",oenv);
 238     }
 239     else
 240     {
 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);
 246     }
 247
 248     clear_rvecs(3,unitaxes);
 249     unitaxes[0][0]=1;
 250     unitaxes[1][1]=1;
 251     unitaxes[2][2]=1;
 252
 253     gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
 254
 255   do
 256   {
 257       /* initialisation for correct distance calculations */
 258     set_pbc(&pbc,ePBC,box);
 259       /* make molecules whole again */
 260     gmx_rmpbc(gpbc,natoms,box,x);
 261
 262       /* copy coords to our smaller arrays */
 263       for(i=0;i<iCA;i++)
 264       {
 265           copy_rvec(x[ind_CA[i]],x_CA[i]);
 266           if(bSC)
 267           {
 268               copy_rvec(x[ind_SC[i]],x_SC[i]);
 269           }
 270       }
 271
 272       for(i=0;i<iCA-3;i++)
 273       {
 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]);
 282
 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]));
 287
 288           svmul(radius[i]/norm(diff13[i]),diff13[i],v1);
 289           svmul(radius[i]/norm(diff24[i]),diff24[i],v2);
 290
 291           rvec_sub(x_CA[i+1],v1,residueorigin[i+1]);
 292           rvec_sub(x_CA[i+2],v2,residueorigin[i+2]);
 293       }
 294       residueradius[0]=residuetwist[0]=residuerise[0]=0;
 295
 296       residueradius[1]=radius[0];
 297       residuetwist[1]=twist[0];
 298       residuerise[1]=rise[0];
 299
 300       residuebending[0]=residuebending[1]=0;
 301       for(i=2;i<iCA-2;i++)
 302         {
 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]) );
 307         }
 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;
 313
 314       clear_rvec(residueorigin[0]);
 315       clear_rvec(residueorigin[iCA-1]);
 316
 317       /* average helix axes to define them on the residues.
 318        * Just extrapolate second first/list atom.
 319        */
 320       copy_rvec(helixaxis[0],residuehelixaxis[0]);
 321       copy_rvec(helixaxis[0],residuehelixaxis[1]);
 322
 323       for(i=2;i<iCA-2;i++)
 324       {
 325           rvec_add(helixaxis[i-2],helixaxis[i-1],residuehelixaxis[i]);
 326           svmul(0.5,residuehelixaxis[i],residuehelixaxis[i]);
 327       }
 328       copy_rvec(helixaxis[iCA-4],residuehelixaxis[iCA-2]);
 329       copy_rvec(helixaxis[iCA-4],residuehelixaxis[iCA-1]);
 330
 331       /* Normalize the axis */
 332       for(i=0;i<iCA;i++)
 333       {
 334           svmul(1.0/norm(residuehelixaxis[i]),residuehelixaxis[i],residuehelixaxis[i]);
 335       }
 336
 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);
 344
 345       for(i=0;i<iCA;i++)
 346       {
 347         if(i==0 || i==iCA-1)
 348           {
 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);
 355           }
 356         else
 357           {
 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]);
 363
 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]);
 368
 369           /* angle with local vector? */
 370           /*
 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],
 375                  residueorigin[i][0],
 376                  residueorigin[i][1],
 377                  residueorigin[i][2],
 378                  residuevector[i][0],
 379                  residuevector[i][1],
 380                  residuevector[i][2],
 381                  180.0/M_PI*acos( cos_angle(residuevector[i],residuehelixaxis[i]) ));
 382            */
 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],
 387                   residuevector[i][0],
 388                   residuevector[i][1],
 389                   residuevector[i][2]);
 390      */
 391           }
 392       }
 393       fprintf(fprise,"\n");
 394       fprintf(fpradius,"\n");
 395       fprintf(fpaxis,"\n");
 396       fprintf(fpcenter,"\n");
 397       fprintf(fptwist,"\n");
 398       fprintf(fpbending,"\n");
 399
 400       if(teller==0)
 401       {
 402           for(i=0;i<iCA;i++)
 403           {
 404               copy_rvec(residuehelixaxis[i],residuehelixaxis_t0[i]);
 405               copy_rvec(residuevector[i],residuevector_t0[i]);
 406               copy_rvec(axis3[i],axis3_t0[i]);
 407           }
 408       }
 409       else
 410       {
 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);
 416
 417           for(i=0;i<iCA;i++)
 418           {
 419             if(i==0 || i==iCA-1)
 420             {
 421                 tilt=rotation=0;
 422             }
 423             else
 424             {
 425               if(!bIncremental)
 426               {
 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]);
 431                }
 432               else
 433               {
 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]);
 438               }
 439               copy_rvec(residuehelixaxis[i],newaxes[0]);
 440               copy_rvec(residuevector[i],newaxes[1]);
 441               copy_rvec(axis3[i],newaxes[2]);
 442
 443               /*
 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",
 445                      teller,i,
 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]);
 452               */
 453
 454               /* rotate reference frame onto unit axes */
 455               calc_fit_R(3,3,weight,unitaxes,refaxes,A);
 456               for(j=0;j<3;j++)
 457                 {
 458                   mvmul(A,refaxes[j],rot_refaxes[j]);
 459                   mvmul(A,newaxes[j],rot_newaxes[j]);
 460                 }
 461
 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.
 466                *
 467                * A contains rotation column vectors.
 468                */
 469
 470               /*
 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]);
 473               */
 474
 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]);
 478
 479               tilt = sqrt(theta1*theta1+theta2*theta2);
 480               rotation = theta3;
 481               fprintf(fptheta1,"%15.12g  ",theta1);
 482               fprintf(fptheta2,"%15.12g  ",theta2);
 483               fprintf(fptheta3,"%15.12g  ",theta3);
 484
 485             }
 486             fprintf(fptilt,"%15.12g  ",tilt);
 487             fprintf(fprotation,"%15.12g  ",rotation);
 488           }
 489           fprintf(fptilt,"\n");
 490           fprintf(fprotation,"\n");
 491           fprintf(fptheta1,"\n");
 492           fprintf(fptheta2,"\n");
 493           fprintf(fptheta3,"\n");
 494       }
 495
 496       for(i=0;i<iCA;i++)
 497       {
 498           copy_rvec(residuehelixaxis[i],residuehelixaxis_tlast[i]);
 499           copy_rvec(residuevector[i],residuevector_tlast[i]);
 500           copy_rvec(axis3[i],axis3_tlast[i]);
 501       }
 502
 503       teller++;
 504   } while (read_next_x(oenv,status,&t,natoms,x,box));
 505
 506     gmx_rmpbc_done(gpbc);
 507
 508     ffclose(fpaxis);
 509     ffclose(fpcenter);
 510     ffclose(fptilt);
 511     ffclose(fprotation);
 512     ffclose(fprise);
 513     ffclose(fpradius);
 514     ffclose(fptwist);
 515     ffclose(fpbending);
 516     ffclose(fptheta1);
 517     ffclose(fptheta2);
 518     ffclose(fptheta3);
 519
 520     close_trj(status);
 521
 522     thanx(stderr);
 523     return 0;
 524 }