src/gromacs/gmxana/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 "statutil.h"
  45 #include "string2.h"
  46 #include "vec.h"
  47 #include "index.h"
  48 #include "pbc.h"
  49 #include "gmx_fatal.h"
  50 #include "futil.h"
  51 #include "gstat.h"
  52 #include "pbc.h"
  53 #include "do_fit.h"
  54 #include "gmx_ana.h"
  55
  56
  57 int gmx_helixorient(int argc, char *argv[])
  58 {
  59     const char      *desc[] = {
  60         "[TT]g_helixorient[tt] calculates the coordinates and direction of the average",
  61         "axis inside an alpha helix, and the direction/vectors of both the",
  62         "C[GRK]alpha[grk] and (optionally) a sidechain atom relative to the axis.[PAR]",
  63         "As input, you need to specify an index group with C[GRK]alpha[grk] atoms",
  64         "corresponding to an [GRK]alpha[grk]-helix of continuous residues. Sidechain",
  65         "directions require a second index group of the same size, containing",
  66         "the heavy atom in each residue that should represent the sidechain.[PAR]",
  67         "[BB]Note[bb] that this program does not do any fitting of structures.[PAR]",
  68         "We need four C[GRK]alpha[grk] coordinates to define the local direction of the helix",
  69         "axis.[PAR]",
  70         "The tilt/rotation is calculated from Euler rotations, where we define",
  71         "the helix axis as the local [IT]x[it]-axis, the residues/C[GRK]alpha[grk] vector as [IT]y[it], and the",
  72         "[IT]z[it]-axis from their cross product. We use the Euler Y-Z-X rotation, meaning",
  73         "we first tilt the helix axis (1) around and (2) orthogonal to the residues",
  74         "vector, and finally apply the (3) rotation around it. For debugging or other",
  75         "purposes, we also write out the actual Euler rotation angles as [TT]theta[1-3].xvg[tt]"
  76     };
  77
  78     t_topology      *top = NULL;
  79     real             t;
  80     rvec            *x = NULL, dx;
  81     matrix           box;
  82     t_trxstatus     *status;
  83     int              natoms;
  84     real             theta1, theta2, theta3;
  85
  86     int              d, i, j, teller = 0;
  87     int              iCA, iSC;
  88     atom_id         *ind_CA;
  89     atom_id         *ind_SC;
  90     char            *gn_CA;
  91     char            *gn_SC;
  92     rvec             averageaxis;
  93     rvec             v1, v2, p1, p2, vtmp, vproj;
  94     rvec            *x_CA, *x_SC;
  95     rvec            *r12;
  96     rvec            *r23;
  97     rvec            *r34;
  98     rvec            *diff13;
  99     rvec            *diff24;
 100     rvec            *helixaxis;
 101     rvec            *residuehelixaxis;
 102     rvec            *residueorigin;
 103     rvec            *residuevector;
 104     rvec            *sidechainvector;
 105
 106     rvec             axes_t0[3];
 107     rvec             axes[3];
 108     rvec            *residuehelixaxis_t0;
 109     rvec            *residuevector_t0;
 110     rvec            *axis3_t0;
 111     rvec            *residuehelixaxis_tlast;
 112     rvec            *residuevector_tlast;
 113     rvec            *axis3_tlast;
 114     rvec             refaxes[3], newaxes[3];
 115     rvec             unitaxes[3];
 116     rvec             rot_refaxes[3], rot_newaxes[3];
 117
 118     real             tilt, rotation;
 119     rvec            *axis3;
 120     real            *twist, *residuetwist;
 121     real            *radius, *residueradius;
 122     real            *rise, *residuerise;
 123     real            *residuebending;
 124
 125     real             tmp, rotangle;
 126     real             weight[3];
 127     t_pbc            pbc;
 128     matrix           A;
 129
 130     FILE            *fpaxis, *fpcenter, *fptilt, *fprotation;
 131     FILE            *fpradius, *fprise, *fptwist;
 132     FILE            *fptheta1, *fptheta2, *fptheta3;
 133     FILE            *fpbending;
 134     int              ePBC;
 135
 136     output_env_t     oenv;
 137     gmx_rmpbc_t      gpbc = NULL;
 138
 139     static  gmx_bool bSC          = FALSE;
 140     static gmx_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     parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_BE_NICE,
 166                       NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv);
 167
 168     top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC);
 169
 170     for (i = 0; i < 3; i++)
 171     {
 172         weight[i] = 1.0;
 173     }
 174
 175     /* read index files */
 176     printf("Select a group of Calpha atoms corresponding to a single continuous helix:\n");
 177     get_index(&(top->atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &iCA, &ind_CA, &gn_CA);
 178     snew(x_CA, iCA);
 179     snew(x_SC, iCA); /* sic! */
 180
 181     snew(r12, iCA-3);
 182     snew(r23, iCA-3);
 183     snew(r34, iCA-3);
 184     snew(diff13, iCA-3);
 185     snew(diff24, iCA-3);
 186     snew(helixaxis, iCA-3);
 187     snew(twist, iCA);
 188     snew(residuetwist, iCA);
 189     snew(radius, iCA);
 190     snew(residueradius, iCA);
 191     snew(rise, iCA);
 192     snew(residuerise, iCA);
 193     snew(residueorigin, iCA);
 194     snew(residuehelixaxis, iCA);
 195     snew(residuevector, iCA);
 196     snew(sidechainvector, iCA);
 197     snew(residuebending, iCA);
 198     snew(residuehelixaxis_t0, iCA);
 199     snew(residuevector_t0, iCA);
 200     snew(axis3_t0, iCA);
 201     snew(residuehelixaxis_tlast, iCA);
 202     snew(residuevector_tlast, iCA);
 203     snew(axis3_tlast, iCA);
 204     snew(axis3, iCA);
 205
 206     if (bSC)
 207     {
 208         printf("Select a group of atoms defining the sidechain direction (1/residue):\n");
 209         get_index(&(top->atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &iSC, &ind_SC, &gn_SC);
 210         if (iSC != iCA)
 211         {
 212             gmx_fatal(FARGS, "Number of sidechain atoms (%d) != number of CA atoms (%d)", iSC, iCA);
 213         }
 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);
 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     }
 505     while (read_next_x(oenv, status, &t, x, box));
 506
 507     gmx_rmpbc_done(gpbc);
 508
 509     ffclose(fpaxis);
 510     ffclose(fpcenter);
 511     ffclose(fptilt);
 512     ffclose(fprotation);
 513     ffclose(fprise);
 514     ffclose(fpradius);
 515     ffclose(fptwist);
 516     ffclose(fpbending);
 517     ffclose(fptheta1);
 518     ffclose(fptheta2);
 519     ffclose(fptheta3);
 520
 521     close_trj(status);
 522
 523     return 0;
 524 }