| 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 #include "EmDNAChemistry.hh"
27
28 #include "G4DNAChemistryManager.hh"
29 #include "G4DNAWaterDissociationDisplacer.hh"
30 #include "G4ProcessManager.hh"
31 #include "G4SystemOfUnits.hh"
32
33 // *** Processes and models for Geant4-DNA
34
35 #include "BoundedBrownianAction.hh"
36
37 #include "G4DNABrownianTransportation.hh"
38 #include "G4DNAElectronHoleRecombination.hh"
39 #include "G4DNAElectronSolvation.hh"
40 #include "G4DNAMolecularDissociation.hh"
41 #include "G4DNAMolecularReactionTable.hh"
42 #include "G4DNAMolecularStepByStepModel.hh"
43 #include "G4DNASancheExcitationModel.hh"
44 #include "G4DNASmoluchowskiReactionModel.hh"
45 #include "G4DNAVibExcitation.hh"
46 // particles
47
48 #include "G4Electron.hh"
49 #include "G4Electron_aq.hh"
50 #include "G4H2O.hh"
51 #include "G4H2O2.hh"
52 #include "G4H3O.hh"
53 #include "G4HO2.hh"
54 #include "G4Hydrogen.hh"
55 #include "G4MoleculeTable.hh"
56 #include "G4O2.hh"
57 #include "G4O3.hh"
58 #include "G4OH.hh"
59 #include "G4Oxygen.hh"
60 #include "G4PhysicsListHelper.hh"
61 /****/
62 #include "G4DNAMoleculeEncounterStepper.hh"
63 #include "G4DNAScavengerProcess.hh"
64 #include "G4MolecularConfiguration.hh"
65 #include "G4ProcessTable.hh"
66 #include "G4VChemistryWorld.hh"
67 /****/
68 #include "G4ChemicalMoleculeFinder.hh"
69 // factory
70 #include "ChemOxygenWaterBuilder.hh"
71 #include "ChemPureWaterBuilder.hh"
72
73 #include "G4ChemDissociationChannels_option1.hh"
74 #include "G4PhysicsConstructorFactory.hh"
75
76 G4_DECLARE_PHYSCONSTR_FACTORY(EmDNAChemistry);
77
78 EmDNAChemistry::EmDNAChemistry() : G4VUserChemistryList(true)
79 {
80 G4DNAChemistryManager::Instance()->SetChemistryList(this);
81 }
82
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
84
85 EmDNAChemistry::~EmDNAChemistry() = default;
86
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88
89 void EmDNAChemistry::ConstructMolecule()
90 {
91 G4ChemDissociationChannels_option1::ConstructMolecule();
92 auto table = G4MoleculeTable::Instance();
93
94 auto H3OpB = table->GetConfiguration("H3Op(B)");
95 H3OpB->SetDiffusionCoefficient(9.46e-9 * (m2 / s));
96
97 auto OHm = table->GetConfiguration("OHm(B)");
98 OHm->SetDiffusionCoefficient(5.3e-9 * (m2 / s));
99 table->CreateConfiguration("H2O", G4H2O::Definition());
100
101 auto G4NO2 = new G4MoleculeDefinition("NO_2", /*mass*/ 30,
102 /*D*/ 0 * (m * m / s),
103 /*charge*/ 0,
104 /*electronL*/ 0,
105 /*radius*/ 0.17 * nm); // should be corrected
106
107 auto G4NO3 = new G4MoleculeDefinition("NO_3", /*mass*/ 38,
108 /*D*/ 0 * (m * m / s),
109 /*charge*/ 0,
110 /*electronL*/ 0,
111 /*radius*/ 0.17 * nm); // should be corrected
112
113 table->CreateConfiguration("NO2", G4NO2);
114 table->CreateConfiguration("NO2m", G4NO2,
115 -1, // charge
116 0 * (m2 / s));
117 table->CreateConfiguration("NO2mm", G4NO2,
118 -2, // charge
119 0 * (m2 / s));
120
121 table->CreateConfiguration("NO3m", G4NO3,
122 -1, // charge
123 0 * (m2 / s));
124
125 table->CreateConfiguration("NO3mm", G4NO3,
126 -2, // charge
127 0 * (m2 / s));
128
129 // FrickeDosimeter
130 auto G4Fe = new G4MoleculeDefinition("Fe",
131 /*mass*/ 55.84 * g / Avogadro * c_squared,
132 /*D*/ 0 * (m * m / s),
133 /*charge*/ 0,
134 /*electronL*/ 0,
135 /*radius*/ 0.35 * nm); // can be adjusted
136
137 table->CreateConfiguration("Fe0", G4Fe);
138
139 table->CreateConfiguration("Feppp", G4Fe,
140 3, // charge
141 4.86e-10 * (m2 / s)); // Michael Spiro* and Andrew M. Creeth
142
143 table->CreateConfiguration("Fepp", G4Fe,
144 2, // charge
145 5.78e-10 * (m2 / s));
146 // HSO4-
147 auto G4HSO4 = new G4MoleculeDefinition("HSO4",
148 /*mass*/ 55.84 * g / Avogadro * c_squared,
149 /*D*/ 0 * (m * m / s),
150 /*charge*/ 0,
151 /*electronL*/ 0,
152 /*radius*/ 0.35 * nm); // can be adjusted
153 table->CreateConfiguration("HSO4m", G4HSO4,
154 -1, // charge
155 0 * (m2 / s));
156
157 // SO4-
158 auto G4SO4 = new G4MoleculeDefinition("SO4",
159 /*mass*/ 55.84 * g / Avogadro * c_squared,
160 /*D*/ 0 * (m * m / s),
161 /*charge*/ 0,
162 /*electronL*/ 0,
163 /*radius*/ 0.35 * nm); // can be adjusted
164 table->CreateConfiguration("SO4m", G4SO4,
165 -1, // charge
166 0 * (m2 / s));
167 }
168
169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
170
171 void EmDNAChemistry::ConstructDissociationChannels()
172 {
173 G4ChemDissociationChannels_option1::ConstructDissociationChannels();
174 }
175
176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
177
178 void EmDNAChemistry::ConstructReactionTable(G4DNAMolecularReactionTable* pReactionTable)
179 {
180 ChemOxygenWaterBuilder::OxygenScavengerReaction(pReactionTable);
181 ChemOxygenWaterBuilder::SecondOrderReactionExtended(pReactionTable);
182 ChemPureWaterBuilder::WaterScavengerReaction(pReactionTable);
183 }
184
185 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
186
187 void EmDNAChemistry::ConstructProcess()
188 {
189 auto table = G4MoleculeTable::Instance();
190 auto O2 = table->GetConfiguration("O2");
191 auto O2m = table->GetConfiguration("O2m");
192 auto HO2 = table->GetConfiguration("HO2°");
193
194 auto e_aq = table->GetConfiguration("e_aq");
195 auto OH = table->GetConfiguration("°OH");
196 auto OHm = table->GetConfiguration("OHm");
197
198 auto NO2 = table->GetConfiguration("NO2");
199 auto NO2m = table->GetConfiguration("NO2m");
200 auto NO2mm = table->GetConfiguration("NO2mm");
201 auto NO3m = table->GetConfiguration("NO3m");
202 auto NO3mm = table->GetConfiguration("NO3mm");
203
204 auto H2O2 = table->GetConfiguration("H2O2");
205 auto H = table->GetConfiguration("H");
206
207 auto* H3OpB = table->GetConfiguration("H3Op(B)");
208 auto* OHmB = table->GetConfiguration("OHm(B)");
209 auto* HO2m = table->GetConfiguration("HO2m");
210 auto* Om = table->GetConfiguration("Om");
211 auto* O3m = table->GetConfiguration("O3m");
212 auto* H3Op = table->GetConfiguration("H3Op");
213
214 fpChemistryWorld->ConstructChemistryComponents();
215 auto confinedBox = fpChemistryWorld->GetChemistryBoundary();
216
217 auto* ph = G4PhysicsListHelper::GetPhysicsListHelper();
218
219 //===============================================================
220 // Extend vibrational to low energy
221 // Anyway, solvation of electrons is taken into account from 7.4 eV
222 // So below this threshold, for now, no accurate modeling is done
223 //
224 G4VProcess* process =
225 G4ProcessTable::GetProcessTable()->FindProcess("e-_G4DNAVibExcitation", "e-");
226
227 if (process) {
228 auto vibExcitation = (G4DNAVibExcitation*)process;
229 G4VEmModel* model = vibExcitation->EmModel();
230 auto sancheExcitationMod = dynamic_cast<G4DNASancheExcitationModel*>(model);
231 if (sancheExcitationMod) {
232 sancheExcitationMod->ExtendLowEnergyLimit(0.025 * eV);
233 }
234 }
235
236 //===============================================================
237 // *** Electron Solvatation ***
238 //
239 process = G4ProcessTable::GetProcessTable()->FindProcess("e-_G4DNAElectronSolvation", "e-");
240
241 if (process == nullptr) {
242 ph->RegisterProcess(new G4DNAElectronSolvation("e-_G4DNAElectronSolvation"),
243 G4Electron::Definition());
244 }
245
246 //===============================================================
247 // Define processes for molecules
248 //
249 auto* theMoleculeTable = G4MoleculeTable::Instance();
250 auto iterator = theMoleculeTable->GetDefintionIterator();
251 iterator.reset();
252
253 while (iterator()) {
254 auto* moleculeDef = iterator.value();
255
256 if (moleculeDef != G4H2O::Definition()) {
257 auto brown = new G4DNABrownianTransportation("BrowianTransportation");
258 // hoang exp
259 auto brownTransport = new BoundedBrownianAction();
260 brownTransport->SetBoundary(*confinedBox);
261 brown->SetUserBrownianAction(brownTransport);
262 // hoang exp
263
264 ph->RegisterProcess(brown, moleculeDef);
265 }
266 else {
267 moleculeDef->GetProcessManager()->AddRestProcess(new G4DNAElectronHoleRecombination(), 2);
268 auto brownTransport = new BoundedBrownianAction();
269 brownTransport->SetBoundary(*confinedBox);
270 auto dissociationProcess = new G4DNAMolecularDissociation("H2O_DNAMolecularDecay", fDecay);
271 dissociationProcess->SetUserBrownianAction(brownTransport);
272 dissociationProcess->SetDisplacer(moleculeDef, new G4DNAWaterDissociationDisplacer);
273 moleculeDef->GetProcessManager()->AddRestProcess(dissociationProcess, 1);
274 }
275
276 if (moleculeDef == G4Hydrogen::Definition()) {
277 // O2
278 auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
279 //------------------------------------------------------------------
280 // H + O2(B) -> HO2
281 auto reactionData = new G4DNAMolecularReactionData(1.3e10 * (1e-3 * m3 / (mole * s)), H, O2);
282 reactionData->AddProduct(HO2);
283 scanvergerProcess->SetReaction(H, reactionData);
284 //------------------------------------------------------------------
285 // H + OH-(B) -> H2O + eaq- 2.49e3 / s
286 reactionData = new G4DNAMolecularReactionData(2.49e7 * (1e-3 * m3 / (mole * s)), H,
287 OHmB); // 2.51e7 (H + OH-)* 1e-7 (pH) = 2.48e0
288 reactionData->AddProduct(e_aq);
289 scanvergerProcess->SetReaction(H, reactionData);
290
291 // H2O2
292 //------------------------------------------------------------------
293 // H + H202 -> OH + H20
294 // reactionData = new G4DNAMolecularReactionData(
295 // 9.0e7 * (1e-3 * m3 / (mole * s)), H,H2O2);
296 // reactionData->AddProduct(OH);
297 // scanvergerProcess->SetReaction(H,reactionData);
298 ph->RegisterProcess(scanvergerProcess, moleculeDef);
299 }
300 if (moleculeDef == G4Electron_aq::Definition()) {
301 auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
302 G4DNAMolecularReactionData* reactionData = nullptr;
303 //------------------------------------------------------------------
304 // e_aq + O2(B) -> O2-
305 reactionData = new G4DNAMolecularReactionData(2.3e10 * (1e-3 * m3 / (mole * s)), e_aq, O2);
306 reactionData->AddProduct(O2m);
307 scanvergerProcess->SetReaction(e_aq, reactionData);
308 //------------------------------------------------------------------
309 // eaq- + H3O+(B) -> H + H2O 2.09e3 / s
310 reactionData =
311 new G4DNAMolecularReactionData(2.25e10 * (1e-3 * m3 / (mole * s)), e_aq,
312 H3OpB); // 2.11e10 (e_aq + H3O+) * 1.0e-7 (Ph=7) = 2.09e3
313 reactionData->AddProduct(H);
314 scanvergerProcess->SetReaction(e_aq, reactionData);
315 //------------------------------------------------------------------
316 // e_aq + NO2- -> NO2--
317 reactionData = new G4DNAMolecularReactionData(3.5e9 * (1e-3 * m3 / (mole * s)), e_aq, NO2m);
318 reactionData->AddProduct(NO2mm);
319 scanvergerProcess->SetReaction(e_aq, reactionData);
320 //------------------------------------------------------------------
321 // e_aq + NO3- -> NO3--
322 reactionData = new G4DNAMolecularReactionData(9.7e9 * (1e-3 * m3 / (mole * s)), e_aq, NO3m);
323 reactionData->AddProduct(NO3mm);
324 scanvergerProcess->SetReaction(e_aq, reactionData);
325 //------------------------------------------------------------------
326
327 // H2O2 + e aq → OHm + OH
328 // reactionData = new G4DNAMolecularReactionData(
329 // 1.1e10 * (1e-3 * m3 / (mole * s)), e_aq, H2O2);//or
330 // reactionData->AddProduct(OHm);
331 // reactionData->AddProduct(OH);
332 // scanvergerProcess->SetReaction(e_aq,reactionData);
333
334 ph->RegisterProcess(scanvergerProcess, moleculeDef);
335 }
336 if (moleculeDef == G4O2::Definition()) {
337 auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
338 G4DNAMolecularReactionData* reactionData = nullptr;
339 //------------------------------------------------------------------
340 // O2- + H3O+(B) -> HO2 + H2O 4.73e3 / s
341 reactionData =
342 new G4DNAMolecularReactionData(4.78e10 * (1e-3 * m3 / (mole * s)), O2m,
343 H3OpB); // 4.78e10(O2- + H3O+) * 1e-7(pH7) = 4.73e3
344 reactionData->AddProduct(HO2);
345 scanvergerProcess->SetReaction(O2m, reactionData);
346 ph->RegisterProcess(scanvergerProcess, moleculeDef);
347 }
348 if (moleculeDef == G4ParticleTable::GetParticleTable()->FindParticle("OHm")) {
349 auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
350 G4DNAMolecularReactionData* reactionData = nullptr;
351 //------------------------------------------------------------------
352 // OH- + H3O+(B) -> 2H2O 1.11e4 / s
353 reactionData =
354 new G4DNAMolecularReactionData(1.13e11 * (1e-3 * m3 / (mole * s)), OHm,
355 H3OpB); // 1.13e11 (H3O+ + OH-) * 1e-7 (pH=7) =1.12e4
356 scanvergerProcess->SetReaction(OHm, reactionData);
357 ph->RegisterProcess(scanvergerProcess, moleculeDef);
358 }
359 if (moleculeDef == G4OH::Definition()) {
360 auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
361 G4DNAMolecularReactionData* reactionData = nullptr;
362
363 //------------------------------------------------------------------
364 // OH + OH-(B) -> O- + H2O 6.24e2 / s
365 reactionData =
366 new G4DNAMolecularReactionData(1.27e10 * (1e-3 * m3 / (mole * s)), OH,
367 OHmB); // 6.30e9 (OH + OH-) * 1e-7 (pH) = 6.24e2
368 reactionData->AddProduct(Om);
369 scanvergerProcess->SetReaction(OH, reactionData);
370
371 //------------------------------------------------------------------
372 // OH + NO2- -> NO2 + OH-
373 reactionData = new G4DNAMolecularReactionData(8e9 * (1e-3 * m3 / (mole * s)), OH, NO2m);
374 reactionData->AddProduct(NO2);
375 reactionData->AddProduct(OHm);
376 scanvergerProcess->SetReaction(OH, reactionData);
377 ph->RegisterProcess(scanvergerProcess, moleculeDef);
378 }
379 if (moleculeDef == G4ParticleTable::GetParticleTable()->FindParticle("HO_2m")) {
380 auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
381 G4DNAMolecularReactionData* reactionData = nullptr;
382 //------------------------------------------------------------------
383 // HO2- + H3O+(B) -> H2O2 + H2O 4.98e3 / s
384 reactionData =
385 new G4DNAMolecularReactionData(4.78e10 * (1e-3 * m3 / (mole * s)), HO2m,
386 H3OpB); // 5.00e10 (H3O+ + HO2-) * 1e-7(pH) = 4.95e3
387 reactionData->AddProduct(H2O2);
388 scanvergerProcess->SetReaction(HO2m, reactionData);
389 ph->RegisterProcess(scanvergerProcess, moleculeDef);
390 }
391
392 if (moleculeDef == G4HO2::Definition()) {
393 auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
394 G4DNAMolecularReactionData* reactionData = nullptr;
395 //------------------------------------------------------------------
396 // HO2 + OH-(B) -> O2- + H2O 6.24e2 / s
397 reactionData = new G4DNAMolecularReactionData(1.27e10 * (1e-3 * m3 / (mole * s)), HO2,
398 OHmB); // 6.30e9(HO2 + OH-)*1e-7 (pH) = 6.24e2
399 reactionData->AddProduct(O2m);
400 scanvergerProcess->SetReaction(HO2, reactionData);
401 //------------------------------------------------------------------
402 ph->RegisterProcess(scanvergerProcess, moleculeDef);
403 }
404 if (moleculeDef == G4Oxygen::Definition()) {
405 auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
406 G4DNAMolecularReactionData* reactionData = nullptr;
407 //------------------------------------------------------------------
408 // O- + H3O+(B) -> OH + H2O 4.73e3 / s
409 reactionData =
410 new G4DNAMolecularReactionData(4.78e10 * (1e-3 * m3 / (mole * s)), Om,
411 H3OpB); // 4.78e10 (H3O+ + O2-) * 1e-7(pH) = 4.73e3
412 reactionData->AddProduct(OH);
413 scanvergerProcess->SetReaction(Om, reactionData);
414 ph->RegisterProcess(scanvergerProcess, moleculeDef);
415 }
416 if (moleculeDef == G4O3::Definition()) {
417 auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
418 G4DNAMolecularReactionData* reactionData = nullptr;
419 //------------------------------------------------------------------
420 // O3- + H3O+(B) -> OH + O2 + H2O 8.91e3 / s
421 reactionData =
422 new G4DNAMolecularReactionData(9.0e10 * (1e-3 * m3 / (mole * s)), O3m,
423 H3OpB); // 9.0e10 (O3- + H3O+) * 1e-7(pH) = 8.91e3
424 reactionData->AddProduct(OH);
425 reactionData->AddProduct(O2);
426 //------------------------------------------------------------------
427 scanvergerProcess->SetReaction(O3m, reactionData);
428 ph->RegisterProcess(scanvergerProcess, moleculeDef);
429 }
430 if (moleculeDef == G4H3O::Definition()) {
431 auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
432 G4DNAMolecularReactionData* reactionData = nullptr;
433 //------------------------------------------------------------------
434 // H3O+ + OH-(B) -> 2H2O 1.11e4 / s
435 reactionData =
436 new G4DNAMolecularReactionData(1.13e11 * (1e-3 * m3 / (mole * s)), H3Op,
437 OHmB); // 1.13e11 (H3O+ + OH-) * 1e-7 (pH=7) = 1.12e4
438 scanvergerProcess->SetReaction(H3Op, reactionData);
439 ph->RegisterProcess(scanvergerProcess, moleculeDef);
440 }
441 if (moleculeDef == G4H2O2::Definition()) {
442 auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
443 G4DNAMolecularReactionData* reactionData = nullptr;
444 //------------------------------------------------------------------
445 // H2O2 + OH-(B) -> HO2- + H2O 4.66e2 / s
446 reactionData =
447 new G4DNAMolecularReactionData(1.27e10 * (1e-3 * m3 / (mole * s)), H2O2,
448 OHmB); // 4.71e8 (H2O2 + OH-) * 1e-7 (pH) = 4.66e1
449 reactionData->AddProduct(HO2m);
450 scanvergerProcess->SetReaction(H2O2, reactionData);
451 ph->RegisterProcess(scanvergerProcess, moleculeDef);
452 }
453 }
454 G4DNAChemistryManager::Instance()->Initialize();
455 }
456
457 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
458
459 void EmDNAChemistry::ConstructTimeStepModel(G4DNAMolecularReactionTable* reactionTable)
460 {
461 auto reactionRadiusComputer = new G4DNASmoluchowskiReactionModel();
462 reactionTable->PrintTable(reactionRadiusComputer);
463 auto stepByStep = new G4DNAMolecularStepByStepModel();
464 stepByStep->SetReactionModel(reactionRadiusComputer);
465 RegisterTimeStepModel(stepByStep, 0);
466 }
467
468 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
469