src/tools/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  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <math.h>
  43 #include <string.h>
  44
  45 #include "statutil.h"
  46 #include "sysstuff.h"
  47 #include "typedefs.h"
  48 #include "smalloc.h"
  49 #include "macros.h"
  50 #include "gmx_fatal.h"
  51 #include "vec.h"
  52 #include "copyrite.h"
  53 #include "futil.h"
  54 #include "statutil.h"
  55 #include "index.h"
  56 #include "pdbio.h"
  57 #include "tpxio.h"
  58 #include "txtdump.h"
  59 #include "physics.h"
  60 #include "random.h"
  61 #include "eigio.h"
  62 #include "gmx_ana.h"
  63
  64
  65 int gmx_nmtraj(int argc,char *argv[])
  66 {
  67     const char *desc[] =
  68     {
  69         "[TT]g_nmtraj[tt] generates an virtual trajectory from an eigenvector, ",
  70         "corresponding to a harmonic Cartesian oscillation around the average ",
  71         "structure. The eigenvectors should normally be mass-weighted, but you can ",
  72         "use non-weighted eigenvectors to generate orthogonal motions. ",
  73         "The output frames are written as a trajectory file covering an entire period, and ",
  74         "the first frame is the average structure. If you write the trajectory in (or convert to) ",
  75         "PDB format you can view it directly in PyMol and also render a photorealistic movie. ",
  76         "Motion amplitudes are calculated from the eigenvalues and a preset temperature, ",
  77         "assuming equipartition of the energy over all modes. To make the motion clearly visible ",
  78         "in PyMol you might want to amplify it by setting an unrealistically high temperature. ",
  79         "However, be aware that both the linear Cartesian displacements and mass weighting will ",
  80         "lead to serious structure deformation for high amplitudes - this is is simply a limitation ",
  81         "of the Cartesian normal mode model. By default the selected eigenvector is set to 7, since ",
  82         " the first six normal modes are the translational and rotational degrees of freedom."
  83     };
  84
  85     static real refamplitude=0.25;
  86     static int  nframes=30;
  87     static real temp=300.0;
  88     static const char *eignrvec = "7";
  89     static const char *phasevec  = "0.0";
  90
  91     t_pargs pa[] =
  92     {
  93         { "-eignr",     FALSE, etSTR,  {&eignrvec}, "String of eigenvectors to use (first is 1)" },
  94         { "-phases",    FALSE, etSTR,  {&phasevec}, "String of phases (default is 0.0)" },
  95         { "-temp",      FALSE, etREAL, {&temp},      "Temperature (K)" },
  96         { "-amplitude", FALSE, etREAL, {&refamplitude}, "Amplitude for modes with eigenvalue<=0" },
  97         { "-nframes",   FALSE, etINT,  {&nframes},   "Number of frames to generate" }
  98     };
  99
 100 #define NPA asize(pa)
 101
 102   t_trxstatus *out;
 103   t_topology top;
 104   int        ePBC;
 105   t_atoms    *atoms;
 106   rvec       *xtop,*xref,*xav,*xout;
 107   int        nvec,*eignr=NULL;
 108   int        *eigvalnr;
 109   rvec       **eigvec=NULL;
 110   matrix     box;
 111   int        natoms;
 112   int        i,j,k,kmode,d,s,v;
 113   gmx_bool       bDMR,bDMA,bFit;
 114   char *     indexfile;
 115
 116   char *     grpname;
 117   real *     eigval;
 118   int        neigval;
 119   int *      dummy;
 120   real *     invsqrtm;
 121   char       title[STRLEN];
 122   real       fraction;
 123   int        *out_eigidx;
 124   real       *out_eigval;
 125   rvec *     this_eigvec;
 126   real       omega,Ekin,sum,m,vel;
 127   gmx_bool       found;
 128   int        nmodes,nphases;
 129   int        *imodes;
 130   real       *amplitude;
 131   real       *phases;
 132   real       dum;
 133   const char       *p;
 134   char *pe;
 135   output_env_t oenv;
 136
 137   t_filenm fnm[] =
 138   {
 139       { efTPS, NULL,    NULL,          ffREAD },
 140       { efTRN, "-v",    "eigenvec",    ffREAD  },
 141       { efTRO, "-o",    "nmtraj",      ffWRITE }
 142   };
 143
 144 #define NFILE asize(fnm)
 145
 146   CopyRight(stderr,argv[0]);
 147   parse_common_args(&argc,argv,PCA_BE_NICE,
 148                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
 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       dummy[i]=i;
 219
 220   /* Find the eigenvalue/vector to match our select one */
 221   snew(out_eigidx,nmodes);
 222   for(i=0;i<nmodes;i++)
 223     out_eigidx[i]=-1;
 224
 225   for(i=0;i<nvec;i++)
 226   {
 227       for(j=0;j<nmodes;j++)
 228       {
 229           if(imodes[j]==eignr[i])
 230               out_eigidx[j]=i;
 231       }
 232   }
 233     for(i=0;i<nmodes;i++)
 234         if(out_eigidx[i]==-1)
 235             gmx_fatal(FARGS,"Could not find mode %d in eigenvector file.\n",imodes[i]);
 236
 237
 238   snew(invsqrtm,natoms);
 239
 240   if (bDMA)
 241   {
 242       for(i=0; (i<natoms); i++)
 243           invsqrtm[i] = gmx_invsqrt(atoms->atom[i].m);
 244   }
 245   else
 246   {
 247       for(i=0; (i<natoms); i++)
 248           invsqrtm[i]=1.0;
 249   }
 250
 251   snew(xout,natoms);
 252   snew(amplitude,nmodes);
 253
 254         printf("mode phases: %g %g\n",phases[0],phases[1]);
 255
 256   for(i=0;i<nmodes;i++)
 257   {
 258       kmode = out_eigidx[i];
 259       this_eigvec=eigvec[kmode];
 260
 261       if( (kmode >= 6) && (eigval[kmode] > 0))
 262       {
 263           /* Derive amplitude from temperature and eigenvalue if we can */
 264
 265           /* Convert eigenvalue to angular frequency, in units s^(-1) */
 266           omega = sqrt(eigval[kmode]*1.0E21/(AVOGADRO*AMU));
 267           /* Harmonic motion will be x=x0 + A*sin(omega*t)*eigenvec.
 268            * The velocity is thus:
 269            *
 270            * v = A*omega*cos(omega*t)*eigenvec.
 271            *
 272            * And the average kinetic energy the integral of mass*v*v/2 over a
 273            * period:
 274            *
 275            * (1/4)*mass*A*omega*eigenvec
 276            *
 277            * For t =2*pi*n, all energy will be kinetic, and v=A*omega*eigenvec.
 278            * The kinetic energy will be sum(0.5*mass*v*v) if we temporarily set A to 1,
 279                    * and the average over a period half of this.
 280            */
 281
 282           Ekin = 0;
 283           for(k=0;k<natoms;k++)
 284           {
 285               m = atoms->atom[k].m;
 286               for(d=0;d<DIM;d++)
 287               {
 288                   vel   = omega*this_eigvec[k][d];
 289                   Ekin += 0.5*0.5*m*vel*vel;
 290               }
 291           }
 292
 293           /* Convert Ekin from amu*(nm/s)^2 to J, i.e., kg*(m/s)^2
 294            * This will also be proportional to A^2
 295            */
 296           Ekin *= AMU*1E-18;
 297
 298           /* Set the amplitude so the energy is kT/2 */
 299           amplitude[i] = sqrt(0.5*BOLTZMANN*temp/Ekin);
 300           }
 301           else
 302           {
 303                   amplitude[i] = refamplitude;
 304           }
 305   }
 306
 307   out=open_trx(ftp2fn(efTRO,NFILE,fnm),"w");
 308
 309     /* Write a sine oscillation around the average structure,
 310      * modulated by the eigenvector with selected amplitude.
 311      */
 312
 313     for(i=0;i<nframes;i++)
 314     {
 315         fraction = (real)i/(real)nframes;
 316                 for(j=0;j<natoms;j++)
 317                 {
 318                         copy_rvec(xav[j],xout[j]);
 319                 }
 320
 321         for(k=0;k<nmodes;k++)
 322         {
 323             kmode=out_eigidx[k];
 324             this_eigvec=eigvec[kmode];
 325
 326             for(j=0;j<natoms;j++)
 327             {
 328                 for(d=0;d<DIM;d++)
 329                 {
 330                                         xout[j][d] += amplitude[k]*sin(2*M_PI*(fraction+phases[k]/360.0))*this_eigvec[j][d];
 331                 }
 332             }
 333         }
 334         write_trx(out,natoms,dummy,atoms,i,(real)i/(real)nframes,box,xout,NULL,NULL);
 335     }
 336
 337     fprintf(stderr,"\n");
 338     close_trx(out);
 339
 340     return 0;
 341 }
 342
 343