Create dynamics interface to combine motion and dist

README.md

@@ -1,32 +1,37 @@
-# Download
+# NMR Random Walk Simulation
 
-1. Clone repository with `git clone https://gitea.pkm.physik.tu-darmstadt.de/NMRRWSims/cpp.git` 
-2. After download, change permissions of build.sh in terminal with `chmod a+x build.sh`
+Simulate solid-state NMR spectra and stimulated echo (STE) decays using random walk models.
 
 # Build
 
-# Running
-
-## Command line
-
-To run a random walk simulation Execute `rwsim` in the command line with 
+Requires CMake 3.18+, a C++17 compiler, and OpenMP.
 
 ```bash
-  ./rwsim /path/to/config.txt MotionModel DistributionModel --ste --spectrum -ARG1 1 -ARG2 2
+cmake -B build
+cmake --build build
 ```
 
-It needs three positional arguments: the path to your basic configuration file (see below) and the names of the motional model and of the distribution of correlation times.
-Set the optional arguments `--ste` and/or `--spectrum` to create Stimulated Echos or normal echo spectra, respectively.
-You can overwrite any parameter given in the configuration file by adding it as optional argument with a numerical value, e.g. `-TAU 1e-3` for a correlation time of 1 ms.
+The executable is `build/src/rwsim`.
 
+# Running
 
-## Configuration file `config.txt`
+```bash
+./build/src/rwsim /path/to/config.txt MotionModel DistributionModel --ste --spectrum -ARG1 1 -ARG2 2
+```
 
-Change the values of delta, eta, mixing times, echo times,... in this file. If a paramter is defined multiple times, only the last one is used.  
+Three positional arguments are required: the path to a configuration file, the motion model name, and the distribution model name.
+Set `--ste` and/or `--spectrum` to run stimulated echo or spectrum simulations.
+Any parameter from the configuration file can be overridden on the command line, e.g. `-tau 1e-3`.
 
-### General parameter 
+Use `-seed 42` for reproducible results. Without a seed, a random one is used.
 
-This list of general parameter are necessary for all simulations:
+**Output files are overwritten if a simulation with the same parameters is run again.**
+
+## Configuration file
+
+Change delta, eta, mixing times, echo times, etc. in this file. If a parameter appears multiple times, only the last value is used.
+
+### General parameters
 
 | Parameter  | Description                |
 |------------|----------------------------|
@@ -36,48 +41,88 @@ This list of general parameter are necessary for all simulations:
 
 ### Distribution of correlation times
 
-Two distributions are available: A delta distribution `Delta`, i.e. the same correlation time for every walker, and a log-normal distribution `LogNormal`.
-
-+ Parameters for delta distribution `Delta`
-
-  | Parameter | Description           |
-  |-----------|-----------------------|
-  | tau       | Jump time in s        |
-
-+ Parameters for log-normal distribution `LogNormal`
-
-  | Parameter | Description                                |
-  |-----------|--------------------------------------------|
-  | tau       | Maximum jump time of the distribution in s |
-  | sigma     | Standard deviation of the distribution     |
+| Model | CLI name | Parameters |
+|-------|----------|------------|
+| Delta distribution (same tau for all walkers) | `Delta` | `tau` (jump time in s) |
+| Log-normal distribution | `LogNormal` | `tau` (peak time in s), `sigma` (width) |
 
 ### Motion models
 
-Four different jump models are implemented: Isotropic random jump `RandomJump`, isotropic jump with a certain angle, isotropic jump with a bimodal distribution of angles, and a tetrahedral jump `TetrahedralJump`.
+| Model | CLI name | Parameters |
+|-------|----------|------------|
+| Isotropic random jump | `RandomJump` | — |
+| Isotropic jump by fixed angle | `IsotropicAngle` | `angle` (degrees) |
+| Bimodal angle distribution | `BimodalAngle` | `angle1`, `angle2` (degrees), `probability1` (0–1) |
+| Tetrahedral 4-site jump | `FourSiteTetrahedral` | — |
+| Octahedral 6-site jump (C3) | `SixSiteOctahedralC3` | — |
+| N-site jump on cone | `NSiteConeJump` | `num_sites`, `chi` (cone angle, degrees) |
+| Random jump on cone | `RandomJumpOnCone` | `angle` (cone angle) |
+| Wobbling in cone | `ConeWobble` | `angle` (cone angle) |
 
