24小时热门版块排行榜    

CyRhmU.jpeg
查看: 6532  |  回复: 101
本帖产生 3 个 模拟EPI ,点击这里进行查看

ghcacj

荣誉版主 (著名写手)

阿超

优秀版主

★ ★ ★
小木虫(金币+0.5):给个红包,谢谢回帖交流
zh1987hs(金币+2):鼓励交流 2010-11-16 16:31:30
设置全一样了,可是结果还是16.5左右,怎么回事,但是我将Ewald参数改成KMAX=7,KAPPA=3时,结果反而和文献比较吻合了,请问这个参数如何取比较合适?
61楼2010-11-11 14:32:15
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ghcacj

荣誉版主 (著名写手)

阿超

优秀版主


小木虫(金币+0.5):给个红包,谢谢回帖交流
换成你提供的CuBTC.mol计算结果居然就成了,找到原因了。呵呵

[ Last edited by ghcacj on 2010-11-11 at 16:51 ]
62楼2010-11-11 15:36:04
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)

引用回帖:
Originally posted by ghcacj at 2010-11-11 15:36:04:
换成你提供的CuBTC.mol计算结果居然就成了,找到原因了。呵呵

[ Last edited by ghcacj on 2010-11-11 at 16:51 ]

是什么原因?

还有, 如果有兴趣的话, 一起测试一下NH3在MOF里面的吸附, 从 CuBTC开始.
Life, Love, Laugh.
63楼2010-11-12 01:06:29
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
ghcacj(金币+10):谢谢 2010-11-12 13:04:00
Update  Total Energy vs. No. of Iteration

因过竹园逢僧话, 偷得浮生半日闲.

上图...





有几个说明:

(1) 该图是测试参数时候所得, 在该例子中CO2的VDW参数已经修改得面目全非, 所以图中的能量绝对值, 没啥意义, 不要试图重现.

(2) GCMC的控制文件中, 第一部分如下

------ General Information ------------------------------------------  
CO2 in CuBTC
20000000   
100   
1000000  
100  
1
91283
3                    
CO2.CuBTC.res
CO2.CuBTC.con

大家可以看到总步数为20000000, 输出频率为 100.

(3) POST的Control文件中, 两个地方需要注意:

50, 0        # Percentages of data to skipped at start and end

这是说我将总步数的前50%扔掉.

2000     # Number of blocks into which data should be divided for stats

这是说, 我将要统计的所有输出, 平均分成2000个Block.

那么就是说:
统计包括了 2E7*50%=1E7 步GCMC;
输出文件(.con文件) 记录了1E7/100=1E5 个 image;
能量图中记录了 2000个点 (这个数值就是ctr.post文件中的那个2000);
每两个点之间的间隔是: 1E5/2000=50 个 image, 也就是 50*100=5000 步GCMC;

这是为什么, 如果你看能量图的横坐标, 第一点是10005000, 第二点是 10010000, 以此累加, 直到最后一点, 20000000.

Clean and Clear ?

^-^
Life, Love, Laugh.
64楼2010-11-12 04:43:25
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ghcacj

荣誉版主 (著名写手)

阿超

优秀版主


小木虫(金币+0.5):给个红包,谢谢回帖交流
原因是那个结构是很久之前建立的,用了FORCITE中的UFF力场优化了下,没想到悲剧了
65楼2010-11-12 08:26:19
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
ghcacj(金币+15, 模拟EPI+1):very very good 2010-11-14 08:43:48
引用回帖:
Originally posted by ghcacj at 2010-11-11 14:32:15:
设置全一样了,可是结果还是16.5左右,怎么回事,但是我将Ewald参数改成KMAX=7,KAPPA=3时,结果反而和文献比较吻合了,请问这个参数如何取比较合适?

针对这个问题, 做一些update.

先来看一下emap的mol_mol_file:


Lattice Lattice NCOUL OFF
Lattice Lattice COUL OFF
Probe Lattice NCOUL OFF
Probe Lattice COUL SUM FAST FIXED EWALD SFACTOR KMAX@15 KAPPA@6.7 LOCUT@1e-10
Probe Probe NCOUL OFF
Probe Probe COUL OFF

