diff --git a/Core/Simulation/ISimulation.cpp b/Core/Simulation/ISimulation.cpp
index a60df0cc46c2bb5963719c4119bc50441522b164..0cc6a25d0e916434acf6b8cfea8a3198227f3abf 100644
--- a/Core/Simulation/ISimulation.cpp
+++ b/Core/Simulation/ISimulation.cpp
@@ -19,6 +19,7 @@
 #include "Core/Computation/IComputation.h"
 #include "Core/Simulation/MPISimulation.h"
 #include "Device/Beam/Beam.h"
+#include "Resample/Options/SimulationOptions.h"
 #include "Resample/Processed/ProcessedSample.h"
 #include "Sample/Multilayer/MultilayerUtils.h"
 #include "Sample/SampleBuilderEngine/ISampleBuilder.h"
@@ -105,20 +106,24 @@ void runComputations(std::vector<std::unique_ptr<IComputation>>& computations)
 //  ************************************************************************************************
 
 ISimulation::ISimulation(const Beam& beam, const MultiLayer& sample, const IDetector& detector)
-    : m_instrument(beam, detector)
+    : m_options(std::make_unique<SimulationOptions>())
+    , m_instrument(beam, detector)
 {
     setSample(sample);
     initialize();
 }
 
 #ifndef SWIG
-ISimulation::ISimulation(const Beam& beam, const IDetector& detector) : m_instrument(beam, detector)
+ISimulation::ISimulation(const Beam& beam, const IDetector& detector)
+    : m_options(std::make_unique<SimulationOptions>())
+    , m_instrument(beam, detector)
 {
     initialize();
 }
 #endif // SWIG
 
 ISimulation::ISimulation()
+    : m_options(std::make_unique<SimulationOptions>())
 {
     initialize();
 }
@@ -131,6 +136,19 @@ void ISimulation::initialize()
     registerChild(&m_sample_provider);
 }
 
