2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2011,2012,2013,2014,2015, 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/correlationfunctions/autocorr.h"
44 #include "gromacs/correlationfunctions/integrate.h"
45 #include "gromacs/fft/fft.h"
46 #include "gromacs/fileio/confio.h"
47 #include "gromacs/fileio/gmxfio.h"
48 #include "gromacs/fileio/trxio.h"
49 #include "gromacs/fileio/xvgr.h"
50 #include "gromacs/gmxana/gmx_ana.h"
51 #include "gromacs/legacyheaders/copyrite.h"
52 #include "gromacs/legacyheaders/macros.h"
53 #include "gromacs/legacyheaders/txtdump.h"
54 #include "gromacs/legacyheaders/typedefs.h"
55 #include "gromacs/legacyheaders/viewit.h"
56 #include "gromacs/math/units.h"
57 #include "gromacs/math/utilities.h"
58 #include "gromacs/math/vec.h"
59 #include "gromacs/topology/index.h"
60 #include "gromacs/utility/fatalerror.h"
61 #include "gromacs/utility/futil.h"
62 #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(t_block *mols, int natoms, atom_id *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*pow(Delta, -4.5)*pow(f, 7.5) -
106 6*pow(Delta, -3)*pow(f, 5) -
107 pow(Delta, -1.5)*pow(f, 3.5) +
108 6*pow(Delta, -1.5)*pow(f, 2.5) +
112 static double YYY(double f, double y)
114 return (2*pow(y*f, 3) - sqr(f)*y*(1+6*y) +
118 static double calc_compress(double y)
124 return ((1+y+sqr(y)-pow(y, 3))/(pow(1-y, 3)));
127 static double bisector(double Delta, double tol,
128 double ff0, double ff1,
129 double ff(double, double))
131 double fd0, fd, fd1, f, f0, f1;
132 double tolmin = 1e-8;
138 fprintf(stderr, "Unrealistic tolerance %g for bisector. Setting it to %g\n", tol, tolmin);
161 while ((f1-f0) > tol);
166 static double calc_fluidicity(double Delta, double tol)
168 return bisector(Delta, tol, 0, 1, FD);
171 static double calc_y(double f, double Delta, double toler)
175 y1 = pow(f/Delta, 1.5);
176 y2 = bisector(f, toler, 0, 10000, YYY);
177 if (fabs((y1-y2)/(y1+y2)) > 100*toler)
179 fprintf(stderr, "Inconsistency computing y: y1 = %f, y2 = %f, using y1.\n",
186 static double calc_Shs(double f, double y)
190 return BOLTZ*(log(calc_compress(fy)) + fy*(3*fy-4)/sqr(1-fy));
193 static real wCsolid(real nu, real beta)
195 real bhn = beta*PLANCK*nu;
206 return sqr(bhn)*ebn/koko;
210 static real wSsolid(real nu, real beta)
212 real bhn = beta*PLANCK*nu;
220 return bhn/gmx_expm1(bhn) - gmx_log1p(-exp(-bhn));
224 static real wAsolid(real nu, real beta)
226 real bhn = beta*PLANCK*nu;
234 return log((1-exp(-bhn))/(exp(-bhn/2))) - log(bhn);
238 static real wEsolid(real nu, real beta)
240 real bhn = beta*PLANCK*nu;
248 return bhn/2 + bhn/gmx_expm1(bhn)-1;
252 int gmx_dos(int argc, char *argv[])
254 const char *desc[] = {
255 "[THISMODULE] computes the Density of States from a simulations.",
256 "In order for this to be meaningful the velocities must be saved",
257 "in the trajecotry with sufficiently high frequency such as to cover",
258 "all vibrations. For flexible systems that would be around a few fs",
259 "between saving. Properties based on the DoS are printed on the",
261 "Note that the density of states is calculated from the mass-weighted",
262 "autocorrelation, and by default only from the square of the real",
263 "component rather than absolute value. This means the shape can differ",
264 "substantially from the plain vibrational power spectrum you can",
265 "calculate with gmx velacc."
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;
290 gmx_bool normalizeAutocorrelation;
292 static gmx_bool bVerbose = TRUE, bAbsolute = FALSE, bNormalizeDos = FALSE;
293 static gmx_bool bRecip = FALSE;
294 static real Temp = 298.15, toler = 1e-6;
297 { "-v", FALSE, etBOOL, {&bVerbose},
298 "Be loud and noisy." },
299 { "-recip", FALSE, etBOOL, {&bRecip},
300 "Use cm^-1 on X-axis instead of 1/ps for DoS plots." },
301 { "-abs", FALSE, etBOOL, {&bAbsolute},
302 "Use the absolute value of the Fourier transform of the VACF as the Density of States. Default is to use the real component only" },
303 { "-normdos", FALSE, etBOOL, {&bNormalizeDos},
304 "Normalize the DoS such that it adds up to 3N. This should usually not be necessary." },
305 { "-T", FALSE, etREAL, {&Temp},
306 "Temperature in the simulation" },
307 { "-toler", FALSE, etREAL, {&toler},
308 "[HIDDEN]Tolerance when computing the fluidicity using bisection algorithm" }
312 { efTRN, "-f", NULL, ffREAD },
313 { efTPR, "-s", NULL, ffREAD },
314 { efNDX, NULL, NULL, ffOPTRD },
315 { efXVG, "-vacf", "vacf", ffWRITE },
316 { efXVG, "-mvacf", "mvacf", ffWRITE },
317 { efXVG, "-dos", "dos", ffWRITE },
318 { efLOG, "-g", "dos", ffWRITE },
320 #define NFILE asize(fnm)
323 const char *DoSlegend[] = {
324 "DoS(v)", "DoS(v)[Solid]", "DoS(v)[Diff]"
328 ppa = add_acf_pargs(&npargs, pa);
329 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
330 NFILE, fnm, npargs, ppa, asize(desc), desc,
331 asize(bugs), bugs, &oenv))
336 beta = 1/(Temp*BOLTZ);
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(efTPR, NFILE, fnm), title, &top, &ePBC, NULL, NULL, box, TRUE);
345 /* Handle index groups */
346 get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &grpNatoms, &index, &grpname);
350 for (i = 0; i < grpNatoms; i++)
352 tmass += top.atoms.atom[index[i]].m;
356 Nmol = calcMoleculesInIndexGroup(&top.mols, top.atoms.nr, index, grpNatoms);
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[index[i/DIM]][XX];
393 c1[i+YY][nframes] = fr.v[index[i/DIM]][YY];
394 c1[i+ZZ][nframes] = fr.v[index[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 /* Unfortunately the -normalize program option for the autocorrelation
416 * function calculation is added as a hack with a static variable in the
417 * autocorrelation.c source. That would work if we called the normal
418 * do_autocorr(), but this routine overrides that by directly calling
419 * the low-level functionality. That unfortunately leads to ignoring the
420 * default value for the option (which is to normalize).
421 * Since the absolute value seems to be important for the subsequent
422 * analysis below, we detect the value directly from the option, calculate
423 * the autocorrelation without normalization, and then apply the
424 * normalization just to the autocorrelation output
425 * (or not, if the user asked for a non-normalized autocorrelation).
427 normalizeAutocorrelation = opt2parg_bool("-normalize", npargs, ppa);
429 /* Note that we always disable normalization here, regardless of user settings */
430 low_do_autocorr(NULL, oenv, NULL, nframes, gnx, nframes, c1, dt, eacNormal, 0, FALSE,
431 FALSE, FALSE, -1, -1, 0);
433 for (j = 0; (j < DOS_NR); j++)
435 snew(dos[j], nframes+4);
440 printf("Going to merge the ACFs into the mass-weighted and plain ACF\n");
442 for (i = 0; (i < gnx); i += DIM)
444 mi = top.atoms.atom[index[i/DIM]].m;
445 for (j = 0; (j < nframes/2); j++)
447 c1j = (c1[i+XX][j] + c1[i+YY][j] + c1[i+ZZ][j]);
448 dos[VACF][j] += c1j/Natom;
449 dos[MVACF][j] += mi*c1j;
453 fp = xvgropen(opt2fn("-vacf", NFILE, fnm), "Velocity autocorrelation function",
454 "Time (ps)", "C(t)", oenv);
457 invNormalize = normalizeAutocorrelation ? 1.0/dos[VACF][0] : 1.0;
459 for (j = 0; (j < nframes/2); j++)
462 fprintf(fp, "%10g %10g\n", tt[j], dos[VACF][j] * invNormalize);
466 fp = xvgropen(opt2fn("-mvacf", NFILE, fnm), "Mass-weighted velocity autocorrelation function",
467 "Time (ps)", "C(t)", oenv);
469 invNormalize = normalizeAutocorrelation ? 1.0/dos[VACF][0] : 1.0;
471 for (j = 0; (j < nframes/2); j++)
473 fprintf(fp, "%10g %10g\n", tt[j], dos[MVACF][j] * invNormalize);
477 if ((fftcode = gmx_fft_init_1d_real(&fft, nframes/2,
478 GMX_FFT_FLAG_NONE)) != 0)
480 gmx_fatal(FARGS, "gmx_fft_init_1d_real returned %d", fftcode);
482 if ((fftcode = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX,
483 (void *)dos[MVACF], (void *)dos[DOS])) != 0)
485 gmx_fatal(FARGS, "gmx_fft_1d_real returned %d", fftcode);
488 /* First compute the DoS */
489 /* Magic factor of 8 included now. */
493 for (j = 0; (j < nframes/4); j++)
496 dos2 += sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]);
499 dos[DOS][j] = bfac*sqrt(sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]));
503 dos[DOS][j] = bfac*dos[DOS][2*j];
507 dostot = evaluate_integral(nframes/4, nu, dos[DOS], NULL, nframes/4, &stddev);
510 for (j = 0; (j < nframes/4); j++)
512 dos[DOS][j] *= 3*Natom/dostot;
519 /* Note this eqn. is incorrect in Pascal2011a! */
520 Delta = ((2*DoS0/(9*Natom))*sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*
521 pow((Natom/V), 1.0/3.0)*pow(6/M_PI, 2.0/3.0));
522 f = calc_fluidicity(Delta, toler);
523 y = calc_y(f, Delta, toler);
524 z = calc_compress(y);
525 Sig = BOLTZ*(5.0/2.0+log(2*M_PI*BOLTZ*Temp/(sqr(PLANCK))*V/(f*Natom)));
526 Shs = Sig+calc_Shs(f, y);
527 rho = (tmass*AMU)/(V*NANO*NANO*NANO);
528 sigHS = pow(6*y*V/(M_PI*Natom), 1.0/3.0);
530 fprintf(fplog, "System = \"%s\"\n", title);
531 fprintf(fplog, "Nmol = %d\n", Nmol);
532 fprintf(fplog, "Natom = %d\n", Natom);
533 fprintf(fplog, "dt = %g ps\n", dt);
534 fprintf(fplog, "tmass = %g amu\n", tmass);
535 fprintf(fplog, "V = %g nm^3\n", V);
536 fprintf(fplog, "rho = %g g/l\n", rho);
537 fprintf(fplog, "T = %g K\n", Temp);
538 fprintf(fplog, "beta = %g mol/kJ\n", beta);
540 fprintf(fplog, "\nDoS parameters\n");
541 fprintf(fplog, "Delta = %g\n", Delta);
542 fprintf(fplog, "fluidicity = %g\n", f);
543 fprintf(fplog, "hard sphere packing fraction = %g\n", y);
544 fprintf(fplog, "hard sphere compressibility = %g\n", z);
545 fprintf(fplog, "ideal gas entropy = %g\n", Sig);
546 fprintf(fplog, "hard sphere entropy = %g\n", Shs);
547 fprintf(fplog, "sigma_HS = %g nm\n", sigHS);
548 fprintf(fplog, "DoS0 = %g\n", DoS0);
549 fprintf(fplog, "Dos2 = %g\n", dos2);
550 fprintf(fplog, "DoSTot = %g\n", dostot);
552 /* Now compute solid (2) and diffusive (3) components */
553 fp = xvgropen(opt2fn("-dos", NFILE, fnm), "Density of states",
554 bRecip ? "E (cm\\S-1\\N)" : "\\f{12}n\\f{4} (1/ps)",
555 "\\f{4}S(\\f{12}n\\f{4})", oenv);
556 xvgr_legend(fp, asize(DoSlegend), DoSlegend, oenv);
557 recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;
558 for (j = 0; (j < nframes/4); j++)
560 dos[DOS_DIFF][j] = DoS0/(1+sqr(DoS0*M_PI*nu[j]/(6*f*Natom)));
561 dos[DOS_SOLID][j] = dos[DOS][j]-dos[DOS_DIFF][j];
562 fprintf(fp, "%10g %10g %10g %10g\n",
564 dos[DOS][j]/recip_fac,
565 dos[DOS_SOLID][j]/recip_fac,
566 dos[DOS_DIFF][j]/recip_fac);
570 /* Finally analyze the results! */
572 wSdiff = Shs/(3*BOLTZ); /* Is this correct? */
574 wAdiff = wEdiff-wSdiff;
575 for (j = 0; (j < nframes/4); j++)
577 dos[DOS_CP][j] = (dos[DOS_DIFF][j]*wCdiff +
578 dos[DOS_SOLID][j]*wCsolid(nu[j], beta));
579 dos[DOS_S][j] = (dos[DOS_DIFF][j]*wSdiff +
580 dos[DOS_SOLID][j]*wSsolid(nu[j], beta));
581 dos[DOS_A][j] = (dos[DOS_DIFF][j]*wAdiff +
582 dos[DOS_SOLID][j]*wAsolid(nu[j], beta));
583 dos[DOS_E][j] = (dos[DOS_DIFF][j]*wEdiff +
584 dos[DOS_SOLID][j]*wEsolid(nu[j], beta));
586 DiffCoeff = evaluate_integral(nframes/2, tt, dos[VACF], NULL, nframes/2, &stddev);
587 DiffCoeff = 1000*DiffCoeff/3.0;
588 fprintf(fplog, "Diffusion coefficient from VACF %g 10^-5 cm^2/s\n",
590 fprintf(fplog, "Diffusion coefficient from DoS %g 10^-5 cm^2/s\n",
591 1000*DoS0/(12*tmass*beta));
593 cP = BOLTZ * evaluate_integral(nframes/4, nu, dos[DOS_CP], NULL,
595 fprintf(fplog, "Heat capacity %g J/mol K\n", 1000*cP/Nmol);
596 fprintf(fplog, "\nArrivederci!\n");
597 gmx_fio_fclose(fplog);
599 do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");