Geant4 Version: 11.4.0
Operating System: Ubuntu 24.04.3 LTS
Compiler/Version: gcc 13.3.0
CMake Version: 3.28.3
Dear Geant4 users,
I’m looking for a way to move forward with debugging my Geant4-DNA code but I feel stuck, maybe one of you would have an idea to help me.
I am experiencing recurrent crashes during the Geant4-DNA chemical stage. My simulation runs correctly for chemical times below 1 microsecond. However, when I set:
/scheduler/endTime 1 microsecond
in my .in file, the simulation eventually aborts. The crash time is not constant and varies from run to run. I have observed different types of fatal exceptions :
- Molecule counter becoming negative :
*** End of step N°1868 T_i= 13.6486 ns dt= 1.37975 ps T_f= 13.6499 ns eCollisionBetweenTracks
*** End of step N°1869 T_i= 13.6499 ns dt= 1.5679 ps T_f= 13.6515 ns eCollisionBetweenTracks
*** End of step N°1870 T_i= 13.6515 ns dt= 1.03904 ps T_f= 13.6526 ns eCollisionBetweenTracks
At time : 13.654 ns Reaction : H3O^1 (-126) + e_aq^-1 (-52) -> H^0 (-240)
-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : N_INF_0
issued by : G4VUserMoleculeCounter<20MoleculeCounterIndex>::RemoveMolecule
After removal of 1 species of e_aq^-1 the final number at time 13.6536 ns is less than zero and so not valid.
Index was :(e_aq^-1, Vesicle)
Global time is 13.6536 ns
Previous selected time is 0 ps
*** Fatal Exception *** core dump ***
G4Track (0x624d260c1b40) - track ID = -239, parent ID = 0
Particle type : H - creator process : not available
Kinetic energy : 38.5195 meV - Momentum direction : (0.763822,-0.478059,0.43363)
Step length : 9.28572e-297 fm - total energy deposit : 1.54417 eV
Pre-step point : (-0.00530331,-0.0024723,-0.000647882) - Physical volume : NanoCluster (G4_WATER)
- defined by : not available
Post-step point : (-0.00530331,-0.0024723,-0.000647882) - Physical volume : NanoCluster (G4_WATER)
- defined by : e-_G4DNAElectronSolvation - step status : 4
*** Note: Step information might not be properly updated.
-------- EEEE -------- G4Exception-END --------- EEEE -------
*** G4Exception: Aborting execution ***
Abandon (core dumped)
- Time comparison error in the molecule counter :
*** End of step N°2906 T_i= 73.7587 ns dt= 1 ps T_f= 73.7597 ns eInteractionWithMedium
*** End of step N°2907 T_i= 73.7597 ns dt= 1.07453 ps T_f= 73.7607 ns eCollisionBetweenTracks
At time : 73.762 ns Reaction : e_aq^-1 (-74) + H3O^1 (-129) -> H^0 (-234)
-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : G4MoleculeCounterTimeComparer
issued by : G4MoleculeCounterTimeComparer::operator()
Unknown comparison type
*** Fatal Exception *** core dump ***
G4Track (0x5b79075617c8) - track ID = -105, parent ID = -40
Particle type : H3O - creator process : H2O_DNAMolecularDecay, creator model : Undefined
Kinetic energy : 38.5195 meV - Momentum direction : (-0.0180335,-0.308808,0.950953)
Step length : 1.10043e-297 fm - total energy deposit : 7.1964 eV
Pre-step point : (-0.00412302,-0.00121582,0.00339967) - Physical volume : NanoCluster (G4_WATER)
- defined by : e-_G4DNAVibExcitation - step status : 4
Post-step point : (-0.00412302,-0.00121582,0.00339967) - Physical volume : NanoCluster (G4_WATER)
- defined by : e-_G4DNAElectronSolvation - step status : 4
*** Note: Step information might not be properly updated.
-------- EEEE -------- G4Exception-END --------- EEEE -------
*** G4Exception: Aborting execution ***
Abandon (core dumped)
The crashes always occur during the chemical stage, after a reaction. However, the crash time, reaction index, and affected volume are not consistent across runs.
My geometry is entirely made of G4_WATER:
- World : 50 µm cube
- Cell : 7 µm radius sphere
- Vesicle : 0.5 µm radius sphere, randomly placed inside the cell
- NanoCluster : 100 nm radius sphere, placed close to the vesicle boundary with a 2 nm margin
The primary electron source is randomly generated inside the NanoCluster.
This source is quite close to the border between two volumes, could proximity to the geometry boundary cause inconsistencies in the molecule counter, even if all volumes use the same material (G4_WATER)?
I would be grateful for any advice on what to investigate next!
Thank you very much for your help,
Regards,
Sarah
Below my DetectorConstruction.cc & .hh files and my PrimaryGeneratorAction.cc file. I don’t want to burden you with too much of my code, but I can provide additional files if needed.
DetectorConstruction.cc :
G4VPhysicalVolume* DetectorConstruction::Construct()
{
auto man = G4NistManager::Instance();
auto water = man->FindOrBuildMaterial("G4_WATER");
fAGuIX = Materials::CreateAGuIX();
const G4double worldXYZ = 50 * um;
auto solidWorld = new G4Box("World", 0.5 * worldXYZ, 0.5 * worldXYZ, 0.5 * worldXYZ);
auto lvWorld = new G4LogicalVolume(solidWorld, water, "World");
auto pvWorld =
new G4PVPlacement(nullptr, G4ThreeVector(), lvWorld, "World", nullptr, false, 0, true);
const auto worldVisAtt = new G4VisAttributes(G4Colour(.5, 1.0, .5));
worldVisAtt->SetVisibility(true);
worldVisAtt->SetForceWireframe(true);
lvWorld->SetVisAttributes(worldVisAtt);
auto lvCell = ConstructCell(pvWorld);
auto lvVesicle = ConstructVesicle(lvCell);
ConstructNanoCluster(lvVesicle);
return pvWorld;
}
G4LogicalVolume* DetectorConstruction::ConstructCell(G4VPhysicalVolume* pvWorld)
{
const G4double cellRadius = 7 * um;
const G4double nucleusRadius = 4 * um;
const G4int nMitochondria = 100;
const G4double mitoA = 0.55 * micrometer;
const G4double mitoB = 0.25 * micrometer;
const G4double mitoC = 0.90 * micrometer;
auto solidCell = new G4Orb("Cell", cellRadius);
auto lvCell = new G4LogicalVolume(
solidCell, G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER"), "Cell");
auto pvCell =
new G4PVPlacement(nullptr, G4ThreeVector(), "Cell", lvCell, pvWorld, false, 0, true);
auto cellVisAtt = new G4VisAttributes(G4Colour(0.0, 0.0, 1.0, 0.5));
cellVisAtt->SetVisibility(false);
lvCell->SetVisAttributes(cellVisAtt);
auto solidNucleus = new G4Orb("Nucleus", nucleusRadius);
auto lvNucleus = new G4LogicalVolume(
solidNucleus, G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER"), "Nucleus");
new G4PVPlacement(nullptr, G4ThreeVector(), "Nucleus", lvNucleus, pvCell, false, 0, true);
auto solidMito = new G4Ellipsoid("Mitochondria", mitoA, mitoB, mitoC);
auto lvMito = new G4LogicalVolume(
solidMito, G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER"), "Mitochondria");
for (auto i = 0; i < nMitochondria; ++i) {
G4bool overlap = true;
do {
auto u = twopi * G4UniformRand();
auto v = std::acos(2 * G4UniformRand() - 1);
auto dr = G4UniformRand() * cellRadius;
auto x = dr * std::cos(u) * std::sin(v);
auto y = dr * std::sin(u) * std::sin(v);
auto z = dr * std::cos(v);
auto pos = G4ThreeVector(x, y, z);
auto phi = G4UniformRand() * 2 * pi;
auto psi = G4UniformRand() * 2 * pi;
auto rot = new G4RotationMatrix();
rot->rotateX(psi);
rot->rotateY(phi);
auto pvMito = new G4PVPlacement(rot, pos, "Mitochondria", lvMito, pvCell, false, i, false);
overlap = pvMito->CheckOverlaps(1000, 0, false);
if (overlap) {
G4PhysicalVolumeStore::DeRegister(pvMito);
}
} while (overlap);
}
return lvCell;
}
G4LogicalVolume* DetectorConstruction::ConstructVesicle(G4LogicalVolume* lvCell)
{
const G4double vesicleRadius = 0.5 * um;
fVesicleRadius = vesicleRadius;
G4double cellRadius = 7*um;
auto solidVesicle = new G4Orb("Vesicle", vesicleRadius);
auto lvVesicle = new G4LogicalVolume(
solidVesicle,
G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER"),
"Vesicle"
);
auto vesicleVisAtt = new G4VisAttributes(G4Colour(1.0, 1.0, 0.0,0.4));
vesicleVisAtt->SetVisibility(true);
lvVesicle->SetVisAttributes(vesicleVisAtt);
G4bool overlap = true;
G4ThreeVector vesicle_pos;
while (overlap)
{
auto u = twopi * G4UniformRand();
auto v = std::acos(2*G4UniformRand() - 1);
G4double maxR = cellRadius - vesicleRadius;
auto r = maxR * std::cbrt(G4UniformRand());
vesicle_pos = G4ThreeVector(
r * std::cos(u) * std::sin(v),
r * std::sin(u) * std::sin(v),
r * std::cos(v)
);
G4PVPlacement* pvVesicle = new G4PVPlacement(nullptr, vesicle_pos, lvVesicle, "Vesicle", lvCell, false, 0, true);
overlap = pvVesicle->CheckOverlaps(1000, 0, false);
if (overlap)
{
G4PhysicalVolumeStore::DeRegister(pvVesicle);
}
else
{
fVesiclePosition = vesicle_pos;
}
}
return lvVesicle;
}
void DetectorConstruction::ConstructNanoCluster(G4LogicalVolume* lvVesicle)
{
const G4double nanoClusterRadius = 100 * nm;
fNanoClusterRadius = nanoClusterRadius;
auto solidNanoCluster = new G4Orb("NanoCluster", nanoClusterRadius);
auto lvNanoCluster = new G4LogicalVolume(solidNanoCluster, G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER"), "NanoCluster");
auto nanoClusterVisAtt = new G4VisAttributes(G4Colour(1.0, 0, 1.0));
nanoClusterVisAtt->SetForceSolid(true);
lvNanoCluster->SetVisAttributes(nanoClusterVisAtt);
double u = G4UniformRand();
double v = G4UniformRand();
double theta = std::acos(2.0*v - 1.0);
double phi = 2.0 * CLHEP::pi * u;
G4double margin = 2 * nm;
double r = fVesicleRadius - nanoClusterRadius - margin;
G4ThreeVector localPos(
r * std::sin(theta) * std::cos(phi),
r * std::sin(theta) * std::sin(phi),
r * std::cos(theta)
);
new G4PVPlacement(
nullptr,
localPos,
lvNanoCluster,
"NanoCluster",
lvVesicle,
false,
0,
true
);
fNanoClusterLocalPosition = localPos;
}
PrimaryGeneratorAction.cc :
PrimaryGeneratorAction::PrimaryGeneratorAction() : G4VUserPrimaryGeneratorAction()
{
G4int n_particle = 1;
fParticleGun = std::make_unique<G4ParticleGun>(n_particle);
// default particle kinematic
G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
G4ParticleDefinition* particle = particleTable->FindParticle("e-");
fParticleGun->SetParticleDefinition(particle);
fParticleGun->SetParticlePosition(G4ThreeVector(0., 0., 0.));
fParticleGun->SetParticleEnergy(100 * keV);
fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., 1.));
}
void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
auto det = static_cast<const DetectorConstruction*>(
G4RunManager::GetRunManager()->GetUserDetectorConstruction());
G4ThreeVector vesicleCenter = det->GetVesiclePosition();
G4ThreeVector nanoClusterLocal = det->GetNanoClusterLocalPosition();
G4ThreeVector clusterGlobal = vesicleCenter + nanoClusterLocal;
G4double R = det->GetNanoClusterRadius();
G4double u = G4UniformRand();
G4double v = G4UniformRand();
G4double w = G4UniformRand();
G4double r = R * std::cbrt(u);
G4double theta = std::acos(2.0*v - 1.0);
G4double phi = CLHEP::twopi * w;
G4ThreeVector localPos(
r * std::sin(theta) * std::cos(phi),
r * std::sin(theta) * std::sin(phi),
r * std::cos(theta)
);
G4ThreeVector globalPos = clusterGlobal + localPos;
G4cout << "=== Event " << anEvent->GetEventID() << " ===" << G4endl;
G4cout << "Vesicle center: " << vesicleCenter << G4endl;
G4cout << "NanoCluster center: " << clusterGlobal << G4endl;
G4cout << "Source local pos : " << localPos << G4endl;
G4cout << "Source global pos: " << globalPos << G4endl;
static bool firstCall = true;
std::ofstream outFile("SourcePositions.csv", std::ios::app);
if (outFile.is_open()) {
if (firstCall) {
outFile << "EventID,X_um,Y_um,Z_um\n";
firstCall = false;
}
outFile << anEvent->GetEventID() << ","
<< globalPos.x()/um << ","
<< globalPos.y()/um << ","
<< globalPos.z()/um << "\n";
outFile.close();
}
G4double phiDir = CLHEP::twopi * G4UniformRand();
G4double cosT = 2.0 * G4UniformRand() - 1.0;
G4double sinT = std::sqrt(1.0 - cosT*cosT);
G4ThreeVector dir(
sinT * std::cos(phiDir),
sinT * std::sin(phiDir),
cosT
);
fParticleGun->SetParticlePosition(globalPos);
fParticleGun->SetParticleMomentumDirection(dir);
fParticleGun->GeneratePrimaryVertex(anEvent);
}
DetectorConstruction.hh :
class G4VPhysicalVolume;
class DetectorConstruction : public G4VUserDetectorConstruction
{
public:
G4VPhysicalVolume* Construct() override;
const G4ThreeVector& GetVesiclePosition() const { return fVesiclePosition; }
G4double GetVesicleRadius() const { return fVesicleRadius; }
G4double GetNanoClusterRadius() const { return fNanoClusterRadius; }
const G4ThreeVector& GetNanoClusterLocalPosition() const { return fNanoClusterLocalPosition; }
protected:
G4LogicalVolume* ConstructCell(G4VPhysicalVolume*);
G4LogicalVolume* ConstructVesicle(G4LogicalVolume*);
void ConstructNanoCluster(G4LogicalVolume*);
private:
G4ThreeVector fVesiclePosition;
G4double fVesicleRadius;
G4ThreeVector fNanoClusterLocalPosition;
G4double fNanoClusterRadius;
G4Material* fAGuIX = nullptr;
};