PARSER中的反应方程式配平

之所以把反应方程式配平拿出来写是因为之前在这段时间写完solver的过程中发现自己还是对动力学的质量作用定理理解的不够深刻，在之前的setup file中的反应方程式我默认的是以已经配平好的基元反应作为输入的，例如

1
2
3
4
5
6
7
8
9
10
11
12
rxn_expressions = [

                'CO_g + *_s -> CO_s',
                '3H2_g + 6*_s -> 6H_s',
                'CO_s + *_s <-> C-O_s + *_s -> C_s + O_s',
                'O_s + H_s -> HO_s + *_s',
                'HO_s + H_s <-> H-OH_s + *_s -> H2O_g + 2*_s',
                'C_s + H_s <-> H-C_s + *_s -> CH_s + *_s',
                'CH_s + H_s -> CH2_s + *_s',
                'CH2_s + H_s -> CH3_s + *_s',
                'CH3_s + H_s -> CH4_g + 2*_s'
]

然后通过正则表达式解析反应方程式后通过集合做差集运算获得总的反应方程式。具体方法在之前的博客中有写：

获取总的反应方程式并检查是否守恒:

主要是通过get_total_rxn_equation()方法，获取总的反应方程式的方法是利用字典和集合数据类型，通过正则表达式获取初态和终态的所有物种和相应总数的字典，然后遍历所有的基元反应获取两边的字典(包含所有物种的类型和数目)，然后将这两个字典转化成集合，进行相互两次差集的运算获取总反应的两边的物种以及每个物种的数目，最终整合成字符串以化学方程式的形式返回。如果最终的总反应方程不守恒的话则需要用户检查基元反应计量数是否正确。

但后来通过牛顿法解稳态的时候发现解出来的覆盖度解有问题，然后一点点向回退，检查Jacobian Matrix，在检查覆盖度变化率动力学方程组，然后检查到反应速率表达式的时候发现和我手推的时候不一样！这才意识到化学计量数被我之前一直忽略了，这里我曾经有想过化学计量数和质量作用定律的关系，看来要使用质量作用定律还是要按照基元反应来写，因为我在后面计算反应所率的时候会按照前面的化学计量数来统一吸附物的覆盖率变化的。反而在最开始写基元反应表达式的时候引入化学计量数会使后面的反应速率表达式错误。所以输入文件中的基元反应表达式必须写成如下这样：

1
2
3
4
5
6
7
8
9
10
11
12
rxn_expressions = [

                'CO_g + *_s -> CO_s',
                'H2_g + 2*_s -> 2H_s',
                'CO_s + *_s <-> C-O_s + *_s -> C_s + O_s',
                'O_s + H_s -> HO_s + *_s',
                'HO_s + H_s <-> H-OH_s + *_s -> H2O_g + 2*_s',
                'C_s + H_s <-> H-C_s + *_s -> CH_s + *_s',
                'CH_s + H_s -> CH2_s + *_s',
                'CH2_s + H_s -> CH3_s + *_s',
                'CH3_s + H_s -> CH4_g + 2*_s'
]

好，那么问题来了。怎么配平？？这里我搁置了两天，一直在想该怎么写能让程序自动把这些基元反应配平，然后求出总的反应表达式，自己也尝试了很多方法，我尝试修改之前的通过集合运算的方法还是解决不了。正在想放弃这一功能的时候，和师兄交流，师兄是用matlab写代码的，由于矩阵是matlab的基本数据类型，必然很多时候会用矩阵的思想，他说貌似可以用类似矩阵消元的方式来配平。忽然就感觉这个问题有戏，就在纸上随便画了两个矩阵，看了看，我去，可以通过求解系数矩阵的零空间的基的方式把配平系数求出来，系数出来了，配平什么的就是自然而然的了。然后就着手把这个问题写完了，实现的方法(以CO氧化的反应为例)：

基元反应方程式：

1
2
3
4
5
6
7
8
rxn_expressions = [

               '*_s + CO_g -> CO*', 
               '2*_s + O2_g <-> O-O* + *_s -> 2O*',
               'CO* +  O* <-> O-CO* + * -> CO2_g + 2*',

                   ]


