2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2011-2018, The GROMACS development team.
5 * Copyright (c) 2019,2020, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
43 #include "gromacs/commandline/pargs.h"
44 #include "gromacs/commandline/viewit.h"
45 #include "gromacs/correlationfunctions/autocorr.h"
46 #include "gromacs/correlationfunctions/integrate.h"
47 #include "gromacs/fft/fft.h"
48 #include "gromacs/fileio/confio.h"
49 #include "gromacs/fileio/gmxfio.h"
50 #include "gromacs/fileio/trxio.h"
51 #include "gromacs/fileio/xvgr.h"
52 #include "gromacs/gmxana/gmx_ana.h"
53 #include "gromacs/math/functions.h"
54 #include "gromacs/math/units.h"
55 #include "gromacs/math/utilities.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/pbcutil/pbc.h"
58 #include "gromacs/topology/index.h"
59 #include "gromacs/topology/topology.h"
60 #include "gromacs/trajectory/trajectoryframe.h"
61 #include "gromacs/utility/arraysize.h"
62 #include "gromacs/utility/fatalerror.h"
63 #include "gromacs/utility/futil.h"
64 #include "gromacs/utility/pleasecite.h"
65 #include "gromacs/utility/smalloc.h"
81 static int calcMoleculesInIndexGroup(const t_block* mols, int natoms, const int* index, int nindex)
90 while (index[i] > mols->index[mol])
95 gmx_fatal(FARGS, "Atom index out of range: %d", index[i] + 1);
98 for (j = mols->index[mol]; j < mols->index[mol + 1]; j++)
102 gmx_fatal(FARGS, "The index group does not consist of whole molecules");
107 gmx_fatal(FARGS, "Index contains atom numbers larger than the topology");
115 static double FD(double Delta, double f)
117 return (2 * std::pow(Delta, -4.5) * std::pow(f, 7.5)
118 - 6 * std::pow(Delta, -3.0) * std::pow(f, 5.0) - std::pow(Delta, -1.5) * std::pow(f, 3.5)
119 + 6 * std::pow(Delta, -1.5) * std::pow(f, 2.5) + 2 * f - 2);
122 static double YYY(double f, double y)
124 return (2 * gmx::power3(y * f) - gmx::square(f) * y * (1 + 6 * y) + (2 + 6 * y) * f - 2);
127 static double calc_compress(double y)
133 return ((1 + y + gmx::square(y) - gmx::power3(y)) / (gmx::power3(1 - y)));
136 static double bisector(double Delta, double tol, double ff0, double ff1, double ff(double, double))
138 double fd, f, f0, f1;
139 double tolmin = 1e-8;
145 fprintf(stderr, "Unrealistic tolerance %g for bisector. Setting it to %g\n", tol, tolmin);
165 } while ((f1 - f0) > tol);
170 static double calc_fluidicity(double Delta, double tol)
172 return bisector(Delta, tol, 0, 1, FD);
175 static double calc_y(double f, double Delta, double toler)
179 y1 = std::pow(f / Delta, 1.5);
180 y2 = bisector(f, toler, 0, 10000, YYY);
181 if (std::abs((y1 - y2) / (y1 + y2)) > 100 * toler)
183 fprintf(stderr, "Inconsistency computing y: y1 = %f, y2 = %f, using y1.\n", y1, y2);
189 static double calc_Shs(double f, double y)
193 return BOLTZ * (std::log(calc_compress(fy)) + fy * (3 * fy - 4) / gmx::square(1 - fy));
196 static real wCsolid(real nu, real beta)
198 real bhn = beta * PLANCK * nu;
208 koko = gmx::square(1 - ebn);
209 return gmx::square(bhn) * ebn / koko;
213 static real wSsolid(real nu, real beta)
215 real bhn = beta * PLANCK * nu;
223 return bhn / std::expm1(bhn) - std::log1p(-std::exp(-bhn));
227 static real wAsolid(real nu, real beta)
229 real bhn = beta * PLANCK * nu;
237 return std::log((1 - std::exp(-bhn)) / (std::exp(-bhn / 2))) - std::log(bhn);
241 static real wEsolid(real nu, real beta)
243 real bhn = beta * PLANCK * nu;
251 return bhn / 2 + bhn / std::expm1(bhn) - 1;
255 int gmx_dos(int argc, char* argv[])
257 const char* desc[] = { "[THISMODULE] 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",
263 "Note that the density of states is calculated from the mass-weighted",
264 "autocorrelation, and by default only from the square of the real",
265 "component rather than absolute value. This means the shape can differ",
266 "substantially from the plain vibrational power spectrum you can",
267 "calculate with gmx velacc." };
268 const char* bugs[] = {
269 "This program needs a lot of memory: total usage equals the number of atoms times "
270 "3 times number of frames times 4 (or 8 when run in double precision)."
274 PbcType pbcType = PbcType::Unset;
280 int nV, nframes, n_alloc, i, j, fftcode, Nmol, Natom;
281 double rho, dt, Vsum, V, tmass, dostot, dos2;
282 real ** c1, **dos, mi, beta, bfac, *nu, *tt, stddev, c1j;
283 gmx_output_env_t* oenv;
285 double cP, DiffCoeff, Delta, f, y, z, sigHS, Shs, Sig, DoS0, recip_fac;
286 double wCdiff, wSdiff, wAdiff, wEdiff;
291 gmx_bool normalizeAutocorrelation;
293 static gmx_bool bVerbose = TRUE, bAbsolute = FALSE, bNormalizeDos = FALSE;
294 static gmx_bool bRecip = FALSE;
295 static real Temp = 298.15, toler = 1e-6;
296 int min_frames = 100;
299 { "-v", FALSE, etBOOL, { &bVerbose }, "Be loud and noisy." },
304 "Use cm^-1 on X-axis instead of 1/ps for DoS plots." },
309 "Use the absolute value of the Fourier transform of the VACF as the Density of States. "
310 "Default is to use the real component only" },
315 "Normalize the DoS such that it adds up to 3N. This should usually not be necessary." },
316 { "-T", FALSE, etREAL, { &Temp }, "Temperature in the simulation" },
321 "[HIDDEN]Tolerance when computing the fluidicity using bisection algorithm" }
325 { efTRN, "-f", nullptr, ffREAD }, { efTPR, "-s", nullptr, ffREAD },
326 { efNDX, nullptr, nullptr, ffOPTRD }, { efXVG, "-vacf", "vacf", ffWRITE },
327 { efXVG, "-mvacf", "mvacf", ffWRITE }, { efXVG, "-dos", "dos", ffWRITE },
328 { efLOG, "-g", "dos", ffWRITE },
330 #define NFILE asize(fnm)
333 const char* DoSlegend[] = { "DoS(v)", "DoS(v)[Solid]", "DoS(v)[Diff]" };
336 ppa = add_acf_pargs(&npargs, pa);
337 if (!parse_common_args(
338 &argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs, &oenv))
344 beta = 1 / (Temp * BOLTZ);
346 fplog = gmx_fio_fopen(ftp2fn(efLOG, NFILE, fnm), "w");
347 fprintf(fplog, "Doing density of states analysis based on trajectory.\n");
348 please_cite(fplog, "Pascal2011a");
349 please_cite(fplog, "Caleman2011b");
351 read_tps_conf(ftp2fn(efTPR, NFILE, fnm), &top, &pbcType, nullptr, nullptr, box, TRUE);
353 /* Handle index groups */
354 get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &grpNatoms, &index, &grpname);
358 for (i = 0; i < grpNatoms; i++)
360 tmass += top.atoms.atom[index[i]].m;
364 Nmol = calcMoleculesInIndexGroup(&top.mols, top.atoms.nr, index, grpNatoms);
367 /* Correlation stuff */
369 for (i = 0; (i < gnx); i++)
374 read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE, fnm), &fr, TRX_NEED_V);
389 if (nframes >= n_alloc)
392 for (i = 0; i < gnx; i++)
394 srenew(c1[i], n_alloc);
397 for (i = 0; i < gnx; i += DIM)
399 c1[i + XX][nframes] = fr.v[index[i / DIM]][XX];
400 c1[i + YY][nframes] = fr.v[index[i / DIM]][YY];
401 c1[i + ZZ][nframes] = fr.v[index[i / DIM]][ZZ];
407 } while (read_next_frame(oenv, status, &fr));
411 if (nframes < min_frames)
413 gmx_fatal(FARGS, "You need at least %d frames in the trajectory and you only have %d.", min_frames, nframes);
415 dt = (t1 - t0) / (nframes - 1);
422 printf("Going to do %d fourier transforms of length %d. Hang on.\n", gnx, nframes);
424 /* Unfortunately the -normalize program option for the autocorrelation
425 * function calculation is added as a hack with a static variable in the
426 * autocorrelation.c source. That would work if we called the normal
427 * do_autocorr(), but this routine overrides that by directly calling
428 * the low-level functionality. That unfortunately leads to ignoring the
429 * default value for the option (which is to normalize).
430 * Since the absolute value seems to be important for the subsequent
431 * analysis below, we detect the value directly from the option, calculate
432 * the autocorrelation without normalization, and then apply the
433 * normalization just to the autocorrelation output
434 * (or not, if the user asked for a non-normalized autocorrelation).
436 normalizeAutocorrelation = opt2parg_bool("-normalize", npargs, ppa);
438 /* Note that we always disable normalization here, regardless of user settings */
440 nullptr, oenv, nullptr, nframes, gnx, nframes, c1, dt, eacNormal, 0, FALSE, FALSE, FALSE, -1, -1, 0);
442 for (j = 0; (j < DOS_NR); j++)
444 snew(dos[j], nframes + 4);
449 printf("Going to merge the ACFs into the mass-weighted and plain ACF\n");
451 for (i = 0; (i < gnx); i += DIM)
453 mi = top.atoms.atom[index[i / DIM]].m;
454 for (j = 0; (j < nframes / 2); j++)
456 c1j = (c1[i + XX][j] + c1[i + YY][j] + c1[i + ZZ][j]);
457 dos[VACF][j] += c1j / Natom;
458 dos[MVACF][j] += mi * c1j;
463 opt2fn("-vacf", NFILE, fnm), "Velocity autocorrelation function", "Time (ps)", "C(t)", oenv);
464 snew(tt, nframes / 2);
466 invNormalize = normalizeAutocorrelation ? 1.0 / dos[VACF][0] : 1.0;
468 for (j = 0; (j < nframes / 2); j++)
471 fprintf(fp, "%10g %10g\n", tt[j], dos[VACF][j] * invNormalize);
475 fp = xvgropen(opt2fn("-mvacf", NFILE, fnm),
476 "Mass-weighted velocity autocorrelation function",
481 invNormalize = normalizeAutocorrelation ? 1.0 / dos[VACF][0] : 1.0;
483 for (j = 0; (j < nframes / 2); j++)
485 fprintf(fp, "%10g %10g\n", tt[j], dos[MVACF][j] * invNormalize);
489 if ((fftcode = gmx_fft_init_1d_real(&fft, nframes / 2, GMX_FFT_FLAG_NONE)) != 0)
491 gmx_fatal(FARGS, "gmx_fft_init_1d_real returned %d", fftcode);
493 if ((fftcode = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX, dos[MVACF], dos[DOS])) != 0)
495 gmx_fatal(FARGS, "gmx_fft_1d_real returned %d", fftcode);
498 /* First compute the DoS */
499 /* Magic factor of 8 included now. */
500 bfac = 8 * dt * beta / 2;
502 snew(nu, nframes / 4);
503 for (j = 0; (j < nframes / 4); j++)
505 nu[j] = 2 * j / (t1 - t0);
506 dos2 += gmx::square(dos[DOS][2 * j]) + gmx::square(dos[DOS][2 * j + 1]);
509 dos[DOS][j] = bfac * std::hypot(dos[DOS][2 * j], dos[DOS][2 * j + 1]);
513 dos[DOS][j] = bfac * dos[DOS][2 * j];
517 dostot = evaluate_integral(nframes / 4, nu, dos[DOS], nullptr, int{ nframes / 4 }, &stddev);
520 for (j = 0; (j < nframes / 4); j++)
522 dos[DOS][j] *= 3 * Natom / dostot;
529 /* Note this eqn. is incorrect in Pascal2011a! */
530 Delta = ((2 * DoS0 / (9 * Natom)) * std::sqrt(M_PI * BOLTZ * Temp * Natom / tmass)
531 * std::pow((Natom / V), 1.0 / 3.0) * std::pow(6.0 / M_PI, 2.0 / 3.0));
532 f = calc_fluidicity(Delta, toler);
533 y = calc_y(f, Delta, toler);
534 z = calc_compress(y);
536 * (5.0 / 2.0 + std::log(2 * M_PI * BOLTZ * Temp / (gmx::square(PLANCK)) * V / (f * Natom)));
537 Shs = Sig + calc_Shs(f, y);
538 rho = (tmass * AMU) / (V * NANO * NANO * NANO);
539 sigHS = std::cbrt(6 * y * V / (M_PI * Natom));
541 fprintf(fplog, "System = \"%s\"\n", *top.name);
542 fprintf(fplog, "Nmol = %d\n", Nmol);
543 fprintf(fplog, "Natom = %d\n", Natom);
544 fprintf(fplog, "dt = %g ps\n", dt);
545 fprintf(fplog, "tmass = %g amu\n", tmass);
546 fprintf(fplog, "V = %g nm^3\n", V);
547 fprintf(fplog, "rho = %g g/l\n", rho);
548 fprintf(fplog, "T = %g K\n", Temp);
549 fprintf(fplog, "beta = %g mol/kJ\n", beta);
551 fprintf(fplog, "\nDoS parameters\n");
552 fprintf(fplog, "Delta = %g\n", Delta);
553 fprintf(fplog, "fluidicity = %g\n", f);
554 fprintf(fplog, "hard sphere packing fraction = %g\n", y);
555 fprintf(fplog, "hard sphere compressibility = %g\n", z);
556 fprintf(fplog, "ideal gas entropy = %g\n", Sig);
557 fprintf(fplog, "hard sphere entropy = %g\n", Shs);
558 fprintf(fplog, "sigma_HS = %g nm\n", sigHS);
559 fprintf(fplog, "DoS0 = %g\n", DoS0);
560 fprintf(fplog, "Dos2 = %g\n", dos2);
561 fprintf(fplog, "DoSTot = %g\n", dostot);
563 /* Now compute solid (2) and diffusive (3) components */
564 fp = xvgropen(opt2fn("-dos", NFILE, fnm),
566 bRecip ? "E (cm\\S-1\\N)" : "\\f{12}n\\f{4} (1/ps)",
567 "\\f{4}S(\\f{12}n\\f{4})",
569 xvgr_legend(fp, asize(DoSlegend), DoSlegend, oenv);
570 recip_fac = bRecip ? (1e7 / SPEED_OF_LIGHT) : 1.0;
571 for (j = 0; (j < nframes / 4); j++)
573 dos[DOS_DIFF][j] = DoS0 / (1 + gmx::square(DoS0 * M_PI * nu[j] / (6 * f * Natom)));
574 dos[DOS_SOLID][j] = dos[DOS][j] - dos[DOS_DIFF][j];
576 "%10g %10g %10g %10g\n",
578 dos[DOS][j] / recip_fac,
579 dos[DOS_SOLID][j] / recip_fac,
580 dos[DOS_DIFF][j] / recip_fac);
584 /* Finally analyze the results! */
586 wSdiff = Shs / (3 * BOLTZ); /* Is this correct? */
588 wAdiff = wEdiff - wSdiff;
589 for (j = 0; (j < nframes / 4); j++)
591 dos[DOS_CP][j] = (dos[DOS_DIFF][j] * wCdiff + dos[DOS_SOLID][j] * wCsolid(nu[j], beta));
592 dos[DOS_S][j] = (dos[DOS_DIFF][j] * wSdiff + dos[DOS_SOLID][j] * wSsolid(nu[j], beta));
593 dos[DOS_A][j] = (dos[DOS_DIFF][j] * wAdiff + dos[DOS_SOLID][j] * wAsolid(nu[j], beta));
594 dos[DOS_E][j] = (dos[DOS_DIFF][j] * wEdiff + dos[DOS_SOLID][j] * wEsolid(nu[j], beta));
596 DiffCoeff = evaluate_integral(nframes / 2, tt, dos[VACF], nullptr, nframes / 2., &stddev);
597 DiffCoeff = 1000 * DiffCoeff / 3.0;
598 fprintf(fplog, "Diffusion coefficient from VACF %g 10^-5 cm^2/s\n", DiffCoeff);
599 fprintf(fplog, "Diffusion coefficient from DoS %g 10^-5 cm^2/s\n", 1000 * DoS0 / (12 * tmass * beta));
601 cP = BOLTZ * evaluate_integral(nframes / 4, nu, dos[DOS_CP], nullptr, int{ nframes / 4 }, &stddev);
602 fprintf(fplog, "Heat capacity %g J/mol K\n", 1000 * cP / Nmol);
603 fprintf(fplog, "\nArrivederci!\n");
604 gmx_fio_fclose(fplog);
606 do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");