【求助】CuBTC MOF + CO2 MuSiC 问题 - 分子模拟 - Monte Carlo - 第4页小木虫论坛-学术科研互动平台

31楼2010-11-04 10:59:01

qphll

金虫 (正式写手)

模拟EPI: 10
应助: 18 (小学生)
金币: 2708.7
散金: 2294
红花: 8
帖子: 554
在线: 226.6小时
虫号: 89654
注册: 2005-08-29
性别: GG
专业: 碳素材料与超硬材料

引用回帖:

Originally posted by ghcacj at 2010-11-04 10:51:55:

不是Map需要改进，而是你制作map的时候，迭代步数可能会有影响，另外一个影响计算结果微弱差别的地方，应该是计算过程中输出的构型数，输出的数目越多，统计的效果就会越好。在同样的步数下，输出的构型数目也会 ...

说得有道理. VDW cutoff的调整, 也是我其中的一个尝试. 回头再来update, 先做会正事....

赞一下

Life, Love, Laugh.

32楼2010-11-04 11:02:01

qphll

金虫 (正式写手)

模拟EPI: 10
应助: 18 (小学生)
金币: 2708.7
散金: 2294
红花: 8
帖子: 554
在线: 226.6小时
虫号: 89654
注册: 2005-08-29
性别: GG
专业: 碳素材料与超硬材料

★
zh1987hs(金币+1):谢谢 2010-11-05 22:37:04

引用回帖:

Originally posted by ghcacj at 2010-11-04 10:59:01:
提议一下，可否将这个计算所涉及到的文件打包发在网盘上呢（map不需要，可以自己制作）让大家都学习下。

好的, 等我将这个例子琢磨透了, 会将所有文件打包上传.

Life, Love, Laugh.

33楼2010-11-04 11:06:33

xiaokongyu

木虫 (著名写手)

应助: 2 (幼儿园)
金币: 2233.7
散金: 35
红花: 6
帖子: 1294
在线: 253.3小时
虫号: 1118520
注册: 2010-10-10
性别: MM
专业: 污染生态化学

终于。。。抢到沙发

新手，请多多关照！谢谢啦！

34楼2010-11-04 11:39:39

qphll

金虫 (正式写手)

模拟EPI: 10
应助: 18 (小学生)
金币: 2708.7
散金: 2294
红花: 8
帖子: 554
在线: 226.6小时
虫号: 89654
注册: 2005-08-29
性别: GG
专业: 碳素材料与超硬材料

引用回帖:

Originally posted by xiaokongyu at 2010-11-04 11:39:39:
终于。。。抢到沙发

楼上来打酱油, 我来再update一下~~

Life, Love, Laugh.

35楼2010-11-05 00:13:28

qphll

金虫 (正式写手)

模拟EPI: 10
应助: 18 (小学生)
金币: 2708.7
散金: 2294
红花: 8
帖子: 554
在线: 226.6小时
虫号: 89654
注册: 2005-08-29
性别: GG
专业: 碳素材料与超硬材料

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
ghcacj(金币+10):谢谢 2010-11-05 09:17:58

update 4

上回说到在生成PMAP和EMAP的时候, 我尝试了不同的High Energy CutOFF和迭代步数.

严格意义上讲, 在生成PAMP和EMAP的时候, 是不需要做迭代的, MAP的意义就是pre-define Energy Grid, 然后在计算的时候, 可以查表, 加速计算. 可怜的MC由于方法上的原因, 并行计算不那么容易, 所以只有想着办法来加速计算. 当然, MD计算的时候, 也可以生成这样的map文件来加速的, 但是这样带来的提速和并行优化比起来, 简直就是你饿着肚子骑自行车和我国和谐号提速的差别 ^-^

回到这边来, 我尝试了四组计算, 由于前面我发现EMAP的变化对最后结果的影响比较大. 这里主要考虑了生成EMAP的参数影响, namely, 迭代步数 1 vs. 100W; HighEnergeCutOFF 100KJ/mol vs. 200KJ/mol.

上图...

大家可以看到, 迭代步数对结果几乎没有影响 (Set1-Set4 and Set2-Set3), 说几乎是因为Set2-Set3最后一个数据有些差别, 当然, 这也有可能是计算步数不够, 或者统计的image数不多造成的.

CUTOFF对结果还是有影响的, 而且蛮大.

从这个测试得到两个信息:

(1) 小心对待map文件中的CUTOFF
(2) 生成map时候的迭代步数不那么重要.
从生成map文件的速度来看, 迭代步数为1和100W是一样的, 我还没有查源程序, 但是有可能迭代步数的增加并没带来任何影响, 源程序可以将1和100W相同处理.

洗洗睡了...

Life, Love, Laugh.

36楼2010-11-05 00:48:12

qphll

金虫 (正式写手)

