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, 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/stringutil.h"
  60
  61 int gmx_nmtraj(int argc, char *argv[])
  62 {
  63     const char *desc[] =
  64     {
  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 limitation ",
  77         "of the Cartesian normal mode model. By default the selected eigenvector is set to 7, since ",
  78         "the first six normal modes are the translational and rotational degrees of freedom."
  79     };
  80
  81     static real        refamplitude = 0.25;
  82     static int         nframes      = 30;
  83     static real        temp         = 300.0;
  84     static const char *eignrvec     = "7";
  85     static const char *phasevec     = "0.0";
  86
  87     t_pargs            pa[] =
  88     {
  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     int               ePBC;
 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     real              fraction;
 114     int              *out_eigidx;
 115     rvec        *     this_eigvec;
 116     real              omega, Ekin, m, vel;
 117     int               nmodes, nphases;
 118     int              *imodes;
 119     real             *amplitude;
 120     real             *phases;
 121     const char       *p;
 122     char             *pe;
 123     gmx_output_env_t *oenv;
 124
 125     t_filenm          fnm[] =
 126     {
 127         { efTPS, nullptr,    nullptr,          ffREAD },
 128         { efTRN, "-v",    "eigenvec",    ffREAD  },
 129         { efTRO, "-o",    "nmtraj",      ffWRITE }
 130     };
 131
 132 #define NFILE asize(fnm)
 133
 134     if (!parse_common_args(&argc, argv, 0,
 135                            NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
 136     {
 137         return 0;
 138     }
 139
 140     read_eigenvectors(opt2fn("-v", NFILE, fnm), &natoms, &bFit,
 141                       &xref, &bDMR, &xav, &bDMA, &nvec, &eignr, &eigvec, &eigval);
 142
 143     read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &xtop, nullptr, box, bDMA);
 144
 145     /* Find vectors and phases */
 146
 147     /* first find number of args in string */
 148     nmodes = gmx::countWords(eignrvec);
 149
 150     snew(imodes, nmodes);
 151     p = eignrvec;
 152     for (i = 0; i < nmodes; i++)
 153     {
 154         /* C indices start on 0 */
 155         imodes[i] = std::strtol(p, &pe, 10)-1;
 156         p         = pe;
 157     }
 158
 159     /* Now read phases */
 160     nphases = gmx::countWords(phasevec);
 161
 162     if (nphases > nmodes)
 163     {
 164         gmx_fatal(FARGS, "More phases than eigenvector indices specified.\n");
 165     }
 166
 167     snew(phases, nmodes);
 168     p = phasevec;
 169
 170     for (i = 0; i < nphases; i++)
 171     {
 172         phases[i] = strtod(p, &pe);
 173         p         = pe;
 174     }
 175
 176     if (nmodes > nphases)
 177     {
 178         printf("Warning: Setting phase of last %d modes to zero...\n", nmodes-nphases);
 179     }
 180
 181     for (i = nphases; i < nmodes; i++)
 182     {
 183         phases[i] = 0;
 184     }
 185
 186     atoms = &top.atoms;
 187
 188     if (atoms->nr != natoms)
 189     {
 190         gmx_fatal(FARGS, "Different number of atoms in topology and eigenvectors.\n");
 191     }
 192
 193     snew(dummy, natoms);
 194     for (i = 0; i < natoms; i++)
 195     {
 196         dummy[i] = i;
 197     }
 198
 199     /* Find the eigenvalue/vector to match our select one */
 200     snew(out_eigidx, nmodes);
 201     for (i = 0; i < nmodes; i++)
 202     {
 203         out_eigidx[i] = -1;
 204     }
 205
 206     for (i = 0; i < nvec; i++)
 207     {
 208         for (j = 0; j < nmodes; j++)
 209         {
 210             if (imodes[j] == eignr[i])
 211             {
 212                 out_eigidx[j] = i;
 213             }
 214         }
 215     }
 216     for (i = 0; i < nmodes; i++)
 217     {
 218         if (out_eigidx[i] == -1)
 219         {
 220             gmx_fatal(FARGS, "Could not find mode %d in eigenvector file.\n", imodes[i]);
 221         }
 222     }
 223
 224
 225     snew(invsqrtm, natoms);
 226
 227     if (bDMA)
 228     {
 229         for (i = 0; (i < natoms); i++)
 230         {
 231             invsqrtm[i] = gmx::invsqrt(atoms->atom[i].m);
 232         }
 233     }
 234     else
 235     {
 236         for (i = 0; (i < natoms); i++)
 237         {
 238             invsqrtm[i] = 1.0;
 239         }
 240     }
 241
 242     snew(xout, natoms);
 243     snew(amplitude, nmodes);
 244
 245     printf("mode phases: %g %g\n", phases[0], phases[1]);
 246
 247     for (i = 0; i < nmodes; i++)
 248     {
 249         kmode       = out_eigidx[i];
 250         this_eigvec = eigvec[kmode];
 251
 252         if ( (kmode >= 6) && (eigval[kmode] > 0))
 253         {
 254             /* Derive amplitude from temperature and eigenvalue if we can */
 255
 256             /* Convert eigenvalue to angular frequency, in units s^(-1) */
 257             omega = std::sqrt(eigval[kmode]*1.0E21/(AVOGADRO*AMU));
 258             /* Harmonic motion will be x=x0 + A*sin(omega*t)*eigenvec.
 259              * The velocity is thus:
 260              *
 261              * v = A*omega*cos(omega*t)*eigenvec.
 262              *
 263              * And the average kinetic energy the integral of mass*v*v/2 over a
 264              * period:
 265              *
 266              * (1/4)*mass*A*omega*eigenvec
 267              *
 268              * For t =2*pi*n, all energy will be kinetic, and v=A*omega*eigenvec.
 269              * The kinetic energy will be sum(0.5*mass*v*v) if we temporarily set A to 1,
 270              * and the average over a period half of this.
 271              */
 272
 273             Ekin = 0;
 274             for (k = 0; k < natoms; k++)
 275             {
 276                 m = atoms->atom[k].m;
 277                 for (d = 0; d < DIM; d++)
 278                 {
 279                     vel   = omega*this_eigvec[k][d];
 280                     Ekin += 0.5*0.5*m*vel*vel;
 281                 }
 282             }
 283
 284             /* Convert Ekin from amu*(nm/s)^2 to J, i.e., kg*(m/s)^2
 285              * This will also be proportional to A^2
 286              */
 287             Ekin *= AMU*1E-18;
 288
 289             /* Set the amplitude so the energy is kT/2 */
 290             amplitude[i] = std::sqrt(0.5*BOLTZMANN*temp/Ekin);
 291         }
 292         else
 293         {
 294             amplitude[i] = refamplitude;
 295         }
 296     }
 297
 298     out = open_trx(ftp2fn(efTRO, NFILE, fnm), "w");
 299
 300     /* Write a sine oscillation around the average structure,
 301      * modulated by the eigenvector with selected amplitude.
 302      */
 303
 304     for (i = 0; i < nframes; i++)
 305     {
 306         fraction = static_cast<real>(i)/nframes;
 307         for (j = 0; j < natoms; j++)
 308         {
 309             copy_rvec(xav[j], xout[j]);
 310         }
 311
 312         for (k = 0; k < nmodes; k++)
 313         {
 314             kmode       = out_eigidx[k];
 315             this_eigvec = eigvec[kmode];
 316
 317             for (j = 0; j < natoms; j++)
 318             {
 319                 for (d = 0; d < DIM; d++)
 320                 {
 321                     xout[j][d] += amplitude[k]*std::sin(2*M_PI*(fraction+phases[k]/360.0))*this_eigvec[j][d];
 322                 }
 323             }
 324         }
 325         write_trx(out, natoms, dummy, atoms, i, static_cast<real>(i)/nframes, box, xout, nullptr, nullptr);
 326     }
 327
 328     fprintf(stderr, "\n");
 329     close_trx(out);
 330
 331     return 0;
 332 }