| 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 //
26 // Author: Hoang TRAN
27
28 #include "G4DiffusionControlledReactionModel.hh"
29 #include "G4Track.hh"
30 #include "G4DNAMolecularReactionTable.hh"
31 #include "G4PhysicalConstants.hh"
32 #include "G4Exp.hh"
33 #include "G4IRTUtils.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4Electron_aq.hh"
36 #include "Randomize.hh"
37 #include "G4Molecule.hh"
38 #include "G4ErrorFunction.hh"
39
40 G4DiffusionControlledReactionModel::G4DiffusionControlledReactionModel() = default;
41
42 G4DiffusionControlledReactionModel::~G4DiffusionControlledReactionModel() = default;
43
44 void G4DiffusionControlledReactionModel::Initialise(
45 const G4MolecularConfiguration* pMolecule, const G4Track&)
46 {
47 fpReactionData = fpReactionTable->GetReactionData(pMolecule);
48 }
49
50 void G4DiffusionControlledReactionModel::InitialiseToPrint(
51 const G4MolecularConfiguration* pMolecule)
52 {
53 fpReactionData = fpReactionTable->GetReactionData(pMolecule);
54 }
55
56 G4double G4DiffusionControlledReactionModel::GetReactionRadius(
57 const G4MolecularConfiguration* pMol1, const G4MolecularConfiguration* pMol2)
58 {
59 auto reactionData = fpReactionTable->GetReactionData(pMol1, pMol2);
60 if(reactionData == nullptr)
61 {
62 G4ExceptionDescription exceptionDescription;
63 exceptionDescription << "No reactionData"
64 << " for : " << pMol1->GetName() << " and "
65 << pMol2->GetName();
66 G4Exception("G4DiffusionControlledReactionModel"
67 "::GetReactionRadius()",
68 "G4DiffusionControlledReactionModel00", FatalException,
69 exceptionDescription);
70 return 0.;
71 }
72 return reactionData->GetEffectiveReactionRadius();
73 }
74
75 G4double G4DiffusionControlledReactionModel::GetReactionRadius(const G4int& i)
76 {
77 auto pMol1 = (*fpReactionData)[i]->GetReactant1();
78 auto pMol2 = (*fpReactionData)[i]->GetReactant2();
79 return GetReactionRadius(pMol1, pMol2);
80 }
81
82 G4double G4DiffusionControlledReactionModel::GetTimeToEncounter(
83 const G4Track& trackA, const G4Track& trackB)
84 {
85 auto pMolConfA = GetMolecule(trackA)->GetMolecularConfiguration();
86 auto pMolConfB = GetMolecule(trackB)->GetMolecularConfiguration();
87
88 G4double D =
89 pMolConfA->GetDiffusionCoefficient() + pMolConfB->GetDiffusionCoefficient();
90
91 if(D == 0)
92 {
93 G4ExceptionDescription exceptionDescription;
94 exceptionDescription << "The total diffusion coefficient for : "
95 << pMolConfA->GetName() << " and "
96 << pMolConfB->GetName() << " is null ";
97 G4Exception("G4DiffusionControlledReactionModel"
98 "::GetTimeToEncounter()",
99 "G4DiffusionControlledReactionModel03", FatalException,
100 exceptionDescription);
101 }
102
103 auto reactionData = G4DNAMolecularReactionTable::Instance()->GetReactionData(
104 pMolConfA, pMolConfB);
105 G4double kobs = reactionData->GetObservedReactionRateConstant();
106 G4double distance = (trackA.GetPosition() - trackB.GetPosition()).mag();
107 G4double SmoluchowskiRadius = reactionData->GetEffectiveReactionRadius();
108
109 if(distance == 0 || distance < SmoluchowskiRadius)
110 {
111 G4ExceptionDescription exceptionDescription;
112 exceptionDescription << "distance = " << distance << " is uncorrected with "
113 << " Reff = " << SmoluchowskiRadius
114 << " for : " << pMolConfA->GetName() << " and "
115 << pMolConfB->GetName();
116 G4Exception("G4DiffusionControlledReactionModel"
117 "::GetTimeToEncounter()",
118 "G4DiffusionControlledReactionModel02", FatalException,
119 exceptionDescription);
120 }
121 else
122 {
123 G4double Winf = SmoluchowskiRadius / distance;
124 G4double U = G4UniformRand();
125 G4double X = 0;
126 G4double irt_1 = -1.0 * ps;
127 if(Winf > 0 && U < Winf)
128 {
129 G4double erfcIn = G4ErrorFunction::erfcInv(U / Winf);
130 if(erfcIn != 0)
131 {
132 G4double d =
133 (distance - SmoluchowskiRadius) / erfcIn;
134 irt_1 = (1.0 / (4 * D)) * d * d;
135 }
136 }
137
138 if(reactionData->GetReactionType() == 0) // Totally diffused contr
139 {
140 return irt_1;
141 }
142
143 if(irt_1 < 0)
144 {
145 return irt_1;
146 }
147
148 G4double kdif = 4 * CLHEP::pi * D * SmoluchowskiRadius * Avogadro;
149
150 if(pMolConfA == pMolConfB)
151 {
152 kdif /= 2;
153 }
154 G4double kact = G4IRTUtils::GetKact(kobs, kdif);
155 G4double sumOfk = kact + kdif;
156 if(sumOfk != 0)
157 {
158 G4double rateFactor = kact / sumOfk;
159 if(G4UniformRand() > rateFactor)
160 {
161 return -1.0 * ps;
162 }
163 G4double Y = std::abs(G4RandGauss::shoot(0.0, std::sqrt(2)));
164
165 if(Y > 0)
166 {
167 X = -(G4Log(G4UniformRand())) / Y;
168 }
169 G4double f = X * SmoluchowskiRadius * kdif / sumOfk;
170 G4double irt_2 = (f * f) / D;
171 return irt_1 + irt_2;
172 }
173 }
174 return -1.0 * ps;
175 }
176