2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2011,2012,2013,2014, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
44 #include "gromacs/fileio/confio.h"
46 #include "gromacs/utility/fatalerror.h"
47 #include "gromacs/utility/futil.h"
50 #include "gromacs/math/utilities.h"
51 #include "gromacs/math/units.h"
52 #include "gromacs/utility/smalloc.h"
53 #include "gromacs/commandline/pargs.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/fileio/xvgr.h"
61 #include "gromacs/fft/fft.h"
62 #include "gromacs/fileio/trxio.h"
65 VACF, MVACF, DOS, DOS_SOLID, DOS_DIFF, DOS_CP, DOS_S, DOS_A, DOS_E, DOS_NR
68 static double FD(double Delta, double f)
70 return (2*pow(Delta, -4.5)*pow(f, 7.5) -
71 6*pow(Delta, -3)*pow(f, 5) -
72 pow(Delta, -1.5)*pow(f, 3.5) +
73 6*pow(Delta, -1.5)*pow(f, 2.5) +
77 static double YYY(double f, double y)
79 return (2*pow(y*f, 3) - sqr(f)*y*(1+6*y) +
83 static double calc_compress(double y)
89 return ((1+y+sqr(y)-pow(y, 3))/(pow(1-y, 3)));
92 static double bisector(double Delta, double tol,
93 double ff0, double ff1,
94 double ff(double, double))
96 double fd0, fd, fd1, f, f0, f1;
103 fprintf(stderr, "Unrealistic tolerance %g for bisector. Setting it to %g\n", tol, tolmin);
126 while ((f1-f0) > tol);
131 static double calc_fluidicity(double Delta, double tol)
133 return bisector(Delta, tol, 0, 1, FD);
136 static double calc_y(double f, double Delta, double toler)
140 y1 = pow(f/Delta, 1.5);
141 y2 = bisector(f, toler, 0, 10000, YYY);
142 if (fabs((y1-y2)/(y1+y2)) > 100*toler)
144 fprintf(stderr, "Inconsistency computing y: y1 = %f, y2 = %f, using y1.\n",
151 static double calc_Shs(double f, double y)
155 return BOLTZ*(log(calc_compress(fy)) + fy*(3*fy-4)/sqr(1-fy));
158 static real wCsolid(real nu, real beta)
160 real bhn = beta*PLANCK*nu;
171 return sqr(bhn)*ebn/koko;
175 static real wSsolid(real nu, real beta)
177 real bhn = beta*PLANCK*nu;
185 return bhn/(exp(bhn)-1) - log(1-exp(-bhn));
189 static real wAsolid(real nu, real beta)
191 real bhn = beta*PLANCK*nu;
199 return log((1-exp(-bhn))/(exp(-bhn/2))) - log(bhn);
203 static real wEsolid(real nu, real beta)
205 real bhn = beta*PLANCK*nu;
213 return bhn/2 + bhn/(exp(bhn)-1)-1;
217 static void dump_fy(output_env_t oenv, real toler)
220 double Delta, f, y, DD;
221 const char *leg[] = { "f", "fy", "y" };
223 DD = pow(10.0, 0.125);
224 fp = xvgropen("fy.xvg", "Fig. 2, Lin2003a", "Delta", "y or fy", oenv);
225 xvgr_legend(fp, asize(leg), leg, oenv);
226 if (output_env_get_print_xvgr_codes(oenv))
228 fprintf(fp, "@ world 1e-05, 0, 1000, 1\n");
229 fprintf(fp, "@ xaxes scale Logarithmic\n");
231 for (Delta = 1e-5; (Delta <= 1000); Delta *= DD)
233 f = calc_fluidicity(Delta, toler);
234 y = calc_y(f, Delta, toler);
235 fprintf(fp, "%10g %10g %10g %10g\n", Delta, f, f*y, y);
240 static void dump_w(output_env_t oenv, real beta)
244 const char *leg[] = { "wCv", "wS", "wA", "wE" };
246 fp = xvgropen("w.xvg", "Fig. 1, Berens1983a", "\\f{12}b\\f{4}h\\f{12}n",
248 xvgr_legend(fp, asize(leg), leg, oenv);
249 for (nu = 1; (nu < 100); nu += 0.05)
251 fprintf(fp, "%10g %10g %10g %10g %10g\n", beta*PLANCK*nu,
252 wCsolid(nu, beta), wSsolid(nu, beta),
253 wAsolid(nu, beta), wEsolid(nu, beta));
258 int gmx_dos(int argc, char *argv[])
260 const char *desc[] = {
261 "[THISMODULE] computes the Density of States from a simulations.",
262 "In order for this to be meaningful the velocities must be saved",
263 "in the trajecotry with sufficiently high frequency such as to cover",
264 "all vibrations. For flexible systems that would be around a few fs",
265 "between saving. Properties based on the DoS are printed on the",
268 const char *bugs[] = {
269 "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)."
280 int nV, nframes, n_alloc, i, j, k, l, fftcode, Nmol, Natom;
281 double rho, dt, V2sum, Vsum, V, tmass, dostot, dos2, dosabs;
282 real **c1, **dos, mi, beta, bfac, *nu, *tt, stddev, c1j;
285 double cP, S, A, E, DiffCoeff, Delta, f, y, z, sigHS, Shs, Sig, DoS0, recip_fac;
286 double wCdiff, wSdiff, wAdiff, wEdiff;
288 static gmx_bool bVerbose = TRUE, bAbsolute = FALSE, bNormalize = FALSE;
289 static gmx_bool bRecip = FALSE, bDump = FALSE;
290 static real Temp = 298.15, toler = 1e-6;
292 { "-v", FALSE, etBOOL, {&bVerbose},
293 "Be loud and noisy." },
294 { "-recip", FALSE, etBOOL, {&bRecip},
295 "Use cm^-1 on X-axis instead of 1/ps for DoS plots." },
296 { "-abs", FALSE, etBOOL, {&bAbsolute},
297 "Use the absolute value of the Fourier transform of the VACF as the Density of States. Default is to use the real component only" },
298 { "-normdos", FALSE, etBOOL, {&bNormalize},
299 "Normalize the DoS such that it adds up to 3N. This is a hack that should not be necessary." },
300 { "-T", FALSE, etREAL, {&Temp},
301 "Temperature in the simulation" },
302 { "-toler", FALSE, etREAL, {&toler},
303 "[HIDDEN]Tolerance when computing the fluidicity using bisection algorithm" },
304 { "-dump", FALSE, etBOOL, {&bDump},
305 "[HIDDEN]Dump the y/fy plot corresponding to Fig. 2 inLin2003a and the and the weighting functions corresponding to Fig. 1 in Berens1983a." }
309 { efTRN, "-f", NULL, ffREAD },
310 { efTPX, "-s", NULL, ffREAD },
311 { efNDX, NULL, NULL, ffOPTRD },
312 { efXVG, "-vacf", "vacf", ffWRITE },
313 { efXVG, "-mvacf", "mvacf", ffWRITE },
314 { efXVG, "-dos", "dos", ffWRITE },
315 { efLOG, "-g", "dos", ffWRITE },
317 #define NFILE asize(fnm)
320 const char *DoSlegend[] = {
321 "DoS(v)", "DoS(v)[Solid]", "DoS(v)[Diff]"
325 ppa = add_acf_pargs(&npargs, pa);
326 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
327 NFILE, fnm, npargs, ppa, asize(desc), desc,
328 asize(bugs), bugs, &oenv))
333 beta = 1/(Temp*BOLTZ);
336 printf("Dumping reference figures. Thanks for your patience.\n");
337 dump_fy(oenv, toler);
342 fplog = gmx_fio_fopen(ftp2fn(efLOG, NFILE, fnm), "w");
343 fprintf(fplog, "Doing density of states analysis based on trajectory.\n");
344 please_cite(fplog, "Pascal2011a");
345 please_cite(fplog, "Caleman2011b");
347 read_tps_conf(ftp2fn(efTPX, NFILE, fnm), title, &top, &ePBC, NULL, NULL, box,
351 for (i = 0; (i < top.atoms.nr); i++)
353 tmass += top.atoms.atom[i].m;
356 Natom = top.atoms.nr;
360 /* Correlation stuff */
362 for (i = 0; (i < gnx); i++)
367 read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE, fnm), &fr, TRX_NEED_V);
383 if (nframes >= n_alloc)
386 for (i = 0; i < gnx; i++)
388 srenew(c1[i], n_alloc);
391 for (i = 0; i < gnx; i += DIM)
393 c1[i+XX][nframes] = fr.v[i/DIM][XX];
394 c1[i+YY][nframes] = fr.v[i/DIM][YY];
395 c1[i+ZZ][nframes] = fr.v[i/DIM][ZZ];
402 while (read_next_frame(oenv, status, &fr));
406 dt = (t1-t0)/(nframes-1);
413 printf("Going to do %d fourier transforms of length %d. Hang on.\n",
416 low_do_autocorr(NULL, oenv, NULL, nframes, gnx, nframes, c1, dt, eacNormal, 0, FALSE,
417 FALSE, FALSE, -1, -1, 0);
419 for (j = 0; (j < DOS_NR); j++)
421 snew(dos[j], nframes+4);
426 printf("Going to merge the ACFs into the mass-weighted and plain ACF\n");
428 for (i = 0; (i < gnx); i += DIM)
430 mi = top.atoms.atom[i/DIM].m;
431 for (j = 0; (j < nframes/2); j++)
433 c1j = (c1[i+XX][j] + c1[i+YY][j] + c1[i+ZZ][j]);
434 dos[VACF][j] += c1j/Natom;
435 dos[MVACF][j] += mi*c1j;
438 fp = xvgropen(opt2fn("-vacf", NFILE, fnm), "Velocity ACF",
439 "Time (ps)", "C(t)", oenv);
441 for (j = 0; (j < nframes/2); j++)
444 fprintf(fp, "%10g %10g\n", tt[j], dos[VACF][j]);
447 fp = xvgropen(opt2fn("-mvacf", NFILE, fnm), "Mass-weighted velocity ACF",
448 "Time (ps)", "C(t)", oenv);
449 for (j = 0; (j < nframes/2); j++)
451 fprintf(fp, "%10g %10g\n", tt[j], dos[MVACF][j]);
455 if ((fftcode = gmx_fft_init_1d_real(&fft, nframes/2,
456 GMX_FFT_FLAG_NONE)) != 0)
458 gmx_fatal(FARGS, "gmx_fft_init_1d_real returned %d", fftcode);
460 if ((fftcode = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX,
461 (void *)dos[MVACF], (void *)dos[DOS])) != 0)
463 gmx_fatal(FARGS, "gmx_fft_1d_real returned %d", fftcode);
466 /* First compute the DoS */
467 /* Magic factor of 8 included now. */
471 for (j = 0; (j < nframes/4); j++)
474 dos2 += sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]);
477 dos[DOS][j] = bfac*sqrt(sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]));
481 dos[DOS][j] = bfac*dos[DOS][2*j];
485 dostot = evaluate_integral(nframes/4, nu, dos[DOS], NULL, nframes/4, &stddev);
488 for (j = 0; (j < nframes/4); j++)
490 dos[DOS][j] *= 3*Natom/dostot;
497 /* Note this eqn. is incorrect in Pascal2011a! */
498 Delta = ((2*DoS0/(9*Natom))*sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*
499 pow((Natom/V), 1.0/3.0)*pow(6/M_PI, 2.0/3.0));
500 f = calc_fluidicity(Delta, toler);
501 y = calc_y(f, Delta, toler);
502 z = calc_compress(y);
503 Sig = BOLTZ*(5.0/2.0+log(2*M_PI*BOLTZ*Temp/(sqr(PLANCK))*V/(f*Natom)));
504 Shs = Sig+calc_Shs(f, y);
505 rho = (tmass*AMU)/(V*NANO*NANO*NANO);
506 sigHS = pow(6*y*V/(M_PI*Natom), 1.0/3.0);
508 fprintf(fplog, "System = \"%s\"\n", title);
509 fprintf(fplog, "Nmol = %d\n", Nmol);
510 fprintf(fplog, "Natom = %d\n", Natom);
511 fprintf(fplog, "dt = %g ps\n", dt);
512 fprintf(fplog, "tmass = %g amu\n", tmass);
513 fprintf(fplog, "V = %g nm^3\n", V);
514 fprintf(fplog, "rho = %g g/l\n", rho);
515 fprintf(fplog, "T = %g K\n", Temp);
516 fprintf(fplog, "beta = %g mol/kJ\n", beta);
518 fprintf(fplog, "\nDoS parameters\n");
519 fprintf(fplog, "Delta = %g\n", Delta);
520 fprintf(fplog, "fluidicity = %g\n", f);
521 fprintf(fplog, "hard sphere packing fraction = %g\n", y);
522 fprintf(fplog, "hard sphere compressibility = %g\n", z);
523 fprintf(fplog, "ideal gas entropy = %g\n", Sig);
524 fprintf(fplog, "hard sphere entropy = %g\n", Shs);
525 fprintf(fplog, "sigma_HS = %g nm\n", sigHS);
526 fprintf(fplog, "DoS0 = %g\n", DoS0);
527 fprintf(fplog, "Dos2 = %g\n", dos2);
528 fprintf(fplog, "DoSTot = %g\n", dostot);
530 /* Now compute solid (2) and diffusive (3) components */
531 fp = xvgropen(opt2fn("-dos", NFILE, fnm), "Density of states",
532 bRecip ? "E (cm\\S-1\\N)" : "\\f{12}n\\f{4} (1/ps)",
533 "\\f{4}S(\\f{12}n\\f{4})", oenv);
534 xvgr_legend(fp, asize(DoSlegend), DoSlegend, oenv);
535 recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;
536 for (j = 0; (j < nframes/4); j++)
538 dos[DOS_DIFF][j] = DoS0/(1+sqr(DoS0*M_PI*nu[j]/(6*f*Natom)));
539 dos[DOS_SOLID][j] = dos[DOS][j]-dos[DOS_DIFF][j];
540 fprintf(fp, "%10g %10g %10g %10g\n",
542 dos[DOS][j]/recip_fac,
543 dos[DOS_SOLID][j]/recip_fac,
544 dos[DOS_DIFF][j]/recip_fac);
548 /* Finally analyze the results! */
550 wSdiff = Shs/(3*BOLTZ); /* Is this correct? */
552 wAdiff = wEdiff-wSdiff;
553 for (j = 0; (j < nframes/4); j++)
555 dos[DOS_CP][j] = (dos[DOS_DIFF][j]*wCdiff +
556 dos[DOS_SOLID][j]*wCsolid(nu[j], beta));
557 dos[DOS_S][j] = (dos[DOS_DIFF][j]*wSdiff +
558 dos[DOS_SOLID][j]*wSsolid(nu[j], beta));
559 dos[DOS_A][j] = (dos[DOS_DIFF][j]*wAdiff +
560 dos[DOS_SOLID][j]*wAsolid(nu[j], beta));
561 dos[DOS_E][j] = (dos[DOS_DIFF][j]*wEdiff +
562 dos[DOS_SOLID][j]*wEsolid(nu[j], beta));
564 DiffCoeff = evaluate_integral(nframes/2, tt, dos[VACF], NULL, nframes/2, &stddev);
565 DiffCoeff = 1000*DiffCoeff/3.0;
566 fprintf(fplog, "Diffusion coefficient from VACF %g 10^-5 cm^2/s\n",
568 fprintf(fplog, "Diffusion coefficient from DoS %g 10^-5 cm^2/s\n",
569 1000*DoS0/(12*tmass*beta));
571 cP = BOLTZ * evaluate_integral(nframes/4, nu, dos[DOS_CP], NULL,
573 fprintf(fplog, "Heat capacity %g J/mol K\n", 1000*cP/Nmol);
576 S = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_S],NULL,
578 fprintf(fplog,"Entropy %g J/mol K\n",1000*S/Nmol);
579 A = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_A],NULL,
581 fprintf(fplog,"Helmholtz energy %g kJ/mol\n",A/Nmol);
582 E = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_E],NULL,
584 fprintf(fplog,"Internal energy %g kJ/mol\n",E/Nmol);
586 fprintf(fplog, "\nArrivederci!\n");
587 gmx_fio_fclose(fplog);
589 do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");