src/gromacs/gmxana/gmx_nmtraj.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 <cctype>
  40 #include <cmath>
  41 #include <cstdlib>
  42 #include <cstring>
  43
  44 #include "gromacs/commandline/pargs.h"
  45 #include "gromacs/fileio/confio.h"
  46 #include "gromacs/fileio/pdbio.h"
  47 #include "gromacs/fileio/trxio.h"
  48 #include "gromacs/gmxana/eigio.h"
  49 #include "gromacs/gmxana/gmx_ana.h"
  50 #include "gromacs/math/functions.h"
  51 #include "gromacs/math/units.h"
  52 #include "gromacs/math/vec.h"
  53 #include "gromacs/topology/topology.h"
  54 #include "gromacs/utility/arraysize.h"
  55 #include "gromacs/utility/cstringutil.h"
  56 #include "gromacs/utility/fatalerror.h"
  57 #include "gromacs/utility/futil.h"
  58 #include "gromacs/utility/smalloc.h"
  59 #include "gromacs/utility/strconvert.h"
  60 #include "gromacs/utility/stringutil.h"
  61
  62 int gmx_nmtraj(int argc, char* argv[])
  63 {
  64     const char* desc[] = {
  65         "[THISMODULE] generates an virtual trajectory from an eigenvector, ",
  66         "corresponding to a harmonic Cartesian oscillation around the average ",
  67         "structure. The eigenvectors should normally be mass-weighted, but you can ",
  68         "use non-weighted eigenvectors to generate orthogonal motions. ",
  69         "The output frames are written as a trajectory file covering an entire period, and ",
  70         "the first frame is the average structure. If you write the trajectory in (or convert to) ",
  71         "PDB format you can view it directly in PyMol and also render a photorealistic movie. ",
  72         "Motion amplitudes are calculated from the eigenvalues and a preset temperature, ",
  73         "assuming equipartition of the energy over all modes. To make the motion clearly visible ",
  74         "in PyMol you might want to amplify it by setting an unrealistically high temperature. ",
  75         "However, be aware that both the linear Cartesian displacements and mass weighting will ",
  76         "lead to serious structure deformation for high amplitudes - this is is simply a ",
  77         "limitation of the Cartesian normal mode model. By default the selected eigenvector ",
  78         "is set to 7, since ",
  79         "the first six normal modes are the translational and rotational degrees of freedom."
  80     };
  81
  82     static real        refamplitude = 0.25;
  83     static int         nframes      = 30;
  84     static real        temp         = 300.0;
  85     static const char* eignrvec     = "7";
  86     static const char* phasevec     = "0.0";
  87
  88     t_pargs pa[] = {
  89         { "-eignr", FALSE, etSTR, { &eignrvec }, "String of eigenvectors to use (first is 1)" },
  90         { "-phases", FALSE, etSTR, { &phasevec }, "String of phases (default is 0.0)" },
  91         { "-temp", FALSE, etREAL, { &temp }, "Temperature (K)" },
  92         { "-amplitude", FALSE, etREAL, { &refamplitude }, "Amplitude for modes with eigenvalue<=0" },
  93         { "-nframes", FALSE, etINT, { &nframes }, "Number of frames to generate" }
  94     };
  95
  96 #define NPA asize(pa)
  97
  98     t_trxstatus* out;
  99     t_topology   top;
 100     PbcType      pbcType;
 101     t_atoms*     atoms;
 102     rvec *       xtop, *xref, *xav, *xout;
 103     int          nvec, *eignr = nullptr;
 104     rvec**       eigvec = nullptr;
 105     matrix       box;
 106     int          natoms;
 107     int          i, j, k, kmode, d;
 108     gmx_bool     bDMR, bDMA, bFit;
 109
 110     real*             eigval;
 111     int*              dummy;
 112     real*             invsqrtm;
 113     int*              out_eigidx;
 114     rvec*             this_eigvec;
 115     real              omega, Ekin, m, vel;
 116     real*             amplitude;
 117     gmx_output_env_t* oenv;
 118
 119     t_filenm fnm[] = { { efTPS, nullptr, nullptr, ffREAD },
 120                        { efTRN, "-v", "eigenvec", ffREAD },
 121                        { efTRO, "-o", "nmtraj", ffWRITE } };
 122
 123 #define NFILE asize(fnm)
 124
 125     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
 126     {
 127         return 0;
 128     }
 129
 130     read_eigenvectors(opt2fn("-v", NFILE, fnm), &natoms, &bFit, &xref, &bDMR, &xav, &bDMA, &nvec,
 131                       &eignr, &eigvec, &eigval);
 132
 133     read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, &xtop, nullptr, box, bDMA);
 134
 135     /* Find vectors and phases */
 136
 137     std::vector<int> imodes;
 138     for (const auto& imodeString : gmx::splitString(eignrvec))
 139     {
 140         imodes.emplace_back(gmx::fromStdString<int>(imodeString));
 141     }
 142     int nmodes = gmx::ssize(imodes);
 143
 144     std::vector<real> phases;
 145     phases.reserve(nmodes);
 146     for (const auto& phaseString : gmx::splitString(phasevec))
 147     {
 148         phases.emplace_back(gmx::fromStdString<int>(phaseString));
 149     }
 150     int nphases = gmx::ssize(phases);
 151
 152     if (nphases > nmodes)
 153     {
 154         gmx_fatal(FARGS, "More phases than eigenvector indices specified.\n");
 155     }
 156
 157     if (nmodes > nphases)
 158     {
 159         printf("Warning: Setting phase of last %d modes to zero...\n", nmodes - nphases);
 160         phases.resize(nmodes, 0);
 161     }
 162
 163     atoms = &top.atoms;
 164
 165     if (atoms->nr != natoms)
 166     {
 167         gmx_fatal(FARGS, "Different number of atoms in topology and eigenvectors.\n");
 168     }
 169
 170     snew(dummy, natoms);
 171     for (i = 0; i < natoms; i++)
 172     {
 173         dummy[i] = i;
 174     }
 175
 176     /* Find the eigenvalue/vector to match our select one */
 177     snew(out_eigidx, nmodes);
 178     for (i = 0; i < nmodes; i++)
 179     {
 180         out_eigidx[i] = -1;
 181     }
 182
 183     for (i = 0; i < nvec; i++)
 184     {
 185         for (j = 0; j < nmodes; j++)
 186         {
 187             if (imodes[j] == eignr[i])
 188             {
 189                 out_eigidx[j] = i;
 190             }
 191         }
 192     }
 193     for (i = 0; i < nmodes; i++)
 194     {
 195         if (out_eigidx[i] == -1)
 196         {
 197             gmx_fatal(FARGS, "Could not find mode %d in eigenvector file.\n", imodes[i]);
 198         }
 199     }
 200
 201
 202     snew(invsqrtm, natoms);
 203
 204     if (bDMA)
 205     {
 206         for (i = 0; (i < natoms); i++)
 207         {
 208             invsqrtm[i] = gmx::invsqrt(atoms->atom[i].m);
 209         }
 210     }
 211     else
 212     {
 213         for (i = 0; (i < natoms); i++)
 214         {
 215             invsqrtm[i] = 1.0;
 216         }
 217     }
 218
 219     snew(xout, natoms);
 220     snew(amplitude, nmodes);
 221
 222     printf("mode phases: %g %g\n", phases[0], phases[1]);
 223
 224     for (i = 0; i < nmodes; i++)
 225     {
 226         kmode       = out_eigidx[i];
 227         this_eigvec = eigvec[kmode];
 228
 229         if ((kmode >= 6) && (eigval[kmode] > 0))
 230         {
 231             /* Derive amplitude from temperature and eigenvalue if we can */
 232
 233             /* Convert eigenvalue to angular frequency, in units s^(-1) */
 234             omega = std::sqrt(eigval[kmode] * 1.0E21 / (AVOGADRO * AMU));
 235             /* Harmonic motion will be x=x0 + A*sin(omega*t)*eigenvec.
 236              * The velocity is thus:
 237              *
 238              * v = A*omega*cos(omega*t)*eigenvec.
 239              *
 240              * And the average kinetic energy the integral of mass*v*v/2 over a
 241              * period:
 242              *
 243              * (1/4)*mass*A*omega*eigenvec
 244              *
 245              * For t =2*pi*n, all energy will be kinetic, and v=A*omega*eigenvec.
 246              * The kinetic energy will be sum(0.5*mass*v*v) if we temporarily set A to 1,
 247              * and the average over a period half of this.
 248              */
 249
 250             Ekin = 0;
 251             for (k = 0; k < natoms; k++)
 252             {
 253                 m = atoms->atom[k].m;
 254                 for (d = 0; d < DIM; d++)
 255                 {
 256                     vel = omega * this_eigvec[k][d];
 257                     Ekin += 0.5 * 0.5 * m * vel * vel;
 258                 }
 259             }
 260
 261             /* Convert Ekin from amu*(nm/s)^2 to J, i.e., kg*(m/s)^2
 262              * This will also be proportional to A^2
 263              */
 264             Ekin *= AMU * 1E-18;
 265
 266             /* Set the amplitude so the energy is kT/2 */
 267             amplitude[i] = std::sqrt(0.5 * BOLTZMANN * temp / Ekin);
 268         }
 269         else
 270         {
 271             amplitude[i] = refamplitude;
 272         }
 273     }
 274
 275     out = open_trx(ftp2fn(efTRO, NFILE, fnm), "w");
 276
 277     /* Write a sine oscillation around the average structure,
 278      * modulated by the eigenvector with selected amplitude.
 279      */
 280
 281     for (i = 0; i < nframes; i++)
 282     {
 283         real fraction = static_cast<real>(i) / static_cast<real>(nframes);
 284         for (j = 0; j < natoms; j++)
 285         {
 286             copy_rvec(xav[j], xout[j]);
 287         }
 288
 289         for (k = 0; k < nmodes; k++)
 290         {
 291             kmode       = out_eigidx[k];
 292             this_eigvec = eigvec[kmode];
 293
 294             for (j = 0; j < natoms; j++)
 295             {
 296                 for (d = 0; d < DIM; d++)
 297                 {
 298                     xout[j][d] += amplitude[k] * std::sin(2 * M_PI * (fraction + phases[k] / 360.0))
 299                                   * this_eigvec[j][d];
 300                 }
 301             }
 302         }
 303         write_trx(out, natoms, dummy, atoms, i, fraction, box, xout, nullptr, nullptr);
 304     }
 305
 306     fprintf(stderr, "\n");
 307     close_trx(out);
 308
 309     return 0;
 310 }