diff --git a/Simulation/G4Atlas/G4AtlasAlg/src/G4AtlasAlg.cxx b/Simulation/G4Atlas/G4AtlasAlg/src/G4AtlasAlg.cxx
index d477355a63a0c8116e9507ccec7489d08fb1f5c8..3553af7ee50c66f13f92d7de06f07384840b8c49 100644
--- a/Simulation/G4Atlas/G4AtlasAlg/src/G4AtlasAlg.cxx
+++ b/Simulation/G4Atlas/G4AtlasAlg/src/G4AtlasAlg.cxx
@@ -25,6 +25,9 @@
 #include "G4StackManager.hh"
 #include "G4UImanager.hh"
 #include "G4ScoringManager.hh"
+#include "G4VUserPhysicsList.hh"
+#include "G4VModularPhysicsList.hh"
+#include "G4ParallelWorldPhysics.hh"
 
 // CLHEP includes
 #include "CLHEP/Random/RandomEngine.h"
@@ -61,6 +64,7 @@ G4AtlasAlg::G4AtlasAlg(const std::string& name, ISvcLocator* pSvcLocator)
   declareProperty("G4Commands", m_g4commands, "Commands to send to the G4UI");
   // Multi-threading specific settings
   declareProperty("MultiThreading", m_useMT, "Multi-threading specific settings");
+  declareProperty("ActivateParallelWorlds",m_activateParallelGeometries,"Toggle on/off the G4 parallel geometry system");
 
   // Verbosities
   declareProperty("Verbosities", m_verbosities);
@@ -88,9 +92,6 @@ StatusCode G4AtlasAlg::initialize()
   ATH_CHECK( m_rndmGenSvc.retrieve() );
   ATH_CHECK( m_userActionSvc.retrieve() );
 
-  // FIXME TOO EARLY???
-  ATH_CHECK(m_g4atlasSvc.retrieve());
-
   ATH_CHECK(m_senDetTool.retrieve());
   ATH_CHECK(m_fastSimTool.retrieve());
 
@@ -179,10 +180,39 @@ void G4AtlasAlg::initializeOnce()
     commandLog(returnCode, g4command);
   }
 
-  // G4 init moved to PyG4AtlasAlg / G4AtlasEngine
-  /// @todo Reinstate or delete?! This can't actually be called from the Py algs
-  //ATH_MSG_INFO("Firing initialization of G4!!!");
-  //initializeG4();
+  // Code from G4AtlasSvc
+  auto* rm = G4RunManager::GetRunManager();
+  if(!rm) {
+    throw std::runtime_error("Run manager retrieval has failed");
+  }
+  rm->Initialize();     // Initialization differs slightly in multi-threading.
+  // TODO: add more details about why this is here.
+  if(!m_useMT && rm->ConfirmBeamOnCondition()) {
+    rm->RunInitialization();
+  }
+
+  ATH_MSG_INFO( "retireving the Detector Geometry Service" );
+  if(m_detGeoSvc.retrieve().isFailure()) {
+    throw std::runtime_error("Could not initialize ATLAS DetectorGeometrySvc!");
+  }
+
+  if(m_userLimitsSvc.retrieve().isFailure()) {
+    throw std::runtime_error("Could not initialize ATLAS UserLimitsSvc!");
+  }
+
+  if (m_activateParallelGeometries) {
+    G4VModularPhysicsList* thePhysicsList=dynamic_cast<G4VModularPhysicsList*>(m_physListSvc->GetPhysicsList());
+    if (!thePhysicsList) {
+      throw std::runtime_error("Failed dynamic_cast!! this is not a G4VModularPhysicsList!");
+    }
+#if G4VERSION_NUMBER >= 1010
+    std::vector<std::string>& parallelWorldNames=m_detGeoSvc->GetParallelWorldNames();
+    for (auto& it: parallelWorldNames) {
+      thePhysicsList->RegisterPhysics(new G4ParallelWorldPhysics(it,true));
+    }
+#endif
+  }
+
   return;
 }
 
diff --git a/Simulation/G4Atlas/G4AtlasAlg/src/G4AtlasAlg.h b/Simulation/G4Atlas/G4AtlasAlg/src/G4AtlasAlg.h
index 416ea22cd7f8d1e8b09bca46cb38e46689152a79..6275e7668450816c06943e318cd1f87a72e93dc3 100644
--- a/Simulation/G4Atlas/G4AtlasAlg/src/G4AtlasAlg.h
+++ b/Simulation/G4Atlas/G4AtlasAlg/src/G4AtlasAlg.h
@@ -25,7 +25,7 @@
 #include "G4AtlasInterfaces/ISensitiveDetectorMasterTool.h"
 #include "G4AtlasInterfaces/IFastSimulationMasterTool.h"
 #include "G4AtlasInterfaces/IPhysicsListSvc.h"
