diff --git a/Base/Util/StringUtil.cpp b/Base/Util/StringUtil.cpp
index 36a48437c9a26255d103e512c5a68e5e4a2fae5e..01e741d7c8a43d071d84007e5a34baa42a88e5b3 100644
--- a/Base/Util/StringUtil.cpp
+++ b/Base/Util/StringUtil.cpp
@@ -17,7 +17,6 @@
 #include <boost/algorithm/string.hpp>
 #include <cerrno>
 #include <charconv>
-#include <regex>
 #include <string>
 
 //! Returns string right-padded with blanks.
diff --git a/Device/IO/IOFactory.cpp b/Device/IO/IOFactory.cpp
index 071c19de807904d6b35e2eaa322f6ec1f4f559e2..12672aa5ec3b4a57add3294ab53268a5308b8ef6 100644
--- a/Device/IO/IOFactory.cpp
+++ b/Device/IO/IOFactory.cpp
@@ -265,9 +265,10 @@ bool Util::RW::dataMatchesFile(const Datafield& data, const std::string& refFile
     std::unique_ptr<Datafield> refDat;
     try {
         refDat.reset(IO::readData2D(refFileName));
-    } catch (...) {
+    } catch (const std::exception& ex) {
         std::cerr << "File comparison: Could not read reference data from file " << refFileName
                   << std::endl;
+        std::cerr << "Error message: " << ex.what() << std::endl;
         return false;
     }
     ASSERT(refDat);
diff --git a/Device/IO/ParseUtil.cpp b/Device/IO/ParseUtil.cpp
index 308df0d6a7938b6c82c1ea92c1b70757532600a6..e932c5579a072967d5b0e1ed11b6fe98f8bb7115 100644
--- a/Device/IO/ParseUtil.cpp
+++ b/Device/IO/ParseUtil.cpp
@@ -21,18 +21,10 @@
 #include <functional>
 #include <iostream>
 #include <iterator>
+#include <regex>
 
 namespace {
 
-std::istringstream getAxisStringRepresentation(std::istream& input_stream)
-{
-    std::string line;
-    std::getline(input_stream, line);
-    const std::vector<std::string> to_replace = {",", "\"", "(", ")", "[", "]"};
-    Base::String::replaceItemsFromString(line, to_replace, " ");
-    return std::istringstream(line);
-}
-
 void readLineOfDoubles(std::vector<double>& buffer, std::istringstream& iss)
 {
     iss.imbue(std::locale::classic());
@@ -40,42 +32,59 @@ void readLineOfDoubles(std::vector<double>& buffer, std::istringstream& iss)
               back_inserter(buffer));
 }
 
+std::vector<double> parse_x_list(std::string text, const std::string& type)
+{
+    if (text.substr(0, 1) == ",") {
+        std::cerr << "Warning: file format from BornAgain <= 20 is obsolete" << std::endl;
+        size_t i = text.find_first_not_of(" ", 1);
+        text.erase(0, i);
+    }
+    if (text.substr(0, 1) == "[" && text.substr(text.size() - 1, 1) == "]")
+        text = text.substr(1, text.size() - 2);
+    std::vector<std::string> arr = Base::String::split(text, ",");
+    std::vector<double> result;
+    for (std::string e : arr) {
+        e = Base::String::trim(e);
+        if (e.empty())
+            continue;
+        double x;
+        if (!(Base::String::to_double(e, &x)))
+            throw std::runtime_error("Reading " + type + ": cannot read entry '" + e + "'");
+        result.push_back(x);
+    }
+    if (!result.size())
+        throw std::runtime_error("Reading " + type + ": found empty list");
+    return result;
+}
+
 } // namespace
 
 //! Creates axis of certain type from input stream
 Scale* Util::Parse::parseScale(std::istream& input_stream)
 {
-    auto iss = ::getAxisStringRepresentation(input_stream);
-    std::string type;
-    if (!(iss >> type))
-        throw std::runtime_error("Cannot read axis type from input");
-
+    std::string line;
+    std::getline(input_stream, line);
+    std::smatch matched;
+    std::regex_match(line, matched, std::regex("(\\w+)\\(\"(.*?)\",\\s*(.+?)\\)"));
+    if (matched.size() != 4)
+        throw std::runtime_error("Cannot read axis from input");
+    const std::string type = matched[1];
+    const std::string name = matched[2];
+    const std::string body = matched[3];
     if (type == "EquiDivision" || type == "FixedBinAxis" /* for compatibility with pre-21 */) {
-        std::string name;
-        size_t nbins;
-        if (!(iss >> name >> nbins))
-            throw std::runtime_error("Reading EquiDivision: cannot read name or size");
-        std::vector<double> boundaries;
-        ::readLineOfDoubles(boundaries, iss);
-        if (boundaries.size() != 2)
-            throw std::runtime_error("Reading EquiDivision: cannot read start or stop");
-        return newEquiDivision(name, nbins, boundaries[0], boundaries[1]);
+        std::vector<std::string> arr = Base::String::split(body, ",");
+        int nbins;
+        double xmi, xma;
+        if (!(Base::String::to_int(arr[0], &nbins) && Base::String::to_double(arr[1], &xmi)
+              && Base::String::to_double(arr[2], &xma)))
+            throw std::runtime_error("Reading EquiDivision: cannot read parameters");
+        return newEquiDivision(name, nbins, xmi, xma);
 
     } else if (type == "ListScan" || type == "DiscreteAxis" /* for compatibility with pre-21 */) {
-        std::string name;
-        if (!(iss >> name))
-            throw std::runtime_error("Reading ListScan: cannot read name");
-        std::vector<double> coordinates;
-        ::readLineOfDoubles(coordinates, iss);
-        return newListScan(name, coordinates);
+        return newListScan(name, parse_x_list(body, type));
 
     } else if (type == "Scale") {
-        std::string name;
-        if (!(iss >> name))
-            throw std::runtime_error("Reading Scale: cannot read name");
-        std::vector<double> coordinates;
-        ::readLineOfDoubles(coordinates, iss);
-        return newGenericScale(name, coordinates);
+        return newGenericScale(name, parse_x_list(body, type));
 
     } else
         throw std::runtime_error("Unknown axis type '" + type + "'");
diff --git a/Device/IO/ReadWriteINT.cpp b/Device/IO/ReadWriteINT.cpp
index 5e58ad25cd9bf2b1b7dc4751167ee686752dcadc..06f4fff3b0f192e99a7b67033393abe70b7db802 100644
--- a/Device/IO/ReadWriteINT.cpp
+++ b/Device/IO/ReadWriteINT.cpp
@@ -60,12 +60,10 @@ void Util::RW::writeBAInt(const Datafield& data, std::ostream& output_stream)
     output_stream << "# BornAgain Intensity Data\n\n";
 
     for (size_t i = 0; i < data.rank(); ++i) {
-        std::string axis_name = std::string("axis") + std::to_string(i);
-        std::unique_ptr<Scale> axis(data.axis(i).clone());
-        axis->setAxisName(axis_name);
+        const Scale& axis = data.axis(i);
         output_stream << std::endl;
         output_stream << "# axis-" << i << "\n";
-        output_stream << (*axis) << "\n";
+        output_stream << axis << "\n";
     }
     size_t n_columns = data.axis(data.rank() - 1).size();
 
diff --git a/Tests/ReferenceData/Suite/PolarizedFeNiBilayerTanhPP.int.gz b/Tests/ReferenceData/Suite/PolarizedFeNiBilayerTanhPP.int.gz
index baff0c472397b9c5371032c6fcde4cfdc9a24d45..d5ce2a7307588a3deae1ab89962fa0c017964e1f 100644
Binary files a/Tests/ReferenceData/Suite/PolarizedFeNiBilayerTanhPP.int.gz and b/Tests/ReferenceData/Suite/PolarizedFeNiBilayerTanhPP.int.gz differ
diff --git a/Tests/ReferenceData/Suite/PolarizedFeNiBilayerTanhPP_Q.int.gz b/Tests/ReferenceData/Suite/PolarizedFeNiBilayerTanhPP_Q.int.gz
index 1902f8abd37dea21791141abafbf45f6dc7b50af..2da5f61237e7c8f5b3c8eed1b819886970de4083 100644
Binary files a/Tests/ReferenceData/Suite/PolarizedFeNiBilayerTanhPP_Q.int.gz and b/Tests/ReferenceData/Suite/PolarizedFeNiBilayerTanhPP_Q.int.gz differ