-+ Isotropic random jump `RandomJump` does not have additional parameters. 
-+ Tetrahedral jump `TetrahedralJump` does not have additional parameters.
-+ Parameters for isotropic jump by an angle `IsotropicAngle`
+### Spectrum parameters
 
-  | Parameter | Description          |
-  |-----------|----------------------|
-  | angle     | Jump angle in degree |
+| Parameter    | Description                     |
+|--------------|---------------------------------|
+| dwell_time   | Acquisition dwell time in s     |
+| num_acq      | Number of acquisition points    |
+| techo_start  | First echo time in s            |
+| techo_stop   | Last echo time in s             |
+| techo_steps  | Number of echo times            |
 
-+ Parameters for isotropic jump with bimodal angle distribution `BimodalAngle`
+### STE parameters
 
-  | Parameter    | Description                          |
-  |--------------|--------------------------------------|
-  | angle1       | First jump angle in degree           |
-  | angle2       | Second jump angle in degree          |
-  | probability1 | Probality that jump has angle1 (0-1) |  
+| Parameter   | Description                     |
+|-------------|---------------------------------|
+| tevo_start  | First evolution time in s       |
+| tevo_stop   | Last evolution time in s        |
+| tevo_steps  | Number of evolution times       |
+| tmix_start  | First mixing time in s          |
+| tmix_stop   | Last mixing time in s           |
+| tmix_steps  | Number of mixing times (log-spaced) |
+| tpulse4     | Delay of 4th pulse in s         |
 
+# Architecture
 
+The simulation is built around three abstractions:
 
+## Dynamics
 
+`Dynamics` (`src/dynamics/base.h`) produces trajectory steps — pairs of `{dt, omega}` representing a waiting time and the NMR frequency after a jump. This is the core interface for defining physical models.
 
+**`DecoupledDynamics`** wraps a separate motion model and time distribution, calling them independently. This covers all existing models.
 
-# Running
+To implement a **coupled model** where motion and time are interdependent (e.g. jump angle affects waiting time, or position-dependent dynamics), implement the `Dynamics` interface directly:
 
-**Because filenames are always the same, previous simulations and results are overwritten!**
+```cpp
+class MyModel : public Dynamics {
+    void initialize(std::mt19937_64& rng) override { /* set initial state */ }
+    Step next(std::mt19937_64& rng) override {
+        // draw angle, compute rate from angle, draw wait time from rate
+        return {dt, omega};
+    }
+    // ... clone(), setParameters(), etc.
+};
+```
 
-1. Execute `build.sh` (in terminal with `./build.sh`). It compiles the source code and executes test.py
+## Experiment
+
+`Experiment` (`src/experiments/base.h`) defines how trajectory data is turned into observables:
+
+- `setup(params)` — configure time axes, allocate accumulators
+- `accumulate(trajectory)` — process one walker's trajectory
+- `save(directory)` — write results to files
+- `clone()` / `merge()` — for thread-parallel accumulation
+
+Existing experiments: `SpectrumExperiment`, `STEExperiment`.
+
+## Simulation runner
+
+`run_simulation()` (`src/sims.h`) connects a `Dynamics` and an `Experiment`:
+
+1. Sets parameters on both
+2. Generates trajectories in parallel (OpenMP, one clone per thread)
+3. Each thread accumulates into its own experiment clone
+4. Merges thread-local results
+5. Saves output
+
+The walker loop is parallelized with OpenMP using static scheduling for deterministic reproducibility with a given seed and thread count.

src/CMakeLists.txt

@@ -3,15 +3,16 @@ add_subdirectory(times)
 add_subdirectory(motions)
 add_subdirectory(utils)
 add_subdirectory(experiments)
+add_subdirectory(dynamics)
 
 add_library(simulation STATIC sims.cpp sims.h)
-target_link_libraries(simulation PRIVATE utils experiments OpenMP::OpenMP_CXX)
+target_link_libraries(simulation PRIVATE utils experiments dynamics OpenMP::OpenMP_CXX)
 
 add_executable(
         rwsim
         main.cpp
 )
 
