Multiplied the DoS by the timestep (between frames).

diff --git a/src/tools/gmx_dos.c b/src/tools/gmx_dos.c
index 8a12d53f88a170e845c9fe5362a3a1b0801a8bf5..11bdaebf9f7a8dae3bbe89c82c6d9e8b114d367c 100644 (file)
--- a/src/tools/gmx_dos.c
+++ b/src/tools/gmx_dos.c
@@ -104,15 +104,15 @@ static real calc_fluidicity(real Delta,real tol)
 static real wCsolid(real nu,real beta)
 {
     real bhn = beta*PLANCK*nu;
-    real ebn;
+    real ebn,koko;
     
     if (bhn == 0)
-        return 1;
+        return 1.0;
     else 
     {
-        ebn = exp(-beta*PLANCK*nu);
-        
-        return ebn/sqr(1-ebn);
+        ebn = exp(bhn);
+        koko = sqr(1-ebn);
+        return sqr(bhn)*ebn/koko;
     }
 }
 
@@ -156,18 +156,25 @@ int gmx_dos(int argc,char *argv[])
     real       t0,t1,m;
     t_trxstatus *status;
     int        nV,teller,n_alloc,i,j,k,l,fftcode,Nmol,Natom;
-    double     dt,V2sum,Vsum,V,tmass;
-    real       **c1,**dos,mi,beta,bfac,*nu,*tt,stddev;
+    double     dt,V2sum,Vsum,V,tmass,dostot,dos2,dosabs;
+    real       **c1,**dos,mi,beta,bfac,*nu,*tt,stddev,c1j;
     output_env_t oenv;
     gmx_fft_t  fft;
     double     cP,S,DiffCoeff,Delta,f,DoS0;
     double     wCdiff,wSdiff;
     
-    static     gmx_bool bVerbose=TRUE;
+    static     gmx_bool bVerbose=TRUE,bAbsolute=FALSE,bNormalize=FALSE;
+    static     gmx_bool bRecip=FALSE;
     static     real Temp=298.15,toler=1e-3;
     t_pargs pa[] = {
         { "-v", FALSE, etBOOL, {&bVerbose},
-          "Be loud and noisy" },
+          "Be loud and noisy." },
+        { "-recip", FALSE, etBOOL, {&bRecip},
+          "Use cm^-1 on X-axis instead of 1/ps for DoS plots." },
+        { "-abs", FALSE, etBOOL, {&bAbsolute},
+          "Use the absolute value of the Fourier transform of the VACF as the Density of States. Default is to use the real component only" },
+        { "-normdos", FALSE, etBOOL, {&bNormalize},
+          "Normalize the DoS such that it adds up to 3N. This is a hack that should not be necessary." },
         { "-T", FALSE, etREAL, {&Temp},
           "Temperature in the simulation" },
         { "-toler", FALSE, etREAL, {&toler},
@@ -262,7 +269,7 @@ int gmx_dos(int argc,char *argv[])
                     FALSE,FALSE,-1,-1,0,0);
     snew(dos,DOS_NR);
     for(j=0; (j<DOS_NR); j++)
-        snew(dos[j],teller);
+        snew(dos[j],teller+4);
 
     if (bVerbose)
         printf("Going to merge the ACFs into the mass-weighted and plain ACF\n");
@@ -271,8 +278,9 @@ int gmx_dos(int argc,char *argv[])
         mi = top.atoms.atom[i/DIM].m;
         for(j=0; (j<teller/2); j++) 
         {
-            dos[VACF][j]  += (c1[i+XX][j] + c1[i+YY][j] + c1[i+ZZ][j])/Natom;
-            dos[MVACF][j] += mi * (c1[i+XX][j] + c1[i+YY][j] + c1[i+ZZ][j]);
+            c1j = (c1[i+XX][j] + c1[i+YY][j] + c1[i+ZZ][j]);
+            dos[VACF][j]  += c1j/Natom;
+            dos[MVACF][j] += mi*c1j;
         }
     }
     fp = xvgropen(opt2fn("-vacf",NFILE,fnm),"Velocity ACF",
@@ -287,10 +295,13 @@ int gmx_dos(int argc,char *argv[])
     fp = xvgropen(opt2fn("-mvacf",NFILE,fnm),"Mass-weighted velocity ACF",
                   "Time (ps)","C(t)",oenv);
     for(j=0; (j<teller/2); j++) 
+    {
         fprintf(fp,"%10g  %10g\n",tt[j],dos[MVACF][j]);
+    }
     fclose(fp);
     
-    if ((fftcode = gmx_fft_init_1d_real(&fft,teller,GMX_FFT_FLAG_NONE)) != 0) 
+    if ((fftcode = gmx_fft_init_1d_real(&fft,teller/2,
+                                        GMX_FFT_FLAG_NONE)) != 0) 
     {
         gmx_fatal(FARGS,"gmx_fft_init_1d_real returned %d",fftcode);
     }
@@ -302,32 +313,49 @@ int gmx_dos(int argc,char *argv[])
 
     /* First compute the DoS */
     beta = 1/(Temp*BOLTZ);
-    bfac = (beta/2);
+    bfac = dt*beta/2;
+    dos2 = 0;
+    snew(nu,teller/2);
     for(j=0; (j<teller/2); j++) 
     {
-        dos[DOS][j] = bfac*dos[DOS][2*j];
+        nu[j] = j/(t1-t0);
+        dos2 += sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]);
+        if (bAbsolute)
+            dos[DOS][j] = bfac*sqrt(sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]));
+        else
+            dos[DOS][j] = bfac*dos[DOS][2*j];
+    }
+    /* Normalize it */
+    dostot = evaluate_integral(teller/2,nu,dos[DOS],NULL,teller/2,&stddev);
+    if (bNormalize) 
+    {
+        for(j=0; (j<teller/2); j++) 
+            dos[DOS][j] *= 3*Natom/dostot;
     }
