src/tools/gmx_nmeig.c

   1 /*  -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   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 -qcorr. See 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:",
 275     "g_nmeig -s topol.tpr -f nm.mtx -first 7 -last 10000 -T 300 -qc [-constr]",
 276     "The constr 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 analyse the quant_corr.xvg file yourself.[PAR]",
 279     "To make things more flexible, the program can also take vsites 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 eigenfreq.xvg",
 283     "output."
 284   };
 285
 286   static gmx_bool bM=TRUE,bCons=FALSE;
 287   static int  begin=1,end=50;
 288   static real T=298.15;
 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     { "-T",     FALSE, etREAL, {&T},
 300       "Temperature for computing quantum heat capacity and enthalpy when using normal mode calculations to correct classical simulations" },
 301     { "-constr", FALSE, etBOOL, {&bCons},
 302       "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." },
 303   };
 304   FILE       *out,*qc;
 305   int        status,trjout;
 306   t_topology top;
 307   gmx_mtop_t mtop;
 308   int        ePBC;
 309   rvec       *top_x;
 310   matrix     box;
 311   real       *eigenvalues;
 312   real       *eigenvectors;
 313   real       rdum,mass_fac,qcvtot,qutot,qcv,qu;
 314   int        natoms,ndim,nrow,ncol,count,nharm,nvsite;
 315   char       *grpname;
 316   int        i,j,k,l,d,gnx;
 317   gmx_bool   bSuck;
 318   atom_id    *index;
 319   t_tpxheader tpx;
 320   int        version,generation;
 321   real       value,omega,nu;
 322   real       factor_gmx_to_omega2;
 323   real       factor_omega_to_wavenumber;
 324   t_commrec  *cr;
 325   output_env_t oenv;
 326   const char *qcleg[] = { "Heat Capacity cV (J/mol K)",
 327                           "Enthalpy H (kJ/mol)" };
 328   real *                 full_hessian   = NULL;
 329   gmx_sparsematrix_t *   sparse_hessian = NULL;
 330
 331   t_filenm fnm[] = {
 332     { efMTX, "-f", "hessian",    ffREAD  },
 333     { efTPX, NULL, NULL,         ffREAD  },
 334     { efXVG, "-of", "eigenfreq", ffWRITE },
 335     { efXVG, "-ol", "eigenval",  ffWRITE },
 336     { efXVG, "-qc", "quant_corr",  ffOPTWR },
 337     { efTRN, "-v", "eigenvec",  ffWRITE }
 338   };
 339 #define NFILE asize(fnm)
 340
 341   cr = init_par(&argc,&argv);
 342
 343   if (MASTER(cr))
 344     CopyRight(stderr,argv[0]);
 345
 346   parse_common_args(&argc,argv,PCA_BE_NICE | (MASTER(cr) ? 0 : PCA_QUIET),
 347                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 348
 349   /* Read tpr file for volume and number of harmonic terms */
 350   read_tpxheader(ftp2fn(efTPX,NFILE,fnm),&tpx,TRUE,&version,&generation);
 351   snew(top_x,tpx.natoms);
 352
 353   read_tpx(ftp2fn(efTPX,NFILE,fnm),NULL,box,&natoms,
 354            top_x,NULL,NULL,&mtop);
 355   if (bCons)
 356   {
 357       nharm = get_nharm(&mtop,&nvsite);
 358   }
 359   else
 360   {
 361       nharm = 0;
 362       nvsite = 0;
 363   }
 364   top = gmx_mtop_t_to_t_topology(&mtop);
 365
 366   bM = TRUE;
 367   ndim = DIM*natoms;
 368
 369   if (opt2bSet("-qc",NFILE,fnm))
 370   {
 371       begin = 7+DIM*nvsite;
 372       end = DIM*natoms;
 373   }
 374   if (begin < 1)
 375       begin = 1;
 376   if (end > ndim)
 377       end = ndim;
 378   printf("Using begin = %d and end = %d\n",begin,end);
 379
 380   /*open Hessian matrix */
 381   gmx_mtxio_read(ftp2fn(efMTX,NFILE,fnm),&nrow,&ncol,&full_hessian,&sparse_hessian);
 382
 383   /* Memory for eigenvalues and eigenvectors (begin..end) */
 384   snew(eigenvalues,nrow);
 385   snew(eigenvectors,nrow*(end-begin+1));
 386
 387   /* If the Hessian is in sparse format we can calculate max (ndim-1) eigenvectors,
 388    * and they must start at the first one. If this is not valid we convert to full matrix
 389    * storage, but warn the user that we might run out of memory...
 390    */
 391   if((sparse_hessian != NULL) && (begin!=1 || end==ndim))
 392   {
 393       if(begin!=1)
 394       {
 395           fprintf(stderr,"Cannot use sparse Hessian with first eigenvector != 1.\n");
 396       }
 397       else if(end==ndim)
 398       {
 399           fprintf(stderr,"Cannot use sparse Hessian to calculate all eigenvectors.\n");
 400       }
 401
 402       fprintf(stderr,"Will try to allocate memory and convert to full matrix representation...\n");
 403
 404       snew(full_hessian,nrow*ncol);
 405       for(i=0;i<nrow*ncol;i++)
 406           full_hessian[i] = 0;
 407
 408       for(i=0;i<sparse_hessian->nrow;i++)
 409       {
 410           for(j=0;j<sparse_hessian->ndata[i];j++)
 411           {
 412               k     = sparse_hessian->data[i][j].col;
 413               value = sparse_hessian->data[i][j].value;
 414               full_hessian[i*ndim+k] = value;
 415               full_hessian[k*ndim+i] = value;
 416           }
 417       }
 418       gmx_sparsematrix_destroy(sparse_hessian);
 419       sparse_hessian = NULL;
 420       fprintf(stderr,"Converted sparse to full matrix storage.\n");
 421   }
 422
 423   if (full_hessian != NULL)
 424   {
 425       /* Using full matrix storage */
 426       nma_full_hessian(full_hessian,nrow,bM,&top,begin,end,
 427                        eigenvalues,eigenvectors);
 428   }
 429   else
 430   {
 431       /* Sparse memory storage, allocate memory for eigenvectors */
 432       snew(eigenvectors,ncol*end);
 433       nma_sparse_hessian(sparse_hessian,bM,&top,end,eigenvalues,eigenvectors);
 434   }
 435
 436   /* check the output, first 6 eigenvalues should be reasonably small */
 437   bSuck=FALSE;
 438   for (i=begin-1; (i<6); i++)
 439   {
 440       if (fabs(eigenvalues[i]) > 1.0e-3)
 441           bSuck=TRUE;
 442   }
 443   if (bSuck)
 444   {
 445       fprintf(stderr,"\nOne of the lowest 6 eigenvalues has a non-zero value.\n");
 446       fprintf(stderr,"This could mean that the reference structure was not\n");
 447       fprintf(stderr,"properly energy minimized.\n");
 448   }
 449
 450   /* now write the output */
 451   fprintf (stderr,"Writing eigenvalues...\n");
 452   out=xvgropen(opt2fn("-ol",NFILE,fnm),
 453                "Eigenvalues","Eigenvalue index","Eigenvalue [Gromacs units]",
 454                oenv);
 455   if (output_env_get_print_xvgr_codes(oenv)) {
 456     if (bM)
 457       fprintf(out,"@ subtitle \"mass weighted\"\n");
 458     else
 459       fprintf(out,"@ subtitle \"not mass weighted\"\n");
 460   }
 461
 462   for (i=0; i<=(end-begin); i++)
 463       fprintf (out,"%6d %15g\n",begin+i,eigenvalues[i]);
 464   ffclose(out);
 465
 466
 467   if (opt2bSet("-qc",NFILE,fnm)) {
 468     qc = xvgropen(opt2fn("-qc",NFILE,fnm),"Quantum Corrections","Eigenvector index","",oenv);
 469     xvgr_legend(qc,asize(qcleg),qcleg,oenv);
 470     qcvtot = qutot = 0;
 471   }
 472   else
 473     qc = NULL;
 474   printf("Writing eigenfrequencies - negative eigenvalues will be set to zero.\n");
 475
 476   out=xvgropen(opt2fn("-of",NFILE,fnm),
 477                "Eigenfrequencies","Eigenvector index","Wavenumber [cm\\S-1\\N]",
 478                oenv);
 479   if (output_env_get_print_xvgr_codes(oenv)) {
 480     if (bM)
 481       fprintf(out,"@ subtitle \"mass weighted\"\n");
 482     else
 483       fprintf(out,"@ subtitle \"not mass weighted\"\n");
 484   }
 485
 486   /* Gromacs units are kJ/(mol*nm*nm*amu),
 487    * where amu is the atomic mass unit.
 488    *
 489    * For the eigenfrequencies we want to convert this to spectroscopic absorption
 490    * wavenumbers given in cm^(-1), which is the frequency divided by the speed of
 491    * light. Do this by first converting to omega^2 (units 1/s), take the square
 492    * root, and finally divide by the speed of light (nm/ps in gromacs).
 493    */
 494   factor_gmx_to_omega2       = 1.0E21/(AVOGADRO*AMU);
 495   factor_omega_to_wavenumber = 1.0E-5/(2.0*M_PI*SPEED_OF_LIGHT);
 496
 497   for (i=begin; (i<=end); i++)
 498   {
 499       value = eigenvalues[i-begin];
 500       if (value < 0)
 501           value = 0;
 502       omega = sqrt(value*factor_gmx_to_omega2);
 503       nu    = 1e-12*omega/(2*M_PI);
 504       value = omega*factor_omega_to_wavenumber;
 505       fprintf (out,"%6d %15g\n",i,value);
 506       if (NULL != qc) {
 507           qcv = cv_corr(nu,T);
 508           qu  = u_corr(nu,T);
 509           if (i > end-nharm)
 510           {
 511               qcv += BOLTZ*KILO;
 512               qu  += BOLTZ*T;
 513           }
 514           fprintf (qc,"%6d %15g %15g\n",i,qcv,qu);
 515           qcvtot += qcv;
 516           qutot += qu;
 517       }
 518   }
 519   ffclose(out);
 520   if (NULL != qc) {
 521     printf("Quantum corrections for harmonic degrees of freedom\n");
 522     printf("Use appropriate -first and -last options to get reliable results.\n");
 523     printf("There were %d constraints and %d vsites in the simulation\n",
 524            nharm,nvsite);
 525     printf("Total correction to cV = %g J/mol K\n",qcvtot);
 526     printf("Total correction to  H = %g kJ/mol\n",qutot);
 527     ffclose(qc);
 528     please_cite(stdout,"Caleman2011b");
 529   }
 530   /* Writing eigenvectors. Note that if mass scaling was used, the eigenvectors
 531    * were scaled back from mass weighted cartesian to plain cartesian in the
 532    * nma_full_hessian() or nma_sparse_hessian() routines. Mass scaled vectors
 533    * will not be strictly orthogonal in plain cartesian scalar products.
 534    */
 535   write_eigenvectors(opt2fn("-v",NFILE,fnm),natoms,eigenvectors,FALSE,begin,end,
 536                      eWXR_NO,NULL,FALSE,top_x,bM,eigenvalues);
 537
 538   thanx(stderr);
 539
 540   return 0;
 541 }
 542