-target_link_libraries(rwsim PUBLIC RWMotion RWTime utils experiments simulation)
+target_link_libraries(rwsim PUBLIC RWMotion RWTime utils experiments dynamics simulation)
 
 target_compile_options(rwsim PUBLIC -Werror -Wall -Wextra -Wconversion -O2)

src/dynamics/CMakeLists.txt

@@ -0,0 +1,6 @@
+add_library(
+        dynamics STATIC
+        base.h
+        decoupled.cpp
+        decoupled.h
+)

src/dynamics/base.h

@@ -0,0 +1,26 @@
+#ifndef RWSIM_DYNAMICSBASE_H
+#define RWSIM_DYNAMICSBASE_H
+
+#include <memory>
+#include <random>
+#include <string>
+#include <unordered_map>
+
+struct Step {
+    double dt;
+    double omega;
+};
+
+class Dynamics {
+public:
+    virtual ~Dynamics() = default;
+
+    virtual void setParameters(const std::unordered_map<std::string, double>& parameters) = 0;
+    virtual void initialize(std::mt19937_64& rng) = 0;
+    virtual Step next(std::mt19937_64& rng) = 0;
+    [[nodiscard]] virtual double getInitOmega() const = 0;
+    [[nodiscard]] virtual std::unique_ptr<Dynamics> clone() const = 0;
+    [[nodiscard]] virtual std::string toString() const = 0;
+};
+
+#endif //RWSIM_DYNAMICSBASE_H

src/dynamics/decoupled.cpp

@@ -0,0 +1,35 @@
+#include "decoupled.h"
+
+DecoupledDynamics::DecoupledDynamics(
+    std::unique_ptr<motions::BaseMotion> motion,
+    std::unique_ptr<times::BaseDistribution> dist)
+    : m_motion(std::move(motion)), m_dist(std::move(dist)) {}
+
+void DecoupledDynamics::setParameters(const std::unordered_map<std::string, double>& parameters) {
+    m_motion->setParameters(parameters);
+    m_dist->setParameters(parameters);
+}
+
+void DecoupledDynamics::initialize(std::mt19937_64& rng) {
+    m_motion->initialize(rng);
+    m_dist->initialize(rng);
+}
+
+Step DecoupledDynamics::next(std::mt19937_64& rng) {
+    double dt = m_dist->tau_wait(rng);
+    double omega = m_motion->jump(rng);
+    return {dt, omega};
+}
+
+double DecoupledDynamics::getInitOmega() const {
+    return m_motion->getInitOmega();
+}
+
+std::unique_ptr<Dynamics> DecoupledDynamics::clone() const {
+    return std::make_unique<DecoupledDynamics>(
+        m_motion->clone(), m_dist->clone());
+}
+
+std::string DecoupledDynamics::toString() const {
+    return m_motion->toString() + "/" + m_dist->toString();
+}

src/dynamics/decoupled.h

@@ -0,0 +1,27 @@
+#ifndef RWSIM_DECOUPLED_H
+#define RWSIM_DECOUPLED_H
+
+#include "base.h"
+#include "../motions/base.h"
+#include "../times/base.h"
+
+#include <memory>
+
+class DecoupledDynamics final : public Dynamics {
+public:
+    DecoupledDynamics(std::unique_ptr<motions::BaseMotion> motion,
+                      std::unique_ptr<times::BaseDistribution> dist);
+
+    void setParameters(const std::unordered_map<std::string, double>& parameters) override;
+    void initialize(std::mt19937_64& rng) override;
+    Step next(std::mt19937_64& rng) override;
+    [[nodiscard]] double getInitOmega() const override;
+    [[nodiscard]] std::unique_ptr<Dynamics> clone() const override;
+    [[nodiscard]] std::string toString() const override;
+
+private:
+    std::unique_ptr<motions::BaseMotion> m_motion;
+    std::unique_ptr<times::BaseDistribution> m_dist;
+};
+
+#endif //RWSIM_DECOUPLED_H

src/experiments/base.h

