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