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 * G4DNAIRT.cc 28 * 29 * Created on: Jul 23, 2019 30 * Author: W. G. Shin 31 * J. Ramos-Mendez and B. Faddegon 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 G4DNAMolecularReactionTable*>(fpReactionTable)), 52 fpReactionModel(nullptr), 53 fTrackHolder(G4ITTrackHolder::Instance()), 54 fReactionSet(nullptr) 55 { 56 timeMin = G4Scheduler::Instance()->GetStartTime(); 57 timeMax = G4Scheduler::Instance()->GetEndTime(); 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 * timeMax); // 95% confidence level 76 77 erfc = new G4ErrorFunction(); 78 } 79 80 81 G4DNAIRT::G4DNAIRT(G4VDNAReactionModel* pReactionModel) 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()->GetStartTime(); 103 timeMax = G4Scheduler::Instance()->GetEndTime(); 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 give any reactions, end the simu. 128 if(fReactionSet->Empty()) 129 { 130 for (auto pTrack : *fTrackHolder->GetMainList()) 131 { 132 pTrack->SetGlobalTime(G4Scheduler::Instance()->GetEndTime()); 133 } 134 } 135 } 136 137 void G4DNAIRT::SpaceBinning(){ 138 auto it_begin = fTrackHolder->GetMainList()->begin(); 139 while(it_begin != fTrackHolder->GetMainList()->end()){ 140 141 G4ThreeVector position = it_begin->GetPosition(); 142 143 if ( fXMin > position.x() ) fXMin = position.x(); 144 if ( fYMin > position.y() ) fYMin = position.y(); 145 if ( fZMin > position.z() ) fZMin = position.z(); 146 147 if ( fXMax < position.x() ) fXMax = position.x(); 148 if ( fYMax < position.y() ) fYMax = position.y(); 149 if ( fZMax < position.z() ) fZMax = position.z(); 150 151 ++it_begin; 152 } 153 154 fNx = G4int((fXMax-fXMin)/fRCutOff) == 0 ? 1 : G4int((fXMax-fXMin)/fRCutOff); 155 fNy = G4int((fYMax-fYMin)/fRCutOff) == 0 ? 1 : G4int((fYMax-fYMin)/fRCutOff); 156 fNz = G4int((fZMax-fZMin)/fRCutOff) == 0 ? 1 : G4int((fZMax-fZMin)/fRCutOff); 157 158 } 159 160 void G4DNAIRT::IRTSampling(){ 161 162 auto it_begin = fTrackHolder->GetMainList()->begin(); 163 while(it_begin != fTrackHolder->GetMainList()->end()){ 164 G4int I = FindBin(fNx, fXMin, fXMax, it_begin->GetPosition().x()); 165 G4int J = FindBin(fNy, fYMin, fYMax, it_begin->GetPosition().y()); 166 G4int K = FindBin(fNz, fZMin, fZMax, it_begin->GetPosition().z()); 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(track); 177 const G4MolecularConfiguration* molConfA = molA->GetMolecularConfiguration(); 178 if(molConfA->GetDiffusionCoefficient() == 0) return; 179 180 const vector<const G4MolecularConfiguration*>* reactivesVector = 181 fMolReactionTable->CanReactWith(molConfA); 182 183 if(reactivesVector == nullptr) return; 184 185 G4double globalTime = G4Scheduler::Instance()->GetGlobalTime(); 186 G4double minTime = timeMax; 187 188 xiniIndex = FindBin(fNx, fXMin, fXMax, track->GetPosition().x()-fRCutOff); 189 xendIndex = FindBin(fNx, fXMin, fXMax, track->GetPosition().x()+fRCutOff); 190 yiniIndex = FindBin(fNy, fYMin, fYMax, track->GetPosition().y()-fRCutOff); 191 yendIndex = FindBin(fNy, fYMin, fYMax, track->GetPosition().y()+fRCutOff); 192 ziniIndex = FindBin(fNz, fZMin, fZMax, track->GetPosition().z()-fRCutOff); 193 zendIndex = FindBin(fNz, fZMin, fZMax, track->GetPosition().z()+fRCutOff); 194 195 for ( int ii = xiniIndex; ii <= xendIndex; ii++ ) { 196 for ( int jj = yiniIndex; jj <= yendIndex; jj++ ) { 197 for ( int kk = ziniIndex; kk <= zendIndex; kk++ ) { 198 199 std::vector<G4Track*> spaceBin = spaceBinned[ii][jj][kk]; 200 for ( int n = 0; n < (int)spaceBinned[ii][jj][kk].size(); n++ ) { 201 if((spaceBin[n] == nullptr) || track == spaceBin[n]) continue; 202 if(spaceBin[n]->GetTrackStatus() == fStopButAlive) continue; 203 204 G4Molecule* molB = G4Molecule::GetMolecule(spaceBin[n]); 205 if(molB == nullptr) continue; 206 207 const G4MolecularConfiguration* molConfB = molB->GetMolecularConfiguration(); 208 if(molConfB->GetDiffusionCoefficient() == 0) continue; 209 210 auto it = std::find(reactivesVector->begin(), reactivesVector->end(), molConfB); 211 if(it == reactivesVector->end()) continue; 212 213 G4ThreeVector orgPosB = spaceBin[n]->GetPosition(); 214 G4double dt = track->GetGlobalTime() - spaceBin[n]->GetGlobalTime(); 215 G4ThreeVector newPosB = orgPosB; 216 217 if(dt > 0){ 218 G4double sigma, x, y, z; 219 G4double diffusionCoefficient = G4Molecule::GetMolecule(spaceBin[n])->GetDiffusionCoefficient(); 220 221 sigma = std::sqrt(2.0 * diffusionCoefficient * dt); 222 223 x = G4RandGauss::shoot(0., 1.0)*sigma; 224 y = G4RandGauss::shoot(0., 1.0)*sigma; 225 z = G4RandGauss::shoot(0., 1.0)*sigma; 226 227 newPosB = orgPosB + G4ThreeVector(x,y,z); 228 }else if(dt < 0) continue; 229 230 G4double r0 = (newPosB - track->GetPosition()).mag(); 231 G4double irt = GetIndependentReactionTime(molConfA, 232 molConfB, 233 r0); 234 if(irt>=0 && irt<timeMax - globalTime) 235 { 236 irt += globalTime; 237 if(irt < minTime) minTime = irt; 238 #ifdef DEBUG 239 G4cout<<irt<<'\t'<<molConfA->GetName()<<" "<<track->GetTrackID()<<'\t'<<molConfB->GetName()<<" "<<spaceBin[n]->GetTrackID()<<'\n'; 240 #endif 241 fReactionSet->AddReaction(irt,track,spaceBin[n]); 242 } 243 } 244 spaceBin.clear(); 245 } 246 } 247 } 248 249 // Scavenging & first order reactions 250 251 auto fReactionDatas = fMolReactionTable->GetReactionData(molConfA); 252 G4double index = -1; 253 //change the scavenging filter of the IRT beyond 1 us proposed by Naoki and Jose 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()->GetDiffusionCoefficient() == 0){ 262 G4double kObs = (*fReactionDatas)[u]->GetObservedReactionRateConstant(); 263 if(kObs == 0) continue; 264 G4double time = -(std::log(1.0 - G4UniformRand())/kObs) + globalTime; 265 if( time < minTime && time >= globalTime && time < timeMax){ 266 minTime = time; 267 index = (G4int)u; 268 } 269 } 270 } 271 272 if(index != -1){ 273 #ifdef DEBUG 274 G4cout<<"scavenged: "<<minTime<<'\t'<<molConfA->GetName()<<it_begin->GetTrackID()<<'\n'; 275 #endif 276 auto fakeMol = new G4Molecule((*fReactionDatas)[index]->GetReactant2()); 277 G4Track* fakeTrack = fakeMol->BuildTrack(globalTime,track->GetPosition()); 278 fTrackHolder->Push(fakeTrack); 279 fReactionSet->AddReaction(minTime, track, fakeTrack); 280 } 281 } 282 283 284 G4double G4DNAIRT::GetIndependentReactionTime(const G4MolecularConfiguration* molA, const G4MolecularConfiguration* molB, G4double distance) { 285 const auto pMoleculeA = molA; 286 const auto pMoleculeB = molB; 287 auto fReactionData = fMolReactionTable->GetReactionData(pMoleculeA, pMoleculeB); 288 G4int reactionType = fReactionData->GetReactionType(); 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->GetOnsagerRadius(); 296 297 if ( reactionType == 0){ 298 G4double sigma = fReactionData->GetEffectiveReactionRadius(); 299 300 if(sigma > r0) return 0; // contact reaction 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) * std::pow( (r0-sigma)/erfc->erfcInv(r0*W/sigma), 2 ); 307 308 return irt; 309 } 310 if ( reactionType == 1 ){ 311 G4double sigma = fReactionData->GetReactionRadius(); 312 G4double kact = fReactionData->GetActivationRateConstant(); 313 G4double kdif = fReactionData->GetDiffusionRateConstant(); 314 G4double kobs = fReactionData->GetObservedReactionRateConstant(); 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*pow(sigma,2) * exp(-rc / sigma)); 323 G4double alpha = v+rc*D/(pow(sigma,2)*(1-exp(-rc/sigma))); 324 a = 4*pow(sigma,2)*alpha/(D*pow(rc,2))*pow(sinh(rc/(2*sigma)),2); 325 b = rc/4*(cosh(rc/(2*r0))/sinh(rc/(2*r0))-cosh(rc/(2*sigma))/sinh(rc/(2*sigma))); 326 r0 = -rc/(1-std::exp(rc/r0)); 327 sigma = fReactionData->GetEffectiveReactionRadius(); 328 } 329 330 if(sigma > r0){ 331 if(fReactionData->GetProbability() > G4UniformRand()) return 0; 332 return irt; 333 } 334 Winf = sigma / r0 * kobs / kdif; 335 336 if(Winf > G4UniformRand()) irt = SamplePDC(a,b)/D; 337 return irt; 338 } 339 340 return -1 * ps; 341 } 342 343 G4int G4DNAIRT::FindBin(G4int n, G4double xmin, G4double xmax, G4double value) { 344 345 G4int bin = -1; 346 if ( value <= xmin ) 347 bin = 0; //1; 348 else if ( value >= xmax) //!(xmax < value) ) //value >= xmax ) 349 bin = n-1; //n; 350 else 351 bin = G4int( n * ( value - xmin )/( xmax - xmin ) ); //bin = 1 + G4int( n * ( value - xmin )/( xmax - xmin ) ); 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, G4double b) { 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 + q * M) / 2, 2); 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 * std::sqrt(CLHEP::pi * X) * erfc->erfcx(b/std::sqrt(X) + a*std::sqrt(X))); 378 379 if ((X <= 2.0*b/a && U <= lambdax) || 380 (X >= 2.0*b/a && U*M/X <= lambdax)) break; 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::MakeReaction(const G4Track& trackA, 393 const G4Track& trackB) 394 { 395 396 std::unique_ptr<G4ITReactionChange> pChanges(new G4ITReactionChange()); 397 pChanges->Initialize(trackA, trackB); 398 399 const auto pMoleculeA = GetMolecule(trackA)->GetMolecularConfiguration(); 400 const auto pMoleculeB = GetMolecule(trackB)->GetMolecularConfiguration(); 401 const auto pReactionData = fMolReactionTable->GetReactionData(pMoleculeA, pMoleculeB); 402 403 G4double globalTime = G4Scheduler::Instance()->GetGlobalTime(); 404 G4double effectiveReactionRadius = pReactionData->GetEffectiveReactionRadius(); 405 406 const G4double D1 = pMoleculeA->GetDiffusionCoefficient(); 407 const G4double D2 = pMoleculeB->GetDiffusionCoefficient(); 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.GetGlobalTime(); 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 * r0 / (2*(D1 + D2)*dt); 429 G4ThreeVector S2 = (r1 + (s12 / s22)*r2) + G4ThreeVector(G4RandGauss::shoot(0, s12 + s22 * s22 / s12), 430 G4RandGauss::shoot(0, s12 + s22 * s22 / s12), 431 G4RandGauss::shoot(0, s12 + s22 * s22 / s12)); 432 433 if(alpha == 0){ 434 return pChanges; 435 } 436 S1.setPhi(rad * G4UniformRand() * 2.0 * CLHEP::pi); 437 S1.setTheta(rad * std::acos(1.0 + 1./alpha * std::log(1.0 - G4UniformRand() * (1 - std::exp(-2.0 * alpha))))); 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->GetTrackA()); 445 auto pTrackB = const_cast<G4Track*>(pChanges->GetTrackB()); 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->GetNbProducts(); 457 458 if(nbProducts != 0){ 459 460 const G4double sqrD1 = D1 == 0. ? 0. : std::sqrt(D1); 461 const G4double sqrD2 = D2 == 0. ? 0. : std::sqrt(D2); 462 if((sqrD1 + sqrD2) == 0){ 463 return pChanges; 464 } 465 const G4double inv_numerator = 1./(sqrD1 + sqrD2); 466 const G4ThreeVector reactionSite = sqrD2 * inv_numerator * trackA.GetPosition() 467 + sqrD1 * inv_numerator * trackB.GetPosition(); 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(pReactionData->GetProduct(u)); 485 auto productTrack = product->BuildTrack(globalTime, 486 position[u]); 487 488 productTrack->SetTrackStatus(fAlive); 489 fTrackHolder->Push(productTrack); 490 491 pChanges->AddSecondary(productTrack); 492 493 G4int I = FindBin(fNx, fXMin, fXMax, position[u].x()); 494 G4int J = FindBin(fNy, fYMin, fYMax, position[u].y()); 495 G4int K = FindBin(fNz, fZMin, fZMax, position[u].z()); 496 497 spaceBinned[I][J][K].push_back(productTrack); 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>> G4DNAIRT::FindReaction( 510 G4ITReactionSet* pReactionSet, 511 const G4double /*currentStepTime*/, 512 const G4double fGlobalTime, 513 const G4bool /*reachedUserStepTimeLimit*/) 514 { 515 std::vector<std::unique_ptr<G4ITReactionChange>> fReactionInfo; 516 fReactionInfo.clear(); 517 518 if (pReactionSet == nullptr) 519 { 520 return fReactionInfo; 521 } 522 523 auto fReactionsetInTime = pReactionSet->GetReactionsPerTime(); 524 assert(fReactionsetInTime.begin() != fReactionsetInTime.end()); 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()->GetReactants().first; 536 G4Track* pTrackB = it_begin->get()->GetReactants().second; 537 auto pReactionChange = MakeReaction(*pTrackA, *pTrackB); 538 539 if(pReactionChange){ 540 fReactionInfo.push_back(std::move(pReactionChange)); 541 } 542 543 fReactionsetInTime = pReactionSet->GetReactionsPerTime(); 544 it_begin = fReactionsetInTime.begin(); 545 } 546 547 return fReactionInfo; 548 } 549 550 G4bool G4DNAIRT::TestReactibility(const G4Track& /*trackA*/, 551 const G4Track& /*trackB*/, 552 G4double /*currentStepTime*/, 553 G4bool /*userStepTimeLimit*/) /*const*/ 554 { 555 return true; 556 } 557 558 void G4DNAIRT::SetReactionModel(G4VDNAReactionModel* model) 559 { 560 fpReactionModel = model; 561 } 562