diff --git a/Framework/API/test/MatrixWorkspaceTest.h b/Framework/API/test/MatrixWorkspaceTest.h
index fb15400a8baa7682b00e3eb56acf99b5a86ff8ce..492ef16ff0ac90bc79ca8ebfe045de51c36cf16b 100644
--- a/Framework/API/test/MatrixWorkspaceTest.h
+++ b/Framework/API/test/MatrixWorkspaceTest.h
@@ -326,7 +326,7 @@ public:
     auto workspace = makeWorkspaceWithDetectors(numHist, 1);
     std::atomic<bool> parallelException{false};
     std::atomic<int> threadCount{1};
-    PARALLEL_FOR1(workspace)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*workspace))
     for (int i = 0; i < numHist; ++i) {
       // Note: Cannot use INTERUPT_REGION macros since not inside an Algorithm.
       threadCount = PARALLEL_NUMBER_OF_THREADS;
diff --git a/Framework/API/test/SpectrumInfoTest.h b/Framework/API/test/SpectrumInfoTest.h
index 606620680f84a0d85ac130ad59e63a5335727af0..6cebe7ba734a73b663a367871320c8b7c9bb5a9b 100644
--- a/Framework/API/test/SpectrumInfoTest.h
+++ b/Framework/API/test/SpectrumInfoTest.h
@@ -97,7 +97,7 @@ public:
     SpectrumInfo info(*ws);
     // This attempts to test threading, but probably it is not really exercising
     // much.
-    PARALLEL_FOR1(ws)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*ws))
     for (int i = 0; i < count; ++i)
       TS_ASSERT_EQUALS(info.isMasked(static_cast<size_t>(i)), i % 2 == 0);
   }
diff --git a/Framework/Algorithms/src/AbsorptionCorrection.cpp b/Framework/Algorithms/src/AbsorptionCorrection.cpp
index d57c910a7a18625c75faa1de6f51c40eed044694..d8291eebacbdcca09310136021a112ae3da312b0 100644
--- a/Framework/Algorithms/src/AbsorptionCorrection.cpp
+++ b/Framework/Algorithms/src/AbsorptionCorrection.cpp
@@ -140,7 +140,7 @@ void AbsorptionCorrection::exec() {
 
   Progress prog(this, 0.0, 1.0, numHists);
   // Loop over the spectra
-  PARALLEL_FOR2(m_inputWS, correctionFactors)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*m_inputWS, *correctionFactors))
   for (int64_t i = 0; i < int64_t(numHists); ++i) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/ApplyTransmissionCorrection.cpp b/Framework/Algorithms/src/ApplyTransmissionCorrection.cpp
index 1dbaf00f3df6a498ad9765bfd1f32ad37962bbbf..469c8d3e8facc26b39593aca19863eb8c08b8a9a 100644
--- a/Framework/Algorithms/src/ApplyTransmissionCorrection.cpp
+++ b/Framework/Algorithms/src/ApplyTransmissionCorrection.cpp
@@ -89,7 +89,7 @@ void ApplyTransmissionCorrection::exec() {
   const auto &spectrumInfo = inputWS->spectrumInfo();
 
   // Loop through the spectra and apply correction
-  PARALLEL_FOR2(inputWS, corrWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *corrWS))
   for (int i = 0; i < numHists; i++) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/CalculateFlatBackground.cpp b/Framework/Algorithms/src/CalculateFlatBackground.cpp
