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