src/gromacs/gmxana/gmx_nmtraj.c

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