src/tools/gmx_velacc.c

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