模拟EPI: 10
应助: 18 (小学生)
金币: 2708.7
散金: 2294
红花: 8
帖子: 554
在线: 226.6小时
虫号: 89654
注册: 2005-08-29
性别: GG
专业: 碳素材料与超硬材料

早起的虫子被鸟吃.

赞一下

Life, Love, Laugh.

37楼2010-11-05 08:06:24

qphll

金虫 (正式写手)

模拟EPI: 10
应助: 18 (小学生)
金币: 2708.7
散金: 2294
红花: 8
帖子: 554
在线: 226.6小时
虫号: 89654
注册: 2005-08-29
性别: GG
专业: 碳素材料与超硬材料

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
ghcacj(金币+20):谢谢 2010-11-05 09:18:10

update 5

前面update 4中说到, 在生成map文件的时候, 迭代次数似乎没有起到什么作用, 从计算时间上来看, 1和100W次的差别, 如果需要统计平均, 那么应该还是能感觉得到的. 但是事实是: 没有差别.

那么为什么会这样? 打一个砂锅看看......

哦, 看到纹路了. 原来是这样的:

首先, 我们看到在生成map的控制文件中, 有可能修改的是这样的两个section:

------ General Information ------------------------------------------
CO2 in CuBTC, to generate C_CO2 Pamp@CuBTC
1 # No. of iteration
1
1
1
1
1
1
C_CO2.CuBTC.pmap.res
C_CO2.CuBTC.pmap.con

或者

------ Mapmaker Information  -----------------------------------------------
1       # Number of maps to make
      #  PLEASE KEEP THIS 'BLANK LINE'!!
CuBTC    # Name of the sorbent
CO2_C_prob_pmap    # Name of prob atoms of the sorbate
NCOUL LJ  # Typer of interaction to be mapped
0.1    # Grid of spacing, Angstrom
200.0    # High end potential cutoff, KJ/mol
AUTO    # Map filename or AUTO

我前面做的修改迭代次数是将General Information里面的第二行, (第一个“1”) 修改为1000000

这里我做了这样的测试:

(1) 将General Information修改为如下:

------ General Information ------------------------------------------
CO2 in CuBTC, to generate C_CO2 Pamp@CuBTC
TEST # No. of iteration
TEST
TEST
TEST
TEST
TEST
TEST
C_CO2.CuBTC.pmap.res
C_CO2.CuBTC.pmap.con

文件的其他部分不做任何改变. 你猜怎地?

mapmake 照样跑得很欢:

Map progress:    1000 points calculated  0.01% of total
Map progress:    2000 points calculated  0.01% of total
Map progress:    3000 points calculated  0.02% of total
Map progress:    4000 points calculated  0.02% of total
Map progress:    5000 points calculated  0.03% of total
Map progress:    6000 points calculated  0.03% of total
Map progress:    7000 points calculated  0.04% of total
Map progress:    8000 points calculated  0.04% of total
Map progress:    9000 points calculated  0.05% of total
Map progress:    10000 points calculated  0.05% of total
Map progress:    11000 points calculated  0.06% of total
Map progress:    12000 points calculated  0.07% of total

半个苹果下去了, 还没有任何错误地继续爬格子下去... 赶紧将任务停了.

这说明General Information里面的数值设置, 对生成map文件没有任何作用.

(2) 来个狠的. 将控制文件的其他section全部删除, 只留下这个:

------ Mapmaker Information  -----------------------------------------------
1       # Number of maps to make
      #  PLEASE KEEP THIS 'BLANK LINE'!!
CuBTC    # Name of the sorbent
CO2_C_prob_pmap    # Name of prob atoms of the sorbate
NCOUL LJ  # Typer of interaction to be mapped
0.1    # Grid of spacing, Angstrom
200.0    # High end potential cutoff, KJ/mol
AUTO    # Map filename or AUTO

毫无意外, 程序报错, 说缺少某某Section之类的. 看来某些Section起不到实际作用, 但是没了它们, 还是不行的.
(等等, 这怎么听起来像是导师和学生的关系呢? ^-^)

(3) 刨根一下, 看看源文件 mapmaker.F90

哦, 原来这样....

那个文件中但是还真有一个DO循环, 也确实是针对 number of mapfile的. 莫非这个mapmake倒真可以并行起来,同时生成几个map文件?

试试吧.

(3)-1

Search一下这个Do循环的ID:

mmparams%nmaps

还好, 没有多少行, 看起来也是在读入control文件的一些信息, 并且做一些预判断.

(3)-2

修改一下试试看能不能给我同时生成两个map文件. 修改的Section如下, 其他未提及的保持原样.

------ Mapmaker Information  -----------------------------------------------
2       # Number of maps to make
      #  PLEASE KEEP THIS 'BLANK LINE'!!
