Version-6.py: Difference between revisions
Jump to navigation
Jump to search
New page: ==CODE== #!/usr/bin/env python import operator import csv import subprocess dict_l = {'1312' : 'HOCH2', '1300' : 'CH3', '1212' : 'HOCH', '1211' : 'OCH2', '1200' : 'CH2', '1112... |
No edit summary |
||
| Line 1: | Line 1: | ||
==CODE== | ==CODE== | ||
#!/usr/bin/env python | #!/usr/bin/env python | ||
import operator | import operator | ||
Revision as of 13:10, 10 February 2017
CODE
- !/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/\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/\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/\s/\//g; s/\'//g; s/ //g" filtered_reactions.txt""", shell = True) # remove '[' and ']' in the output file