@@ -1,9 +1,6 @@
 #ifndef RWSIM_EXPERIMENTBASE_H
 #define RWSIM_EXPERIMENTBASE_H
 
-#include "../motions/base.h"
-#include "../times/base.h"
-
 #include <memory>
 #include <unordered_map>
 #include <string>
@@ -23,7 +20,7 @@ public:
                        const std::unordered_map<std::string, double>& optional) = 0;
     [[nodiscard]] virtual double tmax() const = 0;
     virtual void accumulate(const Trajectory& traj, double init_omega, int num_walkers) = 0;
-    virtual void save(const motions::BaseMotion& motion, const times::BaseDistribution& dist) = 0;
+    virtual void save(const std::string& directory) = 0;
     [[nodiscard]] virtual std::unique_ptr<Experiment> clone() const = 0;
     virtual void merge(const Experiment& other) = 0;
 };

src/experiments/spectrum.cpp

@@ -3,6 +3,7 @@
 #include "../utils/ranges.h"
 #include "../utils/io.h"
 
+
 #include <algorithm>
 #include <cmath>
 
@@ -44,10 +45,9 @@ void SpectrumExperiment::accumulate(const Trajectory& traj, double, int num_walk
     }
 }
 
-void SpectrumExperiment::save(const motions::BaseMotion& motion, const times::BaseDistribution& dist) {
-    const auto path = make_directory(motion, dist);
-    save_parameter_to_file(std::string("timesignal"), path, m_parameter, m_optional);
-    save_data_to_file(std::string("timesignal"), path, m_t_fid, m_fid_dict, m_optional);
+void SpectrumExperiment::save(const std::string& directory) {
+    save_parameter_to_file(std::string("timesignal"), directory, m_parameter, m_optional);
+    save_data_to_file(std::string("timesignal"), directory, m_t_fid, m_fid_dict, m_optional);
 }
 
 std::unique_ptr<Experiment> SpectrumExperiment::clone() const {

src/experiments/spectrum.h

@@ -13,7 +13,7 @@ public:
                const std::unordered_map<std::string, double>& optional) override;
     [[nodiscard]] double tmax() const override;
     void accumulate(const Trajectory& traj, double init_omega, int num_walkers) override;
-    void save(const motions::BaseMotion& motion, const times::BaseDistribution& dist) override;
+    void save(const std::string& directory) override;
     [[nodiscard]] std::unique_ptr<Experiment> clone() const override;
     void merge(const Experiment& other) override;
 

src/experiments/ste.cpp

@@ -70,12 +70,11 @@ void STEExperiment::accumulate(const Trajectory& traj, double init_omega, int nu
     }
 }
 
