src/tools/gmx_nmeig.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 "vec.h"
  48 #include "pbc.h"
  49 #include "copyrite.h"
  50 #include "futil.h"
  51 #include "statutil.h"
  52 #include "index.h"
  53 #include "mshift.h"
  54 #include "xvgr.h"
  55 #include "gstat.h"
  56 #include "txtdump.h"
  57 #include "eigensolver.h"
  58 #include "eigio.h"
  59 #include "mtxio.h"
  60 #include "sparsematrix.h"
  61 #include "physics.h"
  62 #include "main.h"
  63 #include "gmx_ana.h"
  64
  65
  66 static void
  67 nma_full_hessian(real *           hess,
  68                  int              ndim,
  69                  gmx_bool             bM,
  70                  t_topology *     top,
  71                  int              begin,
  72                  int              end,
  73                  real *           eigenvalues,
  74                  real *           eigenvectors)
  75 {
  76     int i,j,k,l;
  77     real mass_fac,rdum;
  78     int natoms;
  79
  80     natoms = top->atoms.nr;
  81
  82     /* divide elements hess[i][j] by sqrt(mas[i])*sqrt(mas[j]) when required */
  83
  84     if (bM)
  85     {
  86         for (i=0; (i<natoms); i++)
  87         {
  88             for (j=0; (j<DIM); j++)
  89             {
  90                 for (k=0; (k<natoms); k++)
  91                 {
  92                     mass_fac=gmx_invsqrt(top->atoms.atom[i].m*top->atoms.atom[k].m);
  93                     for (l=0; (l<DIM); l++)
  94                         hess[(i*DIM+j)*ndim+k*DIM+l]*=mass_fac;
  95                 }
  96             }
  97         }
  98     }
  99
 100     /* call diagonalization routine. */
 101
 102     fprintf(stderr,"\nDiagonalizing to find vectors %d through %d...\n",begin,end);
 103     fflush(stderr);
 104
 105     eigensolver(hess,ndim,begin-1,end-1,eigenvalues,eigenvectors);
 106
 107     /* And scale the output eigenvectors */
 108     if (bM && eigenvectors!=NULL)
 109     {
 110         for(i=0;i<(end-begin+1);i++)
 111         {
 112             for(j=0;j<natoms;j++)
 113             {
 114                 mass_fac = gmx_invsqrt(top->atoms.atom[j].m);
 115                 for (k=0; (k<DIM); k++)
 116                 {
 117                     eigenvectors[i*ndim+j*DIM+k] *= mass_fac;
 118                 }
 119             }
 120         }
 121     }
 122 }
 123
 124
 125
 126 static void
 127 nma_sparse_hessian(gmx_sparsematrix_t *     sparse_hessian,
 128                    gmx_bool                     bM,
 129                    t_topology *             top,
 130                    int                      neig,
 131                    real *                   eigenvalues,
 132                    real *                   eigenvectors)
 133 {
 134     int i,j,k;
 135     int row,col;
 136     real mass_fac;
 137     int iatom,katom;
 138     int natoms;
 139     int ndim;
 140
 141     natoms = top->atoms.nr;
 142     ndim   = DIM*natoms;
 143
 144     /* Cannot check symmetry since we only store half matrix */
 145     /* divide elements hess[i][j] by sqrt(mas[i])*sqrt(mas[j]) when required */
 146
 147     if (bM)
 148     {
 149         for (iatom=0; (iatom<natoms); iatom++)
 150         {
 151             for (j=0; (j<DIM); j++)
 152             {
 153                 row = DIM*iatom+j;
 154                 for(k=0;k<sparse_hessian->ndata[row];k++)
 155                 {
 156                     col = sparse_hessian->data[row][k].col;
 157                     katom = col/3;
 158                     mass_fac=gmx_invsqrt(top->atoms.atom[iatom].m*top->atoms.atom[katom].m);
 159                     sparse_hessian->data[row][k].value *=mass_fac;
 160                 }
 161             }
 162         }
 163     }
 164     fprintf(stderr,"\nDiagonalizing to find eigenvectors 1 through %d...\n",neig);
 165     fflush(stderr);
 166
 167     sparse_eigensolver(sparse_hessian,neig,eigenvalues,eigenvectors,10000000);
 168
 169     /* Scale output eigenvectors */
 170     if (bM && eigenvectors!=NULL)
 171     {
 172         for(i=0;i<neig;i++)
 173         {
 174             for(j=0;j<natoms;j++)
 175             {
 176                 mass_fac = gmx_invsqrt(top->atoms.atom[j].m);
 177                 for (k=0; (k<DIM); k++)
 178                 {
 179                     eigenvectors[i*ndim+j*DIM+k] *= mass_fac;
 180                 }
 181             }
 182         }
 183     }
 184 }
 185
 186
 187
 188 int gmx_nmeig(int argc,char *argv[])
 189 {
 190   const char *desc[] = {
 191     "g_nmeig calculates the eigenvectors/values of a (Hessian) matrix,",
 192     "which can be calculated with [TT]mdrun[tt].",
 193     "The eigenvectors are written to a trajectory file ([TT]-v[tt]).",
 194     "The structure is written first with t=0. The eigenvectors",
 195     "are written as frames with the eigenvector number as timestamp.",
 196     "The eigenvectors can be analyzed with [TT]g_anaeig[tt].",
 197     "An ensemble of structures can be generated from the eigenvectors with",
 198     "[TT]g_nmens[tt]. When mass weighting is used, the generated eigenvectors",
 199     "will be scaled back to plain cartesian coordinates before generating the",
 200     "output - in this case they will no longer be exactly orthogonal in the",
 201     "standard cartesian norm (But in the mass weighted norm they would be)."
 202   };
 203
 204   static gmx_bool bM=TRUE;
 205   static int  begin=1,end=50;
 206   t_pargs pa[] =
 207   {
 208     { "-m",  FALSE, etBOOL, {&bM},
 209       "Divide elements of Hessian by product of sqrt(mass) of involved "
 210       "atoms prior to diagonalization. This should be used for 'Normal Modes' "
 211       "analysis" },
 212     { "-first", FALSE, etINT, {&begin},
 213       "First eigenvector to write away" },
 214     { "-last",  FALSE, etINT, {&end},
 215       "Last eigenvector to write away" }
 216   };
 217   FILE       *out;
 218   int        status,trjout;
 219   t_topology top;
 220   int        ePBC;
 221   rvec       *top_x;
 222   matrix     box;
 223   real       *eigenvalues;
 224   real       *eigenvectors;
 225   real       rdum,mass_fac;
 226   int        natoms,ndim,nrow,ncol,count;
 227   char       *grpname,title[256];
 228   int        i,j,k,l,d,gnx;
 229   gmx_bool       bSuck;
 230   atom_id    *index;
 231   real       value;
 232   real       factor_gmx_to_omega2;
 233   real       factor_omega_to_wavenumber;
 234   t_commrec  *cr;
 235   output_env_t oenv;
 236
 237   real *                 full_hessian   = NULL;
 238   gmx_sparsematrix_t *   sparse_hessian = NULL;
 239
 240   t_filenm fnm[] = {
 241     { efMTX, "-f", "hessian",    ffREAD  },
 242     { efTPS, NULL, NULL,         ffREAD  },
 243     { efXVG, "-of", "eigenfreq", ffWRITE },
 244     { efXVG, "-ol", "eigenval",  ffWRITE },
 245     { efTRN, "-v", "eigenvec",  ffWRITE }
 246   };
 247 #define NFILE asize(fnm)
 248
 249         cr = init_par(&argc,&argv);
 250
 251         if(MASTER(cr))
 252                 CopyRight(stderr,argv[0]);
 253
 254   parse_common_args(&argc,argv,PCA_BE_NICE | (MASTER(cr) ? 0 : PCA_QUIET),
 255                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 256
 257   read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&top_x,NULL,box,bM);
 258
 259   natoms = top.atoms.nr;
 260   ndim = DIM*natoms;
 261
 262   if(begin<1)
 263       begin = 1;
 264   if(end>ndim)
 265       end = ndim;
 266
 267   /*open Hessian matrix */
 268   gmx_mtxio_read(ftp2fn(efMTX,NFILE,fnm),&nrow,&ncol,&full_hessian,&sparse_hessian);
 269
 270   /* Memory for eigenvalues and eigenvectors (begin..end) */
 271   snew(eigenvalues,nrow);
 272   snew(eigenvectors,nrow*(end-begin+1));
 273
 274   /* If the Hessian is in sparse format we can calculate max (ndim-1) eigenvectors,
 275    * and they must start at the first one. If this is not valid we convert to full matrix
 276    * storage, but warn the user that we might run out of memory...
 277    */
 278   if((sparse_hessian != NULL) && (begin!=1 || end==ndim))
 279   {
 280       if(begin!=1)
 281       {
 282           fprintf(stderr,"Cannot use sparse Hessian with first eigenvector != 1.\n");
 283       }
 284       else if(end==ndim)
 285       {
 286           fprintf(stderr,"Cannot use sparse Hessian to calculate all eigenvectors.\n");
 287       }
 288
 289       fprintf(stderr,"Will try to allocate memory and convert to full matrix representation...\n");
 290
 291       snew(full_hessian,nrow*ncol);
 292       for(i=0;i<nrow*ncol;i++)
 293           full_hessian[i] = 0;
 294
 295       for(i=0;i<sparse_hessian->nrow;i++)
 296       {
 297           for(j=0;j<sparse_hessian->ndata[i];j++)
 298           {
 299               k     = sparse_hessian->data[i][j].col;
 300               value = sparse_hessian->data[i][j].value;
 301               full_hessian[i*ndim+k] = value;
 302               full_hessian[k*ndim+i] = value;
 303           }
 304       }
 305       gmx_sparsematrix_destroy(sparse_hessian);
 306       sparse_hessian = NULL;
 307       fprintf(stderr,"Converted sparse to full matrix storage.\n");
 308   }
 309
 310   if(full_hessian != NULL)
 311   {
 312       /* Using full matrix storage */
 313       nma_full_hessian(full_hessian,nrow,bM,&top,begin,end,eigenvalues,eigenvectors);
 314   }
 315   else
 316   {
 317       /* Sparse memory storage, allocate memory for eigenvectors */
 318       snew(eigenvectors,ncol*end);
 319       nma_sparse_hessian(sparse_hessian,bM,&top,end,eigenvalues,eigenvectors);
 320   }
 321
 322
 323   /* check the output, first 6 eigenvalues should be reasonably small */
 324   bSuck=FALSE;
 325   for (i=begin-1; (i<6); i++)
 326   {
 327       if (fabs(eigenvalues[i]) > 1.0e-3)
 328           bSuck=TRUE;
 329   }
 330   if (bSuck)
 331   {
 332       fprintf(stderr,"\nOne of the lowest 6 eigenvalues has a non-zero value.\n");
 333       fprintf(stderr,"This could mean that the reference structure was not\n");
 334       fprintf(stderr,"properly energy minimized.\n");
 335   }
 336
 337
 338   /* now write the output */
 339   fprintf (stderr,"Writing eigenvalues...\n");
 340   out=xvgropen(opt2fn("-ol",NFILE,fnm),
 341                "Eigenvalues","Eigenvalue index","Eigenvalue [Gromacs units]",
 342                oenv);
 343   if (output_env_get_print_xvgr_codes(oenv)) {
 344     if (bM)
 345       fprintf(out,"@ subtitle \"mass weighted\"\n");
 346     else
 347       fprintf(out,"@ subtitle \"not mass weighted\"\n");
 348   }
 349
 350   for (i=0; i<=(end-begin); i++)
 351       fprintf (out,"%6d %15g\n",begin+i,eigenvalues[i]);
 352   ffclose(out);
 353
 354
 355
 356   fprintf(stderr,"Writing eigenfrequencies - negative eigenvalues will be set to zero.\n");
 357
 358   out=xvgropen(opt2fn("-of",NFILE,fnm),
 359                "Eigenfrequencies","Eigenvector index","Wavenumber [cm\\S-1\\N]",
 360                oenv);
 361   if (output_env_get_print_xvgr_codes(oenv)) {
 362     if (bM)
 363       fprintf(out,"@ subtitle \"mass weighted\"\n");
 364     else
 365       fprintf(out,"@ subtitle \"not mass weighted\"\n");
 366   }
 367
 368   /* Gromacs units are kJ/(mol*nm*nm*amu),
 369    * where amu is the atomic mass unit.
 370    *
 371    * For the eigenfrequencies we want to convert this to spectroscopic absorption
 372    * wavenumbers given in cm^(-1), which is the frequency divided by the speed of
 373    * light. Do this by first converting to omega^2 (units 1/s), take the square
 374    * root, and finally divide by the speed of light (nm/ps in gromacs).
 375    */
 376   factor_gmx_to_omega2       = 1.0E21/(AVOGADRO*AMU);
 377   factor_omega_to_wavenumber = 1.0E-5/(2.0*M_PI*SPEED_OF_LIGHT);
 378
 379   for (i=0; i<=(end-begin); i++)
 380   {
 381       value = eigenvalues[i];
 382       if(value < 0)
 383           value = 0;
 384       value=sqrt(value*factor_gmx_to_omega2)*factor_omega_to_wavenumber;
 385       fprintf (out,"%6d %15g\n",begin+i,value);
 386   }
 387   ffclose(out);
 388
 389   /* Writing eigenvectors. Note that if mass scaling was used, the eigenvectors
 390    * were scaled back from mass weighted cartesian to plain cartesian in the
 391    * nma_full_hessian() or nma_sparse_hessian() routines. Mass scaled vectors
 392    * will not be strictly orthogonal in plain cartesian scalar products.
 393    */
 394   write_eigenvectors(opt2fn("-v",NFILE,fnm),natoms,eigenvectors,FALSE,begin,end,
 395                      eWXR_NO,NULL,FALSE,top_x,bM,eigenvalues);
 396
 397   thanx(stderr);
 398
 399   return 0;
 400 }
 401