这几天做了一下速度相关函数，很简单的一个东西，居然搞了3天。真是很郁闷。
不过要特别感谢剑心同学的大力和无私的帮助。同时也感谢mitbbs网上的一些讨论。
他用的是python，实在不懂，自己就编了一个c程序，通过fftw的快速傅立叶变换，也可以得到相同的结果。
既然得到了大家的帮助，自己做成后，也应该分享一下。
附近是速度相关函数的程序。
我把dump出来的格式限制了，dump 1 all custom 10 chain id vx vy vz
“chain” 是文件名，大家可以修改源程序里面的文件名
就是原子的序号，然后就是3个方向的速度，读起来比较快，文件也能小一些。
如果用直接计算方法，还可以并行，但是加入了fftw后，并行没通过，所以放弃了。
但是计算后发现，似乎主要的时间都在读文件上面，计算的时间相对一般，所以单核运算也可以接受。
我一般喜欢现在windows下面写好程序，编译什么的，但是需要到fftw的官网上面下载win版本的fftw，很容易，解压后做成lib，就可以直接用了。
在Linux下，首先要编译fftw，这个很容易，我最后把路径放在了/opt/fftw里。
Linux下编译程序的时候需要给出fftw的路径
icc Corr.c -o Vel -O3 -xS -axS -I/opt/fftw/include -L/opt/fftw/lib -lfftw3
Vel就是生成的可执行文件。
最后输出的就是一个文件，第2，3，4列是Vx，Vy，Vz的自相关函数，5列是速度自相关函数
接下来的4列是一样的，只不过是用fftw算得。
其实我一点都不懂快速傅立叶变换什么的，当年学习的时候就特别不爱学它。
不过幸好陈正隆的书上介绍的还算详细，就按照他说的思路直接编了，结果居然和直接计算的结果一样，那就OK了，也没有做过多的检查。
其实检查了我也不懂。^_^
忘了说一件事，当时刚写完直接计算方法后，一直想验证写的是否正确，因为实在是太简单了。
我就想用扩散系数来判断，速度相关函数的积分是扩散系数的3倍，即3D，而MSD的斜率是扩散系数的6倍。
如果积分的结果正确的话，我就可以肯定的说，程序没错。
首先，我自己算了MSD的结果，和文献一样，所以没有问题。
但是这个速度相关函数的积分却不理想，我尝试过增加计算step（5w-100w），增加相关函数的点（1000-10000）。
但是结果都很不一样，有的时候差很远，而且sum（y）和trapz（y）的结果有时候也完全不一样。
如果用什么样条插值，结果又不一样。最后只能放弃。
下面说的只是我的猜测，我并没有做过深入的研究。
我做的这个系统的D是0.0011，那么速度相关函数的积分应该就是0.0033。
如果看图的话，可以很容易发现，后面大部分都是0，但是放大的话，就发现在0附近震荡，并且大部分值都是10-3,-4,-5的。
这些值虽然很小，但是对于正确的积分结果0.0033来说，可能就比较大了，那么这样的累加

可能

会对数值积分产生影响。

上面的这个是我的想法，测试了一天也没什么特别好的结果，如果大家以后碰到这个问题了，欢迎讨论。
或者有人知道这方面的知识，欢迎提出来，互相学习一下。
程序验证2，就是用FFT的结果来判断，最后发现结果基本上一样，假装说明程序编对了。
最新更新，算了个其他材料，用MSD算得扩散系数D是1.21×10-7，用VACF算得D是1.71×10-7。
再次说明了用VACF算扩散系数不靠谱。

1# cuiloveyang
这个程序有问题吧，下面是里面的一段代码
for i1=1:1:e
a1=a{i1};
for i2=i1:1:e
   m=0;
   n=0;
   a2=a{i2};
   for i3=1:1:number
     for i4=1:1:number
     if a1(i3,1)==a1(i4,1)
       m=ma1(i3,3)*a2(i4,3)a1(i3,4)*a2(i4,4)a1(i3,5)*a2(i4,5);
       i4=number;
     end
     end
   end
end
J(i1,i2)=m/number;
end
#===========
m=0写在了，for i2=i1:1:e的上面，那就是说i2的前e-1个循环，m都为0，这样的话，和
for i2=e:1:e一个效果