src/tools/gmx_velacc.c

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