X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=blobdiff_plain;f=src%2Fgromacs%2Fgmxana%2Fgmx_dos.c;h=d4423adc47ca1e50d943d761b0b704f1f4d2d1c6;hb=6efa33ebea3ba1629d555e0e9e9646ad7e8cc9ad;hp=7592091e9b5253518470e14026979b8945646846;hpb=a751ebf72b76e72acc733b63faaa0215d660e7de;p=alexxy%2Fgromacs.git diff --git a/src/gromacs/gmxana/gmx_dos.c b/src/gromacs/gmxana/gmx_dos.c index 7592091e9b..d4423adc47 100644 --- a/src/gromacs/gmxana/gmx_dos.c +++ b/src/gromacs/gmxana/gmx_dos.c @@ -56,6 +56,7 @@ #include "gromacs/math/units.h" #include "gromacs/math/utilities.h" #include "gromacs/math/vec.h" +#include "gromacs/topology/index.h" #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/futil.h" #include "gromacs/utility/smalloc.h" @@ -65,6 +66,40 @@ enum { VACF, MVACF, DOS, DOS_SOLID, DOS_DIFF, DOS_CP, DOS_S, DOS_A, DOS_E, DOS_NR }; +static int calcMoleculesInIndexGroup(t_block *mols, int natoms, atom_id *index, int nindex) +{ + int i = 0; + int mol = 0; + int nMol = 0; + int j; + + while (i < nindex) + { + while (index[i] > mols->index[mol]) + { + mol++; + if (mol >= mols->nr) + { + gmx_fatal(FARGS, "Atom index out of range: %d", index[i]+1); + } + } + for (j = mols->index[mol]; j < mols->index[mol+1]; j++) + { + if (index[i] != j) + { + gmx_fatal(FARGS, "The index group does not consist of whole molecules"); + } + i++; + if (i == natoms) + { + gmx_fatal(FARGS, "Index contains atom numbers larger than the topology"); + } + } + nMol++; + } + return nMol; +} + static double FD(double Delta, double f) { return (2*pow(Delta, -4.5)*pow(f, 7.5) - @@ -214,47 +249,6 @@ static real wEsolid(real nu, real beta) } } -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); - if (output_env_get_print_xvgr_codes(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, toler); - fprintf(fp, "%10g %10g %10g %10g\n", Delta, f, f*y, y); - } - xvgrclose(fp); -} - -static void dump_w(output_env_t oenv, real beta) -{ - FILE *fp; - double nu; - const char *leg[] = { "wCv", "wS", "wA", "wE" }; - - 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 %10g %10g\n", beta*PLANCK*nu, - wCsolid(nu, beta), wSsolid(nu, beta), - wAsolid(nu, beta), wEsolid(nu, beta)); - } - xvgrclose(fp); -} - int gmx_dos(int argc, char *argv[]) { const char *desc[] = { @@ -264,6 +258,11 @@ int gmx_dos(int argc, char *argv[]) "all vibrations. For flexible systems that would be around a few fs", "between saving. Properties based on the DoS are printed on the", "standard output." + "Note that the density of states is calculated from the mass-weighted", + "autocorrelation, and by default only from the square of the real", + "component rather than absolute value. This means the shape can differ", + "substantially from the plain vibrational power spectrum you can", + "calculate with gmx velacc." }; const char *bugs[] = { "This program needs a lot of memory: total usage equals the number of atoms times 3 times number of frames times 4 (or 8 when run in double precision)." @@ -284,10 +283,16 @@ int gmx_dos(int argc, char *argv[]) gmx_fft_t fft; double cP, S, A, E, DiffCoeff, Delta, f, y, z, sigHS, Shs, Sig, DoS0, recip_fac; double wCdiff, wSdiff, wAdiff, wEdiff; - - static gmx_bool bVerbose = TRUE, bAbsolute = FALSE, bNormalize = FALSE; - static gmx_bool bRecip = FALSE, bDump = FALSE; + int grpNatoms; + atom_id *index; + char *grpname; + double invNormalize; + gmx_bool normalizeAutocorrelation; + + static gmx_bool bVerbose = TRUE, bAbsolute = FALSE, bNormalizeDos = FALSE; + static gmx_bool bRecip = FALSE; static real Temp = 298.15, toler = 1e-6; + t_pargs pa[] = { { "-v", FALSE, etBOOL, {&bVerbose}, "Be loud and noisy." }, @@ -295,14 +300,12 @@ int gmx_dos(int argc, char *argv[]) "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." }, + { "-normdos", FALSE, etBOOL, {&bNormalizeDos}, + "Normalize the DoS such that it adds up to 3N. This should usually not be necessary." }, { "-T", FALSE, etREAL, {&Temp}, "Temperature in the simulation" }, { "-toler", FALSE, etREAL, {&toler}, - "[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." } + "[HIDDEN]Tolerance when computing the fluidicity using bisection algorithm" } }; t_filenm fnm[] = { @@ -331,30 +334,26 @@ int gmx_dos(int argc, char *argv[]) } 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"); please_cite(fplog, "Caleman2011b"); - read_tps_conf(ftp2fn(efTPR, NFILE, fnm), title, &top, &ePBC, NULL, NULL, box, - TRUE); + read_tps_conf(ftp2fn(efTPR, NFILE, fnm), title, &top, &ePBC, NULL, NULL, box, TRUE); + + /* Handle index groups */ + get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &grpNatoms, &index, &grpname); + V = det(box); tmass = 0; - for (i = 0; (i < top.atoms.nr); i++) + for (i = 0; i < grpNatoms; i++) { - tmass += top.atoms.atom[i].m; + tmass += top.atoms.atom[index[i]].m; } - Natom = top.atoms.nr; - Nmol = top.mols.nr; + Natom = grpNatoms; + Nmol = calcMoleculesInIndexGroup(&top.mols, top.atoms.nr, index, grpNatoms); gnx = Natom*DIM; /* Correlation stuff */ @@ -390,9 +389,9 @@ int gmx_dos(int argc, char *argv[]) } for (i = 0; i < gnx; i += DIM) { - 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]; + c1[i+XX][nframes] = fr.v[index[i/DIM]][XX]; + c1[i+YY][nframes] = fr.v[index[i/DIM]][YY]; + c1[i+ZZ][nframes] = fr.v[index[i/DIM]][ZZ]; } t1 = fr.time; @@ -413,6 +412,21 @@ int gmx_dos(int argc, char *argv[]) printf("Going to do %d fourier transforms of length %d. Hang on.\n", gnx, nframes); } + /* Unfortunately the -normalize program option for the autocorrelation + * function calculation is added as a hack with a static variable in the + * autocorrelation.c source. That would work if we called the normal + * do_autocorr(), but this routine overrides that by directly calling + * the low-level functionality. That unfortunately leads to ignoring the + * default value for the option (which is to normalize). + * Since the absolute value seems to be important for the subsequent + * analysis below, we detect the value directly from the option, calculate + * the autocorrelation without normalization, and then apply the + * normalization just to the autocorrelation output + * (or not, if the user asked for a non-normalized autocorrelation). + */ + normalizeAutocorrelation = opt2parg_bool("-normalize", npargs, ppa); + + /* Note that we always disable normalization here, regardless of user settings */ low_do_autocorr(NULL, oenv, NULL, nframes, gnx, nframes, c1, dt, eacNormal, 0, FALSE, FALSE, FALSE, -1, -1, 0); snew(dos, DOS_NR); @@ -427,7 +441,7 @@ int gmx_dos(int argc, char *argv[]) } for (i = 0; (i < gnx); i += DIM) { - mi = top.atoms.atom[i/DIM].m; + mi = top.atoms.atom[index[i/DIM]].m; for (j = 0; (j < nframes/2); j++) { c1j = (c1[i+XX][j] + c1[i+YY][j] + c1[i+ZZ][j]); @@ -435,20 +449,28 @@ int gmx_dos(int argc, char *argv[]) dos[MVACF][j] += mi*c1j; } } - fp = xvgropen(opt2fn("-vacf", NFILE, fnm), "Velocity ACF", + + fp = xvgropen(opt2fn("-vacf", NFILE, fnm), "Velocity autocorrelation function", "Time (ps)", "C(t)", oenv); snew(tt, nframes/2); + + invNormalize = normalizeAutocorrelation ? 1.0/dos[VACF][0] : 1.0; + for (j = 0; (j < nframes/2); j++) { tt[j] = j*dt; - fprintf(fp, "%10g %10g\n", tt[j], dos[VACF][j]); + fprintf(fp, "%10g %10g\n", tt[j], dos[VACF][j] * invNormalize); } xvgrclose(fp); - fp = xvgropen(opt2fn("-mvacf", NFILE, fnm), "Mass-weighted velocity ACF", + + fp = xvgropen(opt2fn("-mvacf", NFILE, fnm), "Mass-weighted velocity autocorrelation function", "Time (ps)", "C(t)", oenv); + + invNormalize = normalizeAutocorrelation ? 1.0/dos[VACF][0] : 1.0; + for (j = 0; (j < nframes/2); j++) { - fprintf(fp, "%10g %10g\n", tt[j], dos[MVACF][j]); + fprintf(fp, "%10g %10g\n", tt[j], dos[MVACF][j] * invNormalize); } xvgrclose(fp); @@ -483,7 +505,7 @@ int gmx_dos(int argc, char *argv[]) } /* Normalize it */ dostot = evaluate_integral(nframes/4, nu, dos[DOS], NULL, nframes/4, &stddev); - if (bNormalize) + if (bNormalizeDos) { for (j = 0; (j < nframes/4); j++) { @@ -571,18 +593,6 @@ int gmx_dos(int argc, char *argv[]) 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(nframes/4,nu,dos[DOS_S],NULL, - nframes/4,&stddev); - fprintf(fplog,"Entropy %g J/mol K\n",1000*S/Nmol); - A = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_A],NULL, - nframes/4,&stddev); - fprintf(fplog,"Helmholtz energy %g kJ/mol\n",A/Nmol); - E = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_E],NULL, - nframes/4,&stddev); - fprintf(fplog,"Internal energy %g kJ/mol\n",E/Nmol); - */ fprintf(fplog, "\nArrivederci!\n"); gmx_fio_fclose(fplog);