G4DNAGillespieDirectMethod Class Reference

#include <G4DNAGillespieDirectMethod.hh>

Public Types
using	EventIt = G4DNAEventSet::EventSet::iterator
using	Index = G4Voxel::Index
using	JumpingData = std::pair< MolType, Index >
using	Key = unsigned int
using	MolType = const G4MolecularConfiguration *
using	ReactionData = const G4DNAMolecularReactionData

Public Member Functions
G4double	ComputeNumberInNode (const Index &index, MolType type)
void	CreateEvent (unsigned int key)
G4double	DiffusiveJumping (const Index &index)
	G4DNAGillespieDirectMethod ()
void	Initialize ()
G4double	PropensityFunction (const Index &index, MolType moleType)
G4double	PropensityFunction (const Index &index, ReactionData *data)
G4double	Reaction (const Index &index)
void	SetEventSet (G4DNAEventSet *)
void	SetTimeStep (const G4double &stepTime)
void	SetVoxelMesh (G4DNAMesh &mesh)
G4double	VolumeOfNode (const Index &index)
	~G4DNAGillespieDirectMethod ()

Private Member Functions
G4bool	FindScavenging (const Index &index, MolType, G4double &)

Private Attributes
std::map< G4double, JumpingData >	fJumpingDataMap
G4DNAMolecularReactionTable *	fMolecularReactions
G4DNAEventSet *	fpEventSet
G4DNAMesh *	fpMesh
G4DNAScavengerMaterial *	fpScavengerMaterial
std::map< G4double, ReactionData * >	fReactionDataMap
G4double	fTimeStep
G4double	fVerbose

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 = G4Voxel::Index

Definition at line 48 of file G4DNAGillespieDirectMethod.hh.

JumpingData

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

Definition at line 49 of file G4DNAGillespieDirectMethod.hh.

Key

using G4DNAGillespieDirectMethod::Key = unsigned int

Definition at line 47 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.

Constructor & Destructor Documentation

G4DNAGillespieDirectMethod()

G4DNAGillespieDirectMethod::G4DNAGillespieDirectMethod

(

)

Definition at line 38 of file G4DNAGillespieDirectMethod.cc.

  : fMolecularReactions(G4DNAMolecularReactionTable::Instance())
  , fpMesh(nullptr)
  , fTimeStep(0)
  , fpEventSet(nullptr)
  , fVerbose(0)
  , fpScavengerMaterial(nullptr)
{}

~G4DNAGillespieDirectMethod()

G4DNAGillespieDirectMethod::~G4DNAGillespieDirectMethod ( )

default

Member Function Documentation

ComputeNumberInNode()

G4double G4DNAGillespieDirectMethod::ComputeNumberInNode	(	const Index &	index,
		MolType	type
	)

Definition at line 288 of file G4DNAGillespieDirectMethod.cc.

{
  if(type->GetDiffusionCoefficient() != 0)
  {
    const auto& node = fpMesh->GetVoxelMapList(index);
    const auto& it   = node.find(type);
    return (it != node.end()) ? (it->second) : 0;
  }
  else
  {
    return 0;
  }
}

References fpMesh, G4MolecularConfiguration::GetDiffusionCoefficient(), and G4DNAMesh::GetVoxelMapList().

Referenced by DiffusiveJumping(), and PropensityFunction().

CreateEvent()

void G4DNAGillespieDirectMethod::CreateEvent

(

unsigned int

key

)

Definition at line 176 of file G4DNAGillespieDirectMethod.cc.

