src/tools/gmx_helixorient.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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, 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.
  11  *
  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.
  16  *
  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.
  21  *
  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.
  26  *
  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.
  34  *
  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.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41 #include <typedefs.h>
  42
  43 #include "smalloc.h"
  44 #include "macros.h"
  45 #include <math.h>
  46 #include "xvgr.h"
  47 #include "copyrite.h"
  48 #include "statutil.h"
  49 #include "string2.h"
  50 #include "vec.h"
  51 #include "index.h"
  52 #include "pbc.h"
  53 #include "gmx_fatal.h"
  54 #include "futil.h"
  55 #include "gstat.h"
  56 #include "pbc.h"
  57 #include "do_fit.h"
  58 #include "gmx_ana.h"
  59
  60
  61 int gmx_helixorient(int argc,char *argv[])
  62 {
  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",
  73     "axis.[PAR]",
  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]"
  80     };
  81
  82     t_topology *top=NULL;
  83     real t;
  84     rvec *x=NULL,dx;
  85     matrix box;
  86     t_trxstatus *status;
  87     int natoms;
  88     real theta1,theta2,theta3;
  89
  90     int d,i,j,teller=0;
  91     int iCA,iSC;
  92     atom_id *ind_CA;
  93     atom_id *ind_SC;
  94     char *gn_CA;
  95     char *gn_SC;
  96     rvec averageaxis;
  97     rvec v1,v2,p1,p2,vtmp,vproj;
  98     rvec *x_CA,*x_SC;
  99     rvec *r12;
 100     rvec *r23;
 101     rvec *r34;
 102     rvec *diff13;
 103     rvec *diff24;
 104     rvec *helixaxis;
 105     rvec *residuehelixaxis;
 106     rvec *residueorigin;
 107     rvec *residuevector;
 108     rvec *sidechainvector;
 109
 110     rvec axes_t0[3];
 111     rvec axes[3];
 112     rvec *residuehelixaxis_t0;
 113     rvec *residuevector_t0;
 114     rvec *axis3_t0;
 115     rvec *residuehelixaxis_tlast;
 116     rvec *residuevector_tlast;
 117     rvec *axis3_tlast;
 118     rvec refaxes[3],newaxes[3];
 119     rvec unitaxes[3];
 120     rvec rot_refaxes[3],rot_newaxes[3];
 121
 122     real tilt,rotation;
 123     rvec *axis3;
 124     real *twist,*residuetwist;
 125     real *radius,*residueradius;
 126     real *rise,*residuerise;
 127     real *residuebending;
 128
 129     real tmp,rotangle;
 130     real weight[3];
 131     t_pbc   pbc;
 132     matrix A;
 133
 134     FILE *fpaxis,*fpcenter,*fptilt,*fprotation;
 135     FILE *fpradius,*fprise,*fptwist;
 136     FILE *fptheta1,*fptheta2,*fptheta3;
 137     FILE *fpbending;
 138     int ePBC;
 139
 140     output_env_t oenv;
 141     gmx_rmpbc_t  gpbc=NULL;
 142
 143     static  gmx_bool bSC=FALSE;
 144     static gmx_bool bIncremental = FALSE;
 145
 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." },
 151     };
 152 #define NPA asize(pa)
 153
 154   t_filenm fnm[] = {
 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 }
 166   };
 167 #define NFILE asize(fnm)
 168
 169
 170   CopyRight(stderr,argv[0]);
 171
 172   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,
 173                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
 174
 175   top=read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);
 176
 177     for(i=0;i<3;i++)
 178         weight[i]=1.0;
 179
 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);
 183     snew(x_CA,iCA);
 184     snew(x_SC,iCA); /* sic! */
 185
 186     snew(r12,iCA-3);
 187     snew(r23,iCA-3);
 188     snew(r34,iCA-3);
 189     snew(diff13,iCA-3);
 190     snew(diff24,iCA-3);
 191     snew(helixaxis,iCA-3);
 192     snew(twist,iCA);
 193     snew(residuetwist,iCA);
 194     snew(radius,iCA);
 195     snew(residueradius,iCA);
 196     snew(rise,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);
 205     snew(axis3_t0,iCA);
 206     snew(residuehelixaxis_tlast,iCA);
 207     snew(residuevector_tlast,iCA);
 208     snew(axis3_tlast,iCA);
 209     snew(axis3,iCA);
 210
 211     if(bSC)
 212     {
 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);
 215         if(iSC!=iCA)
 216             gmx_fatal(FARGS,"Number of sidechain atoms (%d) != number of CA atoms (%d)",iSC,iCA);
 217
 218     }
 219
 220     natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 221
 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");
 228
 229     fptheta1=ffopen("theta1.xvg","w");
 230     fptheta2=ffopen("theta2.xvg","w");
 231     fptheta3=ffopen("theta3.xvg","w");
 232
 233     if(bIncremental)
 234     {
 235         fptilt=xvgropen(opt2fn("-otilt",NFILE,fnm),
 236                 "Incremental local helix tilt","Time(ps)","Tilt (degrees)",
 237                 oenv);
 238         fprotation=xvgropen(opt2fn("-orot",NFILE,fnm),
 239                 "Incremental local helix rotation","Time(ps)",
 240                 "Rotation (degrees)",oenv);
 241     }
 242     else
 243     {
 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);
 249     }
 250
 251     clear_rvecs(3,unitaxes);
 252     unitaxes[0][0]=1;
 253     unitaxes[1][1]=1;
 254     unitaxes[2][2]=1;
 255
 256     gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
 257
 258   do
 259   {
 260       /* initialisation for correct distance calculations */
 261     set_pbc(&pbc,ePBC,box);
 262       /* make molecules whole again */
 263     gmx_rmpbc(gpbc,natoms,box,x);
 264
 265       /* copy coords to our smaller arrays */
 266       for(i=0;i<iCA;i++)
 267       {
 268           copy_rvec(x[ind_CA[i]],x_CA[i]);
 269           if(bSC)
 270           {
 271               copy_rvec(x[ind_SC[i]],x_SC[i]);
 272           }
 273       }
 274
 275       for(i=0;i<iCA-3;i++)
 276       {
 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]);
 285
 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]));
 290
 291           svmul(radius[i]/norm(diff13[i]),diff13[i],v1);
 292           svmul(radius[i]/norm(diff24[i]),diff24[i],v2);
 293
 294           rvec_sub(x_CA[i+1],v1,residueorigin[i+1]);
 295           rvec_sub(x_CA[i+2],v2,residueorigin[i+2]);
 296       }
 297       residueradius[0]=residuetwist[0]=residuerise[0]=0;
 298
 299       residueradius[1]=radius[0];
 300       residuetwist[1]=twist[0];
 301       residuerise[1]=rise[0];
 302
 303       residuebending[0]=residuebending[1]=0;
 304       for(i=2;i<iCA-2;i++)
 305         {
 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]) );
 310         }
 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;
 316
 317       clear_rvec(residueorigin[0]);
 318       clear_rvec(residueorigin[iCA-1]);
 319
 320       /* average helix axes to define them on the residues.
 321        * Just extrapolate second first/list atom.
 322        */
 323       copy_rvec(helixaxis[0],residuehelixaxis[0]);
 324       copy_rvec(helixaxis[0],residuehelixaxis[1]);
 325
 326       for(i=2;i<iCA-2;i++)
 327       {
 328           rvec_add(helixaxis[i-2],helixaxis[i-1],residuehelixaxis[i]);
 329           svmul(0.5,residuehelixaxis[i],residuehelixaxis[i]);
 330       }
 331       copy_rvec(helixaxis[iCA-4],residuehelixaxis[iCA-2]);
 332       copy_rvec(helixaxis[iCA-4],residuehelixaxis[iCA-1]);
 333
 334       /* Normalize the axis */
 335       for(i=0;i<iCA;i++)
 336       {
 337           svmul(1.0/norm(residuehelixaxis[i]),residuehelixaxis[i],residuehelixaxis[i]);
 338       }
 339
 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);
 347
 348       for(i=0;i<iCA;i++)
 349       {
 350         if(i==0 || i==iCA-1)
 351           {
 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);
 358           }
 359         else
 360           {
 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]);
 366
 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]);
 371
 372           /* angle with local vector? */
 373           /*
 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],
 378                  residueorigin[i][0],
 379                  residueorigin[i][1],
 380                  residueorigin[i][2],
 381                  residuevector[i][0],
 382                  residuevector[i][1],
 383                  residuevector[i][2],
 384                  180.0/M_PI*acos( cos_angle(residuevector[i],residuehelixaxis[i]) ));
 385            */
 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],
 390                   residuevector[i][0],
 391                   residuevector[i][1],
 392                   residuevector[i][2]);
 393      */
 394           }
 395       }
 396       fprintf(fprise,"\n");
 397       fprintf(fpradius,"\n");
 398       fprintf(fpaxis,"\n");
 399       fprintf(fpcenter,"\n");
 400       fprintf(fptwist,"\n");
 401       fprintf(fpbending,"\n");
 402
 403       if(teller==0)
 404       {
 405           for(i=0;i<iCA;i++)
 406           {
 407               copy_rvec(residuehelixaxis[i],residuehelixaxis_t0[i]);
 408               copy_rvec(residuevector[i],residuevector_t0[i]);
 409               copy_rvec(axis3[i],axis3_t0[i]);
 410           }
 411       }
 412       else
 413       {
 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);
 419
 420           for(i=0;i<iCA;i++)
 421           {
 422             if(i==0 || i==iCA-1)
 423             {
 424                 tilt=rotation=0;
 425             }
 426             else
 427             {
 428               if(!bIncremental)
 429               {
 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]);
 434                }
 435               else
 436               {
 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]);
 441               }
 442               copy_rvec(residuehelixaxis[i],newaxes[0]);
 443               copy_rvec(residuevector[i],newaxes[1]);
 444               copy_rvec(axis3[i],newaxes[2]);
 445
 446               /*
 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",
 448                      teller,i,
 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]);
 455               */
 456
 457               /* rotate reference frame onto unit axes */
 458               calc_fit_R(3,3,weight,unitaxes,refaxes,A);
 459               for(j=0;j<3;j++)
 460                 {
 461                   mvmul(A,refaxes[j],rot_refaxes[j]);
 462                   mvmul(A,newaxes[j],rot_newaxes[j]);
 463                 }
 464
 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.
 469                *
 470                * A contains rotation column vectors.
 471                */
 472
 473               /*
 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]);
 476               */
 477
 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]);
 481
 482               tilt = sqrt(theta1*theta1+theta2*theta2);
 483               rotation = theta3;
 484               fprintf(fptheta1,"%15.12g  ",theta1);
 485               fprintf(fptheta2,"%15.12g  ",theta2);
 486               fprintf(fptheta3,"%15.12g  ",theta3);
 487
 488             }
 489             fprintf(fptilt,"%15.12g  ",tilt);
 490             fprintf(fprotation,"%15.12g  ",rotation);
 491           }
 492           fprintf(fptilt,"\n");
 493           fprintf(fprotation,"\n");
 494           fprintf(fptheta1,"\n");
 495           fprintf(fptheta2,"\n");
 496           fprintf(fptheta3,"\n");
 497       }
 498
 499       for(i=0;i<iCA;i++)
 500       {
 501           copy_rvec(residuehelixaxis[i],residuehelixaxis_tlast[i]);
 502           copy_rvec(residuevector[i],residuevector_tlast[i]);
 503           copy_rvec(axis3[i],axis3_tlast[i]);
 504       }
 505
 506       teller++;
 507   } while (read_next_x(oenv,status,&t,natoms,x,box));
 508
 509     gmx_rmpbc_done(gpbc);
 510
 511     ffclose(fpaxis);
 512     ffclose(fpcenter);
 513     ffclose(fptilt);
 514     ffclose(fprotation);
 515     ffclose(fprise);
 516     ffclose(fpradius);
 517     ffclose(fptwist);
 518     ffclose(fpbending);
 519     ffclose(fptheta1);
 520     ffclose(fptheta2);
 521     ffclose(fptheta3);
 522
 523     close_trj(status);
 524
 525     thanx(stderr);
 526     return 0;
 527 }