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, "Vibrational Power Spectrum",
 147                   bRecip ? "\\f{12}w\\f{4} (cm\\S-1\\N)" : "\\f{12}n\\f{4} (ps\\S-1\\N)", "a.u.", oenv);
 148     /* This is difficult.
 149      * The length of the ACF is dt (as passed to this routine).
 150      * We pass the vacf with N time steps from 0 to dt.
 151      * That means that after FFT we have lowest frequency = 1/dt
 152      * then 1/(2 dt) etc. (this is the X-axis of the data after FFT).
 153      * To convert to 1/cm we need to have to realize that
 154      * E = hbar w = h nu = h c/lambda. We want to have reciprokal cm
 155      * on the x-axis, that is 1/lambda, so we then have
 156      * 1/lambda = nu/c. Since nu has units of 1/ps and c has gromacs units
 157      * of nm/ps, we need to multiply by 1e7.
 158      * The timestep between saving the trajectory is
 159      * 1e7 is to convert nanometer to cm
 160      */
 161     recip_fac = bRecip ? (1e7 / SPEED_OF_LIGHT) : 1.0;
 162     for (i = 0; (i < n); i += 2)
 163     {
 164         nu    = i / (2 * dt);
 165         omega = nu * recip_fac;
 166         /* Computing the square magnitude of a complex number, since this is a power
 167          * spectrum.
 168          */
 169         fprintf(fp, "%10g  %10g\n", omega, gmx::square(data[i]) + gmx::square(data[i + 1]));
 170     }
 171     xvgrclose(fp);
 172     gmx_fft_destroy(fft);
 173     sfree(data);
 174 }
 175
 176 int gmx_velacc(int argc, char* argv[])
 177 {
 178     const char* desc[] = { "[THISMODULE] computes the velocity autocorrelation function.",
 179                            "When the [TT]-m[tt] option is used, the momentum autocorrelation",
 180                            "function is calculated.[PAR]",
 181                            "With option [TT]-mol[tt] the velocity autocorrelation function of",
 182                            "molecules is calculated. In this case the index group should consist",
 183                            "of molecule numbers instead of atom numbers.[PAR]",
 184                            "By using option [TT]-os[tt] you can also extract the estimated",
 185                            "(vibrational) power spectrum, which is the Fourier transform of the",
 186                            "velocity autocorrelation function.",
 187                            "Be sure that your trajectory contains frames with velocity information",
 188                            "(i.e. [TT]nstvout[tt] was set in your original [REF].mdp[ref] file),",
 189                            "and that the time interval between data collection points is",
 190                            "much shorter than the time scale of the autocorrelation." };
 191
 192     static gmx_bool bMass = FALSE, bMol = FALSE, bRecip = TRUE;
 193     t_pargs         pa[] = {
 194         { "-m", FALSE, etBOOL, { &bMass }, "Calculate the momentum autocorrelation function" },
 195         { "-recip", FALSE, etBOOL, { &bRecip }, "Use cm^-1 on X-axis instead of 1/ps for spectra." },
 196         { "-mol", FALSE, etBOOL, { &bMol }, "Calculate the velocity acf of molecules" }
 197     };
 198
 199     t_topology top;
 200     PbcType    pbcType = PbcType::Unset;
 201     t_trxframe fr;
 202     matrix     box;
 203     gmx_bool   bTPS = FALSE, bTop = FALSE;
 204     int        gnx;
 205     int*       index;
 206     char*      grpname;
 207     /* t0, t1 are the beginning and end time respectively.
 208      * dt is the time step, mass is temp variable for atomic mass.
 209      */
 210     real         t0, t1, dt, mass;
 211     t_trxstatus* status;
 212     int          counter, n_alloc, i, j, counter_dim, k, l;
 213     rvec         mv_mol;
 214     /* Array for the correlation function */
 215     real**            c1;
 216     real*             normm = nullptr;
 217     gmx_output_env_t* oenv;
 218
 219 #define NHISTO 360
 220
 221     t_filenm fnm[] = { { efTRN, "-f", nullptr, ffREAD },
 222                        { efTPS, nullptr, nullptr, ffOPTRD },
 223                        { efNDX, nullptr, nullptr, ffOPTRD },
 224                        { efXVG, "-o", "vac", ffWRITE },
 225                        { efXVG, "-os", "spectrum", ffOPTWR } };
 226 #define NFILE asize(fnm)
 227     int      npargs;
 228     t_pargs* ppa;
 229
 230     npargs = asize(pa);
 231     ppa    = add_acf_pargs(&npargs, pa);
 232     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, npargs, ppa,
 233                            asize(desc), desc, 0, nullptr, &oenv))
 234     {
 235         sfree(ppa);
 236         return 0;
 237     }
 238
 239     if (bMol || bMass)
 240     {
 241         bTPS = ftp2bSet(efTPS, NFILE, fnm) || !ftp2bSet(efNDX, NFILE, fnm);
 242     }
 243
 244     if (bTPS)
 245     {
 246         bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, nullptr, nullptr, box, TRUE);
 247         get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &gnx, &index, &grpname);
 248     }
 249     else
 250     {
 251         rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &gnx, &index, &grpname);
 252     }
 253
 254     if (bMol)
 255     {
 256         if (!bTop)
 257         {
 258             gmx_fatal(FARGS, "Need a topology to determine the molecules");
 259         }
 260         snew(normm, top.atoms.nr);
 261         precalc(top, normm);
 262         index_atom2mol(&gnx, index, &top.mols);
 263     }
 264
 265     /* Correlation stuff */
 266     snew(c1, gnx);
 267     for (i = 0; (i < gnx); i++)
 268     {
 269         c1[i] = nullptr;
 270     }
 271
 272     read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE, fnm), &fr, TRX_NEED_V);
 273     t0 = fr.time;
 274
 275     n_alloc = 0;
 276     counter = 0;
 277     do
 278     {
 279         if (counter >= n_alloc)
 280         {
 281             n_alloc += 100;
 282             for (i = 0; i < gnx; i++)
 283             {
 284                 srenew(c1[i], DIM * n_alloc);
 285             }
 286         }
 287         counter_dim = DIM * counter;
 288         if (bMol)
 289         {
 290             for (i = 0; i < gnx; i++)
 291             {
 292                 clear_rvec(mv_mol);
 293                 k = top.mols.index[index[i]];
 294                 l = top.mols.index[index[i] + 1];
 295                 for (j = k; j < l; j++)
 296                 {
 297                     if (bMass)
 298                     {
 299                         mass = top.atoms.atom[j].m;
 300                     }
 301                     else
 302                     {
 303                         mass = normm[j];
 304                     }
 305                     mv_mol[XX] += mass * fr.v[j][XX];
 306                     mv_mol[YY] += mass * fr.v[j][YY];
 307                     mv_mol[ZZ] += mass * fr.v[j][ZZ];
 308                 }
 309                 c1[i][counter_dim + XX] = mv_mol[XX];
 310                 c1[i][counter_dim + YY] = mv_mol[YY];
 311                 c1[i][counter_dim + ZZ] = mv_mol[ZZ];
 312             }
 313         }
 314         else
 315         {
 316             for (i = 0; i < gnx; i++)
 317             {
 318                 if (bMass)
 319                 {
 320                     mass = top.atoms.atom[index[i]].m;
 321                 }
 322                 else
 323                 {
 324                     mass = 1;
 325                 }
 326                 c1[i][counter_dim + XX] = mass * fr.v[index[i]][XX];
 327                 c1[i][counter_dim + YY] = mass * fr.v[index[i]][YY];
 328                 c1[i][counter_dim + ZZ] = mass * fr.v[index[i]][ZZ];
 329             }
 330         }
 331
 332         t1 = fr.time;
 333
 334         counter++;
 335     } while (read_next_frame(oenv, status, &fr));
 336
 337     close_trx(status);
 338
 339     if (counter >= 4)
 340     {
 341         /* Compute time step between frames */
 342         dt = (t1 - t0) / (counter - 1);
 343         do_autocorr(opt2fn("-o", NFILE, fnm), oenv,
 344                     bMass ? "Momentum Autocorrelation Function" : "Velocity Autocorrelation Function",
 345                     counter, gnx, c1, dt, eacVector, TRUE);
 346
 347         do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
 348
 349         if (opt2bSet("-os", NFILE, fnm))
 350         {
 351             calc_spectrum(counter / 2, (c1[0]), (t1 - t0) / 2, opt2fn("-os", NFILE, fnm), oenv, bRecip);
 352             do_view(oenv, opt2fn("-os", NFILE, fnm), "-nxy");
 353         }
 354     }
 355     else
 356     {
 357         fprintf(stderr, "Not enough frames in trajectory - no output generated.\n");
 358     }
 359
 360     return 0;
 361 }