{
  G4double r1         = G4UniformRand();
  G4double r2         = G4UniformRand();
  auto index          = fpMesh->GetIndex(key);
  G4double dAlpha0    = DiffusiveJumping(index);
  G4double rAlpha0    = Reaction(index);
  G4double alphaTotal = dAlpha0 + rAlpha0;
 
  if(alphaTotal == 0)
  {
    return;
  }
  auto timeStep = ((1.0 / (alphaTotal)) * std::log(1.0 / r1)) + fTimeStep;
 
#ifdef DEBUG
  G4cout << "r2 : " << r2 << " rAlpha0 : " << rAlpha0
         << " dAlpha0 : " << dAlpha0 << "    rAlpha0 / (dAlpha0 + rAlpha0) : "
         << rAlpha0 / (dAlpha0 + rAlpha0) << G4endl;
#endif
  if(r2 < rAlpha0 / alphaTotal)
  {
    if(fVerbose > 1)
    {
      G4cout << "=>>>>reaction at : " << timeStep << " timeStep : "
             << G4BestUnit(((1.0 / alphaTotal) * std::log(1.0 / r1)), "Time")
             << G4endl;
    }
    auto rSelectedIter = fReactionDataMap.upper_bound(r2 * alphaTotal);
    fpEventSet->CreateEvent(timeStep, key, rSelectedIter->second);
  }
  else if(dAlpha0 > 0)
  {
    if(fVerbose > 1)
    {
      G4cout << "=>>>>jumping at : " << timeStep << " timeStep : "
             << G4BestUnit(((1.0 / alphaTotal) * std::log(1.0 / r1)), "Time")
             << G4endl;
    }
 
    auto dSelectedIter = fJumpingDataMap.upper_bound(r2 * alphaTotal - rAlpha0);
    auto pDSelected =
      std::make_unique<std::pair<MolType, Index>>(dSelectedIter->second);
    fpEventSet->CreateEvent(timeStep, key, std::move(pDSelected));
  }
#ifdef DEBUG
  G4cout << G4endl;
#endif
}

References G4DNAEventSet::CreateEvent(), DiffusiveJumping(), fJumpingDataMap, fpEventSet, fpMesh, fReactionDataMap, fTimeStep, fVerbose, G4BestUnit, G4cout, G4endl, G4UniformRand, G4DNAMesh::GetIndex(), and Reaction().

Referenced by Initialize().

DiffusiveJumping()

G4double G4DNAGillespieDirectMethod::DiffusiveJumping

(

const Index &

index

)

Definition at line 252 of file G4DNAGillespieDirectMethod.cc.

{
  fJumpingDataMap.clear();
  G4double alpha0        = 0;
  auto NeighboringVoxels = fpMesh->FindNeighboringVoxels(index);
  if(NeighboringVoxels.empty())
  {
    return 0;
  }
  auto iter = G4MoleculeTable::Instance()->GetConfigurationIterator();
  while(iter())
  {
    const auto conf = iter.value();
    auto propensity = PropensityFunction(index, conf);
    if(propensity == 0)
    {
      continue;
    }
    for(const auto& it_Neighbor : NeighboringVoxels)
    {
      alpha0 += propensity;
      fJumpingDataMap[alpha0] = std::make_pair(conf, it_Neighbor);
#ifdef DEBUG
      G4cout << "mole : " << conf->GetName()
             << " number : " << ComputeNumberInNode(index, conf)
             << " propensity : " << propensity << "  alpha0 : " << alpha0
             << G4endl;
#endif
    }
  }
#ifdef DEBUG
  G4cout << "DiffusiveJumping :alpha0 :  " << alpha0 << G4endl;
#endif
  return alpha0;
}

References ComputeNumberInNode(), G4DNAMesh::FindNeighboringVoxels(), fJumpingDataMap, fpMesh, G4cout, G4endl, G4MoleculeTable::GetConfigurationIterator(), G4MoleculeTable::Instance(), PropensityFunction(), and G4MoleculeIterator< MOLECULE >::value().

Referenced by CreateEvent().

FindScavenging()

G4bool G4DNAGillespieDirectMethod::FindScavenging ( const Index &  index, MolType  moletype, G4double &  numberOfScavenger  )

private

Definition at line 303 of file G4DNAGillespieDirectMethod.cc.

