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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 分子模拟 » Monte Carlo » 【求助】CuBTC MOF + CO2 MuSiC 问题
555/6返回列表上一页123456下一页
查看: 6651  |  回复: 101
查看全部回帖 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
本帖产生 3 个 模拟EPI ，点击这里进行查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

qphll

金虫 (正式写手)

[交流] 【求助】CuBTC MOF + CO2 MuSiC 问题 已有12人参与

在用MuSiC测试体系: CO2+CuBTC的时候, 发现了两个问题:
(1) 采用与文献相同模拟, 相同参数, 相同软件的情况下, 吸附量差别很大;
(2) 改变压力做吸附等温线, 不同压力下的吸附量竟然差别不大.

应该是哪里的设置有些问题. 让我细细道来, 大家评论一下.

拿来测试的文献是: Molecular Simulations and Experimental Studies of CO2, CO and N2 adsorption in Metal-Organic Frameworks, J. Phys. Chem. C, 114, 15735, 2010

更具体一点, 试图重复的数据是该文献Figure-2中, 298K下CO2在CuBTC中的吸附. 在GetData的帮助下, Figure-2中CO2在CuBTC的吸附量大致如下:

Pressure, KPa              Loading, mol/kg
41.28                          2.22
120.67                        5.40
321.05                       10.93
530.09                       12.94
720.96                       14.03
1149.34                      14.95

以41.28KPa这个压力点为例, 贴出所有我计算中用到的文件.
回复此楼

» 猜你喜欢
Life, Love, Laugh.
1楼 2010-11-03 14:13:26
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)
★ ★ ★ ★ ★ ★ ★ ★ ★ ★
ghcacj(金币+10):谢谢 2010-11-11 14:11:03
我应该还有一个update没有做, 就是有关能量的波动曲线. 今天是没有时间做那个事情了. 但是现在如果回头看看做的这些的尝试, 有多少事实上是不需要的呢, 或者有哪些是在计算之前就需要弄清楚的呢?

(1) 体系的大小
(2) VDW参数
(3) CUTOFF数值, VDW, COUL (WFCOUL or EWALD ?)
(4) 计算步数, 能量曲线, 接受概率
(5) MC MOVE的步长
(6) 文献 (实验, 计算)

注意:
(a) 如果要修改CuBTC 为 2*2*2, 很有可能会遇到报错, 说有too many spec-spec interaction. 需要修改ssbasic.F90 Line 1532, 将那里的体系限制从10000改大.

(b) 在 CuBTC 为 2*2*2的情况下, map文件的生成是一个问题, 不管是emap, 还是pamp, 除了需要修改一些参数, 例如, ewald.F90中是hard coded to 50000. 你需要增大这个阀值. 但是, 增大了也没有用. 因为你的内存很有可能是不够的, 我大概估算了一下, 在 0.1 angstrom 为步长的情况下, map文件的生成, 似乎需要大于5G的内存分配. 或者我疏漏了什么, 但是有谁做过这个尝试嘛???

不管怎样, 将这些问题放在心里, 那么你就可以在重现文献数据, 练手熟练的情况下, 放手尝试新的体系了, 而且会对你自己的计算结果有信心!

(欠一个能量图的分析......)
赞一下(2人) 回复此楼
Life, Love, Laugh.
60楼2010-11-11 13:58:15
已阅   回复此楼   关注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 ?

^-^
赞一下(1人) 回复此楼
Life, Love, Laugh.
64楼2010-11-12 04:43:25
已阅   回复此楼   关注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的时候, 就需要注意这一点!
赞一下(2人) 回复此楼
Life, Love, Laugh.
66楼2010-11-13 22:11:52
已阅   回复此楼   关注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里面处理. 估计需要系统地跑一些任务来仔细研究.
赞一下(1人) 回复此楼
Life, Love, Laugh.
68楼2010-11-14 13:02:16
已阅   回复此楼   关注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 和吸附剂的结构优化等等.
赞一下(1人) 回复此楼
Life, Love, Laugh.
70楼2010-11-14 14:45:54
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)
引用回帖:
Originally posted by ghcacj at 2010-11-22 15:39:10:
是啊，没有下文了，可能LZ最近比较忙

