24小时热门版块排行榜    

查看: 2442  |  回复: 5
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

Waxwane

新虫 (初入文坛)

[求助] 关于Ewald sum算法的相关问题 已有1人参与

自己作死地用C编写Monte Carlo的代码,模拟的具体问题是粘土矿物对Ar的吸附。
然后需要加上Ewald Sum算法的问题算矿物羟基的长程力。
但是这个矿物是三斜晶系的,书上给的Ewald Sum是立方格子。
我查过LAMMPS里面,前几年有个Ewald/n方法专门算三斜格子,现在和Ewald本身合并了。
我决定用三斜的格子算,因为矿物的周期性是三斜的,用立方重复不了。但是实验室其他人不搞算法,他们说关于三斜粘土矿物的文献,没有听过把格子划分为三斜格子的说法。
我认为算短程力的格子可以是立方的,但长程力不行,那些文献里没有写清楚。
不知道我的观点正不正确,因为我只能找到立方格子的源码和伪代码,没见过三斜的代码。
还有,如果有编程简单一点的,被大家比较公认的算长程力的算法,求推荐。我已经要被Ewald sum坑死了。
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

王慧1021

金虫 (小有名气)

粘土矿物没有立方的构型么
5楼2015-01-13 11:48:52
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 6 个回答

lsloneil

专家顾问 (正式写手)

【答案】应助回帖

★
感谢参与,应助指数 +1
Waxwane(月只蓝代发): 金币+1, 鼓励交流 2015-01-11 12:37:23
Ewald的k空间求和里的傅立叶变换还是在直角坐标系下进行,理论上来讲跟你用什么晶系关系不大。也就是说你依然可以把盒子里的电荷密度做一个傅立叶变换然后求和。

不过有一点技术上我没太想明白的是在三斜晶系里做傅立叶变换,倒格子的格矢如何选取。立方盒子格矢很容易选取,直接2pi/L就可以了。三斜晶系下可能要做点额外处理来获得一组正交的倒格子格矢?

实空间部分求和跟你的晶系有一定关系,但这一部分很容易处理,因为只是个短程相互作用+考虑minimum image convention。
2楼2015-01-11 03:35:12
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

Waxwane

新虫 (初入文坛)

引用回帖:
2楼: Originally posted by lsloneil at 2015-01-11 03:35:12
Ewald的k空间求和里的傅立叶变换还是在直角坐标系下进行,理论上来讲跟你用什么晶系关系不大。也就是说你依然可以把盒子里的电荷密度做一个傅立叶变换然后求和。

不过有一点技术上我没太想明白的是在三斜晶系里做 ...

实空间部分确实很好做,我个人认为那个虚构的格子是三斜格子,也也就是说那三个矢量是按晶胞构建的,因为他叠加的部分只算的距离,不管是正交的3矢量,还是不正交的三矢量,最后算的距离都是一样的。而用和晶胞一样的矢量算,周期边界方面好操作,正交的周期边界构建不了。
3楼2015-01-11 18:31:08
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lsloneil

专家顾问 (正式写手)

【答案】应助回帖

引用回帖:
3楼: Originally posted by Waxwane at 2015-01-10 22:31:08
实空间部分确实很好做,我个人认为那个虚构的格子是三斜格子,也也就是说那三个矢量是按晶胞构建的,因为他叠加的部分只算的距离,不管是正交的3矢量,还是不正交的三矢量,最后算的距离都是一样的。而用和晶胞一样 ...

看不懂你在说什么。

我又看了一下DL_POLY中Ewald sum的代码,对于所有晶格,它们k空间部分的处理是这样的。

1.先从晶格的格矢得到倒格子(reciprocal lattice)的格矢
2. 由倒格子的格矢再构建一组正交的k空间矢量
3. 有了这组正交的k空间矢量,就可以获得kmax并进行下一步k空间的计算了。

