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