获取此反应的系数矩阵，在这里要把吸附物和气体分子分开，因为最后是要把吸附物的系数全部变为0的。
吸附物顺序：$(*,CO^{*}, O^{*}), 气体顺序：( CO, O_{2}, CO_{2})

设$A$为吸附物系数矩阵，$B$为气体系数矩阵，

吸附物系数矩阵
$$
A = \left[\begin{matrix}
1 & -1 & 0 \\
2 & 0 & -2 \\
-2 & 1 & 1
\end{matrix} \right]
$$
气体分子系数矩阵
$$
B = \left[\begin{matrix}
1 & -1 & 0 \\
2 & 0 & -2 \\
-2 & 1 & 1
\end{matrix} \right]
$$

求矩阵$A$的转置的零空间的基，即$A^{T}\bullet x = 0$
将获取的零空间的基转化成整数形式后再左乘矩阵$B$，也就是( $x\bullet B$ )变得到了配平系数向量了.
由于python没有像matlab那样的直接有null函数求零空间的基的函数，利用numpy的函数有写了个：

1
2
3
4
5
def null(A, eps=1e-10):
    "get null space of transposition of site_matrix"
    u,s,vh = np.linalg.svd(A,full_matrices=1,compute_uv=1)
    null_space = np.compress(s <= eps, vh, axis=0)
    return null_space.T

代买实现分别是写了两个方法先获取矩阵，然后在进行矩阵运算，在最后也进行了守恒的检测，贴上来：
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
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
def get_stoichiometry_matrices(self):
    """
    Go through elementary_rxns_list,
    return sites stoichiometry matrix and gas stoichiometry matrix.
    """
    sites_names = ['*_'+site_name 
                   for site_name in self._owner.site_names] + \
                   list(self._owner.adsorbate_names)
    gas_names = list(self._owner.gas_names)
    #initialize matrices
    m = len(self._owner.elementary_rxns_list)
    n_s, n_g = len(sites_names), len(gas_names)
    site_matrix, gas_matrix = np.matrix(np.zeros((m, n_s))),\
                              np.matrix(np.zeros((m, n_g)))
    #go through all elementary equations
    for i in xrange(m):
        states_list = self._owner.elementary_rxns_list[i]
        for sp in states_list[0]:  #for initial state
            stoichiometry, sp_name = self.split_species(sp)
            if sp_name in sites_names:
                j = sites_names.index(sp_name)
                site_matrix[i, j] += stoichiometry
            if sp_name in gas_names:
                j = gas_names.index(sp_name)
                gas_matrix[i, j] += stoichiometry
        for sp in states_list[-1]:  #for final state
            stoichiometry, sp_name = self.split_species(sp)
            if sp_name in sites_names:
                j = sites_names.index(sp_name)
                site_matrix[i, j] -= stoichiometry
            if sp_name in gas_names:
                j = gas_names.index(sp_name)
                gas_matrix[i, j] -= stoichiometry

    return site_matrix, gas_matrix

def get_total_rxn_equation(self):
    "Get total reaction expression of the kinetic model."
    site_matrix, gas_matrix = self.get_stoichiometry_matrices()

    def null(A, eps=1e-10):
        "get null space of transposition of site_matrix"
        u,s,vh = np.linalg.svd(A,full_matrices=1,compute_uv=1)
        null_space = np.compress(s <= eps, vh, axis=0)
        return null_space.T
    x = null(site_matrix.T)  #basis of null space
    x = map(abs, x.T.tolist()[0])
    #convert entries of x to integer
    min_x = min(x)
    x = [round(i/min_x, 1) for i in x]
    setattr(self._owner, 'trim_coeffients', x)
    x = np.matrix(x)
    total_coefficients = (x*gas_matrix).tolist()[0]
    #create total rxn expression
    reactants_list, products_list = [], []
    for sp_name in self._owner.gas_names:
        idx = self._owner.gas_names.index(sp_name)
        coefficient = total_coefficients[idx]
        if coefficient < 0: #for products
            coefficient = abs(int(coefficient))
            if coefficient == 1:
                coefficient = ''
            else:
                coefficient = str(coefficient)
            products_list.append(coefficient + sp_name)
        else: #for reactants
            coefficient = int(coefficient)
            if coefficient == 1:
                coefficient = ''
            else:
                coefficient = str(coefficient)
            reactants_list.append(coefficient + sp_name)
    reactants_expr = ' + '.join(reactants_list)
    products_expr = ' + '.join(products_list)
    total_rxn_equation = reactants_expr + ' -> ' + products_expr

    #check conservation
    states_dict = self.parse_single_elementary_rxn(total_rxn_equation)[0]
    check_result = self.check_conservation(states_dict)
    if not check_result:
        setattr(self._owner, 'total_rxn_equation', total_rxn_equation)
    else:
        if check_result == 'mass_nonconservative':
            raise ValueError('Mass of total equation \''+
                        total_rxn_equation+'\' is not conservative!')
        if check_result == 'site_nonconservative':
            raise ValueError('Site of total equation \''+
                        total_rxn_equation+'\' is not conservative!')

    return total_rxn_equation

