diff --git a/BEAMS3D/Sources/beams3d_init.f90 b/BEAMS3D/Sources/beams3d_init.f90
index 57d84800..604ed87b 100644
--- a/BEAMS3D/Sources/beams3d_init.f90
+++ b/BEAMS3D/Sources/beams3d_init.f90
@@ -762,10 +762,21 @@ SUBROUTINE beams3d_init
          charge = charge_in(1:nparticles)
          mu_start = mu_start_in(1:nparticles)
          t_end = t_end_in(1:nparticles)
-         beam  = 1
-         nbeams = 1
-         charge_beams(1) = charge_in(1)
-         mass_beams(1)   = mass_in(1)
+         IF (nparticles <= MAXBEAMS) THEN
+            beam(1:nparticles)  = Dex_beams(1:nparticles)
+            charge_beams(1:nparticles) = charge(1:nparticles)
+            mass_beams(1:nparticles) = mass(1:nparticles)
+         ELSE
+            beam(1:MAXBEAMS)  = Dex_beams(1:MAXBEAMS)
+            beam(MAXBEAMS+1:nparticles)  = Dex_beams(MAXBEAMS)
+            charge_beams(1:MAXBEAMS) = charge(1:MAXBEAMS)
+            charge_beams(1:nparticles) = charge(MAXBEAMS)
+            mass_beams(1:MAXBEAMS) = mass(1:MAXBEAMS)
+            mass_beams(1:nparticles) = mass(MAXBEAMS)
+         END IF
+         nbeams = MAXVAL(beam)
+         !charge_beams(1) = charge_in(1)
+         !mass_beams(1)   = mass_in(1)
          lgc2fo_start = .FALSE.
          WHERE ((vr_start == 0) .and. (vphi_start == 0) .and. (vz_start == 0))
             lgc2fo_start = .TRUE.
diff --git a/LIBSTELL/Sources/Modules/beams3d_input_mod.f90 b/LIBSTELL/Sources/Modules/beams3d_input_mod.f90
index 145863f5..943ed0f0 100644
--- a/LIBSTELL/Sources/Modules/beams3d_input_mod.f90
+++ b/LIBSTELL/Sources/Modules/beams3d_input_mod.f90
@@ -394,6 +394,9 @@ SUBROUTINE read_beams3d_input(filename, istat)
             IF (r_start_in(ik) >= 0.0) nparticles = nparticles + 1
          END DO
 
+         ! Assume three is one population if dex_beams is not set.
+         IF (MAXVAL(Dex_beams) < 0) FORALL(ik=1:MAXBEAMS) Dex_beams(ik) = 1
+
 #if !defined(NAG)
       IF (int_type=='NAG') THEN
          int_type = 'LSODE'
@@ -540,6 +543,10 @@ SUBROUTINE write_beams3d_namelist(iunit_out, istat)
          IF (ANY(weight_in /= 1)) WRITE(iunit_out,"(2X,A,1X,'=',10(1X,ES22.12E3))") 'WEIGHT_IN',(weight_in(ik), ik=1,n)
          n = COUNT(t_end_in > -1)
          WRITE(iunit_out,"(2X,A,1X,'=',I6,'*',ES19.12E3)") 'T_END_IN',n,MAXVAL(t_end_in)
+         IF (MAXVAL(dex_beams)>0) THEN
+            n = COUNT(dex_beams>0)
+            WRITE(iunit_out,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'DEX_BEAMS',(dex_beams(ik), ik=1,NION)
+         END IF
       END IF
       WRITE(iunit_out,'(A)') '/'
 
diff --git a/pySTEL/pyproject.toml b/pySTEL/pyproject.toml
index 126dd034..cdbb823e 100644
--- a/pySTEL/pyproject.toml
+++ b/pySTEL/pyproject.toml
@@ -20,7 +20,8 @@ dependencies = [
     "pyvtk",
     "scipy",
     "vtk",
-    "pyyaml"
+    "pyyaml",
+    "h5py"
 ]
 
 [project.urls]