-#include "G4AtlasInterfaces/IG4AtlasSvc.h"
+#include "G4AtlasInterfaces/IUserLimitsSvc.h"
 #include "GeneratorObjects/McEventCollection.h"
 
 // ISF includes
@@ -119,10 +119,11 @@ private:
   std::vector<std::string> m_g4commands;
   /// Activate multi-threading configuration
   bool m_useMT{false};
+  bool m_activateParallelGeometries{false};
   /// Random number service
   ServiceHandle<IAthRNGSvc> m_rndmGenSvc{this, "AtRndmGenSvc", "AthRNGSvc", ""}; // TODO rename property
-  /// G4Atlas Service - handles G4 initialization
-  ServiceHandle<IG4AtlasSvc> m_g4atlasSvc{this, "G4AtlasSvc", "G4AtlasSvc", ""};
+  ///
+  ServiceHandle<IUserLimitsSvc> m_userLimitsSvc{this, "UserLimitsSvc", "UserLimitsSvc", ""};
   /// User Action Service
   ServiceHandle<G4UA::IUserActionSvc> m_userActionSvc{this, "UserActionSvc", "G4UA::UserActionSvc", ""};
   /// Detector Geometry Service (builds G4 Geometry)
diff --git a/Simulation/G4Atlas/G4AtlasApps/test/test_AtlasG4_cosmics_configuration.py b/Simulation/G4Atlas/G4AtlasApps/test/test_AtlasG4_cosmics_configuration.py
index 55260b5f47ef60be690bfd5ab6b60a19204a0c0c..37484ce4240568633e0704a1a20346b30ba8d4ff 100755
--- a/Simulation/G4Atlas/G4AtlasApps/test/test_AtlasG4_cosmics_configuration.py
+++ b/Simulation/G4Atlas/G4AtlasApps/test/test_AtlasG4_cosmics_configuration.py
@@ -238,7 +238,7 @@ class TestAtlasG4Cosmics(unittest.TestCase):
 
 
     def test___G4AtlasAlg_ListOfSetProperties(self):
-        expected_list = ['AtRndmGenSvc', 'DetGeoSvc', 'DetStore', 'EvtStore', 'ExtraInputs', 'ExtraOutputs', 'FastSimMasterTool', 'FlagAbortedEvents', 'G4AtlasSvc', 'G4Commands', 'GeoIDSvc', 'InputConverter', 'InputTruthCollection', 'KillAbortedEvents', 'MultiThreading', 'NeededResources', 'OutputTruthCollection', 'PhysicsListSvc', 'RandomGenerator', 'RecordFlux', 'ReleaseGeoModel', 'SenDetMasterTool', 'TruthRecordService', 'UserActionSvc', 'Verbosities']
+        expected_list = ['AtRndmGenSvc', 'DetGeoSvc', 'DetStore', 'EvtStore', 'ExtraInputs', 'ExtraOutputs', 'FastSimMasterTool', 'FlagAbortedEvents', 'G4Commands', 'GeoIDSvc', 'InputConverter', 'InputTruthCollection', 'KillAbortedEvents', 'MultiThreading', 'NeededResources', 'OutputTruthCollection', 'PhysicsListSvc', 'RandomGenerator', 'RecordFlux', 'ReleaseGeoModel', 'SenDetMasterTool', 'TruthRecordService', 'UserActionSvc', 'UserLimitsSvc', 'Verbosities']
         g4atlasalg = self._job_config_dict['G4AtlasAlg']
         actual_list = g4atlasalg.keys()
         expected_property_value_sorted = sorted(expected_list)
@@ -282,9 +282,9 @@ class TestAtlasG4Cosmics(unittest.TestCase):
         self._assert_Algorithm_property_equal('G4AtlasAlg', 'SenDetMasterTool', expected_tool_name)
 
 