CuBTC    # Name of the sorbent
CO2_C_prob_pmap    # Name of prob atoms of the sorbate
NCOUL LJ  # Typer of interaction to be mapped
0.1    # Grid of spacing, Angstrom
200.0    # High end potential cutoff, KJ/mol
AUTO    # Map filename or AUTO
      #  PLEASE KEEP THIS 'BLANK LINE'!!
CuBTC    # Name of the sorbent
CO2_C_prob_pmap    # Name of prob atoms of the sorbate
NCOUL LJ  # Typer of interaction to be mapped
0.1    # Grid of spacing, Angstrom
100.0    # High end potential cutoff, KJ/mol
AUTO    # Map filename or AUTO

这里的意图是让程序给我生成不同CUTOFF的两个PMAP.

首先得到的报错是:

config.F90: Error: Could not allocate memory for grid%coords

好, 大意是内存不能Allocate, 也许格点太多, 内存不够.

Double check了一下config.F90, 也没有啥其他的发现. 于是修改了一下Grid of sapceing, 将原来的数值从0.1改为10.0
以此来减少生成的格点数目.

修改后的样子如下:

------ Mapmaker Information  -----------------------------------------------
2       # Number of maps to make
      #  PLEASE KEEP THIS 'BLANK LINE'!!
CuBTC    # Name of the sorbent
CO2_C_prob_pmap    # Name of prob atoms of the sorbate
NCOUL LJ  # Typer of interaction to be mapped
10.0    # Grid of spacing, Angstrom
200.0    # High end potential cutoff, KJ/mol
AUTO    # Map filename or AUTO
      #  PLEASE KEEP THIS 'BLANK LINE'!!
CuBTC    # Name of the sorbent
CO2_C_prob_pmap    # Name of prob atoms of the sorbate
NCOUL LJ  # Typer of interaction to be mapped
10.0    # Grid of spacing, Angstrom
100.0    # High end potential cutoff, KJ/mol
AUTO    # Map filename or AUTO

程序倒真能运行了. 但是会“发现”map文件的文件名已经存在.
不管怎样修改, 我都没有解决这个’误会’就只好作罢. 我想大概程序作者只是有这样的规划, 但是也没有能真正在程序中实现吧.

总结一下, 这个update 5 告诉我们, 生成map的时候, 不管是PAMP还是EMAP, 迭代的次数是不需要担心的.

[ Last edited by qphll on 2010-11-5 at 08:17 ]

Life, Love, Laugh.

38楼2010-11-05 08:06:49

ghcacj

荣誉版主 (著名写手)

阿超

模拟EPI: 16
应助: 26 (小学生)
贵宾: 2.764
金币: 2543.5
散金: 8899
红花: 35
沙发: 3
帖子: 1376
在线: 1198.9小时
虫号: 513281
注册: 2008-02-27
性别: GG
专业: 理论和计算化学
管辖: 分子模拟

★ ★ ★
小木虫(金币+0.5):给个红包，谢谢回帖交流
zh1987hs(金币+2):好精彩的讨论啊~ 2010-11-05 22:39:31

好吧，算我认识不足，不过在测试中，你忽略了一个因素，这个因素可能是你计算结果和文献中结果差异的原因，那就是逸度，你用的是压力作为输入参数，而实际上应该用逸度，对于室温下1MPa甚至更高压力下的CO2，逸度和压力的差别已经不可忽略了。

赞一下(2人)

39楼2010-11-05 09:20:02

qphll

金虫 (正式写手)

模拟EPI: 10
应助: 18 (小学生)
金币: 2708.7
散金: 2294
红花: 8
帖子: 554
在线: 226.6小时
虫号: 89654
注册: 2005-08-29
性别: GG
专业: 碳素材料与超硬材料

★ ★ ★ ★ ★
zh1987hs(金币+5):谢谢 2010-11-05 22:39:08

引用回帖:

Originally posted by ghcacj at 2010-11-05 09:20:02:
好吧，算我认识不足，不过在测试中，你忽略了一个因素，这个因素可能是你计算结果和文献中结果差异的原因，那就是逸度，你用的是压力作为输入参数，而实际上应该用逸度，对于室温下1MPa甚至更高压力下的CO2，逸度 ...

这个不是问题, 我测试的压力都不算高. 事实上我注意到这个问题,
从文献读出的是压力值, 在该温度下, 压力和逸度的转换(RP方程) 如下, 单位均为KPa

压力       逸度
11.55    11.54
41.28    41.19
120.67    119.92
321.05    315.43
530.09    514.82
720.96    692.86
1149.34 1078.07

也许我没有看源代码看清楚, 似乎如果在ControlFile中指定从文件中读取的话, 比如指定pressure.dat, 那么在该文件里面放置逸度数值; 但是如果是直接在ControlFile中提供, 那么就提供压力. 似乎源码里面读取ControlFile的值以后, 也会做一个换算.

我有些记不得了, 得重新确认一下.

Life, Love, Laugh.

40楼2010-11-05 09:40:53