simul/DESIGN.md

Simul - A Unified Simulation Framework in Rust

A modular Rust framework for molecular dynamics, stochastic simulation, and Monte Carlo methods, supporting:

• Variable dimensions (1D, 2D, 3D, ...) via const generics
• Continuous and discrete state spaces
• Deterministic and stochastic dynamics
• Molecular dynamics, Brownian motion, and lattice models

---

Table of Contents

1. Workspace Structure
2. Crate Responsibilities
3. Core Abstractions
4. State Space System
5. Dynamics System
6. Interaction System
7. Simulation API
8. Dynamic Dispatch
9. Error Handling
10. Platform Abstraction

---

Workspace Structure

simul/                             # Workspace root (renamed from rust-openmm)
├── Cargo.toml                     # Workspace manifest
├── DESIGN.md                      # This document
├── README.md
│
├── simul-core/                    # Core traits and abstractions
│   ├── Cargo.toml
│   └── src/
│       ├── lib.rs
│       ├── space/                 # StateSpace trait
│       │   ├── mod.rs
│       │   └── traits.rs
│       ├── dynamics/              # Dynamics trait
│       │   └── mod.rs
│       ├── interaction/           # Interaction trait
│       │   └── mod.rs
│       ├── state.rs               # SystemState<S>
│       ├── units.rs               # Physical units
│       ├── error.rs               # Error types
│       └── prelude.rs
│
├── simul-euclidean/               # Continuous Euclidean space simulations
│   ├── Cargo.toml                 # depends on simul-core
│   └── src/
│       ├── lib.rs
│       ├── space.rs               # Euclidean<D> implementation
│       ├── dynamics/
│       │   ├── mod.rs
│       │   ├── verlet.rs          # Velocity Verlet (MD)
│       │   ├── langevin.rs        # Langevin equation
│       │   └── brownian.rs        # Overdamped Brownian
│       ├── interaction/
│       │   ├── mod.rs
│       │   ├── nonbonded.rs       # L-J, Coulomb
│       │   ├── harmonic_bond.rs   # Bonds
│       │   ├── harmonic_angle.rs  # Angles
│       │   └── external.rs        # External fields
│       └── prelude.rs
│
├── simul-lattice/                 # Discrete lattice space simulations
│   ├── Cargo.toml                 # depends on simul-core
│   └── src/
│       ├── lib.rs
│       ├── space.rs               # Lattice<D>, SpinLattice<D>
│       ├── dynamics/
│       │   ├── mod.rs
│       │   ├── metropolis.rs      # Metropolis MC
│       │   ├── glauber.rs         # Glauber dynamics
│       │   ├── gillespie.rs       # Continuous-time MC
│       │   └── kawasaki.rs        # Kawasaki dynamics (conserved)
│       ├── interaction/
│       │   ├── mod.rs
│       │   ├── ising.rs           # Ising model
│       │   ├── potts.rs           # Potts model
│       │   ├── xy.rs              # XY model (future)
│       │   └── lattice_gas.rs     # Lattice gas
│       └── prelude.rs
│
├── simul-sampling/                # Advanced sampling methods (future)
│   ├── Cargo.toml
│   └── src/
│       ├── replica_exchange.rs    # Parallel tempering
│       ├── umbrella.rs            # Umbrella sampling
│       ├── wang_landau.rs         # Wang-Landau
│       └── metadynamics.rs        # Metadynamics
│
├── simul/                         # Facade crate - re-exports all
│   ├── Cargo.toml                 # depends on all simul-* crates
│   └── src/
│       ├── lib.rs
│       └── prelude.rs             # use simul::prelude::*;
│
└── examples/
    ├── hello_argon.rs             # 3D MD with L-J
    ├── brownian_1d.rs             # 1D Brownian motion
    ├── brownian_2d.rs             # 2D Brownian motion
    ├── ising_2d.rs                # 2D Ising model
    ├── random_walk_2d.rs          # 2D lattice random walk
    └── langevin_particle.rs       # Langevin dynamics

---

Crate Responsibilities

simul-core

Purpose: Core traits and abstractions shared by all simulation types.

Contains:

• StateSpace trait - defines the geometric structure
• Dynamics trait - defines time evolution rules
• Interaction trait - defines driving forces/rates
• SystemState<S> - generic system state
• Unit types (Dalton, Nanometer, Picosecond, etc.)
• Error types

