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     "g_helixorient calculates coordinates and direction of the average",
  62     "axis inside an alpha helix, and the direction/vectors of both the",
  63     "alpha carbon and (optionally) a sidechain atom relative to the axis.[PAR]",
  64     "As input, you need to specify an index group with alpha carbon atoms",
  65     "corresponding to an alpha 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     "Note that this program does not do any fitting of structures.[PAR]",
  69     "We need four Calpha 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 X axis, the residues/CA-vector as Y, and the",
  73     "Z 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 theta1-3.xvg"
  77     };
  78
  79     t_topology *top=NULL;
  80     real t;
  81     rvec *x=NULL,dx;
  82     matrix box;
  83     int 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
 139     static  bool bSC=FALSE;
 140     static bool bIncremental = FALSE;
 141
 142     static t_pargs pa[] = {
 143         { "-sidechain",      FALSE, etBOOL, {&bSC},
 144         "Calculate sidechain directions relative to helix axis too." },
 145         { "-incremental",        FALSE, etBOOL, {&bIncremental},
 146         "Calculate incremental rather than total rotation/tilt." },
 147     };
 148 #define NPA asize(pa)
 149
 150   t_filenm fnm[] = {
 151     { efTPX, NULL, NULL, ffREAD },
 152     { efTRX, "-f", NULL, ffREAD },
 153     { efNDX, NULL, NULL, ffOPTRD },
 154     { efDAT, "-oaxis",    "helixaxis", ffWRITE },
 155     { efDAT, "-ocenter",  "center", ffWRITE },
 156     { efXVG, "-orise",    "rise",ffWRITE },
 157     { efXVG, "-oradius",  "radius",ffWRITE },
 158     { efXVG, "-otwist",   "twist",ffWRITE },
 159     { efXVG, "-obending", "bending",ffWRITE },
 160     { efXVG, "-otilt",    "tilt", ffWRITE },
 161     { efXVG, "-orot",     "rotation",ffWRITE }
 162   };
 163 #define NFILE asize(fnm)
 164
 165
 166   CopyRight(stderr,argv[0]);
 167
 168   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,
 169                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
 170
 171   top=read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);
 172
 173     for(i=0;i<3;i++)
 174         weight[i]=1.0;
 175
 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);
 179     snew(x_CA,iCA);
 180     snew(x_SC,iCA); /* sic! */
 181
 182     snew(r12,iCA-3);
 183     snew(r23,iCA-3);
 184     snew(r34,iCA-3);
 185     snew(diff13,iCA-3);
 186     snew(diff24,iCA-3);
 187     snew(helixaxis,iCA-3);
 188     snew(twist,iCA);
 189     snew(residuetwist,iCA);
 190     snew(radius,iCA);
 191     snew(residueradius,iCA);
 192     snew(rise,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);
 201     snew(axis3_t0,iCA);
 202     snew(residuehelixaxis_tlast,iCA);
 203     snew(residuevector_tlast,iCA);
 204     snew(axis3_tlast,iCA);
 205     snew(axis3,iCA);
 206
 207     if(bSC)
 208     {
 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);
 211         if(iSC!=iCA)
 212             gmx_fatal(FARGS,"Number of sidechain atoms (%d) != number of CA atoms (%d)",iSC,iCA);
 213
 214     }
 215
 216     natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 217
 218     fpaxis=ffopen(opt2fn("-oaxis",NFILE,fnm),"w");
 219     fpcenter=ffopen(opt2fn("-ocenter",NFILE,fnm),"w");
 220     fprise=ffopen(opt2fn("-orise",NFILE,fnm),"w");
 221     fpradius=ffopen(opt2fn("-oradius",NFILE,fnm),"w");
 222     fptwist=ffopen(opt2fn("-otwist",NFILE,fnm),"w");
 223     fpbending=ffopen(opt2fn("-obending",NFILE,fnm),"w");
 224
 225     fptheta1=ffopen("theta1.xvg","w");
 226     fptheta2=ffopen("theta2.xvg","w");
 227     fptheta3=ffopen("theta3.xvg","w");
 228
 229     if(bIncremental)
 230     {
 231         fptilt=xvgropen(opt2fn("-otilt",NFILE,fnm),
 232                 "Incremental local helix tilt","Time(ps)","Tilt (degrees)",
 233                 oenv);
 234         fprotation=xvgropen(opt2fn("-orot",NFILE,fnm),
 235                 "Incremental local helix rotation","Time(ps)",
 236                 "Rotation (degrees)",oenv);
 237     }
 238     else
 239     {
 240         fptilt=xvgropen(opt2fn("-otilt",NFILE,fnm),
 241                 "Cumulative local helix tilt","Time(ps)","Tilt (degrees)",oenv);
 242         fprotation=xvgropen(opt2fn("-orot",NFILE,fnm),
 243                 "Cumulative local helix rotation","Time(ps)",
 244                 "Rotation (degrees)",oenv);
 245     }
 246
 247     clear_rvecs(3,unitaxes);
 248     unitaxes[0][0]=1;
 249     unitaxes[1][1]=1;
 250     unitaxes[2][2]=1;
 251
 252   do
 253   {
 254       /* initialisation for correct distance calculations */
 255     set_pbc(&pbc,ePBC,box);
 256       /* make molecules whole again */
 257     rm_pbc(&top->idef,ePBC,natoms,box,x,x);
 258
 259       /* copy coords to our smaller arrays */
 260       for(i=0;i<iCA;i++)
 261       {
 262           copy_rvec(x[ind_CA[i]],x_CA[i]);
 263           if(bSC)
 264           {
 265               copy_rvec(x[ind_SC[i]],x_SC[i]);
 266           }
 267       }
 268
 269       for(i=0;i<iCA-3;i++)
 270       {
 271           rvec_sub(x_CA[i+1],x_CA[i],r12[i]);
 272           rvec_sub(x_CA[i+2],x_CA[i+1],r23[i]);
 273           rvec_sub(x_CA[i+3],x_CA[i+2],r34[i]);
 274           rvec_sub(r12[i],r23[i],diff13[i]);
 275           rvec_sub(r23[i],r34[i],diff24[i]);
 276           /* calculate helix axis */
 277           cprod(diff13[i],diff24[i],helixaxis[i]);
 278           svmul(1.0/norm(helixaxis[i]),helixaxis[i],helixaxis[i]);
 279
 280           tmp = cos_angle(diff13[i],diff24[i]);
 281           twist[i] = 180.0/M_PI * acos( tmp );
 282           radius[i] = sqrt( norm(diff13[i])*norm(diff24[i]) ) / (2.0* (1.0-tmp) );
 283           rise[i]=fabs(iprod(r23[i],helixaxis[i]));
 284
 285           svmul(radius[i]/norm(diff13[i]),diff13[i],v1);
 286           svmul(radius[i]/norm(diff24[i]),diff24[i],v2);
 287
 288           rvec_sub(x_CA[i+1],v1,residueorigin[i+1]);
 289           rvec_sub(x_CA[i+2],v2,residueorigin[i+2]);
 290       }
 291       residueradius[0]=residuetwist[0]=residuerise[0]=0;
 292
 293       residueradius[1]=radius[0];
 294       residuetwist[1]=twist[0];
 295       residuerise[1]=rise[0];
 296
 297       residuebending[0]=residuebending[1]=0;
 298       for(i=2;i<iCA-2;i++)
 299         {
 300           residueradius[i]=0.5*(radius[i-2]+radius[i-1]);
 301           residuetwist[i]=0.5*(twist[i-2]+twist[i-1]);
 302           residuerise[i]=0.5*(rise[i-2]+rise[i-1]);
 303           residuebending[i] = 180.0/M_PI*acos( cos_angle(helixaxis[i-2],helixaxis[i-1]) );
 304         }
 305       residueradius[iCA-2]=radius[iCA-4];
 306       residuetwist[iCA-2]=twist[iCA-4];
 307       residuerise[iCA-2]=rise[iCA-4];
 308       residueradius[iCA-1]=residuetwist[iCA-1]=residuerise[iCA-1]=0;
 309       residuebending[iCA-2]=residuebending[iCA-1]=0;
 310
 311       clear_rvec(residueorigin[0]);
 312       clear_rvec(residueorigin[iCA-1]);
 313
 314       /* average helix axes to define them on the residues.
 315        * Just extrapolate second first/list atom.
 316        */
 317       copy_rvec(helixaxis[0],residuehelixaxis[0]);
 318       copy_rvec(helixaxis[0],residuehelixaxis[1]);
 319
 320       for(i=2;i<iCA-2;i++)
 321       {
 322           rvec_add(helixaxis[i-2],helixaxis[i-1],residuehelixaxis[i]);
 323           svmul(0.5,residuehelixaxis[i],residuehelixaxis[i]);
 324       }
 325       copy_rvec(helixaxis[iCA-4],residuehelixaxis[iCA-2]);
 326       copy_rvec(helixaxis[iCA-4],residuehelixaxis[iCA-1]);
 327
 328       /* Normalize the axis */
 329       for(i=0;i<iCA;i++)
 330       {
 331           svmul(1.0/norm(residuehelixaxis[i]),residuehelixaxis[i],residuehelixaxis[i]);
 332       }
 333
 334       /* calculate vector from origin to residue CA */
 335       fprintf(fpaxis,"%15.12g  ",t);
 336       fprintf(fpcenter,"%15.12g  ",t);
 337       fprintf(fprise,"%15.12g  ",t);
 338       fprintf(fpradius,"%15.12g  ",t);
 339       fprintf(fptwist,"%15.12g  ",t);
 340       fprintf(fpbending,"%15.12g  ",t);
 341
 342       for(i=0;i<iCA;i++)
 343       {
 344         if(i==0 || i==iCA-1)
 345           {
 346             fprintf(fpaxis,"%15.12g %15.12g %15.12g       ",0.0,0.0,0.0);
 347             fprintf(fpcenter,"%15.12g %15.12g %15.12g       ",0.0,0.0,0.0);
 348             fprintf(fprise,"%15.12g  ",0.0);
 349             fprintf(fpradius,"%15.12g  ",0.0);
 350             fprintf(fptwist,"%15.12g  ",0.0);
 351             fprintf(fpbending,"%15.12g  ",0.0);
 352           }
 353         else
 354           {
 355           rvec_sub( bSC ? x_SC[i] : x_CA[i] ,residueorigin[i], residuevector[i]);
 356           svmul(1.0/norm(residuevector[i]),residuevector[i],residuevector[i]);
 357           cprod(residuehelixaxis[i],residuevector[i],axis3[i]);
 358           fprintf(fpaxis,"%15.12g %15.12g %15.12g       ",residuehelixaxis[i][0],residuehelixaxis[i][1],residuehelixaxis[i][2]);
 359           fprintf(fpcenter,"%15.12g %15.12g %15.12g       ",residueorigin[i][0],residueorigin[i][1],residueorigin[i][2]);
 360
 361           fprintf(fprise,"%15.12g  ",residuerise[i]);
 362           fprintf(fpradius,"%15.12g  ",residueradius[i]);
 363           fprintf(fptwist,"%15.12g  ",residuetwist[i]);
 364           fprintf(fpbending,"%15.12g  ",residuebending[i]);
 365
 366           /* angle with local vector? */
 367           /*
 368            printf("res[%2d]:  axis: %g %g %g    origin: %g %g %g   vector: %g %g %g   angle: %g\n",i,
 369                  residuehelixaxis[i][0],
 370                  residuehelixaxis[i][1],
 371                  residuehelixaxis[i][2],
 372                  residueorigin[i][0],
 373                  residueorigin[i][1],
 374                  residueorigin[i][2],
 375                  residuevector[i][0],
 376                  residuevector[i][1],
 377                  residuevector[i][2],
 378                  180.0/M_PI*acos( cos_angle(residuevector[i],residuehelixaxis[i]) ));
 379            */
 380     /*      fprintf(fp,"%15.12g %15.12g %15.12g   %15.12g %15.12g %15.12g\n",
 381                   residuehelixaxis[i][0],
 382                   residuehelixaxis[i][1],
 383                   residuehelixaxis[i][2],
 384                   residuevector[i][0],
 385                   residuevector[i][1],
 386                   residuevector[i][2]);
 387      */
 388           }
 389       }
 390       fprintf(fprise,"\n");
 391       fprintf(fpradius,"\n");
 392       fprintf(fpaxis,"\n");
 393       fprintf(fpcenter,"\n");
 394       fprintf(fptwist,"\n");
 395       fprintf(fpbending,"\n");
 396
 397       if(teller==0)
 398       {
 399           for(i=0;i<iCA;i++)
 400           {
 401               copy_rvec(residuehelixaxis[i],residuehelixaxis_t0[i]);
 402               copy_rvec(residuevector[i],residuevector_t0[i]);
 403               copy_rvec(axis3[i],axis3_t0[i]);
 404           }
 405       }
 406       else
 407       {
 408           fprintf(fptilt,"%15.12g       ",t);
 409           fprintf(fprotation,"%15.12g       ",t);
 410           fprintf(fptheta1,"%15.12g      ",t);
 411           fprintf(fptheta2,"%15.12g      ",t);
 412           fprintf(fptheta3,"%15.12g      ",t);
 413
 414           for(i=0;i<iCA;i++)
 415           {
 416             if(i==0 || i==iCA-1)
 417             {
 418                 tilt=rotation=0;
 419             }
 420             else
 421             {
 422               if(!bIncremental)
 423               {
 424                   /* Total rotation & tilt */
 425                   copy_rvec(residuehelixaxis_t0[i],refaxes[0]);
 426                   copy_rvec(residuevector_t0[i],refaxes[1]);
 427                   copy_rvec(axis3_t0[i],refaxes[2]);
 428                }
 429               else
 430               {
 431                   /* Rotation/tilt since last step */
 432                   copy_rvec(residuehelixaxis_tlast[i],refaxes[0]);
 433                   copy_rvec(residuevector_tlast[i],refaxes[1]);
 434                   copy_rvec(axis3_tlast[i],refaxes[2]);
 435               }
 436               copy_rvec(residuehelixaxis[i],newaxes[0]);
 437               copy_rvec(residuevector[i],newaxes[1]);
 438               copy_rvec(axis3[i],newaxes[2]);
 439
 440               /*
 441               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",
 442                      teller,i,
 443                      refaxes[0][0],refaxes[0][1],refaxes[0][2],
 444                      refaxes[1][0],refaxes[1][1],refaxes[1][2],
 445                      refaxes[2][0],refaxes[2][1],refaxes[2][2],
 446                      newaxes[0][0],newaxes[0][1],newaxes[0][2],
 447                      newaxes[1][0],newaxes[1][1],newaxes[1][2],
 448                      newaxes[2][0],newaxes[2][1],newaxes[2][2]);
 449               */
 450
 451               /* rotate reference frame onto unit axes */
 452               calc_fit_R(3,3,weight,unitaxes,refaxes,A);
 453               for(j=0;j<3;j++)
 454                 {
 455                   mvmul(A,refaxes[j],rot_refaxes[j]);
 456                   mvmul(A,newaxes[j],rot_newaxes[j]);
 457                 }
 458
 459               /* Determine local rotation matrix A */
 460               calc_fit_R(3,3,weight,rot_newaxes,rot_refaxes,A);
 461               /* Calculate euler angles, from rotation order y-z-x, where
 462                * x is helixaxis, y residuevector, and z axis3.
 463                *
 464                * A contains rotation column vectors.
 465                */
 466
 467               /*
 468               printf("frame %d, i=%d, A: %g %g %g  , %g %g %g , %g %g %g\n",
 469                      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]);
 470               */
 471
 472               theta1 = 180.0/M_PI*atan2(A[0][2],A[0][0]);
 473               theta2 = 180.0/M_PI*asin(-A[0][1]);
 474               theta3 = 180.0/M_PI*atan2(A[2][1],A[1][1]);
 475
 476               tilt = sqrt(theta1*theta1+theta2*theta2);
 477               rotation = theta3;
 478               fprintf(fptheta1,"%15.12g  ",theta1);
 479               fprintf(fptheta2,"%15.12g  ",theta2);
 480               fprintf(fptheta3,"%15.12g  ",theta3);
 481
 482             }
 483             fprintf(fptilt,"%15.12g  ",tilt);
 484             fprintf(fprotation,"%15.12g  ",rotation);
 485           }
 486           fprintf(fptilt,"\n");
 487           fprintf(fprotation,"\n");
 488           fprintf(fptheta1,"\n");
 489           fprintf(fptheta2,"\n");
 490           fprintf(fptheta3,"\n");
 491       }
 492
 493       for(i=0;i<iCA;i++)
 494       {
 495           copy_rvec(residuehelixaxis[i],residuehelixaxis_tlast[i]);
 496           copy_rvec(residuevector[i],residuevector_tlast[i]);
 497           copy_rvec(axis3[i],axis3_tlast[i]);
 498       }
 499
 500       teller++;
 501   } while (read_next_x(oenv,status,&t,natoms,x,box));
 502
 503     ffclose(fpaxis);
 504     ffclose(fpcenter);
 505     ffclose(fptilt);
 506     ffclose(fprotation);
 507     ffclose(fprise);
 508     ffclose(fpradius);
 509     ffclose(fptwist);
 510     ffclose(fpbending);
 511     ffclose(fptheta1);
 512     ffclose(fptheta2);
 513     ffclose(fptheta3);
 514
 515     close_trj(status);
 516
 517     thanx(stderr);
 518     return 0;
 519 }