动力学模型工具类--Parser

这部分其实在放假前就开始写了，后面在写的过程中又添加了前面的table_maker和后面的parser的一些修修补补。在这里整体介绍下这个类的一些特性和简要实现过程。
这个类最初的目的还是很简单，就是解析setup file 和 input file把所有的信息全部解析到模型对象中，其中包括模型所必需的一些属性，例如吸附物，过渡态，气体，表面等信息，能量数据，形成能数据等。。。整个类的大概属性和操作在UML图中可以看到，其中包含一个基类ParserBase类，以及子类CsvParser类，在这里我只写了CsvParser，可以再扩展写用于其它表格的Parser，例如利用xlrd模块直接操作excel。

解析基元反应方程式:
在KineticModel中已经通过load()方法中遍历执行setup file后的局部变量的方式将setup file加载到了model中成为model对象的属性。在setup file中定义了一个species_definitions的字典类型变量，用来存储model中所有物种和位点的信息，例如：

'O_s': {'elements': {'O': 1},
        'formation_energy': 0.597,
        'frequencies': [359.5, 393.3, 507.0],
        'information': 'None',
        'name': 'O',
        'site': 's',
        'type': 'adsorbate'}

在前面table_maker中就是利用的这个变量的遍历生成的表格。species_definitions的生成我主要用了dict类型数据，分别写了几个方法来避免3重以上的循环，不然3重以上循环的代码给别人看估计要被抽了。。。
大致思路还是分治啦，大概介绍下分开的这几个方法的作用：

parse_site_expression()和parse_species_expression()
这两个方法类似，都是解析单个物种，例如'2CH3_s'，生成关于这个物种的dict，例如{'CH3_s': {'number': 2, 'site': 's', 'elements': {'C': 1, 'H':3}}}, 其中化学式的匹配是通过正则表达式的方法，在model模型中分别用compile()将三个常用的正则表达式编译，并设为model对象的属性。

#set elementary parse regex(compiled)
self.regex_dict = {}

states_regex = re.compile(r'([^\<\>]*)(?:\<?\-\>)'+
                      r'(?:([^\<\>]*)(?:\<?\-\>))?([^\<\>]*)')
self.regex_dict['IS_TS_FS'] = [states_regex, ['IS','TS','FS']]

species_regex = re.compile(r'(\d*)([^\_\+\*\<\>]+)_(\w+)')
self.regex_dict['species'] = \
            [species_regex, ['stoichiometry','name','site']]

site_regex = re.compile(r'(\d*)(?:\*\_)(\w+)')
self.regex_dict['empty_site'] = \
            [site_regex, ['stoichiometry', 'site']]

parse_state_expression()
这个就感觉像剥洋葱啦，就是分析包含多个物种的一个状态，例如反应中的一个过渡态`’CH2-H_s + _s’,生成相应的字典和物种列表，也就是个tuple`，例如

{'sp_dict': {'CH-H_s': {'number': 1, 
                        'site': 's', 
                        'elements': {'C': 1, 'H': 2}}}},
             {},
             ['CH2-H_s', '*_s']

parse_single_elementary_rxn()
继续剥洋葱，这次是解析整个基元反应方程式，将一个基元反应
‘C_s + H_s <-> H-C_s + _s -> CH_s + _s’
解析成state_dict 和 rxn equation list的tuple:
<!–hexoPostRenderEscape:

({'FS': {'empty_sites_dict': {'s': {'number': 1, 'type': 's'}},
   'species_dict': {'CH_s': {'elements': {'C': 1, 'H': 1},
     'number': 1,
     'site': 's'}},
   'state_expression': 'CH_s + _s'},
  'IS': {'empty_sites_dict': {},
   'species_dict': {'C_s': {'elements': {'C': 1}, 'number': 1, 'site': 's'},
    'H_s': {'elements': {'H': 1}, 'number': 1, 'site': 's'}},
   'state_expression': 'C_s + H_s'},
  'TS': {'empty_sites_dict': {'s': {'number': 1, 'type': 's'}},
   'species_dict': {'H-C_s': {'elements': {'C': 1, 'H': 1},
     'number': 1,
     'site': 's'}},
   'state_expression': 'H-C_s + _s'}},
 [['C_s', 'H_s'], ['H-C_s', '_s'], ['CH_s', '*_s']])

:hexoPostRenderEscape–>

* parse_elementary_rxns()
最后循环分析所有的基元反应，并将gas_names, adsorbate_names, transition_state_names, elementary_rxns_list赋值到model的属性中。

检验基元反应方程守恒:
这里主要是针对每一个基元反应进行质量守恒和位点守恒进行检测。主要用的了两个方法，一个是获取反应的元素个数的dict，一个是获取位点数的dict：

对于分析元素质量守恒

