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