Updated. Fixed computation of frequency.

diff --git a/src/tools/gmx_dos.c b/src/tools/gmx_dos.c
index fcbdcabc19d766bb953637a8e13f80a99a1eb550..504af0a614ff6d784419f7eee6a93472e845c500 100644 (file)
--- a/src/tools/gmx_dos.c
+++ b/src/tools/gmx_dos.c
@@ -119,7 +119,6 @@ static real wCsolid(real nu,real beta)
 static real wSsolid(real nu,real beta)
 {
     real bhn = beta*PLANCK*nu;
-    real ebn;
     
     if (bhn == 0) 
     {
@@ -127,8 +126,6 @@ static real wSsolid(real nu,real beta)
     }
     else 
     {
-        ebn = exp(-bhn);
-        
         return bhn/(exp(bhn)-1) - log(1-exp(-bhn));
     }
 }
@@ -155,8 +152,8 @@ int gmx_dos(int argc,char *argv[])
     char       title[256];
     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,dostot,dos2,dosabs;
+    int        nV,nframes,n_alloc,i,j,k,l,fftcode,Nmol,Natom;
+    double     rho,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;
@@ -188,7 +185,7 @@ int gmx_dos(int argc,char *argv[])
         { efXVG, "-vacf", "vacf",  ffWRITE },
         { efXVG, "-mvacf","mvacf", ffWRITE },
         { efXVG, "-dos",  "dos",   ffWRITE },
-        { efLOG, "-dos",  "dos",   ffWRITE },
+        { efLOG, "-g",    "dos",   ffWRITE },
     };
 #define NFILE asize(fnm)
     int     npargs;
@@ -228,7 +225,7 @@ int gmx_dos(int argc,char *argv[])
     t0=fr.time;
       
     n_alloc=0;
-    teller=0;
+    nframes=0;
     Vsum = V2sum = 0;
     nV = 0;
     do {
@@ -239,7 +236,7 @@ int gmx_dos(int argc,char *argv[])
             Vsum += V;
             nV++;
         }
-        if (teller >= n_alloc) 
+        if (nframes >= n_alloc) 
         {
             n_alloc+=100;
             for(i=0; i<gnx; i++)
@@ -247,38 +244,38 @@ int gmx_dos(int argc,char *argv[])
         }
         for(i=0; i<gnx; i+=DIM) 
         {
-            c1[i+XX][teller] = fr.v[i/DIM][XX];
-            c1[i+YY][teller] = fr.v[i/DIM][YY];
-            c1[i+ZZ][teller] = fr.v[i/DIM][ZZ];
+            c1[i+XX][nframes] = fr.v[i/DIM][XX];
+            c1[i+YY][nframes] = fr.v[i/DIM][YY];
+            c1[i+ZZ][nframes] = fr.v[i/DIM][ZZ];
         }
 
         t1=fr.time;
 
-        teller ++;
+        nframes++;
     } while (read_next_frame(oenv,status,&fr));
   
     close_trj(status);
 
-    dt = (t1-t0)/(teller-1);
+    dt = (t1-t0)/(nframes-1);
     if (nV > 0)
     {
         V = Vsum/nV;
     }
     if (bVerbose)
         printf("Going to do %d fourier transforms of length %d. Hang on.\n",
-               gnx,teller);
-    low_do_autocorr(NULL,oenv,NULL,teller,gnx,teller,c1,dt,eacNormal,0,FALSE,
+               gnx,nframes);
+    low_do_autocorr(NULL,oenv,NULL,nframes,gnx,nframes,c1,dt,eacNormal,0,FALSE,
                     FALSE,FALSE,-1,-1,0,0);
     snew(dos,DOS_NR);
     for(j=0; (j<DOS_NR); j++)
