3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
43 #include "gmx_fatal.h"
62 static void index_atom2mol(int *n
,atom_id
*index
,t_block
*mols
)
71 while (index
[i
] > mols
->index
[mol
]) {
74 gmx_fatal(FARGS
,"Atom index out of range: %d",index
[i
]+1);
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");
84 fprintf(stderr
,"\nSplit group of %d atoms into %d molecules\n",nat
,nmol
);
89 static void precalc(t_topology top
,real normm
[]){
94 for(i
=0;i
<top
.mols
.nr
;i
++){
96 l
=top
.mols
.index
[i
+1];
100 mtot
+=top
.atoms
.atom
[j
].m
;
103 normm
[j
]=top
.atoms
.atom
[j
].m
/mtot
;
109 static void calc_spectrum(int n
,real c
[],real dt
,const char *fn
,
110 output_env_t oenv
,gmx_bool bRecip
)
116 real nu
,omega
,recip_fac
;
122 if ((status
= gmx_fft_init_1d_real(&fft
,n
,GMX_FFT_FLAG_NONE
)) != 0)
124 gmx_fatal(FARGS
,"Invalid fft return status %d",status
);
126 if ((status
= gmx_fft_1d_real(fft
, GMX_FFT_REAL_TO_COMPLEX
,data
,data
)) != 0)
128 gmx_fatal(FARGS
,"Invalid fft return status %d",status
);
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)",
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
147 recip_fac
= bRecip
? (1e7
/SPEED_OF_LIGHT
) : 1.0;
148 for(i
=0; (i
<n
); i
+=2)
151 omega
= nu
*recip_fac
;
152 /* Computing the square magnitude of a complex number, since this is a power
155 fprintf(fp
,"%10g %10g\n",omega
,sqr(data
[i
])+sqr(data
[i
+1]));
158 gmx_fft_destroy(fft
);
162 int gmx_velacc(int argc
,char *argv
[])
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."
177 static gmx_bool bMass
=FALSE
,bMol
=FALSE
,bRecip
=TRUE
;
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" }
191 gmx_bool bTPS
=FALSE
,bTop
=FALSE
;
196 /* t0, t1 are the beginning and end time respectively.
197 * dt is the time step, mass is temp variable for atomic mass.
201 int counter
,n_alloc
,i
,j
,counter_dim
,k
,l
;
203 /* Array for the correlation function */
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
}
217 #define NFILE asize(fnm)
221 CopyRight(stderr
,argv
[0]);
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
);
228 bTPS
= ftp2bSet(efTPS
,NFILE
,fnm
) || !ftp2bSet(efNDX
,NFILE
,fnm
);
232 bTop
=read_tps_conf(ftp2fn(efTPS
,NFILE
,fnm
),title
,&top
,&ePBC
,NULL
,NULL
,box
,
234 get_index(&top
.atoms
,ftp2fn_null(efNDX
,NFILE
,fnm
),1,&gnx
,&index
,&grpname
);
236 rd_index(ftp2fn(efNDX
,NFILE
,fnm
),1,&gnx
,&index
,&grpname
);
240 gmx_fatal(FARGS
,"Need a topology to determine the molecules");
241 snew(normm
,top
.atoms
.nr
);
243 index_atom2mol(&gnx
,index
,&top
.mols
);
246 /* Correlation stuff */
248 for(i
=0; (i
<gnx
); i
++)
251 read_first_frame(oenv
,&status
,ftp2fn(efTRN
,NFILE
,fnm
),&fr
,TRX_NEED_V
);
257 if (counter
>= n_alloc
) {
260 srenew(c1
[i
],DIM
*n_alloc
);
262 counter_dim
=DIM
*counter
;
264 for(i
=0; i
<gnx
; i
++) {
266 k
=top
.mols
.index
[index
[i
]];
267 l
=top
.mols
.index
[index
[i
]+1];
270 mass
= top
.atoms
.atom
[j
].m
;
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
];
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
];
282 for(i
=0; i
<gnx
; i
++) {
284 mass
= top
.atoms
.atom
[index
[i
]].m
;
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
];
295 } while (read_next_frame(oenv
,status
,&fr
));
301 /* Compute time step between frames */
302 dt
= (t1
-t0
)/(counter
-1);
303 do_autocorr(opt2fn("-o",NFILE
,fnm
), oenv
,
305 "Momentum Autocorrelation Function" :
306 "Velocity Autocorrelation Function",
307 counter
,gnx
,c1
,dt
,eacVector
,TRUE
);
309 do_view(oenv
,opt2fn("-o",NFILE
,fnm
),"-nxy");
311 if (opt2bSet("-os",NFILE
,fnm
)) {
312 calc_spectrum(counter
/2,(real
*) (c1
[0]),(t1
-t0
)/2,opt2fn("-os",NFILE
,fnm
),
314 do_view(oenv
,opt2fn("-os",NFILE
,fnm
),"-nxy");
318 fprintf(stderr
,"Not enough frames in trajectory - no output generated.\n");