src/gromacs/gmxana/gmx_helixorient.cpp

   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  * Copyright (c) 2013,2014,2015,2017, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include <cmath>
  40
  41 #include "gromacs/commandline/pargs.h"
  42 #include "gromacs/fileio/trxio.h"
  43 #include "gromacs/fileio/xvgr.h"
  44 #include "gromacs/gmxana/gmx_ana.h"
  45 #include "gromacs/gmxana/gstat.h"
  46 #include "gromacs/math/do_fit.h"
  47 #include "gromacs/math/utilities.h"
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/pbcutil/pbc.h"
  50 #include "gromacs/pbcutil/rmpbc.h"
  51 #include "gromacs/topology/index.h"
  52 #include "gromacs/topology/topology.h"
  53 #include "gromacs/utility/arraysize.h"
  54 #include "gromacs/utility/fatalerror.h"
  55 #include "gromacs/utility/futil.h"
  56 #include "gromacs/utility/smalloc.h"
  57
  58 int gmx_helixorient(int argc, char *argv[])
  59 {
  60     const char       *desc[] = {
  61         "[THISMODULE] 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 = nullptr;
  80     real              t;
  81     rvec             *x = nullptr;
  82     matrix            box;
  83     t_trxstatus      *status;
  84     int               natoms;
  85     real              theta1, theta2, theta3;
  86
  87     int               i, j, teller = 0;
  88     int               iCA, iSC;
  89     int              *ind_CA;
  90     int              *ind_SC;
  91     char             *gn_CA;
  92     char             *gn_SC;
  93     rvec              v1, v2;
  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             *residuehelixaxis_t0;
 107     rvec             *residuevector_t0;
 108     rvec             *axis3_t0;
 109     rvec             *residuehelixaxis_tlast;
 110     rvec             *residuevector_tlast;
 111     rvec             *axis3_tlast;
 112     rvec              refaxes[3], newaxes[3];
 113     rvec              unitaxes[3];
 114     rvec              rot_refaxes[3], rot_newaxes[3];
 115
 116     real              tilt, rotation;
 117     rvec             *axis3;
 118     real             *twist, *residuetwist;
 119     real             *radius, *residueradius;
 120     real             *rise, *residuerise;
 121     real             *residuebending;
 122
 123     real              tmp;
 124     real              weight[3];
 125     t_pbc             pbc;
 126     matrix            A;
 127
 128     FILE             *fpaxis, *fpcenter, *fptilt, *fprotation;
 129     FILE             *fpradius, *fprise, *fptwist;
 130     FILE             *fptheta1, *fptheta2, *fptheta3;
 131     FILE             *fpbending;
 132     int               ePBC;
 133
 134     gmx_output_env_t *oenv;
 135     gmx_rmpbc_t       gpbc = nullptr;
 136
 137     static  gmx_bool  bSC          = FALSE;
 138     static gmx_bool   bIncremental = FALSE;
 139
 140     static t_pargs    pa[] = {
 141         { "-sidechain",      FALSE, etBOOL, {&bSC},
 142           "Calculate sidechain directions relative to helix axis too." },
 143         { "-incremental",        FALSE, etBOOL, {&bIncremental},
 144           "Calculate incremental rather than total rotation/tilt." },
 145     };
 146 #define NPA asize(pa)
 147
 148     t_filenm fnm[] = {
 149         { efTPR, nullptr, nullptr, ffREAD },
 150         { efTRX, "-f", nullptr, ffREAD },
 151         { efNDX, nullptr, nullptr, ffOPTRD },
 152         { efDAT, "-oaxis",    "helixaxis", ffWRITE },
 153         { efDAT, "-ocenter",  "center", ffWRITE },
 154         { efXVG, "-orise",    "rise", ffWRITE },
 155         { efXVG, "-oradius",  "radius", ffWRITE },
 156         { efXVG, "-otwist",   "twist", ffWRITE },
 157         { efXVG, "-obending", "bending", ffWRITE },
 158         { efXVG, "-otilt",    "tilt", ffWRITE },
 159         { efXVG, "-orot",     "rotation", ffWRITE }
 160     };
 161 #define NFILE asize(fnm)
 162
 163     if (!parse_common_args(&argc, argv, PCA_CAN_TIME,
 164                            NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
 165     {
 166         return 0;
 167     }
 168
 169     top = read_top(ftp2fn(efTPR, 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    = gmx_ffopen(opt2fn("-oaxis", NFILE, fnm), "w");
 221     fpcenter  = gmx_ffopen(opt2fn("-ocenter", NFILE, fnm), "w");
 222     fprise    = gmx_ffopen(opt2fn("-orise", NFILE, fnm), "w");
 223     fpradius  = gmx_ffopen(opt2fn("-oradius", NFILE, fnm), "w");
 224     fptwist   = gmx_ffopen(opt2fn("-otwist", NFILE, fnm), "w");
 225     fpbending = gmx_ffopen(opt2fn("-obending", NFILE, fnm), "w");
 226
 227     fptheta1 = gmx_ffopen("theta1.xvg", "w");
 228     fptheta2 = gmx_ffopen("theta2.xvg", "w");
 229     fptheta3 = gmx_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);
 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 * std::acos( tmp );
 286             radius[i] = std::sqrt( norm(diff13[i])*norm(diff24[i]) ) / (2.0* (1.0-tmp) );
 287             rise[i]   = std::abs(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*std::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         }
 371         fprintf(fprise, "\n");
 372         fprintf(fpradius, "\n");
 373         fprintf(fpaxis, "\n");
 374         fprintf(fpcenter, "\n");
 375         fprintf(fptwist, "\n");
 376         fprintf(fpbending, "\n");
 377
 378         if (teller == 0)
 379         {
 380             for (i = 0; i < iCA; i++)
 381             {
 382                 copy_rvec(residuehelixaxis[i], residuehelixaxis_t0[i]);
 383                 copy_rvec(residuevector[i], residuevector_t0[i]);
 384                 copy_rvec(axis3[i], axis3_t0[i]);
 385             }
 386         }
 387         else
 388         {
 389             fprintf(fptilt, "%15.12g       ", t);
 390             fprintf(fprotation, "%15.12g       ", t);
 391             fprintf(fptheta1, "%15.12g      ", t);
 392             fprintf(fptheta2, "%15.12g      ", t);
 393             fprintf(fptheta3, "%15.12g      ", t);
 394
 395             for (i = 0; i < iCA; i++)
 396             {
 397                 if (i == 0 || i == iCA-1)
 398                 {
 399                     tilt = rotation = 0;
 400                 }
 401                 else
 402                 {
 403                     if (!bIncremental)
 404                     {
 405                         /* Total rotation & tilt */
 406                         copy_rvec(residuehelixaxis_t0[i], refaxes[0]);
 407                         copy_rvec(residuevector_t0[i], refaxes[1]);
 408                         copy_rvec(axis3_t0[i], refaxes[2]);
 409                     }
 410                     else
 411                     {
 412                         /* Rotation/tilt since last step */
 413                         copy_rvec(residuehelixaxis_tlast[i], refaxes[0]);
 414                         copy_rvec(residuevector_tlast[i], refaxes[1]);
 415                         copy_rvec(axis3_tlast[i], refaxes[2]);
 416                     }
 417                     copy_rvec(residuehelixaxis[i], newaxes[0]);
 418                     copy_rvec(residuevector[i], newaxes[1]);
 419                     copy_rvec(axis3[i], newaxes[2]);
 420
 421                     /* rotate reference frame onto unit axes */
 422                     calc_fit_R(3, 3, weight, unitaxes, refaxes, A);
 423                     for (j = 0; j < 3; j++)
 424                     {
 425                         mvmul(A, refaxes[j], rot_refaxes[j]);
 426                         mvmul(A, newaxes[j], rot_newaxes[j]);
 427                     }
 428
 429                     /* Determine local rotation matrix A */
 430                     calc_fit_R(3, 3, weight, rot_newaxes, rot_refaxes, A);
 431                     /* Calculate euler angles, from rotation order y-z-x, where
 432                      * x is helixaxis, y residuevector, and z axis3.
 433                      *
 434                      * A contains rotation column vectors.
 435                      */
 436
 437                     theta1 = 180.0/M_PI*std::atan2(A[0][2], A[0][0]);
 438                     theta2 = 180.0/M_PI*std::asin(-A[0][1]);
 439                     theta3 = 180.0/M_PI*std::atan2(A[2][1], A[1][1]);
 440
 441                     tilt     = std::sqrt(theta1*theta1+theta2*theta2);
 442                     rotation = theta3;
 443                     fprintf(fptheta1, "%15.12g  ", theta1);
 444                     fprintf(fptheta2, "%15.12g  ", theta2);
 445                     fprintf(fptheta3, "%15.12g  ", theta3);
 446
 447                 }
 448                 fprintf(fptilt, "%15.12g  ", tilt);
 449                 fprintf(fprotation, "%15.12g  ", rotation);
 450             }
 451             fprintf(fptilt, "\n");
 452             fprintf(fprotation, "\n");
 453             fprintf(fptheta1, "\n");
 454             fprintf(fptheta2, "\n");
 455             fprintf(fptheta3, "\n");
 456         }
 457
 458         for (i = 0; i < iCA; i++)
 459         {
 460             copy_rvec(residuehelixaxis[i], residuehelixaxis_tlast[i]);
 461             copy_rvec(residuevector[i], residuevector_tlast[i]);
 462             copy_rvec(axis3[i], axis3_tlast[i]);
 463         }
 464
 465         teller++;
 466     }
 467     while (read_next_x(oenv, status, &t, x, box));
 468
 469     gmx_rmpbc_done(gpbc);
 470
 471     gmx_ffclose(fpaxis);
 472     gmx_ffclose(fpcenter);
 473     xvgrclose(fptilt);
 474     xvgrclose(fprotation);
 475     gmx_ffclose(fprise);
 476     gmx_ffclose(fpradius);
 477     gmx_ffclose(fptwist);
 478     gmx_ffclose(fpbending);
 479     gmx_ffclose(fptheta1);
 480     gmx_ffclose(fptheta2);
 481     gmx_ffclose(fptheta3);
 482
 483     close_trj(status);
 484
 485     return 0;
 486 }