Dependencies:

• nalgebra - linear algebra
• thiserror - error handling
• rand - random number generation

simul-euclidean

Purpose: Simulations in continuous Euclidean space.

Contains:

• Euclidean<D> space implementation
• Dynamics:
  • VerletDynamics<D> - velocity Verlet (deterministic MD)
  • LangevinDynamics<D> - Langevin equation (stochastic MD)
  • BrownianDynamics<D> - overdamped/Smoluchowski (no inertia)
• Interactions:
  • NonbondedInteraction<D> - Lennard-Jones, Coulomb
  • HarmonicBondInteraction<D> - harmonic bonds
  • HarmonicAngleInteraction<D> - harmonic angles
  • ExternalFieldInteraction<D> - external potentials

Use Cases:

• Molecular dynamics (proteins, materials)
• Brownian motion (colloids, polymers)
• Langevin dynamics (thermostated systems)

Dependencies:

• simul-core
• cubecl (optional, for all backends: CPU SIMD, CUDA, wgpu, Vulkan)

simul-lattice

Purpose: Simulations on discrete lattice spaces.

Contains:

• Lattice<D> space implementation
• SpinLattice<D> for spin models
• Dynamics:
  • MetropolisDynamics<D> - Metropolis-Hastings MC
  • GlauberDynamics<D> - Glauber (heat bath) dynamics
  • GillespieDynamics<D> - continuous-time MC (kinetic MC)
  • KawasakiDynamics<D> - conserved order parameter
• Interactions:
  • IsingInteraction<D> - Ising model
  • PottsInteraction<D> - q-state Potts model
  • LatticeGasInteraction<D> - lattice gas model

Use Cases:

• Ising model (magnetism, phase transitions)
• Lattice random walks
• Lattice gas simulations
• Kinetic Monte Carlo

Dependencies:

• simul-core

simul-sampling (Future)

Purpose: Advanced sampling and enhanced sampling methods.

Contains:

• ReplicaExchange - parallel tempering
• UmbrellaSampling - biased sampling
• WangLandau - flat histogram methods
• Metadynamics - adaptive bias

Dependencies:

• simul-core
• May depend on simul-euclidean or simul-lattice

simul (Facade)

Purpose: Convenience crate that re-exports everything.

// Users can simply:
use simul::prelude::*;

// Instead of:
use simul_core::prelude::*;
use simul_euclidean::prelude::*;
use simul_lattice::prelude::*;

---

Core Abstractions

StateSpace Trait (simul-core)

/// Defines the geometric/topological structure of the state space
pub trait StateSpace: Clone + Send + Sync + 'static {
    /// Compile-time dimension
    const DIM: usize;
    
    /// Type representing a point in this space
    type Point: Clone + Send + Sync + Default + Debug;
    
    /// Type representing momentum/velocity (use () for discrete systems)
    type Momentum: Clone + Send + Sync + Default + Debug;
    
    /// Whether this is a continuous space
    fn is_continuous() -> bool;
    
    /// Human-readable name
    fn name() -> &'static str;
}

Dynamics Trait (simul-core)

/// Defines how the system evolves in time
pub trait Dynamics<S: StateSpace>: Send + Sync + Debug {
    /// Advance system by one step, return time increment
    fn step<R: Rng>(&self, state: &mut SystemState<S>, rng: &mut R) -> f64;
    
    /// Name of the dynamics
    fn name(&self) -> &'static str;
    
    /// Whether dynamics is deterministic
    fn is_deterministic(&self) -> bool;
    
    /// Time type (continuous or discrete)
    fn time_type(&self) -> TimeType;
}

pub enum TimeType {
    Continuous,  // Real-valued time (dt for MD, waiting time for CTMC)
    Discrete,    // Integer steps (MC sweeps)
}

Interaction Trait (simul-core)

/// Defines what drives state changes
pub trait Interaction<S: StateSpace>: Send + Sync + Debug {
    /// Output type of the computation
    type Output;
    
    /// Compute interaction from current state
    fn compute(&self, state: &SystemState<S>) -> Self::Output;
    
    /// Name of the interaction
    fn name(&self) -> &'static str;
}

---

State Space System

Euclidean (simul-euclidean)

use nalgebra::SVector;

