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