Struct ModelRealization 

pub struct ModelRealization<S: Spin> {
    pub couplings: Vec<S::Value>,
    pub spins: Vec<S>,
    pub temperatures: Vec<S::Value>,
    pub system_ids: Vec<usize>,
    pub rngs: Vec<Xoshiro256StarStar>,
    pub pair_rngs: Vec<Xoshiro256StarStar>,
    pub energies: Vec<S::Value>,
    /* private fields */
}

Mutable state for one disorder realization.

Holds the coupling array (fixed after construction), spin configurations for every replica at every temperature, and bookkeeping for parallel tempering.

With n_replicas replicas and n_temps temperatures there are n_systems = n_replicas * n_temps independent spin configurations. Spins are stored in a single flat Vec of length n_systems * n_spins, where system i occupies spins[i*n_spins .. (i+1)*n_spins].

Fields

couplings: Vec<S::Value>

Forward couplings, length n_spins * n_neighbors.

spins: Vec<S>

All spin configurations, length n_systems * n_spins (±1 for Ising, unit vectors for XY).

temperatures: Vec<S::Value>

Temperature assigned to each system slot, length n_systems.

system_ids: Vec<usize>

Parallel-tempering permutation: system_ids[slot] is the system index currently occupying temperature slot slot.

rngs: Vec<Xoshiro256StarStar>

One PRNG per system.

pair_rngs: Vec<Xoshiro256StarStar>

One PRNG per overlap-update pair slot, length n_temps * (n_replicas / 2).

energies: Vec<S::Value>

Scratch cache of the interaction energy per system (-H/N), length n_systems. It is current only right after a measurement or tempering step.

Implementations

impl<S: Spin> ModelRealization<S>

pub fn new( lattice: &Lattice, couplings: Vec<S::Value>, temps: &[S::Value], n_replicas: usize, base_seed: u64, ) -> Self

Initialize a realization with random spins (±1, or uniform unit vectors for XY).

Seeds independent system and pair streams from domain-separated child seeds.

Examples found in repository

examples/bench.rs (lines 93-95)

