CPMD中Fourier工具包的使用：计算Power spectrum或IR spectrum
军
IR spectrum是dipole自相关函数的fourier变换， 两者均可使用CPMD提供的Fourier工具包。
CPMD提供的fourier.tar.gz, 里面提供两个工具：
1. trajc2atom.x, 把你的轨迹文件分解成n个文件，每个文件包含单个原子的速度的轨迹；
2. fourier.x，读入某个原子的速度文件，或者整个体系的total dipole的轨迹文件DIPOLE，执行fourier变换，出来的
结果第一列是频率，第二列是density，plot一下就行了。
注意，这个得到的是fourier transform of autocorrelation function, 不是
autocorrelation function。
fourier.tar.gz包里面提供的是firtran源文件，使用前要先编译. 假设你的电脑使用gfortran编译器。编译方法：
1.解压 > tar zxfv fourier.tar.gz
2. 进入其目录，修改Makefile， 把gfortran那一个block前面的#去掉，如下：
######################################################
# gfortran compiler v4.1.x and up
#OPENMP=-fopenmp
######################################################
F90=gfortran $(OPENMP)
F77=gfortran
LD=$(F90)
F90FLAGS=-O2 -funroll-loops -g -Wall -ffree-form -std=f95 -pedantic -
pg
F77FLAGS=-O2 -funroll-loops -g -Wall -ffixed-form -std=gnu
LDFLAGS=$(F90FLAGS)
######################################################
保存退出后输入> make开始编译
完成后make clean你就会发现多出了两个可执行文件： fourier.x trajec2atom.x
把上述两个文件复制到任何地方，比方你的工作目录，运行
> ./trajec2atom.x 它会提示你输入trajectory文件名input，依提示进行即可；
fourier.x也是这样，根据提示输入。
如果处理很多文件，你可以依次把输入写入一个文本文件比如: job.in， 然后
> ./fourier.x < job.in 即可自动执行批量作业；
fourier.x要求输入 inputfile name, output file name, 采样时间间隔
（timestep*sampling frequency）, 温度（e.g. 300）， 最大频率(e.g. 3000)，分辨率（e.g. 1
或100，数值越大，得到的图越coase，但是有更少的noise）. 下面是我的一个输入文件d.in
DIPOLE_d //输入文件名，我算的是dipole自相关，可以得到IR spectrum,
IR_d_1.dat //输出文件名
40. // 我的步长是4a.u., sampling frequency 10 step
300. // 温度
3000. // 最大频率
1 // 分辨率
其它提示
1. 一般的数据处理并不要很长时间，这个可以在本地电脑上执行；
2. 本地电脑可以方便的安装编译器，比如在ubuntu下，
> sudo apt-get install gfortran
即可安装好gfortran.
3. fourier.x 的执行如下：
它读入输入文件的第4-6列，依次是dipole或速度的x, y, z分量；获取总的记录数和平均值；
然后根据你的采样频率计算频率步长（wave_factor，以波数表示）
然后对每个频率（依波数表示cm-1）在整个轨迹长度上积分<A(t)*A(0)>*exp(-i*w*t) , 得到IR intensity或者
power density
<A(t)*A(0)>就是自相关函数
参考文献：
1. CPMD3.13.2 manual;
2. fourier工具包中的README文件和源代码；
3. Wiki: autocorrelation function
4. CPMD计算IR spectrum, M.P. Gaigeot, JCTC, 1, 772-789, 2005