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, 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     int   ePBC;
 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(&argc, argv, PCA_CAN_TIME, NFILE, fnm, NPA, pa, asize(desc), desc, 0,
 166                            nullptr, &oenv))
 167     {
 168         return 0;
 169     }
 170
 171     top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC);
 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), "Incremental local helix tilt", "Time(ps)",
 235                           "Tilt (degrees)", oenv);
 236         fprotation = xvgropen(opt2fn("-orot", NFILE, fnm), "Incremental local helix rotation",
 237                               "Time(ps)", "Rotation (degrees)", oenv);
 238     }
 239     else
 240     {
 241         fptilt = xvgropen(opt2fn("-otilt", NFILE, fnm), "Cumulative local helix tilt", "Time(ps)",
 242                           "Tilt (degrees)", oenv);
 243         fprotation = xvgropen(opt2fn("-orot", NFILE, fnm), "Cumulative local helix rotation",
 244                               "Time(ps)", "Rotation (degrees)", oenv);
 245     }
 246
 247     clear_rvecs(3, unitaxes);
 248     unitaxes[0][0] = 1;
 249     unitaxes[1][1] = 1;
 250     unitaxes[2][2] = 1;
 251
 252     gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
 253
 254     do
 255     {
 256         /* initialisation for correct distance calculations */
 257         set_pbc(&pbc, ePBC, box);
 258         /* make molecules whole again */
 259         gmx_rmpbc(gpbc, natoms, box, x);
 260
 261         /* copy coords to our smaller arrays */
 262         for (i = 0; i < iCA; i++)
 263         {
 264             copy_rvec(x[ind_CA[i]], x_CA[i]);
 265             if (bSC)
 266             {
 267                 copy_rvec(x[ind_SC[i]], x_SC[i]);
 268             }
 269         }
 270
 271         for (i = 0; i < iCA - 3; i++)
 272         {
 273             rvec_sub(x_CA[i + 1], x_CA[i], r12[i]);
 274             rvec_sub(x_CA[i + 2], x_CA[i + 1], r23[i]);
 275             rvec_sub(x_CA[i + 3], x_CA[i + 2], r34[i]);
 276             rvec_sub(r12[i], r23[i], diff13[i]);
 277             rvec_sub(r23[i], r34[i], diff24[i]);
 278             /* calculate helix axis */
 279             cprod(diff13[i], diff24[i], helixaxis[i]);
 280             svmul(1.0 / norm(helixaxis[i]), helixaxis[i], helixaxis[i]);
 281
 282             tmp       = cos_angle(diff13[i], diff24[i]);
 283             twist[i]  = 180.0 / M_PI * std::acos(tmp);
 284             radius[i] = std::sqrt(norm(diff13[i]) * norm(diff24[i])) / (2.0 * (1.0 - tmp));
 285             rise[i]   = std::abs(iprod(r23[i], helixaxis[i]));
 286
 287             svmul(radius[i] / norm(diff13[i]), diff13[i], v1);
 288             svmul(radius[i] / norm(diff24[i]), diff24[i], v2);
 289
 290             rvec_sub(x_CA[i + 1], v1, residueorigin[i + 1]);
 291             rvec_sub(x_CA[i + 2], v2, residueorigin[i + 2]);
 292         }
 293         residueradius[0] = residuetwist[0] = residuerise[0] = 0;
 294
 295         residueradius[1] = radius[0];
 296         residuetwist[1]  = twist[0];
 297         residuerise[1]   = rise[0];
 298
 299         residuebending[0] = residuebending[1] = 0;
 300         for (i = 2; i < iCA - 2; i++)
 301         {
 302             residueradius[i] = 0.5 * (radius[i - 2] + radius[i - 1]);
 303             residuetwist[i]  = 0.5 * (twist[i - 2] + twist[i - 1]);
 304             residuerise[i]   = 0.5 * (rise[i - 2] + rise[i - 1]);
 305             residuebending[i] = 180.0 / M_PI * std::acos(cos_angle(helixaxis[i - 2], helixaxis[i - 1]));
 306         }
 307         residueradius[iCA - 2] = radius[iCA - 4];
 308         residuetwist[iCA - 2]  = twist[iCA - 4];
 309         residuerise[iCA - 2]   = rise[iCA - 4];
 310         residueradius[iCA - 1] = residuetwist[iCA - 1] = residuerise[iCA - 1] = 0;
 311         residuebending[iCA - 2] = residuebending[iCA - 1] = 0;
 312
 313         clear_rvec(residueorigin[0]);
 314         clear_rvec(residueorigin[iCA - 1]);
 315
 316         /* average helix axes to define them on the residues.
 317          * Just extrapolate second first/list atom.
 318          */
 319         copy_rvec(helixaxis[0], residuehelixaxis[0]);
 320         copy_rvec(helixaxis[0], residuehelixaxis[1]);
 321
 322         for (i = 2; i < iCA - 2; i++)
 323         {
 324             rvec_add(helixaxis[i - 2], helixaxis[i - 1], residuehelixaxis[i]);
 325             svmul(0.5, residuehelixaxis[i], residuehelixaxis[i]);
 326         }
 327         copy_rvec(helixaxis[iCA - 4], residuehelixaxis[iCA - 2]);
 328         copy_rvec(helixaxis[iCA - 4], residuehelixaxis[iCA - 1]);
 329
 330         /* Normalize the axis */
 331         for (i = 0; i < iCA; i++)
 332         {
 333             svmul(1.0 / norm(residuehelixaxis[i]), residuehelixaxis[i], residuehelixaxis[i]);
 334         }
 335
 336         /* calculate vector from origin to residue CA */
 337         fprintf(fpaxis, "%15.12g  ", t);
 338         fprintf(fpcenter, "%15.12g  ", t);
 339         fprintf(fprise, "%15.12g  ", t);
 340         fprintf(fpradius, "%15.12g  ", t);
 341         fprintf(fptwist, "%15.12g  ", t);
 342         fprintf(fpbending, "%15.12g  ", t);
 343
 344         for (i = 0; i < iCA; i++)
 345         {
 346             if (i == 0 || i == iCA - 1)
 347             {
 348                 fprintf(fpaxis, "%15.12g %15.12g %15.12g       ", 0.0, 0.0, 0.0);
 349                 fprintf(fpcenter, "%15.12g %15.12g %15.12g       ", 0.0, 0.0, 0.0);
 350                 fprintf(fprise, "%15.12g  ", 0.0);
 351                 fprintf(fpradius, "%15.12g  ", 0.0);
 352                 fprintf(fptwist, "%15.12g  ", 0.0);
 353                 fprintf(fpbending, "%15.12g  ", 0.0);
 354             }
 355             else
 356             {
 357                 rvec_sub(bSC ? x_SC[i] : x_CA[i], residueorigin[i], residuevector[i]);
 358                 svmul(1.0 / norm(residuevector[i]), residuevector[i], residuevector[i]);
 359                 cprod(residuehelixaxis[i], residuevector[i], axis3[i]);
 360                 fprintf(fpaxis, "%15.12g %15.12g %15.12g       ", residuehelixaxis[i][0],
 361                         residuehelixaxis[i][1], residuehelixaxis[i][2]);
 362                 fprintf(fpcenter, "%15.12g %15.12g %15.12g       ", residueorigin[i][0],
 363                         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                 fprintf(fptilt, "%15.12g  ", tilt);
 448                 fprintf(fprotation, "%15.12g  ", rotation);
 449             }
 450             fprintf(fptilt, "\n");
 451             fprintf(fprotation, "\n");
 452             fprintf(fptheta1, "\n");
 453             fprintf(fptheta2, "\n");
 454             fprintf(fptheta3, "\n");
 455         }
 456
 457         for (i = 0; i < iCA; i++)
 458         {
 459             copy_rvec(residuehelixaxis[i], residuehelixaxis_tlast[i]);
 460             copy_rvec(residuevector[i], residuevector_tlast[i]);
 461             copy_rvec(axis3[i], axis3_tlast[i]);
 462         }
 463
 464         teller++;
 465     } while (read_next_x(oenv, status, &t, x, box));
 466
 467     gmx_rmpbc_done(gpbc);
 468
 469     gmx_ffclose(fpaxis);
 470     gmx_ffclose(fpcenter);
 471     xvgrclose(fptilt);
 472     xvgrclose(fprotation);
 473     gmx_ffclose(fprise);
 474     gmx_ffclose(fpradius);
 475     gmx_ffclose(fptwist);
 476     gmx_ffclose(fpbending);
 477     gmx_ffclose(fptheta1);
 478     gmx_ffclose(fptheta2);
 479     gmx_ffclose(fptheta3);
 480
 481     close_trx(status);
 482
 483     return 0;
 484 }