-    def test___G4AtlasAlg_G4AtlasSvc_setCorrectly(self):
-        expected_service_name = 'G4AtlasSvc'
-        self._assert_Algorithm_property_equal('G4AtlasAlg', 'G4AtlasSvc', expected_service_name)
+    def test___G4AtlasAlg_UserLimitsSvc_setCorrectly(self):
+        expected_service_name = 'UserLimitsSvc'
+        self._assert_Algorithm_property_equal('G4AtlasAlg', 'UserLimitsSvc', expected_service_name)
 
 
     def test___G4AtlasAlg_UserActionSvc_setCorrectly(self):
diff --git a/Simulation/G4Atlas/G4AtlasApps/test/test_AtlasG4_tf_configuration.py b/Simulation/G4Atlas/G4AtlasApps/test/test_AtlasG4_tf_configuration.py
index 6fc74c7b3c67af02231c506473f4dd8e4b9f296a..70be030e872a2c94a41f25a6e4f86e192b2399fe 100755
--- a/Simulation/G4Atlas/G4AtlasApps/test/test_AtlasG4_tf_configuration.py
+++ b/Simulation/G4Atlas/G4AtlasApps/test/test_AtlasG4_tf_configuration.py
@@ -154,7 +154,7 @@ class TestAtlasG4(unittest.TestCase):
 
 
     def test___G4AtlasAlg_ListOfSetProperties(self):
-        expected_list = ['AtRndmGenSvc', 'DetGeoSvc', 'DetStore', 'EvtStore', 'ExtraInputs', 'ExtraOutputs', 'FastSimMasterTool', 'FlagAbortedEvents', 'G4AtlasSvc', 'G4Commands', 'GeoIDSvc', 'InputConverter', 'InputTruthCollection', 'KillAbortedEvents', 'MultiThreading', 'NeededResources', 'OutputTruthCollection', 'PhysicsListSvc', 'RandomGenerator', 'RecordFlux', 'ReleaseGeoModel', 'SenDetMasterTool', 'TruthRecordService', 'UserActionSvc', 'Verbosities']
+        expected_list = ['AtRndmGenSvc', 'DetGeoSvc', 'DetStore', 'EvtStore', 'ExtraInputs', 'ExtraOutputs', 'FastSimMasterTool', 'FlagAbortedEvents', 'G4Commands', 'GeoIDSvc', 'InputConverter', 'InputTruthCollection', 'KillAbortedEvents', 'MultiThreading', 'NeededResources', 'OutputTruthCollection', 'PhysicsListSvc', 'RandomGenerator', 'RecordFlux', 'ReleaseGeoModel', 'SenDetMasterTool', 'TruthRecordService', 'UserActionSvc', 'UserLimitsSvc', 'Verbosities']
         g4atlasalg = self._job_config_dict['G4AtlasAlg']
         actual_list = g4atlasalg.keys()
         expected_property_value_sorted = sorted(expected_list)
@@ -198,9 +198,9 @@ class TestAtlasG4(unittest.TestCase):
         self._assert_Algorithm_property_equal('G4AtlasAlg', 'SenDetMasterTool', expected_tool_name)
 
 
-    def test___G4AtlasAlg_G4AtlasSvc_setCorrectly(self):
-        expected_service_name = 'G4AtlasSvc'
-        self._assert_Algorithm_property_equal('G4AtlasAlg', 'G4AtlasSvc', expected_service_name)
+    def test___G4AtlasAlg_UserLimitsSvc_setCorrectly(self):
+        expected_service_name = 'UserLimitsSvc'
+        self._assert_Algorithm_property_equal('G4AtlasAlg', 'UserLimitsSvc', expected_service_name)
 
 
     def test___G4AtlasAlg_UserActionSvc_setCorrectly(self):
diff --git a/Simulation/G4Atlas/G4AtlasApps/test/test_TestBeam_tf_configuration.py b/Simulation/G4Atlas/G4AtlasApps/test/test_TestBeam_tf_configuration.py
index 7f50108358f495dee92189a93fd0cb93d01ea643..2e759db16d67f38bb29cce937d4c6ef4b05dc472 100755
--- a/Simulation/G4Atlas/G4AtlasApps/test/test_TestBeam_tf_configuration.py
+++ b/Simulation/G4Atlas/G4AtlasApps/test/test_TestBeam_tf_configuration.py
@@ -174,7 +174,7 @@ class TestTestBeam(unittest.TestCase):
 
 
     def test___G4AtlasAlg_ListOfSetProperties(self):
