+ snew(data,n*2);
+ for(i=0; (i<n); i++)
+ data[i] = c[i];
+
+ if ((status = gmx_fft_init_1d_real(&fft,n,GMX_FFT_FLAG_NONE)) != 0)
+ {
+ gmx_fatal(FARGS,"Invalid fft return status %d",status);
+ }
+ if ((status = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX,data,data)) != 0)
+ {
+ gmx_fatal(FARGS,"Invalid fft return status %d",status);
+ }
+ fp = xvgropen(fn,"Vibrational Power Spectrum",
+ bRecip ? "\\f{12}w\\f{4} (cm\\S-1\\N)" :
+ "\\f{12}n\\f{4} (ps\\S-1\\N)",
+ "a.u.",oenv);
+ /* This is difficult.
+ * The length of the ACF is dt (as passed to this routine).
+ * We pass the vacf with N time steps from 0 to dt.
+ * That means that after FFT we have lowest frequency = 1/dt
+ * then 1/(2 dt) etc. (this is the X-axis of the data after FFT).
+ * To convert to 1/cm we need to have to realize that
+ * E = hbar w = h nu = h c/lambda. We want to have reciprokal cm
+ * on the x-axis, that is 1/lambda, so we then have
+ * 1/lambda = nu/c. Since nu has units of 1/ps and c has gromacs units
+ * of nm/ps, we need to multiply by 1e7.
+ * The timestep between saving the trajectory is
+ * 1e7 is to convert nanometer to cm
+ */
+ recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;
+ for(i=0; (i<n); i+=2)
+ {
+ nu = i/(2*dt);
+ omega = nu*recip_fac;
+ /* Computing the square magnitude of a complex number, since this is a power
+ * spectrum.
+ */
+ fprintf(fp,"%10g %10g\n",omega,sqr(data[i])+sqr(data[i+1]));
+ }
+ xvgrclose(fp);
+ gmx_fft_destroy(fft);
+ sfree(data);
+}