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,2019,2020,2021, 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/units.h"
  48 #include "gromacs/math/utilities.h"
  49 #include "gromacs/math/vec.h"
  50 #include "gromacs/pbcutil/pbc.h"
  51 #include "gromacs/pbcutil/rmpbc.h"
  52 #include "gromacs/topology/index.h"
  53 #include "gromacs/topology/topology.h"
  54 #include "gromacs/utility/arraysize.h"
  55 #include "gromacs/utility/fatalerror.h"
  56 #include "gromacs/utility/futil.h"
  57 #include "gromacs/utility/smalloc.h"
  58
  59 int gmx_helixorient(int argc, char* argv[])
  60 {
  61     const char* desc[] = {
  62         "[THISMODULE] calculates the coordinates and direction of the average",
  63         "axis inside an alpha helix, and the direction/vectors of both the",
  64         "C[GRK]alpha[grk] and (optionally) a sidechain atom relative to the axis.[PAR]",
  65         "As input, you need to specify an index group with C[GRK]alpha[grk] atoms",
  66         "corresponding to an [GRK]alpha[grk]-helix of continuous residues. Sidechain",
  67         "directions require a second index group of the same size, containing",
  68         "the heavy atom in each residue that should represent the sidechain.[PAR]",
  69         "[BB]Note[bb] that this program does not do any fitting of structures.[PAR]",
  70         "We need four C[GRK]alpha[grk] coordinates to define the local direction of the helix",
  71         "axis.[PAR]",
  72         "The tilt/rotation is calculated from Euler rotations, where we define",
  73         "the helix axis as the local [IT]x[it]-axis, the residues/C[GRK]alpha[grk] vector as ",
  74         "[IT]y[it], and the",
  75         "[IT]z[it]-axis from their cross product. We use the Euler Y-Z-X rotation, meaning",
  76         "we first tilt the helix axis (1) around and (2) orthogonal to the residues",
  77         "vector, and finally apply the (3) rotation around it. For debugging or other",
  78         "purposes, we also write out the actual Euler rotation angles as [TT]theta[1-3].xvg[tt]"
  79     };
  80
  81     t_topology*  top = nullptr;
  82     real         t;
  83     rvec*        x = nullptr;
  84     matrix       box;
  85     t_trxstatus* status;
  86     int          natoms;
  87     real         theta1, theta2, theta3;
  88
  89     int   i, j, teller = 0;
  90     int   iCA, iSC;
  91     int*  ind_CA;
  92     int*  ind_SC;
  93     char* gn_CA;
  94     char* gn_SC;
  95     rvec  v1, v2;
  96     rvec *x_CA, *x_SC;
  97     rvec* r12;
  98     rvec* r23;
  99     rvec* r34;
 100     rvec* diff13;
 101     rvec* diff24;
 102     rvec* helixaxis;
 103     rvec* residuehelixaxis;
 104     rvec* residueorigin;
 105     rvec* residuevector;
 106     rvec* sidechainvector;
 107
 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;
 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     PbcType pbcType;
 135
 136     gmx_output_env_t* oenv;
 137     gmx_rmpbc_t       gpbc = nullptr;
 138
 139     static gmx_bool bSC          = FALSE;
 140     static gmx_bool bIncremental = FALSE;
 141
 142     static t_pargs pa[] = {
 143         { "-sidechain",
 144           FALSE,
 145           etBOOL,
 146           { &bSC },
 147           "Calculate sidechain directions relative to helix axis too." },
 148         { "-incremental",
 149           FALSE,
 150           etBOOL,
 151           { &bIncremental },
 152           "Calculate incremental rather than total rotation/tilt." },
 153     };
 154 #define NPA asize(pa)
 155
 156     t_filenm fnm[] = {
 157         { efTPR, nullptr, nullptr, ffREAD },        { efTRX, "-f", nullptr, ffREAD },
 158         { efNDX, nullptr, nullptr, ffOPTRD },       { efDAT, "-oaxis", "helixaxis", ffWRITE },
 159         { efDAT, "-ocenter", "center", ffWRITE },   { efXVG, "-orise", "rise", ffWRITE },
 160         { efXVG, "-oradius", "radius", ffWRITE },   { efXVG, "-otwist", "twist", ffWRITE },
 161         { efXVG, "-obending", "bending", ffWRITE }, { efXVG, "-otilt", "tilt", ffWRITE },
 162         { efXVG, "-orot", "rotation", ffWRITE }
 163     };
 164 #define NFILE asize(fnm)
 165
 166     if (!parse_common_args(
 167                 &argc, argv, PCA_CAN_TIME, NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
 168     {
 169         return 0;
 170     }
 171
 172     top = read_top(ftp2fn(efTPR, NFILE, fnm), &pbcType);
 173
 174     for (i = 0; i < 3; i++)
 175     {
 176         weight[i] = 1.0;
 177     }
 178
 179     /* read index files */
 180     printf("Select a group of Calpha atoms corresponding to a single continuous helix:\n");
 181     get_index(&(top->atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &iCA, &ind_CA, &gn_CA);
 182     snew(x_CA, iCA);
 183     snew(x_SC, iCA); /* sic! */
 184
 185     snew(r12, iCA - 3);
 186     snew(r23, iCA - 3);
 187     snew(r34, iCA - 3);
 188     snew(diff13, iCA - 3);
 189     snew(diff24, iCA - 3);
 190     snew(helixaxis, iCA - 3);
 191     snew(twist, iCA);
 192     snew(residuetwist, iCA);
 193     snew(radius, iCA);
 194     snew(residueradius, iCA);
 195     snew(rise, iCA);
 196     snew(residuerise, iCA);
 197     snew(residueorigin, iCA);
 198     snew(residuehelixaxis, iCA);
 199     snew(residuevector, iCA);
 200     snew(sidechainvector, iCA);
 201     snew(residuebending, iCA);
 202     snew(residuehelixaxis_t0, iCA);
 203     snew(residuevector_t0, iCA);
 204     snew(axis3_t0, iCA);
 205     snew(residuehelixaxis_tlast, iCA);
 206     snew(residuevector_tlast, iCA);
 207     snew(axis3_tlast, iCA);
 208     snew(axis3, iCA);
 209
 210     if (bSC)
 211     {
 212         printf("Select a group of atoms defining the sidechain direction (1/residue):\n");
 213         get_index(&(top->atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &iSC, &ind_SC, &gn_SC);
 214         if (iSC != iCA)
 215         {
 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    = gmx_ffopen(opt2fn("-oaxis", NFILE, fnm), "w");
 223     fpcenter  = gmx_ffopen(opt2fn("-ocenter", NFILE, fnm), "w");
 224     fprise    = gmx_ffopen(opt2fn("-orise", NFILE, fnm), "w");
 225     fpradius  = gmx_ffopen(opt2fn("-oradius", NFILE, fnm), "w");
 226     fptwist   = gmx_ffopen(opt2fn("-otwist", NFILE, fnm), "w");
 227     fpbending = gmx_ffopen(opt2fn("-obending", NFILE, fnm), "w");
 228
 229     fptheta1 = gmx_ffopen("theta1.xvg", "w");
 230     fptheta2 = gmx_ffopen("theta2.xvg", "w");
 231     fptheta3 = gmx_ffopen("theta3.xvg", "w");
 232
 233     if (bIncremental)
 234     {
 235         fptilt     = xvgropen(opt2fn("-otilt", NFILE, fnm),
 236                           "Incremental local helix tilt",
 237                           "Time(ps)",
 238                           "Tilt (degrees)",
 239                           oenv);
 240         fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
 241                               "Incremental local helix rotation",
 242                               "Time(ps)",
 243                               "Rotation (degrees)",
 244                               oenv);
 245     }
 246     else
 247     {
 248         fptilt     = xvgropen(opt2fn("-otilt", NFILE, fnm),
 249                           "Cumulative local helix tilt",
 250                           "Time(ps)",
 251                           "Tilt (degrees)",
 252                           oenv);
 253         fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
 254                               "Cumulative local helix rotation",
 255                               "Time(ps)",
 256                               "Rotation (degrees)",
 257                               oenv);
 258     }
 259
 260     clear_rvecs(3, unitaxes);
 261     unitaxes[0][0] = 1;
 262     unitaxes[1][1] = 1;
 263     unitaxes[2][2] = 1;
 264
 265     gpbc = gmx_rmpbc_init(&top->idef, pbcType, natoms);
 266
 267     do
 268     {
 269         /* initialisation for correct distance calculations */
 270         set_pbc(&pbc, pbcType, box);
 271         /* make molecules whole again */
 272         gmx_rmpbc(gpbc, natoms, box, x);
 273
 274         /* copy coords to our smaller arrays */
 275         for (i = 0; i < iCA; i++)
 276         {
 277             copy_rvec(x[ind_CA[i]], x_CA[i]);
 278             if (bSC)
 279             {
 280                 copy_rvec(x[ind_SC[i]], x_SC[i]);
 281             }
 282         }
 283
 284         for (i = 0; i < iCA - 3; i++)
 285         {
 286             rvec_sub(x_CA[i + 1], x_CA[i], r12[i]);
 287             rvec_sub(x_CA[i + 2], x_CA[i + 1], r23[i]);
 288             rvec_sub(x_CA[i + 3], x_CA[i + 2], r34[i]);
 289             rvec_sub(r12[i], r23[i], diff13[i]);
 290             rvec_sub(r23[i], r34[i], diff24[i]);
 291             /* calculate helix axis */
 292             cprod(diff13[i], diff24[i], helixaxis[i]);
 293             svmul(1.0 / norm(helixaxis[i]), helixaxis[i], helixaxis[i]);
 294
 295             tmp       = cos_angle(diff13[i], diff24[i]);
 296             twist[i]  = 180.0 / M_PI * std::acos(tmp);
 297             radius[i] = std::sqrt(norm(diff13[i]) * norm(diff24[i])) / (2.0 * (1.0 - tmp));
 298             rise[i]   = std::abs(iprod(r23[i], helixaxis[i]));
 299
 300             svmul(radius[i] / norm(diff13[i]), diff13[i], v1);
 301             svmul(radius[i] / norm(diff24[i]), diff24[i], v2);
 302
 303             rvec_sub(x_CA[i + 1], v1, residueorigin[i + 1]);
 304             rvec_sub(x_CA[i + 2], v2, residueorigin[i + 2]);
 305         }
 306         residueradius[0] = residuetwist[0] = residuerise[0] = 0;
 307
 308         residueradius[1] = radius[0];
 309         residuetwist[1]  = twist[0];
 310         residuerise[1]   = rise[0];
 311
 312         residuebending[0] = residuebending[1] = 0;
 313         for (i = 2; i < iCA - 2; i++)
 314         {
 315             residueradius[i] = 0.5 * (radius[i - 2] + radius[i - 1]);
 316             residuetwist[i]  = 0.5 * (twist[i - 2] + twist[i - 1]);
 317             residuerise[i]   = 0.5 * (rise[i - 2] + rise[i - 1]);
 318             residuebending[i] = 180.0 / M_PI * std::acos(cos_angle(helixaxis[i - 2], helixaxis[i - 1]));
 319         }
 320         residueradius[iCA - 2] = radius[iCA - 4];
 321         residuetwist[iCA - 2]  = twist[iCA - 4];
 322         residuerise[iCA - 2]   = rise[iCA - 4];
 323         residueradius[iCA - 1] = residuetwist[iCA - 1] = residuerise[iCA - 1] = 0;
 324         residuebending[iCA - 2] = residuebending[iCA - 1] = 0;
 325
 326         clear_rvec(residueorigin[0]);
 327         clear_rvec(residueorigin[iCA - 1]);
 328
 329         /* average helix axes to define them on the residues.
 330          * Just extrapolate second first/list atom.
 331          */
 332         copy_rvec(helixaxis[0], residuehelixaxis[0]);
 333         copy_rvec(helixaxis[0], residuehelixaxis[1]);
 334
 335         for (i = 2; i < iCA - 2; i++)
 336         {
 337             rvec_add(helixaxis[i - 2], helixaxis[i - 1], residuehelixaxis[i]);
 338             svmul(0.5, residuehelixaxis[i], residuehelixaxis[i]);
 339         }
 340         copy_rvec(helixaxis[iCA - 4], residuehelixaxis[iCA - 2]);
 341         copy_rvec(helixaxis[iCA - 4], residuehelixaxis[iCA - 1]);
 342
 343         /* Normalize the axis */
 344         for (i = 0; i < iCA; i++)
 345         {
 346             svmul(1.0 / norm(residuehelixaxis[i]), residuehelixaxis[i], residuehelixaxis[i]);
 347         }
 348
 349         /* calculate vector from origin to residue CA */
 350         fprintf(fpaxis, "%15.12g  ", t);
 351         fprintf(fpcenter, "%15.12g  ", t);
 352         fprintf(fprise, "%15.12g  ", t);
 353         fprintf(fpradius, "%15.12g  ", t);
 354         fprintf(fptwist, "%15.12g  ", t);
 355         fprintf(fpbending, "%15.12g  ", t);
 356
 357         for (i = 0; i < iCA; i++)
 358         {
 359             if (i == 0 || i == iCA - 1)
 360             {
 361                 fprintf(fpaxis, "%15.12g %15.12g %15.12g       ", 0.0, 0.0, 0.0);
 362                 fprintf(fpcenter, "%15.12g %15.12g %15.12g       ", 0.0, 0.0, 0.0);
 363                 fprintf(fprise, "%15.12g  ", 0.0);
 364                 fprintf(fpradius, "%15.12g  ", 0.0);
 365                 fprintf(fptwist, "%15.12g  ", 0.0);
 366                 fprintf(fpbending, "%15.12g  ", 0.0);
 367             }
 368             else
 369             {
 370                 rvec_sub(bSC ? x_SC[i] : x_CA[i], residueorigin[i], residuevector[i]);
 371                 svmul(1.0 / norm(residuevector[i]), residuevector[i], residuevector[i]);
 372                 cprod(residuehelixaxis[i], residuevector[i], axis3[i]);
 373                 fprintf(fpaxis,
 374                         "%15.12g %15.12g %15.12g       ",
 375                         residuehelixaxis[i][0],
 376                         residuehelixaxis[i][1],
 377                         residuehelixaxis[i][2]);
 378                 fprintf(fpcenter,
 379                         "%15.12g %15.12g %15.12g       ",
 380                         residueorigin[i][0],
 381                         residueorigin[i][1],
 382                         residueorigin[i][2]);
 383
 384                 fprintf(fprise, "%15.12g  ", residuerise[i]);
 385                 fprintf(fpradius, "%15.12g  ", residueradius[i]);
 386                 fprintf(fptwist, "%15.12g  ", residuetwist[i]);
 387                 fprintf(fpbending, "%15.12g  ", residuebending[i]);
 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                     /* rotate reference frame onto unit axes */
 441                     calc_fit_R(3, 3, weight, unitaxes, refaxes, A);
 442                     for (j = 0; j < 3; j++)
 443                     {
 444                         mvmul(A, refaxes[j], rot_refaxes[j]);
 445                         mvmul(A, newaxes[j], rot_newaxes[j]);
 446                     }
 447
 448                     /* Determine local rotation matrix A */
 449                     calc_fit_R(3, 3, weight, rot_newaxes, rot_refaxes, A);
 450                     /* Calculate euler angles, from rotation order y-z-x, where
 451                      * x is helixaxis, y residuevector, and z axis3.
 452                      *
 453                      * A contains rotation column vectors.
 454                      */
 455
 456                     theta1 = 180.0 / M_PI * std::atan2(A[0][2], A[0][0]);
 457                     theta2 = 180.0 / M_PI * std::asin(-A[0][1]);
 458                     theta3 = 180.0 / M_PI * std::atan2(A[2][1], A[1][1]);
 459
 460                     tilt     = std::sqrt(theta1 * theta1 + theta2 * theta2);
 461                     rotation = theta3;
 462                     fprintf(fptheta1, "%15.12g  ", theta1);
 463                     fprintf(fptheta2, "%15.12g  ", theta2);
 464                     fprintf(fptheta3, "%15.12g  ", theta3);
 465                 }
 466                 fprintf(fptilt, "%15.12g  ", tilt);
 467                 fprintf(fprotation, "%15.12g  ", rotation);
 468             }
 469             fprintf(fptilt, "\n");
 470             fprintf(fprotation, "\n");
 471             fprintf(fptheta1, "\n");
 472             fprintf(fptheta2, "\n");
 473             fprintf(fptheta3, "\n");
 474         }
 475
 476         for (i = 0; i < iCA; i++)
 477         {
 478             copy_rvec(residuehelixaxis[i], residuehelixaxis_tlast[i]);
 479             copy_rvec(residuevector[i], residuevector_tlast[i]);
 480             copy_rvec(axis3[i], axis3_tlast[i]);
 481         }
 482
 483         teller++;
 484     } while (read_next_x(oenv, status, &t, x, box));
 485
 486     gmx_rmpbc_done(gpbc);
 487
 488     gmx_ffclose(fpaxis);
 489     gmx_ffclose(fpcenter);
 490     xvgrclose(fptilt);
 491     xvgrclose(fprotation);
 492     gmx_ffclose(fprise);
 493     gmx_ffclose(fpradius);
 494     gmx_ffclose(fptwist);
 495     gmx_ffclose(fpbending);
 496     gmx_ffclose(fptheta1);
 497     gmx_ffclose(fptheta2);
 498     gmx_ffclose(fptheta3);
 499
 500     close_trx(status);
 501
 502     return 0;
 503 }