2 * This file is part of the GROMACS molecular simulation package.
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,2013, 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.
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.
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.
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.
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.
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.
46 #include "gmx_fatal.h"
65 static void index_atom2mol(int *n,atom_id *index,t_block *mols)
74 while (index[i] > mols->index[mol]) {
77 gmx_fatal(FARGS,"Atom index out of range: %d",index[i]+1);
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");
87 fprintf(stderr,"\nSplit group of %d atoms into %d molecules\n",nat,nmol);
92 static void precalc(t_topology top,real normm[]){
97 for(i=0;i<top.mols.nr;i++){
99 l=top.mols.index[i+1];
103 mtot+=top.atoms.atom[j].m;
106 normm[j]=top.atoms.atom[j].m/mtot;
112 static void calc_spectrum(int n,real c[],real dt,const char *fn,
113 output_env_t oenv,gmx_bool bRecip)
119 real nu,omega,recip_fac;
125 if ((status = gmx_fft_init_1d_real(&fft,n,GMX_FFT_FLAG_NONE)) != 0)
127 gmx_fatal(FARGS,"Invalid fft return status %d",status);
129 if ((status = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX,data,data)) != 0)
131 gmx_fatal(FARGS,"Invalid fft return status %d",status);
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)",
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
150 recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;
151 for(i=0; (i<n); i+=2)
154 omega = nu*recip_fac;
155 /* Computing the square magnitude of a complex number, since this is a power
158 fprintf(fp,"%10g %10g\n",omega,sqr(data[i])+sqr(data[i+1]));
161 gmx_fft_destroy(fft);
165 int gmx_velacc(int argc,char *argv[])
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."
180 static gmx_bool bMass=FALSE,bMol=FALSE,bRecip=TRUE;
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" }
194 gmx_bool bTPS=FALSE,bTop=FALSE;
199 /* t0, t1 are the beginning and end time respectively.
200 * dt is the time step, mass is temp variable for atomic mass.
204 int counter,n_alloc,i,j,counter_dim,k,l;
206 /* Array for the correlation function */
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 }
220 #define NFILE asize(fnm)
224 CopyRight(stderr,argv[0]);
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);
231 bTPS = ftp2bSet(efTPS,NFILE,fnm) || !ftp2bSet(efNDX,NFILE,fnm);
235 bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,NULL,NULL,box,
237 get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
239 rd_index(ftp2fn(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
243 gmx_fatal(FARGS,"Need a topology to determine the molecules");
244 snew(normm,top.atoms.nr);
246 index_atom2mol(&gnx,index,&top.mols);
249 /* Correlation stuff */
251 for(i=0; (i<gnx); i++)
254 read_first_frame(oenv,&status,ftp2fn(efTRN,NFILE,fnm),&fr,TRX_NEED_V);
260 if (counter >= n_alloc) {
263 srenew(c1[i],DIM*n_alloc);
265 counter_dim=DIM*counter;
267 for(i=0; i<gnx; i++) {
269 k=top.mols.index[index[i]];
270 l=top.mols.index[index[i]+1];
273 mass = top.atoms.atom[j].m;
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];
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];
285 for(i=0; i<gnx; i++) {
287 mass = top.atoms.atom[index[i]].m;
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];
298 } while (read_next_frame(oenv,status,&fr));
304 /* Compute time step between frames */
305 dt = (t1-t0)/(counter-1);
306 do_autocorr(opt2fn("-o",NFILE,fnm), oenv,
308 "Momentum Autocorrelation Function" :
309 "Velocity Autocorrelation Function",
310 counter,gnx,c1,dt,eacVector,TRUE);
312 do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy");
314 if (opt2bSet("-os",NFILE,fnm)) {
315 calc_spectrum(counter/2,(real *) (c1[0]),(t1-t0)/2,opt2fn("-os",NFILE,fnm),
317 do_view(oenv,opt2fn("-os",NFILE,fnm),"-nxy");
321 fprintf(stderr,"Not enough frames in trajectory - no output generated.\n");