-        expected_list = ['AtRndmGenSvc', 'DetGeoSvc', 'DetStore', 'EvtStore', 'ExtraInputs', 'ExtraOutputs', 'FastSimMasterTool', 'G4AtlasSvc', 'G4Commands', 'GeoIDSvc', 'InputConverter', 'InputTruthCollection', 'MultiThreading', 'NeededResources', 'OutputTruthCollection', 'PhysicsListSvc', 'RandomGenerator', 'RecordFlux', 'ReleaseGeoModel', 'SenDetMasterTool', 'TruthRecordService', 'UserActionSvc', 'Verbosities']
+        expected_list = ['AtRndmGenSvc', 'DetGeoSvc', 'DetStore', 'EvtStore', 'ExtraInputs', 'ExtraOutputs', 'FastSimMasterTool', 'G4Commands', 'GeoIDSvc', 'InputConverter', 'InputTruthCollection', 'MultiThreading', 'NeededResources', 'OutputTruthCollection', 'PhysicsListSvc', 'RandomGenerator', 'RecordFlux', 'ReleaseGeoModel', 'SenDetMasterTool', 'TruthRecordService', 'UserActionSvc', 'UserLimitsSvc', 'Verbosities']
         g4atlasalg = self._job_config_dict['G4AtlasAlg']
         actual_list = g4atlasalg.keys()
         expected_property_value_sorted = sorted(expected_list)
@@ -218,9 +218,9 @@ class TestTestBeam(unittest.TestCase):
         self._assert_Algorithm_property_equal('G4AtlasAlg', 'SenDetMasterTool', expected_tool_name)
 
 
-    def test___G4AtlasAlg_G4AtlasSvc_setCorrectly(self):
-        expected_service_name = 'G4AtlasSvc'
-        self._assert_Algorithm_property_equal('G4AtlasAlg', 'G4AtlasSvc', expected_service_name)
+    def test___G4AtlasAlg_UserLimitsSvc_setCorrectly(self):
+        expected_service_name = 'UserLimitsSvc'
+        self._assert_Algorithm_property_equal('G4AtlasAlg', 'UserLimitsSvc', expected_service_name)
 
 
     def test___G4AtlasAlg_UserActionSvc_setCorrectly(self):
diff --git a/Simulation/ISF/ISF_Geant4/ISF_Geant4Tools/src/TransportTool.cxx b/Simulation/ISF/ISF_Geant4/ISF_Geant4Tools/src/TransportTool.cxx
index 45b164a2b6d64d94ba5b4c582361b01e502760ee..66b0e0bdb16cb63844caa55906c220dd209a716c 100644
--- a/Simulation/ISF/ISF_Geant4/ISF_Geant4Tools/src/TransportTool.cxx
+++ b/Simulation/ISF/ISF_Geant4/ISF_Geant4Tools/src/TransportTool.cxx
@@ -38,6 +38,9 @@
 #include "G4ScoringManager.hh"
 #include "G4Timer.hh"
 #include "G4SDManager.hh"
+#include "G4VUserPhysicsList.hh"
+#include "G4VModularPhysicsList.hh"
+#include "G4ParallelWorldPhysics.hh"
 
 #include "AtlasDetDescr/AtlasRegionHelper.h"
 
@@ -63,6 +66,7 @@ iGeant4::G4TransportTool::G4TransportTool(const std::string& type,
   //declareProperty("KillAllNeutrinos",      m_KillAllNeutrinos=true);
   //declareProperty("KillLowEPhotons",       m_KillLowEPhotons=-1.);
   declareProperty("PrintTimingInfo",      m_doTiming       );
+  declareProperty("ActivateParallelWorlds",m_activateParallelGeometries,"Toggle on/off the G4 parallel geometry system");
 
 }
 
@@ -84,9 +88,10 @@ StatusCode iGeant4::G4TransportTool::initialize()
     m_runTimer->Start();
   }
 
-  ATH_CHECK(m_inputConverter.retrieve());
+  // Create the scoring manager if requested
+  if (m_recordFlux) G4ScoringManager::GetScoringManager();
 
- // One-time initialization
+  // One-time initialization
   try {
     std::call_once(initializeOnceFlag, &iGeant4::G4TransportTool::initializeOnce, this);
   }
@@ -98,15 +103,11 @@ StatusCode iGeant4::G4TransportTool::initialize()
   ATH_CHECK( m_rndmGenSvc.retrieve() );
   ATH_CHECK( m_userActionSvc.retrieve() );
 