/// D-dimensional Euclidean space
#[derive(Clone, Debug, Default)]
pub struct Euclidean<const D: usize>;

impl<const D: usize> StateSpace for Euclidean<D> {
    const DIM: usize = D;
    type Point = SVector<f64, D>;
    type Momentum = SVector<f64, D>;
    
    fn is_continuous() -> bool { true }
    fn name() -> &'static str { 
        match D {
            1 => "Euclidean 1D",
            2 => "Euclidean 2D", 
            3 => "Euclidean 3D",
            _ => "Euclidean ND",
        }
    }
}

// Convenient aliases
pub type Space1D = Euclidean<1>;
pub type Space2D = Euclidean<2>;
pub type Space3D = Euclidean<3>;

Lattice (simul-lattice)

/// D-dimensional integer lattice
#[derive(Clone, Debug)]
pub struct Lattice<const D: usize> {
    pub size: SVector<usize, D>,
    pub periodic: bool,
}

impl<const D: usize> StateSpace for Lattice<D> {
    const DIM: usize = D;
    type Point = SVector<i64, D>;
    type Momentum = ();  // No momentum in discrete space
    
    fn is_continuous() -> bool { false }
    fn name() -> &'static str { "Lattice" }
}

---

Dynamics System

Continuous Dynamics (simul-euclidean)

Dynamics	Equation	Deterministic?	Time
VerletDynamics<D>	Newton's equations	Yes	Continuous
LangevinDynamics<D>	mdv = (-∇U - γv)dt + √(2γkT)dW	No	Continuous
BrownianDynamics<D>	dx = (F/γ)dt + √(2kT/γ)dW	No	Continuous

Discrete Dynamics (simul-lattice)

Dynamics	Algorithm	Time
MetropolisDynamics<D>	Metropolis-Hastings	Discrete (sweeps)
GlauberDynamics<D>	Heat bath	Discrete
GillespieDynamics<D>	Gillespie SSA	Continuous (CTMC)
KawasakiDynamics<D>	Conserved exchange	Discrete

---

Simulation API

Builder Pattern

use simul::prelude::*;

// 3D Molecular Dynamics
let md_sim = Simulation::<Euclidean<3>>::builder()
    .state(initial_state)
    .dynamics(VerletDynamics::new(Picosecond(0.002)))
    .add_interaction(NonbondedInteraction::lj())
    .add_interaction(HarmonicBondInteraction::new())
    .seed(42)
    .build()?;

// 1D Brownian Motion
let brownian_sim = Simulation::<Euclidean<1>>::builder()
    .state(initial_state)
    .dynamics(BrownianDynamics::new(temperature, friction, dt))
    .add_interaction(HarmonicTrap::new(spring_constant))
    .seed(42)
    .build()?;

// 2D Ising Model
let ising_sim = Simulation::<Lattice<2>>::builder()
    .state(spin_state)
    .dynamics(MetropolisDynamics::new(temperature))
    .add_interaction(IsingInteraction::new(coupling))
    .build()?;

Iterator Pattern

// Collect trajectory
let trajectory: Vec<_> = simulation
    .run()
    .take(10000)
    .step_by(100)
    .collect();

// On-the-fly analysis
for state in simulation.run().take(1000) {
    println!("t={:.3}, E={:.2}", state.time, state.energy());
}

---

Dynamic Dispatch

For runtime polymorphism with const generics, use enum-based dispatch:

/// Type-erased state for heterogeneous collections
pub enum AnyState {
    Euclidean1D(SystemState<Euclidean<1>>),
    Euclidean2D(SystemState<Euclidean<2>>),
    Euclidean3D(SystemState<Euclidean<3>>),
    Lattice2D(SystemState<Lattice<2>>),
    Lattice3D(SystemState<Lattice<3>>),
}

/// Type-erased dynamics
pub enum AnyDynamics {
    // Euclidean
    Verlet3D(VerletDynamics<3>),
    Langevin3D(LangevinDynamics<3>),
    Brownian1D(BrownianDynamics<1>),
    Brownian2D(BrownianDynamics<2>),
    // Lattice
    Metropolis2D(MetropolisDynamics<2>),
    Glauber2D(GlauberDynamics<2>),
    Gillespie2D(GillespieDynamics<2>),
}

