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