| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // Geant4 class G4ChemDissociationChannels_option1
28 //
29 // H. Tran 16.12.2022
30 //
31
32 #include "G4ChemDissociationChannels_option1.hh"
33
34 #include "G4DNAWaterDissociationDisplacer.hh"
35 #include "G4DNAWaterExcitationStructure.hh"
36 #include "G4Electron_aq.hh"
37 #include "G4FakeMolecule.hh"
38 #include "G4H2.hh"
39 #include "G4H2O.hh"
40 #include "G4H2O2.hh"
41 #include "G4H3O.hh"
42 #include "G4HO2.hh"
43 #include "G4Hydrogen.hh"
44 #include "G4MolecularConfiguration.hh"
45 #include "G4MoleculeTable.hh"
46 #include "G4O2.hh"
47 #include "G4O3.hh"
48 #include "G4OH.hh"
49 #include "G4Oxygen.hh"
50 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4Scheduler.hh"
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54 void G4ChemDissociationChannels_option1::ConstructMolecule()
55 {
56 //-----------------------------------
57 // G4Electron::Definition(); // safety
58
59 //-----------------------------------
60 // Create the definition
61
62 G4H2O::Definition();
63 G4Hydrogen::Definition();
64 G4H3O::Definition();
65 G4OH::Definition();
66 G4Electron_aq::Definition();
67 G4H2O2::Definition();
68 G4H2::Definition();
69
70 G4O2::Definition();
71 G4HO2::Definition();
72 G4Oxygen::Definition();
73 G4O3::Definition();
74
75 auto G4OHm = new G4MoleculeDefinition("OH",/*mass*/ 17.00734 * g / Avogadro * c_squared,
76 2.8e-9 * (m * m / s), -1,
77 5, 0.958 * angstrom, // radius
78 2 // number of atoms
79 );
80
81 auto G4HO2m = new G4MoleculeDefinition("HO_2", 33.0034 * g / Avogadro * c_squared,
82 2.3e-9 * (m * m / s), -1, 0,
83 2.1 * angstrom, 3);
84 auto G4Om = new G4MoleculeDefinition("O", 15.99773 * g / Avogadro * c_squared,
85 2.0e-9 * (m * m / s), 0, 0,
86 2.0 * angstrom, 1);
87 //____________________________________________________________________________
88 auto molTable = G4MoleculeTable::Instance();
89 molTable->CreateConfiguration("H3Op", G4H3O::Definition());
90 molTable->GetConfiguration("H3Op")->SetDiffusionCoefficient(9.46e-9
91 * (m2 / s));
92 molTable->GetConfiguration("H3Op")->SetVanDerVaalsRadius(0.25 * nm);
93
94 molTable->CreateConfiguration("°OH", G4OH::Definition());
95 molTable->GetConfiguration("°OH")->SetDiffusionCoefficient(2.2e-9 * (m2 / s));
96 molTable->GetConfiguration("°OH")->SetVanDerVaalsRadius(0.22 * nm);
97
98 G4MolecularConfiguration* OHm =
99 molTable->CreateConfiguration("OHm", // just a tag to store and retrieve
100 // from G4MoleculeTable
101 G4OHm,
102 -1, // charge
103 5.3e-9 * (m2 / s));
104 OHm->SetMass(17.0079 * g / Avogadro * c_squared);
105 OHm->SetVanDerVaalsRadius(0.33 * nm);
106
107 molTable->CreateConfiguration("e_aq", G4Electron_aq::Definition());
108 molTable->GetConfiguration("e_aq")->SetVanDerVaalsRadius(0.50 * nm);
109
110 molTable->CreateConfiguration("H", G4Hydrogen::Definition());
111 molTable->GetConfiguration("H")->SetVanDerVaalsRadius(0.19 * nm);
112
113 molTable->CreateConfiguration("H2", G4H2::Definition());
114 molTable->GetConfiguration("H2")->SetDiffusionCoefficient(4.8e-9 * (m2 / s));
115 molTable->GetConfiguration("H2")->SetVanDerVaalsRadius(0.14 * nm);
116
117 molTable->CreateConfiguration("H2O2", G4H2O2::Definition());
118 molTable->GetConfiguration("H2O2")->SetDiffusionCoefficient(2.3e-9 * (m2 / s));
119 molTable->GetConfiguration("H2O2")->SetVanDerVaalsRadius(0.21 * nm);
120
121 // molecules extension (RITRACKS)
122
123 molTable->CreateConfiguration("HO2°", G4HO2::Definition());
124 molTable->GetConfiguration("HO2°")->SetVanDerVaalsRadius(0.21 * nm);
125
126 G4MolecularConfiguration* HO2m =
127 molTable->CreateConfiguration("HO2m", // just a tag to store and retrieve
128 // from G4MoleculeTable
129 G4HO2m,
130 -1, // charge
131 1.4e-9 * (m2 / s));
132 HO2m->SetMass(33.00396 * g / Avogadro * c_squared);
133 HO2m->SetVanDerVaalsRadius(0.25 * nm);
134
135 molTable->CreateConfiguration("Oxy", G4Oxygen::Definition());
136 molTable->GetConfiguration("Oxy")->SetVanDerVaalsRadius(0.20 * nm);
137
138 G4MolecularConfiguration* Om =
139 molTable->CreateConfiguration("Om", // just a tag to store and retrieve from
140 // G4MoleculeTable
141 G4Om,
142 -1, // charge
143 2.0e-9 * (m2 / s));
144 Om->SetMass(15.99829 * g / Avogadro * c_squared);
145 Om->SetVanDerVaalsRadius(0.25 * nm);
146
147 molTable->CreateConfiguration("O2", G4O2::Definition());
148 molTable->GetConfiguration("O2")->SetVanDerVaalsRadius(0.17 * nm);
149
150 G4MolecularConfiguration* O2m =
151 molTable->CreateConfiguration("O2m", // just a tag to store and retrieve
152 // from G4MoleculeTable
153 G4O2::Definition(),
154 -1, // charge
155 1.75e-9 * (m2 / s));
156 O2m->SetMass(31.99602 * g / Avogadro * c_squared);
157 O2m->SetVanDerVaalsRadius(0.22 * nm);
158
159 molTable->CreateConfiguration("O3", G4O3::Definition());
160 molTable->GetConfiguration("O3")->SetVanDerVaalsRadius(0.20 * nm);
161
162 G4MolecularConfiguration* O3m =
163 molTable->CreateConfiguration("O3m", // just a tag to store and retrieve
164 // from G4MoleculeTable
165 G4O3::Definition(),
166 -1, // charge
167 2.0e-9 * (m2 / s));
168 O3m->SetMass(47.99375 * g / Avogadro * c_squared);
169 O3m->SetVanDerVaalsRadius(0.20 * nm);
170
171 molTable->CreateConfiguration("H2O(B)", // just a tag to store and retrieve
172 // from G4MoleculeTable
173 G4H2O::Definition(),
174 0, // charge
175 0 * (m2 / s));
176
177 molTable->CreateConfiguration("H3Op(B)", // just a tag to store and retrieve
178 // from G4MoleculeTable
179 G4H3O::Definition(),
180 1, // charge
181 0 * (m2 / s));
182
183 molTable->CreateConfiguration("OHm(B)", // just a tag to store and retrieve
184 // from G4MoleculeTable
185 G4OHm,
186 -1, // charge
187 0 * (m2 / s));
188
189 molTable->CreateConfiguration("NoneM", G4FakeMolecule::Definition());
190 }
191
192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
193
194 void G4ChemDissociationChannels_option1::ConstructDissociationChannels()
195 {
196 //-----------------------------------
197 // Get the molecular configuration
198 auto molTable = G4MoleculeTable::Instance();
199 G4MolecularConfiguration* OH = molTable->GetConfiguration("°OH");
200 G4MolecularConfiguration* OHm = molTable->GetConfiguration("OHm");
201 G4MolecularConfiguration* e_aq = molTable->GetConfiguration("e_aq");
202 G4MolecularConfiguration* H2 = molTable->GetConfiguration("H2");
203 G4MolecularConfiguration* H3O = molTable->GetConfiguration("H3Op");
204 G4MolecularConfiguration* H = molTable->GetConfiguration("H");
205 G4MolecularConfiguration* O = molTable->GetConfiguration("Oxy");
206
207 //-------------------------------------
208 // Define the decay channels
209 G4MoleculeDefinition* water = G4H2O::Definition();
210 G4MolecularDissociationChannel* decCh1;
211 G4MolecularDissociationChannel* decCh2;
212 G4MolecularDissociationChannel* decCh3;
213 G4MolecularDissociationChannel* decCh4;
214 G4MolecularDissociationChannel* decCh5;
215
216 G4ElectronOccupancy* occ = new G4ElectronOccupancy(*(water->GetGroundStateElectronOccupancy()));
217
218 //////////////////////////////////////////////////////////
219 // EXCITATIONS //
220 //////////////////////////////////////////////////////////
221 G4DNAWaterExcitationStructure waterExcitation;
222 //--------------------------------------------------------
223 //---------------Excitation on the fifth layer------------
224
225 decCh1 = new G4MolecularDissociationChannel("A^1B_1_Relax");
226 decCh2 = new G4MolecularDissociationChannel("A^1B_1_DissociDecay");
227 // Decay 1 : OH + H
228 decCh1->SetEnergy(waterExcitation.ExcitationEnergy(0));
229 decCh1->SetProbability(0.35);
230 decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::NoDisplacement);
231
232 decCh2->AddProduct(OH);
233 decCh2->AddProduct(H);
234 decCh2->SetProbability(0.65);
235 decCh2->SetDisplacementType(G4DNAWaterDissociationDisplacer::A1B1_DissociationDecay);
236
237 occ->RemoveElectron(4, 1); // this is the transition form ground state to
238 occ->AddElectron(5, 1); // the first unoccupied orbital: A^1B_1
239
240 water->NewConfigurationWithElectronOccupancy("A^1B_1", *occ);
241 water->AddDecayChannel("A^1B_1", decCh1);
242 water->AddDecayChannel("A^1B_1", decCh2);
243
244 //--------------------------------------------------------
245 //---------------Excitation on the fourth layer-----------
246 decCh1 = new G4MolecularDissociationChannel("B^1A_1_Relax_Channel");
247 decCh2 = new G4MolecularDissociationChannel("B^1A_1_DissociDecay");
248 decCh3 = new G4MolecularDissociationChannel("B^1A_1_AutoIoni_Channel");
249 decCh4 = new G4MolecularDissociationChannel("A^1B_1_DissociDecay");
250 decCh5 = new G4MolecularDissociationChannel("B^1A_1_DissociDecay2");
251
252 // Decay 1 : energy
253 decCh1->SetEnergy(waterExcitation.ExcitationEnergy(1));
254 decCh1->SetProbability(0.175);
255
256 // Decay 2 : 2OH + H_2
257 decCh2->AddProduct(H2);
258 decCh2->AddProduct(OH);
259 decCh2->AddProduct(OH);
260 decCh2->SetProbability(0.0325);
261 decCh2->SetDisplacementType(G4DNAWaterDissociationDisplacer::B1A1_DissociationDecay);
262
263 // Decay 3 : OH + H_3Op + e_aq
264 decCh3->AddProduct(OH);
265 decCh3->AddProduct(H3O);
266 decCh3->AddProduct(e_aq);
267 decCh3->SetProbability(0.50);
268 decCh3->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
269
270 // Decay 4 : H + OH
271 decCh4->AddProduct(H);
272 decCh4->AddProduct(OH);
273 decCh4->SetProbability(0.2535);
274 decCh4->SetDisplacementType(G4DNAWaterDissociationDisplacer::A1B1_DissociationDecay);
275
276 // Decay 5 : 2H + O
277 decCh5->AddProduct(O);
278 decCh5->AddProduct(H);
279 decCh5->AddProduct(H);
280 decCh5->SetProbability(0.039);
281 decCh5->SetDisplacementType(G4DNAWaterDissociationDisplacer::B1A1_DissociationDecay2);
282
283 *occ = *(water->GetGroundStateElectronOccupancy());
284 occ->RemoveElectron(3); // this is the transition form ground state to
285 occ->AddElectron(5, 1); // the first unoccupied orbital: B^1A_1
286
287 water->NewConfigurationWithElectronOccupancy("B^1A_1", *occ);
288 water->AddDecayChannel("B^1A_1", decCh1);
289 water->AddDecayChannel("B^1A_1", decCh2);
290 water->AddDecayChannel("B^1A_1", decCh3);
291 water->AddDecayChannel("B^1A_1", decCh4);
292 water->AddDecayChannel("B^1A_1", decCh5);
293
294 //-------------------------------------------------------
295 //-------------------Excitation of 3rd layer-----------------
296 decCh1 = new G4MolecularDissociationChannel("Exci3rdLayer_AutoIoni_Channel");
297 decCh2 = new G4MolecularDissociationChannel("Exci3rdLayer_Relax_Channel");
298
299 // Decay channel 1 : : OH + H_3Op + e_aq
300 decCh1->AddProduct(OH);
301 decCh1->AddProduct(H3O);
302 decCh1->AddProduct(e_aq);
303
304 decCh1->SetProbability(0.5);
305 decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
306
307 // Decay channel 2 : energy
308 decCh2->SetEnergy(waterExcitation.ExcitationEnergy(2));
309 decCh2->SetProbability(0.5);
310
311 // Electronic configuration of this decay
312 *occ = *(water->GetGroundStateElectronOccupancy());
313 occ->RemoveElectron(2, 1);
314 occ->AddElectron(5, 1);
315
316 // Configure the water molecule
317 water->NewConfigurationWithElectronOccupancy("Exci3rdLayer", *occ);
318 water->AddDecayChannel("Exci3rdLayer", decCh1);
319 water->AddDecayChannel("Exci3rdLayer", decCh2);
320
321 //-------------------------------------------------------
322 //-------------------Excitation of 2nd layer-----------------
323 decCh1 = new G4MolecularDissociationChannel("Exci2ndLayer_AutoIoni_Channel");
324 decCh2 = new G4MolecularDissociationChannel("Exci2ndLayer_Relax_Channel");
325
326 // Decay Channel 1 : : OH + H_3Op + e_aq
327 decCh1->AddProduct(OH);
328 decCh1->AddProduct(H3O);
329 decCh1->AddProduct(e_aq);
330
331 decCh1->SetProbability(0.5);
332 decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
333
334 // Decay channel 2 : energy
335 decCh2->SetEnergy(waterExcitation.ExcitationEnergy(3));
336 decCh2->SetProbability(0.5);
337
338 *occ = *(water->GetGroundStateElectronOccupancy());
339 occ->RemoveElectron(1, 1);
340 occ->AddElectron(5, 1);
341
342 water->NewConfigurationWithElectronOccupancy("Exci2ndLayer", *occ);
343 water->AddDecayChannel("Exci2ndLayer", decCh1);
344 water->AddDecayChannel("Exci2ndLayer", decCh2);
345
346 //-------------------------------------------------------
347 //-------------------Excitation of 1st layer-----------------
348 decCh1 = new G4MolecularDissociationChannel("Exci1stLayer_AutoIoni_Channel");
349 decCh2 = new G4MolecularDissociationChannel("Exci1stLayer_Relax_Channel");
350
351 *occ = *(water->GetGroundStateElectronOccupancy());
352 occ->RemoveElectron(0, 1);
353 occ->AddElectron(5, 1);
354
355 // Decay Channel 1 : : OH + H_3Op + e_aq
356 decCh1->AddProduct(OH);
357 decCh1->AddProduct(H3O);
358 decCh1->AddProduct(e_aq);
359 decCh1->SetProbability(0.5);
360 decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
361
362 // Decay channel 2 : energy
363 decCh2->SetEnergy(waterExcitation.ExcitationEnergy(4));
364 decCh2->SetProbability(0.5);
365
366 water->NewConfigurationWithElectronOccupancy("Exci1stLayer", *occ);
367 water->AddDecayChannel("Exci1stLayer", decCh1);
368 water->AddDecayChannel("Exci1stLayer", decCh2);
369
370 /////////////////////////////////////////////////////////
371 // IONISATION //
372 /////////////////////////////////////////////////////////
373 //--------------------------------------------------------
374 //------------------- Ionisation -------------------------
375
376 decCh1 = new G4MolecularDissociationChannel("Ioni_Channel");
377
378 // Decay Channel 1 : : OH + H_3Op
379 decCh1->AddProduct(H3O);
380 decCh1->AddProduct(OH);
381 decCh1->SetProbability(1);
382 decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::Ionisation_DissociationDecay);
383
384 *occ = *(water->GetGroundStateElectronOccupancy());
385 occ->RemoveElectron(4, 1);
386 // this is a ionized h2O with a hole in its last orbital
387 water->NewConfigurationWithElectronOccupancy("Ioni5", *occ);
388 water->AddDecayChannel("Ioni5", decCh1);
389
390 *occ = *(water->GetGroundStateElectronOccupancy());
391 occ->RemoveElectron(3, 1);
392 water->NewConfigurationWithElectronOccupancy("Ioni4", *occ);
393 water->AddDecayChannel("Ioni4", new G4MolecularDissociationChannel(*decCh1));
394
395 *occ = *(water->GetGroundStateElectronOccupancy());
396 occ->RemoveElectron(2, 1);
397 water->NewConfigurationWithElectronOccupancy("Ioni3", *occ);
398 water->AddDecayChannel("Ioni3", new G4MolecularDissociationChannel(*decCh1));
399
400 *occ = *(water->GetGroundStateElectronOccupancy());
401 occ->RemoveElectron(1, 1);
402 water->NewConfigurationWithElectronOccupancy("Ioni2", *occ);
403 water->AddDecayChannel("Ioni2", new G4MolecularDissociationChannel(*decCh1));
404
405 *occ = *(water->GetGroundStateElectronOccupancy());
406 occ->RemoveElectron(0, 1);
407 water->NewConfigurationWithElectronOccupancy("Ioni1", *occ);
408 water->AddDecayChannel("Ioni1", new G4MolecularDissociationChannel(*decCh1));
409
410 //////////////////////////////////////////////////////////
411 // Dissociative Attachment //
412 //////////////////////////////////////////////////////////
413 decCh1 = new G4MolecularDissociationChannel("DissociAttachment_ch1");
414
415 // Decay 1 : OHm + H
416 decCh1->AddProduct(H2);
417 decCh1->AddProduct(OHm);
418 decCh1->AddProduct(OH);
419 decCh1->SetProbability(1);
420 decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::DissociativeAttachment);
421
422 *occ = *(water->GetGroundStateElectronOccupancy());
423 occ->AddElectron(5, 1); // H_2O^-
424
425 water->NewConfigurationWithElectronOccupancy("DissociAttachment_ch1", *occ);
426 water->AddDecayChannel("DissociAttachment_ch1", decCh1);
427
428 //////////////////////////////////////////////////////////
429 // Electron-hole recombination //
430 //////////////////////////////////////////////////////////
431 decCh1 = new G4MolecularDissociationChannel("H2Ovib_DissociDecay1");
432 decCh2 = new G4MolecularDissociationChannel("H2Ovib_DissociDecay2");
433 decCh3 = new G4MolecularDissociationChannel("H2Ovib_DissociDecay3");
434 decCh4 = new G4MolecularDissociationChannel("H2Ovib_DissociDecay4");
435
436 // Decay 1 : 2OH + H_2
437 decCh1->AddProduct(H2);
438 decCh1->AddProduct(OH);
439 decCh1->AddProduct(OH);
440 decCh1->SetProbability(0.1365);
441 decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::B1A1_DissociationDecay);
442
443 // Decay 2 : OH + H
444 decCh2->AddProduct(OH);
445 decCh2->AddProduct(H);
446 decCh2->SetProbability(0.3575);
447 decCh2->SetDisplacementType(G4DNAWaterDissociationDisplacer::A1B1_DissociationDecay);
448
449 // Decay 3 : 2H + O(3p)
450 decCh3->AddProduct(O);
451 decCh3->AddProduct(H);
452 decCh3->AddProduct(H);
453 decCh3->SetProbability(0.156);
454 decCh3->SetDisplacementType(G4DNAWaterDissociationDisplacer::B1A1_DissociationDecay2);
455
456 // Decay 4 : relaxation
457 decCh4->SetProbability(0.35);
458
459 const auto pH2Ovib = G4H2O::Definition()->NewConfiguration("H2Ovib");
460 assert(pH2Ovib != nullptr);
461
462 water->AddDecayChannel(pH2Ovib, decCh1);
463 water->AddDecayChannel(pH2Ovib, decCh2);
464 water->AddDecayChannel(pH2Ovib, decCh3);
465 water->AddDecayChannel(pH2Ovib, decCh4);
466
467 delete occ;
468 }