| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 /*
27 * G4DNAIRT.cc
28 *
29 * Created on: Jul 23, 2019
30 * Author: W. G. Shin
31 * J. Ramos-Mendez and B. Faddego
32 */
33
34
35 #include "G4DNAIRT.hh"
36 #include "G4ErrorFunction.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4PhysicalConstants.hh"
39 #include "Randomize.hh"
40 #include "G4DNAMolecularReactionTable.hh"
41 #include "G4MolecularConfiguration.hh"
42 #include "G4Molecule.hh"
43 #include "G4ITReactionChange.hh"
44 #include "G4ITTrackHolder.hh"
45 #include "G4ITReaction.hh"
46 #include "G4Scheduler.hh"
47
48 using namespace std;
49
50 G4DNAIRT::G4DNAIRT() :
51 fMolReactionTable(reference_cast<const G4DNAMo
52 fpReactionModel(nullptr),
53 fTrackHolder(G4ITTrackHolder::Instance()),
54 fReactionSet(nullptr)
55 {
56 timeMin = G4Scheduler::Instance()->GetStartT
57 timeMax = G4Scheduler::Instance()->GetEndTim
58
59 fXMin = 1e9*nm;
60 fYMin = 1e9*nm;
61 fZMin = 1e9*nm;
62
63 fXMax = 0e0*nm;
64 fYMax = 0e0*nm;
65 fZMax = 0e0*nm;
66
67 fNx = 0;
68 fNy = 0;
69 fNz = 0;
70
71 xiniIndex = 0, yiniIndex = 0, ziniIndex = 0;
72 xendIndex = 0, yendIndex = 0, zendIndex = 0;
73
74 fRCutOff =
75 1.45 * nm + 2 * std::sqrt(8*9.46e9*nm*nm/s
76
77 erfc = new G4ErrorFunction();
78 }
79
80
81 G4DNAIRT::G4DNAIRT(G4VDNAReactionModel* pReact
82 : G4DNAIRT()
83 {
84 fpReactionModel = pReactionModel;
85 }
86
87 G4DNAIRT::~G4DNAIRT()
88 {
89 delete erfc;
90 }
91
92 void G4DNAIRT::Initialize(){
93
94 fTrackHolder = G4ITTrackHolder::Instance();
95
96 fReactionSet = G4ITReactionSet::Instance();
97 fReactionSet->CleanAllReaction();
98 fReactionSet->SortByTime();
99
100 spaceBinned.clear();
101
102 timeMin = G4Scheduler::Instance()->GetStartT
103 timeMax = G4Scheduler::Instance()->GetEndTim
104
105 xiniIndex = 0;
106 yiniIndex = 0;
107 ziniIndex = 0;
108 xendIndex = 0;
109 yendIndex = 0;
110 zendIndex = 0;
111
112 fXMin = 1e9*nm;
113 fYMin = 1e9*nm;
114 fZMin = 1e9*nm;
115
116 fXMax = 0e0*nm;
117 fYMax = 0e0*nm;
118 fZMax = 0e0*nm;
119
120 fNx = 0;
121 fNy = 0;
122 fNz = 0;
123
124 SpaceBinning(); // 1. binning the space
125 IRTSampling(); // 2. Sampling of the IRT
126
127 //hoang : if the first IRTSampling won't giv
128 if(fReactionSet->Empty())
129 {
130 for (auto pTrack : *fTrackHolder->GetMainL
131 {
132 pTrack->SetGlobalTime(G4Scheduler::Insta
133 }
134 }
135 }
136
137 void G4DNAIRT::SpaceBinning(){
138 auto it_begin = fTrackHolder->GetMainList()-
139 while(it_begin != fTrackHolder->GetMainList(
140
141 G4ThreeVector position = it_begin->GetPosi
142
143 if ( fXMin > position.x() ) fXMin = positi
144 if ( fYMin > position.y() ) fYMin = positi
145 if ( fZMin > position.z() ) fZMin = positi
146
147 if ( fXMax < position.x() ) fXMax = positi
148 if ( fYMax < position.y() ) fYMax = positi
149 if ( fZMax < position.z() ) fZMax = positi
150
151 ++it_begin;
152 }
153
154 fNx = G4int((fXMax-fXMin)/fRCutOff) == 0 ? 1
155 fNy = G4int((fYMax-fYMin)/fRCutOff) == 0 ? 1
156 fNz = G4int((fZMax-fZMin)/fRCutOff) == 0 ? 1
157
158 }
159
160 void G4DNAIRT::IRTSampling(){
161
162 auto it_begin = fTrackHolder->GetMainList()-
163 while(it_begin != fTrackHolder->GetMainList(
164 G4int I = FindBin(fNx, fXMin, fXMax, it_be
165 G4int J = FindBin(fNy, fYMin, fYMax, it_be
166 G4int K = FindBin(fNz, fZMin, fZMax, it_be
167
168 spaceBinned[I][J][K].push_back(*it_begin);
169
170 Sampling(*it_begin);
171 ++it_begin;
172 }
173 }
174
175 void G4DNAIRT::Sampling(G4Track* track){
176 G4Molecule* molA = G4Molecule::GetMolecule(t
177 const G4MolecularConfiguration* molConfA = m
178 if(molConfA->GetDiffusionCoefficient() == 0)
179
180 const vector<const G4MolecularConfiguration*
181 fMolReactionTable->CanReactWith(molConfA);
182
183 if(reactivesVector == nullptr) return;
184
185 G4double globalTime = G4Scheduler::Instance(
186 G4double minTime = timeMax;
187
188 xiniIndex = FindBin(fNx, fXMin, fXMax, track
189 xendIndex = FindBin(fNx, fXMin, fXMax, track
190 yiniIndex = FindBin(fNy, fYMin, fYMax, track
191 yendIndex = FindBin(fNy, fYMin, fYMax, track
192 ziniIndex = FindBin(fNz, fZMin, fZMax, track
193 zendIndex = FindBin(fNz, fZMin, fZMax, track
194
195 for ( int ii = xiniIndex; ii <= xendIndex; i
196 for ( int jj = yiniIndex; jj <= yendIndex;
197 for ( int kk = ziniIndex; kk <= zendInde
198
199 std::vector<G4Track*> spaceBin = space
200 for ( int n = 0; n < (int)spaceBinned[
201 if((spaceBin[n] == nullptr) || track
202 if(spaceBin[n]->GetTrackStatus() ==
203
204 G4Molecule* molB = G4Molecule::GetMo
205 if(molB == nullptr) continue;
206
207 const G4MolecularConfiguration* molC
208 if(molConfB->GetDiffusionCoefficient
209
210 auto it = std::find(reactivesVector-
211 if(it == reactivesVector->end()) con
212
213 G4ThreeVector orgPosB = spaceBin[n]-
214 G4double dt = track->GetGlobalTime()
215 G4ThreeVector newPosB = orgPosB;
216
217 if(dt > 0){
218 G4double sigma, x, y, z;
219 G4double diffusionCoefficient = G4
220
221 sigma = std::sqrt(2.0 * diffusionC
222
223 x = G4RandGauss::shoot(0., 1.0)*si
224 y = G4RandGauss::shoot(0., 1.0)*si
225 z = G4RandGauss::shoot(0., 1.0)*si
226
227 newPosB = orgPosB + G4ThreeVector(
228 }else if(dt < 0) continue;
229
230 G4double r0 = (newPosB - track->GetP
231 G4double irt = GetIndependentReactio
232
233
234 if(irt>=0 && irt<timeMax - globalTim
235 {
236 irt += globalTime;
237 if(irt < minTime) minTime = irt;
238 #ifdef DEBUG
239 G4cout<<irt<<'\t'<<molConfA->GetNa
240 #endif
241 fReactionSet->AddReaction(irt,trac
242 }
243 }
244 spaceBin.clear();
245 }
246 }
247 }
248
249 // Scavenging & first order reactions
250
251 auto fReactionDatas = fMolReactionTable->Get
252 G4double index = -1;
253 //change the scavenging filter of the IRT be
254 if(timeMax > 1*us)
255 {
256 minTime = timeMax;
257 }
258 //
259
260 for(size_t u=0; u<fReactionDatas->size();u++
261 if((*fReactionDatas)[u]->GetReactant2()->G
262 G4double kObs = (*fReactionDatas)[u]->Ge
263 if(kObs == 0) continue;
264 G4double time = -(std::log(1.0 - G4Unifo
265 if( time < minTime && time >= globalTime
266 minTime = time;
267 index = (G4int)u;
268 }
269 }
270 }
271
272 if(index != -1){
273 #ifdef DEBUG
274 G4cout<<"scavenged: "<<minTime<<'\t'<<molC
275 #endif
276 auto fakeMol = new G4Molecule((*fReaction
277 G4Track* fakeTrack = fakeMol->BuildTrack(g
278 fTrackHolder->Push(fakeTrack);
279 fReactionSet->AddReaction(minTime, track,
280 }
281 }
282
283
284 G4double G4DNAIRT::GetIndependentReactionTime(
285 const auto pMoleculeA = molA;
286 const auto pMoleculeB = molB;
287 auto fReactionData = fMolReactionTable->GetR
288 G4int reactionType = fReactionData->GetReact
289 G4double r0 = distance;
290 if(r0 == 0) r0 += 1e-3*nm;
291 G4double irt = -1 * ps;
292 G4double D = molA->GetDiffusionCoefficient()
293 molB->GetDiffusionCoefficient()
294 if(D == 0) D += 1e-20*(m2/s);
295 G4double rc = fReactionData->GetOnsagerRadiu
296
297 if ( reactionType == 0){
298 G4double sigma = fReactionData->GetEffecti
299
300 if(sigma > r0) return 0; // contact reacti
301 if( rc != 0) r0 = -rc / (1-std::exp(rc/r0)
302
303 G4double Winf = sigma/r0;
304 G4double W = G4UniformRand();
305
306 if ( W > 0 && W < Winf ) irt = (0.25/D) *
307
308 return irt;
309 }
310 if ( reactionType == 1 ){
311 G4double sigma = fReactionData->GetReactio
312 G4double kact = fReactionData->GetActivati
313 G4double kdif = fReactionData->GetDiffusio
314 G4double kobs = fReactionData->GetObserved
315
316 G4double a, b, Winf;
317
318 if ( rc == 0 ) {
319 a = 1/sigma * kact / kobs;
320 b = (r0 - sigma) / 2;
321 } else {
322 G4double v = kact/Avogadro/(4*CLHEP::pi*
323 G4double alpha = v+rc*D/(pow(sigma,2)*(1
324 a = 4*pow(sigma,2)*alpha/(D*pow(rc,2))*p
325 b = rc/4*(cosh(rc/(2*r0))/sinh(rc/(2*r0)
326 r0 = -rc/(1-std::exp(rc/r0));
327 sigma = fReactionData->GetEffectiveReact
328 }
329
330 if(sigma > r0){
331 if(fReactionData->GetProbability() > G4U
332 return irt;
333 }
334 Winf = sigma / r0 * kobs / kdif;
335
336 if(Winf > G4UniformRand()) irt = SamplePDC
337 return irt;
338 }
339
340 return -1 * ps;
341 }
342
343 G4int G4DNAIRT::FindBin(G4int n, G4double xmin
344
345 G4int bin = -1;
346 if ( value <= xmin )
347 bin = 0; //1;
348 else if ( value >= xmax) //!(xmax < value)
349 bin = n-1; //n;
350 else
351 bin = G4int( n * ( value - xmin )/( xmax -
352
353 if ( bin < 0 ) bin = 0;
354 if ( bin >= n ) bin = n-1;
355
356 return bin;
357 }
358
359 G4double G4DNAIRT::SamplePDC(G4double a, G4dou
360
361 G4double p = 2.0 * std::sqrt(2.0*b/a);
362 G4double q = 2.0 / std::sqrt(2.0*b/a);
363 G4double M = max(1.0/(a*a),3.0*b/a);
364
365 G4double X, U, lambdax;
366
367 G4int ntrials = 0;
368 while(true) {
369
370 // Generate X
371 U = G4UniformRand();
372 if ( U < p/(p + q * M) ) X = pow(U * (p +
373 else X = pow(2/((1-U)*(p+q*M)/M),2);
374
375 U = G4UniformRand();
376
377 lambdax = std::exp(-b*b/X) * ( 1.0 - a * s
378
379 if ((X <= 2.0*b/a && U <= lambdax) ||
380 (X >= 2.0*b/a && U*M/X <= lambdax)) br
381
382 ntrials++;
383
384 if ( ntrials > 10000 ){
385 G4cout<<"Totally rejected"<<'\n';
386 return -1.0;
387 }
388 }
389 return X;
390 }
391
392 std::unique_ptr<G4ITReactionChange> G4DNAIRT::
393
394 {
395
396 std::unique_ptr<G4ITReactionChange> pChanges
397 pChanges->Initialize(trackA, trackB);
398
399 const auto pMoleculeA = GetMolecule(trackA)-
400 const auto pMoleculeB = GetMolecule(trackB)-
401 const auto pReactionData = fMolReactionTable
402
403 G4double globalTime = G4Scheduler::Instance(
404 G4double effectiveReactionRadius = pReaction
405
406 const G4double D1 = pMoleculeA->GetDiffusion
407 const G4double D2 = pMoleculeB->GetDiffusion
408
409 G4ThreeVector r1 = trackA.GetPosition();
410 G4ThreeVector r2 = trackB.GetPosition();
411
412 if(r1 == r2) r2 += G4ThreeVector(0,0,1e-3*nm
413
414 G4ThreeVector S1 = r1 - r2;
415
416 G4double r0 = S1.mag();
417
418 S1.setMag(effectiveReactionRadius);
419
420 G4double dt = globalTime - trackA.GetGlobalT
421
422 if(dt != 0 && (D1 + D2) != 0 && r0 != 0){
423 G4double s12 = 2.0 * D1 * dt;
424 G4double s22 = 2.0 * D2 * dt;
425 if(s12 == 0) r2 = r1;
426 else if(s22 == 0) r1 = r2;
427 else{
428 G4double alpha = effectiveReactionRadius
429 G4ThreeVector S2 = (r1 + (s12 / s22)*r2)
430
431
432
433 if(alpha == 0){
434 return pChanges;
435 }
436 S1.setPhi(rad * G4UniformRand() * 2.0 *
437 S1.setTheta(rad * std::acos(1.0 + 1./alp
438
439 r1 = (D1 * S1 + D2 * S2) / (D1 + D2);
440 r2 = D2 * (S2 - S1) / (D1 + D2);
441 }
442 }
443
444 auto pTrackA = const_cast<G4Track*>(pChanges
445 auto pTrackB = const_cast<G4Track*>(pChanges
446
447 pTrackA->SetPosition(r1);
448 pTrackB->SetPosition(r2);
449
450 pTrackA->SetGlobalTime(globalTime);
451 pTrackB->SetGlobalTime(globalTime);
452
453 pTrackA->SetTrackStatus(fStopButAlive);
454 pTrackB->SetTrackStatus(fStopButAlive);
455
456 const G4int nbProducts = pReactionData->GetN
457
458 if(nbProducts != 0){
459
460 const G4double sqrD1 = D1 == 0. ? 0. : std
461 const G4double sqrD2 = D2 == 0. ? 0. : std
462 if((sqrD1 + sqrD2) == 0){
463 return pChanges;
464 }
465 const G4double inv_numerator = 1./(sqrD1 +
466 const G4ThreeVector reactionSite = sqrD2 *
467 + sqrD1 *
468
469 std::vector<G4ThreeVector> position;
470
471 if(nbProducts == 1){
472 position.push_back(reactionSite);
473 }else if(nbProducts == 2){
474 position.push_back(trackA.GetPosition())
475 position.push_back(trackB.GetPosition())
476 }else if (nbProducts == 3){
477 position.push_back(reactionSite);
478 position.push_back(trackA.GetPosition())
479 position.push_back(trackB.GetPosition())
480 }
481
482 for(G4int u = 0; u < nbProducts; u++){
483
484 auto product = new G4Molecule(pReactionD
485 auto productTrack = product->BuildTrack(
486
487
488 productTrack->SetTrackStatus(fAlive);
489 fTrackHolder->Push(productTrack);
490
491 pChanges->AddSecondary(productTrack);
492
493 G4int I = FindBin(fNx, fXMin, fXMax, pos
494 G4int J = FindBin(fNy, fYMin, fYMax, pos
495 G4int K = FindBin(fNz, fZMin, fZMax, pos
496
497 spaceBinned[I][J][K].push_back(productTr
498
499 Sampling(productTrack);
500 }
501 }
502
503 fTrackHolder->MergeSecondariesWithMainList()
504 pChanges->KillParents(true);
505 return pChanges;
506 }
507
508
509 std::vector<std::unique_ptr<G4ITReactionChange
510 G4ITReactionSet* pReactionSet,
511 const G4double /*currentStepTime*/,
512 const G4double fGlobalTime,
513 const G4bool /*reachedUserStepTimeLimit*/)
514 {
515 std::vector<std::unique_ptr<G4ITReactionChan
516 fReactionInfo.clear();
517
518 if (pReactionSet == nullptr)
519 {
520 return fReactionInfo;
521 }
522
523 auto fReactionsetInTime = pReactionSet->GetR
524 assert(fReactionsetInTime.begin() != fReacti
525
526 auto it_begin = fReactionsetInTime.begin();
527 while(it_begin != fReactionsetInTime.end())
528 {
529 G4double irt = it_begin->get()->GetTime();
530
531 if(fGlobalTime < irt) break;
532
533 pReactionSet->SelectThisReaction(*it_begin
534
535 G4Track* pTrackA = it_begin->get()->GetRea
536 G4Track* pTrackB = it_begin->get()->GetRea
537 auto pReactionChange = MakeReaction(*pTrac
538
539 if(pReactionChange){
540 fReactionInfo.push_back(std::move(pReact
541 }
542
543 fReactionsetInTime = pReactionSet->GetReac
544 it_begin = fReactionsetInTime.begin();
545 }
546
547 return fReactionInfo;
548 }
549
550 G4bool G4DNAIRT::TestReactibility(const G4Trac
551 const G4Track& /*trackB*/,
552 G4double /*currentStepTime*/,
553 G4bool /*userStepTimeLimit*/) /*const*/
554 {
555 return true;
556 }
557
558 void G4DNAIRT::SetReactionModel(G4VDNAReaction
559 {
560 fpReactionModel = model;
561 }
562