{
  numberOfScavenger = 0;
  if(fpScavengerMaterial == nullptr)
  {
    return false;
  }
  auto volumeOfNode = VolumeOfNode(index);
  if(G4MoleculeTable::Instance()->GetConfiguration("H2O") == moletype)
  {
    auto factor       = Avogadro * volumeOfNode;
    numberOfScavenger = factor;
    return true;
  }
 
  G4double totalNumber =
    fpScavengerMaterial->GetNumberMoleculePerVolumeUnitForMaterialConf(
      moletype);
  if(totalNumber == 0)
  {
    return false;
  }
  else
  {
    G4double numberInDouble =
      volumeOfNode * std::floor(totalNumber) / fpMesh->GetBoundingBox().Volume();
    auto numberInInterg = (int) (std::floor(numberInDouble));
    G4double ram          = G4UniformRand();
    G4double change       = numberInDouble - numberInInterg;
    if(ram > change)
    {
      numberOfScavenger = numberInInterg;
    }
    else
    {
      numberOfScavenger = numberInInterg + 1;
    }
    return true;
  }
}

References source.hepunit::Avogadro, fpMesh, fpScavengerMaterial, G4UniformRand, G4DNAMesh::GetBoundingBox(), G4DNAScavengerMaterial::GetNumberMoleculePerVolumeUnitForMaterialConf(), G4MoleculeTable::Instance(), G4DNABoundingBox::Volume(), and VolumeOfNode().

Referenced by PropensityFunction().

Initialize()

void G4DNAGillespieDirectMethod::Initialize

(

)

Definition at line 154 of file G4DNAGillespieDirectMethod.cc.

{
  // for Scavenger
  fpScavengerMaterial = dynamic_cast<G4DNAScavengerMaterial*>(
    G4Scheduler::Instance()->GetScavengerMaterial());
 
  auto begin = fpMesh->begin();
  auto end   = fpMesh->end();
  for(; begin != end; begin++)
  {
    auto key = begin->first;
#ifdef DEBUG
    fpMesh->PrintVoxel(fpMesh->GetIndex(key));
#endif
    CreateEvent(key);
  }
}

References G4DNAMesh::begin(), CreateEvent(), G4DNAMesh::end(), fpMesh, fpScavengerMaterial, G4DNAMesh::GetIndex(), G4Scheduler::GetScavengerMaterial(), G4Scheduler::Instance(), and G4DNAMesh::PrintVoxel().

PropensityFunction() [1/2]

G4double G4DNAGillespieDirectMethod::PropensityFunction	(	const Index &	index,
		MolType	moleType
	)

Definition at line 68 of file G4DNAGillespieDirectMethod.cc.

{
  if(moleType->GetDiffusionCoefficient() == 0)
  {
    return 0.;
  }
  const auto& node = fpMesh->GetVoxelMapList(index);
  G4double alpha   = 0;
  auto it          = node.find(moleType);
  if(it != node.end())
  {
    auto LengthY = fpMesh->GetBoundingBox(index).Getyhi() -
                   fpMesh->GetBoundingBox(index).Getylo();
    G4double d = it->first->GetDiffusionCoefficient() / std::pow(LengthY, 2);
    alpha      = d * it->second;
 
#ifdef DEBUG
    G4cout << it->first->GetName() << " " << it->second
           << " D : " << it->first->GetDiffusionCoefficient()
           << " LengthY : " << LengthY << " PropensityFunction : " << alpha
           << G4endl;
#endif
  }
  return alpha;
}

References alpha, fpMesh, G4cout, G4endl, G4DNAMesh::GetBoundingBox(), G4MolecularConfiguration::GetDiffusionCoefficient(), G4DNAMesh::GetVoxelMapList(), G4DNABoundingBox::Getyhi(), and G4DNABoundingBox::Getylo().

PropensityFunction() [2/2]

G4double G4DNAGillespieDirectMethod::PropensityFunction	(	const Index &	index,
		ReactionData *	data
	)