impl AnyDynamics {
    pub fn step(&self, state: &mut AnyState, rng: &mut impl Rng) -> Result<f64> {
        match (self, state) {
            (Self::Verlet3D(d), AnyState::Euclidean3D(s)) => Ok(d.step(s, rng)),
            (Self::Brownian2D(d), AnyState::Euclidean2D(s)) => Ok(d.step(s, rng)),
            (Self::Metropolis2D(d), AnyState::Lattice2D(s)) => Ok(d.step(s, rng)),
            _ => Err(SimulationError::TypeMismatch),
        }
    }
}

---

Error Handling

// simul-core/src/error.rs
use thiserror::Error;

#[derive(Debug, Error)]
pub enum SimulationError {
    #[error("Dimension mismatch: expected {expected}, got {actual}")]
    DimensionMismatch { expected: usize, actual: usize },
    
    #[error("Space type mismatch: {message}")]
    TypeMismatch { message: String },
    
    #[error("Energy diverged: {value:.2e}")]
    EnergyDiverged { value: f64 },
    
    #[error("Invalid parameter: {name} = {value}")]
    InvalidParameter { name: &'static str, value: String },
    
    #[error("Builder incomplete: {field} not set")]
    BuilderIncomplete { field: &'static str },
}

---

Platform Abstraction (CubeCL All-In)

Design Philosophy

All compute kernels target CubeCL as the single runtime abstraction. A kernel written once with #[cube(launch)] runs on every supported backend: cubecl-cpu, CUDA, wgpu, and Vulkan. The user selects a backend via a Cargo feature flag; no application code changes are needed.

Using cubecl 0.10.0-pre.2 which includes cubecl-cpu for CPU LLVM SIMD vectorization. The cpu path currently uses a native Rust serial fallback while CubeCL kernels are being implemented.

Feature Flags

# simul-euclidean/Cargo.toml
[features]
default = ["cpu"]
cpu = ["dep:cubecl", "cubecl/cpu"]                    # CubeCL → CPU (LLVM SIMD)
cuda = ["dep:cubecl", "cubecl/cuda"]                  # CubeCL → NVIDIA PTX
wgpu = ["dep:cubecl", "cubecl/wgpu"]                  # CubeCL → WGSL (Metal/DX12)
vulkan = ["dep:cubecl", "cubecl/vulkan"]              # CubeCL → SPIR-V (Vulkan)

cargo build                          # default: cubecl-cpu (LLVM SIMD)
cargo build --features cuda          # NVIDIA GPU
cargo build --features wgpu          # wgpu (WGSL)
cargo build --features vulkan        # Vulkan (SPIR-V)

Architecture Layers

┌─────────────────────────────────────────────────────────────┐
│  User code / Dynamics (Verlet, Langevin, Brownian)          │
│  Calls ForceInteraction::compute_forces()                   │
│  *** Knows nothing about which backend is running ***       │
├─────────────────────────────────────────────────────────────┤
│  ForceInteraction trait                                      │
│  fn compute_forces(&[SVector], &mut [SVector]) -> f64       │
│  Single impl per interaction — internal #[cfg] dispatch     │
├─────────────────────────────────────────────────────────────┤
│  #[cfg] dispatch inside compute_forces()                    │
│                                                              │
│  ┌──────────────────────────────────────────────────────┐   │
│  │  CubeCL kernels — #[cube(launch)]                    │   │
│  │  Single kernel runs on ALL backends:                  │   │
│  │  ┌────────┬──────┬──────────┬──────────────┐         │   │
│  │  │ CPU    │ CUDA │ wgpu     │ Vulkan       │         │   │
│  │  │ (SIMD) │(PTX) │ (WGSL)  │ (SPIR-V)     │         │   │
│  │  └────────┴──────┴──────────┴──────────────┘         │   │
│  └──────────────────────────────────────────────────────┘   │
│                                                              │
│  ┌──────────────┐  (fallback when no feature is active)     │
│  │  Native Rust  │                                           │
│  │  serial loops │                                           │
│  └──────────────┘                                            │
└─────────────────────────────────────────────────────────────┘

CubeCL Backend Matrix

Feature	CubeCL Backend	Target	f64	Status
cpu	cubecl-cpu	CPU LLVM SIMD	yes	ready
cuda	cubecl-cuda	NVIDIA PTX	yes	planned
wgpu	cubecl-wgpu	WGSL (Metal/DX12)	no	planned
vulkan	cubecl-wgpu	SPIR-V (Vulkan)	no	planned

Using cubecl 0.10.0-pre.2 (pre-release). For MD simulations requiring f64, cpu and cuda backends are recommended. The wgpu/vulkan backends may be limited to f32 depending on the WGSL/SPIR-V target.

Kernel Design: Write Once, Run Everywhere

Each force computation is written once as a CubeCL kernel. The same kernel runs on CUDA, wgpu, or Vulkan without modification:

// interaction/kernels.rs — compiled only with cuda/wgpu features
use cubecl::prelude::*;

#[cube(launch)]
fn nonbonded_lj_coulomb<F: Float>(
    positions: &Array<F>,
    forces: &mut Array<F>,
    sigma: &Array<F>,
    epsilon: &Array<F>,
    charge: &Array<F>,
    coulomb_k: F,
    n_particles: u32,
    dim: u32,
) {
    let i = ABSOLUTE_POS;
    if i >= n_particles { return; }
    // ... pairwise LJ + Coulomb — identical code for all backends
}

Interaction Pattern: Single Impl, Internal Dispatch

Each interaction has one ForceInteraction impl. The #[cfg] dispatch is hidden inside the method body:

impl<const D: usize> ForceInteraction<D> for NonbondedInteraction<D> {
    fn compute_forces(&self, positions: &[SVector<f64, D>],
                      forces: &mut [SVector<f64, D>]) -> f64 {
        #[cfg(any(feature = "cuda", feature = "wgpu", feature = "vulkan"))]
        { return self.compute_cubecl(positions, forces); }

        #[cfg(not(any(feature = "cuda", feature = "wgpu", feature = "vulkan")))]
        { self.compute_cpu(positions, forces) }
    }
}

The public API (ForceInteraction) never exposes platform details. Dynamics, tests, and user code all call the same compute_forces().

Physics Parameters vs Computation

Each interaction separates platform-independent parameters from platform-specific computation:

// Platform-independent — shared by all backends
pub struct NonbondedParticleParams {
    pub sigma: f64,
    pub epsilon: f64,
    pub charge: f64,
}

// Interaction holds params; compute is dispatched internally
pub struct NonbondedInteraction<const D: usize> {
    particles: Vec<NonbondedParticleParams>,
    coulomb_constant: f64,
    cutoff: Option<f64>,
}

Roadmap: Full CubeCL Kernel Unification

The cubecl-cpu dependency is in place (v0.10.0-pre.2). Once the #[cube(launch)] kernels are implemented, a single kernel will run on all backends (CPU, CUDA, wgpu, Vulkan) — zero code duplication. The native Rust serial impl will remain as a no-dependency fallback for testing and debugging.

---

Dependency Graph

                    ┌─────────────┐
                    │  simul-core │
                    │  (traits)   │
                    └──────┬──────┘
                           │
           ┌───────────────┼───────────────┐
           │               │               │
           ▼               ▼               ▼
    ┌─────────────┐ ┌─────────────┐ ┌─────────────┐
    │   simul-    │ │   simul-    │ │   simul-    │
    │  euclidean  │ │   lattice   │ │  sampling   │
    │  (MD, BD)   │ │(Ising, RW)  │ │  (RE, WL)   │
    └──────┬──────┘ └──────┬──────┘ └──────┬──────┘
           │               │               │
           └───────────────┼───────────────┘
                           │
                           ▼
                    ┌─────────────┐
                    │    simul    │
                    │  (facade)   │
                    └─────────────┘

---

Summary of Design Decisions

Aspect	Decision	Rationale
Project Name	simul	Short, memorable, general
Workspace	Multi-crate	Modularity, selective deps
Core Crate	simul-core	Shared traits, no physics
Euclidean	simul-euclidean	MD, Brownian, Langevin
Lattice	simul-lattice	Ising, MC, random walk
Dimension	Const generic	Zero-cost, type-safe
Dispatch	Enum-based	Runtime polymorphism
API	Builder + Iterator	Ergonomic, Rust-idiomatic
Errors	thiserror	Typed, descriptive
Platform	CubeCL all-in	Single kernel → all backends (CPU SIMD/CUDA/wgpu/Vulkan)