return buf;
}
+static void remove_drift(int nset,int nbmin,int nbmax,real dt,enerdata_t *edat)
+{
+/* Remove the drift by performing a fit to y = ax+b.
+ Uses 5 iterations. */
+ int i,j,k;
+ double delta,d,sd,sd2;
+
+ edat->npoints = edat->nframes;
+ edat->nsteps = edat->nframes;
+
+ for(k=0; (k<5); k++)
+ {
+ for(i=0; (i<nset); i++)
+ {
+ delta = edat->s[i].slope*dt;
+
+ if (NULL != debug)
+ fprintf(debug,"slope for set %d is %g\n",i,edat->s[i].slope);
+
+ for(j=0; (j<edat->nframes); j++)
+ {
+ edat->s[i].ener[j] -= j*delta;
+ edat->s[i].es[j].sum = 0;
+ edat->s[i].es[j].sum2 = 0;
+ }
+ }
+ calc_averages(nset,edat,nbmin,nbmax);
+ }
+}
+
+static void calc_fluctuation_props(FILE *fp,
+ gmx_bool bDriftCorr,real dt,
+ int nset,int set[],int nmol,
+ char **leg,enerdata_t *edat,
+ int nbmin,int nbmax)
+{
+ int i,j;
+ double vvhh,vv,v,h,hh2,vv2,varv,hh,varh,tt,cv,cp,alpha,kappa,dcp,et,varet;
+ double NANO3;
+ enum { eVol, eEnth, eTemp, eEtot, eNR };
+ const char *my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" };
+ int ii[eNR];
+
+ NANO3 = NANO*NANO*NANO;
+ if (nmol < 2)
+ fprintf(fp,"\nWARNING: nmol = %d, this may not be what you want.\n",
+ nmol);
+ if (!bDriftCorr)
+ {
+ fprintf(fp,"\nYou may want to use the -driftcorr flag in order to correct\n"
+ "for spurious drift in the graphs. Note that this is not\n"
+ "a substitute for proper equilibration and sampling!\n");
+ }
+ else
+ {
+ remove_drift(nset,nbmin,nbmax,dt,edat);
+ }
+ for(i=0; (i<eNR); i++)
+ {
+ for(ii[i]=0; (ii[i]<nset &&
+ (gmx_strcasecmp(leg[ii[i]],my_ener[i]) != 0)); ii[i]++)
+ ;
+/* if (ii[i] < nset)
+ fprintf(fp,"Found %s data.\n",my_ener[i]);
+*/ }
+ /* Compute it all! */
+ vvhh = alpha = kappa = cp = dcp = cv = NOTSET;
+
+ /* Temperature */
+ tt = 0;
+ if (ii[eTemp] < nset)
+ {
+ tt = edat->s[ii[eTemp]].av;
+ }
+ /* Volume */
+ vv = varv = NOTSET;
+ if ((ii[eVol] < nset) && (ii[eTemp] < nset))
+ {
+ vv = edat->s[ii[eVol]].av*NANO3;
+ varv = dsqr(edat->s[ii[eVol]].rmsd*NANO3);
+ kappa = (varv/vv)/(BOLTZMANN*tt);
+ }
+ /* Enthalpy */
+ hh = varh = NOTSET;
+ if ((ii[eEnth] < nset) && (ii[eTemp] < nset))
+ {
+ hh = KILO*edat->s[ii[eEnth]].av/AVOGADRO;
+ varh = dsqr(KILO*edat->s[ii[eEnth]].rmsd/AVOGADRO);
+ cp = AVOGADRO*((varh/nmol)/(BOLTZMANN*tt*tt));
+ }
+ /* Total energy */
+ et = varet = NOTSET;
+ if ((ii[eEtot] < nset) && (hh == NOTSET))
+ {
+ /* Only compute cv in constant volume runs, which we can test
+ by checking whether the enthalpy was computed.
+ */
+ et = edat->s[ii[eEtot]].av;
+ varet = sqr(edat->s[ii[eEtot]].rmsd);
+ cv = KILO*((varet/nmol)/(BOLTZ*tt*tt));
+ }
+ /* Alpha, dcp */
+ if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
+ {
+ vvhh = 0;
+ for(j=0; (j<edat->nframes); j++)
+ {
+ v = edat->s[ii[eVol]].ener[j]*NANO3;
+ h = KILO*edat->s[ii[eEnth]].ener[j]/AVOGADRO;
+ vvhh += (v*h);
+ }
+ vvhh /= edat->nframes;
+ alpha = (vvhh-vv*hh)/(vv*BOLTZMANN*tt*tt);
+ dcp = (vv*AVOGADRO/nmol)*tt*sqr(alpha)/(kappa);
+ }
+
+ fprintf(fp,"\nTemperature dependent fluctuation properties at T = %g.\n",tt);
+ fprintf(fp,"\nHeat capacities obtained from fluctuations do *not* include\n");
+ fprintf(fp,"quantum corrections. If you want to get a more accurate estimate\n");
+ fprintf(fp,"please use the g_dos program.\n\n");
+
+ if (debug != NULL)
+ {
+ if (varv != NOTSET)
+ fprintf(fp,"varv = %10g (m^6)\n",varv*AVOGADRO/nmol);
+ if (vvhh != NOTSET)
+ fprintf(fp,"vvhh = %10g (m^3 J)\n",vvhh);
+ }
+ if (vv != NOTSET)
+ fprintf(fp,"Volume = %10g m^3/mol\n",
+ vv*AVOGADRO/nmol);
+ if (varh != NOTSET)
+ fprintf(fp,"Enthalpy = %10g kJ/mol\n",
+ hh*AVOGADRO/(KILO*nmol));
+ if (alpha != NOTSET)
+ fprintf(fp,"Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n",
+ alpha);
+ if (kappa != NOTSET)
+ {
+ fprintf(fp,"Isothermal Compressibility Kappa = %10g (J/m^3)\n",
+ kappa);
+ fprintf(fp,"Adiabatic bulk modulus = %10g (m^3/J)\n",
+ 1.0/kappa);
+ }
+ if (cp != NOTSET)
+ fprintf(fp,"Heat capacity at constant pressure Cp = %10g J/mol K\n",
+ cp);
+ if (cv != NOTSET)
+ fprintf(fp,"Heat capacity at constant volume Cv = %10g J/mol K\n",
+ cv);
+ if (dcp != NOTSET)
+ fprintf(fp,"Cp-Cv = %10g J/mol K\n",
+ dcp);
+ please_cite(fp,"Allen1987a");
+}
+
static void analyse_ener(gmx_bool bCorr,const char *corrfn,
- gmx_bool bFee,gmx_bool bSum,gmx_bool bFluct,
+ gmx_bool bFee,gmx_bool bSum,gmx_bool bFluct,gmx_bool bTempFluct,
gmx_bool bVisco,const char *visfn,int nmol,
gmx_large_int_t start_step,double start_t,
gmx_large_int_t step,double t,
/* Check out the printed manual for equations! */
double Dt,aver,stddev,errest,delta_t,totaldrift;
enerdata_t *esum=NULL;
- real xxx,integral,intBulk;
+ real xxx,integral,intBulk,Temp,Pres;
real sfrac,oldfrac,diffsum,diffav,fstep,pr_aver,pr_stddev,pr_errest;
double beta=0,expE,expEtot,*fee=NULL;
gmx_large_int_t nsteps;
int nexact,nnotexact;
double x1m,x1mk;
- real Temp=-1,Pres=-1,VarV=-1,VarT=-1,VarEtot=-1,AvEtot=0;
- int i,j,nout;
+ int i,j,k,nout;
real chi2;
char buf[256],eebuf[100];
fee[i] = log(expE/edat->nframes)/beta + aver/nmol;
}
if (strstr(leg[i],"empera") != NULL) {
- VarT = sqr(stddev);
Temp = aver;
} else if (strstr(leg[i],"olum") != NULL) {
- VarV = sqr(stddev);
Vaver= aver;
} else if (strstr(leg[i],"essure") != NULL) {
Pres = aver;
- } else if (strstr(leg[i],"otal") != NULL) {
- VarEtot = sqr(stddev);
- AvEtot = aver;
}
if (bIsEner[i]) {
pr_aver = aver/nmol-ezero;
else
fprintf(stdout,"\n");
}
+
/* Do correlation function */
if (edat->nframes > 1)
{
"The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
+ "Some fluctuation-dependent properties can be calculated provided",
+ "the correct energy terms are selected. The following properties",
+ "will be computed:[BR]",
+ "Property Energy terms needed[BR]",
+ "---------------------------------------------------[BR]",
+ "Heat capacity Cp (NPT sims): Enthalpy, Temp [BR]",
+ "Heat capacity Cv (NVT sims): Etot, Temp [BR]",
+ "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
+ "Isothermal compressibility: Vol, Temp [BR]",
+ "Adiabatic bulk modulus: Vol, Temp [BR]",
+ "---------------------------------------------------[BR]",
+ "You always need to set the number of molecules [TT]-nmol[tt].",
+ "The Cp/Cv computations do [BB]not[bb] include any corrections",
+ "for quantum effects. Use the [TT]g_dos[tt] program if you need that (and you do).[PAR]"
"When the [TT]-viol[tt] option is set, the time averaged",
"violations are plotted and the running time-averaged and",
"instantaneous sum of violations are recalculated. Additionally",
"[BB]Note[bb] that the energies must both be calculated from the same trajectory."
};
- static gmx_bool bSum=FALSE,bFee=FALSE,bPrAll=FALSE,bFluct=FALSE;
+ static gmx_bool bSum=FALSE,bFee=FALSE,bPrAll=FALSE,bFluct=FALSE,bDriftCorr=FALSE;
static gmx_bool bDp=FALSE,bMutot=FALSE,bOrinst=FALSE,bOvec=FALSE;
static int skip=0,nmol=1,nbmin=5,nbmax=5;
static real reftemp=300.0,ezero=0;
"Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
{ "-nmol", FALSE, etINT, {&nmol},
"Number of molecules in your sample: the energies are divided by this number" },
+ { "-driftcorr", FALSE, etBOOL, {&bDriftCorr},
+ "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."},
{ "-fluc", FALSE, etBOOL, {&bFluct},
"Calculate autocorrelation of energy fluctuations rather than energy itself" },
{ "-orinst", FALSE, etBOOL, {&bOrinst},
}
else
{
+ double dt = (frame[cur].t-start_t)/(edat.nframes-1);
analyse_ener(opt2bSet("-corr",NFILE,fnm),opt2fn("-corr",NFILE,fnm),
- bFee,bSum,bFluct,
+ bFee,bSum,bFluct,opt2parg_bSet("-nmol",npargs,ppa),
bVisco,opt2fn("-vis",NFILE,fnm),
- nmol,start_step,start_t,frame[cur].step,frame[cur].t,
+ nmol,
+ start_step,start_t,frame[cur].step,frame[cur].t,
time,reftemp,&edat,
nset,set,bIsEner,leg,enm,Vaver,ezero,nbmin,nbmax,
oenv);
+ calc_fluctuation_props(stdout,bDriftCorr,dt,nset,set,nmol,leg,&edat,
+ nbmin,nbmax);
}
if (opt2bSet("-f2",NFILE,fnm)) {
fec(opt2fn("-f2",NFILE,fnm), opt2fn("-ravg",NFILE,fnm),
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff