1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
43 #include "gmx_fatal.h"
61 #include "gromacs/fft/fft.h"
64 VACF, MVACF, DOS, DOS_SOLID, DOS_DIFF, DOS_CP, DOS_S, DOS_A, DOS_E, DOS_NR
67 static double FD(double Delta, double f)
69 return (2*pow(Delta, -4.5)*pow(f, 7.5) -
70 6*pow(Delta, -3)*pow(f, 5) -
71 pow(Delta, -1.5)*pow(f, 3.5) +
72 6*pow(Delta, -1.5)*pow(f, 2.5) +
76 static double YYY(double f, double y)
78 return (2*pow(y*f, 3) - sqr(f)*y*(1+6*y) +
82 static double calc_compress(double y)
88 return ((1+y+sqr(y)-pow(y, 3))/(pow(1-y, 3)));
91 static double bisector(double Delta, double tol,
92 double ff0, double ff1,
93 double ff(double, double))
95 double fd0, fd, fd1, f, f0, f1;
102 fprintf(stderr, "Unrealistic tolerance %g for bisector. Setting it to %g\n", tol, tolmin);
125 while ((f1-f0) > tol);
130 static double calc_fluidicity(double Delta, double tol)
132 return bisector(Delta, tol, 0, 1, FD);
135 static double calc_y(double f, double Delta, double toler)
139 y1 = pow(f/Delta, 1.5);
140 y2 = bisector(f, toler, 0, 10000, YYY);
141 if (fabs((y1-y2)/(y1+y2)) > 100*toler)
143 fprintf(stderr, "Inconsistency computing y: y1 = %f, y2 = %f, using y1.\n",
150 static double calc_Shs(double f, double y)
154 return BOLTZ*(log(calc_compress(fy)) + fy*(3*fy-4)/sqr(1-fy));
157 static real wCsolid(real nu, real beta)
159 real bhn = beta*PLANCK*nu;
170 return sqr(bhn)*ebn/koko;
174 static real wSsolid(real nu, real beta)
176 real bhn = beta*PLANCK*nu;
184 return bhn/(exp(bhn)-1) - log(1-exp(-bhn));
188 static real wAsolid(real nu, real beta)
190 real bhn = beta*PLANCK*nu;
198 return log((1-exp(-bhn))/(exp(-bhn/2))) - log(bhn);
202 static real wEsolid(real nu, real beta)
204 real bhn = beta*PLANCK*nu;
212 return bhn/2 + bhn/(exp(bhn)-1)-1;
216 static void dump_fy(output_env_t oenv, real toler)
219 double Delta, f, y, DD;
220 const char *leg[] = { "f", "fy", "y" };
222 DD = pow(10.0, 0.125);
223 fp = xvgropen("fy.xvg", "Fig. 2, Lin2003a", "Delta", "y or fy", oenv);
224 xvgr_legend(fp, asize(leg), leg, oenv);
225 fprintf(fp, "@ world 1e-05, 0, 1000, 1\n");
226 fprintf(fp, "@ xaxes scale Logarithmic\n");
227 for (Delta = 1e-5; (Delta <= 1000); Delta *= DD)
229 f = calc_fluidicity(Delta, toler);
230 y = calc_y(f, Delta, toler);
231 fprintf(fp, "%10g %10g %10g %10g\n", Delta, f, f*y, y);
236 static void dump_w(output_env_t oenv, real beta)
240 const char *leg[] = { "wCv", "wS", "wA", "wE" };
242 fp = xvgropen("w.xvg", "Fig. 1, Berens1983a", "\\f{12}b\\f{4}h\\f{12}n",
244 xvgr_legend(fp, asize(leg), leg, oenv);
245 for (nu = 1; (nu < 100); nu += 0.05)
247 fprintf(fp, "%10g %10g %10g %10g %10g\n", beta*PLANCK*nu,
248 wCsolid(nu, beta), wSsolid(nu, beta),
249 wAsolid(nu, beta), wEsolid(nu, beta));
254 int gmx_dos(int argc, char *argv[])
256 const char *desc[] = {
257 "[TT]g_dos[tt] computes the Density of States from a simulations.",
258 "In order for this to be meaningful the velocities must be saved",
259 "in the trajecotry with sufficiently high frequency such as to cover",
260 "all vibrations. For flexible systems that would be around a few fs",
261 "between saving. Properties based on the DoS are printed on the",
264 const char *bugs[] = {
265 "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)."
276 int nV, nframes, n_alloc, i, j, k, l, fftcode, Nmol, Natom;
277 double rho, dt, V2sum, Vsum, V, tmass, dostot, dos2, dosabs;
278 real **c1, **dos, mi, beta, bfac, *nu, *tt, stddev, c1j;
281 double cP, S, A, E, DiffCoeff, Delta, f, y, z, sigHS, Shs, Sig, DoS0, recip_fac;
282 double wCdiff, wSdiff, wAdiff, wEdiff;
284 static gmx_bool bVerbose = TRUE, bAbsolute = FALSE, bNormalize = FALSE;
285 static gmx_bool bRecip = FALSE, bDump = FALSE;
286 static real Temp = 298.15, toler = 1e-6;
288 { "-v", FALSE, etBOOL, {&bVerbose},
289 "Be loud and noisy." },
290 { "-recip", FALSE, etBOOL, {&bRecip},
291 "Use cm^-1 on X-axis instead of 1/ps for DoS plots." },
292 { "-abs", FALSE, etBOOL, {&bAbsolute},
293 "Use the absolute value of the Fourier transform of the VACF as the Density of States. Default is to use the real component only" },
294 { "-normdos", FALSE, etBOOL, {&bNormalize},
295 "Normalize the DoS such that it adds up to 3N. This is a hack that should not be necessary." },
296 { "-T", FALSE, etREAL, {&Temp},
297 "Temperature in the simulation" },
298 { "-toler", FALSE, etREAL, {&toler},
299 "[HIDDEN]Tolerance when computing the fluidicity using bisection algorithm" },
300 { "-dump", FALSE, etBOOL, {&bDump},
301 "[HIDDEN]Dump the y/fy plot corresponding to Fig. 2 inLin2003a and the and the weighting functions corresponding to Fig. 1 in Berens1983a." }
305 { efTRN, "-f", NULL, ffREAD },
306 { efTPX, "-s", NULL, ffREAD },
307 { efNDX, NULL, NULL, ffOPTRD },
308 { efXVG, "-vacf", "vacf", ffWRITE },
309 { efXVG, "-mvacf", "mvacf", ffWRITE },
310 { efXVG, "-dos", "dos", ffWRITE },
311 { efLOG, "-g", "dos", ffWRITE },
313 #define NFILE asize(fnm)
316 const char *DoSlegend[] = {
317 "DoS(v)", "DoS(v)[Solid]", "DoS(v)[Diff]"
321 ppa = add_acf_pargs(&npargs, pa);
322 parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
323 NFILE, fnm, npargs, ppa, asize(desc), desc,
324 asize(bugs), bugs, &oenv);
326 beta = 1/(Temp*BOLTZ);
329 printf("Dumping reference figures. Thanks for your patience.\n");
330 dump_fy(oenv, toler);
335 fplog = gmx_fio_fopen(ftp2fn(efLOG, NFILE, fnm), "w");
336 fprintf(fplog, "Doing density of states analysis based on trajectory.\n");
337 please_cite(fplog, "Pascal2011a");
338 please_cite(fplog, "Caleman2011b");
340 read_tps_conf(ftp2fn(efTPX, NFILE, fnm), title, &top, &ePBC, NULL, NULL, box,
344 for (i = 0; (i < top.atoms.nr); i++)
346 tmass += top.atoms.atom[i].m;
349 Natom = top.atoms.nr;
353 /* Correlation stuff */
355 for (i = 0; (i < gnx); i++)
360 read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE, fnm), &fr, TRX_NEED_V);
376 if (nframes >= n_alloc)
379 for (i = 0; i < gnx; i++)
381 srenew(c1[i], n_alloc);
384 for (i = 0; i < gnx; i += DIM)
386 c1[i+XX][nframes] = fr.v[i/DIM][XX];
387 c1[i+YY][nframes] = fr.v[i/DIM][YY];
388 c1[i+ZZ][nframes] = fr.v[i/DIM][ZZ];
395 while (read_next_frame(oenv, status, &fr));
399 dt = (t1-t0)/(nframes-1);
406 printf("Going to do %d fourier transforms of length %d. Hang on.\n",
409 low_do_autocorr(NULL, oenv, NULL, nframes, gnx, nframes, c1, dt, eacNormal, 0, FALSE,
410 FALSE, FALSE, -1, -1, 0);
412 for (j = 0; (j < DOS_NR); j++)
414 snew(dos[j], nframes+4);
419 printf("Going to merge the ACFs into the mass-weighted and plain ACF\n");
421 for (i = 0; (i < gnx); i += DIM)
423 mi = top.atoms.atom[i/DIM].m;
424 for (j = 0; (j < nframes/2); j++)
426 c1j = (c1[i+XX][j] + c1[i+YY][j] + c1[i+ZZ][j]);
427 dos[VACF][j] += c1j/Natom;
428 dos[MVACF][j] += mi*c1j;
431 fp = xvgropen(opt2fn("-vacf", NFILE, fnm), "Velocity ACF",
432 "Time (ps)", "C(t)", oenv);
434 for (j = 0; (j < nframes/2); j++)
437 fprintf(fp, "%10g %10g\n", tt[j], dos[VACF][j]);
440 fp = xvgropen(opt2fn("-mvacf", NFILE, fnm), "Mass-weighted velocity ACF",
441 "Time (ps)", "C(t)", oenv);
442 for (j = 0; (j < nframes/2); j++)
444 fprintf(fp, "%10g %10g\n", tt[j], dos[MVACF][j]);
448 if ((fftcode = gmx_fft_init_1d_real(&fft, nframes/2,
449 GMX_FFT_FLAG_NONE)) != 0)
451 gmx_fatal(FARGS, "gmx_fft_init_1d_real returned %d", fftcode);
453 if ((fftcode = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX,
454 (void *)dos[MVACF], (void *)dos[DOS])) != 0)
456 gmx_fatal(FARGS, "gmx_fft_1d_real returned %d", fftcode);
459 /* First compute the DoS */
460 /* Magic factor of 8 included now. */
464 for (j = 0; (j < nframes/4); j++)
467 dos2 += sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]);
470 dos[DOS][j] = bfac*sqrt(sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]));
474 dos[DOS][j] = bfac*dos[DOS][2*j];
478 dostot = evaluate_integral(nframes/4, nu, dos[DOS], NULL, nframes/4, &stddev);
481 for (j = 0; (j < nframes/4); j++)
483 dos[DOS][j] *= 3*Natom/dostot;
490 /* Note this eqn. is incorrect in Pascal2011a! */
491 Delta = ((2*DoS0/(9*Natom))*sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*
492 pow((Natom/V), 1.0/3.0)*pow(6/M_PI, 2.0/3.0));
493 f = calc_fluidicity(Delta, toler);
494 y = calc_y(f, Delta, toler);
495 z = calc_compress(y);
496 Sig = BOLTZ*(5.0/2.0+log(2*M_PI*BOLTZ*Temp/(sqr(PLANCK))*V/(f*Natom)));
497 Shs = Sig+calc_Shs(f, y);
498 rho = (tmass*AMU)/(V*NANO*NANO*NANO);
499 sigHS = pow(6*y*V/(M_PI*Natom), 1.0/3.0);
501 fprintf(fplog, "System = \"%s\"\n", title);
502 fprintf(fplog, "Nmol = %d\n", Nmol);
503 fprintf(fplog, "Natom = %d\n", Natom);
504 fprintf(fplog, "dt = %g ps\n", dt);
505 fprintf(fplog, "tmass = %g amu\n", tmass);
506 fprintf(fplog, "V = %g nm^3\n", V);
507 fprintf(fplog, "rho = %g g/l\n", rho);
508 fprintf(fplog, "T = %g K\n", Temp);
509 fprintf(fplog, "beta = %g mol/kJ\n", beta);
511 fprintf(fplog, "\nDoS parameters\n");
512 fprintf(fplog, "Delta = %g\n", Delta);
513 fprintf(fplog, "fluidicity = %g\n", f);
514 fprintf(fplog, "hard sphere packing fraction = %g\n", y);
515 fprintf(fplog, "hard sphere compressibility = %g\n", z);
516 fprintf(fplog, "ideal gas entropy = %g\n", Sig);
517 fprintf(fplog, "hard sphere entropy = %g\n", Shs);
518 fprintf(fplog, "sigma_HS = %g nm\n", sigHS);
519 fprintf(fplog, "DoS0 = %g\n", DoS0);
520 fprintf(fplog, "Dos2 = %g\n", dos2);
521 fprintf(fplog, "DoSTot = %g\n", dostot);
523 /* Now compute solid (2) and diffusive (3) components */
524 fp = xvgropen(opt2fn("-dos", NFILE, fnm), "Density of states",
525 bRecip ? "E (cm\\S-1\\N)" : "\\f{12}n\\f{4} (1/ps)",
526 "\\f{4}S(\\f{12}n\\f{4})", oenv);
527 xvgr_legend(fp, asize(DoSlegend), DoSlegend, oenv);
528 recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;
529 for (j = 0; (j < nframes/4); j++)
531 dos[DOS_DIFF][j] = DoS0/(1+sqr(DoS0*M_PI*nu[j]/(6*f*Natom)));
532 dos[DOS_SOLID][j] = dos[DOS][j]-dos[DOS_DIFF][j];
533 fprintf(fp, "%10g %10g %10g %10g\n",
535 dos[DOS][j]/recip_fac,
536 dos[DOS_SOLID][j]/recip_fac,
537 dos[DOS_DIFF][j]/recip_fac);
541 /* Finally analyze the results! */
543 wSdiff = Shs/(3*BOLTZ); /* Is this correct? */
545 wAdiff = wEdiff-wSdiff;
546 for (j = 0; (j < nframes/4); j++)
548 dos[DOS_CP][j] = (dos[DOS_DIFF][j]*wCdiff +
549 dos[DOS_SOLID][j]*wCsolid(nu[j], beta));
550 dos[DOS_S][j] = (dos[DOS_DIFF][j]*wSdiff +
551 dos[DOS_SOLID][j]*wSsolid(nu[j], beta));
552 dos[DOS_A][j] = (dos[DOS_DIFF][j]*wAdiff +
553 dos[DOS_SOLID][j]*wAsolid(nu[j], beta));
554 dos[DOS_E][j] = (dos[DOS_DIFF][j]*wEdiff +
555 dos[DOS_SOLID][j]*wEsolid(nu[j], beta));
557 DiffCoeff = evaluate_integral(nframes/2, tt, dos[VACF], NULL, nframes/2, &stddev);
558 DiffCoeff = 1000*DiffCoeff/3.0;
559 fprintf(fplog, "Diffusion coefficient from VACF %g 10^-5 cm^2/s\n",
561 fprintf(fplog, "Diffusion coefficient from DoS %g 10^-5 cm^2/s\n",
562 1000*DoS0/(12*tmass*beta));
564 cP = BOLTZ * evaluate_integral(nframes/4, nu, dos[DOS_CP], NULL,
566 fprintf(fplog, "Heat capacity %g J/mol K\n", 1000*cP/Nmol);
569 S = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_S],NULL,
571 fprintf(fplog,"Entropy %g J/mol K\n",1000*S/Nmol);
572 A = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_A],NULL,
574 fprintf(fplog,"Helmholtz energy %g kJ/mol\n",A/Nmol);
575 E = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_E],NULL,
577 fprintf(fplog,"Internal energy %g kJ/mol\n",E/Nmol);
579 fprintf(fplog, "\nArrivederci!\n");
580 gmx_fio_fclose(fplog);
582 do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");