2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
44 #include "gromacs/commandline/pargs.h"
45 #include "gromacs/commandline/viewit.h"
46 #include "gromacs/correlationfunctions/autocorr.h"
47 #include "gromacs/fft/fft.h"
48 #include "gromacs/fileio/confio.h"
49 #include "gromacs/fileio/trxio.h"
50 #include "gromacs/fileio/xvgr.h"
51 #include "gromacs/gmxana/gmx_ana.h"
52 #include "gromacs/gmxana/gstat.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/smalloc.h"
66 static void index_atom2mol(int* n, int* index, const t_block* mols)
68 int nat, i, nmol, mol, j;
76 while (index[i] > mols->index[mol])
81 gmx_fatal(FARGS, "Atom index out of range: %d", index[i] + 1);
84 for (j = mols->index[mol]; j < mols->index[mol + 1]; j++)
86 if (i >= nat || index[i] != j)
88 gmx_fatal(FARGS, "The index group does not consist of whole molecules");
95 fprintf(stderr, "\nSplit group of %d atoms into %d molecules\n", nat, nmol);
100 static void precalc(const t_topology& top, real normm[])
106 for (i = 0; i < top.mols.nr; i++)
108 k = top.mols.index[i];
109 l = top.mols.index[i + 1];
112 for (j = k; j < l; j++)
114 mtot += top.atoms.atom[j].m;
117 for (j = k; j < l; j++)
119 normm[j] = top.atoms.atom[j].m / mtot;
124 static void calc_spectrum(int n, const real c[], real dt, const char* fn, gmx_output_env_t* oenv, gmx_bool bRecip)
130 real nu, omega, recip_fac;
133 for (i = 0; (i < n); i++)
138 if ((status = gmx_fft_init_1d_real(&fft, n, GMX_FFT_FLAG_NONE)) != 0)
140 gmx_fatal(FARGS, "Invalid fft return status %d", status);
142 if ((status = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX, data, data)) != 0)
144 gmx_fatal(FARGS, "Invalid fft return status %d", status);
147 "Vibrational Power Spectrum",
148 bRecip ? "\\f{12}w\\f{4} (cm\\S-1\\N)" : "\\f{12}n\\f{4} (ps\\S-1\\N)",
151 /* This is difficult.
152 * The length of the ACF is dt (as passed to this routine).
153 * We pass the vacf with N time steps from 0 to dt.
154 * That means that after FFT we have lowest frequency = 1/dt
155 * then 1/(2 dt) etc. (this is the X-axis of the data after FFT).
156 * To convert to 1/cm we need to have to realize that
157 * E = hbar w = h nu = h c/lambda. We want to have reciprokal cm
158 * on the x-axis, that is 1/lambda, so we then have
159 * 1/lambda = nu/c. Since nu has units of 1/ps and c has gromacs units
160 * of nm/ps, we need to multiply by 1e7.
161 * The timestep between saving the trajectory is
162 * 1e7 is to convert nanometer to cm
164 recip_fac = bRecip ? (1e7 / SPEED_OF_LIGHT) : 1.0;
165 for (i = 0; (i < n); i += 2)
168 omega = nu * recip_fac;
169 /* Computing the square magnitude of a complex number, since this is a power
172 fprintf(fp, "%10g %10g\n", omega, gmx::square(data[i]) + gmx::square(data[i + 1]));
175 gmx_fft_destroy(fft);
179 int gmx_velacc(int argc, char* argv[])
181 const char* desc[] = { "[THISMODULE] computes the velocity autocorrelation function.",
182 "When the [TT]-m[tt] option is used, the momentum autocorrelation",
183 "function is calculated.[PAR]",
184 "With option [TT]-mol[tt] the velocity autocorrelation function of",
185 "molecules is calculated. In this case the index group should consist",
186 "of molecule numbers instead of atom numbers.[PAR]",
187 "By using option [TT]-os[tt] you can also extract the estimated",
188 "(vibrational) power spectrum, which is the Fourier transform of the",
189 "velocity autocorrelation function.",
190 "Be sure that your trajectory contains frames with velocity information",
191 "(i.e. [TT]nstvout[tt] was set in your original [REF].mdp[ref] file),",
192 "and that the time interval between data collection points is",
193 "much shorter than the time scale of the autocorrelation." };
195 static gmx_bool bMass = FALSE, bMol = FALSE, bRecip = TRUE;
197 { "-m", FALSE, etBOOL, { &bMass }, "Calculate the momentum autocorrelation function" },
198 { "-recip", FALSE, etBOOL, { &bRecip }, "Use cm^-1 on X-axis instead of 1/ps for spectra." },
199 { "-mol", FALSE, etBOOL, { &bMol }, "Calculate the velocity acf of molecules" }
203 PbcType pbcType = PbcType::Unset;
206 gmx_bool bTPS = FALSE, bTop = FALSE;
210 /* t0, t1 are the beginning and end time respectively.
211 * dt is the time step, mass is temp variable for atomic mass.
213 real t0, t1, dt, mass;
215 int counter, n_alloc, i, j, counter_dim, k, l;
217 /* Array for the correlation function */
219 real* normm = nullptr;
220 gmx_output_env_t* oenv;
224 t_filenm fnm[] = { { efTRN, "-f", nullptr, ffREAD },
225 { efTPS, nullptr, nullptr, ffOPTRD },
226 { efNDX, nullptr, nullptr, ffOPTRD },
227 { efXVG, "-o", "vac", ffWRITE },
228 { efXVG, "-os", "spectrum", ffOPTWR } };
229 #define NFILE asize(fnm)
234 ppa = add_acf_pargs(&npargs, pa);
235 if (!parse_common_args(
236 &argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, npargs, ppa, asize(desc), desc, 0, nullptr, &oenv))
244 bTPS = ftp2bSet(efTPS, NFILE, fnm) || !ftp2bSet(efNDX, NFILE, fnm);
249 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, nullptr, nullptr, box, TRUE);
250 get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &gnx, &index, &grpname);
254 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &gnx, &index, &grpname);
261 gmx_fatal(FARGS, "Need a topology to determine the molecules");
263 snew(normm, top.atoms.nr);
265 index_atom2mol(&gnx, index, &top.mols);
268 /* Correlation stuff */
270 for (i = 0; (i < gnx); i++)
275 read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE, fnm), &fr, TRX_NEED_V);
282 if (counter >= n_alloc)
285 for (i = 0; i < gnx; i++)
287 srenew(c1[i], DIM * n_alloc);
290 counter_dim = DIM * counter;
293 for (i = 0; i < gnx; i++)
296 k = top.mols.index[index[i]];
297 l = top.mols.index[index[i] + 1];
298 for (j = k; j < l; j++)
302 mass = top.atoms.atom[j].m;
308 mv_mol[XX] += mass * fr.v[j][XX];
309 mv_mol[YY] += mass * fr.v[j][YY];
310 mv_mol[ZZ] += mass * fr.v[j][ZZ];
312 c1[i][counter_dim + XX] = mv_mol[XX];
313 c1[i][counter_dim + YY] = mv_mol[YY];
314 c1[i][counter_dim + ZZ] = mv_mol[ZZ];
319 for (i = 0; i < gnx; i++)
323 mass = top.atoms.atom[index[i]].m;
329 c1[i][counter_dim + XX] = mass * fr.v[index[i]][XX];
330 c1[i][counter_dim + YY] = mass * fr.v[index[i]][YY];
331 c1[i][counter_dim + ZZ] = mass * fr.v[index[i]][ZZ];
338 } while (read_next_frame(oenv, status, &fr));
344 /* Compute time step between frames */
345 dt = (t1 - t0) / (counter - 1);
346 do_autocorr(opt2fn("-o", NFILE, fnm),
348 bMass ? "Momentum Autocorrelation Function" : "Velocity Autocorrelation Function",
356 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
358 if (opt2bSet("-os", NFILE, fnm))
360 calc_spectrum(counter / 2, (c1[0]), (t1 - t0) / 2, opt2fn("-os", NFILE, fnm), oenv, bRecip);
361 do_view(oenv, opt2fn("-os", NFILE, fnm), "-nxy");
366 fprintf(stderr, "Not enough frames in trajectory - no output generated.\n");