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, 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 ",
  73         "[IT]y[it], and the",
  74         "[IT]z[it]-axis from their cross product. We use the Euler Y-Z-X rotation, meaning",
  75         "we first tilt the helix axis (1) around and (2) orthogonal to the residues",
  76         "vector, and finally apply the (3) rotation around it. For debugging or other",
  77         "purposes, we also write out the actual Euler rotation angles as [TT]theta[1-3].xvg[tt]"
  78     };
  79
  80     t_topology*  top = nullptr;
  81     real         t;
  82     rvec*        x = nullptr;
  83     matrix       box;
  84     t_trxstatus* status;
  85     int          natoms;
  86     real         theta1, theta2, theta3;
  87
  88     int   i, j, teller = 0;
  89     int   iCA, iSC;
  90     int*  ind_CA;
  91     int*  ind_SC;
  92     char* gn_CA;
  93     char* gn_SC;
  94     rvec  v1, v2;
  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* residuehelixaxis_t0;
 108     rvec* residuevector_t0;
 109     rvec* axis3_t0;
 110     rvec* residuehelixaxis_tlast;
 111     rvec* residuevector_tlast;
 112     rvec* axis3_tlast;
 113     rvec  refaxes[3], newaxes[3];
 114     rvec  unitaxes[3];
 115     rvec  rot_refaxes[3], rot_newaxes[3];
 116
 117     real  tilt, rotation;
 118     rvec* axis3;
 119     real *twist, *residuetwist;
 120     real *radius, *residueradius;
 121     real *rise, *residuerise;
 122     real* residuebending;
 123
 124     real   tmp;
 125     real   weight[3];
 126     t_pbc  pbc;
 127     matrix A;
 128
 129     FILE *  fpaxis, *fpcenter, *fptilt, *fprotation;
 130     FILE *  fpradius, *fprise, *fptwist;
 131     FILE *  fptheta1, *fptheta2, *fptheta3;
 132     FILE*   fpbending;
 133     PbcType pbcType;
 134
 135     gmx_output_env_t* oenv;
 136     gmx_rmpbc_t       gpbc = nullptr;
 137
 138     static gmx_bool bSC          = FALSE;
 139     static gmx_bool bIncremental = FALSE;
 140
 141     static t_pargs pa[] = {
 142         { "-sidechain",
 143           FALSE,
 144           etBOOL,
 145           { &bSC },
 146           "Calculate sidechain directions relative to helix axis too." },
 147         { "-incremental",
 148           FALSE,
 149           etBOOL,
 150           { &bIncremental },
 151           "Calculate incremental rather than total rotation/tilt." },
 152     };
 153 #define NPA asize(pa)
 154
 155     t_filenm fnm[] = {
 156         { efTPR, nullptr, nullptr, ffREAD },        { efTRX, "-f", nullptr, ffREAD },
 157         { efNDX, nullptr, nullptr, ffOPTRD },       { efDAT, "-oaxis", "helixaxis", ffWRITE },
 158         { efDAT, "-ocenter", "center", ffWRITE },   { efXVG, "-orise", "rise", ffWRITE },
 159         { efXVG, "-oradius", "radius", ffWRITE },   { efXVG, "-otwist", "twist", ffWRITE },
 160         { efXVG, "-obending", "bending", ffWRITE }, { efXVG, "-otilt", "tilt", ffWRITE },
 161         { efXVG, "-orot", "rotation", ffWRITE }
 162     };
 163 #define NFILE asize(fnm)
 164
 165     if (!parse_common_args(
 166                 &argc, argv, PCA_CAN_TIME, NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
 167     {
 168         return 0;
 169     }
 170
 171     top = read_top(ftp2fn(efTPR, NFILE, fnm), &pbcType);
 172
 173     for (i = 0; i < 3; i++)
 174     {
 175         weight[i] = 1.0;
 176     }
 177
 178     /* read index files */
 179     printf("Select a group of Calpha atoms corresponding to a single continuous helix:\n");
 180     get_index(&(top->atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &iCA, &ind_CA, &gn_CA);
 181     snew(x_CA, iCA);
 182     snew(x_SC, iCA); /* sic! */
 183
 184     snew(r12, iCA - 3);
 185     snew(r23, iCA - 3);
 186     snew(r34, iCA - 3);
 187     snew(diff13, iCA - 3);
 188     snew(diff24, iCA - 3);
 189     snew(helixaxis, iCA - 3);
 190     snew(twist, iCA);
 191     snew(residuetwist, iCA);
 192     snew(radius, iCA);
 193     snew(residueradius, iCA);
 194     snew(rise, iCA);
 195     snew(residuerise, iCA);
 196     snew(residueorigin, iCA);
 197     snew(residuehelixaxis, iCA);
 198     snew(residuevector, iCA);
 199     snew(sidechainvector, iCA);
 200     snew(residuebending, iCA);
 201     snew(residuehelixaxis_t0, iCA);
 202     snew(residuevector_t0, iCA);
 203     snew(axis3_t0, iCA);
 204     snew(residuehelixaxis_tlast, iCA);
 205     snew(residuevector_tlast, iCA);
 206     snew(axis3_tlast, iCA);
 207     snew(axis3, iCA);
 208
 209     if (bSC)
 210     {
 211         printf("Select a group of atoms defining the sidechain direction (1/residue):\n");
 212         get_index(&(top->atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &iSC, &ind_SC, &gn_SC);
 213         if (iSC != iCA)
 214         {
 215             gmx_fatal(FARGS, "Number of sidechain atoms (%d) != number of CA atoms (%d)", iSC, iCA);
 216         }
 217     }
 218
 219     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 220
 221     fpaxis    = gmx_ffopen(opt2fn("-oaxis", NFILE, fnm), "w");
 222     fpcenter  = gmx_ffopen(opt2fn("-ocenter", NFILE, fnm), "w");
 223     fprise    = gmx_ffopen(opt2fn("-orise", NFILE, fnm), "w");
 224     fpradius  = gmx_ffopen(opt2fn("-oradius", NFILE, fnm), "w");
 225     fptwist   = gmx_ffopen(opt2fn("-otwist", NFILE, fnm), "w");
 226     fpbending = gmx_ffopen(opt2fn("-obending", NFILE, fnm), "w");
 227
 228     fptheta1 = gmx_ffopen("theta1.xvg", "w");
 229     fptheta2 = gmx_ffopen("theta2.xvg", "w");
 230     fptheta3 = gmx_ffopen("theta3.xvg", "w");
 231
 232     if (bIncremental)
 233     {
 234         fptilt     = xvgropen(opt2fn("-otilt", NFILE, fnm),
 235                           "Incremental local helix tilt",
 236                           "Time(ps)",
 237                           "Tilt (degrees)",
 238                           oenv);
 239         fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
 240                               "Incremental local helix rotation",
 241                               "Time(ps)",
 242                               "Rotation (degrees)",
 243                               oenv);
 244     }
 245     else
 246     {
 247         fptilt     = xvgropen(opt2fn("-otilt", NFILE, fnm),
 248                           "Cumulative local helix tilt",
 249                           "Time(ps)",
 250                           "Tilt (degrees)",
 251                           oenv);
 252         fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
 253                               "Cumulative local helix rotation",
 254                               "Time(ps)",
 255                               "Rotation (degrees)",
 256                               oenv);
 257     }
 258
 259     clear_rvecs(3, unitaxes);
 260     unitaxes[0][0] = 1;
 261     unitaxes[1][1] = 1;
 262     unitaxes[2][2] = 1;
 263
 264     gpbc = gmx_rmpbc_init(&top->idef, pbcType, natoms);
 265
 266     do
 267     {
 268         /* initialisation for correct distance calculations */
 269         set_pbc(&pbc, pbcType, box);
 270         /* make molecules whole again */
 271         gmx_rmpbc(gpbc, natoms, box, x);
 272
 273         /* copy coords to our smaller arrays */
 274         for (i = 0; i < iCA; i++)
 275         {
 276             copy_rvec(x[ind_CA[i]], x_CA[i]);
 277             if (bSC)
 278             {
 279                 copy_rvec(x[ind_SC[i]], x_SC[i]);
 280             }
 281         }
 282
 283         for (i = 0; i < iCA - 3; i++)
 284         {
 285             rvec_sub(x_CA[i + 1], x_CA[i], r12[i]);
 286             rvec_sub(x_CA[i + 2], x_CA[i + 1], r23[i]);
 287             rvec_sub(x_CA[i + 3], x_CA[i + 2], r34[i]);
 288             rvec_sub(r12[i], r23[i], diff13[i]);
 289             rvec_sub(r23[i], r34[i], diff24[i]);
 290             /* calculate helix axis */
 291             cprod(diff13[i], diff24[i], helixaxis[i]);
 292             svmul(1.0 / norm(helixaxis[i]), helixaxis[i], helixaxis[i]);
 293
 294             tmp       = cos_angle(diff13[i], diff24[i]);
 295             twist[i]  = 180.0 / M_PI * std::acos(tmp);
 296             radius[i] = std::sqrt(norm(diff13[i]) * norm(diff24[i])) / (2.0 * (1.0 - tmp));
 297             rise[i]   = std::abs(iprod(r23[i], helixaxis[i]));
 298
 299             svmul(radius[i] / norm(diff13[i]), diff13[i], v1);
 300             svmul(radius[i] / norm(diff24[i]), diff24[i], v2);
 301
 302             rvec_sub(x_CA[i + 1], v1, residueorigin[i + 1]);
 303             rvec_sub(x_CA[i + 2], v2, residueorigin[i + 2]);
 304         }
 305         residueradius[0] = residuetwist[0] = residuerise[0] = 0;
 306
 307         residueradius[1] = radius[0];
 308         residuetwist[1]  = twist[0];
 309         residuerise[1]   = rise[0];
 310
 311         residuebending[0] = residuebending[1] = 0;
 312         for (i = 2; i < iCA - 2; i++)
 313         {
 314             residueradius[i] = 0.5 * (radius[i - 2] + radius[i - 1]);
 315             residuetwist[i]  = 0.5 * (twist[i - 2] + twist[i - 1]);
 316             residuerise[i]   = 0.5 * (rise[i - 2] + rise[i - 1]);
 317             residuebending[i] = 180.0 / M_PI * std::acos(cos_angle(helixaxis[i - 2], helixaxis[i - 1]));
 318         }
 319         residueradius[iCA - 2] = radius[iCA - 4];
 320         residuetwist[iCA - 2]  = twist[iCA - 4];
 321         residuerise[iCA - 2]   = rise[iCA - 4];
 322         residueradius[iCA - 1] = residuetwist[iCA - 1] = residuerise[iCA - 1] = 0;
 323         residuebending[iCA - 2] = residuebending[iCA - 1] = 0;
 324
 325         clear_rvec(residueorigin[0]);
 326         clear_rvec(residueorigin[iCA - 1]);
 327
 328         /* average helix axes to define them on the residues.
 329          * Just extrapolate second first/list atom.
 330          */
 331         copy_rvec(helixaxis[0], residuehelixaxis[0]);
 332         copy_rvec(helixaxis[0], residuehelixaxis[1]);
 333
 334         for (i = 2; i < iCA - 2; i++)
 335         {
 336             rvec_add(helixaxis[i - 2], helixaxis[i - 1], residuehelixaxis[i]);
 337             svmul(0.5, residuehelixaxis[i], residuehelixaxis[i]);
 338         }
 339         copy_rvec(helixaxis[iCA - 4], residuehelixaxis[iCA - 2]);
 340         copy_rvec(helixaxis[iCA - 4], residuehelixaxis[iCA - 1]);
 341
 342         /* Normalize the axis */
 343         for (i = 0; i < iCA; i++)
 344         {
 345             svmul(1.0 / norm(residuehelixaxis[i]), residuehelixaxis[i], residuehelixaxis[i]);
 346         }
 347
 348         /* calculate vector from origin to residue CA */
 349         fprintf(fpaxis, "%15.12g  ", t);
 350         fprintf(fpcenter, "%15.12g  ", t);
 351         fprintf(fprise, "%15.12g  ", t);
 352         fprintf(fpradius, "%15.12g  ", t);
 353         fprintf(fptwist, "%15.12g  ", t);
 354         fprintf(fpbending, "%15.12g  ", t);
 355
 356         for (i = 0; i < iCA; i++)
 357         {
 358             if (i == 0 || i == iCA - 1)
 359             {
 360                 fprintf(fpaxis, "%15.12g %15.12g %15.12g       ", 0.0, 0.0, 0.0);
 361                 fprintf(fpcenter, "%15.12g %15.12g %15.12g       ", 0.0, 0.0, 0.0);
 362                 fprintf(fprise, "%15.12g  ", 0.0);
 363                 fprintf(fpradius, "%15.12g  ", 0.0);
 364                 fprintf(fptwist, "%15.12g  ", 0.0);
 365                 fprintf(fpbending, "%15.12g  ", 0.0);
 366             }
 367             else
 368             {
 369                 rvec_sub(bSC ? x_SC[i] : x_CA[i], residueorigin[i], residuevector[i]);
 370                 svmul(1.0 / norm(residuevector[i]), residuevector[i], residuevector[i]);
 371                 cprod(residuehelixaxis[i], residuevector[i], axis3[i]);
 372                 fprintf(fpaxis,
 373                         "%15.12g %15.12g %15.12g       ",
 374                         residuehelixaxis[i][0],
 375                         residuehelixaxis[i][1],
 376                         residuehelixaxis[i][2]);
 377                 fprintf(fpcenter,
 378                         "%15.12g %15.12g %15.12g       ",
 379                         residueorigin[i][0],
 380                         residueorigin[i][1],
 381                         residueorigin[i][2]);
 382
 383                 fprintf(fprise, "%15.12g  ", residuerise[i]);
 384                 fprintf(fpradius, "%15.12g  ", residueradius[i]);
 385                 fprintf(fptwist, "%15.12g  ", residuetwist[i]);
 386                 fprintf(fpbending, "%15.12g  ", residuebending[i]);
 387             }
 388         }
 389         fprintf(fprise, "\n");
 390         fprintf(fpradius, "\n");
 391         fprintf(fpaxis, "\n");
 392         fprintf(fpcenter, "\n");
 393         fprintf(fptwist, "\n");
 394         fprintf(fpbending, "\n");
 395
 396         if (teller == 0)
 397         {
 398             for (i = 0; i < iCA; i++)
 399             {
 400                 copy_rvec(residuehelixaxis[i], residuehelixaxis_t0[i]);
 401                 copy_rvec(residuevector[i], residuevector_t0[i]);
 402                 copy_rvec(axis3[i], axis3_t0[i]);
 403             }
 404         }
 405         else
 406         {
 407             fprintf(fptilt, "%15.12g       ", t);
 408             fprintf(fprotation, "%15.12g       ", t);
 409             fprintf(fptheta1, "%15.12g      ", t);
 410             fprintf(fptheta2, "%15.12g      ", t);
 411             fprintf(fptheta3, "%15.12g      ", t);
 412
 413             for (i = 0; i < iCA; i++)
 414             {
 415                 if (i == 0 || i == iCA - 1)
 416                 {
 417                     tilt = rotation = 0;
 418                 }
 419                 else
 420                 {
 421                     if (!bIncremental)
 422                     {
 423                         /* Total rotation & tilt */
 424                         copy_rvec(residuehelixaxis_t0[i], refaxes[0]);
 425                         copy_rvec(residuevector_t0[i], refaxes[1]);
 426                         copy_rvec(axis3_t0[i], refaxes[2]);
 427                     }
 428                     else
 429                     {
 430                         /* Rotation/tilt since last step */
 431                         copy_rvec(residuehelixaxis_tlast[i], refaxes[0]);
 432                         copy_rvec(residuevector_tlast[i], refaxes[1]);
 433                         copy_rvec(axis3_tlast[i], refaxes[2]);
 434                     }
 435                     copy_rvec(residuehelixaxis[i], newaxes[0]);
 436                     copy_rvec(residuevector[i], newaxes[1]);
 437                     copy_rvec(axis3[i], newaxes[2]);
 438
 439                     /* rotate reference frame onto unit axes */
 440                     calc_fit_R(3, 3, weight, unitaxes, refaxes, A);
 441                     for (j = 0; j < 3; j++)
 442                     {
 443                         mvmul(A, refaxes[j], rot_refaxes[j]);
 444                         mvmul(A, newaxes[j], rot_newaxes[j]);
 445                     }
 446
 447                     /* Determine local rotation matrix A */
 448                     calc_fit_R(3, 3, weight, rot_newaxes, rot_refaxes, A);
 449                     /* Calculate euler angles, from rotation order y-z-x, where
 450                      * x is helixaxis, y residuevector, and z axis3.
 451                      *
 452                      * A contains rotation column vectors.
 453                      */
 454
 455                     theta1 = 180.0 / M_PI * std::atan2(A[0][2], A[0][0]);
 456                     theta2 = 180.0 / M_PI * std::asin(-A[0][1]);
 457                     theta3 = 180.0 / M_PI * std::atan2(A[2][1], A[1][1]);
 458
 459                     tilt     = std::sqrt(theta1 * theta1 + theta2 * theta2);
 460                     rotation = theta3;
 461                     fprintf(fptheta1, "%15.12g  ", theta1);
 462                     fprintf(fptheta2, "%15.12g  ", theta2);
 463                     fprintf(fptheta3, "%15.12g  ", theta3);
 464                 }
 465                 fprintf(fptilt, "%15.12g  ", tilt);
 466                 fprintf(fprotation, "%15.12g  ", rotation);
 467             }
 468             fprintf(fptilt, "\n");
 469             fprintf(fprotation, "\n");
 470             fprintf(fptheta1, "\n");
 471             fprintf(fptheta2, "\n");
 472             fprintf(fptheta3, "\n");
 473         }
 474
 475         for (i = 0; i < iCA; i++)
 476         {
 477             copy_rvec(residuehelixaxis[i], residuehelixaxis_tlast[i]);
 478             copy_rvec(residuevector[i], residuevector_tlast[i]);
 479             copy_rvec(axis3[i], axis3_tlast[i]);
 480         }
 481
 482         teller++;
 483     } while (read_next_x(oenv, status, &t, x, box));
 484
 485     gmx_rmpbc_done(gpbc);
 486
 487     gmx_ffclose(fpaxis);
 488     gmx_ffclose(fpcenter);
 489     xvgrclose(fptilt);
 490     xvgrclose(fprotation);
 491     gmx_ffclose(fprise);
 492     gmx_ffclose(fpradius);
 493     gmx_ffclose(fptwist);
 494     gmx_ffclose(fpbending);
 495     gmx_ffclose(fptheta1);
 496     gmx_ffclose(fptheta2);
 497     gmx_ffclose(fptheta3);
 498
 499     close_trx(status);
 500
 501     return 0;
 502 }