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
41 #include "gromacs/fileio/confio.h"
43 #include "gmx_fatal.h"
44 #include "gromacs/fileio/futil.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 fprintf(fp, "@ world 1e-05, 0, 1000, 1\n");
227 fprintf(fp, "@ xaxes scale Logarithmic\n");
228 for (Delta = 1e-5; (Delta <= 1000); Delta *= DD)
230 f = calc_fluidicity(Delta, toler);
231 y = calc_y(f, Delta, toler);
232 fprintf(fp, "%10g %10g %10g %10g\n", Delta, f, f*y, y);
237 static void dump_w(output_env_t oenv, real beta)
241 const char *leg[] = { "wCv", "wS", "wA", "wE" };
243 fp = xvgropen("w.xvg", "Fig. 1, Berens1983a", "\\f{12}b\\f{4}h\\f{12}n",
245 xvgr_legend(fp, asize(leg), leg, oenv);
246 for (nu = 1; (nu < 100); nu += 0.05)
248 fprintf(fp, "%10g %10g %10g %10g %10g\n", beta*PLANCK*nu,
249 wCsolid(nu, beta), wSsolid(nu, beta),
250 wAsolid(nu, beta), wEsolid(nu, beta));
255 int gmx_dos(int argc, char *argv[])
257 const char *desc[] = {
258 "[THISMODULE] computes the Density of States from a simulations.",
259 "In order for this to be meaningful the velocities must be saved",
260 "in the trajecotry with sufficiently high frequency such as to cover",
261 "all vibrations. For flexible systems that would be around a few fs",
262 "between saving. Properties based on the DoS are printed on the",
265 const char *bugs[] = {
266 "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)."
277 int nV, nframes, n_alloc, i, j, k, l, fftcode, Nmol, Natom;
278 double rho, dt, V2sum, Vsum, V, tmass, dostot, dos2, dosabs;
279 real **c1, **dos, mi, beta, bfac, *nu, *tt, stddev, c1j;
282 double cP, S, A, E, DiffCoeff, Delta, f, y, z, sigHS, Shs, Sig, DoS0, recip_fac;
283 double wCdiff, wSdiff, wAdiff, wEdiff;
285 static gmx_bool bVerbose = TRUE, bAbsolute = FALSE, bNormalize = FALSE;
286 static gmx_bool bRecip = FALSE, bDump = FALSE;
287 static real Temp = 298.15, toler = 1e-6;
289 { "-v", FALSE, etBOOL, {&bVerbose},
290 "Be loud and noisy." },
291 { "-recip", FALSE, etBOOL, {&bRecip},
292 "Use cm^-1 on X-axis instead of 1/ps for DoS plots." },
293 { "-abs", FALSE, etBOOL, {&bAbsolute},
294 "Use the absolute value of the Fourier transform of the VACF as the Density of States. Default is to use the real component only" },
295 { "-normdos", FALSE, etBOOL, {&bNormalize},
296 "Normalize the DoS such that it adds up to 3N. This is a hack that should not be necessary." },
297 { "-T", FALSE, etREAL, {&Temp},
298 "Temperature in the simulation" },
299 { "-toler", FALSE, etREAL, {&toler},
300 "[HIDDEN]Tolerance when computing the fluidicity using bisection algorithm" },
301 { "-dump", FALSE, etBOOL, {&bDump},
302 "[HIDDEN]Dump the y/fy plot corresponding to Fig. 2 inLin2003a and the and the weighting functions corresponding to Fig. 1 in Berens1983a." }
306 { efTRN, "-f", NULL, ffREAD },
307 { efTPX, "-s", NULL, ffREAD },
308 { efNDX, NULL, NULL, ffOPTRD },
309 { efXVG, "-vacf", "vacf", ffWRITE },
310 { efXVG, "-mvacf", "mvacf", ffWRITE },
311 { efXVG, "-dos", "dos", ffWRITE },
312 { efLOG, "-g", "dos", ffWRITE },
314 #define NFILE asize(fnm)
317 const char *DoSlegend[] = {
318 "DoS(v)", "DoS(v)[Solid]", "DoS(v)[Diff]"
322 ppa = add_acf_pargs(&npargs, pa);
323 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
324 NFILE, fnm, npargs, ppa, asize(desc), desc,
325 asize(bugs), bugs, &oenv))
330 beta = 1/(Temp*BOLTZ);
333 printf("Dumping reference figures. Thanks for your patience.\n");
334 dump_fy(oenv, toler);
339 fplog = gmx_fio_fopen(ftp2fn(efLOG, NFILE, fnm), "w");
340 fprintf(fplog, "Doing density of states analysis based on trajectory.\n");
341 please_cite(fplog, "Pascal2011a");
342 please_cite(fplog, "Caleman2011b");
344 read_tps_conf(ftp2fn(efTPX, NFILE, fnm), title, &top, &ePBC, NULL, NULL, box,
348 for (i = 0; (i < top.atoms.nr); i++)
350 tmass += top.atoms.atom[i].m;
353 Natom = top.atoms.nr;
357 /* Correlation stuff */
359 for (i = 0; (i < gnx); i++)
364 read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE, fnm), &fr, TRX_NEED_V);
380 if (nframes >= n_alloc)
383 for (i = 0; i < gnx; i++)
385 srenew(c1[i], n_alloc);
388 for (i = 0; i < gnx; i += DIM)
390 c1[i+XX][nframes] = fr.v[i/DIM][XX];
391 c1[i+YY][nframes] = fr.v[i/DIM][YY];
392 c1[i+ZZ][nframes] = fr.v[i/DIM][ZZ];
399 while (read_next_frame(oenv, status, &fr));
403 dt = (t1-t0)/(nframes-1);
410 printf("Going to do %d fourier transforms of length %d. Hang on.\n",
413 low_do_autocorr(NULL, oenv, NULL, nframes, gnx, nframes, c1, dt, eacNormal, 0, FALSE,
414 FALSE, FALSE, -1, -1, 0);
416 for (j = 0; (j < DOS_NR); j++)
418 snew(dos[j], nframes+4);
423 printf("Going to merge the ACFs into the mass-weighted and plain ACF\n");
425 for (i = 0; (i < gnx); i += DIM)
427 mi = top.atoms.atom[i/DIM].m;
428 for (j = 0; (j < nframes/2); j++)
430 c1j = (c1[i+XX][j] + c1[i+YY][j] + c1[i+ZZ][j]);
431 dos[VACF][j] += c1j/Natom;
432 dos[MVACF][j] += mi*c1j;
435 fp = xvgropen(opt2fn("-vacf", NFILE, fnm), "Velocity ACF",
436 "Time (ps)", "C(t)", oenv);
438 for (j = 0; (j < nframes/2); j++)
441 fprintf(fp, "%10g %10g\n", tt[j], dos[VACF][j]);
444 fp = xvgropen(opt2fn("-mvacf", NFILE, fnm), "Mass-weighted velocity ACF",
445 "Time (ps)", "C(t)", oenv);
446 for (j = 0; (j < nframes/2); j++)
448 fprintf(fp, "%10g %10g\n", tt[j], dos[MVACF][j]);
452 if ((fftcode = gmx_fft_init_1d_real(&fft, nframes/2,
453 GMX_FFT_FLAG_NONE)) != 0)
455 gmx_fatal(FARGS, "gmx_fft_init_1d_real returned %d", fftcode);
457 if ((fftcode = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX,
458 (void *)dos[MVACF], (void *)dos[DOS])) != 0)
460 gmx_fatal(FARGS, "gmx_fft_1d_real returned %d", fftcode);
463 /* First compute the DoS */
464 /* Magic factor of 8 included now. */
468 for (j = 0; (j < nframes/4); j++)
471 dos2 += sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]);
474 dos[DOS][j] = bfac*sqrt(sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]));
478 dos[DOS][j] = bfac*dos[DOS][2*j];
482 dostot = evaluate_integral(nframes/4, nu, dos[DOS], NULL, nframes/4, &stddev);
485 for (j = 0; (j < nframes/4); j++)
487 dos[DOS][j] *= 3*Natom/dostot;
494 /* Note this eqn. is incorrect in Pascal2011a! */
495 Delta = ((2*DoS0/(9*Natom))*sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*
496 pow((Natom/V), 1.0/3.0)*pow(6/M_PI, 2.0/3.0));
497 f = calc_fluidicity(Delta, toler);
498 y = calc_y(f, Delta, toler);
499 z = calc_compress(y);
500 Sig = BOLTZ*(5.0/2.0+log(2*M_PI*BOLTZ*Temp/(sqr(PLANCK))*V/(f*Natom)));
501 Shs = Sig+calc_Shs(f, y);
502 rho = (tmass*AMU)/(V*NANO*NANO*NANO);
503 sigHS = pow(6*y*V/(M_PI*Natom), 1.0/3.0);
505 fprintf(fplog, "System = \"%s\"\n", title);
506 fprintf(fplog, "Nmol = %d\n", Nmol);
507 fprintf(fplog, "Natom = %d\n", Natom);
508 fprintf(fplog, "dt = %g ps\n", dt);
509 fprintf(fplog, "tmass = %g amu\n", tmass);
510 fprintf(fplog, "V = %g nm^3\n", V);
511 fprintf(fplog, "rho = %g g/l\n", rho);
512 fprintf(fplog, "T = %g K\n", Temp);
513 fprintf(fplog, "beta = %g mol/kJ\n", beta);
515 fprintf(fplog, "\nDoS parameters\n");
516 fprintf(fplog, "Delta = %g\n", Delta);
517 fprintf(fplog, "fluidicity = %g\n", f);
518 fprintf(fplog, "hard sphere packing fraction = %g\n", y);
519 fprintf(fplog, "hard sphere compressibility = %g\n", z);
520 fprintf(fplog, "ideal gas entropy = %g\n", Sig);
521 fprintf(fplog, "hard sphere entropy = %g\n", Shs);
522 fprintf(fplog, "sigma_HS = %g nm\n", sigHS);
523 fprintf(fplog, "DoS0 = %g\n", DoS0);
524 fprintf(fplog, "Dos2 = %g\n", dos2);
525 fprintf(fplog, "DoSTot = %g\n", dostot);
527 /* Now compute solid (2) and diffusive (3) components */
528 fp = xvgropen(opt2fn("-dos", NFILE, fnm), "Density of states",
529 bRecip ? "E (cm\\S-1\\N)" : "\\f{12}n\\f{4} (1/ps)",
530 "\\f{4}S(\\f{12}n\\f{4})", oenv);
531 xvgr_legend(fp, asize(DoSlegend), DoSlegend, oenv);
532 recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;
533 for (j = 0; (j < nframes/4); j++)
535 dos[DOS_DIFF][j] = DoS0/(1+sqr(DoS0*M_PI*nu[j]/(6*f*Natom)));
536 dos[DOS_SOLID][j] = dos[DOS][j]-dos[DOS_DIFF][j];
537 fprintf(fp, "%10g %10g %10g %10g\n",
539 dos[DOS][j]/recip_fac,
540 dos[DOS_SOLID][j]/recip_fac,
541 dos[DOS_DIFF][j]/recip_fac);
545 /* Finally analyze the results! */
547 wSdiff = Shs/(3*BOLTZ); /* Is this correct? */
549 wAdiff = wEdiff-wSdiff;
550 for (j = 0; (j < nframes/4); j++)
552 dos[DOS_CP][j] = (dos[DOS_DIFF][j]*wCdiff +
553 dos[DOS_SOLID][j]*wCsolid(nu[j], beta));
554 dos[DOS_S][j] = (dos[DOS_DIFF][j]*wSdiff +
555 dos[DOS_SOLID][j]*wSsolid(nu[j], beta));
556 dos[DOS_A][j] = (dos[DOS_DIFF][j]*wAdiff +
557 dos[DOS_SOLID][j]*wAsolid(nu[j], beta));
558 dos[DOS_E][j] = (dos[DOS_DIFF][j]*wEdiff +
559 dos[DOS_SOLID][j]*wEsolid(nu[j], beta));
561 DiffCoeff = evaluate_integral(nframes/2, tt, dos[VACF], NULL, nframes/2, &stddev);
562 DiffCoeff = 1000*DiffCoeff/3.0;
563 fprintf(fplog, "Diffusion coefficient from VACF %g 10^-5 cm^2/s\n",
565 fprintf(fplog, "Diffusion coefficient from DoS %g 10^-5 cm^2/s\n",
566 1000*DoS0/(12*tmass*beta));
568 cP = BOLTZ * evaluate_integral(nframes/4, nu, dos[DOS_CP], NULL,
570 fprintf(fplog, "Heat capacity %g J/mol K\n", 1000*cP/Nmol);
573 S = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_S],NULL,
575 fprintf(fplog,"Entropy %g J/mol K\n",1000*S/Nmol);
576 A = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_A],NULL,
578 fprintf(fplog,"Helmholtz energy %g kJ/mol\n",A/Nmol);
579 E = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_E],NULL,
581 fprintf(fplog,"Internal energy %g kJ/mol\n",E/Nmol);
583 fprintf(fplog, "\nArrivederci!\n");
584 gmx_fio_fclose(fplog);
586 do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");