index 8e7d06c8edb29f49d50ebbeac2820663fbdfd371..5de0b7ac98b02f3a54af220db72bd67655709175 100644
--- a/Framework/Algorithms/src/CalculateFlatBackground.cpp
+++ b/Framework/Algorithms/src/CalculateFlatBackground.cpp
@@ -161,7 +161,7 @@ void CalculateFlatBackground::exec() {
   // the output
   if (outputWS != inputWS) {
     outputWS = WorkspaceFactory::Instance().create(inputWS);
-    PARALLEL_FOR2(inputWS, outputWS)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
     for (int i = 0; i < numHists; ++i) {
       PARALLEL_START_INTERUPT_REGION
       outputWS->setHistogram(i, inputWS->histogram(i));
diff --git a/Framework/Algorithms/src/ChopData.cpp b/Framework/Algorithms/src/ChopData.cpp
index a0e9556d92e9f2d8482f3028ec2242412b97502e..3867c52997728d2e30379965465cbccc5e942216 100644
--- a/Framework/Algorithms/src/ChopData.cpp
+++ b/Framework/Algorithms/src/ChopData.cpp
@@ -129,7 +129,7 @@ void ChopData::exec() {
                                                          nbins + 1, nbins);
 
     // Copy over X, Y and E data
-    PARALLEL_FOR2(inputWS, workspace)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *workspace))
     for (int j = 0; j < nHist; j++) {
 
       auto edges = inputWS->binEdges(j);
diff --git a/Framework/Algorithms/src/ConvertAxesToRealSpace.cpp b/Framework/Algorithms/src/ConvertAxesToRealSpace.cpp
index a04e46808a82c0a1261985acedd7d95352e332ef..34064c126408e2105adbe61508039e32c61a972f 100644
--- a/Framework/Algorithms/src/ConvertAxesToRealSpace.cpp
+++ b/Framework/Algorithms/src/ConvertAxesToRealSpace.cpp
@@ -114,7 +114,7 @@ void ConvertAxesToRealSpace::exec() {
 
   const auto &spectrumInfo = summedWs->spectrumInfo();
   // for each spectra
-  PARALLEL_FOR2(summedWs, outputWs)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*summedWs, *outputWs))
   for (int i = 0; i < nHist; ++i) {
     try {
       V3D pos = spectrumInfo.position(i);
@@ -228,7 +228,7 @@ void ConvertAxesToRealSpace::exec() {
 
   // set all the X arrays - share the same vector
   int nOutputHist = static_cast<int>(outputWs->getNumberHistograms());
-  PARALLEL_FOR1(outputWs)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*outputWs))
   for (int i = 0; i < nOutputHist; ++i) {
     outputWs->setPoints(i, x);
   }
@@ -255,7 +255,7 @@ void ConvertAxesToRealSpace::exec() {
   }
 
   // loop over the data and sqrt the errors to complete the error calculation
-  PARALLEL_FOR1(outputWs)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*outputWs))
   for (int i = 0; i < nOutputHist; ++i) {
     auto &errorVec = outputWs->mutableE(i);
     std::transform(errorVec.begin(), errorVec.end(), errorVec.begin(),
diff --git a/Framework/Algorithms/src/ConvertAxisByFormula.cpp b/Framework/Algorithms/src/ConvertAxisByFormula.cpp
index 050fb69fc574bca7cdf9dffe1b80f842e0cbe004..5941af7807e84f817f720ee98027f19a5d52d29c 100644
--- a/Framework/Algorithms/src/ConvertAxisByFormula.cpp
+++ b/Framework/Algorithms/src/ConvertAxisByFormula.cpp
@@ -214,7 +214,7 @@ void ConvertAxisByFormula::exec() {
           outputWs->getNumberHistograms()); // cast to make openmp happy
       auto xVals = outputWs->refX(0);
       Progress prog(this, 0.6, 1.0, numberOfSpectra_i);
-      PARALLEL_FOR1(outputWs)
+      PARALLEL_FOR_IF(Kernel::threadSafe(*outputWs))
       for (int64_t j = 1; j < numberOfSpectra_i; ++j) {
         PARALLEL_START_INTERUPT_REGION
         outputWs->setX(j, xVals);
diff --git a/Framework/Algorithms/src/ConvertToEventWorkspace.cpp b/Framework/Algorithms/src/ConvertToEventWorkspace.cpp
index c98105f6c15e0df7d830c8ce3d4f2ea676aca32b..7aa9ba3508a64f4b71f453850f3697364995366d 100644
--- a/Framework/Algorithms/src/ConvertToEventWorkspace.cpp
+++ b/Framework/Algorithms/src/ConvertToEventWorkspace.cpp
@@ -62,7 +62,7 @@ void ConvertToEventWorkspace::exec() {
   API::WorkspaceFactory::Instance().initializeFromParent(inWS, outWS, false);
 
   Progress prog(this, 0.0, 1.0, inWS->getNumberHistograms());
-  PARALLEL_FOR1(inWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inWS))
   for (int iwi = 0; iwi < int(inWS->getNumberHistograms()); iwi++) {
     PARALLEL_START_INTERUPT_REGION
     size_t wi = size_t(iwi);
diff --git a/Framework/Algorithms/src/ConvertToMatrixWorkspace.cpp b/Framework/Algorithms/src/ConvertToMatrixWorkspace.cpp
index 0b861ef3f2b5bc09be429924d802f1080b731460..41f179a671d1b75ae08e52209b147841d9c50ce1 100644
--- a/Framework/Algorithms/src/ConvertToMatrixWorkspace.cpp
+++ b/Framework/Algorithms/src/ConvertToMatrixWorkspace.cpp
@@ -48,7 +48,7 @@ void ConvertToMatrixWorkspace::exec() {
     outputWorkspace = WorkspaceFactory::Instance().create(inputWorkspace);
 
     // ...but not the data, so do that here.
-    PARALLEL_FOR2(inputWorkspace, outputWorkspace)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*inputWorkspace, *outputWorkspace))
     for (int64_t i = 0; i < static_cast<int64_t>(numHists); ++i) {
       PARALLEL_START_INTERUPT_REGION
       const auto &inSpec = inputWorkspace->getSpectrum(i);
diff --git a/Framework/Algorithms/src/ConvertUnits.cpp b/Framework/Algorithms/src/ConvertUnits.cpp
index 31623fce6eb6b77e609efda8587d0420c8029d03..47cac08c4a46ba48e122dfdd74e4389d87d4af8d 100644
--- a/Framework/Algorithms/src/ConvertUnits.cpp
+++ b/Framework/Algorithms/src/ConvertUnits.cpp
@@ -282,7 +282,7 @@ API::MatrixWorkspace_sptr ConvertUnits::setupOutputWorkspace(
   if (!m_inputEvents && m_distribution) {
     // Loop over the histograms (detector spectra)
     Progress prog(this, 0.0, 0.2, m_numberOfSpectra);
-    PARALLEL_FOR1(outputWS)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS))
     for (int64_t i = 0; i < static_cast<int64_t>(m_numberOfSpectra); ++i) {
       PARALLEL_START_INTERUPT_REGION
       // Take the bin width dependency out of the Y & E data
@@ -348,7 +348,7 @@ ConvertUnits::convertQuickly(API::MatrixWorkspace_const_sptr inputWS,
 
       auto xVals = outputWS->sharedX(0);
 
-      PARALLEL_FOR1(outputWS)
+      PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS))
       for (int64_t j = 1; j < numberOfSpectra_i; ++j) {
         PARALLEL_START_INTERUPT_REGION
         outputWS->setX(j, xVals);
@@ -368,7 +368,7 @@ ConvertUnits::convertQuickly(API::MatrixWorkspace_const_sptr inputWS,
   // If we get to here then the bins weren't aligned and each spectrum is
   // unique
   // Loop over the histograms (detector spectra)
-  PARALLEL_FOR1(outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS))
   for (int64_t k = 0; k < numberOfSpectra_i; ++k) {
     PARALLEL_START_INTERUPT_REGION
     if (!commonBoundaries) {
@@ -670,7 +670,7 @@ void ConvertUnits::reverse(API::MatrixWorkspace_sptr WS) {
   } else {
     // either events or ragged boundaries
     int numberOfSpectra_i = static_cast<int>(numberOfSpectra);
-    PARALLEL_FOR1(WS)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*WS))
     for (int j = 0; j < numberOfSpectra_i; ++j) {
       PARALLEL_START_INTERUPT_REGION
       if (isInputEvents) {
diff --git a/Framework/Algorithms/src/ConvertUnitsUsingDetectorTable.cpp b/Framework/Algorithms/src/ConvertUnitsUsingDetectorTable.cpp
index 8b57d341e7e833955bb9cb70312da437f6b05ab3..9906fb3a512a4d65e02ca9485a32508cb45d996c 100644
--- a/Framework/Algorithms/src/ConvertUnitsUsingDetectorTable.cpp
+++ b/Framework/Algorithms/src/ConvertUnitsUsingDetectorTable.cpp
@@ -176,7 +176,7 @@ MatrixWorkspace_sptr ConvertUnitsUsingDetectorTable::convertViaTOF(
 
   // TODO: Check why this parallel stuff breaks
   // Loop over the histograms (detector spectra)
-  // PARALLEL_FOR1(outputWS)
+  // PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS))
   for (int64_t i = 0; i < numberOfSpectra_i; ++i) {
 
     // Lets find what row this spectrum Number appears in our detector table.
diff --git a/Framework/Algorithms/src/CorelliCrossCorrelate.cpp b/Framework/Algorithms/src/CorelliCrossCorrelate.cpp
index f66241d770cb7cc0effb9ead55caf9dd78c11050..ff4a27633e09540e514db233157aa51ed61540c5 100644
--- a/Framework/Algorithms/src/CorelliCrossCorrelate.cpp
+++ b/Framework/Algorithms/src/CorelliCrossCorrelate.cpp
@@ -185,7 +185,7 @@ void CorelliCrossCorrelate::exec() {
   int64_t numHistograms = static_cast<int64_t>(inputWS->getNumberHistograms());
   API::Progress prog = API::Progress(this, 0.0, 1.0, numHistograms);
   const auto &spectrumInfo = inputWS->spectrumInfo();
-  PARALLEL_FOR1(outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS))
   for (int64_t i = 0; i < numHistograms; ++i) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/CorrectFlightPaths.cpp b/Framework/Algorithms/src/CorrectFlightPaths.cpp
index 61c54ca76ee4bc8266c791454c19f0559b1b5e24..4d286adfdfa93de553fc6504fbd03cb7bf2ca528 100644
--- a/Framework/Algorithms/src/CorrectFlightPaths.cpp
+++ b/Framework/Algorithms/src/CorrectFlightPaths.cpp
@@ -89,7 +89,7 @@ void CorrectFlightPaths::exec() {
   const auto &spectrumInfo = m_inputWS->spectrumInfo();
 
   // Loop over the histograms (detector spectra)
-  PARALLEL_FOR2(m_inputWS, m_outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*m_inputWS, *m_outputWS))
   for (int64_t i = 0; i < numberOfSpectra_i; ++i) {
     PARALLEL_START_INTERUPT_REGION
     m_outputWS->setHistogram(i, m_outputWS->histogram(i));
diff --git a/Framework/Algorithms/src/CorrectKiKf.cpp b/Framework/Algorithms/src/CorrectKiKf.cpp
index 96b920a13483768b884f52ee4358c34b4aabbb01..49611e20d851c8d88125234d0f79abaad8a2db74 100644
--- a/Framework/Algorithms/src/CorrectKiKf.cpp
+++ b/Framework/Algorithms/src/CorrectKiKf.cpp
@@ -96,7 +96,7 @@ void CorrectKiKf::exec() {
   // Get the parameter map
   const ParameterMap &pmap = outputWS->constInstrumentParameters();
 
-  PARALLEL_FOR2(inputWS, outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
   for (int64_t i = 0; i < int64_t(numberOfSpectra); ++i) {
     PARALLEL_START_INTERUPT_REGION
     double Efi = 0;
@@ -222,7 +222,7 @@ void CorrectKiKf::execEvent() {
 
   int64_t numHistograms = static_cast<int64_t>(inputWS->getNumberHistograms());
   API::Progress prog = API::Progress(this, 0.0, 1.0, numHistograms);
-  PARALLEL_FOR1(outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS))
   for (int64_t i = 0; i < numHistograms; ++i) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/CreatePSDBleedMask.cpp b/Framework/Algorithms/src/CreatePSDBleedMask.cpp
index 6de7570375a7d6ebf22ce769975839b51a4add86..3c096be5f472d65eaacaf2f4717970b1749f2556 100644
--- a/Framework/Algorithms/src/CreatePSDBleedMask.cpp
+++ b/Framework/Algorithms/src/CreatePSDBleedMask.cpp
@@ -149,11 +149,10 @@ void CreatePSDBleedMask::exec() {
 
   progress.resetNumSteps(numTubes, 0, 1);
 
-  PARALLEL_FOR2(inputWorkspace, outputWorkspace)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWorkspace, *outputWorkspace))
   for (int i = 0; i < numTubes; ++i) {
     PARALLEL_START_INTERUPT_REGION
-    auto current = tubeMap.begin();
-    std::advance(current, i);
+    auto current = std::next(tubeMap.begin(), i);
     const TubeIndex::mapped_type tubeIndices = current->second;
     bool mask = performBleedTest(tubeIndices, inputWorkspace, maxRate,
                                  numIgnoredPixels);
diff --git a/Framework/Algorithms/src/CreateWorkspace.cpp b/Framework/Algorithms/src/CreateWorkspace.cpp
index 7758b7d839840cce8d2e4d6b84f0c4868f2ee6ae..37d473426418ea0314c4d5535d3a39bd856c8400 100644
--- a/Framework/Algorithms/src/CreateWorkspace.cpp
+++ b/Framework/Algorithms/src/CreateWorkspace.cpp
@@ -158,7 +158,7 @@ void CreateWorkspace::exec() {
 
   Progress progress(this, 0, 1, nSpec);
 
-  PARALLEL_FOR1(outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS))
   for (int i = 0; i < nSpec; i++) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/CrossCorrelate.cpp b/Framework/Algorithms/src/CrossCorrelate.cpp
index 1edbca628905ae2edab0c522b9787f056ecb1371..27a527d38c44776744bb2bd913dc589e443f7115 100644
--- a/Framework/Algorithms/src/CrossCorrelate.cpp
+++ b/Framework/Algorithms/src/CrossCorrelate.cpp
@@ -169,7 +169,7 @@ void CrossCorrelate::exec() {
   // Initialise the progress reporting object
   out->mutableX(0) = XX;
   m_progress = new Progress(this, 0.0, 1.0, nspecs);
-  PARALLEL_FOR2(inputWS, out)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *out))
   for (int i = 0; i < nspecs; ++i) // Now loop on all spectra
   {
     PARALLEL_START_INTERUPT_REGION
diff --git a/Framework/Algorithms/src/DetectorEfficiencyCor.cpp b/Framework/Algorithms/src/DetectorEfficiencyCor.cpp
index e396a47405e4d66de2705c082ac47205f02445ae..f76796d8070a2a7f5c51ade83d80b9d9f863fcd3 100644
--- a/Framework/Algorithms/src/DetectorEfficiencyCor.cpp
+++ b/Framework/Algorithms/src/DetectorEfficiencyCor.cpp
@@ -121,7 +121,7 @@ void DetectorEfficiencyCor::exec() {
   double numHists_d = static_cast<double>(numHists);
   const int64_t progStep = static_cast<int64_t>(ceil(numHists_d / 100.0));
 
-  PARALLEL_FOR2(m_inputWS, m_outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*m_inputWS, *m_outputWS))
   for (int64_t i = 0; i < numHists; ++i) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/DetectorEfficiencyCorUser.cpp b/Framework/Algorithms/src/DetectorEfficiencyCorUser.cpp
index 233b55d116e72f7d79eb214c67f0657234b6a168..66ec494094395d3a0564e14bdcf33b991cb0bb69 100644
--- a/Framework/Algorithms/src/DetectorEfficiencyCorUser.cpp
+++ b/Framework/Algorithms/src/DetectorEfficiencyCorUser.cpp
@@ -76,7 +76,7 @@ void DetectorEfficiencyCorUser::exec() {
       static_cast<int64_t>(numberOfSpectra); // cast to make openmp happy
 
   // Loop over the histograms (detector spectra)
-  PARALLEL_FOR2(m_outputWS, m_inputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*m_outputWS, *m_inputWS))
   for (int64_t i = 0; i < numberOfSpectra_i; ++i) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/DiffractionFocussing2.cpp b/Framework/Algorithms/src/DiffractionFocussing2.cpp
index b0f6575b90663f62a7d194c09f98be6df3449e85..38cef9356043c7e7e103fcbbd513ca01d881c63e 100644
--- a/Framework/Algorithms/src/DiffractionFocussing2.cpp
+++ b/Framework/Algorithms/src/DiffractionFocussing2.cpp
@@ -170,9 +170,8 @@ void DiffractionFocussing2::exec() {
   Progress *prog;
   prog = new API::Progress(this, 0.2, 1.00,
                            static_cast<int>(totalHistProcess) + nGroups);
-#ifndef __APPLE__
-  PARALLEL_FOR2(m_matrixInputW, out)
-#endif
+
+  PARALLEL_FOR_IF(Kernel::threadSafe(*m_matrixInputW, *out))
   for (int outWorkspaceIndex = 0;
        outWorkspaceIndex < static_cast<int>(m_validGroups.size());
        outWorkspaceIndex++) {
@@ -181,9 +180,8 @@ void DiffractionFocussing2::exec() {
 
     // Get the group
     auto it = group2xvector.find(group);
-    group2vectormap::difference_type dif =
-        std::distance(group2xvector.begin(), it);
-    auto &Xout = (*it).second;
+    auto dif = std::distance(group2xvector.begin(), it);
+    auto &Xout = it->second;
 
     // Assign the new X axis only once (i.e when this group is encountered the
     // first time)
@@ -433,7 +431,7 @@ void DiffractionFocussing2::execEvent() {
     // ------ PARALLELIZE BY GROUPS -------------------------
 
     int nValidGroups = static_cast<int>(this->m_validGroups.size());
-    PARALLEL_FOR1(m_eventW)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*m_eventW))
     for (int iGroup = 0; iGroup < nValidGroups; iGroup++) {
       PARALLEL_START_INTERUPT_REGION
       const int group = this->m_validGroups[iGroup];
diff --git a/Framework/Algorithms/src/EQSANSTofStructure.cpp b/Framework/Algorithms/src/EQSANSTofStructure.cpp
index 08904c774abca9f66f2920acd09139d6ebb38931..9ead21fc227ed6270966e6fc3485160c189c7ce9 100644
--- a/Framework/Algorithms/src/EQSANSTofStructure.cpp
+++ b/Framework/Algorithms/src/EQSANSTofStructure.cpp
@@ -131,7 +131,7 @@ void EQSANSTofStructure::execEvent(
   const auto l1 = spectrumInfo.l1();
 
   // Loop through the spectra and apply correction
-  PARALLEL_FOR1(inputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS))
   for (int64_t ispec = 0; ispec < int64_t(numHists); ++ispec) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/ExtractFFTSpectrum.cpp b/Framework/Algorithms/src/ExtractFFTSpectrum.cpp
index c7c9c5c441845e751745154be265ed91a11dc470..41350451a879537b5bf84a9eae5b85e11d0b156e 100644
--- a/Framework/Algorithms/src/ExtractFFTSpectrum.cpp
+++ b/Framework/Algorithms/src/ExtractFFTSpectrum.cpp
@@ -45,7 +45,7 @@ void ExtractFFTSpectrum::exec() {
 
   Progress prog(this, 0.0, 1.0, numHists);
 
-  PARALLEL_FOR1(outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS))
   for (int i = 0; i < numHists; i++) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/ExtractMask.cpp b/Framework/Algorithms/src/ExtractMask.cpp
index 7bb5ea25d3db88ffd1ce456a4dcde2348dbb2916..4c735bd8db9d69a8b84378a4f35949a62ba7db90 100644
--- a/Framework/Algorithms/src/ExtractMask.cpp
+++ b/Framework/Algorithms/src/ExtractMask.cpp
@@ -61,7 +61,7 @@ void ExtractMask::exec() {
   const int64_t nHist = static_cast<int64_t>(inputWS->getNumberHistograms());
   Progress prog(this, 0.0, 1.0, nHist);
 
-  PARALLEL_FOR2(inputWS, maskWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *maskWS))
   for (int64_t i = 0; i < nHist; ++i) {
     PARALLEL_START_INTERUPT_REGION
     bool inputIsMasked(false);
diff --git a/Framework/Algorithms/src/ExtractSpectra.cpp b/Framework/Algorithms/src/ExtractSpectra.cpp
index 9581e735c9e85696290aaecc7ecaa8fa43e39800..c2871a39620fa2179458f8cc6313d17cd8aeabd7 100644
--- a/Framework/Algorithms/src/ExtractSpectra.cpp
+++ b/Framework/Algorithms/src/ExtractSpectra.cpp
@@ -293,7 +293,7 @@ void ExtractSpectra::execEvent() {
   Progress prog(this, 0.0, 1.0, 2 * m_workspaceIndexList.size());
   eventW->sortAll(Mantid::DataObjects::TOF_SORT, &prog);
   // Loop over the required workspace indices, copying in the desired bins
-  PARALLEL_FOR2(m_inputWorkspace, outputWorkspace)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*m_inputWorkspace, *outputWorkspace))
   for (int j = 0; j < static_cast<int>(m_workspaceIndexList.size()); ++j) {
     PARALLEL_START_INTERUPT_REGION
     auto i = m_workspaceIndexList[j];
diff --git a/Framework/Algorithms/src/FilterByXValue.cpp b/Framework/Algorithms/src/FilterByXValue.cpp
index 9ae5ee3133842468f3159582dde608701e8354e0..8a11b9add4f7157a1056f1ec2839ce6da0a44ce9 100644
--- a/Framework/Algorithms/src/FilterByXValue.cpp
+++ b/Framework/Algorithms/src/FilterByXValue.cpp
@@ -85,7 +85,7 @@ void FilterByXValue::exec() {
 
   Progress prog(this, 0.0, 1.0, numSpec);
   // Loop over the workspace, removing the events that don't pass the filter
-  PARALLEL_FOR1(outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS))
   for (int spec = 0; spec < numSpec; ++spec) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/GetDetectorOffsets.cpp b/Framework/Algorithms/src/GetDetectorOffsets.cpp
index 3f20683ea86634d330b737b1e1dd530d05581fe9..5ce99a549b31ded3c5c9594601522212aa645330 100644
--- a/Framework/Algorithms/src/GetDetectorOffsets.cpp
+++ b/Framework/Algorithms/src/GetDetectorOffsets.cpp
@@ -108,7 +108,7 @@ void GetDetectorOffsets::exec() {
 
   // Fit all the spectra with a gaussian
   Progress prog(this, 0, 1.0, nspec);
-  PARALLEL_FOR1(inputW)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputW))
   for (int wi = 0; wi < nspec; ++wi) {
     PARALLEL_START_INTERUPT_REGION
     // Fit the peak
diff --git a/Framework/Algorithms/src/He3TubeEfficiency.cpp b/Framework/Algorithms/src/He3TubeEfficiency.cpp
index aa3cbf5d2bd85a71de9da557b90c3b9f9f44c1a1..f660924b2d31cee0496c92b43f3123705d5eb1ae 100644
--- a/Framework/Algorithms/src/He3TubeEfficiency.cpp
+++ b/Framework/Algorithms/src/He3TubeEfficiency.cpp
@@ -117,7 +117,7 @@ void He3TubeEfficiency::exec() {
   std::size_t numHists = this->inputWS->getNumberHistograms();
   this->progress = new API::Progress(this, 0.0, 1.0, numHists);
 
-  PARALLEL_FOR2(inputWS, outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
   for (int i = 0; i < static_cast<int>(numHists); ++i) {
     PARALLEL_START_INTERUPT_REGION
 
@@ -422,7 +422,7 @@ void He3TubeEfficiency::execEvent() {
 
   std::size_t numHistograms = outputWS->getNumberHistograms();
   this->progress = new API::Progress(this, 0.0, 1.0, numHistograms);
-  PARALLEL_FOR1(outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS))
   for (int i = 0; i < static_cast<int>(numHistograms); ++i) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/IntegrateByComponent.cpp b/Framework/Algorithms/src/IntegrateByComponent.cpp
index 06b10990bcad9b228fba447e32c2dfed29bf0620..2b0999898b2bf2d7b042f5a45af684d958e4220c 100644
--- a/Framework/Algorithms/src/IntegrateByComponent.cpp
+++ b/Framework/Algorithms/src/IntegrateByComponent.cpp
@@ -78,7 +78,7 @@ void IntegrateByComponent::exec() {
       prog.report();
       std::vector<double> averageYInput, averageEInput;
 
-      PARALLEL_FOR1(integratedWS)
+      PARALLEL_FOR_IF(Kernel::threadSafe(*integratedWS))
       for (int i = 0; i < static_cast<int>(hists.size()); ++i) { // NOLINT
         PARALLEL_START_INTERUPT_REGION
 
@@ -115,7 +115,7 @@ void IntegrateByComponent::exec() {
             gsl_stats_mean(&averageEInput[0], 1, averageYInput.size()));
       }
 
-      PARALLEL_FOR1(integratedWS)
+      PARALLEL_FOR_IF(Kernel::threadSafe(*integratedWS))
       for (int i = 0; i < static_cast<int>(hists.size()); ++i) { // NOLINT
         PARALLEL_START_INTERUPT_REGION
         if (spectrumInfo.isMonitor(hists[i]))
diff --git a/Framework/Algorithms/src/Integration.cpp b/Framework/Algorithms/src/Integration.cpp
index b035bf502efc69fac6a2b829c8584aa3e158b5d1..0e06abbe111e8abc9de45f75ff46614390176744 100644
--- a/Framework/Algorithms/src/Integration.cpp
+++ b/Framework/Algorithms/src/Integration.cpp
@@ -141,7 +141,7 @@ void Integration::exec() {
   const bool axisIsNumeric = localworkspace->getAxis(1)->isNumeric();
 
   // Loop over spectra
-  PARALLEL_FOR2(localworkspace, outputWorkspace)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*localworkspace, *outputWorkspace))
   for (int i = minSpec; i <= maxSpec; ++i) {
     PARALLEL_START_INTERUPT_REGION
     // Workspace index on the output
diff --git a/Framework/Algorithms/src/LorentzCorrection.cpp b/Framework/Algorithms/src/LorentzCorrection.cpp
index 4e7b0d2fbcf6ff7a68bce5e2682a7469317373a1..765c520bc5b33ff9596aaf4299d7d19bcb202188 100644
--- a/Framework/Algorithms/src/LorentzCorrection.cpp
+++ b/Framework/Algorithms/src/LorentzCorrection.cpp
@@ -67,7 +67,7 @@ void LorentzCorrection::exec() {
   const auto &spectrumInfo = inWS->spectrumInfo();
   Progress prog(this, 0, 1, numHistos);
 
-  PARALLEL_FOR1(inWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inWS))
   for (int64_t i = 0; i < int64_t(numHistos); ++i) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/MaskBins.cpp b/Framework/Algorithms/src/MaskBins.cpp
index be0774bc56df61ce24849b13359b293312357b1f..454de81284b4f6471b86a1fc58f24190b39b3300 100644
--- a/Framework/Algorithms/src/MaskBins.cpp
+++ b/Framework/Algorithms/src/MaskBins.cpp
@@ -172,7 +172,7 @@ void MaskBins::execEvent() {
   // Go through all histograms
   if (!this->spectra_list.empty()) {
     // Specific spectra were specified
-    PARALLEL_FOR1(outputWS)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS))
     for (int i = 0; i < static_cast<int>(this->spectra_list.size()); // NOLINT
          ++i) {
       PARALLEL_START_INTERUPT_REGION
@@ -183,7 +183,7 @@ void MaskBins::execEvent() {
     PARALLEL_CHECK_INTERUPT_REGION
   } else {
     // Do all spectra!
-    PARALLEL_FOR1(outputWS)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS))
     for (int64_t i = 0; i < int64_t(numHists); ++i) {
       PARALLEL_START_INTERUPT_REGION
       outputWS->getSpectrum(i).maskTof(m_startX, m_endX);
diff --git a/Framework/Algorithms/src/MaxMin.cpp b/Framework/Algorithms/src/MaxMin.cpp
index 0e8bb126562dcd53cdda890a48a5036491ca9d18..4b8142b6b81858879ee9551b52ff219182355fe6 100644
--- a/Framework/Algorithms/src/MaxMin.cpp
+++ b/Framework/Algorithms/src/MaxMin.cpp
@@ -88,7 +88,7 @@ void MaxMin::exec() {
                                                MaxSpec - MinSpec + 1, 2, 1);
 
   Progress progress(this, 0, 1, (MaxSpec - MinSpec + 1));
-  PARALLEL_FOR2(localworkspace, outputWorkspace)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*localworkspace, *outputWorkspace))
   // Loop over spectra
   for (int i = MinSpec; i <= MaxSpec; ++i) {
     PARALLEL_START_INTERUPT_REGION
diff --git a/Framework/Algorithms/src/MedianDetectorTest.cpp b/Framework/Algorithms/src/MedianDetectorTest.cpp
index 98c72be2d0334abec7f1be8b3b136bf130793a8a..8844f32e172146c517f6e2be5ffcba795eda476d 100644
--- a/Framework/Algorithms/src/MedianDetectorTest.cpp
+++ b/Framework/Algorithms/src/MedianDetectorTest.cpp
@@ -251,7 +251,7 @@ int MedianDetectorTest::maskOutliers(
     std::vector<size_t> &hists = indexmap[i];
     double median = medianvec[i];
 
-    PARALLEL_FOR1(countsWS)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*countsWS))
     for (int j = 0; j < static_cast<int>(hists.size()); ++j) { // NOLINT
       const double value = countsWS->y(hists[j])[0];
       if ((value == 0.) && checkForMask) {
@@ -314,7 +314,7 @@ int MedianDetectorTest::doDetectorTests(
                     (instrument->getSample() != nullptr));
   }
 
-  PARALLEL_FOR2(countsWS, maskWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*countsWS, *maskWS))
   for (int j = 0; j < static_cast<int>(indexmap.size()); ++j) {
     std::vector<size_t> hists = indexmap.at(j);
     double median = medianvec.at(j);
diff --git a/Framework/Algorithms/src/ModeratorTzero.cpp b/Framework/Algorithms/src/ModeratorTzero.cpp
index 97767babcff63f3da5fb963b2aecb9d08e177fb0..488e1099cf704aa023bf61b33ce14a385bcd5b73 100644
--- a/Framework/Algorithms/src/ModeratorTzero.cpp
+++ b/Framework/Algorithms/src/ModeratorTzero.cpp
@@ -127,7 +127,7 @@ void ModeratorTzero::exec() {
 
   const size_t numHists = static_cast<size_t>(inputWS->getNumberHistograms());
   Progress prog(this, 0.0, 1.0, numHists); // report progress of algorithm
-  PARALLEL_FOR2(inputWS, outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
   // iterate over the spectra
   for (int i = 0; i < static_cast<int>(numHists); ++i) {
     PARALLEL_START_INTERUPT_REGION
@@ -260,7 +260,7 @@ void ModeratorTzero::execEvent(const std::string &emode) {
   // Loop over the spectra
   const size_t numHists = static_cast<size_t>(outputWS->getNumberHistograms());
   Progress prog(this, 0.0, 1.0, numHists); // report progress of algorithm
-  PARALLEL_FOR1(outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS))
   for (int i = 0; i < static_cast<int>(numHists); ++i) {
     PARALLEL_START_INTERUPT_REGION
     size_t wsIndex = static_cast<size_t>(i);
diff --git a/Framework/Algorithms/src/ModeratorTzeroLinear.cpp b/Framework/Algorithms/src/ModeratorTzeroLinear.cpp
index 6172b04da5cc1fbbd8dfaec9832f409e28bdcfb9..d13fb6f89e68e43efca0fd6ae96c3f2ccdf0cf21 100644
--- a/Framework/Algorithms/src/ModeratorTzeroLinear.cpp
+++ b/Framework/Algorithms/src/ModeratorTzeroLinear.cpp
@@ -131,7 +131,7 @@ void ModeratorTzeroLinear::exec() {
   // do the shift in X
   const size_t numHists = inputWS->getNumberHistograms();
   Progress prog(this, 0.0, 1.0, numHists); // report progress of algorithm
-  PARALLEL_FOR2(inputWS, outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
   for (int i = 0; i < static_cast<int>(numHists); ++i) {
     PARALLEL_START_INTERUPT_REGION
     double t_f, L_i;
@@ -186,7 +186,7 @@ void ModeratorTzeroLinear::execEvent() {
   // Loop over the spectra
   const size_t numHists = outputWS->getNumberHistograms();
   Progress prog(this, 0.0, 1.0, numHists); // report progress of algorithm
-  PARALLEL_FOR1(outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS))
   for (int i = 0; i < static_cast<int>(numHists); ++i) {
     size_t wsIndex = static_cast<size_t>(i);
     PARALLEL_START_INTERUPT_REGION
diff --git a/Framework/Algorithms/src/MonitorEfficiencyCorUser.cpp b/Framework/Algorithms/src/MonitorEfficiencyCorUser.cpp
index ee44bc929d7d0ca6d2ed0781b875a609b960bfd0..bb91bd7c334fb828df65f1eb2db31b6694ec0ffc 100644
--- a/Framework/Algorithms/src/MonitorEfficiencyCorUser.cpp
+++ b/Framework/Algorithms/src/MonitorEfficiencyCorUser.cpp
@@ -72,7 +72,7 @@ void MonitorEfficiencyCorUser::exec() {
 
   // Loop over the histograms (detector spectra)
   double factor = 1 / eff0;
-  PARALLEL_FOR2(m_outputWS, m_inputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*m_outputWS, *m_inputWS))
   for (int64_t i = 0; i < numberOfSpectra_i; ++i) {
     PARALLEL_START_INTERUPT_REGION
     m_outputWS->setHistogram(i, m_inputWS->histogram(i) * factor);
diff --git a/Framework/Algorithms/src/MonteCarloAbsorption.cpp b/Framework/Algorithms/src/MonteCarloAbsorption.cpp
index f63783f3c4efd98a94ce22c5998bacbacd24d5c3..8162a10ff461403d4d707d177b1483450aa92c34 100644
--- a/Framework/Algorithms/src/MonteCarloAbsorption.cpp
+++ b/Framework/Algorithms/src/MonteCarloAbsorption.cpp
@@ -156,7 +156,7 @@ MonteCarloAbsorption::doSimulation(const MatrixWorkspace &inputWS,
   // Configure strategy
   MCAbsorptionStrategy strategy(*beamProfile, inputWS.sample(), nevents);
 
-  PARALLEL_FOR1(outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS))
   for (int64_t i = 0; i < nhists; ++i) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/MultipleScatteringCylinderAbsorption.cpp b/Framework/Algorithms/src/MultipleScatteringCylinderAbsorption.cpp
index 94eee967db754fd23db014f7a54598e3abe7d4a4..b5b36909d76cdb054925789939bfc6a2e967f700 100644
--- a/Framework/Algorithms/src/MultipleScatteringCylinderAbsorption.cpp
+++ b/Framework/Algorithms/src/MultipleScatteringCylinderAbsorption.cpp
@@ -170,7 +170,7 @@ void MultipleScatteringCylinderAbsorption::exec() {
 
     // now do the correction
     const auto &spectrumInfo = out_WSevent->spectrumInfo();
-    PARALLEL_FOR1(out_WSevent)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*out_WSevent))
     for (int64_t index = 0; index < NUM_HIST; ++index) {
       PARALLEL_START_INTERUPT_REGION
       if (!spectrumInfo.hasDetectors(index))
diff --git a/Framework/Algorithms/src/MultiplyRange.cpp b/Framework/Algorithms/src/MultiplyRange.cpp
index 817fd188d40bcd21709f148bdba2081934edc6fa..877e9bb95cbb1194566007a9f8c3424045d9fe2f 100644
--- a/Framework/Algorithms/src/MultiplyRange.cpp
+++ b/Framework/Algorithms/src/MultiplyRange.cpp
@@ -71,7 +71,7 @@ void MultiplyRange::exec() {
   const int histogramCount = static_cast<int>(inputWS->getNumberHistograms());
   Progress progress(this, 0.0, 1.0, histogramCount);
   // Loop over spectra
-  PARALLEL_FOR2(inputWS, outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
   for (int i = 0; i < histogramCount; ++i) {
     PARALLEL_START_INTERUPT_REGION
     outputWS->setHistogram(i, inputWS->histogram(i));
diff --git a/Framework/Algorithms/src/NormaliseByDetector.cpp b/Framework/Algorithms/src/NormaliseByDetector.cpp
index 72fdd3553c08abc3763af3e045452d67b94062bd..059d36f7dea44f594999a04ac8ec9bebf774139e 100644
--- a/Framework/Algorithms/src/NormaliseByDetector.cpp
+++ b/Framework/Algorithms/src/NormaliseByDetector.cpp
@@ -174,7 +174,7 @@ NormaliseByDetector::processHistograms(MatrixWorkspace_sptr inWS) {
   // Choose between parallel execution and sequential execution then, process
   // histograms accordingly.
   if (m_parallelExecution) {
-    PARALLEL_FOR2(inWS, denominatorWS)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*inWS, *denominatorWS))
     for (int wsIndex = 0; wsIndex < static_cast<int>(nHistograms); ++wsIndex) {
       PARALLEL_START_INTERUPT_REGION
       this->processHistogram(wsIndex, inWS, denominatorWS, prog);
diff --git a/Framework/Algorithms/src/Q1DWeighted.cpp b/Framework/Algorithms/src/Q1DWeighted.cpp
index 383f400fbe506ca99a4e09aefc0385a3e7285243..aeb3ac3fdc18b2477d696af6ef6d1dc15877bf79 100644
--- a/Framework/Algorithms/src/Q1DWeighted.cpp
+++ b/Framework/Algorithms/src/Q1DWeighted.cpp
@@ -160,7 +160,7 @@ void Q1DWeighted::exec() {
 
   const auto &spectrumInfo = inputWS->spectrumInfo();
 
-  PARALLEL_FOR2(inputWS, outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
   // Loop over all xLength-1 detector channels
   // Note: xLength -1, because X is a histogram and has a number of boundaries
   // equal to the number of detector channels + 1.
diff --git a/Framework/Algorithms/src/Qxy.cpp b/Framework/Algorithms/src/Qxy.cpp
index 80b5800f8d4782b55d29c2c29d25b0660ab2c6b5..42d8cc3181becfd31991e61c044d7c0fbd8caa71 100644
--- a/Framework/Algorithms/src/Qxy.cpp
+++ b/Framework/Algorithms/src/Qxy.cpp
@@ -126,9 +126,7 @@ void Qxy::exec() {
   // moved to account for the beam centre
   const V3D samplePos = spectrumInfo.samplePosition();
 
-  //  PARALLEL_FOR2(inputWorkspace,outputWorkspace)
   for (int64_t i = 0; i < int64_t(numSpec); ++i) {
-    //    PARALLEL_START_INTERUPT_REGION
     if (!spectrumInfo.hasDetectors(i)) {
       g_log.warning() << "Workspace index " << i
                       << " has no detector assigned to it - discarding\n";
@@ -298,9 +296,7 @@ void Qxy::exec() {
 
     prog.report("Calculating Q");
 
-    //    PARALLEL_END_INTERUPT_REGION
   } // loop over all spectra
-  //  PARALLEL_CHECK_INTERUPT_REGION
 
   // take sqrt of error weight values
   // left to be executed here for computational efficiency
diff --git a/Framework/Algorithms/src/Rebin.cpp b/Framework/Algorithms/src/Rebin.cpp
index 64388219cf672609ce49b1fa75e91050dc9d1845..b7729e557e175382029736f6cf3648e65c2e97db 100644
--- a/Framework/Algorithms/src/Rebin.cpp
+++ b/Framework/Algorithms/src/Rebin.cpp
@@ -247,7 +247,7 @@ void Rebin::exec() {
       outputWS->replaceAxis(1, inputWS->getAxis(1)->clone(outputWS.get()));
 
     Progress prog(this, 0.0, 1.0, histnumber);
-    PARALLEL_FOR2(inputWS, outputWS)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
     for (int hist = 0; hist < histnumber; ++hist) {
       PARALLEL_START_INTERUPT_REGION
       // get const references to input Workspace arrays (no copying)
diff --git a/Framework/Algorithms/src/Rebin2D.cpp b/Framework/Algorithms/src/Rebin2D.cpp
index 96b916eecd3ee6157cba4ba8de48adbdd64b2bdb..ef9247399c74e18b86d4314f75b1746665b2d7c5 100644
--- a/Framework/Algorithms/src/Rebin2D.cpp
+++ b/Framework/Algorithms/src/Rebin2D.cpp
@@ -119,7 +119,7 @@ void Rebin2D::exec() {
   m_progress = boost::shared_ptr<API::Progress>(
       new API::Progress(this, 0.0, 1.0, nreports));
 
-  PARALLEL_FOR2(inputWS, outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
   for (int64_t i = 0; i < static_cast<int64_t>(numYBins);
        ++i) // signed for openmp
   {
diff --git a/Framework/Algorithms/src/RebinByPulseTimes.cpp b/Framework/Algorithms/src/RebinByPulseTimes.cpp
index c52af972f76e499fc824d228c811eb954eeb1ec2..6067ef19a528fb77dabb56e6f105d6caa690d780 100644
--- a/Framework/Algorithms/src/RebinByPulseTimes.cpp
+++ b/Framework/Algorithms/src/RebinByPulseTimes.cpp
@@ -50,7 +50,7 @@ void RebinByPulseTimes::doHistogramming(IEventWorkspace_sptr inWS,
 
   auto x = Kernel::make_cow<HistogramData::HistogramX>(OutXValues_scaled);
 
-  PARALLEL_FOR2(inWS, outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inWS, *outputWS))
   for (int i = 0; i < histnumber; ++i) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/RebinByTimeAtSample.cpp b/Framework/Algorithms/src/RebinByTimeAtSample.cpp
index e6e0446fb71423397eaf4e00d9ab24dde13efb1e..cb0445ff59620e2e657256eccbde3956477c3e50 100644
--- a/Framework/Algorithms/src/RebinByTimeAtSample.cpp
+++ b/Framework/Algorithms/src/RebinByTimeAtSample.cpp
@@ -61,7 +61,7 @@ void RebinByTimeAtSample::doHistogramming(IEventWorkspace_sptr inWS,
   auto x = Kernel::make_cow<HistogramData::HistogramX>(OutXValues_scaled);
 
   // Go through all the histograms and set the data
-  PARALLEL_FOR2(inWS, outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inWS, *outputWS))
   for (int i = 0; i < histnumber; ++i) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/Rebunch.cpp b/Framework/Algorithms/src/Rebunch.cpp
index 8dfa51318a3a622a9b4b600db4e88fd76644ca66..a7fb024357b83ea2fd0caa1d28b7c8aff511a834 100644
--- a/Framework/Algorithms/src/Rebunch.cpp
+++ b/Framework/Algorithms/src/Rebunch.cpp
@@ -79,7 +79,7 @@ void Rebunch::exec() {
   int progress_step = histnumber / 100;
   if (progress_step == 0)
     progress_step = 1;
-  PARALLEL_FOR2(inputW, outputW)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputW, *outputW))
   for (int hist = 0; hist < histnumber; hist++) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/RemoveBackground.cpp b/Framework/Algorithms/src/RemoveBackground.cpp
index b63c667ae70fb79c03c891c067b3e7b6b739d8df..ffc385ad655973ebfbd0d832acd1b1a07b1fbb0b 100644
--- a/Framework/Algorithms/src/RemoveBackground.cpp
+++ b/Framework/Algorithms/src/RemoveBackground.cpp
@@ -116,7 +116,7 @@ void RemoveBackground::exec() {
                                 inPlace, nullifyNegative);
 
   Progress prog(this, 0.0, 1.0, histnumber);
-  PARALLEL_FOR2(inputWS, outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
   for (int hist = 0; hist < histnumber; ++hist) {
     PARALLEL_START_INTERUPT_REGION
     // get references to output Workspace X-arrays.
diff --git a/Framework/Algorithms/src/RemoveExpDecay.cpp b/Framework/Algorithms/src/RemoveExpDecay.cpp
index f8b53edef4b398246da4f1d31d4f763fd4408d17..c8e7dc7edffda76138828b1b7d2f6d65ca2f7c58 100644
--- a/Framework/Algorithms/src/RemoveExpDecay.cpp
+++ b/Framework/Algorithms/src/RemoveExpDecay.cpp
@@ -77,7 +77,7 @@ void MuonRemoveExpDecay::exec() {
   if (inputWS != outputWS) {
 
     // Copy all the Y and E data
-    PARALLEL_FOR2(inputWS, outputWS)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
     for (int64_t i = 0; i < int64_t(numSpectra); ++i) {
       PARALLEL_START_INTERUPT_REGION
       const auto index = static_cast<size_t>(i);
@@ -91,7 +91,7 @@ void MuonRemoveExpDecay::exec() {
 
   // Do the specified spectra only
   int specLength = static_cast<int>(spectra.size());
-  PARALLEL_FOR2(inputWS, outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
   for (int i = 0; i < specLength; ++i) {
     PARALLEL_START_INTERUPT_REGION
     const auto specNum = static_cast<size_t>(spectra[i]);
diff --git a/Framework/Algorithms/src/ResampleX.cpp b/Framework/Algorithms/src/ResampleX.cpp
index b40eaafee21e290e38ce66b0f60c7772aa4477b6..9e13cdf4d6e483f9aa187ee89ab03b8a6540e016 100644
--- a/Framework/Algorithms/src/ResampleX.cpp
+++ b/Framework/Algorithms/src/ResampleX.cpp
@@ -340,7 +340,7 @@ void ResampleX::exec() {
         Progress prog(this, 0.0, 1.0, numSpectra);
 
         // do the rebinning
-        PARALLEL_FOR2(inputEventWS, outputWS)
+        PARALLEL_FOR_IF(Kernel::threadSafe(*inputEventWS, *outputWS))
         for (int wkspIndex = 0; wkspIndex < numSpectra; ++wkspIndex) {
           PARALLEL_START_INTERUPT_REGION
           BinEdges xValues(0);
@@ -375,7 +375,7 @@ void ResampleX::exec() {
       Progress prog(this, 0.0, 1.0, numSpectra);
 
       // Go through all the histograms and set the data
-      PARALLEL_FOR2(inputEventWS, outputWS)
+      PARALLEL_FOR_IF(Kernel::threadSafe(*inputEventWS, *outputWS))
       for (int wkspIndex = 0; wkspIndex < numSpectra; ++wkspIndex) {
         PARALLEL_START_INTERUPT_REGION
 
@@ -445,7 +445,7 @@ void ResampleX::exec() {
       outputWS->replaceAxis(1, inputWS->getAxis(1)->clone(outputWS.get()));
 
     Progress prog(this, 0.0, 1.0, numSpectra);
-    PARALLEL_FOR2(inputWS, outputWS)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
     for (int wkspIndex = 0; wkspIndex < numSpectra; ++wkspIndex) {
       PARALLEL_START_INTERUPT_REGION
       // get const references to input Workspace arrays (no copying)
diff --git a/Framework/Algorithms/src/ResetNegatives.cpp b/Framework/Algorithms/src/ResetNegatives.cpp
index c1e0ce16578e9f3c04cd3a87c725c51bbb25f2ae..5cb840c07c26989c10f8a930c35c0d44e5636639 100644
--- a/Framework/Algorithms/src/ResetNegatives.cpp
+++ b/Framework/Algorithms/src/ResetNegatives.cpp
@@ -93,7 +93,7 @@ void ResetNegatives::exec() {
 
   // generate output workspace - copy X and dY
   outputWS = API::WorkspaceFactory::Instance().create(inputWS);
-  PARALLEL_FOR2(inputWS, outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
   for (int64_t i = 0; i < nHist; i++) {
     PARALLEL_START_INTERUPT_REGION
     const auto index = static_cast<size_t>(i);
@@ -138,7 +138,7 @@ inline double fixZero(const double value) {
 void ResetNegatives::pushMinimum(MatrixWorkspace_const_sptr minWS,
                                  MatrixWorkspace_sptr wksp, Progress &prog) {
   int64_t nHist = minWS->getNumberHistograms();
-  PARALLEL_FOR2(wksp, minWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*wksp, *minWS))
   for (int64_t i = 0; i < nHist; i++) {
     PARALLEL_START_INTERUPT_REGION
     double minValue = minWS->y(i)[0];
@@ -170,7 +170,7 @@ void ResetNegatives::changeNegatives(MatrixWorkspace_const_sptr minWS,
                                      MatrixWorkspace_sptr wksp,
                                      Progress &prog) {
   int64_t nHist = wksp->getNumberHistograms();
-  PARALLEL_FOR2(minWS, wksp)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*minWS, *wksp))
   for (int64_t i = 0; i < nHist; i++) {
     PARALLEL_START_INTERUPT_REGION
     if (minWS->y(i)[0] <=
diff --git a/Framework/Algorithms/src/SampleCorrections/MayersSampleCorrection.cpp b/Framework/Algorithms/src/SampleCorrections/MayersSampleCorrection.cpp
index 6fc6b35dd61ccec6edcdcc29a73fcebd2f03a9bb..f4e353e0230cc2d8ad3768ee8939f93842c285c1 100644
--- a/Framework/Algorithms/src/SampleCorrections/MayersSampleCorrection.cpp
+++ b/Framework/Algorithms/src/SampleCorrections/MayersSampleCorrection.cpp
@@ -101,7 +101,7 @@ void MayersSampleCorrection::exec() {
   Progress prog(this, 0., 1., nhist);
   prog.setNotifyStep(0.01);
 
-  PARALLEL_FOR2(inputWS, outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
   for (int64_t i = 0; i < static_cast<int64_t>(nhist); ++i) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/ScaleX.cpp b/Framework/Algorithms/src/ScaleX.cpp
index f80e8b9b6b18b7fe2e854ec0661add9efcd0e6b4..a0a1345b621eed69868a93e5650e5a6aabf2e839 100644
--- a/Framework/Algorithms/src/ScaleX.cpp
+++ b/Framework/Algorithms/src/ScaleX.cpp
@@ -123,7 +123,7 @@ void ScaleX::exec() {
   }
 
   // do the shift in X
-  PARALLEL_FOR2(inputW, outputW)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputW, *outputW))
   for (int i = 0; i < histnumber; ++i) {
     PARALLEL_START_INTERUPT_REGION
 
@@ -182,7 +182,7 @@ void ScaleX::execEvent() {
 
   const std::string op = getPropertyValue("Operation");
   int numHistograms = static_cast<int>(outputWS->getNumberHistograms());
-  PARALLEL_FOR1(outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS))
   for (int i = 0; i < numHistograms; ++i) {
     PARALLEL_START_INTERUPT_REGION
     // Do the offsetting
diff --git a/Framework/Algorithms/src/SetUncertainties.cpp b/Framework/Algorithms/src/SetUncertainties.cpp
index 7a4d49ddefb6fb573d35153153be2e2d679e64dc..aa0c7fab054782de81226ef11e4106d694147ab0 100644
--- a/Framework/Algorithms/src/SetUncertainties.cpp
+++ b/Framework/Algorithms/src/SetUncertainties.cpp
@@ -93,7 +93,7 @@ void SetUncertainties::exec() {
   const size_t numHists = inputWorkspace->getNumberHistograms();
   Progress prog(this, 0.0, 1.0, numHists);
 
-  PARALLEL_FOR2(inputWorkspace, outputWorkspace)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWorkspace, *outputWorkspace))
   for (int64_t i = 0; i < int64_t(numHists); ++i) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/SmoothData.cpp b/Framework/Algorithms/src/SmoothData.cpp
index 83c4d93f2fd211d5a70f8c5ee8a853314bb945ee..00c6f94f9793087aabd45fd470a811792f6b6df1 100644
--- a/Framework/Algorithms/src/SmoothData.cpp
+++ b/Framework/Algorithms/src/SmoothData.cpp
@@ -62,7 +62,7 @@ void SmoothData::exec() {
       WorkspaceFactory::Instance().create(inputWorkspace);
 
   Progress progress(this, 0.0, 1.0, inputWorkspace->getNumberHistograms());
-  PARALLEL_FOR2(inputWorkspace, outputWorkspace)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWorkspace, *outputWorkspace))
   // Loop over all the spectra in the workspace
   for (int i = 0; i < static_cast<int>(inputWorkspace->getNumberHistograms());
        ++i) {
diff --git a/Framework/Algorithms/src/SmoothNeighbours.cpp b/Framework/Algorithms/src/SmoothNeighbours.cpp
index 6c2e1a6b4a7be0c12c94064926b9fee1ce6cba7b..d6e4fd4f366357ead20ebd99b0d944d3b8c37b9f 100644
--- a/Framework/Algorithms/src/SmoothNeighbours.cpp
+++ b/Framework/Algorithms/src/SmoothNeighbours.cpp
@@ -615,7 +615,7 @@ void SmoothNeighbours::execWorkspace2D() {
   // API::WorkspaceFactory::Instance().initializeFromParent(inWS, outWS, false);
 
   // Go through all the output workspace
-  PARALLEL_FOR2(inWS, outWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inWS, *outWS))
   for (int outWIi = 0; outWIi < int(numberOfSpectra); outWIi++) {
     PARALLEL_START_INTERUPT_REGION
 
@@ -767,7 +767,7 @@ void SmoothNeighbours::execEvent(Mantid::DataObjects::EventWorkspace_sptr ws) {
                     boost::dynamic_pointer_cast<MatrixWorkspace>(outWS));
 
   // Go through all the output workspace
-  PARALLEL_FOR2(ws, outWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*ws, *outWS))
   for (int outWIi = 0; outWIi < int(numberOfSpectra); outWIi++) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/SofQWNormalisedPolygon.cpp b/Framework/Algorithms/src/SofQWNormalisedPolygon.cpp
index 458e73188cab958fc7e362efac22a97118460809..c3873ef09adceebf2d1214555028dc424c7ca503 100644
--- a/Framework/Algorithms/src/SofQWNormalisedPolygon.cpp
+++ b/Framework/Algorithms/src/SofQWNormalisedPolygon.cpp
@@ -99,7 +99,7 @@ void SofQWNormalisedPolygon::exec() {
   const auto &X = inputWS->x(0);
   int emode = m_EmodeProperties.m_emode;
 
-  PARALLEL_FOR2(inputWS, outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
   for (int64_t i = 0; i < static_cast<int64_t>(nHistos);
        ++i) // signed for openmp
   {
diff --git a/Framework/Algorithms/src/SofQWPolygon.cpp b/Framework/Algorithms/src/SofQWPolygon.cpp
index 1e55da72e902f8f5191a2f815f24a9dec3fdfdba..cceb6adc60c0551f484b3fe66705e59d2dd0ed47 100644
--- a/Framework/Algorithms/src/SofQWPolygon.cpp
+++ b/Framework/Algorithms/src/SofQWPolygon.cpp
@@ -75,7 +75,7 @@ void SofQWPolygon::exec() {
     qCalculator = &SofQWPolygon::calculateIndirectQ;
   }
 
-  PARALLEL_FOR2(inputWS, outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
   for (int64_t i = 0; i < static_cast<int64_t>(nTheta);
        ++i) // signed for openmp
   {
diff --git a/Framework/Algorithms/src/SolidAngle.cpp b/Framework/Algorithms/src/SolidAngle.cpp
index 781a531003b1bc4cd82c9503e1367007f0109ffa..3ac5ef22f01db98cb35735fa0cb0f937de1d6553 100644
--- a/Framework/Algorithms/src/SolidAngle.cpp
+++ b/Framework/Algorithms/src/SolidAngle.cpp
@@ -89,7 +89,7 @@ void SolidAngle::exec() {
   Progress prog(this, 0.0, 1.0, numberOfSpectra);
 
   // Loop over the histograms (detector spectra)
-  PARALLEL_FOR2(outputWS, inputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS, *inputWS))
   for (int j = 0; j <= loopIterations; ++j) {
     PARALLEL_START_INTERUPT_REGION
     int i = j + m_MinSpec;
diff --git a/Framework/Algorithms/src/Stitch1D.cpp b/Framework/Algorithms/src/Stitch1D.cpp
index 247c7dc2166e81c8a14e8439537fc3ff4d465c26..05d0eef3bc2cfa36da72ef8eca1b937247762a4f 100644
--- a/Framework/Algorithms/src/Stitch1D.cpp
+++ b/Framework/Algorithms/src/Stitch1D.cpp
@@ -60,7 +60,7 @@ MatrixWorkspace_sptr Stitch1D::maskAllBut(int a1, int a2,
                                           MatrixWorkspace_sptr &source) {
   MatrixWorkspace_sptr product = WorkspaceFactory::Instance().create(source);
   const int histogramCount = static_cast<int>(source->getNumberHistograms());
-  PARALLEL_FOR2(source, product)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*source, *product))
   for (int i = 0; i < histogramCount; ++i) {
     PARALLEL_START_INTERUPT_REGION
     // Copy over the bin boundaries
@@ -97,7 +97,7 @@ MatrixWorkspace_sptr Stitch1D::maskAllBut(int a1, int a2,
  */
 void Stitch1D::maskInPlace(int a1, int a2, MatrixWorkspace_sptr source) {
   const int histogramCount = static_cast<int>(source->getNumberHistograms());
-  PARALLEL_FOR1(source)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*source))
   for (int i = 0; i < histogramCount; ++i) {
     PARALLEL_START_INTERUPT_REGION
     // Copy over the data
@@ -300,7 +300,7 @@ MatrixWorkspace_sptr Stitch1D::rebin(MatrixWorkspace_sptr &input,
 
   // Record special values and then mask them out as zeros. Special values are
   // remembered and then replaced post processing.
-  PARALLEL_FOR1(outWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*outWS))
   for (int i = 0; i < histogramCount; ++i) {
     PARALLEL_START_INTERUPT_REGION
     std::vector<size_t> &nanYIndexes = m_nanYIndexes[i];
@@ -461,7 +461,7 @@ Stitch1D::findStartEndIndexes(double startOverlap, double endOverlap,
 bool Stitch1D::hasNonzeroErrors(MatrixWorkspace_sptr ws) {
   int64_t ws_size = static_cast<int64_t>(ws->getNumberHistograms());
   bool hasNonZeroErrors = false;
-  PARALLEL_FOR1(ws)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*ws))
   for (int i = 0; i < ws_size; ++i) {
     PARALLEL_START_INTERUPT_REGION
     if (!hasNonZeroErrors) // Keep checking
@@ -625,7 +625,7 @@ void Stitch1D::exec() {
  */
 void Stitch1D::reinsertSpecialValues(MatrixWorkspace_sptr ws) {
   int histogramCount = static_cast<int>(ws->getNumberHistograms());
-  PARALLEL_FOR1(ws)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*ws))
   for (int i = 0; i < histogramCount; ++i) {
     PARALLEL_START_INTERUPT_REGION
     // Copy over the data
diff --git a/Framework/Algorithms/src/SumEventsByLogValue.cpp b/Framework/Algorithms/src/SumEventsByLogValue.cpp
index 455a5f94f0692964227a85a8155d5a3d8d52160d..e344885039b609ef9e7d3c2f168ade6299099383 100644
--- a/Framework/Algorithms/src/SumEventsByLogValue.cpp
+++ b/Framework/Algorithms/src/SumEventsByLogValue.cpp
@@ -148,7 +148,7 @@ void SumEventsByLogValue::createTableOutput(
   std::vector<int> Y(xLength);
   const int numSpec = static_cast<int>(m_inputWorkspace->getNumberHistograms());
   Progress prog(this, 0.0, 1.0, numSpec + xLength);
-  PARALLEL_FOR1(m_inputWorkspace)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*m_inputWorkspace))
   for (int spec = 0; spec < numSpec; ++spec) {
     PARALLEL_START_INTERUPT_REGION
     const IEventList &eventList = m_inputWorkspace->getSpectrum(spec);
@@ -418,7 +418,7 @@ void SumEventsByLogValue::createBinnedOutput(
   auto &Y = outputWorkspace->mutableY(0);
   const int numSpec = static_cast<int>(m_inputWorkspace->getNumberHistograms());
   Progress prog(this, 0.0, 1.0, numSpec);
-  PARALLEL_FOR1(m_inputWorkspace)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*m_inputWorkspace))
   for (int spec = 0; spec < numSpec; ++spec) {
     PARALLEL_START_INTERUPT_REGION
     const IEventList &eventList = m_inputWorkspace->getSpectrum(spec);
diff --git a/Framework/Algorithms/src/TOFSANSResolution.cpp b/Framework/Algorithms/src/TOFSANSResolution.cpp
index b06a7933b76fa8bee33743e1dd1fd373e5b542d0..7352ada947d4bcb26dc0e669bf32f1a847a8ee51 100644
--- a/Framework/Algorithms/src/TOFSANSResolution.cpp
+++ b/Framework/Algorithms/src/TOFSANSResolution.cpp
@@ -131,7 +131,7 @@ void TOFSANSResolution::exec() {
   const auto &spectrumInfo = reducedWS->spectrumInfo();
   const double L1 = spectrumInfo.l1();
 
-  PARALLEL_FOR2(reducedWS, iqWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*reducedWS, *iqWS))
   for (int i = 0; i < numberOfSpectra; i++) {
     PARALLEL_START_INTERUPT_REGION
     if (!spectrumInfo.hasDetectors(i)) {
diff --git a/Framework/Algorithms/src/Transpose.cpp b/Framework/Algorithms/src/Transpose.cpp
index bff41ff5a360e77a008b4a3a9c86a58d1d81c2d5..d1bcfec083481246313e2511249d9a127ddfccb3 100644
--- a/Framework/Algorithms/src/Transpose.cpp
+++ b/Framework/Algorithms/src/Transpose.cpp
@@ -53,8 +53,7 @@ void Transpose::exec() {
 
   Progress progress(this, 0.0, 1.0, newNhist * newYsize);
   progress.report("Swapping data values");
-
-  PARALLEL_FOR2(inputWorkspace, outputWorkspace)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWorkspace, *outputWorkspace))
   for (int64_t i = 0; i < static_cast<int64_t>(newNhist); ++i) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Algorithms/src/UnaryOperation.cpp b/Framework/Algorithms/src/UnaryOperation.cpp
index 966dd3ec4f31055f5cde0e07c0a56e97152ed6df..b86d001ec92df7cef302e5528f8f521656ac33cc 100644
--- a/Framework/Algorithms/src/UnaryOperation.cpp
+++ b/Framework/Algorithms/src/UnaryOperation.cpp
@@ -77,7 +77,7 @@ void UnaryOperation::exec() {
 
   // Loop over every cell in the workspace, calling the abstract correction
   // function
-  PARALLEL_FOR2(in_work, out_work)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*in_work, *out_work))
   for (int64_t i = 0; i < int64_t(numSpec); ++i) {
     PARALLEL_START_INTERUPT_REGION
     // Copy the X values over
@@ -121,7 +121,7 @@ void UnaryOperation::execEvent() {
 
   int64_t numHistograms = static_cast<int64_t>(outputWS->getNumberHistograms());
   API::Progress prog = API::Progress(this, 0.0, 1.0, numHistograms);
-  PARALLEL_FOR1(outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS))
   for (int64_t i = 0; i < numHistograms; ++i) {
     PARALLEL_START_INTERUPT_REGION
     // switch to weighted events if needed, and use the appropriate helper
diff --git a/Framework/Algorithms/src/UnwrapSNS.cpp b/Framework/Algorithms/src/UnwrapSNS.cpp
index 765a29ba2e1939a423a3daaae4f12561dba1d31b..bb222f377302d1bc01bebbade2aa7dfa4eb54c6f 100644
--- a/Framework/Algorithms/src/UnwrapSNS.cpp
+++ b/Framework/Algorithms/src/UnwrapSNS.cpp
@@ -136,7 +136,7 @@ void UnwrapSNS::exec() {
   const auto &spectrumInfo = m_inputWS->spectrumInfo();
   const double L1 = spectrumInfo.l1();
 
-  PARALLEL_FOR2(m_inputWS, outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*m_inputWS, *outputWS))
   for (int workspaceIndex = 0; workspaceIndex < m_numberOfSpectra;
        workspaceIndex++) {
     PARALLEL_START_INTERUPT_REGION
@@ -205,10 +205,8 @@ void UnwrapSNS::execEvent() {
   const double L1 = spectrumInfo.l1();
 
   // do the actual work
-  //  PARALLEL_FOR2(m_inputWS, outW)
   for (int workspaceIndex = 0; workspaceIndex < m_numberOfSpectra;
        workspaceIndex++) {
-    //    PARALLEL_START_INTERUPT_REGION
     std::size_t numEvents = outW->getSpectrum(workspaceIndex).getNumberEvents();
     double Ld = -1.0;
     if (spectrumInfo.hasDetectors(workspaceIndex))
@@ -234,9 +232,7 @@ void UnwrapSNS::execEvent() {
       outW->getSpectrum(workspaceIndex).setTofs(times);
     }
     m_progress->report();
-    //    PARALLEL_END_INTERUPT_REGION
   }
-  //  PARALLEL_CHECK_INTERUPT_REGION
 
   outW->clearMRU();
   this->runMaskDetectors();
diff --git a/Framework/Algorithms/src/WorkspaceJoiners.cpp b/Framework/Algorithms/src/WorkspaceJoiners.cpp
index 77730df8a6526bcf32a0138b5de97b4e2535a6e9..9f63c08f5cb7496e0827d7c2f80cb66583a0db07 100644
--- a/Framework/Algorithms/src/WorkspaceJoiners.cpp
+++ b/Framework/Algorithms/src/WorkspaceJoiners.cpp
@@ -50,7 +50,7 @@ WorkspaceJoiners::execWS2D(API::MatrixWorkspace_const_sptr ws1,
 
   // Loop over the input workspaces in turn copying the data into the output one
   const int64_t &nhist1 = ws1->getNumberHistograms();
-  PARALLEL_FOR2(ws1, output)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*ws1, *output))
   for (int64_t i = 0; i < nhist1; ++i) {
     PARALLEL_START_INTERUPT_REGION
     auto &outSpec = output->getSpectrum(i);
@@ -75,7 +75,7 @@ WorkspaceJoiners::execWS2D(API::MatrixWorkspace_const_sptr ws1,
   // For second loop we use the offset from the first
   const int64_t &nhist2 = ws2->getNumberHistograms();
   const auto &spectrumInfo = ws2->spectrumInfo();
-  PARALLEL_FOR2(ws2, output)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*ws2, *output))
   for (int64_t j = 0; j < nhist2; ++j) {
     PARALLEL_START_INTERUPT_REGION
     // The spectrum in the output workspace
diff --git a/Framework/Algorithms/src/XDataConverter.cpp b/Framework/Algorithms/src/XDataConverter.cpp
index 0e893f24b4234b71befd3205ee5c90fc9778284f..290cc76352d54e4296c88947d73205f3727f8309 100644
--- a/Framework/Algorithms/src/XDataConverter.cpp
+++ b/Framework/Algorithms/src/XDataConverter.cpp
@@ -61,7 +61,7 @@ void XDataConverter::exec() {
     outputWS->replaceAxis(1, inputWS->getAxis(1)->clone(outputWS.get()));
 
   Progress prog(this, 0.0, 1.0, numSpectra);
-  PARALLEL_FOR2(inputWS, outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
   for (int i = 0; i < int(numSpectra); ++i) {
     PARALLEL_START_INTERUPT_REGION
 
@@ -102,4 +102,4 @@ void XDataConverter::setXData(API::MatrixWorkspace_sptr outputWS,
   }
 }
 }
-}
\ No newline at end of file
+}
diff --git a/Framework/Crystal/src/AnvredCorrection.cpp b/Framework/Crystal/src/AnvredCorrection.cpp
index 5d1051413c2fa2168a1817e869065b42a91d30f5..48516938fb85b747bb99fefd5f5c1413c7281c6f 100644
--- a/Framework/Crystal/src/AnvredCorrection.cpp
+++ b/Framework/Crystal/src/AnvredCorrection.cpp
@@ -168,7 +168,7 @@ void AnvredCorrection::exec() {
 
   Progress prog(this, 0.0, 1.0, numHists);
   // Loop over the spectra
-  PARALLEL_FOR2(m_inputWS, correctionFactors)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*m_inputWS, *correctionFactors))
   for (int64_t i = 0; i < int64_t(numHists); ++i) {
     PARALLEL_START_INTERUPT_REGION
 
@@ -288,7 +288,7 @@ void AnvredCorrection::execEvent() {
 
   Progress prog(this, 0.0, 1.0, numHists);
   // Loop over the spectra
-  PARALLEL_FOR2(eventW, correctionFactors)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*eventW, *correctionFactors))
   for (int64_t i = 0; i < int64_t(numHists); ++i) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/Crystal/src/CentroidPeaks.cpp b/Framework/Crystal/src/CentroidPeaks.cpp
index 4f1b59c355f5a739fde98cf74e4a441ea034609d..cc719aa2e16f677f56bc5b3a92bd3a66ee2c763e 100644
--- a/Framework/Crystal/src/CentroidPeaks.cpp
+++ b/Framework/Crystal/src/CentroidPeaks.cpp
@@ -87,7 +87,7 @@ void CentroidPeaks::integrate() {
 
   int Edge = getProperty("EdgePixels");
   Progress prog(this, MinPeaks, 1.0, MaxPeaks);
-  PARALLEL_FOR2(inWS, peakWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inWS, *peakWS))
   for (int i = MinPeaks; i <= MaxPeaks; ++i) {
     PARALLEL_START_INTERUPT_REGION
     // Get a direct ref to that peak.
@@ -224,7 +224,7 @@ void CentroidPeaks::integrateEvent() {
 
   int Edge = getProperty("EdgePixels");
   Progress prog(this, MinPeaks, 1.0, MaxPeaks);
-  PARALLEL_FOR2(inWS, peakWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inWS, *peakWS))
   for (int i = MinPeaks; i <= MaxPeaks; ++i) {
     PARALLEL_START_INTERUPT_REGION
     // Get a direct ref to that peak.
diff --git a/Framework/Crystal/src/FindSXPeaks.cpp b/Framework/Crystal/src/FindSXPeaks.cpp
index e5c0d6ecec38f8d7dfe27ad506b413656d306e10..cb45e50b5d625c9ddd7339ff8760002013c4e27a 100644
--- a/Framework/Crystal/src/FindSXPeaks.cpp
+++ b/Framework/Crystal/src/FindSXPeaks.cpp
@@ -109,7 +109,7 @@ void FindSXPeaks::exec() {
   // unlikely to have more than this.
   entries.reserve(1000);
   // Count the peaks so that we can resize the peakvector at the end.
-  PARALLEL_FOR1(localworkspace)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*localworkspace))
   for (int i = static_cast<int>(m_MinSpec); i <= static_cast<int>(m_MaxSpec);
        ++i) {
     PARALLEL_START_INTERUPT_REGION
diff --git a/Framework/Crystal/src/IntegratePeaksHybrid.cpp b/Framework/Crystal/src/IntegratePeaksHybrid.cpp
index 7ded565d9a703f00a78bcecc81a4efab9a34c9b3..7ccdfd32d0ee4f652fafde3f7edd23b78420b503 100644
--- a/Framework/Crystal/src/IntegratePeaksHybrid.cpp
+++ b/Framework/Crystal/src/IntegratePeaksHybrid.cpp
@@ -173,7 +173,7 @@ void IntegratePeaksHybrid::exec() {
 
   Progress progress(this, 0, 1, peakWS->getNumberPeaks());
 
-  PARALLEL_FOR1(peakWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*peakWS))
   for (int i = 0; i < peakWS->getNumberPeaks(); ++i) {
 
     PARALLEL_START_INTERUPT_REGION
diff --git a/Framework/Crystal/src/IntegratePeaksUsingClusters.cpp b/Framework/Crystal/src/IntegratePeaksUsingClusters.cpp
index 52f214f8939aa11e7707f553ba50f67c2967b927..7836bb9e2d7ac456a48f5812904c22b39a9d3d5e 100644
--- a/Framework/Crystal/src/IntegratePeaksUsingClusters.cpp
+++ b/Framework/Crystal/src/IntegratePeaksUsingClusters.cpp
@@ -150,7 +150,7 @@ void IntegratePeaksUsingClusters::exec() {
   progress.doReport("Performing Peak Integration");
   g_log.information("Starting Integration");
   progress.resetNumSteps(peakWS->getNumberPeaks(), 0.9, 1);
-  PARALLEL_FOR1(peakWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*peakWS))
   for (int i = 0; i < peakWS->getNumberPeaks(); ++i) {
     PARALLEL_START_INTERUPT_REGION
     Geometry::IPeak &peak = peakWS->getPeak(i);
diff --git a/Framework/Crystal/src/NormaliseVanadium.cpp b/Framework/Crystal/src/NormaliseVanadium.cpp
index 08cbfb59c94a774aa0e65298f563e68980e1844d..26423a185c5ccc6cfb9af46a6da47282cb50f6a3 100644
--- a/Framework/Crystal/src/NormaliseVanadium.cpp
+++ b/Framework/Crystal/src/NormaliseVanadium.cpp
@@ -67,7 +67,7 @@ void NormaliseVanadium::exec() {
 
   Progress prog(this, 0.0, 1.0, numHists);
   // Loop over the spectra
-  PARALLEL_FOR2(m_inputWS, correctionFactors)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*m_inputWS, *correctionFactors))
   for (int64_t i = 0; i < int64_t(numHists); ++i) {
     //    PARALLEL_START_INTERUPT_REGION //FIXME: Restore
 
diff --git a/Framework/Crystal/src/PeaksIntersection.cpp b/Framework/Crystal/src/PeaksIntersection.cpp
index 833a6d292a9b35b48db4d4edf4778d0dd4d2cf87..6051f4e2f2eb80051e983c5c6b8efed4fcb66414 100644
--- a/Framework/Crystal/src/PeaksIntersection.cpp
+++ b/Framework/Crystal/src/PeaksIntersection.cpp
@@ -104,7 +104,7 @@ void PeaksIntersection::executePeaksIntersection(const bool checkPeakExtents) {
   }
   Progress prog(this, 0, 1, 100);
 
-  PARALLEL_FOR2(ws, outputWorkspace)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*ws, *outputWorkspace))
   for (int i = 0; i < nPeaks; ++i) {
     PARALLEL_START_INTERUPT_REGION
     IPeak *peak = ws->getPeakPtr(i);
diff --git a/Framework/Crystal/src/SCDCalibratePanels.cpp b/Framework/Crystal/src/SCDCalibratePanels.cpp
index dd4f20a3a277ee2c58a1d9a2178d013ca2bd80ad..fdab7d7c0c83754edb0596d21fb2df1e09346f79 100644
--- a/Framework/Crystal/src/SCDCalibratePanels.cpp
+++ b/Framework/Crystal/src/SCDCalibratePanels.cpp
@@ -173,7 +173,7 @@ void SCDCalibratePanels::exec() {
   std::vector<std::string> fit_workspaces(MyBankNames.size(), "fit_");
   std::vector<std::string> parameter_workspaces(MyBankNames.size(), "params_");
 
-  PARALLEL_FOR1(peaksWs)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*peaksWs))
   for (int i = 0; i < static_cast<int>(MyBankNames.size()); ++i) {
     PARALLEL_START_INTERUPT_REGION
     const std::string &iBank = *std::next(MyBankNames.begin(), i);
@@ -331,7 +331,7 @@ void SCDCalibratePanels::exec() {
       boost::const_pointer_cast<Geometry::Instrument>(peaksWs->getInstrument());
   Geometry::OrientedLattice lattice0 =
       peaksWs->mutableSample().getOrientedLattice();
-  PARALLEL_FOR1(peaksWs)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*peaksWs))
   for (int i = 0; i < nPeaks; i++) {
     PARALLEL_START_INTERUPT_REGION
     DataObjects::Peak &peak = peaksWs->getPeak(i);
diff --git a/Framework/CurveFitting/src/Algorithms/ConvertToYSpace.cpp b/Framework/CurveFitting/src/Algorithms/ConvertToYSpace.cpp
index 51f88a0f47b3cbf797ed44243c000659f93a46e4..fe99f86d374f3b2d78dec702556751f331c83415 100644
--- a/Framework/CurveFitting/src/Algorithms/ConvertToYSpace.cpp
+++ b/Framework/CurveFitting/src/Algorithms/ConvertToYSpace.cpp
@@ -196,7 +196,7 @@ void ConvertToYSpace::exec() {
   const int64_t nreports = nhist;
   auto progress = boost::make_shared<Progress>(this, 0.0, 1.0, nreports);
 
-  PARALLEL_FOR2(m_inputWS, m_outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*m_inputWS, *m_outputWS))
   for (int64_t i = 0; i < nhist; ++i) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/DataHandling/src/CreateSimulationWorkspace.cpp b/Framework/DataHandling/src/CreateSimulationWorkspace.cpp
index 9a724205cf3cc622c2d93864d61be98548eb2369..34b5c1f92848cdc7f770087f09a16ee3e5f47a2a 100644
--- a/Framework/DataHandling/src/CreateSimulationWorkspace.cpp
+++ b/Framework/DataHandling/src/CreateSimulationWorkspace.cpp
@@ -184,7 +184,7 @@ void CreateSimulationWorkspace::createOutputWorkspace() {
 
   m_progress = boost::make_shared<Progress>(this, 0.5, 0.75, nhistograms);
 
-  PARALLEL_FOR1(m_outputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*m_outputWS))
   for (int64_t i = 0; i < static_cast<int64_t>(nhistograms); ++i) {
     m_outputWS->setBinEdges(i, binBoundaries);
     m_outputWS->mutableY(i) = 1.0;
diff --git a/Framework/DataHandling/src/LoadEventNexus.cpp b/Framework/DataHandling/src/LoadEventNexus.cpp
index 9244311039d023f280653b4baed1f61d002b1fff..eb8d8b5b87486c429a91d1537401502391a9f899 100644
--- a/Framework/DataHandling/src/LoadEventNexus.cpp
+++ b/Framework/DataHandling/src/LoadEventNexus.cpp
@@ -1614,7 +1614,7 @@ void LoadEventNexus::loadEvents(API::Progress *const prog,
       if (mT0 != 0.0) {
         int64_t numHistograms =
             static_cast<int64_t>(m_ws->getNumberHistograms());
-        PARALLEL_FOR1(m_ws)
+        PARALLEL_FOR_IF(Kernel::threadSafe(*m_ws))
         for (int64_t i = 0; i < numHistograms; ++i) {
           PARALLEL_START_INTERUPT_REGION
           // Do the offsetting
diff --git a/Framework/DataHandling/src/LoadTOFRawNexus.cpp b/Framework/DataHandling/src/LoadTOFRawNexus.cpp
index 52fa1ad85ecd08081e66cf81c7e5099e7c2c95eb..064855fe9d38f8c30bc7c2ee75a5918ac977e1fd 100644
--- a/Framework/DataHandling/src/LoadTOFRawNexus.cpp
+++ b/Framework/DataHandling/src/LoadTOFRawNexus.cpp
@@ -542,15 +542,11 @@ void LoadTOFRawNexus::exec() {
   const auto id_to_wi = WS->getDetectorIDToWorkspaceIndexMap();
 
   // Load each bank sequentially
-  // PARALLEL_FOR1(WS)
   for (const auto &bankName : bankNames) {
-    //    PARALLEL_START_INTERUPT_REGION
     prog->report("Loading bank " + bankName);
     g_log.debug() << "Loading bank " << bankName << '\n';
     loadBank(filename, entry_name, bankName, WS, id_to_wi);
-    //    PARALLEL_END_INTERUPT_REGION
   }
-  //  PARALLEL_CHECK_INTERUPT_REGION
 
   // Set some units
   if (m_xUnits == "Ang")
diff --git a/Framework/DataHandling/src/SaveToSNSHistogramNexus.cpp b/Framework/DataHandling/src/SaveToSNSHistogramNexus.cpp
index 7536831d0449bc94b8dc7b38ef74fdf3f6a25044..67c9dc75f64d03b59e4ebf2609ffffa4499230d1 100644
--- a/Framework/DataHandling/src/SaveToSNSHistogramNexus.cpp
+++ b/Framework/DataHandling/src/SaveToSNSHistogramNexus.cpp
@@ -287,7 +287,7 @@ int SaveToSNSHistogramNexus::WriteOutDataOrErrors(
     Timer tim1;
     int ypixels = static_cast<int>(det->ypixels());
 
-    PARALLEL_FOR1(inputWorkspace)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*inputWorkspace))
     for (int y = 0; y < ypixels; y++) {
       PARALLEL_START_INTERUPT_REGION
       // Get the workspace index for the detector ID at this spot
diff --git a/Framework/DataObjects/src/Workspace2D.cpp b/Framework/DataObjects/src/Workspace2D.cpp
index 8047123d634af7d2701ee40553da588e5e9ab614..03e93c5c9e4781c84ac34d4b677b5d83d9be9162 100644
--- a/Framework/DataObjects/src/Workspace2D.cpp
+++ b/Framework/DataObjects/src/Workspace2D.cpp
@@ -43,7 +43,7 @@ Workspace2D::~Workspace2D() {
 // http://social.msdn.microsoft.com/Forums/en-US/2fe4cfc7-ca5c-4665-8026-42e0ba634214/visual-studio-$
 
 #ifdef _MSC_VER
-  PARALLEL_FOR1(this)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*this))
   for (int64_t i = 0; i < static_cast<int64_t>(m_noVectors); i++) {
 #else
   for (size_t i = 0; i < m_noVectors; ++i) {
diff --git a/Framework/Kernel/inc/MantidKernel/MultiThreaded.h b/Framework/Kernel/inc/MantidKernel/MultiThreaded.h
index 07014cbba03a54ca45dede942876485463bc8b8b..6c6c5334510abce3c167fc8f9047ca6a8720460f 100644
--- a/Framework/Kernel/inc/MantidKernel/MultiThreaded.h
+++ b/Framework/Kernel/inc/MantidKernel/MultiThreaded.h
@@ -17,9 +17,6 @@ namespace Kernel {
 template <typename Arg>
 inline typename std::enable_if<std::is_pointer<Arg>::value, bool>::type
 threadSafe(Arg workspace) {
-  static_assert(
-      std::is_base_of<DataItem, typename std::remove_pointer<Arg>::type>::value,
-      "Parameter must be derived from Mantid::Kernel::DataItem!");
   return !workspace || workspace->threadSafe();
 }
 
@@ -33,9 +30,6 @@ threadSafe(Arg workspace) {
 template <typename Arg, typename... Args>
 inline typename std::enable_if<std::is_pointer<Arg>::value, bool>::type
 threadSafe(Arg workspace, Args &&... others) {
-  static_assert(
-      std::is_base_of<DataItem, typename std::remove_pointer<Arg>::type>::value,
-      "Parameter must be derived from Mantid::Kernel::DataItem!");
   return (!workspace || workspace->threadSafe()) &&
          threadSafe(std::forward<Args>(others)...);
 }
@@ -48,8 +42,6 @@ threadSafe(Arg workspace, Args &&... others) {
 template <typename Arg>
 inline typename std::enable_if<!std::is_pointer<Arg>::value, bool>::type
 threadSafe(const Arg &workspace) {
-  static_assert(std::is_base_of<DataItem, Arg>::value,
-                "Parameter must be derived from Mantid::Kernel::DataItem!");
   return workspace.threadSafe();
 }
 
@@ -63,8 +55,6 @@ threadSafe(const Arg &workspace) {
 template <typename Arg, typename... Args>
 inline typename std::enable_if<!std::is_pointer<Arg>::value, bool>::type
 threadSafe(const Arg &workspace, Args &&... others) {
-  static_assert(std::is_base_of<DataItem, Arg>::value,
-                "Parameter must be derived from Mantid::Kernel::DataItem!");
   return workspace.threadSafe() && threadSafe(std::forward<Args>(others)...);
 }
 
@@ -145,23 +135,6 @@ threadSafe(const Arg &workspace, Args &&... others) {
 #define PARALLEL_FOR_NO_WSP_CHECK_FIRSTPRIVATE2(variable1, variable2)          \
   PRAGMA(omp parallel for firstprivate(variable1, variable2) )
 
-/** Includes code to add OpenMP commands to run the next for loop in parallel.
-*		The workspace is checked to ensure it is suitable for
-*multithreaded access
-*   NULL workspaces are assumed suitable
-*/
-#define PARALLEL_FOR1(workspace1)                                              \
-    PRAGMA(omp parallel for if ( !workspace1 || workspace1->threadSafe() ) )
-
-/** Includes code to add OpenMP commands to run the next for loop in parallel.
-*	 Both workspaces are checked to ensure they suitable for multithreaded
-*access
-*  or equal to NULL which is also safe
-*/
-#define PARALLEL_FOR2(workspace1, workspace2)                                   \
-    PRAGMA(omp parallel for if ( ( !workspace1 || workspace1->threadSafe() ) && \
-    ( !workspace2 || workspace2->threadSafe() ) ))
-
 /** Ensures that the next execution line or block is only executed if
 * there are multple threads execting in this region
 */
@@ -216,8 +189,6 @@ threadSafe(const Arg &workspace, Args &&... others) {
 #define PARALLEL_FOR_NO_WSP_CHECK()
 #define PARALLEL_FOR_NOWS_CHECK_FIRSTPRIVATE(variable)
 #define PARALLEL_FOR_NO_WSP_CHECK_FIRSTPRIVATE2(variable1, variable2)
-#define PARALLEL_FOR1(workspace1)
-#define PARALLEL_FOR2(workspace1, workspace2)
 #define IF_PARALLEL if (false)
 #define IF_NOT_PARALLEL
 #define PARALLEL_CRITICAL(name)
diff --git a/Framework/MDAlgorithms/src/CalculateCoverageDGS.cpp b/Framework/MDAlgorithms/src/CalculateCoverageDGS.cpp
index 1db1765f7956e41ea08f09473b0be3f9fb2de2a0..4af3ee1f476b02a17f792b2e51dac3eaf8a76abe 100644
--- a/Framework/MDAlgorithms/src/CalculateCoverageDGS.cpp
+++ b/Framework/MDAlgorithms/src/CalculateCoverageDGS.cpp
@@ -366,18 +366,18 @@ void CalculateCoverageDGS::exec() {
   Mantid::Geometry::GeneralFrame frame2("Q2", "");
   Mantid::Geometry::GeneralFrame frame3("Q3", "");
   Mantid::Geometry::GeneralFrame frame4("meV", "");
-  MDHistoDimension_sptr out1(
-      new MDHistoDimension("Q1", "Q1", frame1, static_cast<coord_t>(q1min),
-                           static_cast<coord_t>(q1max), q1NumBins));
-  MDHistoDimension_sptr out2(
-      new MDHistoDimension("Q2", "Q2", frame2, static_cast<coord_t>(q2min),
-                           static_cast<coord_t>(q2max), q2NumBins));
-  MDHistoDimension_sptr out3(
-      new MDHistoDimension("Q3", "Q3", frame3, static_cast<coord_t>(q3min),
-                           static_cast<coord_t>(q3max), q3NumBins));
-  MDHistoDimension_sptr out4(new MDHistoDimension(
+  auto out1 = boost::make_shared<MDHistoDimension>(
+      "Q1", "Q1", frame1, static_cast<coord_t>(q1min),
+      static_cast<coord_t>(q1max), q1NumBins);
+  auto out2 = boost::make_shared<MDHistoDimension>(
+      "Q2", "Q2", frame2, static_cast<coord_t>(q2min),
+      static_cast<coord_t>(q2max), q2NumBins);
+  auto out3 = boost::make_shared<MDHistoDimension>(
+      "Q3", "Q3", frame3, static_cast<coord_t>(q3min),
+      static_cast<coord_t>(q3max), q3NumBins);
+  auto out4 = boost::make_shared<MDHistoDimension>(
       "DeltaE", "DeltaE", frame4, static_cast<coord_t>(m_dEmin),
-      static_cast<coord_t>(m_dEmax), dENumBins));
+      static_cast<coord_t>(m_dEmax), dENumBins);
 
   for (size_t row = 0; row <= 3; row++) {
     if (affineMat[row][0] == 1.) {
@@ -394,7 +394,7 @@ void CalculateCoverageDGS::exec() {
     }
   }
 
-  m_normWS = MDHistoWorkspace_sptr(new MDHistoWorkspace(binDimensions));
+  m_normWS = boost::make_shared<MDHistoWorkspace>(binDimensions);
   m_normWS->setTo(0., 0., 0.);
   setProperty("OutputWorkspace",
               boost::dynamic_pointer_cast<Workspace>(m_normWS));
@@ -403,7 +403,7 @@ void CalculateCoverageDGS::exec() {
 
   const int64_t ndets = static_cast<int64_t>(tt.size());
 
-  PARALLEL_FOR1(inputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS))
   for (int64_t i = 0; i < ndets; i++) {
     PARALLEL_START_INTERUPT_REGION
     auto intersections = calculateIntersections(tt[i], phi[i]);
diff --git a/Framework/MDAlgorithms/src/ConvertToDiffractionMDWorkspace.cpp b/Framework/MDAlgorithms/src/ConvertToDiffractionMDWorkspace.cpp
index b28b4c769dec1b03b78a47bbcf633d3f052ba111..67072208957d5369f20b051c98230b8418fd9496 100644
--- a/Framework/MDAlgorithms/src/ConvertToDiffractionMDWorkspace.cpp
+++ b/Framework/MDAlgorithms/src/ConvertToDiffractionMDWorkspace.cpp
@@ -535,7 +535,7 @@ void ConvertToDiffractionMDWorkspace::exec() {
     }
 
     // 2. Process next chunk of spectra (threaded)
-    PARALLEL_FOR1(m_inWS)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*m_inWS))
     for (int i = start; i < static_cast<int>(wi); ++i) {
       PARALLEL_START_INTERUPT_REGION
       this->convertSpectrum(static_cast<int>(i));
diff --git a/Framework/MDAlgorithms/src/IntegrateEllipsoids.cpp b/Framework/MDAlgorithms/src/IntegrateEllipsoids.cpp
index d282d058f97c9b80fd30c6bfafa23eac566a1632..d66b667fd9674db809b6ba27aaf2614bd53db80b 100644
--- a/Framework/MDAlgorithms/src/IntegrateEllipsoids.cpp
+++ b/Framework/MDAlgorithms/src/IntegrateEllipsoids.cpp
@@ -55,7 +55,7 @@ void IntegrateEllipsoids::qListFromEventWS(Integrate3DEvents &integrator,
   // loop through the eventlists
 
   int numSpectra = static_cast<int>(wksp->getNumberHistograms());
-  PARALLEL_FOR1(wksp)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*wksp))
   for (int i = 0; i < numSpectra; ++i) {
     PARALLEL_START_INTERUPT_REGION
 
@@ -129,7 +129,7 @@ void IntegrateEllipsoids::qListFromHistoWS(Integrate3DEvents &integrator,
 
   int numSpectra = static_cast<int>(wksp->getNumberHistograms());
   const bool histogramForm = wksp->isHistogramData();
-  PARALLEL_FOR1(wksp)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*wksp))
   for (int i = 0; i < numSpectra; ++i) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/MDAlgorithms/src/IntegratePeaksMDHKL.cpp b/Framework/MDAlgorithms/src/IntegratePeaksMDHKL.cpp
index 2a0f1fb1deb2d610abd4bb4a39bb58b3400cd5d3..0e34fc2c51c19ec8a2085c0ec9de16e8e63fb2b4 100644
--- a/Framework/MDAlgorithms/src/IntegratePeaksMDHKL.cpp
+++ b/Framework/MDAlgorithms/src/IntegratePeaksMDHKL.cpp
@@ -106,7 +106,7 @@ void IntegratePeaksMDHKL::exec() {
   int npeaks = peakWS->getNumberPeaks();
 
   auto prog = make_unique<Progress>(this, 0.3, 1.0, npeaks);
-  PARALLEL_FOR1(peakWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*peakWS))
   for (int i = 0; i < npeaks; i++) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/MDAlgorithms/src/MDNormSCD.cpp b/Framework/MDAlgorithms/src/MDNormSCD.cpp
index 1370ab0c8211e48e7c08412df0f24d1300aa8c21..99efceaf9c60422b9efa5850798b430cc6af508c 100644
--- a/Framework/MDAlgorithms/src/MDNormSCD.cpp
+++ b/Framework/MDAlgorithms/src/MDNormSCD.cpp
@@ -401,7 +401,7 @@ void MDNormSCD::calculateNormalization(
       solidAngleWS->getDetectorIDToWorkspaceIndexMap();
 
   auto prog = make_unique<API::Progress>(this, 0.3, 1.0, ndets);
-  PARALLEL_FOR1(integrFlux)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*integrFlux))
   for (int64_t i = 0; i < ndets; i++) {
     PARALLEL_START_INTERUPT_REGION
 
diff --git a/Framework/MDAlgorithms/src/ThresholdMD.cpp b/Framework/MDAlgorithms/src/ThresholdMD.cpp
index 392fdf687502c0718f634695abee6449546201e3..d760285d3e022ac0f7d60fd6afff886cf2c7b181 100644
--- a/Framework/MDAlgorithms/src/ThresholdMD.cpp
+++ b/Framework/MDAlgorithms/src/ThresholdMD.cpp
@@ -109,7 +109,7 @@ void ThresholdMD::exec() {
     frequency = nPoints / 100;
   }
 
-  PARALLEL_FOR2(inputWS, outWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outWS))
   for (int64_t i = 0; i < nPoints; ++i) {
     PARALLEL_START_INTERUPT_REGION
     const double signalAt = inputWS->getSignalAt(i);
diff --git a/Framework/WorkflowAlgorithms/src/SANSSolidAngleCorrection.cpp b/Framework/WorkflowAlgorithms/src/SANSSolidAngleCorrection.cpp
index 41accef33b20b2692474fdffb1e745e07bcf7f9f..5d1452e6baba2bb7c09a6e7d66da2c8bb9c10d37 100644
--- a/Framework/WorkflowAlgorithms/src/SANSSolidAngleCorrection.cpp
+++ b/Framework/WorkflowAlgorithms/src/SANSSolidAngleCorrection.cpp
@@ -110,7 +110,7 @@ void SANSSolidAngleCorrection::exec() {
   // Number of X bins
   const int xLength = static_cast<int>(inputWS->readY(0).size());
 
-  PARALLEL_FOR2(outputWS, inputWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS, *inputWS))
   for (int i = 0; i < numHists; ++i) {
     PARALLEL_START_INTERUPT_REGION
     outputWS->dataX(i) = inputWS->readX(i);
@@ -188,7 +188,7 @@ void SANSSolidAngleCorrection::execEvent() {
   Progress progress(this, 0.0, 1.0, numberOfSpectra);
   progress.report("Solid Angle Correction");
 
-  PARALLEL_FOR1(outputEventWS)
+  PARALLEL_FOR_IF(Kernel::threadSafe(*outputEventWS))
   for (int i = 0; i < numberOfSpectra; i++) {
     PARALLEL_START_INTERUPT_REGION
     IDetector_const_sptr det;
diff --git a/MantidPlot/src/Mantid/MantidMatrix.cpp b/MantidPlot/src/Mantid/MantidMatrix.cpp
index 886a0030b2f9c13b48cce97f933d4debde37bb3e..6251323398288f29ec6f4c3403021900200727fe 100644
--- a/MantidPlot/src/Mantid/MantidMatrix.cpp
+++ b/MantidPlot/src/Mantid/MantidMatrix.cpp
@@ -1136,7 +1136,7 @@ void findYRange(MatrixWorkspace_const_sptr ws, double &miny, double &maxy) {
 
   if (ws) {
 
-    PARALLEL_FOR1(ws)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*ws))
     for (int wi = 0; wi < static_cast<int>(ws->getNumberHistograms()); wi++) {
       double local_min, local_max;
       const Mantid::MantidVec &Y = ws->readY(wi);
diff --git a/MantidPlot/src/Mantid/MantidUI.cpp b/MantidPlot/src/Mantid/MantidUI.cpp
index 9d94ea9d412d8db4741b59e96d07497cbca226f8..1b17315abf88f69c71de45cb2b4432584b6355b2 100644
--- a/MantidPlot/src/Mantid/MantidUI.cpp
+++ b/MantidPlot/src/Mantid/MantidUI.cpp
@@ -1270,7 +1270,7 @@ Table *MantidUI::createDetectorTable(
                                  // value should be displayed
   QVector<QList<QVariant>> tableColValues;
   tableColValues.resize(nrows);
-  PARALLEL_FOR1(ws)
+  PARALLEL_FOR_IF(Mantid::Kernel::threadSafe(*ws))
   for (int row = 0; row < nrows; ++row) {
     // Note PARALLEL_START_INTERUPT_REGION & friends apparently not needed (like
     // in algorithms)
diff --git a/MantidQt/MantidWidgets/src/InstrumentView/InstrumentActor.cpp b/MantidQt/MantidWidgets/src/InstrumentView/InstrumentActor.cpp
index 5341ec549a17d401ffb4a23217cede43ac3b2608..7824ad8d701be064879232d01bc4478d04a020c1 100644
--- a/MantidQt/MantidWidgets/src/InstrumentView/InstrumentActor.cpp
+++ b/MantidQt/MantidWidgets/src/InstrumentView/InstrumentActor.cpp
@@ -638,7 +638,6 @@ void InstrumentActor::resetColors() {
   Instrument_const_sptr inst = getInstrument();
   IMaskWorkspace_sptr mask = getMaskWorkspaceIfExists();
 
-  // PARALLEL_FOR1(m_workspace)
   for (int iwi = 0; iwi < int(m_specIntegrs.size()); iwi++) {
     size_t wi = size_t(iwi);
     double integratedValue = m_specIntegrs[wi];
diff --git a/Vates/VatesAPI/src/LoadVTK.cpp b/Vates/VatesAPI/src/LoadVTK.cpp
index 36fcc0ca5ae8a7d3a9683adb993a8554ed18c2fd..25e2570d3a38b6a727bfece221fbe10a557bd181 100644
--- a/Vates/VatesAPI/src/LoadVTK.cpp
+++ b/Vates/VatesAPI/src/LoadVTK.cpp
@@ -258,7 +258,7 @@ void LoadVTK::execMDEvent(vtkDataSet *readDataset,
   ws->initialize();
 
   if (errorsSQ == NULL) {
-    PARALLEL_FOR1(ws)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*ws))
     for (int64_t i = 0; i < nPoints; ++i) {
       PARALLEL_START_INTERUPT_REGION
       double coordinates[3];
@@ -275,7 +275,7 @@ void LoadVTK::execMDEvent(vtkDataSet *readDataset,
     }
     PARALLEL_CHECK_INTERUPT_REGION
   } else {
-    PARALLEL_FOR1(ws)
+    PARALLEL_FOR_IF(Kernel::threadSafe(*ws))
     for (int64_t i = 0; i < nPoints; ++i) {
       PARALLEL_START_INTERUPT_REGION
       double coordinates[3];