Deuterium fusion to create neutrons in Geant4

I’m trying to simulate the Deuterium - Deuterium interaction D(D,n)3He to create a neutron source using the Geant4 code.

However, from 100k events, where deuterons were fired on a deuterium target, not a single neutron was produced. Could this have something to do with the definition of the deuterium material?

My definition of the deuterium material is

    G4Isotope* D  = new G4Isotope("Deuteron", 1, 2, 2.0141018* CLHEP::g / CLHEP::mole);
    G4Element* elD = new G4Element("Deuterium","elD", 1);
    elD->AddIsotope(D, 1);
    G4Material* matD = new G4Material("matD", 0.00018* CLHEP::g / CLHEP::cm3, 1);
    matD->AddElement(elD, 1);

Is this the correct way to define the deuterium?

From earlier tests with deuterons on lithium I tested the deuteron particle beam and it does produce neutrons by means of the 7Li(D,n)8Be reaction, so I think the beam part of the simulation is working correctly.

Thanks in advance,

Koen

Hello,

it seems that your description of Deuteron is correct. You need to check value of cross section. For that you may take

$G4INSTALL/examples/extended/hadronic/Hadr00

add your extra material and check cross sections.

VI

Hey Koen,

I have been struggling with this issue for the longest time. You need to download the TENDL datafiles and set/create the G4PARTICLEHP environment variable to the dataset. The physics list you need to use is G4HadronPhysicsQGSP_BIC_AllHP and I also needed to use G4IonPhysicsPHP in order to get this to work.

List of nuclear reactions: 
                deuteron + Deuterium Isotope --> neutron + He3:    3850   Q =   3.271 MeV
              deuteron + Deuterium Isotope --> proton + triton:    3545   Q =  4.0359 MeV
                deuteron + Tritium Isotope --> neutron + alpha:  992605   Q =  17.631 MeV

List of generated particles:
            He3:    3850  Emean =  871.83 keV	( 500.31 keV --> 1.2552 MeV)
          alpha:  992605  Emean =   3.619 MeV	( 2.898 MeV --> 4.3413 MeV)
        neutron:  996455  Emean =  14.081 MeV	( 2.1298 MeV --> 14.847 MeV)
         proton:    3545  Emean =  3.0796 MeV	( 2.6627 MeV --> 3.4997 MeV)
         triton:    3545  Emean =  1.0703 MeV	( 650.22 keV --> 1.4871 MeV)

Thanks,

  • Brendan

Brendan,

Can you expand a little bit more about how you solved this problem? I am having the same issue and followed your steps such as downloading the TENDL datafiles but am unable to get any sort of neutron production from deuteron + Deuterium interactions. I have set the various particle environments in the following way:

export G4PROTONHPDATA=/path/to/Geant4/data/G4TENDL1.4/Proton
export G4DEUTERONHPDATA=/path/to/Geant4/data/G4TENDL1.4/Deuteron

and have also added the following Physics Lists:

G4int verb = 2;
RegisterPhysics( new G4HadronPhysicsQGSP_BIC_AllHP(verb));
RegisterPhysics( new G4IonPhysicsPHP(verb));

but after simulating 5M deuterons on a 10micron target of deuterium I am not producing any neutrons.

Thanks!

Let me see if I can help, one thing I noticed was that I set the whole G4PARTICLEHPDATA to G4TENDL1.4. I did not set the proton or deuteron data because I believe G4PARTICLEHPDATA includes those. I also set my system environment variables manually, though I’m not sure that will make a difference.

Here was the Physics List that I used:

  G4int verb = 1;  
  SetVerboseLevel(verb);
  
  // mandatory for G4NuclideTable
  const G4double meanLife = 1*picosecond;  
  G4NuclideTable::GetInstance()->SetMeanLifeThreshold(meanLife);  
    
  // Hadron Elastic scattering
  RegisterPhysics( new G4HadronElasticPhysicsHP(verb));  

  // Hadron Inelastic physics
  RegisterPhysics( new G4HadronPhysicsQGSP_BIC_AllHP(verb));
    
  // Ion Elastic scattering
  RegisterPhysics( new G4IonElasticPhysics(verb));
  
  // Ion Inelastic physics
  RegisterPhysics( new G4IonPhysicsPHP(verb));

  // Gamma physics
  RegisterPhysics( new GammaNuclearPhysics("gamma"));
  
  // stopping Particles
  RegisterPhysics( new G4StoppingPhysics(verb));

  // Radioactive decay
  RegisterPhysics(new G4RadioactiveDecayPhysics());

  //decay
  RegisterPhysics(new G4DecayPhysics());

Hopefully that helps! Good luck.