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;
}
}
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},
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");
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",
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);
}
/* 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);
/* 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");