+void ISimulation::setOptions(const SimulationOptions& options) {
+    ASSERT(m_options);
+    *m_options = options;
+}
+
+const SimulationOptions& ISimulation::options() const {
+    ASSERT(m_options);
+    return *m_options; }
+
+SimulationOptions& ISimulation::options() {
+    ASSERT(m_options);
+    return *m_options; }
+
 //! Initializes a progress monitor that prints to stdout.
 void ISimulation::setTerminalProgressMonitor()
 {
@@ -170,8 +188,8 @@ void ISimulation::runSimulation()
     m_progress.setExpectedNTicks(param_combinations * total_size);
 
     // restrict calculation to current batch
-    const size_t n_batches = m_options.getNumberOfBatches();
-    const size_t current_batch = m_options.getCurrentBatch();
+    const size_t n_batches = m_options->getNumberOfBatches();
+    const size_t current_batch = m_options->getCurrentBatch();
 
     const size_t batch_start = startIndex(n_batches, current_batch, total_size);
     const size_t batch_size = batchSize(n_batches, current_batch, total_size);
@@ -293,7 +311,7 @@ void ISimulation::runSingleSimulation(const ProcessedSample& re_sample, size_t b
 {
     initElementVector();
 
-    const size_t n_threads = m_options.getNumberOfThreads();
+    const size_t n_threads = m_options->getNumberOfThreads();
     ASSERT(n_threads > 0);
 
     std::vector<std::unique_ptr<IComputation>> computations;
diff --git a/Core/Simulation/ISimulation.h b/Core/Simulation/ISimulation.h
index 399201f364c94eb8d4844c2bf6d8e8b232b6080d..fbb1b886d13594f347e77441ec2675f15546a7b6 100644
--- a/Core/Simulation/ISimulation.h
+++ b/Core/Simulation/ISimulation.h
@@ -20,7 +20,6 @@
 #include "Device/Detector/IDetector2D.h"
 #include "Device/Instrument/Instrument.h"
 #include "Param/Distrib/DistributionHandler.h"
-#include "Resample/Options/SimulationOptions.h"
 #include "Sample/SampleBuilderEngine/SampleProvider.h"
 
 template <class T> class OutputData;
@@ -30,6 +29,7 @@ class ICoordSystem;
 class ISampleBuilder;
 class MultiLayer;
 class ProcessedSample;
+class SimulationOptions;
 class SimulationResult;
 
 //! Abstract base class of OffSpecularSimulation, GISASSimulation and SpecularSimulation.
@@ -80,9 +80,9 @@ public:
     void addParameterDistribution(const ParameterDistribution& par_distr);
     const std::vector<ParameterDistribution>& getDistributions() const;
 
-    void setOptions(const SimulationOptions& options) { m_options = options; }
-    const SimulationOptions& options() const { return m_options; }
-    SimulationOptions& options() { return m_options; }
+    void setOptions(const SimulationOptions& options);
+    const SimulationOptions& options() const;
+    SimulationOptions& options();
 
     void subscribe(ProgressHandler::Callback_t inform) { m_progress.subscribe(inform); }
     void setTerminalProgressMonitor();
@@ -147,7 +147,7 @@ private:
     virtual std::vector<double> rawResults() const = 0;
     virtual void setRawResults(const std::vector<double>& raw_data) = 0;
 
-    SimulationOptions m_options;
+    std::unique_ptr<SimulationOptions> m_options;
     ProgressHandler m_progress;
     SampleProvider m_sample_provider;
     DistributionHandler m_distribution_handler;
diff --git a/auto/Wrap/doxygenCore.i b/auto/Wrap/doxygenCore.i
index 66a3dffd792180c7d71181ba8758763da7582f09..a34f704bbeb0f2735625181eb192d84b87031c5c 100644
--- a/auto/Wrap/doxygenCore.i
+++ b/auto/Wrap/doxygenCore.i
@@ -964,10 +964,10 @@ Returns the results of the simulation in a format that supports unit conversion
 %feature("docstring")  ISimulation::setOptions "void ISimulation::setOptions(const SimulationOptions &options)
 ";
 
-%feature("docstring")  ISimulation::options "const SimulationOptions& ISimulation::options() const
+%feature("docstring")  ISimulation::options "const SimulationOptions & ISimulation::options() const
 ";
 
-%feature("docstring")  ISimulation::options "SimulationOptions& ISimulation::options()
+%feature("docstring")  ISimulation::options "SimulationOptions & ISimulation::options()
 ";
 
 %feature("docstring")  ISimulation::subscribe "void ISimulation::subscribe(ProgressHandler::Callback_t inform)
@@ -2527,6 +2527,8 @@ GISAS simulation with an extra long wavelength.
 
 
 // File: ISimulation_8cpp.xml
+%feature("docstring")  ASSERT "ASSERT(m_options)
+";
 
 
 // File: ISimulation_8h.xml
diff --git a/auto/Wrap/libBornAgainCore.py b/auto/Wrap/libBornAgainCore.py
index a1bef13a9d7a9260e909383d08470b56536e9cc3..2d56d591bf8d0cefa21e933bd2747b4d92b15c0a 100644
--- a/auto/Wrap/libBornAgainCore.py
+++ b/auto/Wrap/libBornAgainCore.py
@@ -3591,7 +3591,7 @@ class ISimulation(libBornAgainParam.INode):
         r"""
         options(ISimulation self) -> SimulationOptions const
         options(ISimulation self) -> SimulationOptions &
-        SimulationOptions& ISimulation::options()
+        SimulationOptions & ISimulation::options()
 
         """
         return _libBornAgainCore.ISimulation_options(self, *args)
diff --git a/auto/Wrap/libBornAgainCore_wrap.cpp b/auto/Wrap/libBornAgainCore_wrap.cpp
index 2ab4a84528098a15aa2758a7369c163304a0fe01..3d36675495ac94848aa9be7d1821e7ecc85ac305 100644
--- a/auto/Wrap/libBornAgainCore_wrap.cpp
+++ b/auto/Wrap/libBornAgainCore_wrap.cpp
@@ -44282,7 +44282,7 @@ static PyMethodDef SwigMethods[] = {
 	 { "ISimulation_options", _wrap_ISimulation_options, METH_VARARGS, "\n"
 		"ISimulation_options(ISimulation self) -> SimulationOptions const\n"
 		"ISimulation_options(ISimulation self) -> SimulationOptions &\n"
-		"SimulationOptions& ISimulation::options()\n"
+		"SimulationOptions & ISimulation::options()\n"
 		"\n"
 		""},
 	 { "ISimulation_subscribe", _wrap_ISimulation_subscribe, METH_VARARGS, "\n"