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     "[TT]g_nmeig[tt] 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