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.
41 #include "gromacs/fileio/confio.h"
43 #include "gmx_fatal.h"
44 #include "gromacs/fileio/futil.h"
47 #include "gromacs/math/utilities.h"
50 #include "gromacs/utility/smalloc.h"
51 #include "gromacs/commandline/pargs.h"
60 #include "gromacs/fft/fft.h"
61 #include "gromacs/fileio/trxio.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 if (output_env_get_print_xvgr_codes(oenv))
227 fprintf(fp, "@ world 1e-05, 0, 1000, 1\n");
228 fprintf(fp, "@ xaxes scale Logarithmic\n");
230 for (Delta = 1e-5; (Delta <= 1000); Delta *= DD)
232 f = calc_fluidicity(Delta, toler);
233 y = calc_y(f, Delta, toler);
234 fprintf(fp, "%10g %10g %10g %10g\n", Delta, f, f*y, y);
239 static void dump_w(output_env_t oenv, real beta)
243 const char *leg[] = { "wCv", "wS", "wA", "wE" };
245 fp = xvgropen("w.xvg", "Fig. 1, Berens1983a", "\\f{12}b\\f{4}h\\f{12}n",
247 xvgr_legend(fp, asize(leg), leg, oenv);
248 for (nu = 1; (nu < 100); nu += 0.05)
250 fprintf(fp, "%10g %10g %10g %10g %10g\n", beta*PLANCK*nu,
251 wCsolid(nu, beta), wSsolid(nu, beta),
252 wAsolid(nu, beta), wEsolid(nu, beta));
257 int gmx_dos(int argc, char *argv[])
259 const char *desc[] = {
260 "[THISMODULE] computes the Density of States from a simulations.",
261 "In order for this to be meaningful the velocities must be saved",
262 "in the trajecotry with sufficiently high frequency such as to cover",
263 "all vibrations. For flexible systems that would be around a few fs",
264 "between saving. Properties based on the DoS are printed on the",
267 const char *bugs[] = {
268 "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)."
279 int nV, nframes, n_alloc, i, j, k, l, fftcode, Nmol, Natom;
280 double rho, dt, V2sum, Vsum, V, tmass, dostot, dos2, dosabs;
281 real **c1, **dos, mi, beta, bfac, *nu, *tt, stddev, c1j;
284 double cP, S, A, E, DiffCoeff, Delta, f, y, z, sigHS, Shs, Sig, DoS0, recip_fac;
285 double wCdiff, wSdiff, wAdiff, wEdiff;
287 static gmx_bool bVerbose = TRUE, bAbsolute = FALSE, bNormalize = FALSE;
288 static gmx_bool bRecip = FALSE, bDump = FALSE;
289 static real Temp = 298.15, toler = 1e-6;
291 { "-v", FALSE, etBOOL, {&bVerbose},
292 "Be loud and noisy." },
293 { "-recip", FALSE, etBOOL, {&bRecip},
294 "Use cm^-1 on X-axis instead of 1/ps for DoS plots." },
295 { "-abs", FALSE, etBOOL, {&bAbsolute},
296 "Use the absolute value of the Fourier transform of the VACF as the Density of States. Default is to use the real component only" },
297 { "-normdos", FALSE, etBOOL, {&bNormalize},
298 "Normalize the DoS such that it adds up to 3N. This is a hack that should not be necessary." },
299 { "-T", FALSE, etREAL, {&Temp},
300 "Temperature in the simulation" },
301 { "-toler", FALSE, etREAL, {&toler},
302 "[HIDDEN]Tolerance when computing the fluidicity using bisection algorithm" },
303 { "-dump", FALSE, etBOOL, {&bDump},
304 "[HIDDEN]Dump the y/fy plot corresponding to Fig. 2 inLin2003a and the and the weighting functions corresponding to Fig. 1 in Berens1983a." }
308 { efTRN, "-f", NULL, ffREAD },
309 { efTPX, "-s", NULL, ffREAD },
310 { efNDX, NULL, NULL, ffOPTRD },
311 { efXVG, "-vacf", "vacf", ffWRITE },
312 { efXVG, "-mvacf", "mvacf", ffWRITE },
313 { efXVG, "-dos", "dos", ffWRITE },
314 { efLOG, "-g", "dos", ffWRITE },
316 #define NFILE asize(fnm)
319 const char *DoSlegend[] = {
320 "DoS(v)", "DoS(v)[Solid]", "DoS(v)[Diff]"
324 ppa = add_acf_pargs(&npargs, pa);
325 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
326 NFILE, fnm, npargs, ppa, asize(desc), desc,
327 asize(bugs), bugs, &oenv))
332 beta = 1/(Temp*BOLTZ);
335 printf("Dumping reference figures. Thanks for your patience.\n");
336 dump_fy(oenv, toler);
341 fplog = gmx_fio_fopen(ftp2fn(efLOG, NFILE, fnm), "w");
342 fprintf(fplog, "Doing density of states analysis based on trajectory.\n");
343 please_cite(fplog, "Pascal2011a");
344 please_cite(fplog, "Caleman2011b");
346 read_tps_conf(ftp2fn(efTPX, NFILE, fnm), title, &top, &ePBC, NULL, NULL, box,
350 for (i = 0; (i < top.atoms.nr); i++)
352 tmass += top.atoms.atom[i].m;
355 Natom = top.atoms.nr;
359 /* Correlation stuff */
361 for (i = 0; (i < gnx); i++)
366 read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE, fnm), &fr, TRX_NEED_V);
382 if (nframes >= n_alloc)
385 for (i = 0; i < gnx; i++)
387 srenew(c1[i], n_alloc);
390 for (i = 0; i < gnx; i += DIM)
392 c1[i+XX][nframes] = fr.v[i/DIM][XX];
393 c1[i+YY][nframes] = fr.v[i/DIM][YY];
394 c1[i+ZZ][nframes] = fr.v[i/DIM][ZZ];
401 while (read_next_frame(oenv, status, &fr));
405 dt = (t1-t0)/(nframes-1);
412 printf("Going to do %d fourier transforms of length %d. Hang on.\n",
415 low_do_autocorr(NULL, oenv, NULL, nframes, gnx, nframes, c1, dt, eacNormal, 0, FALSE,
416 FALSE, FALSE, -1, -1, 0);
418 for (j = 0; (j < DOS_NR); j++)
420 snew(dos[j], nframes+4);
425 printf("Going to merge the ACFs into the mass-weighted and plain ACF\n");
427 for (i = 0; (i < gnx); i += DIM)
429 mi = top.atoms.atom[i/DIM].m;
430 for (j = 0; (j < nframes/2); j++)
432 c1j = (c1[i+XX][j] + c1[i+YY][j] + c1[i+ZZ][j]);
433 dos[VACF][j] += c1j/Natom;
434 dos[MVACF][j] += mi*c1j;
437 fp = xvgropen(opt2fn("-vacf", NFILE, fnm), "Velocity ACF",
438 "Time (ps)", "C(t)", oenv);
440 for (j = 0; (j < nframes/2); j++)
443 fprintf(fp, "%10g %10g\n", tt[j], dos[VACF][j]);
446 fp = xvgropen(opt2fn("-mvacf", NFILE, fnm), "Mass-weighted velocity ACF",
447 "Time (ps)", "C(t)", oenv);
448 for (j = 0; (j < nframes/2); j++)
450 fprintf(fp, "%10g %10g\n", tt[j], dos[MVACF][j]);
454 if ((fftcode = gmx_fft_init_1d_real(&fft, nframes/2,
455 GMX_FFT_FLAG_NONE)) != 0)
457 gmx_fatal(FARGS, "gmx_fft_init_1d_real returned %d", fftcode);
459 if ((fftcode = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX,
460 (void *)dos[MVACF], (void *)dos[DOS])) != 0)
462 gmx_fatal(FARGS, "gmx_fft_1d_real returned %d", fftcode);
465 /* First compute the DoS */
466 /* Magic factor of 8 included now. */
470 for (j = 0; (j < nframes/4); j++)
473 dos2 += sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]);
476 dos[DOS][j] = bfac*sqrt(sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]));
480 dos[DOS][j] = bfac*dos[DOS][2*j];
484 dostot = evaluate_integral(nframes/4, nu, dos[DOS], NULL, nframes/4, &stddev);
487 for (j = 0; (j < nframes/4); j++)
489 dos[DOS][j] *= 3*Natom/dostot;
496 /* Note this eqn. is incorrect in Pascal2011a! */
497 Delta = ((2*DoS0/(9*Natom))*sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*
498 pow((Natom/V), 1.0/3.0)*pow(6/M_PI, 2.0/3.0));
499 f = calc_fluidicity(Delta, toler);
500 y = calc_y(f, Delta, toler);
501 z = calc_compress(y);
502 Sig = BOLTZ*(5.0/2.0+log(2*M_PI*BOLTZ*Temp/(sqr(PLANCK))*V/(f*Natom)));
503 Shs = Sig+calc_Shs(f, y);
504 rho = (tmass*AMU)/(V*NANO*NANO*NANO);
505 sigHS = pow(6*y*V/(M_PI*Natom), 1.0/3.0);
507 fprintf(fplog, "System = \"%s\"\n", title);
508 fprintf(fplog, "Nmol = %d\n", Nmol);
509 fprintf(fplog, "Natom = %d\n", Natom);
510 fprintf(fplog, "dt = %g ps\n", dt);
511 fprintf(fplog, "tmass = %g amu\n", tmass);
512 fprintf(fplog, "V = %g nm^3\n", V);
513 fprintf(fplog, "rho = %g g/l\n", rho);
514 fprintf(fplog, "T = %g K\n", Temp);
515 fprintf(fplog, "beta = %g mol/kJ\n", beta);
517 fprintf(fplog, "\nDoS parameters\n");
518 fprintf(fplog, "Delta = %g\n", Delta);
519 fprintf(fplog, "fluidicity = %g\n", f);
520 fprintf(fplog, "hard sphere packing fraction = %g\n", y);
521 fprintf(fplog, "hard sphere compressibility = %g\n", z);
522 fprintf(fplog, "ideal gas entropy = %g\n", Sig);
523 fprintf(fplog, "hard sphere entropy = %g\n", Shs);
524 fprintf(fplog, "sigma_HS = %g nm\n", sigHS);
525 fprintf(fplog, "DoS0 = %g\n", DoS0);
526 fprintf(fplog, "Dos2 = %g\n", dos2);
527 fprintf(fplog, "DoSTot = %g\n", dostot);
529 /* Now compute solid (2) and diffusive (3) components */
530 fp = xvgropen(opt2fn("-dos", NFILE, fnm), "Density of states",
531 bRecip ? "E (cm\\S-1\\N)" : "\\f{12}n\\f{4} (1/ps)",
532 "\\f{4}S(\\f{12}n\\f{4})", oenv);
533 xvgr_legend(fp, asize(DoSlegend), DoSlegend, oenv);
534 recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;
535 for (j = 0; (j < nframes/4); j++)
537 dos[DOS_DIFF][j] = DoS0/(1+sqr(DoS0*M_PI*nu[j]/(6*f*Natom)));
538 dos[DOS_SOLID][j] = dos[DOS][j]-dos[DOS_DIFF][j];
539 fprintf(fp, "%10g %10g %10g %10g\n",
541 dos[DOS][j]/recip_fac,
542 dos[DOS_SOLID][j]/recip_fac,
543 dos[DOS_DIFF][j]/recip_fac);
547 /* Finally analyze the results! */
549 wSdiff = Shs/(3*BOLTZ); /* Is this correct? */
551 wAdiff = wEdiff-wSdiff;
552 for (j = 0; (j < nframes/4); j++)
554 dos[DOS_CP][j] = (dos[DOS_DIFF][j]*wCdiff +
555 dos[DOS_SOLID][j]*wCsolid(nu[j], beta));
556 dos[DOS_S][j] = (dos[DOS_DIFF][j]*wSdiff +
557 dos[DOS_SOLID][j]*wSsolid(nu[j], beta));
558 dos[DOS_A][j] = (dos[DOS_DIFF][j]*wAdiff +
559 dos[DOS_SOLID][j]*wAsolid(nu[j], beta));
560 dos[DOS_E][j] = (dos[DOS_DIFF][j]*wEdiff +
561 dos[DOS_SOLID][j]*wEsolid(nu[j], beta));
563 DiffCoeff = evaluate_integral(nframes/2, tt, dos[VACF], NULL, nframes/2, &stddev);
564 DiffCoeff = 1000*DiffCoeff/3.0;
565 fprintf(fplog, "Diffusion coefficient from VACF %g 10^-5 cm^2/s\n",
567 fprintf(fplog, "Diffusion coefficient from DoS %g 10^-5 cm^2/s\n",
568 1000*DoS0/(12*tmass*beta));
570 cP = BOLTZ * evaluate_integral(nframes/4, nu, dos[DOS_CP], NULL,
572 fprintf(fplog, "Heat capacity %g J/mol K\n", 1000*cP/Nmol);
575 S = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_S],NULL,
577 fprintf(fplog,"Entropy %g J/mol K\n",1000*S/Nmol);
578 A = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_A],NULL,
580 fprintf(fplog,"Helmholtz energy %g kJ/mol\n",A/Nmol);
581 E = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_E],NULL,
583 fprintf(fplog,"Internal energy %g kJ/mol\n",E/Nmol);
585 fprintf(fplog, "\nArrivederci!\n");
586 gmx_fio_fclose(fplog);
588 do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");