src/gromacs/gmxana/gmx_nmtraj.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <math.h>
  40 #include <string.h>
  41
  42 #include "statutil.h"
  43 #include "sysstuff.h"
  44 #include "typedefs.h"
  45 #include "smalloc.h"
  46 #include "macros.h"
  47 #include "gmx_fatal.h"
  48 #include "vec.h"
  49 #include "futil.h"
  50 #include "statutil.h"
  51 #include "index.h"
  52 #include "pdbio.h"
  53 #include "tpxio.h"
  54 #include "txtdump.h"
  55 #include "physics.h"
  56 #include "random.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         "[TT]g_nmtraj[tt] 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     parse_common_args(&argc, argv, PCA_BE_NICE,
 143                       NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv);
 144
 145     read_eigenvectors(opt2fn("-v", NFILE, fnm), &natoms, &bFit,
 146                       &xref, &bDMR, &xav, &bDMA, &nvec, &eignr, &eigvec, &eigval);
 147
 148     read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, NULL, box, bDMA);
 149
 150     /* Find vectors and phases */
 151
 152     /* first find number of args in string */
 153     nmodes = 0;
 154     p      = eignrvec;
 155     while (*p != 0)
 156     {
 157         dum = strtod(p, &pe);
 158         p   = pe;
 159         nmodes++;
 160     }
 161
 162     snew(imodes, nmodes);
 163     p = eignrvec;
 164     for (i = 0; i < nmodes; i++)
 165     {
 166         /* C indices start on 0 */
 167         imodes[i] = strtol(p, &pe, 10)-1;
 168         p         = pe;
 169     }
 170
 171     /* Now read phases */
 172     nphases = 0;
 173     p       = phasevec;
 174     while (*p != 0)
 175     {
 176         dum = strtod(p, &pe);
 177         p   = pe;
 178         nphases++;
 179     }
 180     if (nphases > nmodes)
 181     {
 182         gmx_fatal(FARGS, "More phases than eigenvector indices specified.\n");
 183     }
 184
 185     snew(phases, nmodes);
 186     p = phasevec;
 187
 188     for (i = 0; i < nphases; i++)
 189     {
 190         phases[i] = strtod(p, &pe);
 191         p         = pe;
 192     }
 193
 194     if (nmodes > nphases)
 195     {
 196         printf("Warning: Setting phase of last %d modes to zero...\n", nmodes-nphases);
 197     }
 198
 199     for (i = nphases; i < nmodes; i++)
 200     {
 201         phases[i] = 0;
 202     }
 203
 204     atoms = &top.atoms;
 205
 206     if (atoms->nr != natoms)
 207     {
 208         gmx_fatal(FARGS, "Different number of atoms in topology and eigenvectors.\n");
 209     }
 210
 211     snew(dummy, natoms);
 212     for (i = 0; i < natoms; i++)
 213     {
 214         dummy[i] = i;
 215     }
 216
 217     /* Find the eigenvalue/vector to match our select one */
 218     snew(out_eigidx, nmodes);
 219     for (i = 0; i < nmodes; i++)
 220     {
 221         out_eigidx[i] = -1;
 222     }
 223
 224     for (i = 0; i < nvec; i++)
 225     {
 226         for (j = 0; j < nmodes; j++)
 227         {
 228             if (imodes[j] == eignr[i])
 229             {
 230                 out_eigidx[j] = i;
 231             }
 232         }
 233     }
 234     for (i = 0; i < nmodes; i++)
 235     {
 236         if (out_eigidx[i] == -1)
 237         {
 238             gmx_fatal(FARGS, "Could not find mode %d in eigenvector file.\n", imodes[i]);
 239         }
 240     }
 241
 242
 243     snew(invsqrtm, natoms);
 244
 245     if (bDMA)
 246     {
 247         for (i = 0; (i < natoms); i++)
 248         {
 249             invsqrtm[i] = gmx_invsqrt(atoms->atom[i].m);
 250         }
 251     }
 252     else
 253     {
 254         for (i = 0; (i < natoms); i++)
 255         {
 256             invsqrtm[i] = 1.0;
 257         }
 258     }
 259
 260     snew(xout, natoms);
 261     snew(amplitude, nmodes);
 262
 263     printf("mode phases: %g %g\n", phases[0], phases[1]);
 264
 265     for (i = 0; i < nmodes; i++)
 266     {
 267         kmode       = out_eigidx[i];
 268         this_eigvec = eigvec[kmode];
 269
 270         if ( (kmode >= 6) && (eigval[kmode] > 0))
 271         {
 272             /* Derive amplitude from temperature and eigenvalue if we can */
 273
 274             /* Convert eigenvalue to angular frequency, in units s^(-1) */
 275             omega = sqrt(eigval[kmode]*1.0E21/(AVOGADRO*AMU));
 276             /* Harmonic motion will be x=x0 + A*sin(omega*t)*eigenvec.
 277              * The velocity is thus:
 278              *
 279              * v = A*omega*cos(omega*t)*eigenvec.
 280              *
 281              * And the average kinetic energy the integral of mass*v*v/2 over a
 282              * period:
 283              *
 284              * (1/4)*mass*A*omega*eigenvec
 285              *
 286              * For t =2*pi*n, all energy will be kinetic, and v=A*omega*eigenvec.
 287              * The kinetic energy will be sum(0.5*mass*v*v) if we temporarily set A to 1,
 288              * and the average over a period half of this.
 289              */
 290
 291             Ekin = 0;
 292             for (k = 0; k < natoms; k++)
 293             {
 294                 m = atoms->atom[k].m;
 295                 for (d = 0; d < DIM; d++)
 296                 {
 297                     vel   = omega*this_eigvec[k][d];
 298                     Ekin += 0.5*0.5*m*vel*vel;
 299                 }
 300             }
 301
 302             /* Convert Ekin from amu*(nm/s)^2 to J, i.e., kg*(m/s)^2
 303              * This will also be proportional to A^2
 304              */
 305             Ekin *= AMU*1E-18;
 306
 307             /* Set the amplitude so the energy is kT/2 */
 308             amplitude[i] = sqrt(0.5*BOLTZMANN*temp/Ekin);
 309         }
 310         else
 311         {
 312             amplitude[i] = refamplitude;
 313         }
 314     }
 315
 316     out = open_trx(ftp2fn(efTRO, NFILE, fnm), "w");
 317
 318     /* Write a sine oscillation around the average structure,
 319      * modulated by the eigenvector with selected amplitude.
 320      */
 321
 322     for (i = 0; i < nframes; i++)
 323     {
 324         fraction = (real)i/(real)nframes;
 325         for (j = 0; j < natoms; j++)
 326         {
 327             copy_rvec(xav[j], xout[j]);
 328         }
 329
 330         for (k = 0; k < nmodes; k++)
 331         {
 332             kmode       = out_eigidx[k];
 333             this_eigvec = eigvec[kmode];
 334
 335             for (j = 0; j < natoms; j++)
 336             {
 337                 for (d = 0; d < DIM; d++)
 338                 {
 339                     xout[j][d] += amplitude[k]*sin(2*M_PI*(fraction+phases[k]/360.0))*this_eigvec[j][d];
 340                 }
 341             }
 342         }
 343         write_trx(out, natoms, dummy, atoms, i, (real)i/(real)nframes, box, xout, NULL, NULL);
 344     }
 345
 346     fprintf(stderr, "\n");
 347     close_trx(out);
 348
 349     return 0;
 350 }