应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 分子模拟 » Gromacs/Amber/NAMD » 【求助】Ewald长程力中Fourier空间形状问题?
101/1返回列表
查看: 727  |  回复: 9
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前主题已经存档。

liuyalihd

铜虫 (小有名气)

[交流] 【求助】Ewald长程力中Fourier空间形状问题?

我看某文献1中计算Ewald长程力时,分成实空间和Fourier空间.
实空间中其实还是在截断半径内计算两个粒子之间的作用. Fourier空间是在半球形空间(z轴0到nz, y轴从-ny到ny, x轴从-nx到nx)内计算.

但在gromacs4.0.5中,我看了Ewald.c的代码,发现它没有分成实空间和Fourier空间两部分计算.它似乎仅在Fourier空间中计算.
最让我奇怪的是它的Fourier空间不是半球形,而是在半正方体(z轴0到nz, y轴从-ny到ny, x轴从-nx到nx)中计算.

请问为何gromacs4.0.5中计算Ewald的Fourier空间形状不是传统的球形呢?
gromacs4.0.5中计算Ewald是咋回事?为何没有实空间的计算呢?

[ Last edited by zdhlover on 2009-11-15 at 02:29 ]
回复此楼

» 猜你喜欢
1楼2009-09-01 20:34:38
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

bay__gulf

金虫 (著名写手)

刘苏州
★ ★ ★
小木虫(金币+0.5):给个红包，谢谢回帖交流
yjcmwgk(金币+2,VIP+0):感谢交流! 9-2 09:20
gmx 为保证速度使用了大量的数学和编程技巧，如此一来就失去了可读性
你再仔细看看，可能有你没有注意到的地方

仅作参考，我功力没你深，还没有到阅读源代码的级别
赞一下(7人) 回复此楼
2楼2009-09-01 20:52:07
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

liuyalihd

铜虫 (小有名气)
谢谢bay__gulf

很惭愧, 我是写code的, 不太懂MD.
我想问从理论上说, Fourier空间应该是球形的吧! 我没理解错吧?
赞一下 回复此楼
3楼2009-09-01 22:22:09
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

yongleli

木虫 (正式写手)
★ ★ ★ ★
小木虫(金币+0.5):给个红包，谢谢回帖交流
lei0736(金币+3,VIP+0):谢谢 9-2 16:11
没听说过fouirer空间是球形的。。。
Ewald求和涉及一个正格空间到倒格空间的变换。
这两个空间只是单位向量不同，
一个是\vec{a}_i，i=1,2,3
一个是\vec{b}_i=2\pi\vec{a}_{j}\vec{a}_{k}/\vec{a}_i\dot\vec{a}_j\cross\vec{a}_k; i=1,2,3

所以倒空间对称性随正空间特征而定。
赞一下(11人) 回复此楼
4楼2009-09-02 15:01:01
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

liuyalihd

铜虫 (小有名气)
请问: yongleli  
因为我看code中有三个嵌套for循环 (分别是  nz : z轴0到5,  ny: y轴从-5到5,  nx: x轴从-5到5), 内有个条件判断,镜像空间与真实空间之间距离要小于25.

这也就近似相当于从真实模拟空间开始向外(这个外是指无穷无尽的镜像)扩展,x轴(y轴)正负各镜像5份,z轴仅仅往正方向镜像5份(这是为了计算量减半),对这些镜像出来的空间计算Fourier.  
关键是代码中有个if语句,判断这个镜像空间与真实空间(真实空间为0,0,0)的距离小于25才真正做计算.这不就是相当于在真实空间外镜像出近似半球形的空间做计算吗?

我把代码写在下面: 请多指教!


  计算三角函数并做成表;
  for (nz = 0; nz <= 5; nz ++) {
    for (ny = - 5; ny <= 5; ny ++) {
      for (nx = - 5; nx <= 5; nx ++) {

        len=nx*nx+ny*ny+nz*nz;
        if (len > 5*5 ) continue;
        
        查表,并逐个粒子计算Fourier;
     }
    }
   }
赞一下 回复此楼
5楼2009-09-02 22:02:00
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zixi

金虫 (小有名气)
★ ★ ★
小木虫(金币+0.5):给个红包，谢谢回帖交流
yjcmwgk(金币+2,VIP+0): 9-4 13:43
没读过gromacs的代码，不过感觉是不是这样的，你可以读一下代码求证。