+    
     /* Now analyze it */
     DoS0 = dos[DOS][0];
-    Delta = ((2*DoS0/(9*Natom))*pow(M_PI*BOLTZ*Temp/tmass,1.0/3.0)*
+    
+    /* Note this eqn. is incorrect in Pascal2011a! */
+    Delta = ((2*DoS0/(9*Natom))*sqrt(M_PI*BOLTZ*Temp/tmass)*
              pow((Natom/V),1.0/3.0)*pow(6/M_PI,2.0/3.0));
     f = calc_fluidicity(Delta,toler);
-    printf("Nmol = %d, Natom = %d, dt = %g ps\n"
-           "V = %g nm^3, Delta = %g, f = %g\n",
-           Nmol,Natom,dt,V,Delta,f);
+    printf("Nmol = %d, Natom = %d, dt = %g ps, tmass = %g amu\n"
+           "V = %g nm^3, Delta = %g, f = %g beta = %g (mol/kJ)\n"
+           "DoS0 = %g, Dos2 = %g, DoSTot = %g\n",
+           Nmol,Natom,dt,tmass,V,Delta,f,beta,DoS0,dos2,dostot);
            
-    snew(nu,teller/2);
     /* Now compute solid (2) and diffusive (3) components */
     fp = xvgropen(opt2fn("-dos",NFILE,fnm),"Density of states",
-                  "E (cm\\S-1\\N)","\\f{4}S(\\f{12}n\\f{4})",oenv);
+                  bRecip ? "E (cm\\S-1\\N)" : "\\f{12}n\\f{4} (1/ps)",
+                  "\\f{4}S(\\f{12}n\\f{4})",oenv);
     xvgr_legend(fp,asize(DoSlegend),DoSlegend,oenv);
     for(j=0; (j<teller/2); j++) 
     {
-        nu[j] = j/(t1-t0);
         dos[DOS_DIFF][j] = DoS0/(1+sqr(DoS0*M_PI*nu[j]/(6*f*Natom)));
         dos[DOS_SOLID][j] = dos[DOS][j]-dos[DOS_DIFF][j];
         fprintf(fp,"%10g  %10g  %10g  %10g\n",
-                nu[j]/(1e-7*SPEED_OF_LIGHT),
+                bRecip ? nu[j]/(1e-7*SPEED_OF_LIGHT) : nu[j],
                 dos[DOS][j],dos[DOS_SOLID][j],dos[DOS_DIFF][j]);
     }
     fclose(fp);
@@ -335,21 +363,28 @@ int gmx_dos(int argc,char *argv[])
     /* Finally analyze the results! */    
     wCdiff = 0.5;
     wSdiff = DoS0/(3*BOLTZ); /* Is this correct? */
+/*    fp = xvgropen("test.xvg","Cp(v)","v","Cp",oenv);*/
     for(j=0; (j<teller/2); j++) 
     {
-        dos[DOS_CP][j] = dos[DOS_DIFF][j]*wCdiff + dos[DOS_SOLID][j]*wCsolid(nu[j],beta);
-        dos[DOS_S][j] = dos[DOS_DIFF][j]*wSdiff + dos[DOS_SOLID][j]*wSsolid(nu[j],beta);
-    }
+        dos[DOS_CP][j] = (dos[DOS_DIFF][j]*wCdiff + 
+                          dos[DOS_SOLID][j]*wCsolid(nu[j],beta));
+        dos[DOS_S][j]  = (dos[DOS_DIFF][j]*wSdiff + 
+                          dos[DOS_SOLID][j]*wSsolid(nu[j],beta));
+/*        fprintf(fp,"%10g  %10g  %10g\n",nu[j],
+                dos[DOS_CP][j],wCsolid(nu[j],beta));
+*/  }
+    /*  fclose(fp);*/
     DiffCoeff = evaluate_integral(teller/2,tt,dos[VACF],NULL,teller/2,&stddev);
     DiffCoeff = 1000*DiffCoeff/3.0;
     printf("Diffusion coefficient from VACF %g 10^-5 cm^2/s\n",DiffCoeff);
-    printf("Diffusion coefficient from DoS %g 10^-5 cm^2/s\n",1000*dos[DOS][0]/(12*tmass*beta));
+    printf("Diffusion coefficient from DoS %g 10^-5 cm^2/s\n",
+           1000*DoS0/(12*tmass*beta));
     
     cP = BOLTZ * evaluate_integral(teller/2,nu,dos[DOS_CP],NULL,teller/2,&stddev);
-    printf("Heat capacity %g J/mol K\n",1000*cP);
+    printf("Heat capacity %g J/mol K\n",1000*cP/Nmol);
     
     S  = BOLTZ * evaluate_integral(teller/2,nu,dos[DOS_S],NULL,teller/2,&stddev);
-    printf("Entropy %g J/mol K\n",1000*S);
+    printf("Entropy %g J/mol K\n",1000*S/Nmol);
     
     do_view(oenv,ftp2fn(efXVG,NFILE,fnm),"-nxy");