24小时热门版块排行榜    

CyRhmU.jpeg
查看: 6508  |  回复: 101
本帖产生 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这个压力点为例, 贴出所有我计算中用到的文件.
回复此楼

» 收录本帖的淘帖专辑推荐

有机场效应晶体管及有机太阳能电池 MOF相关 第一性原理和电化学 Software Operation
mof 材料

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

Life, Love, Laugh.
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
zh1987hs(金币+4):鼓励 2010-11-04 10:49:21
ghcacj(金币+10, 模拟EPI+1):鼓励,支持 2010-11-04 10:53:28
update 2:

重新生成了CO2的两个pmap文件, 和CuBTC的EMAP文件. 改动如下:

(1) Pmap的control file中 cutoff阀值由原先500KJ/mol 修改为200KJ/mol

(2) EMAP的control file中 cutoff阀值由原先500KJ/mol 修改为100KJ/mol

(3) atm_atm_file中,
CuBTC原子的vdw HICUT为13.1A, LOCUT为0.01A;
CO2原子和CO2-CBTC的vdw HICUT为12.8A, LOCUT为0.01A;

(4) PMAP的迭代次数由原先的1修改为1000000; EMAP的生成方法和原先一样, 还是用一个带单位电荷的探针, 迭代步数依然为1.

在map文件生成的情况下, 单点计算前面提到的那留个压力, 同时在不适用map文件的情况下计算41.28Kpa下的吸附

试验结果如下:

(1) 图一张



尽管我和文献一样, 每点迭代2E7, 取后50%平均, 但是文献是用了3*3*3, 我只是用1*1*1, 这样看起来结果吻合得还算可以.

(2) 同样的压力(41.28KPa), 不用map, 计算结果为2.24kg/mol, 使用PMAP+EMAP后的结果是 2.34kg/mol. 文献报导的改点吸附为 2.22kg/mol.

貌似不用map文件, 计算更加可靠一些? (假设文献报导的值是standard reference)

言下之意, PMAP或者EMAP还有需要改进的地方??? HOW????

(3) MAP文件对计算速度的提升还是很明显的, 依然是(2)中的比较, 不用map文件, 耗时7.1759 hrs, 使用了PMAP+EMAP, 耗时 58.714 min.

注意的是, 不用map文件时, 迭代步数为 1E7; 用Map文件下的计算, 都是 2E7

这样折算下来, 速度差别, 那是相当相当的大........

» 本帖已获得的红花(最新10朵)

Life, Love, Laugh.
28楼2010-11-04 10:40:02
已阅   回复此楼   关注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的回帖
相关版块跳转 我要订阅楼主 qphll 的主题更新
普通表情 高级回复(可上传附件)
信息提示
请填处理意见