1 : 关于只有Ewald空间的部分
     一般将Ewald求和分成Real空间和Ewald空间，感觉这样的好处是，如果需要对长程力做截断处理而不做Ewald求和，那么只要取Ewald的Real空间部分，并取高斯分布电荷的参数alpha为0，就可以了，不需要另外写代码。如果需要Ewald求和，那么直接在Real空间上添加上Ewald空间的能量值就可以了。这样增加的代码的重用性。Real空间的求和可能在其他文件里面，你仔细找找。

2. Ewald空间使用半正方形
    我看你描述的k值的半正方形和半球形的取值范围没有多大区别嘛？！都是z轴0到nz, y轴从-ny到ny, x轴从-nx到nx。
既然Ewald空间的求和随着k值的增加迅速收敛的，那么我只要取足够大的k值，使其包含所必需的球形空间，所得的结果差别多少呢？ 而且实际操作起来，正方形空间比球形空见更容易处理，为什么不用正方形空间呢？
赞一下(7人) 回复此楼
6楼2009-09-04 02:30:07
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

liuyalihd

铜虫 (小有名气)
致  zixi:

关于您的第一条：非常感谢，我再找找real空间的计算。

关于您的第二条：
1）原理部分： 前面是我没有讲清楚。球形空间是Ewald求和法的周期作用区域，这在很多书中都提到了（如：分子模拟-陈正隆，9.2节 图9-5）。
2）代码部分：另外，我给出的代码不是gromacs的。但是我看在gromacs的ewald.c（src/mdlib/ewald.c）中有函数do_ewald()，我看代码非常像是计算Ewald中的Fourier计算。但是它没有if判断镜像是否在球形空间中，我认为应该是球形才计算，在正方体空间计算不符合原理呀！？
3）正方体空间与球形空间相比，去掉if语句的开销不提，球形空间减少了一半的计算量。
赞一下 回复此楼
7楼2009-09-04 12:37:02
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

yongleli

木虫 (正式写手)
★ ★ ★ ★
小木虫(金币+0.5):给个红包，谢谢回帖交流
yjcmwgk(金币+3,VIP+0):感谢您的真知灼见。管理层决定：授予您2010年炸药奖 9-5 09:56
我又看了一下理论， 在Computer Simulations of Liquids里，pp. 157, Ewald求和的范围是在球形区域里的。并不是说福利叶变换是把实空间变成了球形空间。

我编的Ewald求和用的是pp.159页的公式，Gromacs没看过，应该跟书里那个公式差不多吧。建议最好看看原始文献。有一本中文书里详细解释了这个理论和算法。

黄昆，Max Born，晶格动力学，北大版。

很老的书了，但是近年出了新版。
赞一下(7人) 回复此楼
8楼2009-09-04 22:43:26
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

liuyalihd

铜虫 (小有名气)
致: yongleli
非常感谢!

1)Ewald求和的范围是在球形区域里的.
a)它的实空间计算是在模拟系统自身(周期性边界条件)计算的,  截断半径是L/2.   L是模拟系统的边长(假设立方体空间). 这里与球形空间没关系!
b)它的Fourier空间计算是对镜像出的多个镜像系统计算,多个镜像组成近似球形空间的. 对某个镜像空间计算的前提是它必须在极限球内,那么才对其做Fourier计算啦!
上述是我看The Art of Molecular Dynamics Simulation得出的.有错误请指正!

2)我看了Computer Simulations of Liquids中的Ewald求和.它给了Fortran程序F22计算Ewald的.可原代码要给UK的人写信,或者从USA上ftp下载(见P343).那个ftp我上不去,请有源代码的发给我,谢谢了!
在页面http://muchong.com/html/200803/729220.html中
lei0736 给出的地址http://www.fs2you.com/files/69a170a6-eddf-11dc-93a7-0014221b798a/我登不上!

我信箱liuyalihd@163.com
赞一下 回复此楼
9楼2009-09-05 00:02:34
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

yongleli

木虫 (正式写手)
★ ★ ★ ★ ★ ★ ★
小木虫(金币+0.5):给个红包，谢谢回帖交流
yjcmwgk(金币+6,VIP+0):感谢原创交流。管理层决定：授予您加入2010年中国科学院院士增选名单的资格 9-5 13:56
这是我用的那个f22。

不过它比较简陋。
我写的时候修改了一些。
建议当个模版看看好了。
赞一下(7人) 回复此楼
10楼2009-09-05 12:29:58
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 liuyalihd 的主题更新
101/1返回列表
普通表情 龙 兔 虎 猫 高级回复(可上传附件)
信息提示
关闭
请填处理意见
关闭