2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2011,2012,2013,2014,2015,2016, 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.
42 #include "gromacs/commandline/pargs.h"
43 #include "gromacs/commandline/viewit.h"
44 #include "gromacs/correlationfunctions/autocorr.h"
45 #include "gromacs/correlationfunctions/integrate.h"
46 #include "gromacs/fft/fft.h"
47 #include "gromacs/fileio/confio.h"
48 #include "gromacs/fileio/gmxfio.h"
49 #include "gromacs/fileio/trxio.h"
50 #include "gromacs/fileio/xvgr.h"
51 #include "gromacs/gmxana/gmx_ana.h"
52 #include "gromacs/math/functions.h"
53 #include "gromacs/math/units.h"
54 #include "gromacs/math/utilities.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/topology/index.h"
57 #include "gromacs/topology/topology.h"
58 #include "gromacs/trajectory/trajectoryframe.h"
59 #include "gromacs/utility/arraysize.h"
60 #include "gromacs/utility/fatalerror.h"
61 #include "gromacs/utility/futil.h"
62 #include "gromacs/utility/pleasecite.h"
63 #include "gromacs/utility/smalloc.h"
66 VACF, MVACF, DOS, DOS_SOLID, DOS_DIFF, DOS_CP, DOS_S, DOS_A, DOS_E, DOS_NR
69 static int calcMoleculesInIndexGroup(const t_block *mols, int natoms, const int *index, int nindex)
78 while (index[i] > mols->index[mol])
83 gmx_fatal(FARGS, "Atom index out of range: %d", index[i]+1);
86 for (j = mols->index[mol]; j < mols->index[mol+1]; j++)
90 gmx_fatal(FARGS, "The index group does not consist of whole molecules");
95 gmx_fatal(FARGS, "Index contains atom numbers larger than the topology");
103 static double FD(double Delta, double f)
105 return (2*std::pow(Delta, -4.5)*std::pow(f, 7.5) -
106 6*std::pow(Delta, -3.0)*std::pow(f, 5.0) -
107 std::pow(Delta, -1.5)*std::pow(f, 3.5) +
108 6*std::pow(Delta, -1.5)*std::pow(f, 2.5) +
112 static double YYY(double f, double y)
114 return (2*gmx::power3(y*f) - gmx::square(f)*y*(1+6*y) +
118 static double calc_compress(double y)
124 return ((1+y+gmx::square(y)-gmx::power3(y))/(gmx::power3(1-y)));
127 static double bisector(double Delta, double tol,
128 double ff0, double ff1,
129 double ff(double, double))
131 double fd, f, f0, f1;
132 double tolmin = 1e-8;
138 fprintf(stderr, "Unrealistic tolerance %g for bisector. Setting it to %g\n", tol, tolmin);
159 while ((f1-f0) > tol);
164 static double calc_fluidicity(double Delta, double tol)
166 return bisector(Delta, tol, 0, 1, FD);
169 static double calc_y(double f, double Delta, double toler)
173 y1 = std::pow(f/Delta, 1.5);
174 y2 = bisector(f, toler, 0, 10000, YYY);
175 if (std::abs((y1-y2)/(y1+y2)) > 100*toler)
177 fprintf(stderr, "Inconsistency computing y: y1 = %f, y2 = %f, using y1.\n",
184 static double calc_Shs(double f, double y)
188 return BOLTZ*(std::log(calc_compress(fy)) + fy*(3*fy-4)/gmx::square(1-fy));
191 static real wCsolid(real nu, real beta)
193 real bhn = beta*PLANCK*nu;
203 koko = gmx::square(1-ebn);
204 return gmx::square(bhn)*ebn/koko;
208 static real wSsolid(real nu, real beta)
210 real bhn = beta*PLANCK*nu;
218 return bhn/std::expm1(bhn) - std::log1p(-std::exp(-bhn));
222 static real wAsolid(real nu, real beta)
224 real bhn = beta*PLANCK*nu;
232 return std::log((1-std::exp(-bhn))/(std::exp(-bhn/2))) - std::log(bhn);
236 static real wEsolid(real nu, real beta)
238 real bhn = beta*PLANCK*nu;
246 return bhn/2 + bhn/std::expm1(bhn)-1;
250 int gmx_dos(int argc, char *argv[])
252 const char *desc[] = {
253 "[THISMODULE] computes the Density of States from a simulations.",
254 "In order for this to be meaningful the velocities must be saved",
255 "in the trajecotry with sufficiently high frequency such as to cover",
256 "all vibrations. For flexible systems that would be around a few fs",
257 "between saving. Properties based on the DoS are printed on the",
259 "Note that the density of states is calculated from the mass-weighted",
260 "autocorrelation, and by default only from the square of the real",
261 "component rather than absolute value. This means the shape can differ",
262 "substantially from the plain vibrational power spectrum you can",
263 "calculate with gmx velacc."
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)."
276 int nV, nframes, n_alloc, i, j, fftcode, Nmol, Natom;
277 double rho, dt, Vsum, V, tmass, dostot, dos2;
278 real **c1, **dos, mi, beta, bfac, *nu, *tt, stddev, c1j;
279 gmx_output_env_t *oenv;
281 double cP, DiffCoeff, Delta, f, y, z, sigHS, Shs, Sig, DoS0, recip_fac;
282 double wCdiff, wSdiff, wAdiff, wEdiff;
287 gmx_bool normalizeAutocorrelation;
289 static gmx_bool bVerbose = TRUE, bAbsolute = FALSE, bNormalizeDos = FALSE;
290 static gmx_bool bRecip = FALSE;
291 static real Temp = 298.15, toler = 1e-6;
292 int min_frames = 100;
295 { "-v", FALSE, etBOOL, {&bVerbose},
296 "Be loud and noisy." },
297 { "-recip", FALSE, etBOOL, {&bRecip},
298 "Use cm^-1 on X-axis instead of 1/ps for DoS plots." },
299 { "-abs", FALSE, etBOOL, {&bAbsolute},
300 "Use the absolute value of the Fourier transform of the VACF as the Density of States. Default is to use the real component only" },
301 { "-normdos", FALSE, etBOOL, {&bNormalizeDos},
302 "Normalize the DoS such that it adds up to 3N. This should usually not be necessary." },
303 { "-T", FALSE, etREAL, {&Temp},
304 "Temperature in the simulation" },
305 { "-toler", FALSE, etREAL, {&toler},
306 "[HIDDEN]Tolerance when computing the fluidicity using bisection algorithm" }
310 { efTRN, "-f", NULL, ffREAD },
311 { efTPR, "-s", NULL, ffREAD },
312 { efNDX, NULL, NULL, ffOPTRD },
313 { efXVG, "-vacf", "vacf", ffWRITE },
314 { efXVG, "-mvacf", "mvacf", ffWRITE },
315 { efXVG, "-dos", "dos", ffWRITE },
316 { efLOG, "-g", "dos", ffWRITE },
318 #define NFILE asize(fnm)
321 const char *DoSlegend[] = {
322 "DoS(v)", "DoS(v)[Solid]", "DoS(v)[Diff]"
326 ppa = add_acf_pargs(&npargs, pa);
327 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
328 NFILE, fnm, npargs, ppa, asize(desc), desc,
329 asize(bugs), bugs, &oenv))
334 beta = 1/(Temp*BOLTZ);
336 fplog = gmx_fio_fopen(ftp2fn(efLOG, NFILE, fnm), "w");
337 fprintf(fplog, "Doing density of states analysis based on trajectory.\n");
338 please_cite(fplog, "Pascal2011a");
339 please_cite(fplog, "Caleman2011b");
341 read_tps_conf(ftp2fn(efTPR, NFILE, fnm), &top, &ePBC, NULL, NULL, box, TRUE);
343 /* Handle index groups */
344 get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &grpNatoms, &index, &grpname);
348 for (i = 0; i < grpNatoms; i++)
350 tmass += top.atoms.atom[index[i]].m;
354 Nmol = calcMoleculesInIndexGroup(&top.mols, top.atoms.nr, index, grpNatoms);
357 /* Correlation stuff */
359 for (i = 0; (i < gnx); i++)
364 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[index[i/DIM]][XX];
390 c1[i+YY][nframes] = fr.v[index[i/DIM]][YY];
391 c1[i+ZZ][nframes] = fr.v[index[i/DIM]][ZZ];
398 while (read_next_frame(oenv, status, &fr));
402 if (nframes < min_frames)
404 gmx_fatal(FARGS, "You need at least %d frames in the trajectory and you only have %d.", min_frames, nframes);
406 dt = (t1-t0)/(nframes-1);
413 printf("Going to do %d fourier transforms of length %d. Hang on.\n",
416 /* Unfortunately the -normalize program option for the autocorrelation
417 * function calculation is added as a hack with a static variable in the
418 * autocorrelation.c source. That would work if we called the normal
419 * do_autocorr(), but this routine overrides that by directly calling
420 * the low-level functionality. That unfortunately leads to ignoring the
421 * default value for the option (which is to normalize).
422 * Since the absolute value seems to be important for the subsequent
423 * analysis below, we detect the value directly from the option, calculate
424 * the autocorrelation without normalization, and then apply the
425 * normalization just to the autocorrelation output
426 * (or not, if the user asked for a non-normalized autocorrelation).
428 normalizeAutocorrelation = opt2parg_bool("-normalize", npargs, ppa);
430 /* Note that we always disable normalization here, regardless of user settings */
431 low_do_autocorr(NULL, oenv, NULL, nframes, gnx, nframes, c1, dt, eacNormal, 0, FALSE,
432 FALSE, FALSE, -1, -1, 0);
434 for (j = 0; (j < DOS_NR); j++)
436 snew(dos[j], nframes+4);
441 printf("Going to merge the ACFs into the mass-weighted and plain ACF\n");
443 for (i = 0; (i < gnx); i += DIM)
445 mi = top.atoms.atom[index[i/DIM]].m;
446 for (j = 0; (j < nframes/2); j++)
448 c1j = (c1[i+XX][j] + c1[i+YY][j] + c1[i+ZZ][j]);
449 dos[VACF][j] += c1j/Natom;
450 dos[MVACF][j] += mi*c1j;
454 fp = xvgropen(opt2fn("-vacf", NFILE, fnm), "Velocity autocorrelation function",
455 "Time (ps)", "C(t)", oenv);
458 invNormalize = normalizeAutocorrelation ? 1.0/dos[VACF][0] : 1.0;
460 for (j = 0; (j < nframes/2); j++)
463 fprintf(fp, "%10g %10g\n", tt[j], dos[VACF][j] * invNormalize);
467 fp = xvgropen(opt2fn("-mvacf", NFILE, fnm), "Mass-weighted velocity autocorrelation function",
468 "Time (ps)", "C(t)", oenv);
470 invNormalize = normalizeAutocorrelation ? 1.0/dos[VACF][0] : 1.0;
472 for (j = 0; (j < nframes/2); j++)
474 fprintf(fp, "%10g %10g\n", tt[j], dos[MVACF][j] * invNormalize);
478 if ((fftcode = gmx_fft_init_1d_real(&fft, nframes/2,
479 GMX_FFT_FLAG_NONE)) != 0)
481 gmx_fatal(FARGS, "gmx_fft_init_1d_real returned %d", fftcode);
483 if ((fftcode = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX,
484 (void *)dos[MVACF], (void *)dos[DOS])) != 0)
486 gmx_fatal(FARGS, "gmx_fft_1d_real returned %d", fftcode);
489 /* First compute the DoS */
490 /* Magic factor of 8 included now. */
494 for (j = 0; (j < nframes/4); j++)
497 dos2 += gmx::square(dos[DOS][2*j]) + gmx::square(dos[DOS][2*j+1]);
500 dos[DOS][j] = bfac*std::hypot(dos[DOS][2*j], dos[DOS][2*j+1]);
504 dos[DOS][j] = bfac*dos[DOS][2*j];
508 dostot = evaluate_integral(nframes/4, nu, dos[DOS], NULL, nframes/4, &stddev);
511 for (j = 0; (j < nframes/4); j++)
513 dos[DOS][j] *= 3*Natom/dostot;
520 /* Note this eqn. is incorrect in Pascal2011a! */
521 Delta = ((2*DoS0/(9*Natom))*std::sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*
522 std::pow((Natom/V), 1.0/3.0)*std::pow(6.0/M_PI, 2.0/3.0));
523 f = calc_fluidicity(Delta, toler);
524 y = calc_y(f, Delta, toler);
525 z = calc_compress(y);
526 Sig = BOLTZ*(5.0/2.0+std::log(2*M_PI*BOLTZ*Temp/(gmx::square(PLANCK))*V/(f*Natom)));
527 Shs = Sig+calc_Shs(f, y);
528 rho = (tmass*AMU)/(V*NANO*NANO*NANO);
529 sigHS = std::cbrt(6*y*V/(M_PI*Natom));
531 fprintf(fplog, "System = \"%s\"\n", *top.name);
532 fprintf(fplog, "Nmol = %d\n", Nmol);
533 fprintf(fplog, "Natom = %d\n", Natom);
534 fprintf(fplog, "dt = %g ps\n", dt);
535 fprintf(fplog, "tmass = %g amu\n", tmass);
536 fprintf(fplog, "V = %g nm^3\n", V);
537 fprintf(fplog, "rho = %g g/l\n", rho);
538 fprintf(fplog, "T = %g K\n", Temp);
539 fprintf(fplog, "beta = %g mol/kJ\n", beta);
541 fprintf(fplog, "\nDoS parameters\n");
542 fprintf(fplog, "Delta = %g\n", Delta);
543 fprintf(fplog, "fluidicity = %g\n", f);
544 fprintf(fplog, "hard sphere packing fraction = %g\n", y);
545 fprintf(fplog, "hard sphere compressibility = %g\n", z);
546 fprintf(fplog, "ideal gas entropy = %g\n", Sig);
547 fprintf(fplog, "hard sphere entropy = %g\n", Shs);
548 fprintf(fplog, "sigma_HS = %g nm\n", sigHS);
549 fprintf(fplog, "DoS0 = %g\n", DoS0);
550 fprintf(fplog, "Dos2 = %g\n", dos2);
551 fprintf(fplog, "DoSTot = %g\n", dostot);
553 /* Now compute solid (2) and diffusive (3) components */
554 fp = xvgropen(opt2fn("-dos", NFILE, fnm), "Density of states",
555 bRecip ? "E (cm\\S-1\\N)" : "\\f{12}n\\f{4} (1/ps)",
556 "\\f{4}S(\\f{12}n\\f{4})", oenv);
557 xvgr_legend(fp, asize(DoSlegend), DoSlegend, oenv);
558 recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;
559 for (j = 0; (j < nframes/4); j++)
561 dos[DOS_DIFF][j] = DoS0/(1+gmx::square(DoS0*M_PI*nu[j]/(6*f*Natom)));
562 dos[DOS_SOLID][j] = dos[DOS][j]-dos[DOS_DIFF][j];
563 fprintf(fp, "%10g %10g %10g %10g\n",
565 dos[DOS][j]/recip_fac,
566 dos[DOS_SOLID][j]/recip_fac,
567 dos[DOS_DIFF][j]/recip_fac);
571 /* Finally analyze the results! */
573 wSdiff = Shs/(3*BOLTZ); /* Is this correct? */
575 wAdiff = wEdiff-wSdiff;
576 for (j = 0; (j < nframes/4); j++)
578 dos[DOS_CP][j] = (dos[DOS_DIFF][j]*wCdiff +
579 dos[DOS_SOLID][j]*wCsolid(nu[j], beta));
580 dos[DOS_S][j] = (dos[DOS_DIFF][j]*wSdiff +
581 dos[DOS_SOLID][j]*wSsolid(nu[j], beta));
582 dos[DOS_A][j] = (dos[DOS_DIFF][j]*wAdiff +
583 dos[DOS_SOLID][j]*wAsolid(nu[j], beta));
584 dos[DOS_E][j] = (dos[DOS_DIFF][j]*wEdiff +
585 dos[DOS_SOLID][j]*wEsolid(nu[j], beta));
587 DiffCoeff = evaluate_integral(nframes/2, tt, dos[VACF], NULL, nframes/2, &stddev);
588 DiffCoeff = 1000*DiffCoeff/3.0;
589 fprintf(fplog, "Diffusion coefficient from VACF %g 10^-5 cm^2/s\n",
591 fprintf(fplog, "Diffusion coefficient from DoS %g 10^-5 cm^2/s\n",
592 1000*DoS0/(12*tmass*beta));
594 cP = BOLTZ * evaluate_integral(nframes/4, nu, dos[DOS_CP], NULL,
596 fprintf(fplog, "Heat capacity %g J/mol K\n", 1000*cP/Nmol);
597 fprintf(fplog, "\nArrivederci!\n");
598 gmx_fio_fclose(fplog);
600 do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");