首先需要指出的是, 这里除了 Probe 与 Lattice的COUL作用之外, 其他的相互作用全部关闭.
来仔细看看这句:
Probe Lattice COUL SUM FAST FIXED EWALD SFACTOR KMAX@15 KAPPA@6.7 LOCUT@1e10
(1) “Probe” 和 “Lattice”, 这个就像“路人甲”一样, 你想怎么称呼都可以;
(2) “COUL”是从 ff.F90过来的, 在那个源程序里面定义了体系的能量从COUL和NONCOUL计算得到;
(3) “SUM”, “FAST”, “FIXED”,“EWALD”来自于ssdrive.F90, 我们来一一解读一下.
(a) SUM
Line 94, 定于了 “SUM”:  “Type(SSSumParams), Pointer       :: sum”. 很容易理解, 在任意位置上, Probe与Lattice的COUL作用, 需要累加起来.
(b) FAST
Line 495,是注释行, 关于“FAST”是这样说的, “fast -- True (False) => evaluate "Fast" (Slow) interactions”
字面上的理解就是, 用了“FAST”那么计算的时候, 关于相互作用的计算就能“加速”,至于如何加速? 不知道, 也不需要知道.
(c) FIXED
没有从源程序找到很好的对应, 但是这个是和Control文件中对于Lattice的处理对应的, 我们将Lattice文件Fixed处理, 那么这里也是需要告诉mapmake文件这个信息. 除此之外, 没啥用处.
(d) EWALD
这个指定了如果生成EMAP, 前面有提到过, GCMC的EMAP只能通过EWALD得到, 而on-the-fly COUL计算(i.e. 不用EMAP, 在GCMC中直接计算COUL)只能WFCOUL而不能使用EWALD.
Line542再次证明了这一点认识: “on the fly ewald sums are implemented for md only”
也就是说, MD时, on the fly COUL计算, 是可以采用EWALD SUM的, 但是MC(包括GCMC)不行....

(4) SFACTOR, KMAX, KAPPA, LOCUT是EWALD方法中的参数. 具体的原理,参见Frenkel&Smit 的第十二章, 或者MuSiC的相关文档. 你也可以从ewald.F90 的subroutine (Line 152 to Line 269) ewald_init看到程序中是怎样指定这些参数的. 另外, sssum.F90也很值得一看, 如果你想了解更多EWALD是怎样在MuSiC中实现的话.
下面一一分析一下这四个参数.
(a) SFACTOR
Line 205 to Line 207:      
        Case ('sfactor')
        !** We want to use the structure factor.
        eparams%useSF = .TRUE.
它是一个逻辑变量, 默认值是.FALSE., 所以在这里需要加入这个关键词来激活它, 变为.TRUE.
激活这个逻辑变量以后, 程序做了一些什么呢? 参加 Line 284 那边的注释.
      !** Call the structure factor routine. This calculates the k space
      !** vectors we will use later when computing the Ewald summation. If
      !** the coordinates of the ions are fixed, this needs to be calculated
      !** only once. If they change, it must be recalculated.   
(b) KMAX
Line 184: “Number of k "boxes" in each direction (x,y,z)”. 你可能注意到在程序里面, KAMX is hard coded to be 5000, 你可能回想增大这个数值, 来处理更大的体系. 但是, 不幸的是, 你的内存很有可能比KMAX参数更快地达到阀值, 所以, 修改这个参数在实际使用上, 没有意义, 我没有测试过, 但是我觉得这里 KMAX@15, 其实15这个数值, 是没什么意义的, 只是需要一个input而已, 你用15, 10还是8, 没有区别.
(c) KAPPA 是和 EWALD SUM里面的Gaussian function有关的一个参数, 它定义了 Gaussian Function的带宽.
如果你看Frenkel&Smit 的第十二章, 这里程序中的KAPPA参数值等于Frenkel&Smit Chap.12中的alpha^0.5.
(d) LOCUT
这个很好理解, 请看ewald.F90文件中的这一行:
Line 200, “   !** Low end cutoff for the reciprocal space k vectors.”

