|
|
| Line 1: |
Line 1: |
| ==CODE==
| | Generate Glycerol Decomposition Network and Complete Intermediates |
| | |
| #!/usr/bin/env python
| |
| import operator
| |
| import csv
| |
| import subprocess
| |
| | |
| dict_l = {'1312' : 'HOCH2', '1300' : 'CH3', '1212' : 'HOCH', '1211' : 'OCH2',
| |
| '1200' : 'CH2', '1112' : 'HOC', '1111' : 'OCH', '1100' : 'CH',
| |
| '1011' : 'OC', '1000' : 'C', '0000' : '0000'}
| |
| dict_r = {'1312' : 'CH2OH', '1300' : 'CH3', '1212' : 'CHOH', '1211' : 'CH2O',
| |
| '1200' : 'CH2', '1112' : 'COH', '1111' : 'CHO', '1100' : 'CH',
| |
| '1011' : 'CO', '1000' : 'C', '0000' : '0000'}
| |
| | |
| # Function :Get Formulas
| |
| | |
| def Get_formula1(label): # get formula such as HOCH2-CHOH-CH2OH from 1312-1212-1312
| |
| if ':' in label:
| |
| return '%s-%s:%s' %(dict_l.get(label[0:4]), dict_r.get(label[5:9]),dict_r.get(label[10:]) )
| |
| else:
| |
| return '%s-%s-%s' %(dict_l.get(label[0:4]), dict_r.get(label[5:9]),dict_r.get(label[10:]) )
| |
| | |
| def Get_formula2(label): # Get formula such as C3H8O3 from 1312-1212-1312
| |
| if ':' not in label: # ':' is used to define the broken C--C bonds, for example, 1312-1212:1312 means HOCH-CHOH------CH2OH bond breaks
| |
| C_number = int(label[0]) + int(label[5]) + int(label[10])
| |
| H_number = int(label[1]) + int(label[6]) + int(label[11])
| |
| O_number = int(label[2]) + int(label[7]) + int(label[12])
| |
| return 'C%sH%sO%s' %(C_number, H_number, O_number)
| |
| else:
| |
| C_number = int(label[0]) + int(label[5])
| |
| H_number = int(label[1]) + int(label[6])
| |
| O_number = int(label[2]) + int(label[7])
| |
| return 'C%sH%sO%s' %(C_number, H_number, O_number)
| |
| #print Get_formula('1312-1212-1312') # test function
| |
| | |
| def Get_formula3(label): # this function is used to get smaller species (C1) in C--C bond breaking products:
| |
| return '%s' %(dict_r.get(label[10:]))
| |
| | |
| # Keep label sequence unique for C3 Species
| |
| def label_reverse(molecular):
| |
| if float(molecular[10:]) < float(molecular[0:4]):
| |
| molecular = molecular[10:] + molecular[4:10] + molecular[0:4]
| |
| return molecular
| |
| else:
| |
| return molecular
| |
| | |
| # Keep label sequence unique for C2 Species
| |
| def label_reverse_C2(molecular):
| |
| if float(molecular[5:]) < float(molecular[0:4]):
| |
| molecular = molecular[5:] + '-' + molecular[0:4]
| |
| return molecular
| |
| else:
| |
| return molecular
| |
| | |
| # Defination of O--H bond breaking
| |
| def OH_breaking(Start_molecular):
| |
| product_list = []
| |
| if Start_molecular[3] == '2' : # There is a OH here!
| |
| product = label_reverse(str(int(Start_molecular[0:4]) - 101) + Start_molecular[4:])
| |
| product_list.append(product)
| |
| if Start_molecular[8] == '2' :
| |
| product = label_reverse(Start_molecular[0:5] + str(int(Start_molecular[5:9]) - 101) + Start_molecular[9:])
| |
| product_list.append(product)
| |
| if Start_molecular[13] == '2' :
| |
| product = label_reverse(Start_molecular[0:10] + str(int(Start_molecular[10:]) - 101) )
| |
| if product not in product_list:
| |
| product_list.append(product)
| |
| return product_list
| |
| #print OH_breaking('1312-1212-1312')
| |
| | |
| # Defination of C--H bond breaking
| |
| def CH_breaking(Start_molecular):
| |
| product_list = []
| |
| if float(Start_molecular[1]) > 1 : # There is a H attaching to C
| |
| product = label_reverse(str(int(Start_molecular[0:4]) - 100) + Start_molecular[4:])
| |
| product_list.append(product)
| |
| if float(Start_molecular[1]) == 1 : # To judge CH or OH
| |
| if float(Start_molecular[3]) <= 1 : # If OH; Start_molecular[3] == 2
| |
| product = label_reverse(str(int(Start_molecular[0:4]) - 100) + Start_molecular[4:])
| |
| product_list.append(product)
| |
|
| |
| if float(Start_molecular[6]) > 1 :
| |
| product = label_reverse(Start_molecular[0:5] + str(int(Start_molecular[5:9]) - 100) + Start_molecular[9:])
| |
| product_list.append(product)
| |
| if float(Start_molecular[6]) == 1 :
| |
| if float(Start_molecular[8]) <= 1 :
| |
| product = label_reverse(Start_molecular[0:5] + str(int(Start_molecular[5:9]) - 100) + Start_molecular[9:])
| |
| product_list.append(product)
| |
| if float(Start_molecular[11]) > 1 :
| |
| product = label_reverse(Start_molecular[0:10] + str(int(Start_molecular[10:]) - 100) )
| |
| if product not in product_list:
| |
| product_list.append(product)
| |
| if float(Start_molecular[11]) == 1 :
| |
| if float(Start_molecular[13]) <= 1 :
| |
| product = label_reverse(Start_molecular[0:10] + str(int(Start_molecular[10:]) - 100) )
| |
| if product not in product_list:
| |
| product_list.append(product)
| |
| return product_list
| |
| #print CH_breaking('1312-1212-1312')
| |
| | |
| # Defination of C--O bond breaking, C--OH is not included
| |
| def CO_breaking(Start_molecular):
| |
| product_list = []
| |
| if Start_molecular[3] == '1' : # Start_molecular[3] == 1 means only O terminal
| |
| product = label_reverse(str(int(Start_molecular[0:4]) - 11) + Start_molecular[4:])
| |
| product_list.append(product)
| |
| if Start_molecular[8] == '1' :
| |
| product = label_reverse(Start_molecular[0:5] + str(int(Start_molecular[5:9]) - 11) + Start_molecular[9:])
| |
| product_list.append(product)
| |
| if Start_molecular[13] == '1' :
| |
| product = label_reverse(Start_molecular[0:10] + str(int(Start_molecular[10:]) - 11) )
| |
| if product not in product_list:
| |
| product_list.append(product)
| |
| return product_list
| |
| #print CO_breaking('1211-1211-1211')
| |
| | |
| # Defination of C--OH bond breaking
| |
| def COH_breaking(Start_molecular):
| |
| product_list = []
| |
| if Start_molecular[3] == '2' : # Start_molecular[3] == 2 means there exist one OH
| |
| product = label_reverse(str(int(Start_molecular[0:4])-112) + Start_molecular[4:])
| |
| product_list.append(product)
| |
| if Start_molecular[8] == '2' :
| |
| product = label_reverse(Start_molecular[0:5] + str(int(Start_molecular[5:9]) - 112) + Start_molecular[9:])
| |
| product_list.append(product)
| |
| if Start_molecular[13] == '2' :
| |
| product = label_reverse(Start_molecular[0:10] + str(int(Start_molecular[10:]) - 112) )
| |
| if product not in product_list:
| |
| product_list.append(product)
| |
| return product_list
| |
| #print COH_breaking('1312-1212-1312')
| |
| | |
| # Defination of C--C bond breaking
| |
| def CC_breaking(Start_molecular):
| |
| product_list = []
| |
| if Start_molecular[10:] != '0000' :
| |
| product1 = label_reverse_C2(Start_molecular[0:9]) + ':' + Start_molecular[10:]
| |
| product2 = label_reverse_C2(Start_molecular[5:]) + ':' + Start_molecular[0:4]
| |
| product_list.append(product1)
| |
| if product2 != product1:
| |
| product_list.append(product2)
| |
| return product_list
| |
| #print CC_breaking('1312-1212-1312')
| |
| | |
| Start_molecular = '1312-1212-1112'
| |
| | |
| def one_step(Start):
| |
| initial_species = []
| |
| def append_species(products):
| |
| for i in products:
| |
| initial_species.append(i)
| |
| append_species(OH_breaking(Start))
| |
| append_species(CH_breaking(Start))
| |
| append_species(COH_breaking(Start))
| |
| append_species(CO_breaking(Start))
| |
| # append_species(CC_breaking(Start)) # Be cautious!!!!!!!!!!!!!!
| |
| return initial_species
| |
|
| |
| list_product = []
| |
| list_product.append('1312-1212-1312')
| |
| | |
| #print '1312-1212-1312', one_step('1312-1212-1312')
| |
| for i in one_step('1312-1212-1312'):
| |
| list_product.append(i)
| |
| # print i , one_step(i)
| |
| for j in one_step(i):
| |
| if j not in list_product:
| |
| list_product.append(j)
| |
| for k in one_step(j):
| |
| if k not in list_product:
| |
| list_product.append(k)
| |
| for g in one_step(k):
| |
| if g not in list_product:
| |
| list_product.append(g)
| |
| for x in one_step(g):
| |
| if x not in list_product:
| |
| list_product.append(x)
| |
| for y in one_step(x):
| |
| if y not in list_product:
| |
| list_product.append(y)
| |
| for z in one_step(y):
| |
| if z not in list_product:
| |
| list_product.append(z)
| |
| for a in one_step(z):
| |
| if a not in list_product:
| |
| list_product.append(z)
| |
| for b in one_step(a):
| |
| if b not in list_product:
| |
| list_product.append(b)
| |
| for c in one_step(b):
| |
| if c not in list_product:
| |
| list_product.append(c)
| |
| # write species' list
| |
| subprocess.call('rm reaction_list.txt species_list.txt -fr', shell = True) # remove old files
| |
| | |
| final_list_product = sorted(set(list_product))
| |
| | |
| sl = open ('species_list.txt', 'wr')
| |
| for i in final_list_product:
| |
| sl.write('%s,%s,%s,%s,%s\n' %(i, Get_formula1(i), Get_formula2(i), Get_formula2(i)[3], Get_formula2(i)[5] ))
| |
| sl.close()
| |
| | |
| # write reaction list
| |
| rl = open('reaction_list.txt', 'wr')
| |
| # Print all reactions but C--C bond breakings
| |
| for i in final_list_product:
| |
| for j in OH_breaking(i):
| |
| rl.write( '%s,%s,%s,%s,%s,%s,%s,%s,%s\n' %(i, Get_formula1(i), j , Get_formula1(j),'0101', 'H', 'O-H', Get_formula2(i), Get_formula2(j)))
| |
| for k in CH_breaking(i):
| |
| rl.write( '%s,%s,%s,%s,%s,%s,%s,%s,%s\n' %(i, Get_formula1(i), k, Get_formula1(k), '0100','H', 'C-H', Get_formula2(i), Get_formula2(k)))
| |
| for l in CO_breaking(i):
| |
| rl.write( '%s,%s,%s,%s,%s,%s,%s,%s,%s\n' %(i, Get_formula1(i), l, Get_formula1(l), '0011', 'O', 'C-O', Get_formula2(i), Get_formula2(l)))
| |
| for m in COH_breaking(i):
| |
| rl.write( '%s,%s,%s,%s,%s,%s,%s,%s,%s\n' %(i, Get_formula1(i), m, Get_formula1(m), '0112','OH', 'C-OH', Get_formula2(i), Get_formula2(m)))
| |
| for x in CC_breaking(i):
| |
| rl.write( '%s,%s,%s,%s,%s,%s,%s,%s,%s\n' %(i, Get_formula1(i), x , Get_formula1(x), x[10:], Get_formula3(x), 'C-C', Get_formula2(i), Get_formula2(x) ))
| |
| | |
| rl.close()
| |
| # From results above, species_list.txt file can be generated, then the species are given the name such as 731, 732,..(in label_dic.txt)
| |
| # Performing DFT calculations about these structures and get their energies. Add energies to label_dic.txt and get c1c2c3.txt file,
| |
| # label_dic.txt are only used during Calculation.
| |
| # c1c2c3.txt contains C1, C2, C3 species' energies
| |
|
| |
| ####################################################
| |
| #Step2: Get reference energies relative to Glycerol
| |
| ####################################################
| |
| | |
| E_slab = -317.09787996 # Ru(0001) slab
| |
| E_H_ads = -320.9686537675 # H atom on Ru(0001) slab
| |
| E_OH_ads = -327.92243824 # OH on Ru(0001) slab
| |
| E_G_gas = -73.469085999 # Glycerol gas phase energy
| |
| | |
| def dicts_1():
| |
| dict9 = {}
| |
| with open ('c1c2c3.txt', mode = 'r') as data1:
| |
| reader = csv.reader(data1, delimiter = ',')
| |
| for rows in reader:
| |
| dict9.update({rows[0] : rows[2]}) # E_DFT_ZPVE_corrected
| |
| return dict9
| |
| | |
| dict9 = dicts_1()
| |
| | |
| def get_reference_energy(label):
| |
| # get absolute energies for Energy Profile Plotting
| |
| # E_OH, E_H, E_slab, and E_G_gas are used as reference energies for all the species.
| |
| if ':' not in label:
| |
| N_H = float(Get_formula2(label)[3])
| |
| N_O = float(Get_formula2(label)[5])
| |
| E_ref = (3 - N_O ) * E_OH_ads + (5 + N_O - N_H) * E_H_ads - (9 - N_H) * E_slab - E_G_gas
| |
| E_label = float(dict9.get(label)) + E_ref
| |
| else :
| |
| N_H = float(Get_formula2(label.replace(':','-'))[3]) # Change Products to Reactants and use its energies
| |
| N_O = float(Get_formula2(label.replace(':','-'))[5])
| |
| E_ref = (3 - N_O ) * E_OH_ads + (5 + N_O - N_H) * E_H_ads - (9 - N_H) * E_slab - E_G_gas
| |
| E_label = float(dict9.get(label[0:9])) + float(dict9.get(label[10:])) - E_slab + E_ref
| |
| return E_label
| |
| | |
| def reference_energies(IS, FS1, FS2):
| |
| N_H = float(Get_formula2(IS)[3])
| |
| N_O = float(Get_formula2(IS)[5])
| |
| E_ref = (3 - N_O ) * E_OH_ads + (5 + N_O - N_H) * E_H_ads - (9 - N_H) * E_slab - E_G_gas
| |
| E_refered_IS = float(dict9.get(IS)) + E_ref
| |
| #Get E_dft_FS
| |
| if ':' in FS1:
| |
| E_refered_FS = float(dict9.get(FS1[0:9])) + float(dict9.get(FS1[10:])) - E_slab + E_ref
| |
| else:
| |
| E_refered_FS = float(dict9.get(FS1)) + float(dict9.get(FS2)) - E_slab + E_ref
| |
|
| |
| Delta_E = E_refered_FS - E_refered_IS
| |
| | |
| return E_refered_IS, E_refered_FS, Delta_E
| |
| | |
| #1312-1212-1312,HOCH2-CHOH-CH2OH,1212-1312:1312,HOCH-CH2OH:CH2OH,CH2OH,C-C,C3H8O3,C2H5O2
| |
| | |
| #Write referenced energies and transition states' energies with scaling equations from Ethylene glycol
| |
| | |
| re = open('reaction_energy.txt', 'wr')
| |
| with open('reaction_list.txt', mode = 'r') as data2:
| |
| reader = csv.reader(data2, delimiter = ',' )
| |
| for rows in reader:
| |
| E1, E2, E3 = reference_energies(rows[0], rows[2], rows[4]) # E1, E2, E3: E_IS, E_FS, Delta_E
| |
| g0 = get_reference_energy(rows[0])
| |
| g2 = get_reference_energy(rows[2])
| |
| r0 = rows[0]
| |
| r2 = rows[2]
| |
| r4 = rows[4]
| |
| r6 = rows[6]
| |
| if 'O-H' in rows[6]:
| |
| re.write('%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s\n'
| |
| %(r0, r2, r4, r6, E1, E2, 0.857*E1+0.251, 0.857*E1+0.251 - E1, E3, g0, g2))
| |
| elif 'C-C' in rows[6]:
| |
| re.write('%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s\n'
| |
| %(r0, r2, r4, r6, E1, E2, 0.955*E2+1.099, 0.955*E2+1.099 - E1, E3, g0, g2))
| |
| elif 'C-OH' in rows[6]:
| |
| re.write('%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s\n'
| |
| %(r0, r2, r4, r6, E1, E2, 1.062*E2+1.777, 1.062*E2+1.777 - E1, E3, g0, g2))
| |
| elif 'C-H' in rows[6]:
| |
| re.write('%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s\n'
| |
| %(r0, r2, r4, r6, E1, E2, 0.903*E2+0.592, 0.903*E2+0.592 - E1, E3, g0, g2))
| |
| else: # C--O breaking, not C--OH !!!
| |
| re.write('%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s\n'
| |
| %(r0, r2, r4, r6, E1, E2, 0.732*E2+1.227, 0.732*E2+1.227 - E1, E3, g0, g2))
| |
| re.close()
| |
| | |
| subprocess.call("""sed -i "s/\[//g; s/\]//g; s/\'//g; s/ //g" reaction_energy.txt""", shell = True) # remove '[' and ']' in the output file
| |
| | |
| ## reaction_list.txt: this file containes the reactions, barriers, and reaction energies
| |
| #
| |
| #
| |
| ########################################################################################
| |
| ## Step 3 Select and filter Reactions
| |
| ########################################################################################
| |
| ## Select Reactions
| |
| #####################################
| |
| #
| |
| def get_line_number(phrase): # Not used!
| |
| numbers = []
| |
| with open('reaction_energy.txt') as f:
| |
| for i, line in enumerate(f, 0):
| |
| if phrase in line:
| |
| numbers.append(i)
| |
| return numbers
| |
| #print get_line_number('1312-1212-1312')
| |
| | |
| def get_reactions(label):
| |
| lis = []
| |
| Ea = []
| |
| lis_s = []
| |
| with open('reaction_energy.txt', mode = 'r') as vertical:
| |
| reader = csv.reader(vertical, delimiter=',')
| |
| for rows in reader:
| |
| if label in rows[0]:
| |
| lis.append(rows[:])
| |
| # Filter_Step1 : Remove reactions that can not proceed
| |
| for i in lis:
| |
| if i[0][5:9] == '1212' and i[1][5:9] == '1112' : # XXXX-CHOH-XXXX --> XXXX-COH-XXXX
| |
| lis.pop(lis.index(i)) #remove these kinds of reactions
| |
| elif i[0][5:9] == '1111' and i[1][5:9] == '1011' : # XXXX-CHO-XXXX --> XXXX-CO-XXXX
| |
| lis.pop(lis.index(i))
| |
| # Filter_Step2: Remove reactions with high barriers
| |
| for i in lis:
| |
| Ea.append(float(i[7]))
| |
|
| |
| lowest = min(i for i in Ea)
| |
| | |
| for i in lis:
| |
| if float(i[7]) - lowest < 0.300:
| |
| lis_s.append(i)
| |
| return lis_s
| |
| | |
| ## rs: selected reactions
| |
| rs = open ('selected_reactions.txt', 'w')
| |
| f = open ('species_list.txt')
| |
| for line in f:
| |
| for i in get_reactions(line.rstrip().split(',')[0]):
| |
| rs.write('%s\n' %(i))
| |
| rs.close()
| |
| f.close()
| |
| | |
| subprocess.call("""sed -i "s/\[//g; s/\]//g; s/\'//g; s/ //g" selected_reactions.txt""", shell = True) # remove '[' and ']' in the output file
| |
| | |
| ###################################################################################################
| |
| ###Filter Reactions
| |
| ###################################################################################################
| |
| | |
| def get_reactions_2(label):
| |
| lis = []
| |
| with open('selected_reactions.txt', mode = 'r') as vertical:
| |
| reader = csv.reader(vertical, delimiter=',')
| |
| for rows in reader:
| |
| if label in rows[0]:
| |
| lis.append(rows[:])
| |
| return lis
| |
| | |
| def cycle(IS):
| |
| reactions = []
| |
| FS = []
| |
| for ele in get_reactions_2(IS):
| |
| if ele[1] not in FS:
| |
| FS.append(ele[1])
| |
| if ele not in reactions:
| |
| reactions.append(ele)
| |
| return reactions, FS
| |
| | |
| ########Get all Intermediate Species involved in the reaction system ############
| |
| f_list = []
| |
| f_list.append('1312-1212-1312')
| |
| def list_append(i):
| |
| if i not in f_list:
| |
| f_list.append(i)
| |
| | |
| for f1 in cycle('1312-1212-1312')[1]:
| |
| list_append(f1)
| |
| for f2 in cycle(f1)[1]:
| |
| list_append(f2)
| |
| for f3 in cycle(f2)[1]:
| |
| list_append(f3)
| |
| for f4 in cycle(f3)[1]:
| |
| list_append(f4)
| |
| for f5 in cycle(f4)[1]:
| |
| list_append(f5)
| |
| for f6 in cycle(f5)[1]:
| |
| list_append(f6)
| |
| for f7 in cycle(f6)[1]:
| |
| list_append(f7)
| |
| for f8 in cycle(f7)[1]:
| |
| list_append(f8)
| |
| for f9 in cycle(f8)[1]:
| |
| list_append(f9)
| |
| fs = open('FS_list.txt', 'w')
| |
| for i in f_list:
| |
| fs.write(i + '\n')
| |
| fs.close()
| |
| #########################################################
| |
| ### Filter Reactions based on previous FS
| |
| #########################################################
| |
| fr = open('filtered_reactions.txt', 'w')
| |
| fs = open('FS_list.txt')
| |
| for line in fs:
| |
| IS = line.rstrip()
| |
| if ':' not in IS:
| |
| for i in get_reactions(IS):
| |
| fr.write('%s\n' %(i))
| |
| fr.close()
| |
| subprocess.call("""sed -i "s/\[//g; s/\]//g; s/\'//g; s/ //g" filtered_reactions.txt""", shell = True) # remove '[' and ']' in the output file
| |
| #############################################################
| |