scavenger/src/EmDNAChemistry.cc

0001 // ********************************************************************
0002 // * License and Disclaimer                                           *
0003 // *                                                                  *
0004 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0005 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0006 // * conditions of the Geant4 Software License,  included in the file *
0007 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0008 // * include a list of copyright holders.                             *
0009 // *                                                                  *
0010 // * Neither the authors of this software system, nor their employing *
0011 // * institutes,nor the agencies providing financial support for this *
0012 // * work  make  any representation or  warranty, express or implied, *
0013 // * regarding  this  software system or assume any liability for its *
0014 // * use.  Please see the license in the file  LICENSE  and URL above *
0015 // * for the full disclaimer and the limitation of liability.         *
0016 // *                                                                  *
0017 // * This  code  implementation is the result of  the  scientific and *
0018 // * technical work of the GEANT4 collaboration.                      *
0019 // * By using,  copying,  modifying or  distributing the software (or *
0020 // * any work based  on the software)  you  agree  to acknowledge its *
0021 // * use  in  resulting  scientific  publications,  and indicate your *
0022 // * acceptance of all terms of the Geant4 Software license.          *
0023 // ********************************************************************
0024 //
0025 /// \file scavenger/src/EmDNAChemistry.cc
0026 /// \brief Implementation of the scavenger::EmDNAChemistry class
0027
0028 #include "EmDNAChemistry.hh"
0029
0030 #include "G4DNAChemistryManager.hh"
0031 #include "G4DNAElectronHoleRecombination.hh"
0032 #include "G4DNAElectronSolvation.hh"
0033 #include "G4DNAIRT.hh"
0034 #include "G4DNAMolecularDissociation.hh"
0035 #include "G4DNAMolecularIRTModel.hh"
0036 #include "G4DNAMolecularReactionTable.hh"
0037 #include "G4DNASancheExcitationModel.hh"
0038 #include "G4DNAVibExcitation.hh"
0039 #include "G4DNAWaterDissociationDisplacer.hh"
0040 #include "G4DNAWaterExcitationStructure.hh"
0041 #include "G4PhysicalConstants.hh"
0042 #include "G4ProcessManager.hh"
0043 #include "G4SystemOfUnits.hh"
0044 #include "G4VDNAReactionModel.hh"
0045 // particles
0046 #include "G4Electron.hh"
0047 #include "G4Electron_aq.hh"
0048 #include "G4H2O.hh"
0049 #include "G4MoleculeTable.hh"
0050 #include "G4O2.hh"
0051 #include "G4O3.hh"
0052 #include "G4OH.hh"
0053 #include "G4Oxygen.hh"
0054 #include "G4PhysicsListHelper.hh"
0055 /****/
0056 #include "G4MolecularConfiguration.hh"
0057 #include "G4ProcessTable.hh"
0058 /****/
0059 // factory
0060 #include "G4PhysicsConstructorFactory.hh"
0061 #include "G4ChemDissociationChannels_option1.hh"
0062 #include "G4ChemDissociationChannels_option1.hh"
0063 namespace scavenger
0064 {
0065
0066 G4_DECLARE_PHYSCONSTR_FACTORY(EmDNAChemistry);
0067
0068 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0069
0070 EmDNAChemistry::EmDNAChemistry()
0071   : G4VUserChemistryList(true),
0072     G4UImessenger(),
0073     fpParserDir(new G4UIdirectory("/chem/reactionTable/")),
0074     fpReactionTableNameCmd(new G4UIcmdWithAString("/chem/reactionTable/name", this))
0075 {
0076   G4DNAChemistryManager::Instance()->SetChemistryList(this);
0077   fpParserDir->SetGuidance("Chemistry control");
0078   fpReactionTableNameCmd->SetGuidance("Name of the chemical reaction table");
0079 }
0080
0081 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0082
0083 void EmDNAChemistry::SetNewValue(G4UIcommand* command, G4String newValue)
0084 {
0085   if (command == fpReactionTableNameCmd.get()) {
0086     fReactionTableName = newValue;
0087   }
0088 }
0089
0090 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0091
0092 void EmDNAChemistry::ConstructMolecule()
0093 {
0094   // Create the definition
0095   G4ChemDissociationChannels_option1::ConstructMolecule();
0096
0097   auto G4NO2 = new G4MoleculeDefinition("NO_2",
0098                                         /*mass*/ 46.0055 * g / Avogadro * c_squared,
0099                                         /*D*/ 1.7e-9 * (m * m / s),
0100                                         /*charge*/ 0,
0101                                         /*electronL*/ 0,
0102                                         /*radius*/ 0.35 * nm);  // can be adjusted
0103
0104   auto G4NO3 = new G4MoleculeDefinition("NO_3",
0105                                         /*mass*/ 62.0049 * g / Avogadro * c_squared,
0106                                         /*D*/ 1.7e-9 * (m * m / s),
0107                                         /*charge*/ 0,
0108                                         /*electronL*/ 0,
0109                                         /*radius*/ 0.35 * nm);  // can be adjusted
0110   //____________________________________________________________________________
0111   // Note: Parameters Value changed according to Plante Paper
0112
0113   G4MoleculeTable::Instance()->CreateConfiguration("O2(B)",  // just a tag to store and retrieve
0114                                                              // from G4MoleculeTable
0115                                                    G4O2::Definition(),
0116                                                    0,  // charge
0117                                                    0 * (m2 / s));
0118
0119   G4MoleculeTable::Instance()->CreateConfiguration("NO2m(B)",  // just a tag to store and retrieve
0120                                                                // from G4MoleculeTable
0121                                                    G4NO2,
0122                                                    -1,  // charge
0123                                                    0 * (m2 / s));
0124
0125   G4MoleculeTable::Instance()->CreateConfiguration("NO3m(B)",  // just a tag to store and retrieve
0126                                                                // from G4MoleculeTable
0127                                                    G4NO3,
0128                                                    -1,  // charge
0129                                                    0 * (m2 / s));
0130
0131 }
0132
0133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0134
0135 void EmDNAChemistry::ConstructDissociationChannels()
0136 {
0137   //-----------------------------------
0138   // Get the molecular configuration
0139   auto OH = G4MoleculeTable::Instance()->GetConfiguration("°OH");
0140   auto OHm = G4MoleculeTable::Instance()->GetConfiguration("OHm");
0141   auto e_aq = G4MoleculeTable::Instance()->GetConfiguration("e_aq");
0142   auto H2 = G4MoleculeTable::Instance()->GetConfiguration("H2");
0143   auto H3O = G4MoleculeTable::Instance()->GetConfiguration("H3Op");
0144   auto H = G4MoleculeTable::Instance()->GetConfiguration("H");
0145
0146   //-------------------------------------
0147   // Define the decay channels
0148   auto* water = G4H2O::Definition();
0149   G4MolecularDissociationChannel* decCh1;
0150   G4MolecularDissociationChannel* decCh2;
0151
0152   auto occ = new G4ElectronOccupancy(*(water->GetGroundStateElectronOccupancy()));
0153
0154   //////////////////////////////////////////////////////////
0155   //            EXCITATIONS                               //
0156   //////////////////////////////////////////////////////////
0157   G4DNAWaterExcitationStructure waterExcitation;
0158   //--------------------------------------------------------
0159   //---------------Excitation on the fifth layer------------
0160
0161   decCh1 = new G4MolecularDissociationChannel("A^1B_1_Relaxation");
0162   decCh2 = new G4MolecularDissociationChannel("A^1B_1_DissociativeDecay");
0163   // Decay 1 : OH + H
0164   decCh1->SetEnergy(waterExcitation.ExcitationEnergy(0));
0165   decCh1->SetProbability(0.35);
0166   decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::NoDisplacement);
0167
0168   decCh2->AddProduct(OH);
0169   decCh2->AddProduct(H);
0170   decCh2->SetProbability(0.65);
0171   decCh2->SetDisplacementType(G4DNAWaterDissociationDisplacer::A1B1_DissociationDecay);
0172
0173   //  water->AddExcitedState("A^1B_1");
0174   occ->RemoveElectron(4, 1);  // this is the transition form ground state to
0175   occ->AddElectron(5, 1);  // the first unoccupied orbital: A^1B_1
0176
0177   water->NewConfigurationWithElectronOccupancy("A^1B_1", *occ);
0178   water->AddDecayChannel("A^1B_1", decCh1);
0179   water->AddDecayChannel("A^1B_1", decCh2);
0180
0181   //--------------------------------------------------------
0182   //---------------Excitation on the fourth layer-----------
0183   decCh1 = new G4MolecularDissociationChannel("B^1A_1_Relaxation_Channel");
0184   decCh2 = new G4MolecularDissociationChannel("B^1A_1_DissociativeDecay");
0185   auto decCh3 = new G4MolecularDissociationChannel("B^1A_1_AutoIonisation_Channel");
0186
0187   // Decay 1 : energy
0188   decCh1->SetEnergy(waterExcitation.ExcitationEnergy(1));
0189   decCh1->SetProbability(0.3);
0190
0191   // Decay 2 : 2OH + H_2
0192   decCh2->AddProduct(H2);
0193   decCh2->AddProduct(OH);
0194   decCh2->AddProduct(OH);
0195   decCh2->SetProbability(0.15);
0196   decCh2->SetDisplacementType(G4DNAWaterDissociationDisplacer::B1A1_DissociationDecay);
0197
0198   // Decay 3 : OH + H_3Op + e_aq
0199   decCh3->AddProduct(OH);
0200   decCh3->AddProduct(H3O);
0201   decCh3->AddProduct(e_aq);
0202   decCh3->SetProbability(0.55);
0203   decCh3->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
0204
0205   *occ = *(water->GetGroundStateElectronOccupancy());
0206   occ->RemoveElectron(3);  // this is the transition form ground state to
0207   occ->AddElectron(5, 1);  // the first unoccupied orbital: B^1A_1
0208
0209   water->NewConfigurationWithElectronOccupancy("B^1A_1", *occ);
0210   water->AddDecayChannel("B^1A_1", decCh1);
0211   water->AddDecayChannel("B^1A_1", decCh2);
0212   water->AddDecayChannel("B^1A_1", decCh3);
0213
0214   //-------------------------------------------------------
0215   //-------------------Excitation of 3rd layer-----------------
0216   decCh1 = new G4MolecularDissociationChannel("Excitation3rdLayer_AutoIonisation_Channel");
0217   decCh2 = new G4MolecularDissociationChannel("Excitation3rdLayer_Relaxation_Channel");
0218
0219   // Decay channel 1 : : OH + H_3Op + e_aq
0220   decCh1->AddProduct(OH);
0221   decCh1->AddProduct(H3O);
0222   decCh1->AddProduct(e_aq);
0223
0224   decCh1->SetProbability(0.5);
0225   decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
0226
0227   // Decay channel 2 : energy
0228   decCh2->SetEnergy(waterExcitation.ExcitationEnergy(2));
0229   decCh2->SetProbability(0.5);
0230
0231   // Electronic configuration of this decay
0232   *occ = *(water->GetGroundStateElectronOccupancy());
0233   occ->RemoveElectron(2, 1);
0234   occ->AddElectron(5, 1);
0235
0236   // Configure the water molecule
0237   water->NewConfigurationWithElectronOccupancy("Excitation3rdLayer", *occ);
0238   water->AddDecayChannel("Excitation3rdLayer", decCh1);
0239   water->AddDecayChannel("Excitation3rdLayer", decCh2);
0240
0241   //-------------------------------------------------------
0242   //-------------------Excitation of 2nd layer-----------------
0243   decCh1 = new G4MolecularDissociationChannel("Excitation2ndLayer_AutoIonisation_Channel");
0244   decCh2 = new G4MolecularDissociationChannel("Excitation2ndLayer_Relaxation_Channel");
0245
0246   // Decay Channel 1 : : OH + H_3Op + e_aq
0247   decCh1->AddProduct(OH);
0248   decCh1->AddProduct(H3O);
0249   decCh1->AddProduct(e_aq);
0250
0251   decCh1->SetProbability(0.5);
0252   decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
0253
0254   // Decay channel 2 : energy
0255   decCh2->SetEnergy(waterExcitation.ExcitationEnergy(3));
0256   decCh2->SetProbability(0.5);
0257
0258   *occ = *(water->GetGroundStateElectronOccupancy());
0259   occ->RemoveElectron(1, 1);
0260   occ->AddElectron(5, 1);
0261
0262   water->NewConfigurationWithElectronOccupancy("Excitation2ndLayer", *occ);
0263   water->AddDecayChannel("Excitation2ndLayer", decCh1);
0264   water->AddDecayChannel("Excitation2ndLayer", decCh2);
0265
0266   //-------------------------------------------------------
0267   //-------------------Excitation of 1st layer-----------------
0268   decCh1 = new G4MolecularDissociationChannel("Excitation1stLayer_AutoIonisation_Channel");
0269   decCh2 = new G4MolecularDissociationChannel("Excitation1stLayer_Relaxation_Channel");
0270
0271   *occ = *(water->GetGroundStateElectronOccupancy());
0272   occ->RemoveElectron(0, 1);
0273   occ->AddElectron(5, 1);
0274
0275   // Decay Channel 1 : : OH + H_3Op + e_aq
0276   decCh1->AddProduct(OH);
0277   decCh1->AddProduct(H3O);
0278   decCh1->AddProduct(e_aq);
0279   decCh1->SetProbability(0.5);
0280   decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
0281
0282   // Decay channel 2 : energy
0283   decCh2->SetEnergy(waterExcitation.ExcitationEnergy(4));
0284   decCh2->SetProbability(0.5);
0285
0286   water->NewConfigurationWithElectronOccupancy("Excitation1stLayer", *occ);
0287   water->AddDecayChannel("Excitation1stLayer", decCh1);
0288   water->AddDecayChannel("Excitation1stLayer", decCh2);
0289
0290   /////////////////////////////////////////////////////////
0291   //                  IONISATION                         //
0292   /////////////////////////////////////////////////////////
0293   //--------------------------------------------------------
0294   //------------------- Ionisation -------------------------
0295
0296   decCh1 = new G4MolecularDissociationChannel("Ionisation_Channel");
0297
0298   // Decay Channel 1 : : OH + H_3Op
0299   decCh1->AddProduct(H3O);
0300   decCh1->AddProduct(OH);
0301   decCh1->SetProbability(1);
0302   decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::Ionisation_DissociationDecay);
0303
0304   *occ = *(water->GetGroundStateElectronOccupancy());
0305   occ->RemoveElectron(4, 1);
0306   // this is a ionized h2O with a hole in its last orbital
0307   water->NewConfigurationWithElectronOccupancy("Ionisation5", *occ);
0308   water->AddDecayChannel("Ionisation5", decCh1);
0309
0310   *occ = *(water->GetGroundStateElectronOccupancy());
0311   occ->RemoveElectron(3, 1);
0312   water->NewConfigurationWithElectronOccupancy("Ionisation4", *occ);
0313   water->AddDecayChannel("Ionisation4", new G4MolecularDissociationChannel(*decCh1));
0314
0315   *occ = *(water->GetGroundStateElectronOccupancy());
0316   occ->RemoveElectron(2, 1);
0317   water->NewConfigurationWithElectronOccupancy("Ionisation3", *occ);
0318   water->AddDecayChannel("Ionisation3", new G4MolecularDissociationChannel(*decCh1));
0319
0320   *occ = *(water->GetGroundStateElectronOccupancy());
0321   occ->RemoveElectron(1, 1);
0322   water->NewConfigurationWithElectronOccupancy("Ionisation2", *occ);
0323   water->AddDecayChannel("Ionisation2", new G4MolecularDissociationChannel(*decCh1));
0324
0325   *occ = *(water->GetGroundStateElectronOccupancy());
0326   occ->RemoveElectron(0, 1);
0327   water->NewConfigurationWithElectronOccupancy("Ionisation1", *occ);
0328   water->AddDecayChannel("Ionisation1", new G4MolecularDissociationChannel(*decCh1));
0329
0330   //////////////////////////////////////////////////////////
0331   //            Dissociative Attachment                   //
0332   //////////////////////////////////////////////////////////
0333   decCh1 = new G4MolecularDissociationChannel("DissociativeAttachment");
0334
0335   // Decay 1 : 2OH + H_2
0336   decCh1->AddProduct(H2);
0337   decCh1->AddProduct(OHm);
0338   decCh1->AddProduct(OH);
0339   decCh1->SetProbability(1);
0340   decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::DissociativeAttachment);
0341
0342   *occ = *(water->GetGroundStateElectronOccupancy());
0343   occ->AddElectron(5, 1);  // H_2O^-
0344   water->NewConfigurationWithElectronOccupancy("DissociativeAttachment", *occ);
0345   water->AddDecayChannel("DissociativeAttachment", decCh1);
0346
0347   //////////////////////////////////////////////////////////
0348   //            Electron-hole recombination               //
0349   //////////////////////////////////////////////////////////
0350   decCh1 = new G4MolecularDissociationChannel("H2Ovib_DissociationDecay1");
0351   decCh2 = new G4MolecularDissociationChannel("H2Ovib_DissociationDecay2");
0352   decCh3 = new G4MolecularDissociationChannel("H2Ovib_DissociationDecay3");
0353
0354   // Decay 1 : 2OH + H_2
0355   decCh1->AddProduct(H2);
0356   decCh1->AddProduct(OH);
0357   decCh1->AddProduct(OH);
0358   decCh1->SetProbability(0.15);
0359   decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::B1A1_DissociationDecay);
0360
0361   // Decay 2 : OH + H
0362   decCh2->AddProduct(OH);
0363   decCh2->AddProduct(H);
0364   decCh2->SetProbability(0.55);
0365   decCh2->SetDisplacementType(G4DNAWaterDissociationDisplacer::A1B1_DissociationDecay);
0366
0367   // Decay 3 : relaxation
0368   decCh3->SetProbability(0.30);
0369
0370   const auto pH2Ovib = G4H2O::Definition()->NewConfiguration("H2Ovib");
0371   assert(pH2Ovib != nullptr);
0372
0373   water->AddDecayChannel(pH2Ovib, decCh1);
0374   water->AddDecayChannel(pH2Ovib, decCh2);
0375   water->AddDecayChannel(pH2Ovib, decCh3);
0376
0377   delete occ;
0378 }
0379
0380 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0381
0382 void EmDNAChemistry::ConstructReactionTable(G4DNAMolecularReactionTable* theReactionTable)
0383 {
0384   // Construct chemical reactions from user input file
0385   ParserChemReaction parser;
0386   parser.ReadReactionFile(fReactionTableName);
0387
0388   auto ListReactant1 = parser.GetListReactant1();
0389   auto ListReactant2 = parser.GetListReactant2();
0390   auto ListProduct = parser.GetListProduct();
0391   auto ListRate = parser.GetListRate();
0392
0393   const G4int Ntype = 5;
0394   for (size_t i = 0; i < Ntype; i++) {
0395     for (size_t j = 0; j < ListReactant1[i].size(); j++) {
0396       G4MolecularConfiguration* Reactant1 =
0397         G4MoleculeTable::Instance()->GetConfiguration(ListReactant1[i][j]);
0398       G4MolecularConfiguration* Reactant2 =
0399         G4MoleculeTable::Instance()->GetConfiguration(ListReactant2[i][j]);
0400
0401       auto pReactionData = new G4DNAMolecularReactionData(ListRate[i][j], Reactant1, Reactant2);
0402
0403       for (size_t k = 0; k < ListProduct[i][j].size(); k++) {
0404         G4MolecularConfiguration* Product =
0405           G4MoleculeTable::Instance()->GetConfiguration(ListProduct[i][j][k]);
0406         pReactionData->AddProduct(Product);
0407       }
0408       // Reaction type II and IV
0409       if (i == 1 || i == 3) {
0410         pReactionData->SetReactionType(1);
0411       }
0412
0413       theReactionTable->SetReaction(pReactionData);
0414     }
0415   }
0416 }
0417
0418 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0419
0420 void EmDNAChemistry::ConstructProcess()
0421 {
0422   auto ph = G4PhysicsListHelper::GetPhysicsListHelper();
0423   G4VProcess* process =
0424     G4ProcessTable::GetProcessTable()->FindProcess("e-_G4DNAVibExcitation", "e-");
0425   if (process) {
0426     auto vibExcitation = (G4DNAVibExcitation*)process;
0427     G4VEmModel* model = vibExcitation->EmModel();
0428     auto sancheExcitationMod = dynamic_cast<G4DNASancheExcitationModel*>(model);
0429     if (sancheExcitationMod) {
0430       sancheExcitationMod->ExtendLowEnergyLimit(0.025 * eV);
0431     }
0432   }
0433   process = G4ProcessTable::GetProcessTable()->FindProcess("e-_G4DNAElectronSolvation", "e-");
0434
0435   if (process == nullptr) {
0436     ph->RegisterProcess(new G4DNAElectronSolvation("e-_G4DNAElectronSolvation"),
0437                         G4Electron::Definition());
0438   }
0439   auto* theMoleculeTable = G4MoleculeTable::Instance();
0440   G4MoleculeDefinitionIterator iterator = theMoleculeTable->GetDefintionIterator();
0441   iterator.reset();
0442   while (iterator()) {
0443     G4MoleculeDefinition* moleculeDef = iterator.value();
0444     if (moleculeDef == G4H2O::Definition()) {
0445       moleculeDef->GetProcessManager()->AddRestProcess(new G4DNAElectronHoleRecombination(), 2);
0446
0447       auto dissociationProcess = new G4DNAMolecularDissociation("H2O_DNAMolecularDecay");
0448       dissociationProcess->SetDisplacer(moleculeDef, new G4DNAWaterDissociationDisplacer);
0449       // dissociationProcess->SetVerboseLevel(3);
0450       moleculeDef->GetProcessManager()->AddRestProcess(dissociationProcess, 1);
0451     }
0452   }
0453   G4DNAChemistryManager::Instance()->Initialize();
0454 }
0455
0456 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0457
0458 void EmDNAChemistry::ConstructTimeStepModel(G4DNAMolecularReactionTable*
0459                                             /*reactionTable*/)
0460 {
0461   auto irt = new G4DNAMolecularIRTModel();
0462   RegisterTimeStepModel(irt, 0);
0463 }
0464
0465 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0466
0467 }  // namespace scavenger