1
2
3
4
5
6
7
8
9
10
11
12
13
def get_elements_num_dict(self, species_dict):
sum_element_dict = {}
for sp in species_dict:
sp_num = species_dict[sp]['number']
group = {}
#generate a dict e.g. group = {'C': 2, 'O': 2}
for element in species_dict[sp]['elements']:
group.setdefault(element,
species_dict[sp]['elements'][element]*sp_num)
sum_element_dict = \
self.merge_elements_dict(sum_element_dict, group)

return sum_element_dict

这个是分析通过解析以后的得到的state_dict来获取元素数量的dict
例如获取
{'C_s': {'elements': {'C': 1}, 'number': 1, 'site': 's'}, 'CO_s': {'elements': {'C': 1, 'O': 1}, 'number': 2, 'site': 's'}}
得到
{'C': 2, 'O': 1}

对于分析位点守恒

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
def get_total_site_dict(self, state_dict):
total_site_dict = {}
if state_dict['empty_sites_dict']:
for empty_site in state_dict['empty_sites_dict']:
total_site_dict.setdefault(empty_site,
state_dict['empty_sites_dict'][empty_site]['number'])

#get site number from species dict
for species in state_dict['species_dict']:
site = state_dict['species_dict'][species]['site']
if site == 'g': #neglect gas site when check conservation
continue
number = state_dict['species_dict'][species]['number']
if site in total_site_dict:
total_site_dict[site] += number
else:
total_site_dict.setdefault(site, number)
return total_site_dict

这个是同时分析空位和吸附物中来获取总的吸附位点的数量。

获取总的反应方程式并检查是否守恒:
这一部分一开始在写parser的时候没想写，昨天在奶奶家把parser的获取数据并计算形成能的部分写完以后闲着没事就写了，写完了才发现，这个是非常重要的还！主要是通过get_total_rxn_equation()方法，获取总的反应方程式的方法是利用字典和集合数据类型，通过正则表达式获取初态和终态的所有物种和相应总数的字典，然后遍历所有的基元反应获取两边的字典(包含所有物种的类型和数目)，然后将这两个字典转化成集合，进行相互两次差集的运算获取总反应的两边的物种以及每个物种的数目，最终整合成字符串以化学方程式的形式返回。如果最终的总反应方程不守恒的话则需要用户检查基元反应计量数是否正确。
这部分代码贴上来有点占篇幅，有兴趣的可以下载源码。

计算 generalized formation energy 等数据:
这一部分是写在CsvParser的类方法，主要工作是分析由table_maker生成的表格中的DFT_energy这一列数据，然后生成相应的generalized formation energy.最后把formation energy，frequencies， information的数据添加到species_definitions相应的key中。以便后面计算的时候获取能量数据。
generalized formation energy要相对于共同的参考能量:

[latex]E{i} =U{i}-\sum{j}^{}{n{j}R_{j}} [/latex]

其中[latex]E{i} [/latex]是物种i的 generalized formation energy， [latex]U{i}[/latex]是物种i的初始DFT能量，[latex]n{j}[/latex]为是原子j在物种i中的数目，[latex]R{j}[/latex]
为原子j的参考能量。
实现代码：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
def parse_data(self):
"""
Parse input file, add data e.g. formation energy,
frequencies, information, to _owner.species_definition
and set corresponding dict as model's attrs.
"""
file_obj = open('./energy.csv', 'rU')
file_obj.readline()
#formation_energy_dict, frequencies_dict = {}, {}
for line_str in file_obj:
line_list = line_str.strip('\n').split(',')
site_name, species_name, formation_energy,\
frequencies, information =\
line_list[1], line_list[2], float(line_list[4]),\
eval(line_list[5].replace(' ', ',')), line_list[6]
#get full name
#for gas
if site_name == 'gas':
full_name = species_name + '_g'
#for any species on site
else:
#get site symbol
for site_symbol in self._owner.site_names:
if self._owner.species_definitions[site_symbol]\
['site_name'] == site_name:
break
if species_name == 'slab':
full_name = site_symbol
else:
full_name = species_name + '_' + site_symbol
#update _owner.species_definition
for key in ['formation_energy', 'frequencies', 'information']:
#print eval(key)
self._owner.species_definitions[full_name].\
setdefault(key, eval(key))
file_obj.close()

解析以后的甲酸裂解的总反应方程式：

计算formation energy之前的 table:

计算formation energy之后的 table:

parse data之后的 species_definitions:

———————————————————————-我是分割线——————————————————————–
今天是大年初一，祝所有人身体健康，工作顺利。今天红包抢的很开心，嘿嘿～

Prev

Archive

Next

Comments

2015-02-19

学术26

python57

kinetic model10

kinetics9

Kynetix15

catalysis28

chemistry25