| Geant4 Cross Reference |
1 // ********************************************************************
2 // * License and Disclaimer *
3 // * *
4 // * The Geant4 software is copyright of the Copyright Holders of *
5 // * the Geant4 Collaboration. It is provided under the terms and *
6 // * conditions of the Geant4 Software License, included in the file *
7 // * LICENSE and available at http://cern.ch/geant4/license . These *
8 // * include a list of copyright holders. *
9 // * *
10 // * Neither the authors of this software system, nor their employing *
11 // * institutes,nor the agencies providing financial support for this *
12 // * work make any representation or warranty, express or implied, *
13 // * regarding this software system or assume any liability for its *
14 // * use. Please see the license in the file LICENSE and URL above *
15 // * for the full disclaimer and the limitation of liability. *
16 // * *
17 // * This code implementation is the result of the scientific and *
18 // * technical work of the GEANT4 collaboration. *
19 // * By using, copying, modifying or distributing the software (or *
20 // * any work based on the software) you agree to acknowledge its *
21 // * use in resulting scientific publications, and indicate your *
22 // * acceptance of all terms of the Geant4 Software license. *
23 // ********************************************************************
24 //
25 /// \file scavenger/src/EmDNAChemistry.cc
26 /// \brief Implementation of the scavenger::EmDNAChemistry class
27
28 #include "EmDNAChemistry.hh"
29
30 #include "G4DNAChemistryManager.hh"
31 #include "G4DNAElectronHoleRecombination.hh"
32 #include "G4DNAElectronSolvation.hh"
33 #include "G4DNAIRT.hh"
34 #include "G4DNAMolecularDissociation.hh"
35 #include "G4DNAMolecularIRTModel.hh"
36 #include "G4DNAMolecularReactionTable.hh"
37 #include "G4DNASancheExcitationModel.hh"
38 #include "G4DNAVibExcitation.hh"
39 #include "G4DNAWaterDissociationDisplacer.hh"
40 #include "G4DNAWaterExcitationStructure.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4ProcessManager.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4VDNAReactionModel.hh"
45 // particles
46 #include "G4Electron.hh"
47 #include "G4Electron_aq.hh"
48 #include "G4H2O.hh"
49 #include "G4MoleculeTable.hh"
50 #include "G4O2.hh"
51 #include "G4O3.hh"
52 #include "G4OH.hh"
53 #include "G4Oxygen.hh"
54 #include "G4PhysicsListHelper.hh"
55 /****/
56 #include "G4MolecularConfiguration.hh"
57 #include "G4ProcessTable.hh"
58 /****/
59 // factory
60 #include "G4PhysicsConstructorFactory.hh"
61 #include "G4ChemDissociationChannels_option1.hh"
62 #include "G4ChemDissociationChannels_option1.hh"
63 namespace scavenger
64 {
65
66 G4_DECLARE_PHYSCONSTR_FACTORY(EmDNAChemistry);
67
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
69
70 EmDNAChemistry::EmDNAChemistry()
71 : G4VUserChemistryList(true),
72 G4UImessenger(),
73 fpParserDir(new G4UIdirectory("/chem/reactionTable/")),
74 fpReactionTableNameCmd(new G4UIcmdWithAString("/chem/reactionTable/name", this))
75 {
76 G4DNAChemistryManager::Instance()->SetChemistryList(this);
77 fpParserDir->SetGuidance("Chemistry control");
78 fpReactionTableNameCmd->SetGuidance("Name of the chemical reaction table");
79 }
80
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
82
83 void EmDNAChemistry::SetNewValue(G4UIcommand* command, G4String newValue)
84 {
85 if (command == fpReactionTableNameCmd.get()) {
86 fReactionTableName = newValue;
87 }
88 }
89
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
91
92 void EmDNAChemistry::ConstructMolecule()
93 {
94 // Create the definition
95 G4ChemDissociationChannels_option1::ConstructMolecule();
96
97 auto G4NO2 = new G4MoleculeDefinition("NO_2",
98 /*mass*/ 46.0055 * g / Avogadro * c_squared,
99 /*D*/ 1.7e-9 * (m * m / s),
100 /*charge*/ 0,
101 /*electronL*/ 0,
102 /*radius*/ 0.35 * nm); // can be adjusted
103
104 auto G4NO3 = new G4MoleculeDefinition("NO_3",
105 /*mass*/ 62.0049 * g / Avogadro * c_squared,
106 /*D*/ 1.7e-9 * (m * m / s),
107 /*charge*/ 0,
108 /*electronL*/ 0,
109 /*radius*/ 0.35 * nm); // can be adjusted
110 //____________________________________________________________________________
111 // Note: Parameters Value changed according to Plante Paper
112
113 G4MoleculeTable::Instance()->CreateConfiguration("O2(B)", // just a tag to store and retrieve
114 // from G4MoleculeTable
115 G4O2::Definition(),
116 0, // charge
117 0 * (m2 / s));
118
119 G4MoleculeTable::Instance()->CreateConfiguration("NO2m(B)", // just a tag to store and retrieve
120 // from G4MoleculeTable
121 G4NO2,
122 -1, // charge
123 0 * (m2 / s));
124
125 G4MoleculeTable::Instance()->CreateConfiguration("NO3m(B)", // just a tag to store and retrieve
126 // from G4MoleculeTable
127 G4NO3,
128 -1, // charge
129 0 * (m2 / s));
130
131 }
132
133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
134
135 void EmDNAChemistry::ConstructDissociationChannels()
136 {
137 //-----------------------------------
138 // Get the molecular configuration
139 auto OH = G4MoleculeTable::Instance()->GetConfiguration("°OH");
140 auto OHm = G4MoleculeTable::Instance()->GetConfiguration("OHm");
141 auto e_aq = G4MoleculeTable::Instance()->GetConfiguration("e_aq");
142 auto H2 = G4MoleculeTable::Instance()->GetConfiguration("H2");
143 auto H3O = G4MoleculeTable::Instance()->GetConfiguration("H3Op");
144 auto H = G4MoleculeTable::Instance()->GetConfiguration("H");
145
146 //-------------------------------------
147 // Define the decay channels
148 auto* water = G4H2O::Definition();
149 G4MolecularDissociationChannel* decCh1;
150 G4MolecularDissociationChannel* decCh2;
151
152 auto occ = new G4ElectronOccupancy(*(water->GetGroundStateElectronOccupancy()));
153
154 //////////////////////////////////////////////////////////
155 // EXCITATIONS //
156 //////////////////////////////////////////////////////////
157 G4DNAWaterExcitationStructure waterExcitation;
158 //--------------------------------------------------------
159 //---------------Excitation on the fifth layer------------
160
161 decCh1 = new G4MolecularDissociationChannel("A^1B_1_Relaxation");
162 decCh2 = new G4MolecularDissociationChannel("A^1B_1_DissociativeDecay");
163 // Decay 1 : OH + H
164 decCh1->SetEnergy(waterExcitation.ExcitationEnergy(0));
165 decCh1->SetProbability(0.35);
166 decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::NoDisplacement);
167
168 decCh2->AddProduct(OH);
169 decCh2->AddProduct(H);
170 decCh2->SetProbability(0.65);
171 decCh2->SetDisplacementType(G4DNAWaterDissociationDisplacer::A1B1_DissociationDecay);
172
173 // water->AddExcitedState("A^1B_1");
174 occ->RemoveElectron(4, 1); // this is the transition form ground state to
175 occ->AddElectron(5, 1); // the first unoccupied orbital: A^1B_1
176
177 water->NewConfigurationWithElectronOccupancy("A^1B_1", *occ);
178 water->AddDecayChannel("A^1B_1", decCh1);
179 water->AddDecayChannel("A^1B_1", decCh2);
180
181 //--------------------------------------------------------
182 //---------------Excitation on the fourth layer-----------
183 decCh1 = new G4MolecularDissociationChannel("B^1A_1_Relaxation_Channel");
184 decCh2 = new G4MolecularDissociationChannel("B^1A_1_DissociativeDecay");
185 auto decCh3 = new G4MolecularDissociationChannel("B^1A_1_AutoIonisation_Channel");
186
187 // Decay 1 : energy
188 decCh1->SetEnergy(waterExcitation.ExcitationEnergy(1));
189 decCh1->SetProbability(0.3);
190
191 // Decay 2 : 2OH + H_2
192 decCh2->AddProduct(H2);
193 decCh2->AddProduct(OH);
194 decCh2->AddProduct(OH);
195 decCh2->SetProbability(0.15);
196 decCh2->SetDisplacementType(G4DNAWaterDissociationDisplacer::B1A1_DissociationDecay);
197
198 // Decay 3 : OH + H_3Op + e_aq
199 decCh3->AddProduct(OH);
200 decCh3->AddProduct(H3O);
201 decCh3->AddProduct(e_aq);
202 decCh3->SetProbability(0.55);
203 decCh3->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
204
205 *occ = *(water->GetGroundStateElectronOccupancy());
206 occ->RemoveElectron(3); // this is the transition form ground state to
207 occ->AddElectron(5, 1); // the first unoccupied orbital: B^1A_1
208
209 water->NewConfigurationWithElectronOccupancy("B^1A_1", *occ);
210 water->AddDecayChannel("B^1A_1", decCh1);
211 water->AddDecayChannel("B^1A_1", decCh2);
212 water->AddDecayChannel("B^1A_1", decCh3);
213
214 //-------------------------------------------------------
215 //-------------------Excitation of 3rd layer-----------------
216 decCh1 = new G4MolecularDissociationChannel("Excitation3rdLayer_AutoIonisation_Channel");
217 decCh2 = new G4MolecularDissociationChannel("Excitation3rdLayer_Relaxation_Channel");
218
219 // Decay channel 1 : : OH + H_3Op + e_aq
220 decCh1->AddProduct(OH);
221 decCh1->AddProduct(H3O);
222 decCh1->AddProduct(e_aq);
223
224 decCh1->SetProbability(0.5);
225 decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
226
227 // Decay channel 2 : energy
228 decCh2->SetEnergy(waterExcitation.ExcitationEnergy(2));
229 decCh2->SetProbability(0.5);
230
231 // Electronic configuration of this decay
232 *occ = *(water->GetGroundStateElectronOccupancy());
233 occ->RemoveElectron(2, 1);
234 occ->AddElectron(5, 1);
235
236 // Configure the water molecule
237 water->NewConfigurationWithElectronOccupancy("Excitation3rdLayer", *occ);
238 water->AddDecayChannel("Excitation3rdLayer", decCh1);
239 water->AddDecayChannel("Excitation3rdLayer", decCh2);
240
241 //-------------------------------------------------------
242 //-------------------Excitation of 2nd layer-----------------
243 decCh1 = new G4MolecularDissociationChannel("Excitation2ndLayer_AutoIonisation_Channel");
244 decCh2 = new G4MolecularDissociationChannel("Excitation2ndLayer_Relaxation_Channel");
245
246 // Decay Channel 1 : : OH + H_3Op + e_aq
247 decCh1->AddProduct(OH);
248 decCh1->AddProduct(H3O);
249 decCh1->AddProduct(e_aq);
250
251 decCh1->SetProbability(0.5);
252 decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
253
254 // Decay channel 2 : energy
255 decCh2->SetEnergy(waterExcitation.ExcitationEnergy(3));
256 decCh2->SetProbability(0.5);
257
258 *occ = *(water->GetGroundStateElectronOccupancy());
259 occ->RemoveElectron(1, 1);
260 occ->AddElectron(5, 1);
261
262 water->NewConfigurationWithElectronOccupancy("Excitation2ndLayer", *occ);
263 water->AddDecayChannel("Excitation2ndLayer", decCh1);
264 water->AddDecayChannel("Excitation2ndLayer", decCh2);
265
266 //-------------------------------------------------------
267 //-------------------Excitation of 1st layer-----------------
268 decCh1 = new G4MolecularDissociationChannel("Excitation1stLayer_AutoIonisation_Channel");
269 decCh2 = new G4MolecularDissociationChannel("Excitation1stLayer_Relaxation_Channel");
270
271 *occ = *(water->GetGroundStateElectronOccupancy());
272 occ->RemoveElectron(0, 1);
273 occ->AddElectron(5, 1);
274
275 // Decay Channel 1 : : OH + H_3Op + e_aq
276 decCh1->AddProduct(OH);
277 decCh1->AddProduct(H3O);
278 decCh1->AddProduct(e_aq);
279 decCh1->SetProbability(0.5);
280 decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
281
282 // Decay channel 2 : energy
283 decCh2->SetEnergy(waterExcitation.ExcitationEnergy(4));
284 decCh2->SetProbability(0.5);
285
286 water->NewConfigurationWithElectronOccupancy("Excitation1stLayer", *occ);
287 water->AddDecayChannel("Excitation1stLayer", decCh1);
288 water->AddDecayChannel("Excitation1stLayer", decCh2);
289
290 /////////////////////////////////////////////////////////
291 // IONISATION //
292 /////////////////////////////////////////////////////////
293 //--------------------------------------------------------
294 //------------------- Ionisation -------------------------
295
296 decCh1 = new G4MolecularDissociationChannel("Ionisation_Channel");
297
298 // Decay Channel 1 : : OH + H_3Op
299 decCh1->AddProduct(H3O);
300 decCh1->AddProduct(OH);
301 decCh1->SetProbability(1);
302 decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::Ionisation_DissociationDecay);
303
304 *occ = *(water->GetGroundStateElectronOccupancy());
305 occ->RemoveElectron(4, 1);
306 // this is a ionized h2O with a hole in its last orbital
307 water->NewConfigurationWithElectronOccupancy("Ionisation5", *occ);
308 water->AddDecayChannel("Ionisation5", decCh1);
309
310 *occ = *(water->GetGroundStateElectronOccupancy());
311 occ->RemoveElectron(3, 1);
312 water->NewConfigurationWithElectronOccupancy("Ionisation4", *occ);
313 water->AddDecayChannel("Ionisation4", new G4MolecularDissociationChannel(*decCh1));
314
315 *occ = *(water->GetGroundStateElectronOccupancy());
316 occ->RemoveElectron(2, 1);
317 water->NewConfigurationWithElectronOccupancy("Ionisation3", *occ);
318 water->AddDecayChannel("Ionisation3", new G4MolecularDissociationChannel(*decCh1));
319
320 *occ = *(water->GetGroundStateElectronOccupancy());
321 occ->RemoveElectron(1, 1);
322 water->NewConfigurationWithElectronOccupancy("Ionisation2", *occ);
323 water->AddDecayChannel("Ionisation2", new G4MolecularDissociationChannel(*decCh1));
324
325 *occ = *(water->GetGroundStateElectronOccupancy());
326 occ->RemoveElectron(0, 1);
327 water->NewConfigurationWithElectronOccupancy("Ionisation1", *occ);
328 water->AddDecayChannel("Ionisation1", new G4MolecularDissociationChannel(*decCh1));
329
330 //////////////////////////////////////////////////////////
331 // Dissociative Attachment //
332 //////////////////////////////////////////////////////////
333 decCh1 = new G4MolecularDissociationChannel("DissociativeAttachment");
334
335 // Decay 1 : 2OH + H_2
336 decCh1->AddProduct(H2);
337 decCh1->AddProduct(OHm);
338 decCh1->AddProduct(OH);
339 decCh1->SetProbability(1);
340 decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::DissociativeAttachment);
341
342 *occ = *(water->GetGroundStateElectronOccupancy());
343 occ->AddElectron(5, 1); // H_2O^-
344 water->NewConfigurationWithElectronOccupancy("DissociativeAttachment", *occ);
345 water->AddDecayChannel("DissociativeAttachment", decCh1);
346
347 //////////////////////////////////////////////////////////
348 // Electron-hole recombination //
349 //////////////////////////////////////////////////////////
350 decCh1 = new G4MolecularDissociationChannel("H2Ovib_DissociationDecay1");
351 decCh2 = new G4MolecularDissociationChannel("H2Ovib_DissociationDecay2");
352 decCh3 = new G4MolecularDissociationChannel("H2Ovib_DissociationDecay3");
353
354 // Decay 1 : 2OH + H_2
355 decCh1->AddProduct(H2);
356 decCh1->AddProduct(OH);
357 decCh1->AddProduct(OH);
358 decCh1->SetProbability(0.15);
359 decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::B1A1_DissociationDecay);
360
361 // Decay 2 : OH + H
362 decCh2->AddProduct(OH);
363 decCh2->AddProduct(H);
364 decCh2->SetProbability(0.55);
365 decCh2->SetDisplacementType(G4DNAWaterDissociationDisplacer::A1B1_DissociationDecay);
366
367 // Decay 3 : relaxation
368 decCh3->SetProbability(0.30);
369
370 const auto pH2Ovib = G4H2O::Definition()->NewConfiguration("H2Ovib");
371 assert(pH2Ovib != nullptr);
372
373 water->AddDecayChannel(pH2Ovib, decCh1);
374 water->AddDecayChannel(pH2Ovib, decCh2);
375 water->AddDecayChannel(pH2Ovib, decCh3);
376
377 delete occ;
378 }
379
380 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
381
382 void EmDNAChemistry::ConstructReactionTable(G4DNAMolecularReactionTable* theReactionTable)
383 {
384 // Construct chemical reactions from user input file
385 ParserChemReaction parser;
386 parser.ReadReactionFile(fReactionTableName);
387
388 auto ListReactant1 = parser.GetListReactant1();
389 auto ListReactant2 = parser.GetListReactant2();
390 auto ListProduct = parser.GetListProduct();
391 auto ListRate = parser.GetListRate();
392
393 const G4int Ntype = 5;
394 for (size_t i = 0; i < Ntype; i++) {
395 for (size_t j = 0; j < ListReactant1[i].size(); j++) {
396 G4MolecularConfiguration* Reactant1 =
397 G4MoleculeTable::Instance()->GetConfiguration(ListReactant1[i][j]);
398 G4MolecularConfiguration* Reactant2 =
399 G4MoleculeTable::Instance()->GetConfiguration(ListReactant2[i][j]);
400
401 auto pReactionData = new G4DNAMolecularReactionData(ListRate[i][j], Reactant1, Reactant2);
402
403 for (size_t k = 0; k < ListProduct[i][j].size(); k++) {
404 G4MolecularConfiguration* Product =
405 G4MoleculeTable::Instance()->GetConfiguration(ListProduct[i][j][k]);
406 pReactionData->AddProduct(Product);
407 }
408 // Reaction type II and IV
409 if (i == 1 || i == 3) {
410 pReactionData->SetReactionType(1);
411 }
412
413 theReactionTable->SetReaction(pReactionData);
414 }
415 }
416 }
417
418 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
419
420 void EmDNAChemistry::ConstructProcess()
421 {
422 auto ph = G4PhysicsListHelper::GetPhysicsListHelper();
423 G4VProcess* process =
424 G4ProcessTable::GetProcessTable()->FindProcess("e-_G4DNAVibExcitation", "e-");
425 if (process) {
426 auto vibExcitation = (G4DNAVibExcitation*)process;
427 G4VEmModel* model = vibExcitation->EmModel();
428 auto sancheExcitationMod = dynamic_cast<G4DNASancheExcitationModel*>(model);
429 if (sancheExcitationMod) {
430 sancheExcitationMod->ExtendLowEnergyLimit(0.025 * eV);
431 }
432 }
433 process = G4ProcessTable::GetProcessTable()->FindProcess("e-_G4DNAElectronSolvation", "e-");
434
435 if (process == nullptr) {
436 ph->RegisterProcess(new G4DNAElectronSolvation("e-_G4DNAElectronSolvation"),
437 G4Electron::Definition());
438 }
439 auto* theMoleculeTable = G4MoleculeTable::Instance();
440 G4MoleculeDefinitionIterator iterator = theMoleculeTable->GetDefintionIterator();
441 iterator.reset();
442 while (iterator()) {
443 G4MoleculeDefinition* moleculeDef = iterator.value();
444 if (moleculeDef == G4H2O::Definition()) {
445 moleculeDef->GetProcessManager()->AddRestProcess(new G4DNAElectronHoleRecombination(), 2);
446
447 auto dissociationProcess = new G4DNAMolecularDissociation("H2O_DNAMolecularDecay");
448 dissociationProcess->SetDisplacer(moleculeDef, new G4DNAWaterDissociationDisplacer);
449 // dissociationProcess->SetVerboseLevel(3);
450 moleculeDef->GetProcessManager()->AddRestProcess(dissociationProcess, 1);
451 }
452 }
453 G4DNAChemistryManager::Instance()->Initialize();
454 }
455
456 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
457
458 void EmDNAChemistry::ConstructTimeStepModel(G4DNAMolecularReactionTable*
459 /*reactionTable*/)
460 {
461 auto irt = new G4DNAMolecularIRTModel();
462 RegisterTimeStepModel(irt, 0);
463 }
464
465 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
466
467 } // namespace scavenger
468