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