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 if (!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))
329 beta = 1/(Temp*BOLTZ);
332 printf("Dumping reference figures. Thanks for your patience.\n");
333 dump_fy(oenv, toler);
338 fplog = gmx_fio_fopen(ftp2fn(efLOG, NFILE, fnm), "w");
339 fprintf(fplog, "Doing density of states analysis based on trajectory.\n");
340 please_cite(fplog, "Pascal2011a");
341 please_cite(fplog, "Caleman2011b");
343 read_tps_conf(ftp2fn(efTPX, NFILE, fnm), title, &top, &ePBC, NULL, NULL, box,
347 for (i = 0; (i < top.atoms.nr); i++)
349 tmass += top.atoms.atom[i].m;
352 Natom = top.atoms.nr;
356 /* Correlation stuff */
358 for (i = 0; (i < gnx); i++)
363 read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE, fnm), &fr, TRX_NEED_V);
379 if (nframes >= n_alloc)
382 for (i = 0; i < gnx; i++)
384 srenew(c1[i], n_alloc);
387 for (i = 0; i < gnx; i += DIM)
389 c1[i+XX][nframes] = fr.v[i/DIM][XX];
390 c1[i+YY][nframes] = fr.v[i/DIM][YY];
391 c1[i+ZZ][nframes] = fr.v[i/DIM][ZZ];
398 while (read_next_frame(oenv, status, &fr));
402 dt = (t1-t0)/(nframes-1);
409 printf("Going to do %d fourier transforms of length %d. Hang on.\n",
412 low_do_autocorr(NULL, oenv, NULL, nframes, gnx, nframes, c1, dt, eacNormal, 0, FALSE,
413 FALSE, FALSE, -1, -1, 0);
415 for (j = 0; (j < DOS_NR); j++)
417 snew(dos[j], nframes+4);
422 printf("Going to merge the ACFs into the mass-weighted and plain ACF\n");
424 for (i = 0; (i < gnx); i += DIM)
426 mi = top.atoms.atom[i/DIM].m;
427 for (j = 0; (j < nframes/2); j++)
429 c1j = (c1[i+XX][j] + c1[i+YY][j] + c1[i+ZZ][j]);
430 dos[VACF][j] += c1j/Natom;
431 dos[MVACF][j] += mi*c1j;
434 fp = xvgropen(opt2fn("-vacf", NFILE, fnm), "Velocity ACF",
435 "Time (ps)", "C(t)", oenv);
437 for (j = 0; (j < nframes/2); j++)
440 fprintf(fp, "%10g %10g\n", tt[j], dos[VACF][j]);
443 fp = xvgropen(opt2fn("-mvacf", NFILE, fnm), "Mass-weighted velocity ACF",
444 "Time (ps)", "C(t)", oenv);
445 for (j = 0; (j < nframes/2); j++)
447 fprintf(fp, "%10g %10g\n", tt[j], dos[MVACF][j]);
451 if ((fftcode = gmx_fft_init_1d_real(&fft, nframes/2,
452 GMX_FFT_FLAG_NONE)) != 0)
454 gmx_fatal(FARGS, "gmx_fft_init_1d_real returned %d", fftcode);
456 if ((fftcode = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX,
457 (void *)dos[MVACF], (void *)dos[DOS])) != 0)
459 gmx_fatal(FARGS, "gmx_fft_1d_real returned %d", fftcode);
462 /* First compute the DoS */
463 /* Magic factor of 8 included now. */
467 for (j = 0; (j < nframes/4); j++)
470 dos2 += sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]);
473 dos[DOS][j] = bfac*sqrt(sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]));
477 dos[DOS][j] = bfac*dos[DOS][2*j];
481 dostot = evaluate_integral(nframes/4, nu, dos[DOS], NULL, nframes/4, &stddev);
484 for (j = 0; (j < nframes/4); j++)
486 dos[DOS][j] *= 3*Natom/dostot;
493 /* Note this eqn. is incorrect in Pascal2011a! */
494 Delta = ((2*DoS0/(9*Natom))*sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*
495 pow((Natom/V), 1.0/3.0)*pow(6/M_PI, 2.0/3.0));
496 f = calc_fluidicity(Delta, toler);
497 y = calc_y(f, Delta, toler);
498 z = calc_compress(y);
499 Sig = BOLTZ*(5.0/2.0+log(2*M_PI*BOLTZ*Temp/(sqr(PLANCK))*V/(f*Natom)));
500 Shs = Sig+calc_Shs(f, y);
501 rho = (tmass*AMU)/(V*NANO*NANO*NANO);
502 sigHS = pow(6*y*V/(M_PI*Natom), 1.0/3.0);
504 fprintf(fplog, "System = \"%s\"\n", title);
505 fprintf(fplog, "Nmol = %d\n", Nmol);
506 fprintf(fplog, "Natom = %d\n", Natom);
507 fprintf(fplog, "dt = %g ps\n", dt);
508 fprintf(fplog, "tmass = %g amu\n", tmass);
509 fprintf(fplog, "V = %g nm^3\n", V);
510 fprintf(fplog, "rho = %g g/l\n", rho);
511 fprintf(fplog, "T = %g K\n", Temp);
512 fprintf(fplog, "beta = %g mol/kJ\n", beta);
514 fprintf(fplog, "\nDoS parameters\n");
515 fprintf(fplog, "Delta = %g\n", Delta);
516 fprintf(fplog, "fluidicity = %g\n", f);
517 fprintf(fplog, "hard sphere packing fraction = %g\n", y);
518 fprintf(fplog, "hard sphere compressibility = %g\n", z);
519 fprintf(fplog, "ideal gas entropy = %g\n", Sig);
520 fprintf(fplog, "hard sphere entropy = %g\n", Shs);
521 fprintf(fplog, "sigma_HS = %g nm\n", sigHS);
522 fprintf(fplog, "DoS0 = %g\n", DoS0);
523 fprintf(fplog, "Dos2 = %g\n", dos2);
524 fprintf(fplog, "DoSTot = %g\n", dostot);
526 /* Now compute solid (2) and diffusive (3) components */
527 fp = xvgropen(opt2fn("-dos", NFILE, fnm), "Density of states",
528 bRecip ? "E (cm\\S-1\\N)" : "\\f{12}n\\f{4} (1/ps)",
529 "\\f{4}S(\\f{12}n\\f{4})", oenv);
530 xvgr_legend(fp, asize(DoSlegend), DoSlegend, oenv);
531 recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;
532 for (j = 0; (j < nframes/4); j++)
534 dos[DOS_DIFF][j] = DoS0/(1+sqr(DoS0*M_PI*nu[j]/(6*f*Natom)));
535 dos[DOS_SOLID][j] = dos[DOS][j]-dos[DOS_DIFF][j];
536 fprintf(fp, "%10g %10g %10g %10g\n",
538 dos[DOS][j]/recip_fac,
539 dos[DOS_SOLID][j]/recip_fac,
540 dos[DOS_DIFF][j]/recip_fac);
544 /* Finally analyze the results! */
546 wSdiff = Shs/(3*BOLTZ); /* Is this correct? */
548 wAdiff = wEdiff-wSdiff;
549 for (j = 0; (j < nframes/4); j++)
551 dos[DOS_CP][j] = (dos[DOS_DIFF][j]*wCdiff +
552 dos[DOS_SOLID][j]*wCsolid(nu[j], beta));
553 dos[DOS_S][j] = (dos[DOS_DIFF][j]*wSdiff +
554 dos[DOS_SOLID][j]*wSsolid(nu[j], beta));
555 dos[DOS_A][j] = (dos[DOS_DIFF][j]*wAdiff +
556 dos[DOS_SOLID][j]*wAsolid(nu[j], beta));
557 dos[DOS_E][j] = (dos[DOS_DIFF][j]*wEdiff +
558 dos[DOS_SOLID][j]*wEsolid(nu[j], beta));
560 DiffCoeff = evaluate_integral(nframes/2, tt, dos[VACF], NULL, nframes/2, &stddev);
561 DiffCoeff = 1000*DiffCoeff/3.0;
562 fprintf(fplog, "Diffusion coefficient from VACF %g 10^-5 cm^2/s\n",
564 fprintf(fplog, "Diffusion coefficient from DoS %g 10^-5 cm^2/s\n",
565 1000*DoS0/(12*tmass*beta));
567 cP = BOLTZ * evaluate_integral(nframes/4, nu, dos[DOS_CP], NULL,
569 fprintf(fplog, "Heat capacity %g J/mol K\n", 1000*cP/Nmol);
572 S = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_S],NULL,
574 fprintf(fplog,"Entropy %g J/mol K\n",1000*S/Nmol);
575 A = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_A],NULL,
577 fprintf(fplog,"Helmholtz energy %g kJ/mol\n",A/Nmol);
578 E = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_E],NULL,
580 fprintf(fplog,"Internal energy %g kJ/mol\n",E/Nmol);
582 fprintf(fplog, "\nArrivederci!\n");
583 gmx_fio_fclose(fplog);
585 do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");