From 5cd759f24e237f2e90e88895e5f0c4d37d421ea4 Mon Sep 17 00:00:00 2001
From: David van der Spoel <spoel@xray.bmc.uu.se>
Date: Tue, 28 Feb 2012 21:53:28 +0100
Subject: [PATCH] Vibration spectra from velocity ACF or Normal modes.

The velocity autocorrelation function can be converted into
a vibrational spectrum by performing a fourier transform
and plotting the amplitude squared as a function of
frequency or wavenumber. This has now been automated in
g_velacc with the -os flag. Likewise, an option -os
has been added to the g_nmeig program to produce a spectrum
based on the harmonic approximation.

Change-Id: I7e3d669b9108e3a91856ff20688f6216a035bf4a
---
 src/tools/gmx_nmeig.c  |  56 +++++-
 src/tools/gmx_velacc.c | 390 ++++++++++++++++++++++++-----------------
 2 files changed, 278 insertions(+), 168 deletions(-)

diff --git a/src/tools/gmx_nmeig.c b/src/tools/gmx_nmeig.c
index 13793d1809..febc81f5da 100644
--- a/src/tools/gmx_nmeig.c
+++ b/src/tools/gmx_nmeig.c
@@ -1,4 +1,4 @@
-/*  -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+/*
  * 
  *                This source code is part of
  * 
@@ -284,8 +284,8 @@ int gmx_nmeig(int argc,char *argv[])
   };
     
   static gmx_bool bM=TRUE,bCons=FALSE;
-  static int  begin=1,end=50;
-  static real T=298.15;
+  static int  begin=1,end=50,maxspec=4000;
+  static real T=298.15,width=1;
   t_pargs pa[] = 
   {
     { "-m",  FALSE, etBOOL, {&bM},
@@ -296,12 +296,16 @@ int gmx_nmeig(int argc,char *argv[])
       "First eigenvector to write away" },
     { "-last",  FALSE, etINT, {&end}, 
       "Last eigenvector to write away" },
+    { "-maxspec", FALSE, etINT, {&maxspec},
+      "Highest frequency (1/cm) to consider in the spectrum" },
     { "-T",     FALSE, etREAL, {&T},
       "Temperature for computing quantum heat capacity and enthalpy when using normal mode calculations to correct classical simulations" },
     { "-constr", FALSE, etBOOL, {&bCons},
       "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." },
+    { "-width",  FALSE, etREAL, {&width},
+      "Width (sigma) of the gaussian peaks (1/cm) when generating a spectrum" }
   };
-  FILE       *out,*qc;
+  FILE       *out,*qc,*spec;
   int        status,trjout;
   t_topology top;
   gmx_mtop_t mtop;
@@ -321,10 +325,12 @@ int gmx_nmeig(int argc,char *argv[])
   real       value,omega,nu;
   real       factor_gmx_to_omega2;
   real       factor_omega_to_wavenumber;
+  real       *spectrum=NULL;
+  real       wfac;
   t_commrec  *cr;
   output_env_t oenv;
   const char *qcleg[] = { "Heat Capacity cV (J/mol K)", 
-			  "Enthalpy H (kJ/mol)" };
+                          "Enthalpy H (kJ/mol)" };
   real *                 full_hessian   = NULL;
   gmx_sparsematrix_t *   sparse_hessian = NULL;
 
@@ -333,9 +339,10 @@ int gmx_nmeig(int argc,char *argv[])
     { efTPX, NULL, NULL,         ffREAD  },
     { efXVG, "-of", "eigenfreq", ffWRITE },
     { efXVG, "-ol", "eigenval",  ffWRITE },
+    { efXVG, "-os", "spectrum",  ffOPTWR },
     { efXVG, "-qc", "quant_corr",  ffOPTWR },
     { efTRN, "-v", "eigenvec",  ffWRITE }
-  }; 
+  };
 #define NFILE asize(fnm) 
 
   cr = init_par(&argc,&argv);
@@ -449,7 +456,7 @@ int gmx_nmeig(int argc,char *argv[])
                       
   /* now write the output */
   fprintf (stderr,"Writing eigenvalues...\n");
