Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Types | Public Member Functions | List of all members
G4DNAGillespieDirectMethod Class Reference

#include <G4DNAGillespieDirectMethod.hh>

Public Types

using MolType = const G4MolecularConfiguration*
 
using Index = G4VDNAMesh::Index
 
using Voxel = G4DNAMesh::Voxel
 
using JumpingData = std::pair<MolType, Index>
 
using ReactionData = const G4DNAMolecularReactionData
 
using EventIt = G4DNAEventSet::EventSet::iterator
 

Public Member Functions

 G4DNAGillespieDirectMethod ()
 
 ~G4DNAGillespieDirectMethod ()
 
G4double PropensityFunction (const Voxel &voxel, ReactionData *data)
 
G4double PropensityFunction (const Voxel &voxel, MolType moleType)
 
void SetVoxelMesh (G4DNAMesh &mesh)
 
void SetTimeStep (const G4double &stepTime)
 
void Initialize ()
 
void CreateEvent (const Index &index)
 
void CreateEvents ()
 
void SetEventSet (G4DNAEventSet *)
 
void SetVerbose (const G4int &verbose)
 
G4bool SetEquilibrium (const G4DNAMolecularReactionData *pReaction)
 
void ResetEquilibrium ()
 

Detailed Description

Definition at line 41 of file G4DNAGillespieDirectMethod.hh.

Member Typedef Documentation

◆ EventIt

using G4DNAGillespieDirectMethod::EventIt = G4DNAEventSet::EventSet::iterator

Definition at line 51 of file G4DNAGillespieDirectMethod.hh.

◆ Index

using G4DNAGillespieDirectMethod::Index = G4VDNAMesh::Index

Definition at line 47 of file G4DNAGillespieDirectMethod.hh.

◆ JumpingData

using G4DNAGillespieDirectMethod::JumpingData = std::pair<MolType, Index>

Definition at line 49 of file G4DNAGillespieDirectMethod.hh.

◆ MolType

using G4DNAGillespieDirectMethod::MolType = const G4MolecularConfiguration*

Definition at line 46 of file G4DNAGillespieDirectMethod.hh.

◆ ReactionData

using G4DNAGillespieDirectMethod::ReactionData = const G4DNAMolecularReactionData

Definition at line 50 of file G4DNAGillespieDirectMethod.hh.

◆ Voxel

using G4DNAGillespieDirectMethod::Voxel = G4DNAMesh::Voxel

Definition at line 48 of file G4DNAGillespieDirectMethod.hh.

Constructor & Destructor Documentation

◆ G4DNAGillespieDirectMethod()

G4DNAGillespieDirectMethod::G4DNAGillespieDirectMethod ( )

Definition at line 38 of file G4DNAGillespieDirectMethod.cc.

39 : fMolecularReactions(G4DNAMolecularReactionTable::Instance())
40{}
G4DNAMolecularReactionTable::Instance
static G4DNAMolecularReactionTable * Instance()
Definition G4DNAMolecularReactionTable.cc:360

◆ ~G4DNAGillespieDirectMethod()

G4DNAGillespieDirectMethod::~G4DNAGillespieDirectMethod ( )
default

Member Function Documentation

◆ CreateEvent()

void G4DNAGillespieDirectMethod::CreateEvent ( const Index & index)

Definition at line 193 of file G4DNAGillespieDirectMethod.cc.

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}
FatalErrorInArgument
@ FatalErrorInArgument
Definition G4ExceptionSeverity.hh:68
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition G4Exception.hh:40
G4BestUnit
#define G4BestUnit(a, b)
Definition G4SteppingVerbose.cc:46
G4double
double G4double
Definition G4Types.hh:83
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
G4DNAEventSet::CreateEvent
void CreateEvent(const G4double &time, const Index &index, Event::ReactionData *pReactionData)
Definition G4DNAEventSet.cc:75
G4DNAMesh::GetVoxel
Voxel & GetVoxel(const Index &index)
Definition G4DNAMesh.cc:36

Referenced by CreateEvents().

◆ CreateEvents()

void G4DNAGillespieDirectMethod::CreateEvents ( )

Definition at line 175 of file G4DNAGillespieDirectMethod.cc.

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}
G4DNAMesh::PrintVoxel
void PrintVoxel(const Index &index)
Definition G4DNAMesh.cc:100
G4DNAMesh::begin
auto begin()
Definition G4DNAMesh.hh:60
G4DNAMesh::end
auto end()
Definition G4DNAMesh.hh:59

◆ Initialize()

void G4DNAGillespieDirectMethod::Initialize ( )

Definition at line 157 of file G4DNAGillespieDirectMethod.cc.

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}
G4Scheduler::Instance
static G4Scheduler * Instance()
Definition G4Scheduler.cc:101
G4Scheduler::GetScavengerMaterial
G4VScavengerMaterial * GetScavengerMaterial() const
Definition G4Scheduler.hh:197

◆ PropensityFunction() [1/2]

G4double G4DNAGillespieDirectMethod::PropensityFunction ( const Voxel & voxel,
MolType moleType )

Definition at line 68 of file G4DNAGillespieDirectMethod.cc.

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}

◆ PropensityFunction() [2/2]

G4double G4DNAGillespieDirectMethod::PropensityFunction ( const Voxel & voxel,
ReactionData * data )

Definition at line 96 of file G4DNAGillespieDirectMethod.cc.

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}

◆ ResetEquilibrium()

void G4DNAGillespieDirectMethod::ResetEquilibrium ( )

Definition at line 394 of file G4DNAGillespieDirectMethod.cc.

395{
396 for(auto& it : fEquilibriumProcesses)
397 {
398 it.second->Reset();
399 }
400}

◆ SetEquilibrium()

G4bool G4DNAGillespieDirectMethod::SetEquilibrium ( const G4DNAMolecularReactionData * pReaction)

Definition at line 383 of file G4DNAGillespieDirectMethod.cc.

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}

◆ SetEventSet()

void G4DNAGillespieDirectMethod::SetEventSet ( G4DNAEventSet * pEventSet)

Definition at line 44 of file G4DNAGillespieDirectMethod.cc.

45{
46 fpEventSet = pEventSet;
47}

◆ SetTimeStep()

void G4DNAGillespieDirectMethod::SetTimeStep ( const G4double & stepTime)

Definition at line 189 of file G4DNAGillespieDirectMethod.cc.

190{
191 fTimeStep = stepTime;
192}

◆ SetVerbose()

void G4DNAGillespieDirectMethod::SetVerbose ( const G4int & verbose)
inline

Definition at line 61 of file G4DNAGillespieDirectMethod.hh.

62 {
63 fVerbose = verbose;
64 }

◆ SetVoxelMesh()

void G4DNAGillespieDirectMethod::SetVoxelMesh ( G4DNAMesh & mesh)
inline

Definition at line 55 of file G4DNAGillespieDirectMethod.hh.

55{ fpMesh = &mesh; }
The documentation for this class was generated from the following files: