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,2018, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include <cmath>
  40 #include <cstdio>
  41 #include <cstring>
  42
  43 #include "gromacs/commandline/pargs.h"
  44 #include "gromacs/commandline/viewit.h"
  45 #include "gromacs/correlationfunctions/autocorr.h"
  46 #include "gromacs/fft/fft.h"
  47 #include "gromacs/fileio/confio.h"
  48 #include "gromacs/fileio/trxio.h"
  49 #include "gromacs/fileio/xvgr.h"
  50 #include "gromacs/gmxana/gmx_ana.h"
  51 #include "gromacs/gmxana/gstat.h"
  52 #include "gromacs/math/functions.h"
  53 #include "gromacs/math/units.h"
  54 #include "gromacs/math/utilities.h"
  55 #include "gromacs/math/vec.h"
  56 #include "gromacs/topology/index.h"
  57 #include "gromacs/topology/topology.h"
  58 #include "gromacs/trajectory/trajectoryframe.h"
  59 #include "gromacs/utility/arraysize.h"
  60 #include "gromacs/utility/fatalerror.h"
  61 #include "gromacs/utility/futil.h"
  62 #include "gromacs/utility/smalloc.h"
  63
  64 static void index_atom2mol(int *n, int *index, const t_block *mols)
  65 {
  66     int nat, i, nmol, mol, j;
  67
  68     nat  = *n;
  69     i    = 0;
  70     nmol = 0;
  71     mol  = 0;
  72     while (i < nat)
  73     {
  74         while (index[i] > mols->index[mol])
  75         {
  76             mol++;
  77             if (mol >= mols->nr)
  78             {
  79                 gmx_fatal(FARGS, "Atom index out of range: %d", index[i]+1);
  80             }
  81         }
  82         for (j = mols->index[mol]; j < mols->index[mol+1]; j++)
  83         {
  84             if (i >= nat || index[i] != j)
  85             {
  86                 gmx_fatal(FARGS, "The index group does not consist of whole molecules");
  87             }
  88             i++;
  89         }
  90         index[nmol++] = mol;
  91     }
  92
  93     fprintf(stderr, "\nSplit group of %d atoms into %d molecules\n", nat, nmol);
  94
  95     *n = nmol;
  96 }
  97
  98 static void precalc(const t_topology &top, real normm[])
  99 {
 100
 101     real mtot;
 102     int  i, j, k, l;
 103
 104     for (i = 0; i < top.mols.nr; i++)
 105     {
 106         k    = top.mols.index[i];
 107         l    = top.mols.index[i+1];
 108         mtot = 0.0;
 109
 110         for (j = k; j < l; j++)
 111         {
 112             mtot += top.atoms.atom[j].m;
 113         }
 114
 115         for (j = k; j < l; j++)
 116         {
 117             normm[j] = top.atoms.atom[j].m/mtot;
 118         }
 119
 120     }
 121
 122 }
 123
 124 static void calc_spectrum(int n, const real c[], real dt, const char *fn,
 125                           gmx_output_env_t *oenv, gmx_bool bRecip)
 126 {
 127     FILE     *fp;
 128     gmx_fft_t fft;
 129     int       i, status;
 130     real     *data;
 131     real      nu, omega, recip_fac;
 132
 133     snew(data, n*2);
 134     for (i = 0; (i < n); i++)
 135     {
 136         data[i] = c[i];
 137     }
 138
 139     if ((status = gmx_fft_init_1d_real(&fft, n, GMX_FFT_FLAG_NONE)) != 0)
 140     {
 141         gmx_fatal(FARGS, "Invalid fft return status %d", status);
 142     }
 143     if ((status = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX, data, data)) != 0)
 144     {
 145         gmx_fatal(FARGS, "Invalid fft return status %d", status);
 146     }
 147     fp = xvgropen(fn, "Vibrational Power Spectrum",
 148                   bRecip ? "\\f{12}w\\f{4} (cm\\S-1\\N)" :
 149                   "\\f{12}n\\f{4} (ps\\S-1\\N)",
 150                   "a.u.", 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[] = {
 182         "[THISMODULE] computes the velocity autocorrelation function.",
 183         "When the [TT]-m[tt] option is used, the momentum autocorrelation",
 184         "function is calculated.[PAR]",
 185         "With option [TT]-mol[tt] the velocity autocorrelation function of",
 186         "molecules is calculated. In this case the index group should consist",
 187         "of molecule numbers instead of atom numbers.[PAR]",
 188         "By using option [TT]-os[tt] you can also extract the estimated",
 189         "(vibrational) power spectrum, which is the Fourier transform of the",
 190         "velocity autocorrelation function.",
 191         "Be sure that your trajectory contains frames with velocity information",
 192         "(i.e. [TT]nstvout[tt] was set in your original [REF].mdp[ref] file),",
 193         "and that the time interval between data collection points is",
 194         "much shorter than the time scale of the autocorrelation."
 195     };
 196
 197     static gmx_bool bMass = FALSE, bMol = FALSE, bRecip = TRUE;
 198     t_pargs         pa[]  = {
 199         { "-m", FALSE, etBOOL, {&bMass},
 200           "Calculate the momentum autocorrelation function" },
 201         { "-recip", FALSE, etBOOL, {&bRecip},
 202           "Use cm^-1 on X-axis instead of 1/ps for spectra." },
 203         { "-mol", FALSE, etBOOL, {&bMol},
 204           "Calculate the velocity acf of molecules" }
 205     };
 206
 207     t_topology      top;
 208     int             ePBC = -1;
 209     t_trxframe      fr;
 210     matrix          box;
 211     gmx_bool        bTPS = FALSE, bTop = FALSE;
 212     int             gnx;
 213     int            *index;
 214     char           *grpname;
 215     /* t0, t1 are the beginning and end time respectively.
 216      * dt is the time step, mass is temp variable for atomic mass.
 217      */
 218     real              t0, t1, dt, mass;
 219     t_trxstatus      *status;
 220     int               counter, n_alloc, i, j, counter_dim, k, l;
 221     rvec              mv_mol;
 222     /* Array for the correlation function */
 223     real            **c1;
 224     real             *normm = nullptr;
 225     gmx_output_env_t *oenv;
 226
 227 #define NHISTO 360
 228
 229     t_filenm  fnm[] = {
 230         { efTRN, "-f",    nullptr,   ffREAD  },
 231         { efTPS, nullptr,    nullptr,   ffOPTRD },
 232         { efNDX, nullptr,    nullptr,   ffOPTRD },
 233         { efXVG, "-o",    "vac",  ffWRITE },
 234         { efXVG, "-os",   "spectrum", ffOPTWR }
 235     };
 236 #define NFILE asize(fnm)
 237     int       npargs;
 238     t_pargs  *ppa;
 239
 240     npargs = asize(pa);
 241     ppa    = add_acf_pargs(&npargs, pa);
 242     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
 243                            NFILE, fnm, npargs, ppa, asize(desc), desc, 0, nullptr, &oenv))
 244     {
 245         sfree(ppa);
 246         return 0;
 247     }
 248
 249     if (bMol || bMass)
 250     {
 251         bTPS = ftp2bSet(efTPS, NFILE, fnm) || !ftp2bSet(efNDX, NFILE, fnm);
 252     }
 253
 254     if (bTPS)
 255     {
 256         bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, nullptr, nullptr, box,
 257                              TRUE);
 258         get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &gnx, &index, &grpname);
 259     }
 260     else
 261     {
 262         rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &gnx, &index, &grpname);
 263     }
 264
 265     if (bMol)
 266     {
 267         if (!bTop)
 268         {
 269             gmx_fatal(FARGS, "Need a topology to determine the molecules");
 270         }
 271         snew(normm, top.atoms.nr);
 272         precalc(top, normm);
 273         index_atom2mol(&gnx, index, &top.mols);
 274     }
 275
 276     /* Correlation stuff */
 277     snew(c1, gnx);
 278     for (i = 0; (i < gnx); i++)
 279     {
 280         c1[i] = nullptr;
 281     }
 282
 283     read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE, fnm), &fr, TRX_NEED_V);
 284     t0 = fr.time;
 285
 286     n_alloc = 0;
 287     counter = 0;
 288     do
 289     {
 290         if (counter >= n_alloc)
 291         {
 292             n_alloc += 100;
 293             for (i = 0; i < gnx; i++)
 294             {
 295                 srenew(c1[i], DIM*n_alloc);
 296             }
 297         }
 298         counter_dim = DIM*counter;
 299         if (bMol)
 300         {
 301             for (i = 0; i < gnx; i++)
 302             {
 303                 clear_rvec(mv_mol);
 304                 k = top.mols.index[index[i]];
 305                 l = top.mols.index[index[i]+1];
 306                 for (j = k; j < l; j++)
 307                 {
 308                     if (bMass)
 309                     {
 310                         mass = top.atoms.atom[j].m;
 311                     }
 312                     else
 313                     {
 314                         mass = normm[j];
 315                     }
 316                     mv_mol[XX] += mass*fr.v[j][XX];
 317                     mv_mol[YY] += mass*fr.v[j][YY];
 318                     mv_mol[ZZ] += mass*fr.v[j][ZZ];
 319                 }
 320                 c1[i][counter_dim+XX] = mv_mol[XX];
 321                 c1[i][counter_dim+YY] = mv_mol[YY];
 322                 c1[i][counter_dim+ZZ] = mv_mol[ZZ];
 323             }
 324         }
 325         else
 326         {
 327             for (i = 0; i < gnx; i++)
 328             {
 329                 if (bMass)
 330                 {
 331                     mass = top.atoms.atom[index[i]].m;
 332                 }
 333                 else
 334                 {
 335                     mass = 1;
 336                 }
 337                 c1[i][counter_dim+XX] = mass*fr.v[index[i]][XX];
 338                 c1[i][counter_dim+YY] = mass*fr.v[index[i]][YY];
 339                 c1[i][counter_dim+ZZ] = mass*fr.v[index[i]][ZZ];
 340             }
 341         }
 342
 343         t1 = fr.time;
 344
 345         counter++;
 346     }
 347     while (read_next_frame(oenv, status, &fr));
 348
 349     close_trx(status);
 350
 351     if (counter >= 4)
 352     {
 353         /* Compute time step between frames */
 354         dt = (t1-t0)/(counter-1);
 355         do_autocorr(opt2fn("-o", NFILE, fnm), oenv,
 356                     bMass ?
 357                     "Momentum Autocorrelation Function" :
 358                     "Velocity Autocorrelation Function",
 359                     counter, gnx, c1, dt, eacVector, TRUE);
 360
 361         do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
 362
 363         if (opt2bSet("-os", NFILE, fnm))
 364         {
 365             calc_spectrum(counter/2, (c1[0]), (t1-t0)/2, opt2fn("-os", NFILE, fnm),
 366                           oenv, bRecip);
 367             do_view(oenv, opt2fn("-os", NFILE, fnm), "-nxy");
 368         }
 369     }
 370     else
 371     {
 372         fprintf(stderr, "Not enough frames in trajectory - no output generated.\n");
 373     }
 374
 375     return 0;
 376 }