-        snew(dos[j],teller+4);
+        snew(dos[j],nframes+4);
 
     if (bVerbose)
         printf("Going to merge the ACFs into the mass-weighted and plain ACF\n");
     for(i=0; (i<gnx); i+=DIM) 
     {
         mi = top.atoms.atom[i/DIM].m;
-        for(j=0; (j<teller/2); j++) 
+        for(j=0; (j<nframes/2); j++) 
         {
             c1j = (c1[i+XX][j] + c1[i+YY][j] + c1[i+ZZ][j]);
             dos[VACF][j]  += c1j/Natom;
@@ -287,8 +284,8 @@ int gmx_dos(int argc,char *argv[])
     }
     fp = xvgropen(opt2fn("-vacf",NFILE,fnm),"Velocity ACF",
                   "Time (ps)","C(t)",oenv);
-    snew(tt,teller/2);
-    for(j=0; (j<teller/2); j++) 
+    snew(tt,nframes/2);
+    for(j=0; (j<nframes/2); j++) 
     {
         tt[j] = j*dt;
         fprintf(fp,"%10g  %10g\n",tt[j],dos[VACF][j]);
@@ -296,13 +293,13 @@ int gmx_dos(int argc,char *argv[])
     fclose(fp);
     fp = xvgropen(opt2fn("-mvacf",NFILE,fnm),"Mass-weighted velocity ACF",
                   "Time (ps)","C(t)",oenv);
-    for(j=0; (j<teller/2); j++) 
+    for(j=0; (j<nframes/2); j++) 
     {
         fprintf(fp,"%10g  %10g\n",tt[j],dos[MVACF][j]);
     }
     fclose(fp);
     
-    if ((fftcode = gmx_fft_init_1d_real(&fft,teller/2,
+    if ((fftcode = gmx_fft_init_1d_real(&fft,nframes/2,
                                         GMX_FFT_FLAG_NONE)) != 0) 
     {
         gmx_fatal(FARGS,"gmx_fft_init_1d_real returned %d",fftcode);
@@ -318,10 +315,10 @@ int gmx_dos(int argc,char *argv[])
     /* Magic factor of 8 included now. */
     bfac = 8*dt*beta/2;
     dos2 = 0;
-    snew(nu,teller/2);
-    for(j=0; (j<teller/2); j++) 
+    snew(nu,nframes/4);
+    for(j=0; (j<nframes/4); j++) 
     {
-        nu[j] = j/(t1-t0);
+        nu[j] = 2*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]));
@@ -329,10 +326,10 @@ int gmx_dos(int argc,char *argv[])
             dos[DOS][j] = bfac*dos[DOS][2*j];
     }
     /* Normalize it */
-    dostot = evaluate_integral(teller/2,nu,dos[DOS],NULL,teller/2,&stddev);
+    dostot = evaluate_integral(nframes/4,nu,dos[DOS],NULL,nframes/4,&stddev);
     if (bNormalize) 
     {
-        for(j=0; (j<teller/2); j++) 
+        for(j=0; (j<nframes/4); j++) 
             dos[DOS][j] *= 3*Natom/dostot;
     }
     
@@ -343,10 +340,12 @@ int gmx_dos(int argc,char *argv[])
     Delta = ((2*DoS0/(9*Natom))*sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*
              pow((Natom/V),1.0/3.0)*pow(6/M_PI,2.0/3.0));
     f = calc_fluidicity(Delta,toler);
+    rho = (tmass*AMU)/(V*NANO*NANO*NANO);
     fprintf(fplog,"Nmol = %d\nNatom = %d\ndt = %g ps\ntmass = %g amu\n"
-            "V = %g nm^3\nDelta = %g\nf = %g\nbeta = %g (mol/kJ)\n"
+            "V = %g nm^3\nrho = %g g/l\n"
+            "Delta = %g\nf = %g\nbeta = %g (mol/kJ)\n"
             "DoS0 = %g\nDos2 = %g\nDoSTot = %g\n",
-            Nmol,Natom,dt,tmass,V,Delta,f,beta,DoS0,dos2,dostot);
+            Nmol,Natom,dt,tmass,V,rho,Delta,f,beta,DoS0,dos2,dostot);
            
     /* Now compute solid (2) and diffusive (3) components */
     fp = xvgropen(opt2fn("-dos",NFILE,fnm),"Density of states",
@@ -354,7 +353,7 @@ int gmx_dos(int argc,char *argv[])
                   "\\f{4}S(\\f{12}n\\f{4})",oenv);
     xvgr_legend(fp,asize(DoSlegend),DoSlegend,oenv);
     recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;
-    for(j=0; (j<teller/2); j++) 
+    for(j=0; (j<nframes/4); j++) 
     {
         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];
@@ -369,25 +368,28 @@ int gmx_dos(int argc,char *argv[])
     /* Finally analyze the results! */    
     wCdiff = 0.5;
     wSdiff = DoS0/(3*BOLTZ); /* Is this correct? */
-    for(j=0; (j<teller/2); j++) 
+    for(j=0; (j<nframes/4); 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));
     }
-    DiffCoeff = evaluate_integral(teller/2,tt,dos[VACF],NULL,teller/2,&stddev);
+    DiffCoeff = evaluate_integral(nframes/2,tt,dos[VACF],NULL,nframes/2,&stddev);
     DiffCoeff = 1000*DiffCoeff/3.0;
     fprintf(fplog,"Diffusion coefficient from VACF %g 10^-5 cm^2/s\n",
             DiffCoeff);
     fprintf(fplog,"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);
+    cP = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_CP],NULL,
+                                   nframes/4,&stddev);
     fprintf(fplog,"Heat capacity %g J/mol K\n",1000*cP/Nmol);
     
-    S  = BOLTZ * evaluate_integral(teller/2,nu,dos[DOS_S],NULL,teller/2,&stddev);
+    S  = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_S],NULL,
+                                   nframes/4,&stddev);
     fprintf(fplog,"Entropy %g J/mol K\n",1000*S/Nmol);
+    fprintf(fplog,"\nArrivederci!\n");
     fclose(fplog);
     
     do_view(oenv,ftp2fn(efXVG,NFILE,fnm),"-nxy");