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