接下来的问题就是, 我们应该如何选择合适的KAPPA和KMAX参数? 这两个参数的选择是基于EWALD SUM 计算时对于实空间和倒空间的“利用”和平衡. 一般而言, KMAX必须让SUMMATION在倒空间的计算满足随着KAPPA参数收敛, 同时尽量让实空间的计算也能快速收敛. 比方说, 减小KAPPA参数可以使得倒空间的收敛加快, 但是此时实空间的收敛就会减慢. 在非常小的KAPPA参数下, 实空间的EWALD SUM 和WFCOUL方法一样. 在其他MAP MAKE 参数相同的情况下, 测试不同的KAPPA数值, 看生成MAP需要的不同时间, 来找到最优值, 但是一般来讲, 这个对于MC(包括GCMC, 和MC里面相应的EMAP的生成)没有实际意义. 但是对于MD来讲, 由于on the fly COUL计算可以采用EWALD SUM 方法, 这样的测试还是能有明显提速功效的.
好了, 对于生成EMAP时EWALD 方法的介绍差不多了解这些就足够了. 还有值得提醒的就是, MuSiC中支持的EWALD 方法只是对正交盒子有效 (orthorhombic unit cell). 这个对于GCMC问题不大, zeolite, MOF等材料, 大多数还是orthombic的, 但是如果是自己生成的初始结构, 想要算吸附, 然后又想要使用EMAP的时候, 就需要注意这一点!
Life, Love, Laugh.
66楼2010-11-13 22:11:52
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ghcacj

荣誉版主 (著名写手)

阿超

优秀版主

★ ★ ★
小木虫(金币+0.5):给个红包,谢谢回帖交流
zh1987hs(金币+2):鼓励交流 2010-11-16 16:31:08
顺便问一下,如果一个晶体CIF文件打开后,删除杂原子后,结构比较乱,比如苯环都是扁平的,碰到这样的情况怎么办?用Forcite或Dmol3优化?
67楼2010-11-14 08:46:31
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
ghcacj(金币+10):谢谢 2010-11-14 18:07:08
引用回帖:
Originally posted by ghcacj at 2010-11-14 08:46:31:
顺便问一下,如果一个晶体CIF文件打开后,删除杂原子后,结构比较乱,比如苯环都是扁平的,碰到这样的情况怎么办?用Forcite或Dmol3优化?

这个比较难回答. 实事求是地讲, 很多时候模拟是被动的, 因为模拟结果常常'被'拿来和实验结果做比较, 似乎只有在能验证实验结果基础上的预测才是正确, 合理, 有科学价值的.

但是这里要注意的是, 实验也是有很多假设, 不管是实验方法, 仪器, 还是实验结果的处理. 实验的结果因此也是一个多因素的综合. 比如你这里提到的晶体CIF的问题, 这个文件来源于实验测定, 因此很容易包含杂原子. 那么, 在准备模拟初始构型的时候, 删除了杂原子以后, 那个结构已经不是严格意义上CIF报导的结构了.  那么还需要优化嘛? 我觉得答案取决于你怎样用这个结构, 多数的zeolite, MOF, 因为吸附前后结构的变化不是很大, 大家都倾向于将它们rigid处理. 所以在删除杂原子以后, 对结构是不需要优化的.

当然, 如果确实某些关键的吸附/反应点在有没有杂原子前后, 非常不一样. 那么结构的优化是必须的. 对于那样的体系, 我想光光是确定结构与杂质的对应关系就是非常麻烦了, 往往不能简单放在 Dmol3或者Forcite里面处理. 估计需要系统地跑一些任务来仔细研究.
Life, Love, Laugh.
68楼2010-11-14 13:02:16
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ghcacj

荣誉版主 (著名写手)

阿超

优秀版主

★ ★ ★
小木虫(金币+0.5):给个红包,谢谢回帖交流
zh1987hs(金币+2):谢谢 2010-11-16 16:30:56
那也就是说是否对结构优化,没有一个强制规定,取决于个人的选择,是吗?我感觉吸附模拟中有很多地方比较随性,没有一个统一规定,每个课题组都有自己的一套习惯,因为我看过用DMOL3和FORCITE优化的文献

[ Last edited by ghcacj on 2010-11-14 at 14:46 ]
69楼2010-11-14 14:29:32
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)

★ ★ ★ ★ ★ ★
ghcacj(金币+6):谢谢 2010-11-14 18:07:29
引用回帖:
Originally posted by ghcacj at 2010-11-14 14:29:32:
那也就是说是否对结构优化,没有一个强制规定,取决于个人的选择,是吗?我感觉吸附模拟中有很多地方比较随性,没有一个统一规定,每个课题组都有自己的一套习惯

吸附的理想状态是, 吸附质和吸附剂都是flexable model, 但是有谁这样做呢?

所以, 大概这样的原则: 如果有文献, 那么最开始还是顺着文献大家公认的习惯来做; 如果对一个体系吃透了, 第二步做一些尝试, 比如吸附质的 flexable model 和吸附剂的结构优化等等.
Life, Love, Laugh.
70楼2010-11-14 14:45:54
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 qphll 的主题更新
普通表情 高级回复(可上传附件)
信息提示
请填处理意见