动力学模型已兼容相对能量的输入文件

之前的动力学代码的数据文件我是在catmap的基础上进行改进，写了个tablemaker来生成一个csv文件，然后让要计算动力学的人把每个需要的物种的绝对能量输入到表格中，模型会在处理数据的时候自动计算其所有基元反应的能垒以及反应能量变化。之所有这样进行能量输入，是方便在后面能够基于每个物种的能量进行计算，比如在做出是猜测的时候进行Boltzmann分布，在进行敏感度分析时候能够针对每个中间态和过渡态进行分析。
但是之前有一次帮师姐计算它的反应，他只有相对能量。如下，

反应：

1
2
3
4
5
6
7
8
9
10
% 反应
O2 + # <-> O2#
NO + O2# <-> ONOO#
ONOO# <-> NO2#
O# + NO <-> NO2#
NO2# <-> # + NO2

% 能量
Ea = [0.63 0.65 0.22 0.65 0];
G0 = [0.368 -0.455 -1.13 -0.49 -0.19];

于是我就在纸上手动算出来每个物种的所谓’generalized formation energy’。当然首先把O2、N0以及#的能量为0，然后根据相对能量的关系列出来线性方程组，通过矩阵乘法一步计算出所有其他物种的形成能。(顺便在这里吐槽下，所谓generalized formation energy通过这样的方法也能得到，感觉是个很虚的概念诶～)。得到矩阵的具体方法就是下得到一个coefficients matrix，m是从能量数据中得到的等式数，n是除去能量为0的物种以外其他所有物种。贴上获取矩阵的代码：

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
def get_unknown_coeff_vector(self, elementary_rxn_list):
    """
    Expect a elementary_rxn_list,
    e.g. [['O2_s', 'NO_g'], ['*_s', 'ON-OO_s'], ['*_s', 'ONOO_s']]
    return coefficient vectors, Ea, G0.
    e.g. ([[0, 0, -1, 0, 0, 0, 0, 0, 1], [0, 0, -1, 1, 0, 0, 0, 0, 0]], 0.655, -0.455)
    """
    idx = self._owner.elementary_rxns_list.index(elementary_rxn_list)
    Ea, G0 = self.Ea[idx], self.G0[idx]

    if not hasattr(self, 'unknowns'):
        self.get_unknown_species()

    coeff_vects = []
    if Ea != 0 and len(elementary_rxn_list) != 2:  # has barrier
        #get ts coefficient vector
        is_list, ts_list = elementary_rxn_list[0], elementary_rxn_list[1]
        is_dict, ts_dict = self.list2dict(is_list), self.list2dict(ts_list)
        coeff_vect = []
        for unknown in self.unknowns:
            if unknown in is_dict:
                coeff = -is_dict[unknown]
            elif unknown in ts_dict:
                coeff = ts_dict[unknown]
            else:
                coeff = 0
            coeff_vect.append(coeff)
        coeff_vects.append(coeff_vect)

    #coefficient vector for G0
    is_list, fs_list = elementary_rxn_list[0], elementary_rxn_list[-1]
    is_dict, fs_dict = self.list2dict(is_list), self.list2dict(fs_list)
    coeff_vect = []
    for unknown in self.unknowns:
        if unknown in is_dict:
            coeff = -is_dict[unknown]
        elif unknown in fs_dict:
            coeff = fs_dict[unknown]
        else:
            coeff = 0
        coeff_vect.append(coeff)
    coeff_vects.append(coeff_vect)

    if Ea:
        return coeff_vects, [Ea, G0]
    else:
        return coeff_vects, [G0]

后面的处理就是修改了下parser里面的parse_data()方法，将上面计算得到的形成能解析到动力学模型中。
对应于csv_parser，新的rel_energy_parser中也有相应的parse_data()方法来读取能量数据，
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
def parse_data(self):  # correspond with parse_data() in csv_parser.py
    """
    Solve Axb equation to get value of generalized free energies.
    """
    A, b = [], []
    for rxn_list in self._owner.elementary_rxns_list:
        coeff_vects, value = self.get_unknown_coeff_vector(rxn_list)
        A.extend(coeff_vects)
        b.extend(value)

    A, b = np.matrix(A), np.matrix(b).reshape(-1, 1)
    x = A.I*b  # values for unknowns

    #convert column vector to list
    x = x.reshape(1, -1).tolist()[0]

    #put values to G_dict
    for sp_name, G in zip(self.unknowns, x):
        self.G_dict.setdefault(sp_name, G)

    #put generalized formation energy into species_definition
    for sp_name in self.G_dict:
        sp_dict = self._owner.species_definitions
        sp_dict[sp_name].setdefault('formation_energy', self.G_dict[sp_name])

    setattr(self._owner, 'hasdata', True)

    return


为此我又专门写了一个能够读取相对能量数据并计算动力学模型的运行脚本run_relative.py
能量数据放在一个.py结尾的文件中，形如

Comments

▲