对于步骤2,如果是立方盒子,这一步相当于什么都没做,因为原先倒格子的格矢就是正交的。但对于三斜的格子,因为其倒格子格矢不是正交的,这一步就可以获得一组正交的k空间矢量。步骤2的算法的大概意思是,对于三斜晶系,总可以找到一组正交的格矢,它们形成的晶胞体积跟原晶胞相等。

如果你想看dlpoly的代码,可以在网上搜DL_POLY classic的代码(下载需要注册,但免费)。Ewald的代码在文件ewald_module.f里,参看subroutine ewald1.
4楼2015-01-12 01:20:43
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考博] 售SCI一区T0P文章,我:8.O55.1.O.54,科目全,可十急 +5 qR2AdCWTVtkw 2026-10-09 9/450 2026-10-11 22:15 by XQNyDCa1Tk5n
[考博] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +3 qR2AdCWTVtkw 2026-10-09 7/350 2026-10-11 22:11 by XQNyDCa1Tk5n
[基金申请] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +5 qR2AdCWTVtkw 2026-10-09 9/450 2026-10-11 22:06 by XQNyDCa1Tk5n
[教师之家] 售SCI文章,我:8O5.5.1.O.54,科目齐全,可+急 (EPI+-1)(金币-50) +5 qR2AdCWTVtkw 2026-10-09 10/500 2026-10-11 22:06 by XQNyDCa1Tk5n
[有机交流] OpenAI又攻克数学难题了,离AI做有机合成还有多远? +3 wangyl_123 2026-10-10 6/300 2026-10-11 15:04 by Abin0406
[教师之家] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 (EPI+-1)(金币-50) +5 qR2AdCWTVtkw 2026-10-09 10/500 2026-10-11 14:25 by VOMexwnDzOFZ
[公派出国] 售SCI一区文章,我:8O5.5.1.O5.4,科目全,可伽急 +4 qR2AdCWTVtkw 2026-10-09 9/450 2026-10-11 14:13 by VOMexwnDzOFZ
[教师之家] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 (EPI+-1)(金币-50) +4 qR2AdCWTVtkw 2026-10-09 7/350 2026-10-11 14:05 by VOMexwnDzOFZ
[考研] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +5 qR2AdCWTVtkw 2026-10-09 8/400 2026-10-11 13:51 by VOMexwnDzOFZ
[硕博家园] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +5 qR2AdCWTVtkw 2026-10-09 8/400 2026-10-11 13:45 by VOMexwnDzOFZ
[找工作] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +3 qR2AdCWTVtkw 2026-10-09 7/350 2026-10-11 10:54 by VOMexwnDzOFZ
[找工作] 售SCI一区T0P文章,我:8.O55.1.O.54,科目全,可十急 +4 qR2AdCWTVtkw 2026-10-09 7/350 2026-10-11 10:53 by VOMexwnDzOFZ
[公派出国] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +4 qR2AdCWTVtkw 2026-10-09 8/400 2026-10-11 10:53 by VOMexwnDzOFZ
[找工作] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +4 qR2AdCWTVtkw 2026-10-09 7/350 2026-10-11 10:34 by VOMexwnDzOFZ
[考博] 招收博士生一名 +3 mandolin 2026-10-10 3/150 2026-10-11 10:21 by 不是山谷!
[教师之家] 课题组会帮你通过考核么? +3 好运来~9527 2026-10-10 5/250 2026-10-11 09:00 by 好运来~9527
[公派出国] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +4 qR2AdCWTVtkw 2026-10-09 4/200 2026-10-11 03:14 by M0zrYVSAnxpE
[有机交流] 各位大佬,路线设计求助 +6 zapen 2026-10-08 8/400 2026-10-10 15:59 by zapen
[基金申请] 入职后第一次管课题组的专利缴费,对账才发现官费比别人多掏了一倍 +5 13108017953 2026-10-08 5/250 2026-10-10 11:15 by wdwd123
[有机交流] 我应该选哪条合成路线让实验成功率高一些? +5 jing816y 2026-10-08 9/450 2026-10-10 09:05 by zapen
信息提示
请填处理意见