-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathkmc_writer.py
793 lines (624 loc) · 31.4 KB
/
kmc_writer.py
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
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
# Class definitions for KMC writers
import os
import copy
import yaml
import numpy as np
import shutil
import datetime
import abc
import cantera as ct
import rmgpy.kinetics
import rmgpy.reaction
import scipy.optimize
# import matplotlib
# matplotlib.use('TkAgg')
import matplotlib.pyplot as plt
PT111_SDEN = 1.0e-5 # surface site density of Pt(111) surface in mol/m^2
def plot_kinetics(rxns, labels=None, savedir='./'):
"""Function for plotting reaction kinetics
Takes in a list of RMG reactions (rmgpy.reaction.Reaction) or a single reaction
"""
plt.xlabel('1000 / T (K^-1)')
plt.ylabel('ln(k)')
if type(rxns) != list:
rxns = [rxns]
T = np.linspace(300, 3000, 1001)
for rxn in rxns:
k = np.zeros(len(T))
for i in range(0, len(T)):
k[i] = rxn.get_rate_coefficient(T[i], 101325, surface_site_density=PT111_SDEN)
# if type(rxn.kinetics) == rmgpy.kinetics.surface.StickingCoefficient:
# k[i] = rxn.get_rate_coefficient(T[i], 101325, surface_site_density=PT111_SDEN)
# else:
# k[i] = rxn.get_rate_coefficient(T[i], 101325)
plt.plot(1000.0 / T, np.log(k))
if labels:
plt.legend(labels)
plt.show()
plt.savefig(os.path.join(savedir, 'test.png'))
class ZacrosReaction(rmgpy.reaction.Reaction):
"""Reaction with extra attributes for conversion to Zacros file
This is just a place to hold extra conversion functions"""
def __init__(self, rmg_reaction, SDEN=PT111_SDEN):
super().__init__(
reactants=rmg_reaction.reactants,
products=rmg_reaction.products,
kinetics=rmg_reaction.kinetics,
reversible=rmg_reaction.reversible,
)
# super().__init__(*args, **kwargs)
self.KMC_kinetics = None
# self.reactants = None
# self.products = None
self.species_name_dict = {}
def convert_kinetics(self):
pass
def write_mechanism_step(self, T, SDEN=PT111_SDEN, savedir='./'):
print(f'writing a mechanism step for {self}...')
step_lines = []
translation_method = 'pre_exponential'
# TODO add reversibility
step_name = ''.join(str(self).split())
# Figure out the type of reaction this is:
reaction_type = get_reaction_type(self)
gas_reacs_prods_string = ''
for sp in self.reactants + self.products:
if sp.contains_surface_site():
continue
gas_reacs_prods_string += f' {self.species_name_dict[sp]} {self.get_stoichiometric_coefficient(sp)}'
if self.reversible:
step_lines.append(f'reversible_step {step_name}\n\n')
else:
step_lines.append(f'step {step_name}\n\n')
gas_reacs_prods_string = ''
for sp in self.reactants + self.products:
if sp.contains_surface_site():
continue
gas_reacs_prods_string += f' {self.species_name_dict[sp]} {self.get_stoichiometric_coefficient(sp)}'
step_lines.append(f' gas_reacs_prods {gas_reacs_prods_string}\n')
# TODO read the sites from the adjacency list - this is very far down the road
dentate = 1
print(f'Reaction Type: {reaction_type}')
if reaction_type == 'adsorption':
step_lines.append(' sites 1\n')
step_lines.append(' initial\n')
step_lines.append(f' 1 * {dentate}\n')
step_lines.append(' final\n')
step_lines.append(f' 1 {self.products[0].label} {dentate}\n\n')
step_lines.append(' variant top\n')
step_lines.append(' site_types\ttop\n')
elif reaction_type == 'desorption':
step_lines.append(' sites 1\n')
step_lines.append(' initial\n')
step_lines.append(f' 1 {self.reactants[0].label} {dentate}\n\n')
step_lines.append(' final\n')
step_lines.append(f' 1 * {dentate}\n')
step_lines.append(' variant top\n')
step_lines.append(' site_types\ttop\n')
elif reaction_type == 'dissociative_adsorption':
step_lines.append(' sites 2\n')
step_lines.append(' neighboring 1-2\n')
step_lines.append(' initial\n')
step_lines.append(f' 1 * {dentate}\n')
step_lines.append(f' 2 * {dentate}\n')
step_lines.append(' final\n')
step_lines.append(f' 1 {self.products[0].label} {dentate}\n')
step_lines.append(f' 2 {self.products[1].label} {dentate}\n\n')
step_lines.append(' variant top\n')
step_lines.append(' site_types\ttop top\n')
elif reaction_type == 'associative_desorption':
step_lines.append(' sites 2\n')
step_lines.append(' neighboring 1-2\n')
step_lines.append(' initial\n')
step_lines.append(f' 1 {self.reactants[0].label} {dentate}\n')
step_lines.append(f' 2 {self.reactants[1].label} {dentate}\n')
step_lines.append(' final\n')
step_lines.append(f' 1 * {dentate}\n')
step_lines.append(f' 2 * {dentate}\n\n')
step_lines.append(' variant top\n')
step_lines.append(' site_types\ttop top\n')
else:
raise NotImplementedError(f'Cannot handle reaction type {reaction_type}')
# TODO add more site types
R = 8.314 # J/mol/K
if translation_method == 'pre_exponential':
# TODO assert monodentate
# convert whatever it is to a simple arrhenius form (only pre-exponential factor and activation energy)):
simple_fwd_kinetics = other_kinetics2simple_arrhenius(self, T, SDEN=SDEN)
A = simple_fwd_kinetics.A.value_si # m^3/mol/s
n = simple_fwd_kinetics.n.value_si # dimensionless
Ea = simple_fwd_kinetics.Ea.value_si # J/mol
assert n == 0, 'n should equal zero in this format'
fake_reaction = copy.deepcopy(self)
fake_reaction.kinetics = simple_fwd_kinetics
plot_kinetics([self, fake_reaction], labels=['RMG', 'Zacros'], savedir=savedir)
# I should automatically make a report about the rate conversion process. What the units are...
# Assuming only a pre-exponential factor, no activation energy
if type(self.kinetics) == rmgpy.kinetics.surface.StickingCoefficient or \
type(self.kinetics) == rmgpy.kinetics.surface.SurfaceArrhenius or \
type(self.kinetics) == rmgpy.kinetics.arrhenius.Arrhenius:
Ea = J_to_eV(Ea)
if self.reversible:
PE_ratio = np.exp(self.get_entropy_of_reaction(T) / R)
# divide by RT
# now the reaction rate is in units # m^3/mol/s and we need to convert to 1/bar/s
if reaction_type == 'adsorption':
# if reaction.reversible:
# PE_ratio = np.exp(reaction.get_entropy_of_reaction(T) / R)
A = A / R / T * 1e5
elif reaction_type == 'dissociative_adsorption':
coordination_number = 6.0
A = 2.0 * A / R / T * 1e5 * SDEN / coordination_number
# coordination_number = 6.0
# A = A / R / T * 1e5 * SDEN / coordination_number
# if reaction.reversible:
# PE_ratio = np.exp(reaction.get_entropy_of_reaction(T) / R)
# A = A / R / T * 1e5 * site_density
# print(site_density)
# A = A / R / T / R / T * 1e5 / 2.0
elif reaction_type == 'associative_desorption':
A = A / R / T * 1e5
else:
raise NotImplementedError(f'Cannot handle reaction type {reaction_type}')
PE_units = '1/bar/s'
else:
raise NotImplementedError(f'Cannot handle kinetics type {type(self.kinetics)}')
step_lines.append(f' pre_expon\t{A}\t# {PE_units}\n')
step_lines.append(f' activ_eng\t{Ea} # eV\n')
if self.reversible:
step_lines.append(f' pe_ratio\t{PE_ratio}\n')
step_lines.append(' end_variant\n\n')
if self.reversible:
step_lines.append('end_reversible_step\n\n')
else:
step_lines.append('end_step\n\n\n')
return step_lines
def write_diffusion_step(species):
if species.is_surface_site():
return [] # an empty site cannot diffuse
print(f'writing a mechanism step for {species}...')
pre_expon = 10000000.0 # 1/bar/s
barrier = 0.15 # eV
step_lines = []
gas_reacs_prods_string = ''
step_name = f'diffusion_{species.label}'
# For simplicity, make it irreversible
step_lines.append(f'step {step_name}\n\n')
step_lines.append(f' gas_reacs_prods {gas_reacs_prods_string}\n')
# TODO read the sites from the adjacency list - this is very far down the road
dentate = 1
step_lines.append(' sites 2\n')
step_lines.append(' neighboring 1-2\n')
step_lines.append(' initial\n')
step_lines.append(f' 1 {species.label} {dentate}\n')
step_lines.append(f' 2 * {dentate}\n')
step_lines.append(' final\n')
step_lines.append(f' 1 * {dentate}\n')
step_lines.append(f' 2 {species.label} {dentate}\n\n')
step_lines.append(' variant top\n')
step_lines.append(' site_types\ttop top\n')
step_lines.append(f' pre_expon\t{pre_expon}\n')
step_lines.append(f' activ_eng\t{barrier}\n')
step_lines.append(' end_variant\n\n')
step_lines.append('end_step\n\n\n')
return step_lines
def get_reaction_type(reaction):
"""
function to detect the forward reaction type
"""
if len(reaction.reactants) == 2 and len(reaction.products) == 1 and \
any([sp.is_surface_site() for sp in reaction.reactants]) and \
any([not sp.contains_surface_site() for sp in reaction.reactants]):
return 'adsorption'
elif len(reaction.reactants) == 1 and len(reaction.products) == 2 and \
any([sp.is_surface_site() for sp in reaction.products]) and \
any([not sp.contains_surface_site() for sp in reaction.products]):
return 'desorption'
elif len(reaction.reactants) == 3 and len(reaction.products) == 2 and \
np.sum([sp.is_surface_site() for sp in reaction.reactants]) == 2 and \
np.sum([sp.contains_surface_site() for sp in reaction.products]) == 2 and \
np.sum([not sp.contains_surface_site() for sp in reaction.reactants]) == 1:
return 'dissociative_adsorption'
elif len(reaction.reactants) == 2 and len(reaction.products) == 3 and \
np.sum([sp.is_surface_site() for sp in reaction.products]) == 2 and \
np.sum([sp.contains_surface_site() for sp in reaction.reactants]) == 2 and \
np.sum([not sp.contains_surface_site() for sp in reaction.products]) == 1:
return 'associative_desorption'
raise NotImplementedError
def rename_species(name):
"""function to rename species names to strings that can be used as Python variables"""
name = name.replace('(', '_')
name = name.replace(')', '')
return name
class KMCWriter(abc.ABC):
@abc.abstractmethod
def write(self):
raise NotImplementedError
class MonteCoffeeWriter(KMCWriter):
def write(self, output_dir, species_list, reaction_list):
# copy the template files to the output_dir
template_dir = os.path.join(os.path.dirname(__file__), 'montecoffee', 'template_files')
shutil.copy(os.path.join(template_dir, 'run.py'), output_dir)
shutil.copy(os.path.join(template_dir, 'user_kmc.py'), output_dir)
shutil.copy(os.path.join(template_dir, 'user_system.py'), output_dir)
# Make a dictionary mapping species variable names to species
# need to write the user_sites.py classes
user_sites_path = os.path.join(output_dir, 'user_sites.py')
user_sites_lines = self.generate_user_sites(species_list)
with open(user_sites_path, 'w') as f:
f.writelines(user_sites_lines)
# need to write the user_events.py classes
user_events_path = os.path.join(output_dir, 'user_events.py')
user_events_lines = self.generate_user_events(reaction_list)
with open(user_events_path, 'w') as f:
f.writelines(user_events_lines)
def generate_user_sites(self, species_list):
lines = [f'# User sites autogenerated by rmg2kmc {datetime.datetime.now()}\n']
lines.append('import base.sites\n\n\n')
# for now, this only handles 4 site types for 111 facets
lines.append('# Define constant site types\n')
lines.append('SITE_FCC111_TOP = 0\n')
lines.append('SITE_FCC111_BRIDGE = 1\n')
lines.append('SITE_FCC111_FCC_HOLLOW = 2\n')
lines.append('SITE_FCC111_HCP_HOLLOW = 3\n')
lines.append('\n')
lines.append('# Define species constants\n')
# convert cantera names to python variable names
# TODO check that no species was lost in the conversion
# TODO reserve zero for vacant site
for i, species in enumerate(species_list):
lines.append(f'SPECIES_{rename_species(species.name)} = {i}\n')
lines.append('\n')
# add a basic site
lines.append('class Site(base.sites.SiteBase):\n')
lines.append(' def __init__(\n')
lines.append(' self,\n')
lines.append(' stype=SITE_FCC111_TOP,\n')
lines.append(' covered=0, # covered is the species index\n')
lines.append(' ind=[], # ase indices of site-related atoms\n')
lines.append(' lattice_pos=None\n')
lines.append(' ):\n')
lines.append(' base.sites.SiteBase.__init__(\n')
lines.append(' self,\n')
lines.append(' stype=stype,\n')
lines.append(' covered=covered,\n')
lines.append(' ind=ind,\n')
lines.append(' lattice_pos=lattice_pos\n')
lines.append(' )\n')
return lines
def generate_user_events(self, reaction_list):
# TODO remove dependence of MonteCoffee simulations on Cantera
lines = [f'# User events autogenerated by rmg2kmc {datetime.datetime.now()}\n']
lines.append('import cantera as ct\n')
lines.append('import base.events\n')
lines.append('import user_sites\n\n\n')
# need to convert reaction into a class name
# TODO add Diffusion reactions
for i, reaction in enumerate(reaction_list):
if len(reaction.reactants) > 2:
raise NotImplementedError('Trimolecular reactions not yet supported')
# figure out if it's a surface reaction or adsorption
reaction_type = get_reaction_type(reaction)
# define forward reaction
lines.append(f'class Reaction{i}Fwd(base.events.EventBase):\n')
lines.append(f' # {reaction.equation}\n')
lines.append(' def __init__(self, params):\n')
lines.append(' base.events.EventBase.__init__(self, params, name={reaction.equation})\n\n')
lines.append(' def possible(self, system, site, other_site):\n')
# TODO actually fill this out
lines.append(' return True\n\n')
# Define the reaction rate
lines.append(' def get_rate(self, system, site, other_site):\n')
# TODO add temperature dependence
if not isinstance(reaction.rate, ct.Arrhenius):
raise NotImplementedError('Kinetics only implemented for type Arrhenius')
lines.append(' T = 1000.0\n')
A = reaction.rate.pre_exponential_factor
Ea = reaction.rate.activation_energy
b = reaction.rate.temperature_exponent
lines.append(f' kinetics = ct.Arrhenius(A={A}, b={b}, E={Ea})\n')
lines.append(' return kinetics(T)\n\n')
# TODO fill out actual event actions
lines.append(' def do_event(self, system, site, other_site):\n')
lines.append(' pass\n\n')
lines.append(' def get_involve_other(self):\n')
lines.append(' return False\n\n\n')
if reaction.reversible:
# TODO define reverse reaction
pass
return lines
def get_species_name(species, species_list):
for sp in species_list:
if sp.is_isomorphic(species):
return sp.label
def J_to_eV(joules_per_mol):
"""Converts a value in J/mol to eV"""
N_A = 6.02214076e23
e = 1.602176634e-19
return joules_per_mol / N_A / e
def other_kinetics2simple_arrhenius(rxn, T0, SDEN=2.77e-5):
"""converts any kinetic model into a simple arrhenius form
with just a pre-exponential and activation energy
The kinetics will match the original model exactly at Temperature T0
This is expecting a reaction object, not a kinetics
"""
def my_function(x, m):
return m * x
reaction_type = get_reaction_type(rxn)
# handle case where there's just a pre-exponential factor
if type(rxn.kinetics) != rmgpy.kinetics.surface.StickingCoefficient and \
rxn.kinetics.Ea.value_si == 0 and rxn.kinetics.n.value_si == 0:
print(rxn.kinetics)
return rxn.kinetics
A = rxn.kinetics.A.value_si
simple_arrhenius_kinetics = rmgpy.kinetics.surface.SurfaceArrhenius()
if reaction_type == 'adsorption':
simple_arrhenius_kinetics.A = (A, 'm^3/(mol*s)')
elif reaction_type == 'dissociative_adsorption':
simple_arrhenius_kinetics.A = (A, 'm^5/(mol^2*s)')
simple_arrhenius_kinetics.Ea = (0, 'kJ/mol')
simple_arrhenius_kinetics.n = 0
return simple_arrhenius_kinetics
Ts = np.linspace(300, 3000, 1001)
lnks = np.zeros(len(Ts))
needs_site_density = False
if type(rxn.kinetics) == rmgpy.kinetics.surface.StickingCoefficient:
needs_site_density = True
for i in range(0, len(Ts)):
if needs_site_density:
lnks[i] = np.log(rxn.get_rate_coefficient(Ts[i], 101325, surface_site_density=SDEN))
else:
lnks[i] = np.log(rxn.get_rate_coefficient(Ts[i], 101325))
R = 8.314 # J/mol/K
xs = 1.0 / Ts
x0 = 1.0 / T0
if needs_site_density:
y0 = np.log(rxn.get_rate_coefficient(T0, 101325, surface_site_density=SDEN))
else:
y0 = np.log(rxn.get_rate_coefficient(T0, 101325))
ydata = lnks - y0
xdata = xs - x0
popt, pconv = scipy.optimize.curve_fit(my_function, xdata, ydata)
# transform the variables back
m = popt[0]
b = y0 - m * x0
A = np.exp(b)
Ea = -m * R
simple_arrhenius_kinetics = rmgpy.kinetics.surface.SurfaceArrhenius()
simple_arrhenius_kinetics.Ea = (Ea, 'J/mol')
if reaction_type == 'adsorption':
simple_arrhenius_kinetics.A = (A, 'm^3/(mol*s)')
elif reaction_type == 'dissociative_adsorption':
simple_arrhenius_kinetics.A = (A, 'm^3/(mol*s)')
# simple_arrhenius_kinetics.A = (A, 'm^5/(mol^2*s)')
print(simple_arrhenius_kinetics)
return simple_arrhenius_kinetics
class ZacrosWriter(KMCWriter):
def write_lattice_file(self, output_dir):
lattice_path = os.path.join(output_dir, 'lattice_input.dat')
lines = []
lines.append(f'# Lattice input autogenerated by rmg2kmc {datetime.datetime.now()}\n')
lines.append('lattice periodic_cell\n\n')
lines.append('# R = sqrt(2) a / 4\n')
lines.append('# a = 3.9239A\n')
lines.append('# R = 1.3873 A\n')
lines.append('# (2R, 0)\n')
lines.append('# (R, sqrt(3)R)\n')
lines.append('cell_vectors # in row format (Angstroms)\n')
lines.append(' 2.774616298697894 0.000000000000000\n')
lines.append(' 1.387308149348947 2.402888200426728\n\n')
lines.append('repeat_cell 20 20\n\n')
lines.append('n_site_types 1\n')
lines.append('site_type_names top\n\n')
lines.append('n_cell_sites 1\n')
lines.append('site_types top\n\n')
lines.append('site_coordinates # fractional coordinates (x,y) in row format\n')
lines.append(' 0.000000000000000 0.000000000000000\n\n')
lines.append('neighboring_structure\n')
lines.append(' 1-1 north\n')
lines.append(' 1-1 east\n')
lines.append(' 1-1 southeast\n\n')
lines.append('end_neighboring_structure\n\n')
lines.append('end_lattice\n\n')
with open(lattice_path, 'w') as f:
f.writelines(lines)
def write_energetics_file(self, output_dir, species_list, T):
energetics_path = os.path.join(output_dir, 'energetics_input.dat')
lines = []
lines.append(f'# Energetics input autogenerated by rmg2kmc {datetime.datetime.now()}\n')
lines.append('energetics\n\n')
# # Get the energy of an empty site
# for species in species_list:
# if species.is_surface_site():
# h_empty = J_to_eV(species.get_enthalpy(T)) # TODO get the enthalpy of the surface site
# break
# else:
# raise ValueError('No empty surface site found in species list')
# # TODO - perhaps limp on with h_empty = 0?
# Get the energy of the gas version
for species in species_list:
if not species.contains_surface_site():
continue
if species.is_surface_site():
continue
lines.append(f'cluster {species.label}_top\n') # TODO other sites?
lines.append(' sites 1\n')
lines.append(' lattice_state\n')
lines.append(f' 1 {species.label} 1\n')
lines.append(' site_types top\n')
lines.append(f' cluster_eng {np.round(J_to_eV(species.get_enthalpy(T)), 5)} # eV\n\n')
lines.append(f'end_cluster {species.label}_top\n') # TODO other sites?
lines.append('end_energetics\n\n')
with open(energetics_path, 'w') as f:
f.writelines(lines)
def write_mechanism_file(self, output_dir, species_list, reaction_list, T, site_density):
"""
Write the Zacros mechanism file
translate the rate using either a simple pre-exponential or Zacros's 7-param equation
"""
translation_method = 'pre_exponential' # TODO add option for 7-param
# translation_method = '7_param'
# WRITE THE MECHANISM FILE
mechanism_path = os.path.join(output_dir, 'mechanism_input.dat')
lines = []
lines.append(f'# Mechanism file autogenerated by rmg2kmc {datetime.datetime.now()}\n')
lines.append(f'# Temperature: {T}K\n')
lines.append('\n')
lines.append('\n')
lines.append('mechanism\n')
lines.append('\n')
reaction_list = [ZacrosReaction(rxn) for rxn in reaction_list]
diffusion = True
# write each reaction as an irreversible step
for reaction in reaction_list:
for sp in reaction.reactants + reaction.products:
sp_name = get_species_name(sp, species_list)
reaction.species_name_dict[sp] = sp_name
if diffusion and sp.contains_surface_site():
lines += write_diffusion_step(sp)
lines += reaction.write_mechanism_step(T, SDEN=site_density, savedir=output_dir)
lines.append('end_mechanism\n')
with open(mechanism_path, 'w') as f:
f.writelines(lines)
def write_simulation_file(self, output_dir, species_list, T, P, starting_gas_conc, simulation_file=None):
# WRITE THE GENERAL SIMULATION INPUT FILE
simulation_path = os.path.join(output_dir, 'simulation_input.dat')
# TODO get timeframe from Cantera simulation?
# get the simulation settings
if simulation_file:
with open(simulation_file, 'r') as f:
simulation_settings = yaml.safe_load(f)
t_end = simulation_settings['t_end']
try:
random_seed = simulation_settings['random_seed']
except KeyError:
random_seed = 8949321
else:
t_end = 1e-6
random_seed = 8949321
lines = []
lines.append(f'# Simulation input file autogenerated by rmg2kmc {datetime.datetime.now()}\n')
lines.append(f'random_seed\t\t{random_seed}\n\n')
lines.append(f'temperature\t\t{float(T)} # K\n')
lines.append(f'pressure\t\t{P} # bar\n\n')
# get num gas species
num_gas_species = 0
num_surf_species = 0
gas_species = []
surf_species = []
surface_site_sp_name = ''
for species in species_list:
if species.is_surface_site():
surface_site_sp_name = str(species) # what we expect to see in cantera settings
continue
if species.contains_surface_site():
surf_species.append(species)
num_surf_species += 1
else:
gas_species.append(species)
num_gas_species += 1
lines.append(f'n_gas_species\t\t{num_gas_species}\n')
lines.append(f'gas_specs_names\t\t{" ".join([str(sp.label) for sp in gas_species])}\n')
# get the gas energies using the RMG thermo # convert from J/mol to eV per molecule
lines.append(f'gas_energies\t\t{" ".join([str(np.round(J_to_eV(sp.get_enthalpy(T)), 5)) for sp in gas_species])} # eV\n')
lines.append(f'gas_molec_weights\t\t{" ".join([str(np.round(sp.molecular_weight.value, 5)) for sp in gas_species])}\n')
# TODO add something that makes sense
lines.append(f'gas_molar_fracs\t\t{starting_gas_conc}\n\n')
# lines.append(f'gas_molar_fracs\t\t{" ".join([str(np.round(1.0 / num_gas_species, 2)) for sp in gas_species])}\n\n')
lines.append(f'n_surf_species\t\t{num_surf_species}\n')
lines.append(f'surf_specs_names\t\t{" ".join([str(sp.label) for sp in surf_species])}\n')
# TODO add bidentate species
lines.append(f'surf_specs_dent\t\t{" ".join(["1" for sp in surf_species])}\n\n')
lines.append(f'snapshots\t\ton time {t_end / 100.0}\n')
lines.append(f'process_statistics\t\ton time {t_end / 100.0}\n')
lines.append(f'species_numbers\t\ton time {t_end / 100.0}\n\n')
lines.append(f'# event_report\t\ton\n\n')
lines.append(f'max_steps\t\tinfinity # steps\n')
lines.append(f'max_time\t\t{t_end} # seconds\n\n')
lines.append(f'wall_time\t\t3000 # seconds\n\n')
lines.append(f'no_restart\t\t\n\n')
lines.append(f'# debug_report_processes\n')
lines.append(f'# debug_report_global_energetics\n')
lines.append(f'# debug_check_processes\n')
lines.append(f'# debug_check_lattice\n\n')
lines.append(f'finish\n\n')
with open(simulation_path, 'w') as f:
f.writelines(lines)
def write_initial_state_file(self, output_dir, species_list, simulation_file=None):
initial_state_file = os.path.join(output_dir, 'state_input.dat')
if os.path.exists(initial_state_file):
print('Cleaning up old initial state file...')
os.remove(initial_state_file)
if not simulation_file:
return
with open(simulation_file, 'r') as f:
simulation_settings = yaml.safe_load(f)
try:
random_seed = simulation_settings['random_seed']
except KeyError:
random_seed = 8949321
try:
initial_surface_coverage_setting = simulation_settings['surface_coverages_ct']
except KeyError:
# no initial surface coverage setting
return
assert type(initial_surface_coverage_setting) == dict, 'surface_coverages_ct must be a dictionary'
# check the simulation file surface_coverages_ct to see if it's anything other than an empty surface
surf_species = []
surface_site_sp_name = ''
for species in species_list:
if species.is_surface_site():
surface_site_sp_name = str(species) # what we expect to see in cantera settings
surf_species.append(species)
elif species.contains_surface_site():
surf_species.append(species)
if not surface_site_sp_name:
raise ValueError('No empty surface site found in species list')
if len(initial_surface_coverage_setting) == 0:
raise ValueError('Initial surface coverage setting must contain at least one species')
elif len(initial_surface_coverage_setting) == 1:
if list(initial_surface_coverage_setting.keys())[0] == surface_site_sp_name and \
initial_surface_coverage_setting[surface_site_sp_name] == 1.0:
return # because it's the default state
print('Creating initial state file')
total_coverage_sum = 0
surf_spec_names = [str(sp) for sp in surf_species]
# surf_spec_names = [str(sp.label) for sp in surf_species]
num_spec_to_place = {}
N_sites = 20 * 20 # TODO add flexibility in lattice size
for sp_key in initial_surface_coverage_setting.keys():
total_coverage_sum += initial_surface_coverage_setting[sp_key]
# also check that the species name is valid and get the real species name??
if sp_key not in surf_spec_names:
print(surf_spec_names)
raise ValueError(f'Invalid species name {sp_key} in initial surface coverage setting')
if sp_key == surface_site_sp_name:
continue
# get the corresponding species
sp = [sp for sp in surf_species if str(sp) == sp_key][0]
# use the species label for the actual key
num_spec_to_place[str(sp.label)] = int(initial_surface_coverage_setting[sp_key] * N_sites)
if total_coverage_sum != 1.0:
raise ValueError('Initial surface coverage setting must sum to 1.0')
# TODO modify # sites by dentate
with open(initial_state_file, 'w') as f:
f.write(f'# Initial state file autogenerated by rmg2kmc {datetime.datetime.now()}\n\n')
f.write('initial_state\n\n')
for sp_key in num_spec_to_place.keys(): # but not species list, just the
f.write(f' seed_multiple {sp_key} {num_spec_to_place[sp_key]}' + '\n')
f.write(' site_types top\n')
f.write(' end_seed_multiple\n\n')
f.write('end_initial_state\n\n')
def write(self, output_dir, species_list, reaction_list, T, P, starting_gas_conc, site_density=2.72e-5, simulation_file=None):
# make the mechanism file
if not os.path.exists(output_dir):
os.makedirs(output_dir, exist_ok=True)
# WRITE THE SIMULATION FILE
self.write_simulation_file(output_dir, species_list, T, P, starting_gas_conc, simulation_file=simulation_file)
self.write_initial_state_file(output_dir, species_list, simulation_file=simulation_file)
# WRITE THE MECHANISM FILE
self.write_mechanism_file(output_dir, species_list, reaction_list, T, site_density)
# WRITE THE LATTICE FILE
self.write_lattice_file(output_dir)
# WRITE THE ENERGETICS FILE
self.write_energetics_file(output_dir, species_list, T)