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 "mtop_util.h"
  61 #include "sparsematrix.h"
  62 #include "physics.h"
  63 #include "main.h"
  64 #include "gmx_ana.h"
  65
  66 static double cv_corr(double nu,double T)
  67 {
  68     double x = PLANCK*nu/(BOLTZ*T);
  69     double ex = exp(x);
  70
  71     if (nu <= 0)
  72         return BOLTZ*KILO;
  73     else
  74         return BOLTZ*KILO*(ex*sqr(x)/sqr(ex-1) - 1);
  75 }
  76
  77 static double u_corr(double nu,double T)
  78 {
  79     double x = PLANCK*nu/(BOLTZ*T);
  80     double ex = exp(x);
  81
  82     if (nu <= 0)
  83         return BOLTZ*T;
  84     else
  85         return BOLTZ*T*(0.5*x - 1 + x/(ex-1));
  86 }
  87
  88 static int get_nharm_mt(gmx_moltype_t *mt)
  89 {
  90     static int harm_func[] = { F_BONDS };
  91     int    i,ft,nh;
  92
  93     nh = 0;
  94     for(i=0; (i<asize(harm_func)); i++)
  95     {
  96         ft = harm_func[i];
  97         nh += mt->ilist[ft].nr/(interaction_function[ft].nratoms+1);
  98     }
  99     return nh;
 100 }
 101
 102 static int get_nvsite_mt(gmx_moltype_t *mt)
 103 {
 104     static int vs_func[] = { F_VSITE2, F_VSITE3, F_VSITE3FD, F_VSITE3FAD,
 105                              F_VSITE3OUT, F_VSITE4FD, F_VSITE4FDN, F_VSITEN };
 106     int    i,ft,nh;
 107
 108     nh = 0;
 109     for(i=0; (i<asize(vs_func)); i++)
 110     {
 111         ft = vs_func[i];
 112         nh += mt->ilist[ft].nr/(interaction_function[ft].nratoms+1);
 113     }
 114     return nh;
 115 }
 116
 117 static int get_nharm(gmx_mtop_t *mtop,int *nvsites)
 118 {
 119     int j,mt,nh,nv;
 120
 121     nh = 0;
 122     nv = 0;
 123     for(j=0; (j<mtop->nmolblock); j++)
 124     {
 125         mt = mtop->molblock[j].type;
 126         nh += mtop->molblock[j].nmol * get_nharm_mt(&(mtop->moltype[mt]));
 127         nv += mtop->molblock[j].nmol * get_nvsite_mt(&(mtop->moltype[mt]));
 128     }
 129     *nvsites = nv;
 130     return nh;
 131 }
 132
 133 static void
 134 nma_full_hessian(real *           hess,
 135                  int              ndim,
 136                  gmx_bool             bM,
 137                  t_topology *     top,
 138                  int              begin,
 139                  int              end,
 140                  real *           eigenvalues,
 141                  real *           eigenvectors)
 142 {
 143     int i,j,k,l;
 144     real mass_fac,rdum;
 145     int natoms;
 146
 147     natoms = top->atoms.nr;
 148
 149     /* divide elements hess[i][j] by sqrt(mas[i])*sqrt(mas[j]) when required */
 150
 151     if (bM)
 152     {
 153         for (i=0; (i<natoms); i++)
 154         {
 155             for (j=0; (j<DIM); j++)
 156             {
 157                 for (k=0; (k<natoms); k++)
 158                 {
 159                     mass_fac=gmx_invsqrt(top->atoms.atom[i].m*top->atoms.atom[k].m);
 160                     for (l=0; (l<DIM); l++)
 161                         hess[(i*DIM+j)*ndim+k*DIM+l]*=mass_fac;
 162                 }
 163             }
 164         }
 165     }
 166
 167     /* call diagonalization routine. */
 168
 169     fprintf(stderr,"\nDiagonalizing to find vectors %d through %d...\n",begin,end);
 170     fflush(stderr);
 171
 172     eigensolver(hess,ndim,begin-1,end-1,eigenvalues,eigenvectors);
 173
 174     /* And scale the output eigenvectors */
 175     if (bM && eigenvectors!=NULL)
 176     {
 177         for(i=0;i<(end-begin+1);i++)
 178         {
 179             for(j=0;j<natoms;j++)
 180             {
 181                 mass_fac = gmx_invsqrt(top->atoms.atom[j].m);
 182                 for (k=0; (k<DIM); k++)
 183                 {
 184                     eigenvectors[i*ndim+j*DIM+k] *= mass_fac;
 185                 }
 186             }
 187         }
 188     }
 189 }
 190
 191
 192
 193 static void
 194 nma_sparse_hessian(gmx_sparsematrix_t *     sparse_hessian,
 195                    gmx_bool                     bM,
 196                    t_topology *             top,
 197                    int                      neig,
 198                    real *                   eigenvalues,
 199                    real *                   eigenvectors)
 200 {
 201     int i,j,k;
 202     int row,col;
 203     real mass_fac;
 204     int iatom,katom;
 205     int natoms;
 206     int ndim;
 207
 208     natoms = top->atoms.nr;
 209     ndim   = DIM*natoms;
 210
 211     /* Cannot check symmetry since we only store half matrix */
 212     /* divide elements hess[i][j] by sqrt(mas[i])*sqrt(mas[j]) when required */
 213
 214     if (bM)
 215     {
 216         for (iatom=0; (iatom<natoms); iatom++)
 217         {
 218             for (j=0; (j<DIM); j++)
 219             {
 220                 row = DIM*iatom+j;
 221                 for(k=0;k<sparse_hessian->ndata[row];k++)
 222                 {
 223                     col = sparse_hessian->data[row][k].col;
 224                     katom = col/3;
 225                     mass_fac=gmx_invsqrt(top->atoms.atom[iatom].m*top->atoms.atom[katom].m);
 226                     sparse_hessian->data[row][k].value *=mass_fac;
 227                 }
 228             }
 229         }
 230     }
 231     fprintf(stderr,"\nDiagonalizing to find eigenvectors 1 through %d...\n",neig);
 232     fflush(stderr);
 233
 234     sparse_eigensolver(sparse_hessian,neig,eigenvalues,eigenvectors,10000000);
 235
 236     /* Scale output eigenvectors */
 237     if (bM && eigenvectors!=NULL)
 238     {
 239         for(i=0;i<neig;i++)
 240         {
 241             for(j=0;j<natoms;j++)
 242             {
 243                 mass_fac = gmx_invsqrt(top->atoms.atom[j].m);
 244                 for (k=0; (k<DIM); k++)
 245                 {
 246                     eigenvectors[i*ndim+j*DIM+k] *= mass_fac;
 247                 }
 248             }
 249         }
 250     }
 251 }
 252
 253
 254
 255 int gmx_nmeig(int argc,char *argv[])
 256 {
 257   const char *desc[] = {
 258     "[TT]g_nmeig[tt] calculates the eigenvectors/values of a (Hessian) matrix,",
 259     "which can be calculated with [TT]mdrun[tt].",
 260     "The eigenvectors are written to a trajectory file ([TT]-v[tt]).",
 261     "The structure is written first with t=0. The eigenvectors",
 262     "are written as frames with the eigenvector number as timestamp.",
 263     "The eigenvectors can be analyzed with [TT]g_anaeig[tt].",
 264     "An ensemble of structures can be generated from the eigenvectors with",
 265     "[TT]g_nmens[tt]. When mass weighting is used, the generated eigenvectors",
 266     "will be scaled back to plain Cartesian coordinates before generating the",
 267     "output. In this case, they will no longer be exactly orthogonal in the",
 268     "standard Cartesian norm, but in the mass-weighted norm they would be.[PAR]",
 269     "This program can be optionally used to compute quantum corrections to heat capacity",
 270     "and enthalpy by providing an extra file argument [TT]-qcorr[tt]. See the GROMACS",
 271     "manual, Chapter 1, for details. The result includes subtracting a harmonic",
 272     "degree of freedom at the given temperature.",
 273     "The total correction is printed on the terminal screen.",
 274     "The recommended way of getting the corrections out is:[PAR]",
 275     "[TT]g_nmeig -s topol.tpr -f nm.mtx -first 7 -last 10000 -T 300 -qc [-constr][tt][PAR]",
 276     "The [TT]-constr[tt] option should be used when bond constraints were used during the",
 277     "simulation [BB]for all the covalent bonds[bb]. If this is not the case, ",
 278     "you need to analyze the [TT]quant_corr.xvg[tt] file yourself.[PAR]",
 279     "To make things more flexible, the program can also take virtual sites into account",
 280     "when computing quantum corrections. When selecting [TT]-constr[tt] and",
 281     "[TT]-qc[tt], the [TT]-begin[tt] and [TT]-end[tt] options will be set automatically as well.",
 282     "Again, if you think you know it better, please check the [TT]eigenfreq.xvg[tt]",
 283     "output."
 284   };
 285
 286   static gmx_bool bM=TRUE,bCons=FALSE;
 287   static int  begin=1,end=50,maxspec=4000;
 288   static real T=298.15,width=1;
 289   t_pargs pa[] =
 290   {
 291     { "-m",  FALSE, etBOOL, {&bM},
 292       "Divide elements of Hessian by product of sqrt(mass) of involved "
 293       "atoms prior to diagonalization. This should be used for 'Normal Modes' "
 294       "analysis" },
 295     { "-first", FALSE, etINT, {&begin},
 296       "First eigenvector to write away" },
 297     { "-last",  FALSE, etINT, {&end},
 298       "Last eigenvector to write away" },
 299     { "-maxspec", FALSE, etINT, {&maxspec},
 300       "Highest frequency (1/cm) to consider in the spectrum" },
 301     { "-T",     FALSE, etREAL, {&T},
 302       "Temperature for computing quantum heat capacity and enthalpy when using normal mode calculations to correct classical simulations" },
 303     { "-constr", FALSE, etBOOL, {&bCons},
 304       "If constraints were used in the simulation but not in the normal mode analysis (this is the recommended way of doing it) you will need to set this for computing the quantum corrections." },
 305     { "-width",  FALSE, etREAL, {&width},
 306       "Width (sigma) of the gaussian peaks (1/cm) when generating a spectrum" }
 307   };
 308   FILE       *out,*qc,*spec;
 309   int        status,trjout;
 310   t_topology top;
 311   gmx_mtop_t mtop;
 312   int        ePBC;
 313   rvec       *top_x;
 314   matrix     box;
 315   real       *eigenvalues;
 316   real       *eigenvectors;
 317   real       rdum,mass_fac,qcvtot,qutot,qcv,qu;
 318   int        natoms,ndim,nrow,ncol,count,nharm,nvsite;
 319   char       *grpname;
 320   int        i,j,k,l,d,gnx;
 321   gmx_bool   bSuck;
 322   atom_id    *index;
 323   t_tpxheader tpx;
 324   int        version,generation;
 325   real       value,omega,nu;
 326   real       factor_gmx_to_omega2;
 327   real       factor_omega_to_wavenumber;
 328   real       *spectrum=NULL;
 329   real       wfac;
 330   t_commrec  *cr;
 331   output_env_t oenv;
 332   const char *qcleg[] = { "Heat Capacity cV (J/mol K)",
 333                           "Enthalpy H (kJ/mol)" };
 334   real *                 full_hessian   = NULL;
 335   gmx_sparsematrix_t *   sparse_hessian = NULL;
 336
 337   t_filenm fnm[] = {
 338     { efMTX, "-f", "hessian",    ffREAD  },
 339     { efTPX, NULL, NULL,         ffREAD  },
 340     { efXVG, "-of", "eigenfreq", ffWRITE },
 341     { efXVG, "-ol", "eigenval",  ffWRITE },
 342     { efXVG, "-os", "spectrum",  ffOPTWR },
 343     { efXVG, "-qc", "quant_corr",  ffOPTWR },
 344     { efTRN, "-v", "eigenvec",  ffWRITE }
 345   };
 346 #define NFILE asize(fnm)
 347
 348   cr = init_par(&argc,&argv);
 349
 350   if (MASTER(cr))
 351     CopyRight(stderr,argv[0]);
 352
 353   parse_common_args(&argc,argv,PCA_BE_NICE | (MASTER(cr) ? 0 : PCA_QUIET),
 354                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 355
 356   /* Read tpr file for volume and number of harmonic terms */
 357   read_tpxheader(ftp2fn(efTPX,NFILE,fnm),&tpx,TRUE,&version,&generation);
 358   snew(top_x,tpx.natoms);
 359
 360   read_tpx(ftp2fn(efTPX,NFILE,fnm),NULL,box,&natoms,
 361            top_x,NULL,NULL,&mtop);
 362   if (bCons)
 363   {
 364       nharm = get_nharm(&mtop,&nvsite);
 365   }
 366   else
 367   {
 368       nharm = 0;
 369       nvsite = 0;
 370   }
 371   top = gmx_mtop_t_to_t_topology(&mtop);
 372
 373   bM = TRUE;
 374   ndim = DIM*natoms;
 375
 376   if (opt2bSet("-qc",NFILE,fnm))
 377   {
 378       begin = 7+DIM*nvsite;
 379       end = DIM*natoms;
 380   }
 381   if (begin < 1)
 382       begin = 1;
 383   if (end > ndim)
 384       end = ndim;
 385   printf("Using begin = %d and end = %d\n",begin,end);
 386
 387   /*open Hessian matrix */
 388   gmx_mtxio_read(ftp2fn(efMTX,NFILE,fnm),&nrow,&ncol,&full_hessian,&sparse_hessian);
 389
 390   /* Memory for eigenvalues and eigenvectors (begin..end) */
 391   snew(eigenvalues,nrow);
 392   snew(eigenvectors,nrow*(end-begin+1));
 393
 394   /* If the Hessian is in sparse format we can calculate max (ndim-1) eigenvectors,
 395    * and they must start at the first one. If this is not valid we convert to full matrix
 396    * storage, but warn the user that we might run out of memory...
 397    */
 398   if((sparse_hessian != NULL) && (begin!=1 || end==ndim))
 399   {
 400       if(begin!=1)
 401       {
 402           fprintf(stderr,"Cannot use sparse Hessian with first eigenvector != 1.\n");
 403       }
 404       else if(end==ndim)
 405       {
 406           fprintf(stderr,"Cannot use sparse Hessian to calculate all eigenvectors.\n");
 407       }
 408
 409       fprintf(stderr,"Will try to allocate memory and convert to full matrix representation...\n");
 410
 411       snew(full_hessian,nrow*ncol);
 412       for(i=0;i<nrow*ncol;i++)
 413           full_hessian[i] = 0;
 414
 415       for(i=0;i<sparse_hessian->nrow;i++)
 416       {
 417           for(j=0;j<sparse_hessian->ndata[i];j++)
 418           {
 419               k     = sparse_hessian->data[i][j].col;
 420               value = sparse_hessian->data[i][j].value;
 421               full_hessian[i*ndim+k] = value;
 422               full_hessian[k*ndim+i] = value;
 423           }
 424       }
 425       gmx_sparsematrix_destroy(sparse_hessian);
 426       sparse_hessian = NULL;
 427       fprintf(stderr,"Converted sparse to full matrix storage.\n");
 428   }
 429
 430   if (full_hessian != NULL)
 431   {
 432       /* Using full matrix storage */
 433       nma_full_hessian(full_hessian,nrow,bM,&top,begin,end,
 434                        eigenvalues,eigenvectors);
 435   }
 436   else
 437   {
 438       /* Sparse memory storage, allocate memory for eigenvectors */
 439       snew(eigenvectors,ncol*end);
 440       nma_sparse_hessian(sparse_hessian,bM,&top,end,eigenvalues,eigenvectors);
 441   }
 442
 443   /* check the output, first 6 eigenvalues should be reasonably small */
 444   bSuck=FALSE;
 445   for (i=begin-1; (i<6); i++)
 446   {
 447       if (fabs(eigenvalues[i]) > 1.0e-3)
 448           bSuck=TRUE;
 449   }
 450   if (bSuck)
 451   {
 452       fprintf(stderr,"\nOne of the lowest 6 eigenvalues has a non-zero value.\n");
 453       fprintf(stderr,"This could mean that the reference structure was not\n");
 454       fprintf(stderr,"properly energy minimized.\n");
 455   }
 456
 457   /* now write the output */
 458   fprintf (stderr,"Writing eigenvalues...\n");
 459   out=xvgropen(opt2fn("-ol",NFILE,fnm),
 460                "Eigenvalues","Eigenvalue index","Eigenvalue [Gromacs units]",
 461                oenv);
 462   if (output_env_get_print_xvgr_codes(oenv)) {
 463     if (bM)
 464       fprintf(out,"@ subtitle \"mass weighted\"\n");
 465     else
 466       fprintf(out,"@ subtitle \"not mass weighted\"\n");
 467   }
 468
 469   for (i=0; i<=(end-begin); i++)
 470       fprintf (out,"%6d %15g\n",begin+i,eigenvalues[i]);
 471   ffclose(out);
 472
 473
 474   if (opt2bSet("-qc",NFILE,fnm)) {
 475     qc = xvgropen(opt2fn("-qc",NFILE,fnm),"Quantum Corrections","Eigenvector index","",oenv);
 476     xvgr_legend(qc,asize(qcleg),qcleg,oenv);
 477     qcvtot = qutot = 0;
 478   }
 479   else
 480     qc = NULL;
 481   printf("Writing eigenfrequencies - negative eigenvalues will be set to zero.\n");
 482
 483   out=xvgropen(opt2fn("-of",NFILE,fnm),
 484                "Eigenfrequencies","Eigenvector index","Wavenumber [cm\\S-1\\N]",
 485                oenv);
 486   if (output_env_get_print_xvgr_codes(oenv)) {
 487     if (bM)
 488       fprintf(out,"@ subtitle \"mass weighted\"\n");
 489     else
 490       fprintf(out,"@ subtitle \"not mass weighted\"\n");
 491   }
 492   /* Spectrum ? */
 493   spec = NULL;
 494   if (opt2bSet("-os",NFILE,fnm) && (maxspec > 0))
 495   {
 496       snew(spectrum,maxspec);
 497       spec=xvgropen(opt2fn("-os",NFILE,fnm),
 498                     "Vibrational spectrum based on harmonic approximation",
 499                     "\\f{12}w\\f{4} (cm\\S-1\\N)",
 500                     "Intensity [Gromacs units]",
 501                     oenv);
 502       for(i=0; (i<maxspec); i++)
 503       {
 504           spectrum[i] = 0;
 505       }
 506   }
 507
 508   /* Gromacs units are kJ/(mol*nm*nm*amu),
 509    * where amu is the atomic mass unit.
 510    *
 511    * For the eigenfrequencies we want to convert this to spectroscopic absorption
 512    * wavenumbers given in cm^(-1), which is the frequency divided by the speed of
 513    * light. Do this by first converting to omega^2 (units 1/s), take the square
 514    * root, and finally divide by the speed of light (nm/ps in gromacs).
 515    */
 516   factor_gmx_to_omega2       = 1.0E21/(AVOGADRO*AMU);
 517   factor_omega_to_wavenumber = 1.0E-5/(2.0*M_PI*SPEED_OF_LIGHT);
 518
 519   for (i=begin; (i<=end); i++)
 520   {
 521       value = eigenvalues[i-begin];
 522       if (value < 0)
 523           value = 0;
 524       omega = sqrt(value*factor_gmx_to_omega2);
 525       nu    = 1e-12*omega/(2*M_PI);
 526       value = omega*factor_omega_to_wavenumber;
 527       fprintf (out,"%6d %15g\n",i,value);
 528       if (NULL != spec)
 529       {
 530           wfac = eigenvalues[i-begin]/(width*sqrt(2*M_PI));
 531           for(j=0; (j<maxspec); j++)
 532           {
 533               spectrum[j] += wfac*exp(-sqr(j-value)/(2*sqr(width)));
 534           }
 535       }
 536       if (NULL != qc) {
 537           qcv = cv_corr(nu,T);
 538           qu  = u_corr(nu,T);
 539           if (i > end-nharm)
 540           {
 541               qcv += BOLTZ*KILO;
 542               qu  += BOLTZ*T;
 543           }
 544           fprintf (qc,"%6d %15g %15g\n",i,qcv,qu);
 545           qcvtot += qcv;
 546           qutot += qu;
 547       }
 548   }
 549   ffclose(out);
 550   if (NULL != spec)
 551   {
 552       for(j=0; (j<maxspec); j++)
 553       {
 554           fprintf(spec,"%10g  %10g\n",1.0*j,spectrum[j]);
 555       }
 556       ffclose(spec);
 557   }
 558   if (NULL != qc) {
 559     printf("Quantum corrections for harmonic degrees of freedom\n");
 560     printf("Use appropriate -first and -last options to get reliable results.\n");
 561     printf("There were %d constraints and %d vsites in the simulation\n",
 562            nharm,nvsite);
 563     printf("Total correction to cV = %g J/mol K\n",qcvtot);
 564     printf("Total correction to  H = %g kJ/mol\n",qutot);
 565     ffclose(qc);
 566     please_cite(stdout,"Caleman2011b");
 567   }
 568   /* Writing eigenvectors. Note that if mass scaling was used, the eigenvectors
 569    * were scaled back from mass weighted cartesian to plain cartesian in the
 570    * nma_full_hessian() or nma_sparse_hessian() routines. Mass scaled vectors
 571    * will not be strictly orthogonal in plain cartesian scalar products.
 572    */
 573   write_eigenvectors(opt2fn("-v",NFILE,fnm),natoms,eigenvectors,FALSE,begin,end,
 574                      eWXR_NO,NULL,FALSE,top_x,bM,eigenvalues);
 575
 576   thanx(stderr);
 577
 578   return 0;
 579 }
 580