| 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 #include "G4DNAGillespieDirectMethod.hh"
27 #include "Randomize.hh"
28 #include "G4PhysicalConstants.hh"
29 #include <memory>
30 #include <tuple>
31 #include "G4DNAEventSet.hh"
32 #include "G4UnitsTable.hh"
33 #include "G4DNAScavengerMaterial.hh"
34 #include "G4Scheduler.hh"
35 #include "G4DNAMolecularReactionTable.hh"
36 #include "G4ChemEquilibrium.hh"
37
38 G4DNAGillespieDirectMethod::G4DNAGillespieDirectMethod()
39 : fMolecularReactions(G4DNAMolecularReactionTable::Instance())
40 {}
41
42 G4DNAGillespieDirectMethod::~G4DNAGillespieDirectMethod() = default;
43
44 void G4DNAGillespieDirectMethod::SetEventSet(G4DNAEventSet* pEventSet)
45 {
46 fpEventSet = pEventSet;
47 }
48
49 //#define DEBUG 1
50
51 G4double G4DNAGillespieDirectMethod::VolumeOfNode(const Voxel& voxel)
52 {
53 auto box = std::get<1>(voxel);
54 auto LengthY = box.Getyhi() - box.Getylo();
55 auto LengthX = box.Getxhi() - box.Getxlo();
56 auto LengthZ = box.Getzhi() - box.Getzlo();
57 G4double V = LengthY * LengthX * LengthZ;
58 if(V <= 0)
59 {
60 G4ExceptionDescription exceptionDescription;
61 exceptionDescription << "V > 0 !! ";
62 G4Exception("G4DNAGillespieDirectMethod::VolumeOfNode",
63 "G4DNAGillespieDirectMethod03", FatalErrorInArgument,
64 exceptionDescription);
65 }
66 return V;
67 }
68 G4double G4DNAGillespieDirectMethod::PropensityFunction(const Voxel& voxel,
69 MolType moleType)
70 {
71 if(moleType->GetDiffusionCoefficient() == 0)
72 {
73 return 0.;
74 }
75 const auto& node = std::get<2>(voxel);
76 const auto& box = std::get<1>(voxel);
77
78 G4double alpha = 0;
79 auto it = node.find(moleType);
80 if(it != node.end())
81 {
82 auto LengthY = box.Getyhi() - box.Getylo();
83 G4double d = it->first->GetDiffusionCoefficient() / std::pow(LengthY, 2);
84 alpha = d * it->second;
85
86 #ifdef DEBUG
87 G4cout << it->first->GetName() << " " << it->second
88 << " D : " << it->first->GetDiffusionCoefficient()
89 << " LengthY : " << LengthY << " PropensityFunction : " << alpha
90 << G4endl;
91 #endif
92 }
93 return alpha;
94 }
95
96 G4double G4DNAGillespieDirectMethod::PropensityFunction(const Voxel& voxel,
97 ReactionData* data)
98 {
99 G4double value;
100 auto ConfA = data->GetReactant1();
101 auto ConfB = data->GetReactant2();
102 G4double scavengerNumber = 0;
103 auto typeANumber = FindScavenging(voxel, ConfA, scavengerNumber)
104 ? scavengerNumber
105 : ComputeNumberInNode(voxel, ConfA);
106
107 auto typeBNumber = FindScavenging(voxel, ConfB, scavengerNumber)
108 ? scavengerNumber
109 : ComputeNumberInNode(voxel, ConfB);
110
111 if(typeANumber == 0 || typeBNumber == 0)
112 {
113 return 0;
114 }
115
116 auto k =
117 data->GetObservedReactionRateConstant() / (Avogadro * VolumeOfNode(voxel));
118 if(ConfA == ConfB)
119 {
120 value = typeANumber * (typeBNumber - 1) * k;
121 }
122 else
123 {
124 value = typeANumber * typeBNumber * k;
125 }
126
127 if(value < 0)
128 {
129 G4ExceptionDescription exceptionDescription;
130 exceptionDescription
131 << "G4DNAGillespieDirectMethod::PropensityFunction for : "
132 << ConfA->GetName() << "(" << typeANumber << ") + " << ConfB->GetName()
133 << "(" << typeBNumber << ") : propensity : " << value
134 << " GetObservedReactionRateConstant : "
135 << data->GetObservedReactionRateConstant()
136 << " GetEffectiveReactionRadius : "
137 << G4BestUnit(data->GetEffectiveReactionRadius(), "Length")
138 << " k : " << k << " volume : " << VolumeOfNode(voxel) << G4endl;
139 G4Exception("G4DNAGillespieDirectMethod::PropensityFunction",
140 "G4DNAGillespieDirectMethod013", FatalErrorInArgument,
141 exceptionDescription);
142 }
143
144 #ifdef DEBUG
145 if(value > 0)
146 G4cout << "G4DNAGillespieDirectMethod::PropensityFunction for : "
147 << ConfA->GetName() << "(" << typeANumber << ") + "
148 << ConfB->GetName() << "(" << typeBNumber
149 << ") : propensity : " << value
150 << " Time to Reaction : " << G4BestUnit(timeToReaction, "Time")
151 << " k : " << k << " Index : " << index << G4endl;
152 #endif
153
154 return value;
155 }
156
157 void G4DNAGillespieDirectMethod::Initialize()
158 {
159 // for Scavenger
160 fpScavengerMaterial = dynamic_cast<G4DNAScavengerMaterial*>(
161 G4Scheduler::Instance()->GetScavengerMaterial());
162 fEquilibriumProcesses.emplace(
163 std::make_pair(6, std::make_unique<G4ChemEquilibrium>(6, 10 * CLHEP::us)));//reactionType6 and 10 * us
164 fEquilibriumProcesses.emplace(
165 std::make_pair(7, std::make_unique<G4ChemEquilibrium>(7, 10 * CLHEP::us)));//reactionType6 and 10 * us
166 fEquilibriumProcesses.emplace(
167 std::make_pair(8, std::make_unique<G4ChemEquilibrium>(8, 10 * CLHEP::us)));//reactionType6 and 10 * us
168 for(auto& it : fEquilibriumProcesses)
169 {
170 it.second->Initialize();
171 it.second->SetVerbose(fVerbose);
172 }
173 }
174
175 void G4DNAGillespieDirectMethod::CreateEvents()
176 {
177 auto begin = fpMesh->begin();
178 auto end = fpMesh->end();
179 for(; begin != end; begin++)
180 {
181 auto index = std::get<0>(*begin);
182 #ifdef DEBUG
183 fpMesh->PrintVoxel(index);
184 #endif
185 CreateEvent(index);
186 }
187 }
188
189 void G4DNAGillespieDirectMethod::SetTimeStep(const G4double& stepTime)
190 {
191 fTimeStep = stepTime;
192 }
193 void G4DNAGillespieDirectMethod::CreateEvent(const Index& index)
194 {
195 const auto& voxel = fpMesh->GetVoxel(index);
196 if(std::get<2>(voxel).empty())
197 {
198 G4ExceptionDescription exceptionDescription;
199 exceptionDescription << "This voxel : " << index
200 << " is not ready to make event" << G4endl;
201 G4Exception("G4DNAGillespieDirectMethod::CreateEvent",
202 "G4DNAGillespieDirectMethod05", FatalErrorInArgument,
203 exceptionDescription);
204 }
205 G4double r1 = G4UniformRand();
206 G4double r2 = G4UniformRand();
207 G4double dAlpha0 = DiffusiveJumping(voxel);
208 G4double rAlpha0 = Reaction(voxel);
209 G4double alphaTotal = dAlpha0 + rAlpha0;
210
211 if(alphaTotal == 0)
212 {
213 return;
214 }
215 auto timeStep = ((1.0 / (alphaTotal)) * std::log(1.0 / r1)) + fTimeStep;
216
217 #ifdef DEBUG
218 G4cout << "r2 : " << r2 << " rAlpha0 : " << rAlpha0
219 << " dAlpha0 : " << dAlpha0 << " rAlpha0 / (dAlpha0 + rAlpha0) : "
220 << rAlpha0 / (dAlpha0 + rAlpha0) << G4endl;
221 #endif
222 if(r2 < rAlpha0 / alphaTotal)
223 {
224 if(fVerbose > 1)
225 {
226 G4cout << "=>>>>reaction at : " << timeStep << " timeStep : "
227 << G4BestUnit(((1.0 / alphaTotal) * std::log(1.0 / r1)), "Time")
228 << G4endl;
229 }
230 auto rSelectedIter = fReactionDataMap.upper_bound(r2 * alphaTotal);
231 fpEventSet->CreateEvent(timeStep, index, rSelectedIter->second);
232 }
233 else if(dAlpha0 > 0)
234 {
235 if(fVerbose > 1)
236 {
237 G4cout << "=>>>>jumping at : " << timeStep << " timeStep : "
238 << G4BestUnit(((1.0 / alphaTotal) * std::log(1.0 / r1)), "Time")
239 << G4endl;
240 }
241
242 auto dSelectedIter = fJumpingDataMap.upper_bound(r2 * alphaTotal - rAlpha0);
243 auto pDSelected =
244 std::make_unique<std::pair<MolType, Index>>(dSelectedIter->second);
245 fpEventSet->CreateEvent(timeStep, index, std::move(pDSelected));
246 }
247 #ifdef DEBUG
248 G4cout << G4endl;
249 #endif
250 }
251
252 G4double G4DNAGillespieDirectMethod::Reaction(const Voxel& voxel)
253 {
254 fReactionDataMap.clear();
255 G4double alpha0 = 0;
256 const auto& dataList =
257 fMolecularReactions->GetVectorOfReactionData(); // shoud make a member
258 if(dataList.empty())
259 {
260 G4ExceptionDescription exceptionDescription;
261 exceptionDescription << "MolecularReactionTable empty" << G4endl;
262 G4Exception("G4DNAGillespieDirectMethod::Reaction",
263 "G4DNAGillespieDirectMethod01", FatalErrorInArgument,
264 exceptionDescription);
265 }
266
267 for(const auto& it : dataList)
268 {
269 if(!IsEquilibrium(it->GetReactionType()))
270 {
271 continue;
272 }
273 auto propensity = PropensityFunction(voxel, it);
274 if(propensity == 0)
275 {
276 continue;
277 }
278 alpha0 += propensity;
279 fReactionDataMap[alpha0] = it;
280 }
281 #ifdef DEBUG
282 G4cout << "Reaction :alpha0 : " << alpha0 << G4endl;
283 #endif
284 return alpha0;
285 }
286
287 G4double G4DNAGillespieDirectMethod::DiffusiveJumping(const Voxel& voxel)
288 {
289 fJumpingDataMap.clear();
290 G4double alpha0 = 0;
291 auto index = std::get<0>(voxel);
292 auto NeighboringVoxels = fpMesh->FindNeighboringVoxels(index);
293 if(NeighboringVoxels.empty())
294 {
295 return 0;
296 }
297 auto iter = G4MoleculeTable::Instance()->GetConfigurationIterator();
298 while(iter())
299 {
300 const auto* conf = iter.value();
301 auto propensity = PropensityFunction(voxel, conf);
302 if(propensity == 0)
303 {
304 continue;
305 }
306 for(const auto& it_Neighbor : NeighboringVoxels)
307 {
308 alpha0 += propensity;
309 fJumpingDataMap[alpha0] = std::make_pair(conf, it_Neighbor);
310 #ifdef DEBUG
311 G4cout << "mole : " << conf->GetName()
312 << " number : " << ComputeNumberInNode(index, conf)
313 << " propensity : " << propensity << " alpha0 : " << alpha0
314 << G4endl;
315 #endif
316 }
317 }
318 #ifdef DEBUG
319 G4cout << "DiffusiveJumping :alpha0 : " << alpha0 << G4endl;
320 #endif
321 return alpha0;
322 }
323
324 G4double G4DNAGillespieDirectMethod::ComputeNumberInNode(
325 const Voxel& voxel, MolType type) // depend node ?
326 {
327 if(type->GetDiffusionCoefficient() != 0)
328 {
329 const auto& node = std::get<2>(voxel);
330 const auto& it = node.find(type);
331 return (it != node.end()) ? (it->second) : 0;
332 }
333 return 0;
334 }
335
336 G4bool G4DNAGillespieDirectMethod::FindScavenging(const Voxel& voxel,
337 MolType moletype,
338 G4double& numberOfScavenger)
339 {
340 numberOfScavenger = 0;
341 if(fpScavengerMaterial == nullptr)
342 {
343 return false;
344 }
345 auto volumeOfNode = VolumeOfNode(voxel);
346 if(G4MoleculeTable::Instance()->GetConfiguration("H2O") == moletype)
347 {
348 auto factor = Avogadro * volumeOfNode;
349 numberOfScavenger = factor;
350 return true;
351 }
352
353 G4double totalNumber =
354 fpScavengerMaterial->GetNumberMoleculePerVolumeUnitForMaterialConf(
355 moletype);
356 if(totalNumber == 0)
357 {
358 return false;
359 }
360
361 G4double numberInDouble = volumeOfNode * std::floor(totalNumber) /
362 fpMesh->GetBoundingBox().Volume();
363 auto numberInInterg = (int64_t) (std::floor(numberInDouble));
364 G4double change = numberInDouble - numberInInterg;
365 G4UniformRand() > change ? numberOfScavenger = numberInInterg
366 : numberOfScavenger = numberInInterg + 1;
367 return true;
368 }
369
370 G4bool G4DNAGillespieDirectMethod::IsEquilibrium(const G4int& reactionType) const
371 {
372 auto reaction = fEquilibriumProcesses.find(reactionType);
373 if(reaction == fEquilibriumProcesses.end())
374 {
375 return true;
376 }else
377 {
378 return (reaction->second->GetEquilibriumStatus());
379 }
380
381 }
382
383 G4bool G4DNAGillespieDirectMethod::SetEquilibrium(const G4DNAMolecularReactionData* pReaction)
384 {
385 for(auto& it : fEquilibriumProcesses)
386 {
387 it.second->SetGlobalTime(fTimeStep);
388 it.second->SetEquilibrium(pReaction);
389 if(it.second->IsStatusChanged()) return true;
390 }
391 return false;
392 }
393
394 void G4DNAGillespieDirectMethod::ResetEquilibrium()
395 {
396 for(auto& it : fEquilibriumProcesses)
397 {
398 it.second->Reset();
399 }
400 }
401
402
403
404