src/gromacs/gmxana/gmx_velacc.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
  11  *
  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.
  16  *
  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.
  21  *
  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.
  26  *
  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.
  34  *
  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.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include <cmath>
  41 #include <cstdio>
  42 #include <cstring>
  43
  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"
  65
  66 static void index_atom2mol(int* n, int* index, const t_block* mols)
  67 {
  68     int nat, i, nmol, mol, j;
  69
  70     nat  = *n;
  71     i    = 0;
  72     nmol = 0;
  73     mol  = 0;
  74     while (i < nat)
  75     {
  76         while (index[i] > mols->index[mol])
  77         {
  78             mol++;
  79             if (mol >= mols->nr)
  80             {
  81                 gmx_fatal(FARGS, "Atom index out of range: %d", index[i] + 1);
  82             }
  83         }
  84         for (j = mols->index[mol]; j < mols->index[mol + 1]; j++)
  85         {
  86             if (i >= nat || index[i] != j)
  87             {
  88                 gmx_fatal(FARGS, "The index group does not consist of whole molecules");
  89             }
  90             i++;
  91         }
  92         index[nmol++] = mol;
  93     }
  94
  95     fprintf(stderr, "\nSplit group of %d atoms into %d molecules\n", nat, nmol);
  96
  97     *n = nmol;
  98 }
  99
 100 static void precalc(const t_topology& top, real normm[])
 101 {
 102
 103     real mtot;
 104     int  i, j, k, l;
 105
 106     for (i = 0; i < top.mols.nr; i++)
 107     {
 108         k    = top.mols.index[i];
 109         l    = top.mols.index[i + 1];
 110         mtot = 0.0;
 111
 112         for (j = k; j < l; j++)
 113         {
 114             mtot += top.atoms.atom[j].m;
 115         }
 116
 117         for (j = k; j < l; j++)
 118         {
 119             normm[j] = top.atoms.atom[j].m / mtot;
 120         }
 121     }
 122 }
 123
 124 static void calc_spectrum(int n, const real c[], real dt, const char* fn, gmx_output_env_t* oenv, gmx_bool bRecip)
 125 {
 126     FILE*     fp;
 127     gmx_fft_t fft;
 128     int       i, status;
 129     real*     data;
 130     real      nu, omega, recip_fac;
 131
 132     snew(data, n * 2);
 133     for (i = 0; (i < n); i++)
 134     {
 135         data[i] = c[i];
 136     }
 137
 138     if ((status = gmx_fft_init_1d_real(&fft, n, GMX_FFT_FLAG_NONE)) != 0)
 139     {
 140         gmx_fatal(FARGS, "Invalid fft return status %d", status);
 141     }
 142     if ((status = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX, data, data)) != 0)
 143     {
 144         gmx_fatal(FARGS, "Invalid fft return status %d", status);
 145     }
 146     fp = xvgropen(fn,
 147                   "Vibrational Power Spectrum",
 148                   bRecip ? "\\f{12}w\\f{4} (cm\\S-1\\N)" : "\\f{12}n\\f{4} (ps\\S-1\\N)",
 149                   "a.u.",
 150                   oenv);
 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
 163      */
 164     recip_fac = bRecip ? (1e7 / SPEED_OF_LIGHT) : 1.0;
 165     for (i = 0; (i < n); i += 2)
 166     {
 167         nu    = i / (2 * dt);
 168         omega = nu * recip_fac;
 169         /* Computing the square magnitude of a complex number, since this is a power
 170          * spectrum.
 171          */
 172         fprintf(fp, "%10g  %10g\n", omega, gmx::square(data[i]) + gmx::square(data[i + 1]));
 173     }
 174     xvgrclose(fp);
 175     gmx_fft_destroy(fft);
 176     sfree(data);
 177 }
 178
 179 int gmx_velacc(int argc, char* argv[])
 180 {
 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." };
 194
 195     static gmx_bool bMass = FALSE, bMol = FALSE, bRecip = TRUE;
 196     t_pargs         pa[] = {
 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" }
 200     };
 201
 202     t_topology top;
 203     PbcType    pbcType = PbcType::Unset;
 204     t_trxframe fr;
 205     matrix     box;
 206     gmx_bool   bTPS = FALSE, bTop = FALSE;
 207     int        gnx;
 208     int*       index;
 209     char*      grpname;
 210     /* t0, t1 are the beginning and end time respectively.
 211      * dt is the time step, mass is temp variable for atomic mass.
 212      */
 213     real         t0, t1, dt, mass;
 214     t_trxstatus* status;
 215     int          counter, n_alloc, i, j, counter_dim, k, l;
 216     rvec         mv_mol;
 217     /* Array for the correlation function */
 218     real**            c1;
 219     real*             normm = nullptr;
 220     gmx_output_env_t* oenv;
 221
 222 #define NHISTO 360
 223
 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)
 230     int      npargs;
 231     t_pargs* ppa;
 232
 233     npargs = asize(pa);
 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))
 237     {
 238         sfree(ppa);
 239         return 0;
 240     }
 241
 242     if (bMol || bMass)
 243     {
 244         bTPS = ftp2bSet(efTPS, NFILE, fnm) || !ftp2bSet(efNDX, NFILE, fnm);
 245     }
 246
 247     if (bTPS)
 248     {
 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);
 251     }
 252     else
 253     {
 254         rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &gnx, &index, &grpname);
 255     }
 256
 257     if (bMol)
 258     {
 259         if (!bTop)
 260         {
 261             gmx_fatal(FARGS, "Need a topology to determine the molecules");
 262         }
 263         snew(normm, top.atoms.nr);
 264         precalc(top, normm);
 265         index_atom2mol(&gnx, index, &top.mols);
 266     }
 267
 268     /* Correlation stuff */
 269     snew(c1, gnx);
 270     for (i = 0; (i < gnx); i++)
 271     {
 272         c1[i] = nullptr;
 273     }
 274
 275     read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE, fnm), &fr, TRX_NEED_V);
 276     t0 = fr.time;
 277
 278     n_alloc = 0;
 279     counter = 0;
 280     do
 281     {
 282         if (counter >= n_alloc)
 283         {
 284             n_alloc += 100;
 285             for (i = 0; i < gnx; i++)
 286             {
 287                 srenew(c1[i], DIM * n_alloc);
 288             }
 289         }
 290         counter_dim = DIM * counter;
 291         if (bMol)
 292         {
 293             for (i = 0; i < gnx; i++)
 294             {
 295                 clear_rvec(mv_mol);
 296                 k = top.mols.index[index[i]];
 297                 l = top.mols.index[index[i] + 1];
 298                 for (j = k; j < l; j++)
 299                 {
 300                     if (bMass)
 301                     {
 302                         mass = top.atoms.atom[j].m;
 303                     }
 304                     else
 305                     {
 306                         mass = normm[j];
 307                     }
 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];
 311                 }
 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];
 315             }
 316         }
 317         else
 318         {
 319             for (i = 0; i < gnx; i++)
 320             {
 321                 if (bMass)
 322                 {
 323                     mass = top.atoms.atom[index[i]].m;
 324                 }
 325                 else
 326                 {
 327                     mass = 1;
 328                 }
 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];
 332             }
 333         }
 334
 335         t1 = fr.time;
 336
 337         counter++;
 338     } while (read_next_frame(oenv, status, &fr));
 339
 340     close_trx(status);
 341
 342     if (counter >= 4)
 343     {
 344         /* Compute time step between frames */
 345         dt = (t1 - t0) / (counter - 1);
 346         do_autocorr(opt2fn("-o", NFILE, fnm),
 347                     oenv,
 348                     bMass ? "Momentum Autocorrelation Function" : "Velocity Autocorrelation Function",
 349                     counter,
 350                     gnx,
 351                     c1,
 352                     dt,
 353                     eacVector,
 354                     TRUE);
 355
 356         do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
 357
 358         if (opt2bSet("-os", NFILE, fnm))
 359         {
 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");
 362         }
 363     }
 364     else
 365     {
 366         fprintf(stderr, "Not enough frames in trajectory - no output generated.\n");
 367     }
 368
 369     return 0;
 370 }