diff --git a/Simulation/G4Extensions/Monopole/src/G4mplEquationSetup.cxx b/Simulation/G4Extensions/Monopole/src/G4mplEquationSetup.cxx
index 796ebf49b9d8ca46da1ac4da32cf3d9edf0ea9e8..724428a9c654fd0cbc910235cb46c41e2cf3c952 100644
--- a/Simulation/G4Extensions/Monopole/src/G4mplEquationSetup.cxx
+++ b/Simulation/G4Extensions/Monopole/src/G4mplEquationSetup.cxx
@@ -35,6 +35,12 @@
 //
 
 // =======================================================================
+// Modified: 19 May 2019, M. Bandieramonte: introduced MT mode. The class
+//           was a Singleton and it was not thread-safe.
+//           Added the #ifdef G4MULTITHREADED directive to handle
+//           the multithreaded case. One instance of the class will be created
+//           per each thread and stored in a tbb::concurrent_unordered_map that
+//           is hashed with the threadID number.
 // Modified: 28 August 2013, W. Taylor: adapted for ATLAS
 // Created:  23 May 2013,    J. Apostolakis
 //            Adapted from G4MonopoleFieldSetup by B. Bozsogi
@@ -62,7 +68,11 @@
 //    AtlasRK4, NystromRK4 - these have the equation of motion embedded inside
 #include "G4SystemOfUnits.hh"
 
-G4mplEquationSetup* G4mplEquationSetup::fG4mplEquationSetup=0;
+#ifdef G4MULTITHREADED
+G4mplEquationSetup::ESThreadMap_t G4mplEquationSetup::m_ESThreadMap;
+#endif
+
+
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
@@ -84,13 +94,17 @@ G4mplEquationSetup::G4mplEquationSetup():
 
 G4mplEquationSetup* G4mplEquationSetup::GetInstance()
 {
-  if ( fG4mplEquationSetup == 0 )
-    {
-      static G4mplEquationSetup theInstance;
-      fG4mplEquationSetup = &theInstance;
-    }
+#ifdef G4MULTITHREADED
+    auto es = getES();
+    if (!es) //nullpointer if it is not found
+      return setES();
+    else return es;
+#else
+    //Standard implementation of a Singleton Pattern
+    static G4mplEquationSetup instance;
+    return &instance;
+#endif
 
-  return fG4mplEquationSetup;
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -109,6 +123,30 @@ G4mplEquationSetup::~G4mplEquationSetup()
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
+#ifdef G4MULTITHREADED
+  G4mplEquationSetup* G4mplEquationSetup::getES()
+  {
+   // Get current thread-ID
+   const auto tid = std::this_thread::get_id();
+   auto esPair = m_ESThreadMap.find(tid);
+   if(esPair == m_ESThreadMap.end())
+      return nullptr; //if not found return null pointer
+   else return esPair->second;
+  }
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......   
+
+   G4mplEquationSetup* G4mplEquationSetup::setES()
+   {
+      G4mplEquationSetup* instance = new G4mplEquationSetup;
+      const auto tid = std::this_thread::get_id();
+      auto inserted = m_ESThreadMap.insert( std::make_pair(tid, instance)).first;
+      return (G4mplEquationSetup*) inserted->second;
+   }
+#endif
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
 void
 G4mplEquationSetup::InitialiseForField(G4FieldManager* fieldManager )
 {
diff --git a/Simulation/G4Extensions/Monopole/src/G4mplEquationSetup.hh b/Simulation/G4Extensions/Monopole/src/G4mplEquationSetup.hh
index 87961cbd7afbbc1f6002b196683c553516348427..fca58626b2b3febeaf0df8fe4bc905b42ce36349 100644
--- a/Simulation/G4Extensions/Monopole/src/G4mplEquationSetup.hh
+++ b/Simulation/G4Extensions/Monopole/src/G4mplEquationSetup.hh
@@ -30,6 +30,13 @@
 // particles.
 //
 // =======================================================================
+// Modified: 19 May 2019, M. Bandieramonte: introduced MT mode. The class
+// 	     was a Singleton and it was not thread-safe.
+//           Added the #ifdef G4MULTITHREADED directive to handle
+//           the multithreaded case. One instance of the class will be created
+//           per each thread and stored in a tbb::concurrent_unordered_map that
+//           is hashed with the threadID number.
+//
 // Modified: 28 August 2013, W. Taylor: adapted for ATLAS
 // Created:  23 May 2013, J. Apostolakis
 //            Adapted from G4MonopoleFieldSetup by B. Bozsogi
@@ -38,6 +45,11 @@
 #ifndef MONOPOLE_G4mplEquationSetup_H
 #define MONOPOLE_G4mplEquationSetup_H
 
+#include <thread>
+#ifdef G4MULTITHREADED
+#  include "tbb/concurrent_unordered_map.h"
+#endif
+
 // Geant4 headers
 #include "G4MagneticField.hh"
 
@@ -63,7 +75,21 @@ public:
 
   ~G4mplEquationSetup() ;
 
-private:
+private: 
+
+#ifdef G4MULTITHREADED
+     // Thread-to-EquationSetup concurrent map type
+     using ESThreadMap_t = tbb::concurrent_unordered_map< std::thread::id, G4mplEquationSetup*, std::hash<std::thread::id> >;
+     // Concurrent map of EquationsSetup, one for each thread
+     static ESThreadMap_t m_ESThreadMap;
+     //@brief Search inside m_ESThreadMap the element with the current threadID 
+     // and return it or return a null pointer if the element is not found
+     static G4mplEquationSetup* getES();
+     // @brief Insert the current ES in m_ESThreadMap and 
+     // associate it with the current threadID
+     static G4mplEquationSetup* setES();
+ #endif
+     //
 
   G4mplEquationSetup();
 
@@ -83,8 +109,6 @@ private:
   G4MagIntegratorStepper* fStepper ;
   G4bool                  fCreatedOrdinaryStepper; // If set, created stepper.
 
-  // For Singleton
-  static G4mplEquationSetup*  fG4mplEquationSetup;
   G4bool                      fVerbose;
   //
   // State - changed during tracking