三个老师在组里访问. 过去的一周, 没有做成任何事情....

Next week will be better? :-)
赞一下 回复此楼
Life, Love, Laugh.
75楼2010-11-23 00:43:47
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)
★ ★ ★ ★ ★ ★ ★ ★ ★ ★
ghcacj(金币+10):精彩 2010-11-23 09:19:17
引用回帖:
Originally posted by zyj8119 at 2010-11-22 13:58:47:

没有下文了？

小update一下.

常常看到文献上有这样的说法:

"At least 20 million trials were used in the simulations, among which the first half was used for equilibrium, and the last half was used to calculate the ensemble averages."

那么这里的"20 million trials" 和 "50%" 这两个关键词是怎样来的?

同样用我们这里的例子, 同样用类似前面某次update提到的体系能量变化图, 我们再仔细看看:

上图!



大家可以看到, 我这里也是跑了 20 million, 也是将前50%扔掉, 后50%系综平均.

但是很明显, 从大图(黑线)可以看出, 体系的总能量在(10E6~14E6)这个阶段还是有很明显的降低. 从小图1(红线, 18E6~20E6)看出, 就算是将总迭代(2E7)的80%扔掉, 体系的能量在最后的20%那里还是在降低.

从小图2(蓝线, 19E6~20E6) 看起来, 似乎体系还是在波动, 但是考虑到波动的绝对值小于 20KJ (注意, 这里是 20KJ for the whole system, NOT 20KJ/mol), 再结合大图(黑线)最后 19E6~20E6处体现出来的'平坦', 我们有理由说服自己, 在19E6步以后, 体系已经达到平衡态了. 差不多可以安全地做系综平均了.

但是, 为了百分百放心, 最好在这个基础上, 再往后算10E6步, 然后看看这个区间 19E6~30E6 的能量波动, 如果能量的波动像小图2(蓝线)所示范围震荡, 那么OK, 放心大胆系综平均, 得出你的结论吧. 反之, 如果体系能量还是像小图1(红线)所示, 有一路下降的趋势(哪怕斜率平坦), 那么对不起, 体系还没有平衡, 继续跑吧......


以上小小update, 可以看作如何选择迭代总步数, 如果选择系综平均(到底扔掉前面的百分之多少)的一些依据.

[ Last edited by qphll on 2010-11-23 at 07:02 ]
赞一下(2人) 回复此楼
Life, Love, Laugh.
76楼2010-11-23 06:56:08
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)
★ ★ ★
zh1987hs(金币+3): 谢谢 2010-11-23 12:46:03
引用回帖:
Originally posted by ghcacj at 2010-11-23 09:20:27:


不过目前我看到过的文献，计算步数最大一般是4E7了，似乎很难有更大的，这个感觉不光和计算体系有关，也和硬件配置相关，几年前几百万步一样可以被接受，现在似乎一定要1E7以上了。这个就和计算能接受的误差有 ...

几年前, 实验对于表面,界面和受限孔内的现象几乎束手无策, 理论计算怎么说怎么有理. 现在仪器和实验技术的发展, 要求理论计算能提供更加精细的信息, 不单单满足于往左还是往右这样的趋势, 而且要求能说出往左多少度这样的信息. 所以计算精度要求的提高,也是势在必行.
赞一下(1人) 回复此楼
Life, Love, Laugh.
79楼2010-11-23 12:42:57
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)
★ ★ ★
ghcacj(金币+3):谢谢 2010-11-23 14:41:16
引用回帖:
Originally posted by ghcacj at 2010-11-23 13:56:39:
很好，LZ上次说模拟Cu-BTC吸附NH3的例子，有文献吗？

没. 光NH3本身的力场就是问题, 它和H2O的力场一样, 力场让人头疼. 值得做得东西还有就是将MOF的力场精细化, 一般采用的UFF太粗略了. 而且rigid处理总是不很恰当, 尤其是要研究分子在MOF里面的动态行为的时候, 比如扩散.

[ Last edited by qphll on 2010-11-23 at 14:25 ]
赞一下(1人) 回复此楼
Life, Love, Laugh.
81楼2010-11-23 14:21:33
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 qphll 的主题更新
555/6返回列表上一页123456下一页
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
信息提示
关闭
请填处理意见
关闭