效果图：

→_→欲拯救世界必先学线性代数

顺便把之前那个化学式配平的代码也贴上来吧，虽然以后再也用不到了，留个纪念也不错。

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
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
def get_total_rxn_equation_orig(self, elementary_rxns_list):
    """
    Analyse elementary_rxns_list, get total_rxn_equation and 
    check the conservation of it, set total_rxn_equation as attr of model.
    """
    reactants_dict = {}
    products_dict = {}
    for equation_list in elementary_rxns_list:
        #for reactants
        for reactant in equation_list[0]:
            if not '*' in reactant:
                m = self._owner.regex_dict['species'][0].search(reactant)
                if m.group(1):
                    stoichiometry = float(m.group(1))
                else:
                    stoichiometry = 1.0
                species_name = m.group(2)
                site = m.group(3)
                total_name = species_name + '_' + site
            else: #for free site
                m = self._owner.regex_dict['empty_site'][0].search(reactant)
                if m.group(1):
                    stoichiometry = float(m.group(1))
                else:
                    stoichiometry = 1.0
                site_name = m.group(2)
                total_name = '*_' + site_name
            if reactants_dict.has_key(total_name):
                reactants_dict[total_name] += stoichiometry
            else:
                reactants_dict.setdefault(total_name, stoichiometry)
        #for products
        for product in equation_list[-1]:
            if not '*' in product:
                m = self._owner.regex_dict['species'][0].search(product)
                if m.group(1):
                    stoichiometry = float(m.group(1))
                else:
                    stoichiometry = 1.0
                species_name = m.group(2)
                site = m.group(3)
                total_name = species_name + '_' + site
            else: #for free site
                m = self._owner.regex_dict['empty_site'][0].search(product)
                if m.group(1):
                    stoichiometry = float(m.group(1))
                else:
                    stoichiometry = 1.0
                site_name = m.group(2)
                total_name = '*_' + site_name
            if products_dict.has_key(total_name):
                products_dict[total_name] += stoichiometry
            else:
                products_dict.setdefault(total_name, stoichiometry)
    #return sorted(reactants_dict.items()), sorted(products_dict.items())

    reactants_set, products_set = set(reactants_dict.items()), \
                                  set(products_dict.items())
    #return reactants_set, products_set
    def set2str(rxn_set):
        rxn_list = []
        for sp_tuple in rxn_set:
            species_name, number = sp_tuple[0], str(sp_tuple[1])
            if number == '1':
                number = ''
            sp_str = number + species_name
            rxn_list.append(sp_str)
        return ' + '.join(rxn_list)
    total_reactants_str = set2str(reactants_set - products_set)
    total_products_str = set2str(products_set - reactants_set)

    #return total_reactants_str + ' -> ' + total_products_str
    total_rxn_equation = total_reactants_str + ' -> ' + total_products_str

    #check conservation
    states_dict = self.parse_single_elementary_rxn(total_rxn_equation)[0]
    check_result = self.check_conservation(states_dict)
    if not check_result:
        setattr(self._owner, 'total_rxn_equation', total_rxn_equation)
    else:
        if check_result == 'mass_nonconservative':
            raise ValueError('Mass of total equation \''+
                        total_rxn_equation+'\' is not conservative!')
        if check_result == 'site_nonconservative':
            raise ValueError('Site of total equation \''+
                        total_rxn_equation+'\' is not conservative!')

Comments

▲