fn main() {
    let l = env_usize("PEAPODS_L", 128);
    let n_dims = env_usize("PEAPODS_DIM", 2);
    let n_temps = env_usize("PEAPODS_TEMPS", 16);
    let n_replicas = env_usize("PEAPODS_REPLICAS", 2);
    let n_sweeps = env_usize("PEAPODS_SWEEPS", 50);
    let n_realizations = env_usize("PEAPODS_NREAL", 100);
    let sequential = env::var_os("PEAPODS_SEQUENTIAL").is_some();
    let overlap_cluster_mode = if env::var_os("PEAPODS_OVERLAP_WOLFF").is_some() {
        ClusterMode::Wolff
    } else {
        ClusterMode::Sw
    };
    let generic_lattice = env::var_os("PEAPODS_GENERIC_LATTICE").is_some();
    let mode = env::var("PEAPODS_MODE").unwrap_or_else(|_| "cmr".to_string());
    let couplings_kind = env::var("PEAPODS_COUPLINGS").unwrap_or_else(|_| "bimodal".to_string());
    let sweep_mode = match env::var("PEAPODS_SWEEP").as_deref().unwrap_or("metropolis") {
        "metropolis" => SweepMode::Metropolis,
        "gibbs" => SweepMode::Gibbs,
        other => panic!("unknown PEAPODS_SWEEP '{other}'"),
    };
    let shape = vec![l; n_dims];
    let lattice = if generic_lattice {
        Lattice::with_offsets(shape, hypercubic(n_dims))
    } else {
        Lattice::new(shape)
    };
    let n_spins = lattice.n_spins;
    let n_neighbors = lattice.n_neighbors;

    let t_min: f32 = env::var("PEAPODS_TMIN").map_or(0.1, |v| v.parse().unwrap());
    let t_max: f32 = env::var("PEAPODS_TMAX").map_or(5.0, |v| v.parse().unwrap());
    let temps: Vec<f32> = (0..n_temps)
        .map(|i| t_min * (t_max / t_min).powf(i as f32 / (n_temps.max(2) - 1) as f32))
        .collect();

    let n_pairs = n_replicas / 2;
    let n_systems = n_replicas * n_temps;
    let rngs_per_real = n_systems + n_temps * n_pairs;

    let mut rng = Xoshiro256StarStar::seed_from_u64(0x5eed);
    let mut realizations = Vec::with_capacity(n_realizations);
    for r in 0..n_realizations {
        let couplings: Vec<f32> = (0..n_spins * n_neighbors)
            .map(|_| match couplings_kind.as_str() {
                "bimodal" => {
                    if rng.gen::<bool>() {
                        1.0
                    } else {
                        -1.0
                    }
                }
                "ferro" => 1.0,
                // Box-Muller keeps the bench free of extra distribution crates.
                "gaussian" => {
                    let u1: f64 = 1.0 - rng.gen::<f64>();
                    let u2: f64 = rng.gen();
                    ((-2.0 * u1.ln()).sqrt() * (std::f64::consts::TAU * u2).cos()) as f32
                }
                other => panic!("unknown PEAPODS_COUPLINGS '{other}'"),
            })
            .collect();
        let base_seed = 42 + (r * rngs_per_real) as u64;
        realizations.push(Realization::new(
            &lattice, couplings, &temps, n_replicas, base_seed,
        ));
    }

    let interrupted = AtomicBool::new(false);
    let (cluster_update, pt_interval, overlap_cluster) = match mode.as_str() {
        "metropolis" => (None, None, None),
        "pt" => (None, Some(1), None),
        "sw" => (
            Some(ClusterConfig {
                interval: 1,
                mode: ClusterMode::Sw,
                action: ClusterAction::Update,
                collect_stats: false,
            }),
            None,
            None,
        ),
        "wolff" => (
            Some(ClusterConfig {
                interval: 1,
                mode: ClusterMode::Wolff,
                action: ClusterAction::Update,
                collect_stats: false,
            }),
            None,
            None,
        ),
        "cmr" => (
            None,
            Some(1),
            Some(OverlapClusterConfig {
                interval: 1,
                modes: vec![OverlapClusterBuildMode::Cmr],
                cluster_mode: overlap_cluster_mode,
                action: ClusterAction::Update,
                collect_stats: false,
                snapshot_interval: None,
                max_temperature: None,
            }),
        ),
        "sw_pt" => (
            Some(ClusterConfig {
                interval: 1,
                mode: ClusterMode::Sw,
                action: ClusterAction::Update,
                collect_stats: false,
            }),
            Some(1),
            None,
        ),
        "houdayer" | "jorg" => (
            None,
            Some(1),
            Some(OverlapClusterConfig {
                interval: 1,
                modes: vec![if mode == "jorg" {
                    OverlapClusterBuildMode::Jorg(2)
                } else {
                    OverlapClusterBuildMode::Houdayer(2)
                }],
                cluster_mode: overlap_cluster_mode,
                action: ClusterAction::Update,
                collect_stats: false,
                snapshot_interval: None,
                max_temperature: None,
            }),
        ),
        // Any other overlap build mode string, e.g. "pair4", "jorg4" or "rmc", with PT.
        other => (
            None,
            Some(1),
            Some(OverlapClusterConfig {
                interval: 1,
                modes: parse_overlap_modes(other)
                    .unwrap_or_else(|err| panic!("unknown PEAPODS_MODE '{other}': {err}")),
                cluster_mode: overlap_cluster_mode,
                action: ClusterAction::Update,
                collect_stats: false,
                snapshot_interval: None,
                max_temperature: None,
            }),
        ),
    };

    let config = SimConfig {
        n_sweeps,
        warmup_sweeps: 0,
        sweep_mode,
        cluster_update,
        pt_interval,
        pt_schedule: PtSchedule::SingleRandomEdge,
        overlap_cluster,
        autocorrelation_max_lag: None,
        autocorrelation_backend: AutocorrelationBackend::Ring,
        sequential,
        equilibration_diagnostic: false,
        conditional_link_overlap: true,
    };

    println!(
        "Lattice: {l}^{n_dims}  |  Temps: {n_temps}  |  Replicas: {n_replicas}  |  Sweeps: {n_sweeps}  |  Realizations: {n_realizations}"
    );
    println!(
        "Config: {couplings_kind}, mode={mode}, sequential={sequential}, generic_lattice={generic_lattice}"
    );
    println!("{}", "-".repeat(70));

    let t0 = Instant::now();
    let result = run_sweep_parallel(
        &lattice,
        &mut realizations,
        n_replicas,
        n_temps,
        &config,
        &interrupted,
        &|| {},
    )
    .unwrap();
    let elapsed = t0.elapsed().as_secs_f64();

    let mut state_hash = DefaultHasher::new();
    for realization in &realizations {
        realization.spins.hash(&mut state_hash);
        realization.system_ids.hash(&mut state_hash);
        for rng in realization.rngs.iter().chain(&realization.pair_rngs) {
            let mut rng = rng.clone();
            rng.next_u64().hash(&mut state_hash);
        }
    }
    for values in [
        &result.mags,
        &result.mags2,
        &result.mags4,
        &result.energies,
        &result.energies2,
    ] {
        for value in values {
            value.to_bits().hash(&mut state_hash);
        }
    }
    for values in [
        &result.overlap_stats.overlap,
        &result.overlap_stats.overlap2,
        &result.overlap_stats.overlap4,
        &result.overlap_stats.link_overlap,
        &result.overlap_stats.link_overlap2,
        &result.overlap_stats.link_overlap4,
    ] {
        for value in values {
            value.to_bits().hash(&mut state_hash);
        }
    }
    result.overlap_stats.histogram.hash(&mut state_hash);
    result
        .overlap_stats
        .per_sample_histogram
        .hash(&mut state_hash);
    for samples in [
        &result.overlap_stats.ql_at_q_sum,
        &result.overlap_stats.ql2_at_q_sum,
    ] {
        for bins in samples {
            for value in bins {
                value.to_bits().hash(&mut state_hash);
            }
        }
    }
    for samples in [
        &result.overlap_stats.per_sample_ql_at_q_sum,
        &result.overlap_stats.per_sample_ql2_at_q_sum,
    ] {
        for temperatures in samples {
            for bins in temperatures {
                for value in bins {
                    value.to_bits().hash(&mut state_hash);
                }
            }
        }
    }

    let per_sweep = elapsed / n_sweeps as f64 * 1000.0;
    println!("Total: {:.3} s  |  {:.3} ms/sweep", elapsed, per_sweep);
    println!("State checksum: {:016x}", state_hash.finish());
}

pub fn reset( &mut self, lattice: &Lattice, n_replicas: usize, n_temps: usize, base_seed: u64, )

Re-randomize all spins and reset the tempering permutation.

pub fn validate( &self, lattice: &Lattice, n_replicas: usize, n_temps: usize, ) -> Result<(), String>

Validate externally mutable buffers before entering parallel kernels.

pub fn pt_edge_attempts(&self) -> &[u64]

Swap attempts per temperature edge since creation or the last reset. Like the acceptances and round trips, they count warmup sweeps too.

pub fn pt_edge_acceptances(&self) -> &[u64]

pub fn pt_round_trips(&self) -> &[u64]

pub fn pt_flow_up(&self) -> &[u64]

Per temperature slot, visits after each tempering event by systems that last visited the cold end of the ladder.

pub fn pt_flow_labelled(&self) -> &[u64]

Per temperature slot, visits by systems that have visited either end; the denominator of the flow fraction pt_flow_up / pt_flow_labelled.