-  ATH_CHECK(m_g4atlasSvc.retrieve());
-
-  if (m_recordFlux) G4ScoringManager::GetScoringManager();
-
-  ATH_CHECK (m_detGeoSvc.retrieve());
-
   ATH_CHECK(m_senDetTool.retrieve());
   ATH_CHECK(m_fastSimTool.retrieve());
 
+  ATH_CHECK(m_inputConverter.retrieve());
+
   return StatusCode::SUCCESS;
 }
 
@@ -167,6 +168,39 @@ void iGeant4::G4TransportTool::initializeOnce()
     commandLog(returnCode, g4command);
   }
 
+  // Code from G4AtlasSvc
+  auto* rm = G4RunManager::GetRunManager();
+  if(!rm) {
+    throw std::runtime_error("Run manager retrieval has failed");
+  }
+  rm->Initialize();     // Initialization differs slightly in multi-threading.
+  // TODO: add more details about why this is here.
+  if(!m_useMT && rm->ConfirmBeamOnCondition()) {
+    rm->RunInitialization();
+  }
+
+  ATH_MSG_INFO( "retireving the Detector Geometry Service" );
+  if(m_detGeoSvc.retrieve().isFailure()) {
+    throw std::runtime_error("Could not initialize ATLAS DetectorGeometrySvc!");
+  }
+
+  if(m_userLimitsSvc.retrieve().isFailure()) {
+    throw std::runtime_error("Could not initialize ATLAS UserLimitsSvc!");
+  }
+
+  if (m_activateParallelGeometries) {
+    G4VModularPhysicsList* thePhysicsList=dynamic_cast<G4VModularPhysicsList*>(m_physListSvc->GetPhysicsList());
+    if (!thePhysicsList) {
+      throw std::runtime_error("Failed dynamic_cast!! this is not a G4VModularPhysicsList!");
+    }
+#if G4VERSION_NUMBER >= 1010
+    std::vector<std::string>& parallelWorldNames=m_detGeoSvc->GetParallelWorldNames();
+    for (auto& it: parallelWorldNames) {
+      thePhysicsList->RegisterPhysics(new G4ParallelWorldPhysics(it,true));
+    }
+#endif
+  }
+
   return;
 }
 
diff --git a/Simulation/ISF/ISF_Geant4/ISF_Geant4Tools/src/TransportTool.h b/Simulation/ISF/ISF_Geant4/ISF_Geant4Tools/src/TransportTool.h
index df0871742452a50565d491bae26a083051299d59..ee3fcd7668243331600cd73dbdb42016984ab80a 100644
--- a/Simulation/ISF/ISF_Geant4/ISF_Geant4Tools/src/TransportTool.h
+++ b/Simulation/ISF/ISF_Geant4/ISF_Geant4Tools/src/TransportTool.h
@@ -19,12 +19,12 @@
 // Athena headers
 #include "AthenaKernel/IAthRNGSvc.h"
 #include "AthenaKernel/SlotSpecificObj.h"
-#include "G4AtlasInterfaces/IG4AtlasSvc.h"
 #include "G4AtlasInterfaces/IUserActionSvc.h"
 #include "G4AtlasInterfaces/IDetectorGeometrySvc.h"
 #include "G4AtlasInterfaces/ISensitiveDetectorMasterTool.h"
 #include "G4AtlasInterfaces/IFastSimulationMasterTool.h"
 #include "G4AtlasInterfaces/IPhysicsListSvc.h"
+#include "G4AtlasInterfaces/IUserLimitsSvc.h"
 #include "CxxUtils/checker_macros.h"
 
 // ISF includes
@@ -134,10 +134,11 @@ namespace iGeant4
     std::vector<std::string> m_g4commands;
     /// Activate multi-threading configuration
     bool m_useMT{false};
+    bool m_activateParallelGeometries{false};
     // Random number service
     ServiceHandle<IAthRNGSvc> m_rndmGenSvc{this, "RandomNumberService", "AthRNGSvc", ""};
-    /// G4AtlasSvc
-    ServiceHandle<IG4AtlasSvc> m_g4atlasSvc{this, "G4AtlasSvc", "G4AtlasSvc", ""};
+    ///
+    ServiceHandle<IUserLimitsSvc> m_userLimitsSvc{this, "UserLimitsSvc", "UserLimitsSvc", ""};
     /// user action service
     ServiceHandle<G4UA::IUserActionSvc> m_userActionSvc{this, "UserActionSvc", "", ""};
     /// Detector Geometry Service (builds G4 Geometry)