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 "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     parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_BE_NICE,
 167                       NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv);
 168
 169     top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC);
 170
 171     for (i = 0; i < 3; i++)
 172     {
 173         weight[i] = 1.0;
 174     }
 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         {
 213             gmx_fatal(FARGS, "Number of sidechain atoms (%d) != number of CA atoms (%d)", iSC, iCA);
 214         }
 215
 216     }
 217
 218     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 219
 220     fpaxis    = ffopen(opt2fn("-oaxis", NFILE, fnm), "w");
 221     fpcenter  = ffopen(opt2fn("-ocenter", NFILE, fnm), "w");
 222     fprise    = ffopen(opt2fn("-orise", NFILE, fnm), "w");
 223     fpradius  = ffopen(opt2fn("-oradius", NFILE, fnm), "w");
 224     fptwist   = ffopen(opt2fn("-otwist", NFILE, fnm), "w");
 225     fpbending = ffopen(opt2fn("-obending", NFILE, fnm), "w");
 226
 227     fptheta1 = ffopen("theta1.xvg", "w");
 228     fptheta2 = ffopen("theta2.xvg", "w");
 229     fptheta3 = ffopen("theta3.xvg", "w");
 230
 231     if (bIncremental)
 232     {
 233         fptilt = xvgropen(opt2fn("-otilt", NFILE, fnm),
 234                           "Incremental local helix tilt", "Time(ps)", "Tilt (degrees)",
 235                           oenv);
 236         fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
 237                               "Incremental local helix rotation", "Time(ps)",
 238                               "Rotation (degrees)", oenv);
 239     }
 240     else
 241     {
 242         fptilt = xvgropen(opt2fn("-otilt", NFILE, fnm),
 243                           "Cumulative local helix tilt", "Time(ps)", "Tilt (degrees)", oenv);
 244         fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
 245                               "Cumulative local helix rotation", "Time(ps)",
 246                               "Rotation (degrees)", oenv);
 247     }
 248
 249     clear_rvecs(3, unitaxes);
 250     unitaxes[0][0] = 1;
 251     unitaxes[1][1] = 1;
 252     unitaxes[2][2] = 1;
 253
 254     gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms, box);
 255
 256     do
 257     {
 258         /* initialisation for correct distance calculations */
 259         set_pbc(&pbc, ePBC, box);
 260         /* make molecules whole again */
 261         gmx_rmpbc(gpbc, natoms, box, x);
 262
 263         /* copy coords to our smaller arrays */
 264         for (i = 0; i < iCA; i++)
 265         {
 266             copy_rvec(x[ind_CA[i]], x_CA[i]);
 267             if (bSC)
 268             {
 269                 copy_rvec(x[ind_SC[i]], x_SC[i]);
 270             }
 271         }
 272
 273         for (i = 0; i < iCA-3; i++)
 274         {
 275             rvec_sub(x_CA[i+1], x_CA[i], r12[i]);
 276             rvec_sub(x_CA[i+2], x_CA[i+1], r23[i]);
 277             rvec_sub(x_CA[i+3], x_CA[i+2], r34[i]);
 278             rvec_sub(r12[i], r23[i], diff13[i]);
 279             rvec_sub(r23[i], r34[i], diff24[i]);
 280             /* calculate helix axis */
 281             cprod(diff13[i], diff24[i], helixaxis[i]);
 282             svmul(1.0/norm(helixaxis[i]), helixaxis[i], helixaxis[i]);
 283
 284             tmp       = cos_angle(diff13[i], diff24[i]);
 285             twist[i]  = 180.0/M_PI * acos( tmp );
 286             radius[i] = sqrt( norm(diff13[i])*norm(diff24[i]) ) / (2.0* (1.0-tmp) );
 287             rise[i]   = fabs(iprod(r23[i], helixaxis[i]));
 288
 289             svmul(radius[i]/norm(diff13[i]), diff13[i], v1);
 290             svmul(radius[i]/norm(diff24[i]), diff24[i], v2);
 291
 292             rvec_sub(x_CA[i+1], v1, residueorigin[i+1]);
 293             rvec_sub(x_CA[i+2], v2, residueorigin[i+2]);
 294         }
 295         residueradius[0] = residuetwist[0] = residuerise[0] = 0;
 296
 297         residueradius[1] = radius[0];
 298         residuetwist[1]  = twist[0];
 299         residuerise[1]   = rise[0];
 300
 301         residuebending[0] = residuebending[1] = 0;
 302         for (i = 2; i < iCA-2; i++)
 303         {
 304             residueradius[i]  = 0.5*(radius[i-2]+radius[i-1]);
 305             residuetwist[i]   = 0.5*(twist[i-2]+twist[i-1]);
 306             residuerise[i]    = 0.5*(rise[i-2]+rise[i-1]);
 307             residuebending[i] = 180.0/M_PI*acos( cos_angle(helixaxis[i-2], helixaxis[i-1]) );
 308         }
 309         residueradius[iCA-2]  = radius[iCA-4];
 310         residuetwist[iCA-2]   = twist[iCA-4];
 311         residuerise[iCA-2]    = rise[iCA-4];
 312         residueradius[iCA-1]  = residuetwist[iCA-1] = residuerise[iCA-1] = 0;
 313         residuebending[iCA-2] = residuebending[iCA-1] = 0;
 314
 315         clear_rvec(residueorigin[0]);
 316         clear_rvec(residueorigin[iCA-1]);
 317
 318         /* average helix axes to define them on the residues.
 319          * Just extrapolate second first/list atom.
 320          */
 321         copy_rvec(helixaxis[0], residuehelixaxis[0]);
 322         copy_rvec(helixaxis[0], residuehelixaxis[1]);
 323
 324         for (i = 2; i < iCA-2; i++)
 325         {
 326             rvec_add(helixaxis[i-2], helixaxis[i-1], residuehelixaxis[i]);
 327             svmul(0.5, residuehelixaxis[i], residuehelixaxis[i]);
 328         }
 329         copy_rvec(helixaxis[iCA-4], residuehelixaxis[iCA-2]);
 330         copy_rvec(helixaxis[iCA-4], residuehelixaxis[iCA-1]);
 331
 332         /* Normalize the axis */
 333         for (i = 0; i < iCA; i++)
 334         {
 335             svmul(1.0/norm(residuehelixaxis[i]), residuehelixaxis[i], residuehelixaxis[i]);
 336         }
 337
 338         /* calculate vector from origin to residue CA */
 339         fprintf(fpaxis, "%15.12g  ", t);
 340         fprintf(fpcenter, "%15.12g  ", t);
 341         fprintf(fprise, "%15.12g  ", t);
 342         fprintf(fpradius, "%15.12g  ", t);
 343         fprintf(fptwist, "%15.12g  ", t);
 344         fprintf(fpbending, "%15.12g  ", t);
 345
 346         for (i = 0; i < iCA; i++)
 347         {
 348             if (i == 0 || i == iCA-1)
 349             {
 350                 fprintf(fpaxis, "%15.12g %15.12g %15.12g       ", 0.0, 0.0, 0.0);
 351                 fprintf(fpcenter, "%15.12g %15.12g %15.12g       ", 0.0, 0.0, 0.0);
 352                 fprintf(fprise, "%15.12g  ", 0.0);
 353                 fprintf(fpradius, "%15.12g  ", 0.0);
 354                 fprintf(fptwist, "%15.12g  ", 0.0);
 355                 fprintf(fpbending, "%15.12g  ", 0.0);
 356             }
 357             else
 358             {
 359                 rvec_sub( bSC ? x_SC[i] : x_CA[i], residueorigin[i], residuevector[i]);
 360                 svmul(1.0/norm(residuevector[i]), residuevector[i], residuevector[i]);
 361                 cprod(residuehelixaxis[i], residuevector[i], axis3[i]);
 362                 fprintf(fpaxis, "%15.12g %15.12g %15.12g       ", residuehelixaxis[i][0], residuehelixaxis[i][1], residuehelixaxis[i][2]);
 363                 fprintf(fpcenter, "%15.12g %15.12g %15.12g       ", residueorigin[i][0], residueorigin[i][1], residueorigin[i][2]);
 364
 365                 fprintf(fprise, "%15.12g  ", residuerise[i]);
 366                 fprintf(fpradius, "%15.12g  ", residueradius[i]);
 367                 fprintf(fptwist, "%15.12g  ", residuetwist[i]);
 368                 fprintf(fpbending, "%15.12g  ", residuebending[i]);
 369
 370                 /* angle with local vector? */
 371                 /*
 372                    printf("res[%2d]:  axis: %g %g %g    origin: %g %g %g   vector: %g %g %g   angle: %g\n",i,
 373                        residuehelixaxis[i][0],
 374                        residuehelixaxis[i][1],
 375                        residuehelixaxis[i][2],
 376                        residueorigin[i][0],
 377                        residueorigin[i][1],
 378                        residueorigin[i][2],
 379                        residuevector[i][0],
 380                        residuevector[i][1],
 381                        residuevector[i][2],
 382                        180.0/M_PI*acos( cos_angle(residuevector[i],residuehelixaxis[i]) ));
 383                  */
 384                 /*      fprintf(fp,"%15.12g %15.12g %15.12g   %15.12g %15.12g %15.12g\n",
 385                               residuehelixaxis[i][0],
 386                               residuehelixaxis[i][1],
 387                               residuehelixaxis[i][2],
 388                               residuevector[i][0],
 389                               residuevector[i][1],
 390                               residuevector[i][2]);
 391                  */
 392             }
 393         }
 394         fprintf(fprise, "\n");
 395         fprintf(fpradius, "\n");
 396         fprintf(fpaxis, "\n");
 397         fprintf(fpcenter, "\n");
 398         fprintf(fptwist, "\n");
 399         fprintf(fpbending, "\n");
 400
 401         if (teller == 0)
 402         {
 403             for (i = 0; i < iCA; i++)
 404             {
 405                 copy_rvec(residuehelixaxis[i], residuehelixaxis_t0[i]);
 406                 copy_rvec(residuevector[i], residuevector_t0[i]);
 407                 copy_rvec(axis3[i], axis3_t0[i]);
 408             }
 409         }
 410         else
 411         {
 412             fprintf(fptilt, "%15.12g       ", t);
 413             fprintf(fprotation, "%15.12g       ", t);
 414             fprintf(fptheta1, "%15.12g      ", t);
 415             fprintf(fptheta2, "%15.12g      ", t);
 416             fprintf(fptheta3, "%15.12g      ", t);
 417
 418             for (i = 0; i < iCA; i++)
 419             {
 420                 if (i == 0 || i == iCA-1)
 421                 {
 422                     tilt = rotation = 0;
 423                 }
 424                 else
 425                 {
 426                     if (!bIncremental)
 427                     {
 428                         /* Total rotation & tilt */
 429                         copy_rvec(residuehelixaxis_t0[i], refaxes[0]);
 430                         copy_rvec(residuevector_t0[i], refaxes[1]);
 431                         copy_rvec(axis3_t0[i], refaxes[2]);
 432                     }
 433                     else
 434                     {
 435                         /* Rotation/tilt since last step */
 436                         copy_rvec(residuehelixaxis_tlast[i], refaxes[0]);
 437                         copy_rvec(residuevector_tlast[i], refaxes[1]);
 438                         copy_rvec(axis3_tlast[i], refaxes[2]);
 439                     }
 440                     copy_rvec(residuehelixaxis[i], newaxes[0]);
 441                     copy_rvec(residuevector[i], newaxes[1]);
 442                     copy_rvec(axis3[i], newaxes[2]);
 443
 444                     /*
 445                        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                        teller,i,
 447                        refaxes[0][0],refaxes[0][1],refaxes[0][2],
 448                                refaxes[1][0],refaxes[1][1],refaxes[1][2],
 449                                refaxes[2][0],refaxes[2][1],refaxes[2][2],
 450                                newaxes[0][0],newaxes[0][1],newaxes[0][2],
 451                                newaxes[1][0],newaxes[1][1],newaxes[1][2],
 452                                newaxes[2][0],newaxes[2][1],newaxes[2][2]);
 453                      */
 454
 455                     /* rotate reference frame onto unit axes */
 456                     calc_fit_R(3, 3, weight, unitaxes, refaxes, A);
 457                     for (j = 0; j < 3; j++)
 458                     {
 459                         mvmul(A, refaxes[j], rot_refaxes[j]);
 460                         mvmul(A, newaxes[j], rot_newaxes[j]);
 461                     }
 462
 463                     /* Determine local rotation matrix A */
 464                     calc_fit_R(3, 3, weight, rot_newaxes, rot_refaxes, A);
 465                     /* Calculate euler angles, from rotation order y-z-x, where
 466                      * x is helixaxis, y residuevector, and z axis3.
 467                      *
 468                      * A contains rotation column vectors.
 469                      */
 470
 471                     /*
 472                        printf("frame %d, i=%d, A: %g %g %g  , %g %g %g , %g %g %g\n",
 473                        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]);
 474                      */
 475
 476                     theta1 = 180.0/M_PI*atan2(A[0][2], A[0][0]);
 477                     theta2 = 180.0/M_PI*asin(-A[0][1]);
 478                     theta3 = 180.0/M_PI*atan2(A[2][1], A[1][1]);
 479
 480                     tilt     = sqrt(theta1*theta1+theta2*theta2);
 481                     rotation = theta3;
 482                     fprintf(fptheta1, "%15.12g  ", theta1);
 483                     fprintf(fptheta2, "%15.12g  ", theta2);
 484                     fprintf(fptheta3, "%15.12g  ", theta3);
 485
 486                 }
 487                 fprintf(fptilt, "%15.12g  ", tilt);
 488                 fprintf(fprotation, "%15.12g  ", rotation);
 489             }
 490             fprintf(fptilt, "\n");
 491             fprintf(fprotation, "\n");
 492             fprintf(fptheta1, "\n");
 493             fprintf(fptheta2, "\n");
 494             fprintf(fptheta3, "\n");
 495         }
 496
 497         for (i = 0; i < iCA; i++)
 498         {
 499             copy_rvec(residuehelixaxis[i], residuehelixaxis_tlast[i]);
 500             copy_rvec(residuevector[i], residuevector_tlast[i]);
 501             copy_rvec(axis3[i], axis3_tlast[i]);
 502         }
 503
 504         teller++;
 505     }
 506     while (read_next_x(oenv, status, &t, natoms, x, box));
 507
 508     gmx_rmpbc_done(gpbc);
 509
 510     ffclose(fpaxis);
 511     ffclose(fpcenter);
 512     ffclose(fptilt);
 513     ffclose(fprotation);
 514     ffclose(fprise);
 515     ffclose(fpradius);
 516     ffclose(fptwist);
 517     ffclose(fpbending);
 518     ffclose(fptheta1);
 519     ffclose(fptheta2);
 520     ffclose(fptheta3);
 521
 522     close_trj(status);
 523
 524     thanx(stderr);
 525     return 0;
 526 }