Definition at line 95 of file G4DNAGillespieDirectMethod.cc.

{
  G4double value;
  auto ConfA               = data->GetReactant1();
  auto ConfB               = data->GetReactant2();
  G4double scavengerNumber = 0;
  auto typeANumber         = FindScavenging(index, ConfA, scavengerNumber)
                       ? scavengerNumber
                       : ComputeNumberInNode(index, ConfA);
 
  auto typeBNumber = FindScavenging(index, ConfB, scavengerNumber)
                       ? scavengerNumber
                       : ComputeNumberInNode(index, ConfB);
 
  if(typeANumber == 0 || typeBNumber == 0)
  {
    return 0;
  }
 
  auto k =
    data->GetObservedReactionRateConstant() / (Avogadro * VolumeOfNode(index));
  if(ConfA == ConfB)
  {
    value = typeANumber * (typeBNumber - 1) * k;
  }
  else
  {
    value = typeANumber * typeBNumber * k;
  }
 
  if(value < 0)
  {
    G4cout << "G4DNAGillespieDirectMethod::PropensityFunction for : "
           << ConfA->GetName() << "(" << typeANumber << ") + "
           << ConfB->GetName() << "(" << typeBNumber
           << ") : propensity : " << value
           << " GetObservedReactionRateConstant : "
           << data->GetObservedReactionRateConstant()
           << " GetEffectiveReactionRadius : "
           << G4BestUnit(data->GetEffectiveReactionRadius(), "Length")
           << " k : " << k << " volume : " << VolumeOfNode(index)
           << " Index : " << index << G4endl;
    assert(false);
  }
 
#ifdef DEBUG
  if(value > 0)
    G4cout << "G4DNAGillespieDirectMethod::PropensityFunction for : "
           << ConfA->GetName() << "(" << typeANumber << ") + "
           << ConfB->GetName() << "(" << typeBNumber
           << ") : propensity : " << value
           << "  Time to Reaction : " << G4BestUnit(timeToReaction, "Time")
           << " k : " << k << " Index : " << index << G4endl;
#endif
 
  return value;
}

References source.hepunit::Avogadro, ComputeNumberInNode(), FindScavenging(), G4BestUnit, G4cout, G4endl, G4DNAMolecularReactionData::GetEffectiveReactionRadius(), G4DNAMolecularReactionData::GetObservedReactionRateConstant(), G4DNAMolecularReactionData::GetReactant1(), G4DNAMolecularReactionData::GetReactant2(), and VolumeOfNode().

Referenced by DiffusiveJumping(), and Reaction().

Reaction()

G4double G4DNAGillespieDirectMethod::Reaction

(

const Index &

index

)

Definition at line 226 of file G4DNAGillespieDirectMethod.cc.

{
  fReactionDataMap.clear();
  G4double alpha0 = 0;
  auto dataList   = fMolecularReactions->GetVectorOfReactionData();
  if(dataList.empty())
  {
    G4cout << "MolecularReactionTable empty" << G4endl;
    assert(false);
  }
  for(const auto& it : dataList)
  {
    auto propensity = PropensityFunction(index, it);
    if(propensity == 0)
    {
      continue;
    }
    alpha0 += propensity;
    fReactionDataMap[alpha0] = it;
  }
#ifdef DEBUG
  G4cout << "Reaction :alpha0 :  " << alpha0 << G4endl;
#endif
  return alpha0;
}

References fMolecularReactions, fReactionDataMap, G4cout, G4endl, G4DNAMolecularReactionTable::GetVectorOfReactionData(), and PropensityFunction().

Referenced by CreateEvent().

SetEventSet()

void G4DNAGillespieDirectMethod::SetEventSet

(

G4DNAEventSet *

pEventSet

)

Definition at line 49 of file G4DNAGillespieDirectMethod.cc.

{
  fpEventSet = pEventSet;
}

References fpEventSet.

SetTimeStep()

void G4DNAGillespieDirectMethod::SetTimeStep

