| 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 //
27 /*
28 * G4DNAIRTMoleculeEncounterStepper.cc
29 *
30 * Created on: Jul 23, 2019
31 * Author: W. G. Shin
32 * J. Ramos-Mendez and B. Faddegon
33 */
34
35 #include "G4DNAIRTMoleculeEncounterStepper.hh"
36
37 #include "G4DNAMolecularReactionTable.hh"
38 #include "G4H2O.hh"
39 #include "G4ITReaction.hh"
40 #include "G4MolecularConfiguration.hh"
41 #include "G4MoleculeFinder.hh"
42 #include "G4Scheduler.hh"
43 #include "G4UnitsTable.hh"
44 #include "G4VDNAReactionModel.hh"
45 #include "G4memory.hh"
46
47 #include <memory>
48
49 using namespace std;
50 using namespace CLHEP;
51
52 //#define DEBUG_MEM
53
54 #ifdef DEBUG_MEM
55 #include "G4MemStat.hh"
56 using namespace G4MemStat;
57 #endif
58
59 G4DNAIRTMoleculeEncounterStepper::Utils::Utils(const G4Track& tA,
60 const G4MolecularConfiguration* pMoleculeB)
61 : fpTrackA(tA)
62 , fpMoleculeB(pMoleculeB)
63 {
64 fpMoleculeA = GetMolecule(tA);
65 fDA = fpMoleculeA->GetDiffusionCoefficient();
66 fDB = fpMoleculeB->GetDiffusionCoefficient();
67 fConstant = 8 * (fDA + fDB + 2 * sqrt(fDA * fDB));
68 }
69
70 G4DNAIRTMoleculeEncounterStepper::G4DNAIRTMoleculeEncounterStepper()
71 :
72 fMolecularReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
73 {
74 fpTrackContainer = G4ITTrackHolder::Instance();
75 fReactionSet = G4ITReactionSet::Instance();
76 }
77
78 G4DNAIRTMoleculeEncounterStepper::~G4DNAIRTMoleculeEncounterStepper() = default;
79
80 void G4DNAIRTMoleculeEncounterStepper::Prepare()
81 {
82 fSampledMinTimeStep = DBL_MAX;
83 if(G4Scheduler::Instance()->GetGlobalTime() == G4Scheduler::Instance()->GetStartTime()){
84 G4VITTimeStepComputer::Prepare();
85 G4MoleculeFinder::Instance()->UpdatePositionMap();
86 }
87 }
88
89 void G4DNAIRTMoleculeEncounterStepper::InitializeForNewTrack()
90 {
91 if (fReactants)
92 {
93 fReactants.reset();
94 }
95 fSampledMinTimeStep = DBL_MAX;
96 fHasAlreadyReachedNullTime = false;
97 }
98
99 template<typename T>
100 inline bool IsInf(T value)
101 {
102 return std::numeric_limits<T>::has_infinity
103 && value == std::numeric_limits<T>::infinity();
104 }
105
106 G4double
107 G4DNAIRTMoleculeEncounterStepper::CalculateStep(const G4Track& trackA,
108 const G4double& userMinTimeStep)
109 {
110
111 auto pMoleculeA = GetMolecule(trackA);
112 InitializeForNewTrack();
113 fUserMinTimeStep = userMinTimeStep;
114
115 #ifdef G4VERBOSE
116 if (fVerbose != 0)
117 {
118 G4cout
119 << "_______________________________________________________________________"
120 << G4endl;
121 G4cout << "G4DNAMoleculeEncounterStepper::CalculateStep" << G4endl;
122 G4cout << "Check done for molecule : " << pMoleculeA->GetName()
123 << " (" << trackA.GetTrackID() << ") "
124 << G4endl;
125 }
126 #endif
127
128 //__________________________________________________________________
129 // Retrieve general informations for making reactions
130 auto pMolConfA = pMoleculeA->GetMolecularConfiguration();
131
132 const auto pReactantList = fMolecularReactionTable->CanReactWith(pMolConfA);
133
134 if (pReactantList == nullptr)
135 {
136 #ifdef G4VERBOSE
137 // DEBUG
138 if (fVerbose > 1)
139 {
140 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
141 G4cout << "!!! WARNING" << G4endl;
142 G4cout << "G4MoleculeEncounterStepper::CalculateStep will return infinity "
143 "for the reaction because the molecule "
144 << pMoleculeA->GetName()
145 << " does not have any reactants given in the reaction table."
146 << G4endl;
147 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
148 }
149 #endif
150 return DBL_MAX;
151 }
152
153 auto nbReactives = (G4int)pReactantList->size();
154
155 if (nbReactives == 0)
156 {
157 #ifdef G4VERBOSE
158 // DEBUG
159 if (fVerbose != 0)
160 {
161 // TODO replace with the warning mode of G4Exception
162 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
163 G4cout << "!!! WARNING" << G4endl;
164 G4cout << "G4MoleculeEncounterStepper::CalculateStep will return infinity "
165 "for the reaction because the molecule "
166 << pMoleculeA->GetName()
167 << " does not have any reactants given in the reaction table."
168 << "This message can also result from a wrong implementation of the reaction table."
169 << G4endl;
170 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
171 }
172 #endif
173 return DBL_MAX;
174 }
175
176 fReactants = std::make_shared<vector<G4Track*>>();
177 fReactionModel->Initialise(pMolConfA, trackA);
178
179 //__________________________________________________________________
180 // Start looping on possible reactants
181 for (G4int i = 0; i < nbReactives; ++i)
182 {
183 auto pMoleculeB = (*pReactantList)[i];
184
185 //______________________________________________________________
186 // Retrieve reaction range
187 const G4double R = fReactionModel->GetReactionRadius(i);
188
189 //______________________________________________________________
190 // Use KdTree algorithm to find closest reactants
191 G4KDTreeResultHandle resultsNearest(
192 G4MoleculeFinder::Instance()->FindNearest(pMoleculeA,
193 pMoleculeB->GetMoleculeID()));
194
195 if (static_cast<int>(resultsNearest) == 0) continue;
196
197 G4double r2 = resultsNearest->GetDistanceSqr();
198 Utils utils(trackA, pMoleculeB);
199
200 if (r2 <= R * R) // ==> Record in range
201 {
202 // Entering in this condition may due to the fact that molecules are very close
203 // to each other
204 // Therefore, if we only take the nearby reactant into account, it might have already
205 // reacted. Instead, we will take all possible reactants that satisfy the condition r<R
206
207 if (!fHasAlreadyReachedNullTime)
208 {
209 fReactants->clear();
210 fHasAlreadyReachedNullTime = true;
211 }
212
213 fSampledMinTimeStep = 0.;
214 G4KDTreeResultHandle resultsInRange(
215 G4MoleculeFinder::Instance()->FindNearestInRange(pMoleculeA,
216 pMoleculeB->GetMoleculeID(),
217 R));
218 CheckAndRecordResults(utils,
219 #ifdef G4VERBOSE
220 R,
221 #endif
222 resultsInRange);
223 }
224 else
225 {
226 G4double r = sqrt(r2);
227 G4double tempMinET = pow(r - R, 2) / utils.fConstant;
228 // constant = 16 * (fDA + fDB + 2*sqrt(fDA*fDB))
229
230 if (tempMinET <= fSampledMinTimeStep)
231 {
232 if (fUserMinTimeStep < DBL_MAX/*IsInf(fUserMinTimeStep) == false*/
233 && tempMinET <= fUserMinTimeStep) // ==> Record in range
234 {
235 if (fSampledMinTimeStep > fUserMinTimeStep)
236 {
237 fReactants->clear();
238 }
239
240 fSampledMinTimeStep = fUserMinTimeStep;
241
242 G4double range = R + sqrt(fUserMinTimeStep*utils.fConstant);
243
244 G4KDTreeResultHandle resultsInRange(
245 G4MoleculeFinder::Instance()->
246 FindNearestInRange(pMoleculeA,
247 pMoleculeB->GetMoleculeID(),
248 range));
249
250 CheckAndRecordResults(utils,
251 #ifdef G4VERBOSE
252 range,
253 #endif
254 resultsInRange);
255 }
256 else // ==> Record nearest
257 {
258 if (tempMinET < fSampledMinTimeStep)
259 // to avoid cases where fSampledMinTimeStep == tempMinET
260 {
261 fSampledMinTimeStep = tempMinET;
262 fReactants->clear();
263 }
264
265 CheckAndRecordResults(utils,
266 #ifdef G4VERBOSE
267 R,
268 #endif
269 resultsNearest);
270 }
271 }
272 }
273 }
274
275 #ifdef G4VERBOSE
276 if (fVerbose != 0)
277 {
278 G4cout << "G4MoleculeEncounterStepper::CalculateStep will finally return :"
279 << G4BestUnit(fSampledMinTimeStep, "Time") << G4endl;
280
281 if (fVerbose > 1)
282 {
283 G4cout << "Selected reactants for trackA: " << pMoleculeA->GetName()
284 << " (" << trackA.GetTrackID() << ") are: ";
285
286 vector<G4Track*>::iterator it;
287 for (it = fReactants->begin(); it != fReactants->end(); it++)
288 {
289 G4Track* trackB = *it;
290 G4cout << GetMolecule(trackB)->GetName() << " ("
291 << trackB->GetTrackID() << ") \t ";
292 }
293 G4cout << G4endl;
294 }
295 }
296 #endif
297 return fSampledMinTimeStep;
298 }
299
300
301
302
303 void G4DNAIRTMoleculeEncounterStepper::CheckAndRecordResults(const Utils& utils,
304 #ifdef G4VERBOSE
305 const G4double R,
306 #endif
307 G4KDTreeResultHandle& results)
308 {
309 if (static_cast<int>(results) == 0)
310 {
311 #ifdef G4VERBOSE
312 if (fVerbose > 1)
313 {
314 G4cout << "No molecule " << utils.fpMoleculeB->GetName()
315 << " found to react with " << utils.fpMoleculeA->GetName()
316 << G4endl;
317 }
318 #endif
319 return;
320 }
321
322 for (results->Rewind(); !results->End(); results->Next())
323 {
324 auto reactiveB = results->GetItem<G4IT>();
325
326 if (reactiveB == nullptr)
327 {
328 continue;
329 }
330
331 G4Track *trackB = reactiveB->GetTrack();
332
333 if (trackB == nullptr)
334 {
335 G4ExceptionDescription exceptionDescription;
336 exceptionDescription
337 << "The reactant B found using the MoleculeFinder does not have a valid "
338 "track attached to it. If this is done on purpose, please do "
339 "not record this molecule in the MoleculeFinder."
340 << G4endl;
341 G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
342 "MoleculeEncounterStepper001", FatalErrorInArgument,
343 exceptionDescription);
344 continue;
345 }
346
347 if (trackB->GetTrackStatus() != fAlive)
348 {
349 continue;
350 }
351
352 if (trackB == &utils.fpTrackA)
353 {
354 G4ExceptionDescription exceptionDescription;
355 exceptionDescription
356 << "A track is reacting with itself (which is impossible) ie fpTrackA == trackB"
357 << G4endl;
358 exceptionDescription << "Molecule A (and B) is of type : "
359 << utils.fpMoleculeA->GetName() << " with trackID : "
360 << utils.fpTrackA.GetTrackID() << G4endl;
361
362 G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
363 "MoleculeEncounterStepper003", FatalErrorInArgument,
364 exceptionDescription);
365
366 }
367
368 if (fabs(trackB->GetGlobalTime() - utils.fpTrackA.GetGlobalTime())
369 > utils.fpTrackA.GetGlobalTime() * (1 - 1 / 100))
370 {
371 // DEBUG
372 G4ExceptionDescription exceptionDescription;
373 exceptionDescription
374 << "The interacting tracks are not synchronized in time" << G4endl;
375 exceptionDescription
376 << "trackB->GetGlobalTime() != fpTrackA.GetGlobalTime()" << G4endl;
377
378 exceptionDescription << "fpTrackA : trackID : " << utils.fpTrackA.GetTrackID()
379 << "\t Name :" << utils.fpMoleculeA->GetName()
380 << "\t fpTrackA->GetGlobalTime() = "
381 << G4BestUnit(utils.fpTrackA.GetGlobalTime(), "Time") << G4endl;
382
383 exceptionDescription << "trackB : trackID : " << trackB->GetTrackID()
384 << "\t Name :" << utils.fpMoleculeB->GetName()
385 << "\t trackB->GetGlobalTime() = "
386 << G4BestUnit(trackB->GetGlobalTime(), "Time") << G4endl;
387
388 G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
389 "MoleculeEncounterStepper004", FatalErrorInArgument,
390 exceptionDescription);
391 }
392
393 #ifdef G4VERBOSE
394 if (fVerbose > 1)
395 {
396
397 G4double r2 = results->GetDistanceSqr();
398 G4cout << "\t ************************************************** " << G4endl;
399 G4cout << "\t Reaction between "
400 << utils.fpMoleculeA->GetName() << " (" << utils.fpTrackA.GetTrackID() << ") "
401 << " & " << utils.fpMoleculeB->GetName() << " (" << trackB->GetTrackID() << "), "
402 << "Interaction Range = "
403 << G4BestUnit(R, "Length") << G4endl;
404 G4cout << "\t Real distance between reactants = "
405 << G4BestUnit((utils.fpTrackA.GetPosition() - trackB->GetPosition()).mag(), "Length") << G4endl;
406 G4cout << "\t Distance between reactants calculated by nearest neighbor algorithm = "
407 << G4BestUnit(sqrt(r2), "Length") << G4endl;
408
409 }
410 #endif
411
412 fReactants->push_back(trackB);
413 }
414 }
415
416 void G4DNAIRTMoleculeEncounterStepper::SetReactionModel(G4VDNAReactionModel* pReactionModel)
417 {
418 fReactionModel = pReactionModel;
419 }
420
421 G4VDNAReactionModel* G4DNAIRTMoleculeEncounterStepper::GetReactionModel()
422 {
423 return fReactionModel;
424 }
425
426 void G4DNAIRTMoleculeEncounterStepper::SetVerbose(int flag)
427 {
428 fVerbose = flag;
429 }
430
431 G4double G4DNAIRTMoleculeEncounterStepper::CalculateMinTimeStep(G4double currentGlobalTime, G4double definedMinTimeStep){
432
433 G4bool start = true;
434 G4bool active = false;
435
436 fUserMinTimeStep = definedMinTimeStep;
437
438 if(fReactionSet->Empty()){
439 if(currentGlobalTime == G4Scheduler::Instance()->GetStartTime()){
440
441 for (auto pTrack : *fpTrackContainer->GetMainList())
442 {
443 if (pTrack == nullptr)
444 {
445 G4ExceptionDescription exceptionDescription;
446 exceptionDescription << "No track found.";
447 G4Exception("G4Scheduler::CalculateMinStep", "ITScheduler006",
448 FatalErrorInArgument, exceptionDescription);
449 continue;
450 }
451
452 G4TrackStatus trackStatus = pTrack->GetTrackStatus();
453 if (trackStatus == fStopAndKill || trackStatus == fStopButAlive)
454 {
455 start = false;
456 continue;
457 }
458 active = true;
459 }
460
461 if(start){
462 return -1;
463 }if(!active){
464 G4Scheduler::Instance()->Stop();
465 return fSampledMinTimeStep;
466 }
467 return fSampledMinTimeStep;
468 }
469 for (auto pTrack : *fpTrackContainer->GetMainList())
470 {
471 pTrack->SetGlobalTime(G4Scheduler::Instance()->GetEndTime());
472 }
473 return fSampledMinTimeStep;
474 }
475
476 const auto& fReactionSetInTime = fReactionSet->GetReactionsPerTime();
477 fSampledMinTimeStep = fReactionSetInTime.begin()->get()->GetTime() - currentGlobalTime;
478
479 return fSampledMinTimeStep;
480 }
481