-  out=xvgropen(opt2fn("-ol",NFILE,fnm), 
+  out=xvgropen(opt2fn("-ol",NFILE,fnm),
                "Eigenvalues","Eigenvalue index","Eigenvalue [Gromacs units]",
                oenv);
   if (output_env_get_print_xvgr_codes(oenv)) {
@@ -462,7 +469,7 @@ int gmx_nmeig(int argc,char *argv[])
   for (i=0; i<=(end-begin); i++)
       fprintf (out,"%6d %15g\n",begin+i,eigenvalues[i]);
   ffclose(out);
-  
+
 
   if (opt2bSet("-qc",NFILE,fnm)) {
     qc = xvgropen(opt2fn("-qc",NFILE,fnm),"Quantum Corrections","Eigenvector index","",oenv);
@@ -482,7 +489,22 @@ int gmx_nmeig(int argc,char *argv[])
     else 
       fprintf(out,"@ subtitle \"not mass weighted\"\n");
   }
-  
+  /* Spectrum ? */
+  spec = NULL;
+  if (opt2bSet("-os",NFILE,fnm) && (maxspec > 0))
+  {
+      snew(spectrum,maxspec);
+      spec=xvgropen(opt2fn("-os",NFILE,fnm), 
+                    "Vibrational spectrum based on harmonic approximation",
+                    "\\f{12}w\\f{4} (cm\\S-1\\N)",
+                    "Intensity [Gromacs units]",
+                    oenv);
+      for(i=0; (i<maxspec); i++)
+      {
+          spectrum[i] = 0;
+      }
+  }
+
   /* Gromacs units are kJ/(mol*nm*nm*amu),
    * where amu is the atomic mass unit.
    *
@@ -503,6 +525,14 @@ int gmx_nmeig(int argc,char *argv[])
       nu    = 1e-12*omega/(2*M_PI);
       value = omega*factor_omega_to_wavenumber;
       fprintf (out,"%6d %15g\n",i,value);
+      if (NULL != spec)
+      {
+          wfac = eigenvalues[i-begin]/(width*sqrt(2*M_PI));
+          for(j=0; (j<maxspec); j++)
+          {
+              spectrum[j] += wfac*exp(-sqr(j-value)/(2*sqr(width)));
+          }
+      }
       if (NULL != qc) {
           qcv = cv_corr(nu,T);
           qu  = u_corr(nu,T);
@@ -517,6 +547,14 @@ int gmx_nmeig(int argc,char *argv[])
       }
   }
   ffclose(out);
+  if (NULL != spec)
+  {
+      for(j=0; (j<maxspec); j++)
+      {
+          fprintf(spec,"%10g  %10g\n",1.0*j,spectrum[j]);
+      }
+      ffclose(spec);
+  }
   if (NULL != qc) {
     printf("Quantum corrections for harmonic degrees of freedom\n");
     printf("Use appropriate -first and -last options to get reliable results.\n");
diff --git a/src/tools/gmx_velacc.c b/src/tools/gmx_velacc.c
index 0151928643..ee3083a3d8 100644
--- a/src/tools/gmx_velacc.c
+++ b/src/tools/gmx_velacc.c
@@ -57,196 +57,268 @@
 #include "strdb.h"
 #include "xvgr.h"
 #include "gmx_ana.h"
-
+#include "gmx_fft.h"
 
 static void index_atom2mol(int *n,atom_id *index,t_block *mols)
 {
-  int nat,i,nmol,mol,j;
-
-  nat = *n;
-  i = 0;
-  nmol = 0;
-  mol = 0;
-  while (i < nat) {
-    while (index[i] > mols->index[mol]) {
-      mol++;
-      if (mol >= mols->nr)
-	gmx_fatal(FARGS,"Atom index out of range: %d",index[i]+1);
-    }
-    for(j=mols->index[mol]; j<mols->index[mol+1]; j++) {
-      if (i >= nat || index[i] != j)
-      	gmx_fatal(FARGS,"The index group does not consist of whole molecules");
-      i++;
+    int nat,i,nmol,mol,j;
+
+    nat = *n;
+    i = 0;
+    nmol = 0;
+    mol = 0;
+    while (i < nat) {
+        while (index[i] > mols->index[mol]) {
+            mol++;
+            if (mol >= mols->nr)
+                gmx_fatal(FARGS,"Atom index out of range: %d",index[i]+1);
+        }
+        for(j=mols->index[mol]; j<mols->index[mol+1]; j++) {
+            if (i >= nat || index[i] != j)
+                gmx_fatal(FARGS,"The index group does not consist of whole molecules");
+            i++;
+        }
+        index[nmol++] = mol;
     }
-    index[nmol++] = mol;
-  }
 
-  fprintf(stderr,"\nSplit group of %d atoms into %d molecules\n",nat,nmol);
+    fprintf(stderr,"\nSplit group of %d atoms into %d molecules\n",nat,nmol);
 
-  *n = nmol;
+    *n = nmol;
 }
 
 static void precalc(t_topology top,real normm[]){
 
-  real mtot;
-  int i,j,k,l;
+    real mtot;
+    int i,j,k,l;
 
-  for(i=0;i<top.mols.nr;i++){
-    k=top.mols.index[i];
-    l=top.mols.index[i+1];
-    mtot=0.0;
+    for(i=0;i<top.mols.nr;i++){
+        k=top.mols.index[i];
+        l=top.mols.index[i+1];
+        mtot=0.0;
 
-    for(j=k;j<l;j++)
-      mtot+=top.atoms.atom[j].m;
+        for(j=k;j<l;j++)
+            mtot+=top.atoms.atom[j].m;
 
-    for(j=k;j<l;j++)
-      normm[j]=top.atoms.atom[j].m/mtot;
+        for(j=k;j<l;j++)
+            normm[j]=top.atoms.atom[j].m/mtot;
 
-  }
+    }
 
 }
 
+static void calc_spectrum(int n,real c[],real dt,const char *fn,
+                          output_env_t oenv,gmx_bool bRecip)
+{
+    FILE *fp;
+    gmx_fft_t fft;
+    int  i,status;
+    real *data;
+    real nu,omega,recip_fac;
 
+    snew(data,n*2);
+    for(i=0; (i<n); i++)
+        data[i] = c[i];
+
+    if ((status = gmx_fft_init_1d_real(&fft,n,GMX_FFT_FLAG_NONE)) != 0)
+    {
+        gmx_fatal(FARGS,"Invalid fft return status %d",status);
+    }
+    if ((status = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX,data,data)) != 0)
+    {
+        gmx_fatal(FARGS,"Invalid fft return status %d",status);
+    }
+    fp = xvgropen(fn,"Vibrational Power Spectrum",
+                  bRecip ? "\\f{12}w\\f{4} (cm\\S-1\\N)" :
+                  "\\f{12}n\\f{4} (ps\\S-1\\N)",
+                  "a.u.",oenv);
+    /* This is difficult.
+     * The length of the ACF is dt (as passed to this routine).
+     * We pass the vacf with N time steps from 0 to dt.
+     * That means that after FFT we have lowest frequency = 1/dt
+     * then 1/(2 dt) etc. (this is the X-axis of the data after FFT).
+     * To convert to 1/cm we need to have to realize that
+     * E = hbar w = h nu = h c/lambda. We want to have reciprokal cm
+     * on the x-axis, that is 1/lambda, so we then have
+     * 1/lambda = nu/c. Since nu has units of 1/ps and c has gromacs units
+     * of nm/ps, we need to multiply by 1e7.
+     * The timestep between saving the trajectory is
+     * 1e7 is to convert nanometer to cm
+     */
+    recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;
+    for(i=0; (i<n); i+=2) 
+    {
+        nu = i/(2*dt);
+        omega = nu*recip_fac;
+        /* Computing the square magnitude of a complex number, since this is a power
+         * spectrum.
+         */
+        fprintf(fp,"%10g  %10g\n",omega,sqr(data[i])+sqr(data[i+1]));
+    }
+    xvgrclose(fp);
+    gmx_fft_destroy(fft);
+    sfree(data);
+}
 
 int gmx_velacc(int argc,char *argv[])
 {
-  const char *desc[] = {
-    "[TT]g_velacc[tt] computes the velocity autocorrelation function.",
-    "When the [TT]-m[tt] option is used, the momentum autocorrelation",
-    "function is calculated.[PAR]",
-    "With option [TT]-mol[tt] the velocity autocorrelation function of",
-    "molecules is calculated. In this case the index group should consist",
-    "of molecule numbers instead of atom numbers.[PAR]",
-    "Be sure that your trajectory contains frames with velocity information",
-    "(i.e. [TT]nstvout[tt] was set in your original [TT].mdp[tt] file),",
-    "and that the time interval between data collection points is",
-    "much shorter than the time scale of the autocorrelation."
-  };
+    const char *desc[] = {
+        "[TT]g_velacc[tt] computes the velocity autocorrelation function.",
+        "When the [TT]-m[tt] option is used, the momentum autocorrelation",
+        "function is calculated.[PAR]",
+        "With option [TT]-mol[tt] the velocity autocorrelation function of",
+        "molecules is calculated. In this case the index group should consist",
+        "of molecule numbers instead of atom numbers.[PAR]",
+        "Be sure that your trajectory contains frames with velocity information",
+        "(i.e. [TT]nstvout[tt] was set in your original [TT].mdp[tt] file),",
+        "and that the time interval between data collection points is",
+        "much shorter than the time scale of the autocorrelation."
+    };
   
-  static gmx_bool bM=FALSE,bMol=FALSE;
-  t_pargs pa[] = {
-    { "-m", FALSE, etBOOL, {&bM},
-      "Calculate the momentum autocorrelation function" },
-    { "-mol", FALSE, etBOOL, {&bMol},
-      "Calculate the velocity acf of molecules" }
-  };
-
-  t_topology top;
-  int        ePBC=-1;
-  t_trxframe fr;
-  matrix     box;
-  gmx_bool       bTPS=FALSE,bTop=FALSE;
-  int        gnx;
-  atom_id    *index;
-  char       *grpname;
-  char       title[256];
-  real       t0,t1,m;
-  t_trxstatus *status;
-  int        teller,n_alloc,i,j,tel3,k,l;
-  rvec       mv_mol;
-  real       **c1;
-  real	     *normm=NULL;
-  output_env_t oenv;
+    static gmx_bool bMass=FALSE,bMol=FALSE,bRecip=TRUE;
+    t_pargs pa[] = {
+        { "-m", FALSE, etBOOL, {&bMass},
+          "Calculate the momentum autocorrelation function" },
+        { "-recip", FALSE, etBOOL, {&bRecip},
+          "Use cm^-1 on X-axis instead of 1/ps for spectra." },
+        { "-mol", FALSE, etBOOL, {&bMol},
+          "Calculate the velocity acf of molecules" }
+    };
+
+    t_topology top;
+    int        ePBC=-1;
+    t_trxframe fr;
+    matrix     box;
+    gmx_bool       bTPS=FALSE,bTop=FALSE;
+    int        gnx;
+    atom_id    *index;
+    char       *grpname;
+    char       title[256];
+    /* t0, t1 are the beginning and end time respectively. 
+     * dt is the time step, mass is temp variable for atomic mass.
+     */
+    real       t0,t1,dt,mass;
+    t_trxstatus *status;
+    int        counter,n_alloc,i,j,counter_dim,k,l;
+    rvec       mv_mol;
+    /* Array for the correlation function */
+    real       **c1;
+    real             *normm=NULL;
+    output_env_t oenv;
   
 #define NHISTO 360
     
-  t_filenm  fnm[] = {
-    { efTRN, "-f",    NULL,   ffREAD  },
-    { efTPS, NULL,    NULL,   ffOPTRD }, 
-    { efNDX, NULL,    NULL,   ffOPTRD },
-    { efXVG, "-o",    "vac",  ffWRITE }
-  };
+    t_filenm  fnm[] = {
+        { efTRN, "-f",    NULL,   ffREAD  },
+        { efTPS, NULL,    NULL,   ffOPTRD }, 
+        { efNDX, NULL,    NULL,   ffOPTRD },
+        { efXVG, "-o",    "vac",  ffWRITE },
+        { efXVG, "-os",   "spectrum", ffOPTWR }
+    };
 #define NFILE asize(fnm)
-  int     npargs;
-  t_pargs *ppa;
-
-  CopyRight(stderr,argv[0]);
-  npargs = asize(pa);
-  ppa    = add_acf_pargs(&npargs,pa);
-  parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
-		    NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
-
-  if (bMol || bM) {
-    bTPS = ftp2bSet(efTPS,NFILE,fnm) || !ftp2bSet(efNDX,NFILE,fnm);
-  }
-
-  if (bTPS) {
-    bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,NULL,NULL,box,
-		       TRUE);
-    get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
-  } else
-    rd_index(ftp2fn(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
-
-  if (bMol) {
-    if (!bTop)
-      gmx_fatal(FARGS,"Need a topology to determine the molecules");
-    snew(normm,top.atoms.nr);
-    precalc(top,normm);
-    index_atom2mol(&gnx,index,&top.mols);
-  }
+    int     npargs;
+    t_pargs *ppa;
+
+    CopyRight(stderr,argv[0]);
+    npargs = asize(pa);
+    ppa    = add_acf_pargs(&npargs,pa);
+    parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
+                      NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
+
+    if (bMol || bMass) {
+        bTPS = ftp2bSet(efTPS,NFILE,fnm) || !ftp2bSet(efNDX,NFILE,fnm);
+    }
+
+    if (bTPS) {
+        bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,NULL,NULL,box,
+                           TRUE);
+        get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
+    } else
+        rd_index(ftp2fn(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
+
+    if (bMol) {
+        if (!bTop)
+            gmx_fatal(FARGS,"Need a topology to determine the molecules");
+        snew(normm,top.atoms.nr);
+        precalc(top,normm);
+        index_atom2mol(&gnx,index,&top.mols);
+    }
   
-  /* Correlation stuff */
-  snew(c1,gnx);
-  for(i=0; (i<gnx); i++)
-    c1[i]=NULL;
+    /* Correlation stuff */
+    snew(c1,gnx);
+    for(i=0; (i<gnx); i++)
+        c1[i]=NULL;
   
-  read_first_frame(oenv,&status,ftp2fn(efTRN,NFILE,fnm),&fr,TRX_NEED_V);
-  t0=fr.time;
+    read_first_frame(oenv,&status,ftp2fn(efTRN,NFILE,fnm),&fr,TRX_NEED_V);
+    t0=fr.time;
       
-  n_alloc=0;
-  teller=0;
-  do {
-    if (teller >= n_alloc) {
-      n_alloc+=100;
-      for(i=0; i<gnx; i++)
-	srenew(c1[i],DIM*n_alloc);
-    }
-    tel3=3*teller;
-    if (bMol)
-      for(i=0; i<gnx; i++) {
-	clear_rvec(mv_mol);
-	k=top.mols.index[index[i]];
-	l=top.mols.index[index[i]+1];
-	for(j=k; j<l; j++) {
-   	  if (bM)
-	    m = top.atoms.atom[j].m;
-	  else
-	    m = normm[j];
-	  mv_mol[XX] += m*fr.v[j][XX];
-	  mv_mol[YY] += m*fr.v[j][YY];
-	  mv_mol[ZZ] += m*fr.v[j][ZZ];
-	}
-	c1[i][tel3+XX]=mv_mol[XX];
-	c1[i][tel3+YY]=mv_mol[YY];
-	c1[i][tel3+ZZ]=mv_mol[ZZ];
-      }
-     else
-      for(i=0; i<gnx; i++) {
-        if (bM)
-	  m = top.atoms.atom[index[i]].m;
-	else
-	 m = 1;
-	c1[i][tel3+XX]=m*fr.v[index[i]][XX];
-	c1[i][tel3+YY]=m*fr.v[index[i]][YY];
-	c1[i][tel3+ZZ]=m*fr.v[index[i]][ZZ];
-      }
+    n_alloc=0;
+    counter=0;
+    do {
+        if (counter >= n_alloc) {
+            n_alloc+=100;
+            for(i=0; i<gnx; i++)
+                srenew(c1[i],DIM*n_alloc);
+        }
+        counter_dim=DIM*counter;
+        if (bMol)
+            for(i=0; i<gnx; i++) {
+                clear_rvec(mv_mol);
+                k=top.mols.index[index[i]];
+                l=top.mols.index[index[i]+1];
+                for(j=k; j<l; j++) {
+                    if (bMass)
+                        mass = top.atoms.atom[j].m;
+                    else
+                        mass = normm[j];
+                    mv_mol[XX] += mass*fr.v[j][XX];
+                    mv_mol[YY] += mass*fr.v[j][YY];
+                    mv_mol[ZZ] += mass*fr.v[j][ZZ];
+                }
+                c1[i][counter_dim+XX]=mv_mol[XX];
+                c1[i][counter_dim+YY]=mv_mol[YY];
+                c1[i][counter_dim+ZZ]=mv_mol[ZZ];
+            }
+        else
+            for(i=0; i<gnx; i++) {
+                if (bMass)
+                    mass = top.atoms.atom[index[i]].m;
+                else
+                    mass = 1;
+                c1[i][counter_dim+XX]=mass*fr.v[index[i]][XX];
+                c1[i][counter_dim+YY]=mass*fr.v[index[i]][YY];
+                c1[i][counter_dim+ZZ]=mass*fr.v[index[i]][ZZ];
+            }
 
-    t1=fr.time;
+        t1=fr.time;
 
-    teller ++;
-  } while (read_next_frame(oenv,status,&fr));
+        counter ++;
+    } while (read_next_frame(oenv,status,&fr));
   
-	close_trj(status);
+    close_trj(status);
 
-  do_autocorr(ftp2fn(efXVG,NFILE,fnm), oenv,
-	      bM ? 
-	      "Momentum Autocorrelation Function" :
-	      "Velocity Autocorrelation Function",
-	      teller,gnx,c1,(t1-t0)/(teller-1),eacVector,TRUE);
-  
-  do_view(oenv,ftp2fn(efXVG,NFILE,fnm),"-nxy");
-  
-  thanx(stderr);
-  
-  return 0;
+    if (counter >= 4)
+    {
+      /* Compute time step between frames */
+      dt = (t1-t0)/(counter-1);
+      do_autocorr(opt2fn("-o",NFILE,fnm), oenv,
+		  bMass ? 
+		  "Momentum Autocorrelation Function" :
+		  "Velocity Autocorrelation Function",
+		  counter,gnx,c1,dt,eacVector,TRUE);
+
+      do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy");
+
+      if (opt2bSet("-os",NFILE,fnm)) {
+        calc_spectrum(counter/2,(real *) (c1[0]),(t1-t0)/2,opt2fn("-os",NFILE,fnm),
+                      oenv,bRecip);
+        do_view(oenv,opt2fn("-os",NFILE,fnm),"-nxy");
+      }
+    }
+    else {
+      fprintf(stderr,"Not enough frames in trajectory - no output generated.\n");
+    }
+
+    thanx(stderr);
+
+    return 0;
 }
-- 
2.22.0