# 写在前面
在进行动力学求解实操前，请仔细浏览官方网站的学习资料，避免不必要错误发生。源码十分重要，github上有全部源码，尽可能多读源码，有助于debug和自定义实现目标功能和模块。
官网：https://catmap.readthedocs.io/en/latest/installation.html
github: https://github.com/SUNCAT-Center/catmapCatMAP主要需要两个输入文件，分别是**能量文件和机理文件**能量文件
以下是CatMAP官网给的准备能量文件的相关代码
from ase.atoms import string2symbols
abinitio_energies = {
'CO_gas': -626.611970497,
'H2_gas': -32.9625308725,
'CH4_gas': -231.60983421,
'H2O_gas': -496.411394229,
'CO_111': -115390.445596,
'C_111': -114926.212205,
'O_111': -115225.106527,
'H_111': -114779.038569,
'CH_111': -114943.455431,
'OH_111': -115241.861661,
'CH2_111': -114959.776961,
'CH3_111': -114976.7397,
'C-O_111': -115386.76440668429,
'H-OH_111': -115257.78796158083,
'H-C_111': -114942.25042955727,
'slab_111': -114762.254842,
}
ref_dict = {}
ref_dict['H'] = 0.5*abinitio_energies['H2_gas']
ref_dict['O'] = abinitio_energies['H2O_gas'] - 2*ref_dict['H']
ref_dict['C'] = abinitio_energies['CH4_gas'] - 4*ref_dict['H']
ref_dict['111'] = abinitio_energies['slab_111']
def get_formation_energies(energy_dict,ref_dict):
formation_energies = {}
for key in energy_dict.keys(): #iterate through keys
E0 = energy_dict[key] #raw energy
name,site = key.split('_') #split key into name/site
if 'slab' not in name: #do not include empty site energy (0)
if site == '111':
E0 -= ref_dict[site] #subtract slab energy if adsorbed
#remove - from transition-states
formula = name.replace('-','')
#get the composition as a list of atomic species
composition = string2symbols(formula)
#for each atomic species, subtract off the reference energy
for atom in composition:
E0 -= ref_dict[atom]
#round to 3 decimals since this is the accuracy of DFT
E0 = round(E0,3)
formation_energies[key] = E0
return formation_energies
formation_energies = get_formation_energies(abinitio_energies,ref_dict)
for key in formation_energies:
print key, formation_energies[key]然而当你运行复现时可能出现报错，提示无法使用string2symbols，此时可更换成fomular，实现代码复现。若出现其他报错，请根据报错信息进行相应修改和调整。具体如下所示
from ase.formula import Formula
abinitio_energies = {
'CO_gas': -626.611970497,
'H2_gas': -32.9625308725,
'CH4_gas': -231.60983421,
'H2O_gas': -496.411394229,
'CO_111': -115390.445596,
'C_111': -114926.212205,
'O_111': -115225.106527,
'H_111': -114779.038569,
'CH_111': -114943.455431,
'OH_111': -115241.861661,
'CH2_111': -114959.776961,
'CH3_111': -114976.7397,
'C-O_111': -115386.76440668429,
'H-OH_111': -115257.78796158083,
'H-C_111': -114942.25042955727,
'slab_111': -114762.254842,
}
ref_dict = {
'H': 0.5 * abinitio_energies['H2_gas'],
'O': abinitio_energies['H2O_gas'] - 2 * 0.5 * abinitio_energies['H2_gas'],
'C': abinitio_energies['CH4_gas'] - 4 * 0.5 * abinitio_energies['H2_gas'],
'111': abinitio_energies['slab_111']
}
def get_formation_energies(energy_dict, ref_dict):
formation_energies = {}
for key in energy_dict.keys():
E0 = energy_dict[key] # raw energy
name, site = key.split('_') # split key into name/site
if 'slab' not in name: # do not include empty site energy (0)
if site == '111':
E0 -= ref_dict[site] # subtract slab energy if adsorbed
formula = name.replace('-', '') # remove - from transition-states
composition = Formula(formula).count() # get the composition as a dict of atomic species
for atom, count in composition.items():
E0 -= ref_dict[atom] * count # subtract off the reference energy for each atom
formation_energies[key] = round(E0, 3) # round to 3 decimals since this is the accuracy of DFT
return formation_energies
formation_energies = get_formation_energies(abinitio_energies, ref_dict)
for key, value in formation_energies.items():
print(key, value)关于生成能的参考态的选择是任意的，因为参考能量|Rj|在任何相对量中都会抵消，因此对于气相能量的选择可根据自身体系确定。运行以上代码后可检查生成能中参考态的能量是否为0进一步确认。
实操案例，仅作参考
from ase.formula import Formula
abinitio_energies = {
'CO_gas': -14.77871128,
'H2_gas': -6.77185318,
'O2_gas': -8.75496776,
'CO_100': -615.02843792,
'COHH_100': -621.85821680,
'CHO-H_100': -621.76479902,
'CH2O_100': -622.76728347,
'HCH3O_100': -630.7248231,
'CH3OH_100': -630.77961136,
'CH3OH_gas': -30.22753676,
'CO-HH_100': -621.02512236,
'HH-CO_100': -621.74680964,
'H-CHO_100': -621.15969570,
'H-HCH2O_100': -629.29653002,
'CH3O-H_100':-630.71280466,
'slab_100': -599.96338123,
}
ref_dict = {
'H': 0.5 * abinitio_energies['H2_gas'],
'O': 0.5 * abinitio_energies['O2_gas'],
'C': abinitio_energies['CO_gas'] - 0.5 * abinitio_energies['O2_gas'],
'100': abinitio_energies['slab_100']
}
def get_formation_energies(energy_dict, ref_dict):
formation_energies = {}
for key in energy_dict.keys():
E0 = energy_dict[key] # raw energy
parts = key.split('_') # split key into parts
if len(parts) > 1: # check if there is a site identifier
name, site = parts # unpack if there are exactly two parts
else:
continue # skip keys without a site identifier
if 'slab' not in name: # do not include empty site energy (0)
if site == '100':
E0 -= ref_dict[site] # subtract slab energy if adsorbed
formula = name.replace('-', '') # remove - from transition-states
composition = Formula(formula).count() # get the composition as a dict of atomic species
for atom, count in composition.items():
E0 -= ref_dict.get(atom, 0) * count # subtract off the reference energy for each atom
formation_energies[key] = round(E0, 3) # round to 3 decimals since this is the accuracy of DFT
return formation_energies
formation_energies = get_formation_energies(abinitio_energies, ref_dict)
for key, value in formation_energies.items():
print(key, value)继续运行以下代码，生成能量文件
frequency_dict = {
'CO_gas': [2127.2],
'H2_gas': [4336.4],
'O2_gas': [1551.2],
'CO_100': [2071.2, 258.0, 140.8, 91.2, 52.3, 34.2],
'COHH_100': [3550.2, 2056.4, 1724.8, 873.8, 803.3, 565.7, 399.9, 268.4, 172.6, 151.6, 48.4,57.4],
'CHO-H_100': [3529.5, 2505.9, 1687.1, 1267.2, 748.2, 697.8, 541.1, 481.9, 243.6, 168.4, 35.7, 19.4],
'CH2O_100': [2988.1, 2863.2, 1666.7, 1478.4, 1251.9, 1178.3, 335.6, 240.5, 169.4, 164.2, 99.5, 28.8],
'HCH3O_100': [3003.6, 2959.6, 2844.2, 2139.7, 1447.7, 1436.2, 1419.3, 1240.9, 1153.3, 1134.8, 1122.2, 1013.9, 386.5, 224.7, 196.6, 128.5, 87.3, 58.6],
'CH3OH_100': [3074.4, 3022.8, 2919.3, 2520.9, 1490.9, 1447.6, 1438.7, 1418.2, 1139.5, 1134.1, 1017.0, 976.9, 295.3, 232.5, 174.4, 150.7, 124.4, 80.1],
'CH3OH_gas': [3775.4, 3054.5, 2974.6, 2926.2, 1457.6, 1448.0, 1427.1, 1325.0, 1134.2, 1055.9, 997.3, 286.0],
'CO-HH_100': [2036.5, 1710.7, 1379.7, 1225.6, 581.1, 403.5, 255.1, 223.8, 143.2, 44.9, 28.2],
'HH-CO_100': [3535.0, 1862.9, 1017.6, 999.9, 771.0, 613.8, 583.0, 398.3, 219.2, 127.0, 95.7],
'H-CHO_100': [3121.8, 2680.2, 1531.1, 1280.6, 923.2, 774.4, 599.5, 284.4, 169.7, 118.7, 28.7],
'H-HCH2O_100': [3017.6, 2918.8, 2831.7, 1535.4, 1407.0, 1215.2, 1174.8, 1060.9, 1028.5, 891.3, 830.2, 553.0, 316.7, 205.9, 165.1, 103.0, 68.6],
'CH3O-H_100': [3028.1, 2983.9, 2873.1, 1504.7, 1449.2, 1438.5, 1428.2, 1292.8, 1141.4, 1130.1, 1013.6, 504.3, 359.2, 205.7, 157.9, 138.8, 55.4],
}
def make_input_file(file_name, energy_dict, frequency_dict):
# create a header
header = '\t'.join(['surface_name', 'site_name',
'species_name', 'formation_energy',
'frequencies', 'reference'])
lines = [] # list of lines in the output
for key in energy_dict.keys(): # iterate through keys
E = energy_dict[key] # raw energy
name, site = key.split('_') # split key into name/site
if 'slab' not in name: # do not include empty site energy (0)
frequency = frequency_dict[key]
if site == 'gas':
surface = None
else:
surface = 'ZnGa2O4'
outline = [surface, site, name, E, frequency, 'Input File Tutorial.']
line = '\t'.join([str(w) for w in outline])
lines.append(line)
lines.sort() # The file is easier to read if sorted (optional)
lines = [header] + lines # add header to top
input_file = '\n'.join(lines) # Join the lines with a line break
with open(file_name, 'w') as input: # use with statement to handle file opening and closing
input.write(input_file) # write the text
print('Successfully created input file') # correct print statement in Python 3
file_name = 'energies.txt'
make_input_file(file_name, formation_energies, frequency_dict)机理文件
机理文件通常以.mkm尾缀命名。由于体系差异性大，机理文件中可自定义的内容较多，在此仅给出基础必须内容和相关案例，请以官方学习资料为准。
#以温度/压力作为描述符
scaler = 'ThermodynamicScaler'
# 此部分为微观动力学中涉及到的基元反应，其中_s表示催化剂表面位点，若考虑不同吸附位点，可以用其他表达式加以区分，如_t;_h。
rxn_expressions = [
'N2_g + 2*_s <-> N-N* + *_s -> 2N*',
'H2_g + 2*_s -> 2H*',
'N* + H* <-> N-H* + *_s -> NH* + *_s',
'NH* + H* <-> NH-H* + *_s -> NH2* + *_s',
'NH2* + H* <-> NH2-H* + *_s -> NH3* + *_s',
'NH3* -> NH3_g + *_s',
]
#指明催化剂表面，此处需与能量文件中的信息保持一致
surface_names = ['Ni']
descriptor_names= ['temperature','logPressure']
descriptor_ranges = [[400,1000],[-8,3]]
#若需要更多的数据点可适当增大，如29等，相应的计算时长增加
resolution = 15
species_definitions = {}
species_definitions['N2_g'] = {'concentration':1./4.}
species_definitions['H2_g'] = {'concentration':3./4.}
species_definitions['NH3_g'] = {'concentration':0}
#定义反应活性位点，也需与能量文件中的保持一致
species_definitions['s'] = {'site_names': ['111'], 'total':1}
scaling_constraint_dict = {
'N_s':['+',0,None],
'H_s':['+',0,None],
'NH_s':['+',0,None],
'NH2_s':['+',0,None],
'NH3_s':['+',0,None],
'N-N_s':'final_state',
'N-H_s':'final_state',
'NH-H_s':'final_state',
'NH2-H_s':'final_state',
}
data_file = 'furfural_hydro.pkl'
#指明能量文件的路径和名称
input_file = 'energies.txt'
#对气体进行何种校正，具体区别请查官网
#gas_thermo_mode = "shomate_gas"
gas_thermo_mode = "ideal_gas" #Ideal gas approximation
#gas_thermo_mode = "zero_point_gas" #uses zero-point corrections only
#gas_thermo_mode = "fixed_entropy_gas" #assumes entropy of 0.002 eV/K for all gasses except H2 (H2 is 0.00135 eV/K)
#gas_thermo_mode = "frozen_gas" #neglect thermal contributions
#指明对吸附物种进行何种校正
#adsorbate_thermo_mode = "frozen_adsorbate"
adsorbate_thermo_mode = "harmonic_adsorbate"
#adsorbate_thermo_mode = "hindered_adsorbate"
#adsorbate_thermo_mode = "zero_point_adsorbate"
#反应机制，涉及反应步骤，当需输出Free Energy Diagrams，必须指明
rxn_mechanisms = {
"steps": [1, 2, 3, 4, 5, 6],
}
decimal_precision = 500 #precision of numbers involved
tolerance = 1e-50 #all d_theta/d_t's must be less than this at the solution
max_rootfinding_iterations = 500
max_bisections = 5 
