请选择 进入手机版 | 继续访问电脑版

分子模拟论坛 Molecular Simulation Forums

 找回密码
 立即注册

QQ登录

只需一步,快速开始

扫一扫,访问微社区

搜索
查看: 102|回复: 3

速度相关函数

[复制链接]

21

主题

337

帖子

701

积分

高级会员

Rank: 4

积分
701
发表于 2012-6-1 10:32:00 | 显示全部楼层 |阅读模式
本帖最后由 cuiloveyang 于 2012-6-1 23:38 编辑
这几天做了一下速度相关函数,很简单的一个东西,居然搞了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算扩散系数不靠谱。

回复

使用道具 举报

6

主题

19

帖子

48

积分

新手上路

Rank: 1

积分
48
发表于 2012-6-1 18:54:00 | 显示全部楼层
替乐乐兄顶一顶!
回复

使用道具 举报

8

主题

22

帖子

54

积分

注册会员

Rank: 2

积分
54
发表于 2013-9-28 19:47:00 | 显示全部楼层
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一个效果
回复

使用道具 举报

0

主题

1

帖子

14

积分

新手上路

Rank: 1

积分
14
发表于 2017-12-3 08:54:51 | 显示全部楼层
看不到附件?研究了好多天acf也搞不明白。。。
回复

使用道具 举报

Archiver|手机版|小黑屋|分子模拟论坛  

GMT+8, 2017-12-13 05:34 , Processed in 0.305098 second(s), 28 queries .

Powered by Discuz! X3.3

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表