2*f - 2);
}
+static real YYY(real f,real y)
+{
+ return (2*pow(y*f,3) - sqr(f)*y*(1+6*y) +
+ (2+6*y)*f - 2);
+}
+
+static real calc_compress(real y)
+{
+ if (y==1)
+ return 0;
+ return ((1+y+sqr(y)-pow(y,3))/(pow(1-y,3)));
+}
+
+static real bisector(real Delta,real tol,
+ real ff0,real ff1,
+ real ff(real,real))
+{
+ real fd0,fd,fd1,f,f0,f1;
+ real tolmin = 1e-6;
+
+ f0 = ff0;
+ f1 = ff1;
+ if (tol < tolmin)
+ {
+ fprintf(stderr,"Unrealistic tolerance %g for bisector. Setting it to %g\n",tol,tolmin);
+ tol=1e-6;
+ }
+
+ do {
+ fd0 = ff(Delta,f0);
+ fd1 = ff(Delta,f1);
+ f = (f0+f1)*0.5;
+ fd = ff(Delta,f);
+ if (fd < 0)
+ f0 = f;
+ else if (fd > 0)
+ f1 = f;
+ else
+ return f;
+ } while ((f1-f0) > tol);
+
+ return f;
+}
+
static real calc_fluidicity(real Delta,real tol)
+{
+ return bisector(Delta,tol,0,1,FD);
+}
+
+static real calc_y(real f,real Delta)
+{
+ return pow(f/Delta,1.5);
+/* return bisector(f,1e-5,0,10000,YYY);*/
+}
+
+static real calc_Shs(real f,real y)
+{
+ real fy = f*y;
+
+ return BOLTZ*(log(calc_compress(fy)) + fy*(3*fy-4)/sqr(1-fy));
+}
+
+static real old_calc_fluidicity(real Delta,real tol)
{
real fd0,fd,fd1,f,f0=0,f1=1;
real tolmin = 1e-6;
}
}
+static void dump_fy(output_env_t oenv,real toler)
+{
+ FILE *fp;
+ double Delta,f,y,DD;
+ const char *leg[] = { "f", "fy", "y" };
+
+ DD = pow(10.0,0.125);
+ fp=xvgropen("fy.xvg","Fig. 2, Lin2003a","Delta","y or fy",oenv);
+ xvgr_legend(fp,asize(leg),leg,oenv);
+ fprintf(fp,"@ world 1e-05, 0, 1000, 1\n");
+ fprintf(fp,"@ xaxes scale Logarithmic\n");
+ for (Delta=1e-5; (Delta <= 1000); Delta*=DD)
+ {
+ f = calc_fluidicity(Delta,toler);
+ y = calc_y(f,Delta);
+ fprintf(fp,"%10g %10g %10g %10g\n",Delta,f,f*y,y);
+ }
+ fclose(fp);
+}
+
+static void dump_w(output_env_t oenv,real beta)
+{
+ FILE *fp;
+ double nu;
+ const char *leg[] = { "wCv", "wS" };
+
+ fp=xvgropen("w.xvg","Fig. 1, Berens1983a","\\f{12}b\\f{4}h\\f{12}n",
+ "w",oenv);
+ xvgr_legend(fp,asize(leg),leg,oenv);
+ for (nu=1; (nu<100); nu+=0.05)
+ {
+ fprintf(fp,"%10g %10g %10g\n",beta*PLANCK*nu,
+ wCsolid(nu,beta),wSsolid(nu,beta));
+ }
+ fclose(fp);
+}
+
int gmx_dos(int argc,char *argv[])
{
const char *desc[] = {
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,recip_fac;
+ double cP,S,DiffCoeff,Delta,f,y,z,sigHS,Shs,Sig,DoS0,recip_fac;
double wCdiff,wSdiff;
static gmx_bool bVerbose=TRUE,bAbsolute=FALSE,bNormalize=FALSE;
- static gmx_bool bRecip=FALSE;
- static real Temp=298.15,toler=1e-3;
+ static gmx_bool bRecip=FALSE,bDump=FALSE;
+ static real Temp=298.15,toler=1e-6;
t_pargs pa[] = {
{ "-v", FALSE, etBOOL, {&bVerbose},
"Be loud and noisy." },
{ "-T", FALSE, etREAL, {&Temp},
"Temperature in the simulation" },
{ "-toler", FALSE, etREAL, {&toler},
- "[HIDDEN]Tolerance when computing the fluidicity using bisection algorithm" }
+ "[HIDDEN]Tolerance when computing the fluidicity using bisection algorithm" },
+ { "-dump", FALSE, etBOOL, {&bDump},
+ "[HIDDEN]Dump the y/fy plot corresponding to Fig. 2 inLin2003a and the and the weighting functions corresponding to Fig. 1 in Berens1983a." }
};
t_filenm fnm[] = {
parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
NFILE,fnm,npargs,ppa,asize(desc),desc,
asize(bugs),bugs,&oenv);
+
+ beta = 1/(Temp*BOLTZ);
+ if (bDump)
+ {
+ printf("Dumping reference figures. Thanks for your patience.\n");
+ dump_fy(oenv,toler);
+ dump_w(oenv,beta);
+ exit(0);
+ }
+
fplog = gmx_fio_fopen(ftp2fn(efLOG,NFILE,fnm),"w");
fprintf(fplog,"Doing density of states analysis based on trajectory.\n");
please_cite(fplog,"Pascal2011a");
}
/* First compute the DoS */
- beta = 1/(Temp*BOLTZ);
/* Magic factor of 8 included now. */
bfac = 8*dt*beta/2;
dos2 = 0;
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);
+ y = calc_y(f,Delta);
+ z = calc_compress(y);
+ Sig = 5.0/3.0;
+ Shs = Sig+calc_Shs(f,y);
rho = (tmass*AMU)/(V*NANO*NANO*NANO);
- fprintf(fplog,"Nmol = %d\nNatom = %d\ndt = %g ps\ntmass = %g amu\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,rho,Delta,f,beta,DoS0,dos2,dostot);
-
+ sigHS = pow(6*y*V/(M_PI*Natom),1.0/3.0);
+
+ fprintf(fplog,"System = %s\n",title);
+ fprintf(fplog,"Nmol = %d\n",Nmol);
+ fprintf(fplog,"Natom = %d\n",Natom);
+ fprintf(fplog,"dt = %g ps\n",dt);
+ fprintf(fplog,"tmass = %g amu\n",tmass);
+ fprintf(fplog,"V = %g nm^3\n",V);
+ fprintf(fplog,"rho = %g g/l\n",rho);
+ fprintf(fplog,"T = %g K\n",Temp);
+ fprintf(fplog,"beta = %g mol/kJ\n",beta);
+
+ fprintf(fplog,"\nDoS parameters\n");
+ fprintf(fplog,"Delta = %g\n",Delta);
+ fprintf(fplog,"fluidicity = %g\n",f);
+ fprintf(fplog,"hard sphere packing fraction = %g\n",y);
+ fprintf(fplog,"hard sphere compressibility = %g\n",z);
+ fprintf(fplog,"ideal gas entropy = %g\n",Sig);
+ fprintf(fplog,"hard sphere entropy = %g\n",Shs);
+ fprintf(fplog,"sigma_HS = %g nm\n",sigHS);
+ fprintf(fplog,"DoS0 = %g\n",DoS0);
+ fprintf(fplog,"Dos2 = %g\n",dos2);
+ fprintf(fplog,"DoSTot = %g\n",dostot);
+
/* Now compute solid (2) and diffusive (3) components */
fp = xvgropen(opt2fn("-dos",NFILE,fnm),"Density of states",
bRecip ? "E (cm\\S-1\\N)" : "\\f{12}n\\f{4} (1/ps)",
/* Finally analyze the results! */
wCdiff = 0.5;
- wSdiff = DoS0/(3*BOLTZ); /* Is this correct? */
+ wSdiff = Shs/(3*BOLTZ); /* Is this correct? */
for(j=0; (j<nframes/4); j++)
{
dos[DOS_CP][j] = (dos[DOS_DIFF][j]*wCdiff +
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff