UHDR/src/EmDNAChemistry.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 #include "EmDNAChemistry.hh"
0027
0028 #include "DetectorConstruction.hh"
0029
0030 #include "G4DNAChemistryManager.hh"
0031 #include "G4DNAWaterDissociationDisplacer.hh"
0032 #include "G4ProcessManager.hh"
0033 #include "G4RunManager.hh"
0034 #include "G4SystemOfUnits.hh"
0035
0036 // *** Processes and models for Geant4-DNA
0037
0038 #include "BoundedBrownianAction.hh"
0039
0040 #include "G4ChemReboundTransportation.hh"
0041 #include "G4DNAElectronHoleRecombination.hh"
0042 #include "G4DNAElectronSolvation.hh"
0043 #include "G4DNAMolecularDissociation.hh"
0044 #include "G4DNAMolecularReactionTable.hh"
0045 #include "G4DNAMolecularStepByStepModel.hh"
0046 #include "G4DNASancheExcitationModel.hh"
0047 #include "G4DNASmoluchowskiReactionModel.hh"
0048 #include "G4DNAVibExcitation.hh"
0049 // particles
0050
0051 #include "G4Electron.hh"
0052 #include "G4Electron_aq.hh"
0053 #include "G4H2O.hh"
0054 #include "G4H2O2.hh"
0055 #include "G4H3O.hh"
0056 #include "G4HO2.hh"
0057 #include "G4Hydrogen.hh"
0058 #include "G4MoleculeTable.hh"
0059 #include "G4O2.hh"
0060 #include "G4O3.hh"
0061 #include "G4OH.hh"
0062 #include "G4Oxygen.hh"
0063 #include "G4PhysicsListHelper.hh"
0064 /****/
0065 #include "G4DNAMoleculeEncounterStepper.hh"
0066 #include "G4DNAScavengerProcess.hh"
0067 #include "G4MolecularConfiguration.hh"
0068 #include "G4ProcessTable.hh"
0069 #include "G4VChemistryWorld.hh"
0070 /****/
0071 #include "G4ChemicalMoleculeFinder.hh"
0072 // factory
0073 #include "ChemNO2_NO3ScavengerBuilder.hh"
0074 #include "ChemOxygenWaterBuilder.hh"
0075 #include "ChemPureWaterBuilder.hh"
0076
0077 #include "G4ChemDissociationChannels_option1.hh"
0078 #include "G4DNAIndependentReactionTimeModel.hh"
0079 #include "G4GenericMessenger.hh"
0080 #include "G4PhysicsConstructorFactory.hh"
0081 G4_DECLARE_PHYSCONSTR_FACTORY(EmDNAChemistry);
0082
0083 EmDNAChemistry::EmDNAChemistry() : G4VUserChemistryList(true)
0084 {
0085   G4DNAChemistryManager::Instance()->SetChemistryList(this);
0086   //  DefineCommands(); // Le Tuan Anh: create cmds
0087 }
0088
0089 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0090
0091 EmDNAChemistry::~EmDNAChemistry() = default;
0092
0093 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0094
0095 void EmDNAChemistry::ConstructMolecule()
0096 {
0097   G4ChemDissociationChannels_option1::ConstructMolecule();
0098   auto table = G4MoleculeTable::Instance();
0099
0100   auto H3OpB = table->GetConfiguration("H3Op(B)");
0101   H3OpB->SetDiffusionCoefficient(9.46e-9 * (m2 / s));
0102
0103   auto OHm = table->GetConfiguration("OHm(B)");
0104   OHm->SetDiffusionCoefficient(5.3e-9 * (m2 / s));
0105   table->CreateConfiguration("H2O", G4H2O::Definition());
0106
0107   auto G4NO2 = new G4MoleculeDefinition("NO_2", /*mass*/ 30,
0108                                         /*D*/ 5.3e-9 * (m * m / s),
0109                                         /*charge*/ 0,
0110                                         /*electronL*/ 0,
0111                                         /*radius*/ 0.17 * nm);  // should be corrected
0112
0113   auto G4NO3 = new G4MoleculeDefinition("NO_3", /*mass*/ 38,
0114                                         /*D*/ 0 * (m * m / s),
0115                                         /*charge*/ 0,
0116                                         /*electronL*/ 0,
0117                                         /*radius*/ 0.17 * nm);  // should be corrected
0118
0119   table->CreateConfiguration("NO2", G4NO2);
0120   table->CreateConfiguration("NO2m", G4NO2,
0121                              -1,  // charge
0122                              1.9e-9 * (m2 / s));
0123   table->CreateConfiguration("NO2mm", G4NO2,
0124                              -2,  // charge
0125                              1.9e-9 * (m2 / s));
0126
0127   table->CreateConfiguration("NO3m", G4NO3,
0128                              -1,  // charge
0129                              1.9e-9 * (m2 / s));
0130
0131   table->CreateConfiguration("NO3mm", G4NO3,
0132                              -2,  // charge
0133                              1.9e-9 * (m2 / s));
0134
0135   // FrickeDosimeter
0136   auto G4Fe = new G4MoleculeDefinition("Fe",
0137                                        /*mass*/ 55.84 * g / Avogadro * c_squared,
0138                                        /*D*/ 0 * (m * m / s),
0139                                        /*charge*/ 0,
0140                                        /*electronL*/ 0,
0141                                        /*radius*/ 0.35 * nm);  // can be adjusted
0142
0143   table->CreateConfiguration("Fe0", G4Fe);
0144
0145   table->CreateConfiguration("Feppp", G4Fe,
0146                              3,  // charge
0147                              4.86e-10 * (m2 / s));  // Michael Spiro* and Andrew M. Creeth
0148
0149   table->CreateConfiguration("Fepp", G4Fe,
0150                              2,  // charge
0151                              5.78e-10 * (m2 / s));
0152   // HSO4-
0153   auto G4HSO4 = new G4MoleculeDefinition("HSO4",
0154                                          /*mass*/ 55.84 * g / Avogadro * c_squared,
0155                                          /*D*/ 0 * (m * m / s),
0156                                          /*charge*/ 0,
0157                                          /*electronL*/ 0,
0158                                          /*radius*/ 0.35 * nm);  // can be adjusted
0159   table->CreateConfiguration("HSO4m", G4HSO4,
0160                              -1,  // charge
0161                              0 * (m2 / s));
0162
0163   // SO4-
0164   auto G4SO4 = new G4MoleculeDefinition("SO4",
0165                                         /*mass*/ 55.84 * g / Avogadro * c_squared,
0166                                         /*D*/ 0 * (m * m / s),
0167                                         /*charge*/ 0,
0168                                         /*electronL*/ 0,
0169                                         /*radius*/ 0.35 * nm);  // can be adjusted
0170   table->CreateConfiguration("SO4m", G4SO4,
0171                              -1,  // charge
0172                              0 * (m2 / s));
0173
0174   // CO2
0175   auto G4CO2 = new G4MoleculeDefinition("CO_2",
0176                                         /*mass*/ 44.01 * g / Avogadro * c_squared,
0177                                         /*D*/ 1.88e-9 * (m * m / s),
0178                                         /*charge*/ 0,
0179                                         /*electronL*/ 0,
0180                                         /*radius*/ 0.35 * nm);  // can be adjusted
0181   table->CreateConfiguration("CO2", G4CO2,
0182                              0,  // charge
0183                              1.88e-9 * (m2 / s));
0184
0185   table->CreateConfiguration("CO2m", G4CO2,
0186                              -1,  // charge
0187                              1.88e-9 * (m2 / s));
0188
0189   // HCO3-
0190   auto G4HCO3 = new G4MoleculeDefinition("HCO_3",
0191                                          /*mass*/ 61.01 * g / Avogadro * c_squared,
0192                                          /*D*/ 1.88e-9 * (m * m / s),
0193                                          /*charge*/ 0,
0194                                          /*electronL*/ 0,
0195                                          /*radius*/ 0.35 * nm);  // can be adjusted
0196
0197   table->CreateConfiguration("HCO3", G4HCO3,
0198                              0,  // charge
0199                              1.88e-9 * (m2 / s));
0200   table->CreateConfiguration("HCO3m", G4HCO3,
0201                              -1,  // charge
0202                              1.88e-9 * (m2 / s));
0203
0204   // CO3-
0205   auto G4CO3 = new G4MoleculeDefinition("CO_3",
0206                                         /*mass*/ 61.01 * g / Avogadro * c_squared,
0207                                         /*D*/ 0.8e-9 * (m * m / s),
0208                                         /*charge*/ 0,
0209                                         /*electronL*/ 0,
0210                                         /*radius*/ 0.35 * nm);  // can be adjusted
0211   table->CreateConfiguration("CO3m", G4CO3,
0212                              -1,  // charge
0213                              0.8e-9 * (m2 / s));
0214
0215   table->CreateConfiguration("CO3mm", G4CO3,
0216                              -2,  // charge
0217                              0.8e-9 * (m2 / s));
0218
0219   // G4N2O
0220   auto G4N2O = new G4MoleculeDefinition("N_2O",
0221                                         /*mass*/ 61.01 * g / Avogadro * c_squared,
0222                                         /*D*/ 0.8e-9 * (m * m / s),  // not corrected
0223                                         /*charge*/ 0,
0224                                         /*electronL*/ 0,
0225                                         /*radius*/ 0.35 * nm);  // can be adjusted
0226   table->CreateConfiguration("N2O", G4N2O,
0227                              0,  // charge
0228                              0.8e-9 * (m2 / s));  // not corrected
0229
0230   // MeOH
0231   auto G4MeOH = new G4MoleculeDefinition("MeOH",
0232                                          /*mass*/ 61.01 * g / Avogadro * c_squared,
0233                                          /*D*/ 1e-10 * (m * m / s),  // not corrected
0234                                          /*charge*/ 0,
0235                                          /*electronL*/ 0,
0236                                          /*radius*/ 0.35 * nm);  // can be adjusted
0237   table->CreateConfiguration("CH3OH", G4MeOH,
0238                              0,  // charge
0239                              1e-10 * (m2 / s));  // n
0240   table->CreateConfiguration("CH2OH", G4MeOH,
0241                              0,  // charge
0242                              1e-10 * (m2 / s));  // n
0243   // CH2OH
0244 }
0245
0246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0247
0248 void EmDNAChemistry::ConstructDissociationChannels()
0249 {
0250   G4ChemDissociationChannels_option1::ConstructDissociationChannels();
0251 }
0252
0253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0254
0255 void EmDNAChemistry::ConstructReactionTable(G4DNAMolecularReactionTable* pReactionTable)
0256 {
0257   ChemOxygenWaterBuilder::OxygenScavengerReaction(pReactionTable);
0258   ChemOxygenWaterBuilder::CO2ScavengerReaction(pReactionTable);
0259   ChemOxygenWaterBuilder::SecondOrderReactionExtended(pReactionTable);
0260   ChemPureWaterBuilder::WaterScavengerReaction(pReactionTable);
0261   ChemNO2_NO3ScavengerBuilder::NO2_NO3ScavengerReaction(pReactionTable);
0262 }
0263
0264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0265
0266 void EmDNAChemistry::ConstructProcess()
0267 {
0268   auto table = G4MoleculeTable::Instance();
0269   auto O2 = table->GetConfiguration("O2");
0270   auto O2m = table->GetConfiguration("O2m");
0271   auto HO2 = table->GetConfiguration("HO2°");
0272
0273   auto e_aq = table->GetConfiguration("e_aq");
0274   auto OH = table->GetConfiguration("°OH");
0275   auto OHm = table->GetConfiguration("OHm");
0276
0277   auto NO2 = table->GetConfiguration("NO2");
0278   auto NO2m = table->GetConfiguration("NO2m");
0279   auto NO2mm = table->GetConfiguration("NO2mm");
0280   auto NO3m = table->GetConfiguration("NO3m");
0281   auto NO3mm = table->GetConfiguration("NO3mm");
0282
0283   auto HCO3m = table->GetConfiguration("HCO3m");
0284   auto* CO2 = table->GetConfiguration("CO2");
0285   auto* CO2m = table->GetConfiguration("CO2m");
0286   auto H2O2 = table->GetConfiguration("H2O2");
0287   auto H = table->GetConfiguration("H");
0288   auto* HCO3 = table->GetConfiguration("HCO3");
0289   auto* H3OpB = table->GetConfiguration("H3Op(B)");
0290   auto* OHmB = table->GetConfiguration("OHm(B)");
0291   auto* HO2m = table->GetConfiguration("HO2m");
0292   auto* Om = table->GetConfiguration("Om");
0293   auto* O3m = table->GetConfiguration("O3m");
0294   auto* H3Op = table->GetConfiguration("H3Op");
0295   auto H2O = table->GetConfiguration("H2O");
0296   auto* N2O = table->GetConfiguration("N2O");
0297   auto* MeOH = table->GetConfiguration("CH3OH");
0298   auto* CH2OH = table->GetConfiguration("CH2OH");
0299   //auto* CO3m = table->GetConfiguration("CO3m");
0300
0301   const auto* fpDetector = dynamic_cast<const DetectorConstruction*>(
0302     G4RunManager::GetRunManager()->GetUserDetectorConstruction());
0303   auto fpChemistryWorld = fpDetector->GetChemistryWorld();
0304   fpChemistryWorld->ConstructChemistryComponents();
0305   auto confinedBox = fpChemistryWorld->GetChemistryBoundary();
0306
0307   auto* ph = G4PhysicsListHelper::GetPhysicsListHelper();
0308
0309   //===============================================================
0310   // Extend vibrational to low energy
0311   // Anyway, solvation of electrons is taken into account from 7.4 eV
0312   // So below this threshold, for now, no accurate modeling is done
0313   //
0314   G4VProcess* process =
0315     G4ProcessTable::GetProcessTable()->FindProcess("e-_G4DNAVibExcitation", "e-");
0316
0317   if (process) {
0318     auto vibExcitation = (G4DNAVibExcitation*)process;
0319     G4VEmModel* model = vibExcitation->EmModel();
0320     auto sancheExcitationMod = dynamic_cast<G4DNASancheExcitationModel*>(model);
0321     if (sancheExcitationMod) {
0322       sancheExcitationMod->ExtendLowEnergyLimit(0.025 * eV);
0323     }
0324   }
0325
0326   //===============================================================
0327   // *** Electron Solvatation ***
0328   //
0329   process = G4ProcessTable::GetProcessTable()->FindProcess("e-_G4DNAElectronSolvation", "e-");
0330
0331   if (process == nullptr) {
0332     ph->RegisterProcess(new G4DNAElectronSolvation("e-_G4DNAElectronSolvation"),
0333                         G4Electron::Definition());
0334   }
0335
0336   //===============================================================
0337   // Define processes for molecules
0338   //
0339   auto* theMoleculeTable = G4MoleculeTable::Instance();
0340   auto iterator = theMoleculeTable->GetDefintionIterator();
0341   iterator.reset();
0342
0343   while (iterator()) {
0344     auto* moleculeDef = iterator.value();
0345
0346     if (moleculeDef != G4H2O::Definition()) {
0347       auto brown = new G4ChemReboundTransportation("ReboundTransport");
0348       brown->SetBoundary(confinedBox);
0349       ph->RegisterProcess(brown, moleculeDef);
0350     }
0351     else {
0352       moleculeDef->GetProcessManager()->AddRestProcess(new G4DNAElectronHoleRecombination(), 2);
0353       auto brownTransport = new BoundedBrownianAction();
0354       brownTransport->SetBoundary(*confinedBox);
0355       auto dissociationProcess = new G4DNAMolecularDissociation("H2O_DNAMolecularDecay", fDecay);
0356       dissociationProcess->SetUserBrownianAction(brownTransport);
0357       dissociationProcess->SetDisplacer(moleculeDef, new G4DNAWaterDissociationDisplacer);
0358       // dissociationProcess->SetVerbose(1);
0359       moleculeDef->GetProcessManager()->AddRestProcess(dissociationProcess, 1);
0360     }
0361
0362     if (moleculeDef == G4Hydrogen::Definition()) {
0363       // O2
0364       auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
0365       //------------------------------------------------------------------
0366       // H + O2(B) -> HO2
0367       auto reactionData = new G4DNAMolecularReactionData(2.1e10 * (1e-3 * m3 / (mole * s)), H, O2);
0368       reactionData->AddProduct(HO2);
0369       scanvergerProcess->SetReaction(H, reactionData);
0370       //------------------------------------------------------------------
0371       // H + OH-(B) -> H2O + eaq- 2.49e3 / s
0372       reactionData = new G4DNAMolecularReactionData(2.49e7 * (1e-3 * m3 / (mole * s)), H,
0373                                                     OHmB);  // 2.51e7 (H + OH-)* 1e-7 (pH) = 2.48e0
0374       reactionData->AddProduct(e_aq);
0375       scanvergerProcess->SetReaction(H, reactionData);
0376       //------------------------------------------------------------------
0377       // H + H2O -> eaq- + H3O+ 5.94 / s pkA = 9.5515  //2
0378       reactionData = new G4DNAMolecularReactionData(6.32 / s, H, H2O);  // 6.32e0 *
0379       reactionData->AddProduct(e_aq);
0380       reactionData->AddProduct(H3OpB);
0381       scanvergerProcess->SetReaction(H, reactionData);
0382       // H2O2
0383       //------------------------------------------------------------------
0384       // H + H202 -> OH + H20
0385       reactionData = new G4DNAMolecularReactionData(9.0e7 * (1e-3 * m3 / (mole * s)), H, H2O2);
0386       reactionData->AddProduct(OH);
0387       scanvergerProcess->SetReaction(H, reactionData);
0388
0389       ph->RegisterProcess(scanvergerProcess, moleculeDef);
0390     }
0391     if (moleculeDef == G4Electron_aq::Definition()) {
0392       auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
0393       G4DNAMolecularReactionData* reactionData = nullptr;
0394       //------------------------------------------------------------------
0395       // e_aq + O2(B) -> O2-
0396       reactionData = new G4DNAMolecularReactionData(1.74e10 * (1e-3 * m3 / (mole * s)), e_aq, O2);
0397       reactionData->AddProduct(O2m);
0398       scanvergerProcess->SetReaction(e_aq, reactionData);
0399       //------------------------------------------------------------------
0400       // eaq- + H3O+(B) -> H + H2O 2.09e3 / s  //2
0401       reactionData =
0402         new G4DNAMolecularReactionData(2.11e10 * (1e-3 * m3 / (mole * s)), e_aq,
0403                                        H3OpB);  // 2.11e10 (e_aq + H3O+) * 1.0e-7 (Ph=7) = 2.09e3
0404       reactionData->AddProduct(H);
0405       scanvergerProcess->SetReaction(e_aq, reactionData);
0406       //------------------------------------------------------------------
0407       // e_aq + NO2- -> NO2--
0408       reactionData = new G4DNAMolecularReactionData(3.5e9 * (1e-3 * m3 / (mole * s)), e_aq, NO2m);
0409       reactionData->AddProduct(NO2mm);
0410       scanvergerProcess->SetReaction(e_aq, reactionData);
0411       //------------------------------------------------------------------
0412       // e_aq + NO3- -> NO3--
0413       reactionData = new G4DNAMolecularReactionData(9.7e9 * (1e-3 * m3 / (mole * s)), e_aq, NO3m);
0414       reactionData->AddProduct(NO3mm);
0415       scanvergerProcess->SetReaction(e_aq, reactionData);
0416       //------------------------------------------------------------------                                                                                                      // eaq- + H2O -> H + OH- 15.7 / M * s pKa = ???  //3
0417                                                                                                       reactionData =
0418         new G4DNAMolecularReactionData(1.57e1 * 55.3 / s, e_aq, H2O);
0419       reactionData->AddProduct(H);
0420       reactionData->AddProduct(OHmB);
0421       scanvergerProcess->SetReaction(e_aq, reactionData);
0422
0423       //------------------------------------------------------------------
0424       // Oxygen concentration
0425       // e_aq + CO2(B) -> CO2- k = 0.77 × 1010
0426       reactionData = new G4DNAMolecularReactionData(0.77e10 * (1e-3 * m3 / (mole * s)), e_aq, CO2);
0427       reactionData->AddProduct(CO2m);
0428       scanvergerProcess->SetReaction(e_aq, reactionData);
0429       // scanvergerProcess->SetVerboseLevel(1);
0430
0431       //------------------------------------------------------------------
0432       reactionData = new G4DNAMolecularReactionData(0.9e10 * (1e-3 * m3 / (mole * s)), e_aq, N2O);
0433       reactionData->AddProduct(Om);
0434       scanvergerProcess->SetReaction(e_aq, reactionData);
0435
0436       //------------------------------------------------------------------
0437       ph->RegisterProcess(scanvergerProcess, moleculeDef);
0438     }
0439     if (moleculeDef == G4O2::Definition()) {
0440       auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
0441       G4DNAMolecularReactionData* reactionData = nullptr;
0442       //------------------------------------------------------------------
0443       // O2- + H3O+(B) -> HO2 + H2O 4.73e3 / s  //1
0444       reactionData =
0445         new G4DNAMolecularReactionData(4.78e10 * (1e-3 * m3 / (mole * s)), O2m,
0446                                        H3OpB);  // 4.78e10(O2- + H3O+) * 1e-7(pH7) = 4.73e3
0447       reactionData->AddProduct(HO2);
0448       reactionData->AddProduct(H2O);
0449       reactionData->SetReactionType(6);  // Equilibrium 6
0450       scanvergerProcess->SetReaction(O2m, reactionData);
0451       //------------------------------------------------------------------
0452       // O2- + H2O -> HO2 + OH- 0.15 / s  //4
0453       reactionData = new G4DNAMolecularReactionData(0.15 * 55.3 / s, O2m, H2O);
0454       reactionData->AddProduct(HO2);
0455       reactionData->AddProduct(OHmB);
0456       scanvergerProcess->SetReaction(O2m, reactionData);
0457
0458       //------------------------------------------------------------------
0459       ph->RegisterProcess(scanvergerProcess, moleculeDef);
0460     }
0461     if (moleculeDef == G4OH::Definition()) {
0462       auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
0463       // scanvergerProcess->SetVerboseLevel(1);
0464       G4DNAMolecularReactionData* reactionData = nullptr;
0465       //------------------------------------------------------------------
0466       // OH + OH-(B) -> O- + H2O 6.24e2 / s  //6
0467       reactionData =
0468         new G4DNAMolecularReactionData(1.27e10 * (1e-3 * m3 / (mole * s)), OH,
0469                                        OHmB);  // 6.30e9 (OH + OH-) * 1e-7 (pH) = 6.24e2
0470       reactionData->AddProduct(Om);
0471       reactionData->AddProduct(H2O);
0472       reactionData->SetReactionType(8);  // Equilibrium 8
0473       scanvergerProcess->SetReaction(OH, reactionData);
0474       //------------------------------------------------------------------
0475       // OH + NO2- -> NO2 + OH-
0476       reactionData = new G4DNAMolecularReactionData(8e9 * (1e-3 * m3 / (mole * s)), OH, NO2m);
0477       reactionData->AddProduct(NO2);
0478       reactionData->AddProduct(OHm);
0479       scanvergerProcess->SetReaction(OH, reactionData);
0480       //------------------------------------------------------------------
0481       // OH + HCO3m -> NO2 + OH-
0482       reactionData = new G4DNAMolecularReactionData(8.5e6 * (1e-3 * m3 / (mole * s)), OH, HCO3m);
0483       scanvergerProcess->SetReaction(OH, reactionData);
0484       // scanvergerProcess->SetVerboseLevel(1);
0485       //------------------------------------------------------------------
0486       //  OH -> O- + H3O+(B)  //8
0487       reactionData = new G4DNAMolecularReactionData(0.060176635 / s, OH,
0488                                                     H2O);  //
0489       reactionData->AddProduct(Om);
0490       reactionData->AddProduct(H3OpB);
0491       scanvergerProcess->SetReaction(OH, reactionData);
0492       //------------------------------------------------------------------
0493       // OH + CO2(B) -> HCO3
0494       reactionData = new G4DNAMolecularReactionData(1.e7 * (1e-3 * m3 / (mole * s)), OH, CO2);
0495       reactionData->AddProduct(HCO3);
0496       scanvergerProcess->SetReaction(OH, reactionData);
0497       // scanvergerProcess->SetVerboseLevel(1);
0498       //------------------------------------------------------------------
0499       //  OH + CH3OH -> CH2OH + H2O
0500       reactionData = new G4DNAMolecularReactionData(9.7e8 * (1e-3 * m3 / (mole * s)), OH, MeOH);
0501       reactionData->AddProduct(CH2OH);
0502       scanvergerProcess->SetReaction(OH, reactionData);
0503       // scanvergerProcess->SetVerboseLevel(1);
0504       //------------------------------------------------------------------
0505       ph->RegisterProcess(scanvergerProcess, moleculeDef);
0506     }
0507     if (moleculeDef == G4MoleculeTable::Instance()->GetMoleculeDefinition("OH"))  // OH-
0508     {
0509       auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
0510       G4DNAMolecularReactionData* reactionData = nullptr;
0511       //------------------------------------------------------------------
0512       // OH- + H3O+(B) -> 2H2O 1.11e4 / s
0513       reactionData =
0514         new G4DNAMolecularReactionData(1.13e11 * (1e-3 * m3 / (mole * s)), OHm,
0515                                        H3OpB);  // 1.13e11 (H3O+ + OH-) * 1e-7 (pH=7) =1.12e4
0516       scanvergerProcess->SetReaction(OHm, reactionData);
0517       //------------------------------------------------------------------
0518       ph->RegisterProcess(scanvergerProcess, moleculeDef);
0519     }
0520     if (moleculeDef == G4HO2::Definition()) {
0521       auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
0522       // scanvergerProcess->SetVerboseLevel(1);
0523       G4DNAMolecularReactionData* reactionData = nullptr;
0524       //------------------------------------------------------------------
0525       // HO2 + OH-(B) -> O2- + H2O 6.24e2 / s //4
0526       reactionData = new G4DNAMolecularReactionData(1.27e10 * (1e-3 * m3 / (mole * s)), HO2,
0527                                                     OHmB);  // 6.30e9(HO2 + OH-)*1e-7 (pH) = 6.24e2
0528       reactionData->AddProduct(O2m);
0529       scanvergerProcess->SetReaction(HO2, reactionData);
0530       //------------------------------------------------------------------
0531       // HO2 + H2O -> H3O+ + O2- //1
0532       reactionData = new G4DNAMolecularReactionData(7.58e5 / s, HO2, H2O);
0533       reactionData->AddProduct(H3OpB);
0534       reactionData->AddProduct(O2m);
0535       reactionData->SetReactionType(6);  // Equilibrium 6
0536       scanvergerProcess->SetReaction(HO2, reactionData);
0537       //------------------------------------------------------------------
0538       ph->RegisterProcess(scanvergerProcess, moleculeDef);
0539     }
0540
0541     if (moleculeDef == G4MoleculeTable::Instance()->GetMoleculeDefinition("HO_2")) /*HO2-*/ {
0542       auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
0543       // scanvergerProcess->SetVerboseLevel(1);
0544       G4DNAMolecularReactionData* reactionData = nullptr;
0545       //------------------------------------------------------------------
0546       // HO2- + H3O+(B) -> H2O2 + H2O 4.98e3 / s  //7
0547       reactionData =
0548         new G4DNAMolecularReactionData(4.78e10 * (1e-3 * m3 / (mole * s)), HO2m,
0549                                        H3OpB);  // 5.00e10 (H3O+ + HO2-) * 1e-7(pH) = 4.95e3
0550       reactionData->AddProduct(H2O2);
0551       scanvergerProcess->SetReaction(HO2m, reactionData);
0552       //------------------------------------------------------------------
0553       // HO2- + H2O -> H2O2 + OH- 1.36e6 / M * s pka = 11.784  //5
0554       reactionData = new G4DNAMolecularReactionData(1.36e6 * 55.3 / s, HO2m, H2O);  //
0555       reactionData->AddProduct(H2O2);
0556       reactionData->AddProduct(OHmB);
0557       reactionData->SetReactionType(7);  // Equilibrium 7
0558       scanvergerProcess->SetReaction(HO2m, reactionData);
0559       //------------------------------------------------------------------
0560       ph->RegisterProcess(scanvergerProcess, moleculeDef);
0561     }
0562
0563     if (moleculeDef == G4Oxygen::Definition()) {
0564       auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
0565       G4DNAMolecularReactionData* reactionData = nullptr;
0566       //------------------------------------------------------------------
0567       // O- + H3O+(B) -> OH + H2O 4.73e3 / s //8
0568       reactionData =
0569         new G4DNAMolecularReactionData(9.56e10 * (1e-3 * m3 / (mole * s)), Om,
0570                                        H3OpB);  // 4.78e10 (H3O+ + O2-) * 1e-7(pH) = 4.73e3
0571       reactionData->AddProduct(OH);
0572       scanvergerProcess->SetReaction(Om, reactionData);
0573       //------------------------------------------------------------------
0574       // O- + H2O -> OH + OH- 1.8e6 / s pka = 11.9  //6
0575       reactionData = new G4DNAMolecularReactionData(1.8e6 * 55.3 / s, Om, H2O);
0576       reactionData->AddProduct(OH);
0577       reactionData->AddProduct(OHmB);
0578       reactionData->SetReactionType(8);  // Equilibrium 8
0579       scanvergerProcess->SetReaction(Om, reactionData);
0580       //------------------------------------------------------------------
0581       // O- + O2(B) -> O3-
0582       reactionData = new G4DNAMolecularReactionData(3.7e9 * (1e-3 * m3 / (mole * s)), Om, O2);
0583       reactionData->AddProduct(O3m);
0584       scanvergerProcess->SetReaction(Om, reactionData);
0585       //------------------------------------------------------------------
0586       ph->RegisterProcess(scanvergerProcess, moleculeDef);
0587     }
0588     if (moleculeDef == G4O3::Definition()) {
0589       auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
0590       G4DNAMolecularReactionData* reactionData = nullptr;
0591       //------------------------------------------------------------------
0592       // O3- + H3O+(B) -> OH + O2 + H2O 8.91e3 / s
0593       reactionData =
0594         new G4DNAMolecularReactionData(9.0e10 * (1e-3 * m3 / (mole * s)), O3m,
0595                                        H3OpB);  // 9.0e10 (O3- + H3O+) * 1e-7(pH) = 8.91e3
0596       reactionData->AddProduct(OH);
0597       reactionData->AddProduct(O2);
0598       scanvergerProcess->SetReaction(O3m, reactionData);
0599       //------------------------------------------------------------------
0600       // O3- + H2OB -> O- + O2
0601       reactionData = new G4DNAMolecularReactionData(2.66e3 / s, O3m, H2O);
0602       reactionData->AddProduct(Om);
0603       reactionData->AddProduct(O2);
0604       scanvergerProcess->SetReaction(O3m, reactionData);
0605       //------------------------------------------------------------------
0606       ph->RegisterProcess(scanvergerProcess, moleculeDef);
0607     }
0608     if (moleculeDef == G4H3O::Definition()) {
0609       auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
0610       G4DNAMolecularReactionData* reactionData = nullptr;
0611       //------------------------------------------------------------------
0612       // H3O+ + OH-(B) -> 2H2O 1.11e4 / s
0613       reactionData =
0614         new G4DNAMolecularReactionData(1.13e11 * (1e-3 * m3 / (mole * s)), H3Op,
0615                                        OHmB);  // 1.13e11 (H3O+ + OH-) * 1e-7 (pH=7) = 1.12e4
0616       scanvergerProcess->SetReaction(H3Op, reactionData);
0617       //------------------------------------------------------------------
0618       ph->RegisterProcess(scanvergerProcess, moleculeDef);
0619     }
0620     if (moleculeDef == G4H2O2::Definition()) {
0621       auto scanvergerProcess = new G4DNAScavengerProcess("G4DNAScavengerProcess", *confinedBox);
0622       // scanvergerProcess->SetVerboseLevel(1);
0623       G4DNAMolecularReactionData* reactionData = nullptr;
0624       //------------------------------------------------------------------
0625       // H2O2 + OH-(B) -> HO2- + H2O 4.66e2 / s //5
0626       reactionData =
0627         new G4DNAMolecularReactionData(1.27e10 * (1e-3 * m3 / (mole * s)), H2O2,
0628                                        OHmB);  // 4.71e8 (H2O2 + OH-) * 1e-7 (pH) = 4.66e1
0629       reactionData->AddProduct(HO2m);
0630       reactionData->SetReactionType(7);  // Equilibrium 7
0631       scanvergerProcess->SetReaction(H2O2, reactionData);
0632       //------------------------------------------------------------------
0633       // H2O2 + H2O -> H+ + HO2- First order pka = 11.784  //7
0634       reactionData = new G4DNAMolecularReactionData(7.86e-2 / s, H2O2, H2O);
0635       reactionData->AddProduct(HO2m);
0636       reactionData->AddProduct(H3OpB);
0637       scanvergerProcess->SetReaction(H2O2, reactionData);
0638       //------------------------------------------------------------------
0639
0640       ph->RegisterProcess(scanvergerProcess, moleculeDef);
0641     }
0642   }
0643   G4DNAChemistryManager::Instance()->Initialize();
0644 }
0645
0646 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0647
0648 void EmDNAChemistry::ConstructTimeStepModel(G4DNAMolecularReactionTable* reactionTable)
0649 {
0650   // Le Tuan Anh: swtich between time-step-models
0651   G4EmParameters* param = G4EmParameters::Instance();
0652   auto model = param->GetTimeStepModel();
0653   if (model == G4ChemTimeStepModel::SBS) {
0654     auto reactionRadiusComputer = new G4DNASmoluchowskiReactionModel();
0655     reactionTable->PrintTable(reactionRadiusComputer);
0656     auto stepByStep = new G4DNAMolecularStepByStepModel();
0657     stepByStep->SetReactionModel(reactionRadiusComputer);
0658     RegisterTimeStepModel(stepByStep, 0);
0659   }
0660   else if (model == G4ChemTimeStepModel::IRT_syn) {
0661     RegisterTimeStepModel(new G4DNAIndependentReactionTimeModel(), 0);
0662   }
0663 }
0664
0665 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....