-void STEExperiment::save(const motions::BaseMotion& motion, const times::BaseDistribution& dist) {
-    const auto folders = make_directory(motion, dist);
-    save_parameter_to_file(std::string("ste"), folders, m_parameter, m_optional);
-    save_data_to_file(std::string("coscos"), folders, m_mixing_times, m_cc_dict, m_optional);
-    save_data_to_file(std::string("sinsin"), folders, m_mixing_times, m_ss_dict, m_optional);
-    save_data_to_file(std::string("f2"), folders, m_mixing_times, m_f2, m_optional);
+void STEExperiment::save(const std::string& directory) {
+    save_parameter_to_file(std::string("ste"), directory, m_parameter, m_optional);
+    save_data_to_file(std::string("coscos"), directory, m_mixing_times, m_cc_dict, m_optional);
+    save_data_to_file(std::string("sinsin"), directory, m_mixing_times, m_ss_dict, m_optional);
+    save_data_to_file(std::string("f2"), directory, m_mixing_times, m_f2, m_optional);
 }
 
 std::unique_ptr<Experiment> STEExperiment::clone() const {

src/experiments/ste.h

@@ -13,7 +13,7 @@ public:
                const std::unordered_map<std::string, double>& optional) override;
     [[nodiscard]] double tmax() const override;
     void accumulate(const Trajectory& traj, double init_omega, int num_walkers) override;
-    void save(const motions::BaseMotion& motion, const times::BaseDistribution& dist) override;
+    void save(const std::string& directory) override;
     [[nodiscard]] std::unique_ptr<Experiment> clone() const override;
     void merge(const Experiment& other) override;
 

src/main.cpp

@@ -1,5 +1,6 @@
 
 #include "sims.h"
+#include "dynamics/decoupled.h"
 #include "experiments/spectrum.h"
 #include "experiments/ste.h"
 #include "utils/io.h"
@@ -43,16 +44,17 @@ int main (const int argc, char *argv[]) {
         rng.seed(rd());
     }
 
-    auto motion = motions::BaseMotion::createFromInput(args.motion_type);
-    auto dist = times::BaseDistribution::createFromInput(args.distribution_type);
+    DecoupledDynamics dynamics(
+        motions::BaseMotion::createFromInput(args.motion_type),
+        times::BaseDistribution::createFromInput(args.distribution_type));
 
     if (args.spectrum) {
         SpectrumExperiment experiment;
-        run_simulation(experiment, parameter, args.optional, *motion, *dist, rng);
+        run_simulation(experiment, parameter, args.optional, dynamics, rng);
     }
     if (args.ste) {
         STEExperiment experiment;
-        run_simulation(experiment, parameter, args.optional, *motion, *dist, rng);
+        run_simulation(experiment, parameter, args.optional, dynamics, rng);
     }
 
     return 0;

src/sims.cpp

@@ -1,5 +1,6 @@
 #include "sims.h"
 #include "utils/functions.h"
+#include "utils/io.h"
 
 #include <iostream>
 #include <chrono>
@@ -10,14 +11,12 @@ void run_simulation(
     Experiment& experiment,
     std::unordered_map<std::string, double>& parameter,
     std::unordered_map<std::string, double>& optional,
-    motions::BaseMotion& motion,
-    times::BaseDistribution& dist,
+    Dynamics& dynamics,
     std::mt19937_64& rng
     ) {
     const int num_walker = static_cast<int>(parameter["num_walker"]);
 
-    dist.setParameters(parameter);
-    motion.setParameters(parameter);
+    dynamics.setParameters(parameter);
     experiment.setup(parameter, optional);
 
     const auto start = printStart(optional);
@@ -31,13 +30,11 @@ void run_simulation(
         thread_rngs.emplace_back(rng());
     }
 
-    // Create per-thread clones of motion, distribution, and experiment
-    std::vector<std::unique_ptr<motions::BaseMotion>> thread_motions;
-    std::vector<std::unique_ptr<times::BaseDistribution>> thread_dists;
+    // Create per-thread clones of dynamics and experiment
+    std::vector<std::unique_ptr<Dynamics>> thread_dynamics;
     std::vector<std::unique_ptr<Experiment>> thread_experiments;
     for (int i = 0; i < num_threads; i++) {
-        thread_motions.push_back(motion.clone());
-        thread_dists.push_back(dist.clone());
+        thread_dynamics.push_back(dynamics.clone());
         thread_experiments.push_back(experiment.clone());
     }
 
@@ -48,14 +45,13 @@ void run_simulation(
     {
         const int tid = omp_get_thread_num();
         auto& local_rng = thread_rngs[tid];
-        auto& local_motion = *thread_motions[tid];
-        auto& local_dist = *thread_dists[tid];
+        auto& local_dynamics = *thread_dynamics[tid];
         auto& local_experiment = *thread_experiments[tid];
 
         #pragma omp for schedule(static)
         for (int mol_i = 0; mol_i < num_walker; mol_i++) {
-            auto traj = make_trajectory(local_motion, local_dist, experiment.tmax(), local_rng);
-            local_experiment.accumulate(traj, local_motion.getInitOmega(), num_walker);
+            auto traj = make_trajectory(local_dynamics, experiment.tmax(), local_rng);
+            local_experiment.accumulate(traj, local_dynamics.getInitOmega(), num_walker);
 
             if (tid == 0) {
                 #pragma omp atomic
@@ -73,24 +69,23 @@ void run_simulation(
         experiment.merge(*thread_experiments[i]);
     }
 
-    experiment.save(motion, dist);
+    const auto directory = make_directory(dynamics.toString());
+    experiment.save(directory);
     printEnd(start);
 }
 
 
 Trajectory make_trajectory(
-    motions::BaseMotion& motion,
-    times::BaseDistribution& dist,
+    Dynamics& dynamics,
     const double t_max,
     std::mt19937_64& rng
     ) {
     double t_passed = 0;
     double phase = 0;
 
-    motion.initialize(rng);
-    dist.initialize(rng);
+    dynamics.initialize(rng);
 
-    double omega = motion.getInitOmega();
+    double omega = dynamics.getInitOmega();
 
     Trajectory traj;
     traj.time.emplace_back(t_passed);
@@ -98,10 +93,10 @@ Trajectory make_trajectory(
     traj.omega.emplace_back(omega);
 
     while (t_passed < t_max) {
-        const double t = dist.tau_wait(rng);
-        t_passed += t;
-        omega = motion.jump(rng);
-        phase += omega * t;
+        auto [dt, new_omega] = dynamics.next(rng);
+        phase += omega * dt;
+        t_passed += dt;
+        omega = new_omega;
 
         traj.time.emplace_back(t_passed);
         traj.phase.emplace_back(phase);

src/sims.h

@@ -1,9 +1,8 @@
 #ifndef RWSIM_SIMS_H
 #define RWSIM_SIMS_H
 
+#include "dynamics/base.h"
 #include "experiments/base.h"
-#include "motions/base.h"
-#include "times/base.h"
 
 #include <unordered_map>
 #include <string>
@@ -14,11 +13,10 @@ void run_simulation(
     Experiment& experiment,
     std::unordered_map<std::string, double>& parameter,
     std::unordered_map<std::string, double>& optional,
-    motions::BaseMotion& motion,
-    times::BaseDistribution& dist,
+    Dynamics& dynamics,
     std::mt19937_64& rng);
 
-Trajectory make_trajectory(motions::BaseMotion& motion, times::BaseDistribution& dist, double t_max, std::mt19937_64& rng);
+Trajectory make_trajectory(Dynamics& dynamics, double t_max, std::mt19937_64& rng);
 
 std::chrono::system_clock::time_point printStart(std::unordered_map<std::string, double> &optional);
 void printEnd(std::chrono::system_clock::time_point start);

src/utils/io.cpp

@@ -11,11 +11,6 @@
 #include <string>
 #include <filesystem>
 
-#include "../motions/base.h"
-#include "../times/base.h"
-
-
-class Motion;
 
 std::pair<std::string, double> get_optional_parameter(std::vector<std::string>::const_iterator &it) {
     std::string stripped_arg;
@@ -128,19 +123,12 @@ std::unordered_map<std::string, double> read_parameter(const std::filesystem::pa
 }
 
 
-std::string make_directory(
-    const motions::BaseMotion& motion,
-    const times::BaseDistribution& distribution
-    ) {
-    std::ostringstream path_name;
-    path_name << motion.toString() << "/" << distribution.toString();
-
-    if (!std::filesystem::create_directories(path_name.str())) {
-        std::cout << "Created directory " << path_name.str() << std::endl;
+std::string make_directory(const std::string& path_name) {
+    if (!std::filesystem::create_directories(path_name)) {
+        std::cout << "Created directory " << path_name << std::endl;
     }
 
-
-    return path_name.str();
+    return path_name;
 }
 
 

src/utils/io.h

@@ -7,9 +7,6 @@
 #include <filesystem>
 #include <vector>
 
-#include "../motions/base.h"
-#include "../times/base.h"
-
 struct Arguments {
     std::string parameter_file{};
     bool ste = false;
@@ -23,7 +20,7 @@ Arguments parse_args(int argc, char* argv[]);
 std::pair<std::string, double> get_optional_parameter(std::vector<std::string>::const_iterator &it);
 std::unordered_map<std::string, double> read_parameter(const std::filesystem::path&);
 
-std::string make_directory(const motions::BaseMotion&, const times::BaseDistribution&);
+std::string make_directory(const std::string& path_name);
 
 void save_parameter_to_file(const std::string&, const std::string&, const std::unordered_map<std::string, double>&, const std::unordered_map<std::string, double>&);
 void save_data_to_file(const std::string&, const std::string&, const std::vector<double>&, const std::map<double, std::vector<double>>&, const std::unordered_map<std::string, double>&);