(

const G4double &

stepTime

)

Definition at line 172 of file G4DNAGillespieDirectMethod.cc.

{
  fTimeStep = stepTime;
}

References fTimeStep.

SetVoxelMesh()

void G4DNAGillespieDirectMethod::SetVoxelMesh ( G4DNAMesh &  mesh )

inline

Definition at line 55 of file G4DNAGillespieDirectMethod.hh.

{ fpMesh = &mesh; }

References fpMesh.

VolumeOfNode()

G4double G4DNAGillespieDirectMethod::VolumeOfNode

(

const Index &

index

)

Definition at line 56 of file G4DNAGillespieDirectMethod.cc.

{
  auto LengthY = fpMesh->GetBoundingBox(index).Getyhi() -
                 fpMesh->GetBoundingBox(index).Getylo();
  auto LengthX = fpMesh->GetBoundingBox(index).Getxhi() -
                 fpMesh->GetBoundingBox(index).Getxlo();
  auto LengthZ = fpMesh->GetBoundingBox(index).Getzhi() -
                 fpMesh->GetBoundingBox(index).Getzlo();
  G4double V = LengthY * LengthX * LengthZ;
  assert(V > 0);
  return V;
}

References fpMesh, G4DNAMesh::GetBoundingBox(), G4DNABoundingBox::Getxhi(), G4DNABoundingBox::Getxlo(), G4DNABoundingBox::Getyhi(), G4DNABoundingBox::Getylo(), G4DNABoundingBox::Getzhi(), and G4DNABoundingBox::Getzlo().

Referenced by FindScavenging(), and PropensityFunction().

Field Documentation

fJumpingDataMap

std::map<G4double , JumpingData> G4DNAGillespieDirectMethod::fJumpingDataMap

private

Definition at line 72 of file G4DNAGillespieDirectMethod.hh.

Referenced by CreateEvent(), and DiffusiveJumping().

fMolecularReactions

G4DNAMolecularReactionTable* G4DNAGillespieDirectMethod::fMolecularReactions

private

Definition at line 66 of file G4DNAGillespieDirectMethod.hh.

Referenced by Reaction().

fpEventSet

G4DNAEventSet* G4DNAGillespieDirectMethod::fpEventSet

private

Definition at line 69 of file G4DNAGillespieDirectMethod.hh.

Referenced by CreateEvent(), and SetEventSet().

fpMesh

G4DNAMesh* G4DNAGillespieDirectMethod::fpMesh

private

Definition at line 67 of file G4DNAGillespieDirectMethod.hh.

Referenced by ComputeNumberInNode(), CreateEvent(), DiffusiveJumping(), FindScavenging(), Initialize(), PropensityFunction(), SetVoxelMesh(), and VolumeOfNode().

fpScavengerMaterial

G4DNAScavengerMaterial* G4DNAGillespieDirectMethod::fpScavengerMaterial

private

Definition at line 74 of file G4DNAGillespieDirectMethod.hh.

Referenced by FindScavenging(), and Initialize().

fReactionDataMap

std::map<G4double , ReactionData*> G4DNAGillespieDirectMethod::fReactionDataMap

private

Definition at line 71 of file G4DNAGillespieDirectMethod.hh.

Referenced by CreateEvent(), and Reaction().

fTimeStep

G4double G4DNAGillespieDirectMethod::fTimeStep

private

Definition at line 68 of file G4DNAGillespieDirectMethod.hh.

Referenced by CreateEvent(), and SetTimeStep().

fVerbose

G4double G4DNAGillespieDirectMethod::fVerbose

private

Definition at line 70 of file G4DNAGillespieDirectMethod.hh.

Referenced by CreateEvent().

---

The documentation for this class was generated from the following files:

• source/processes/electromagnetic/dna/models/include/G4DNAGillespieDirectMethod.hh
• source/processes/electromagnetic/dna/models/src/G4DNAGillespieDirectMethod.cc