src/tools/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 "copyrite.h"
  50 #include "futil.h"
  51 #include "statutil.h"
  52 #include "index.h"
  53 #include "pdbio.h"
  54 #include "tpxio.h"
  55 #include "txtdump.h"
  56 #include "physics.h"
  57 #include "random.h"
  58 #include "eigio.h"
  59 #include "gmx_ana.h"
  60
  61
  62 int gmx_nmtraj(int argc,char *argv[])
  63 {
  64     const char *desc[] =
  65     {
  66         "[TT]g_nmtraj[tt] generates an virtual trajectory from an eigenvector, ",
  67         "corresponding to a harmonic cartesian oscillation around the average ",
  68         "structure. The eigenvectors should normally be mass-weighted, but you can ",
  69         "use non-weighted eigenvectors to generate orthogonal motions. ",
  70         "The output frames are written as a trajectory file covering an entire period, and ",
  71         "the first frame is the average structure. If you write the trajectory in (or convert to) ",
  72         "PDB format you can view it directly in PyMol and also render a photorealistic movie. ",
  73         "Motion amplitudes are calculated from the eigenvalues and a preset temperature, ",
  74         "assuming equipartition of the energy over all modes. To make the motion clearly visible ",
  75         "in PyMol you might want to amplify it by setting an unrealistic high temperature. ",
  76         "However, be aware that both the linear cartesian displacements and mass weighting will ",
  77         "lead to serious structure deformation for high amplitudes - this is is simply a limitation ",
  78         "of the cartesian normal mode model. By default the selected eigenvector is set to 7, since ",
  79         " the first six normal modes are the translational and rotational degrees of freedom."
  80     };
  81
  82     static real refamplitude=0.25;
  83     static int  nframes=30;
  84     static real temp=300.0;
  85     static const char *eignrvec = "7";
  86     static const char *phasevec  = "0.0";
  87
  88     t_pargs pa[] =
  89     {
  90         { "-eignr",     FALSE, etSTR,  {&eignrvec}, "String of eigenvectors to use (first is 1)" },
  91         { "-phases",    FALSE, etSTR,  {&phasevec}, "String of phases (default is 0.0)" },
  92         { "-temp",      FALSE, etREAL, {&temp},      "Temperature in Kelvin" },
  93         { "-amplitude", FALSE, etREAL, {&refamplitude}, "Amplitude for modes with eigenvalue<=0" },
  94         { "-nframes",   FALSE, etINT,  {&nframes},   "Number of frames to generate" }
  95     };
  96
  97 #define NPA asize(pa)
  98
  99   int        out;
 100   t_topology top;
 101   int        ePBC;
 102   t_atoms    *atoms;
 103   rvec       *xtop,*xref,*xav,*xout;
 104   int        nvec,*eignr=NULL;
 105   int        *eigvalnr;
 106   rvec       **eigvec=NULL;
 107   matrix     box;
 108   int        natoms;
 109   int        i,j,k,kmode,d,s,v;
 110   bool       bDMR,bDMA,bFit;
 111   char *     indexfile;
 112
 113   char *     grpname;
 114   real *     eigval;
 115   int        neigval;
 116   int *      dummy;
 117   real *     invsqrtm;
 118   char       title[STRLEN];
 119   real       fraction;
 120   int        *out_eigidx;
 121   real       *out_eigval;
 122   rvec *     this_eigvec;
 123   real       omega,Ekin,sum,m,vel;
 124   bool       found;
 125   int        nmodes,nphases;
 126   int        *imodes;
 127   real       *amplitude;
 128   real       *phases;
 129   real       dum;
 130   const char       *p;
 131   char *pe;
 132   output_env_t oenv;
 133
 134   t_filenm fnm[] =
 135   {
 136       { efTPS, NULL,    NULL,          ffREAD },
 137       { efTRN, "-v",    "eigenvec",    ffREAD  },
 138       { efTRO, "-o",    "nmtraj",      ffWRITE }
 139   };
 140
 141 #define NFILE asize(fnm)
 142
 143   CopyRight(stderr,argv[0]);
 144   parse_common_args(&argc,argv,PCA_BE_NICE,
 145                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
 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       dummy[i]=i;
 216
 217   /* Find the eigenvalue/vector to match our select one */
 218   snew(out_eigidx,nmodes);
 219   for(i=0;i<nmodes;i++)
 220     out_eigidx[i]=-1;
 221
 222   for(i=0;i<nvec;i++)
 223   {
 224       for(j=0;j<nmodes;j++)
 225       {
 226           if(imodes[j]==eignr[i])
 227               out_eigidx[j]=i;
 228       }
 229   }
 230     for(i=0;i<nmodes;i++)
 231         if(out_eigidx[i]==-1)
 232             gmx_fatal(FARGS,"Could not find mode %d in eigenvector file.\n",imodes[i]);
 233
 234
 235   snew(invsqrtm,natoms);
 236
 237   if (bDMA)
 238   {
 239       for(i=0; (i<natoms); i++)
 240           invsqrtm[i] = gmx_invsqrt(atoms->atom[i].m);
 241   }
 242   else
 243   {
 244       for(i=0; (i<natoms); i++)
 245           invsqrtm[i]=1.0;
 246   }
 247
 248   snew(xout,natoms);
 249   snew(amplitude,nmodes);
 250
 251         printf("mode phases: %g %g\n",phases[0],phases[1]);
 252
 253   for(i=0;i<nmodes;i++)
 254   {
 255       kmode = out_eigidx[i];
 256       this_eigvec=eigvec[kmode];
 257
 258       if( (kmode >= 6) && (eigval[kmode] > 0))
 259       {
 260           /* Derive amplitude from temperature and eigenvalue if we can */
 261
 262           /* Convert eigenvalue to angular frequency, in units s^(-1) */
 263           omega = sqrt(eigval[kmode]*1.0E21/(AVOGADRO*AMU));
 264           /* Harmonic motion will be x=x0 + A*sin(omega*t)*eigenvec.
 265            * The velocity is thus:
 266            *
 267            * v = A*omega*cos(omega*t)*eigenvec.
 268            *
 269            * And the average kinetic energy the integral of mass*v*v/2 over a
 270            * period:
 271            *
 272            * (1/4)*mass*A*omega*eigenvec
 273            *
 274            * For t =2*pi*n, all energy will be kinetic, and v=A*omega*eigenvec.
 275            * The kinetic energy will be sum(0.5*mass*v*v) if we temporarily set A to 1,
 276                    * and the average over a period half of this.
 277            */
 278
 279           Ekin = 0;
 280           for(k=0;k<natoms;k++)
 281           {
 282               m = atoms->atom[k].m;
 283               for(d=0;d<DIM;d++)
 284               {
 285                   vel   = omega*this_eigvec[k][d];
 286                   Ekin += 0.5*0.5*m*vel*vel;
 287               }
 288           }
 289
 290           /* Convert Ekin from amu*(nm/s)^2 to J, i.e., kg*(m/s)^2
 291            * This will also be proportional to A^2
 292            */
 293           Ekin *= AMU*1E-18;
 294
 295           /* Set the amplitude so the energy is kT/2 */
 296           amplitude[i] = sqrt(0.5*BOLTZMANN*temp/Ekin);
 297           }
 298           else
 299           {
 300                   amplitude[i] = refamplitude;
 301           }
 302   }
 303
 304   out=open_trx(ftp2fn(efTRO,NFILE,fnm),"w");
 305
 306     /* Write a sine oscillation around the average structure,
 307      * modulated by the eigenvector with selected amplitude.
 308      */
 309
 310     for(i=0;i<nframes;i++)
 311     {
 312         fraction = (real)i/(real)nframes;
 313                 for(j=0;j<natoms;j++)
 314                 {
 315                         copy_rvec(xav[j],xout[j]);
 316                 }
 317
 318         for(k=0;k<nmodes;k++)
 319         {
 320             kmode=out_eigidx[k];
 321             this_eigvec=eigvec[kmode];
 322
 323             for(j=0;j<natoms;j++)
 324             {
 325                 for(d=0;d<DIM;d++)
 326                 {
 327                                         xout[j][d] += amplitude[k]*sin(2*M_PI*(fraction+phases[k]/360.0))*this_eigvec[j][d];
 328                 }
 329             }
 330         }
 331         write_trx(out,natoms,dummy,atoms,i,(real)i/(real)nframes,box,xout,NULL,NULL);
 332     }
 333
 334     fprintf(stderr,"\n");
 335     close_trx(out);
 336
 337     return 0;
 338 }
 339
 340