From 6134638c28e5b8dd3632fe29f55ad48cce1870c7 Mon Sep 17 00:00:00 2001 From: Dimitri Bonanni-Surprenant Date: Thu, 27 Mar 2025 09:41:10 -0400 Subject: [PATCH 1/4] (fix) MonteCarlo systematic error. A systematic error was made because we did not count the first element of the Markov chain into the computation of the energy. --- src/bin/dvmc_pairwf.rs | 2 +- src/bin/dvmc_pairwf2.rs | 195 ++++++++++++++++++++++++++++++++ src/constants.rs | 2 +- src/dvmc.rs | 243 ++++++++++++++++++++++++++++++++++++++++ src/fock_state.rs | 3 +- src/lib.rs | 11 +- src/monte_carlo.rs | 10 +- 7 files changed, 455 insertions(+), 11 deletions(-) create mode 100644 src/bin/dvmc_pairwf2.rs create mode 100644 src/dvmc.rs diff --git a/src/bin/dvmc_pairwf.rs b/src/bin/dvmc_pairwf.rs index 3e0af27..5e2371c 100644 --- a/src/bin/dvmc_pairwf.rs +++ b/src/bin/dvmc_pairwf.rs @@ -16,7 +16,7 @@ const NVIJ: usize = SIZE*(SIZE - 1) / 2; const NGI: usize = SIZE; const NPARAMS: usize = NFIJ + NGI + NVIJ; const NELEC: usize = SIZE; -const NMCSAMP: usize = 1_000; +const NMCSAMP: usize = 1000; const NBOOTSTRAP: usize = 1; const NMCWARMUP: usize = 1000; const MCSAMPLE_INTERVAL: usize = 1; diff --git a/src/bin/dvmc_pairwf2.rs b/src/bin/dvmc_pairwf2.rs new file mode 100644 index 0000000..d3e3196 --- /dev/null +++ b/src/bin/dvmc_pairwf2.rs @@ -0,0 +1,195 @@ +//use blas::dcopy; +//use log::{debug, info}; +use rand_mt::Mt64; +use rand::Rng; +//use indicatif::{ProgressBar, ProgressStyle}; +use std::fs::File; +use std::io::Write; + +use impurity::{VarParams, SysParams, generate_bitmask, FockState, RandomStateGeneration}; +use impurity::dvmc::{variationnal_monte_carlo, EnergyOptimisationMethod, EnergyComputationMethod, VMCParams}; + +const SEED: u64 = 14; +const SIZE: usize = 16; +const NFIJ: usize = 4*SIZE*SIZE; +const NVIJ: usize = SIZE*(SIZE - 1) / 2; +const NGI: usize = SIZE; +const NPARAMS: usize = NFIJ + NGI + NVIJ; +const NELEC: usize = SIZE; +const NMCSAMP: usize = 1_000; +const NBOOTSTRAP: usize = 1; +const NMCWARMUP: usize = 1000; +const MCSAMPLE_INTERVAL: usize = 4; +const _NTHREADS: usize = 1; +const CLEAN_UPDATE_FREQUENCY: usize = 32; +const TOLERENCE_SHERMAN_MORRISSON: f64 = 1e-12; +const TOLERENCE_SINGULARITY: f64 = 1e-12; +const CONS_U: f64 = 1.0; +const CONS_T: f64 = 1.0; +const INITIAL_RATIO_UT: f64 = 1.0; +const FINAL_RATIO_UT: f64 = 32.0; +const NRATIO_POINTS: usize = 1; +const EPSILON_CG: f64 = 1e-16; +const EPSILON_SHIFT: f64 = 1e-2; +const OPTIMISATION_TIME_STEP: f64 = 1e-2; +const OPTIMISATION_DECAY: f64 = 0.0; +const NOPTITER: usize = 1000; +const KMAX: usize = NPARAMS; +const PARAM_THRESHOLD: f64 = ::MIN; +//const PARAM_THRESHOLD: f64 = 0.0; +const OPTIMISE: bool = true; +const OPTIMISE_GUTZ: bool = true; +const OPTIMISE_JAST: bool = true; +const OPTIMISE_ORB: bool = true; +const SET_EXPVALO_ZERO: bool = false; +const COMPUTE_ENERGY_METHOD: EnergyComputationMethod = EnergyComputationMethod::MonteCarlo; +const OPTIMISE_ENERGY_METHOD: EnergyOptimisationMethod = EnergyOptimisationMethod::ExactInverse; +const ENERGY_CONV_AVERAGE_SAMPLE: usize = 20; + +// 4 et 2 Sites +//pub const HOPPINGS: [f64; SIZE*SIZE] = [ +// //0.0, 1.0, 1.0, 0.0, +// 0.0, 1.0, 1.0, 0.0, +// 1.0, 0.0, 0.0, 1.0, +// 1.0, 0.0, 0.0, 1.0, +// 0.0, 1.0, 1.0, 0.0 +//]; + +// 8 Sites +//pub const HOPPINGS: [f64; SIZE*SIZE] = [ +// 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, +// 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, +// 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, +// 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, +// 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, +// 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, +// 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, +// 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, +//]; + +// 8 Sites +pub const HOPPINGS: [f64; SIZE*SIZE] = [ + 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, + 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, + 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, + 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, + 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, + 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, + 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, +]; + + +fn sq(a: f64) -> f64 { + ::abs(a) * ::abs(a) +} + +fn write_energy(fp: &mut File, e: &[f64]) { + let mut mean_energy = 0.0; + let mut mean_error = 0.0; + let mut mean_corr = 0.0; + for i in 0..ENERGY_CONV_AVERAGE_SAMPLE { + mean_energy += e[NOPTITER*3 - 3*i - 3]; + mean_error += e[NOPTITER*3 - 3*i - 2]; + mean_corr += e[NOPTITER*3 - 3*i - 1]; + } + mean_energy *= 1.0 / ENERGY_CONV_AVERAGE_SAMPLE as f64; + mean_error *= 1.0 / ENERGY_CONV_AVERAGE_SAMPLE as f64; + mean_corr *= 1.0 / ENERGY_CONV_AVERAGE_SAMPLE as f64; + + writeln!(fp, "{}", format!("{} {} {}", mean_energy, mean_error, mean_corr)).unwrap(); +} + +fn log_energy_convs(e: &[f64], fp: &mut File) { + for i in 0..NOPTITER { + writeln!(fp, "{}", format!("{} {} {}", e[3*i], e[3*i +1], e[3*i+2])).unwrap(); + } +} + +fn main() { + let mut fp = File::create("u_t_sweep").unwrap(); + writeln!(fp, "{}", format!("# {} {} {}", SIZE, NMCSAMP, NOPTITER)).unwrap(); + let mut paramsfp = File::create("params").unwrap(); + writeln!(paramsfp, "{}", format!("# {} {} {}", SIZE, NMCSAMP, NOPTITER)).unwrap(); + let mut _save: bool = true; + // Initialize logger + env_logger::init(); + let bitmask = generate_bitmask(&HOPPINGS, SIZE); + let mut rng = Mt64::new(SEED); + + for nsweep_iter in 0..NRATIO_POINTS { + let system_params = SysParams { + size: SIZE, + nelec: NELEC, + array_size: (SIZE + 7) / 8, + cons_t: -CONS_T, + cons_u: CONS_T * INITIAL_RATIO_UT + (nsweep_iter as f64 / NRATIO_POINTS as f64) * CONS_T * (FINAL_RATIO_UT - INITIAL_RATIO_UT), + nfij: NFIJ, + nvij: NVIJ, + ngi: NGI, + transfert_matrix: &HOPPINGS, + hopping_bitmask: &bitmask, + clean_update_frequency: CLEAN_UPDATE_FREQUENCY, + nmcsample: NMCSAMP, + nbootstrap: NBOOTSTRAP, + nmcwarmup: NMCWARMUP, + mcsample_interval: MCSAMPLE_INTERVAL, + tolerance_sherman_morrison: TOLERENCE_SHERMAN_MORRISSON, + tolerance_singularity: TOLERENCE_SINGULARITY, + pair_wavefunction: true, + }; + println!("U = {}, T = {}", system_params.cons_u, system_params.cons_t); + + let mut all_params: Vec = Vec::with_capacity(NGI + NVIJ + NFIJ); + for _ in 0..(NGI + NVIJ + NFIJ) { + all_params.push(rng.gen()); + } + let (gi, params) = all_params.split_at_mut(NGI); + let (vij, fij) = params.split_at_mut(NVIJ); + let mut parameters = VarParams { + fij, + gi, + vij, + size: SIZE + }; + + let vmcparams = VMCParams { + dt: OPTIMISATION_TIME_STEP, + optimisation_decay: OPTIMISATION_DECAY, + threshold: PARAM_THRESHOLD, + kmax: KMAX, + epsilon: EPSILON_SHIFT, + epsilon_cg: EPSILON_CG, + noptiter: NOPTITER, + nparams: NGI + NVIJ + NFIJ, + optimise: OPTIMISE, + optimise_gutzwiller: OPTIMISE_GUTZ, + optimise_jastrow: OPTIMISE_JAST, + optimise_orbital: OPTIMISE_ORB, + compute_energy_method: COMPUTE_ENERGY_METHOD, + optimise_energy_method: OPTIMISE_ENERGY_METHOD, + }; + + let state: FockState = { + let mut tmp: FockState = FockState::generate_from_nelec(&mut rng, NELEC, SIZE); + while tmp.spin_up.count_ones() != tmp.spin_down.count_ones() { + tmp = FockState::generate_from_nelec(&mut rng, NELEC, SIZE); + } + tmp + }; + + let e_array = variationnal_monte_carlo(&mut rng, state, &mut parameters, &system_params, &vmcparams); + write_energy(&mut fp, &e_array); + + log_energy_convs(&e_array, &mut paramsfp); + + } +} diff --git a/src/constants.rs b/src/constants.rs index 08256ad..a3eb0cc 100644 --- a/src/constants.rs +++ b/src/constants.rs @@ -1,4 +1,4 @@ -const SIZE: usize = 4; +const SIZE: usize = 9; const ARRAY_SIZE: usize = (SIZE + 7) / 8; /// System parameters /// TODOC diff --git a/src/dvmc.rs b/src/dvmc.rs new file mode 100644 index 0000000..c3bbf90 --- /dev/null +++ b/src/dvmc.rs @@ -0,0 +1,243 @@ +use log::{info, error}; +use rand::distributions::{Distribution, Standard}; +use rand::Rng; +use std::fmt::{Display, Debug}; +use blas::{idamax, daxpy, dcopy, dscal}; + +use crate::{VarParams, DerivativeOperator, SysParams, FockState, BitOps}; +use crate::monte_carlo::{compute_mean_energy, compute_mean_energy_exact}; +use crate::optimisation::{conjugate_gradiant, exact_overlap_inverse}; + +#[derive(Debug)] +pub enum EnergyOptimisationMethod { + ExactInverse, + ConjugateGradiant, +} + +#[derive(Debug)] +pub enum EnergyComputationMethod { + MonteCarlo, + ExactSum, +} + +#[derive(Debug)] +pub struct VMCParams { + pub noptiter: usize, + pub compute_energy_method: EnergyComputationMethod, + pub optimise_energy_method: EnergyOptimisationMethod, + pub epsilon: f64, + pub epsilon_cg: f64, + pub kmax: usize, + pub dt: f64, + pub optimisation_decay: f64, + pub nparams: usize, + pub threshold: f64, + pub optimise: bool, + pub optimise_gutzwiller: bool, + pub optimise_jastrow: bool, + pub optimise_orbital: bool, +} + +fn zero_out_derivatives(der: &mut DerivativeOperator, sys: &SysParams) { + for i in 0.. (sys.nfij + sys.nvij + sys.ngi) * sys.nmcsample { + der.o_tilde[i] = 0.0; + } + for i in 0..sys.nfij + sys.nvij + sys.ngi { + der.expval_o[i] = 0.0; + der.ho[i] = 0.0; + } + for i in 0..sys.nmcsample { + der.visited[i] = 0; + } + der.mu = -1; +} + +pub fn variationnal_monte_carlo( + rng: &mut R, + initial_state: FockState, + params: &mut VarParams, + sys: &SysParams, + vmcparams: &VMCParams +) -> Vec +where T: BitOps + From + Display + Debug, Standard: Distribution +{ + let mut output_energy_array = vec![0.0; vmcparams.noptiter * 3]; + let mut otilde: Vec = vec![0.0; (4*sys.nfij + sys.nvij + sys.ngi) * (sys.nmcsample + 1)]; + //let mut work_otilde: Vec = vec![0.0; (sys.nfij + sys.nvij + sys.ngi) * (sys.nmcsample + 1)]; + let mut expvalo: Vec = vec![0.0; 4*sys.nfij + sys.nvij + sys.ngi]; + //let mut work_expvalo: Vec = vec![0.0; sys.nfij + sys.nvij + sys.ngi]; + let mut expval_ho: Vec = vec![0.0; 4*sys.nfij + sys.nvij + sys.ngi]; + //let mut work_expval_ho: Vec = vec![0.0; sys.nfij + sys.nvij + sys.ngi]; + let mut visited: Vec = vec![0; sys.nmcsample + 1]; + //let mut work_visited: Vec = vec![0; sys.nmcsample + 1]; + let mut der = DerivativeOperator { + o_tilde: &mut otilde, + expval_o: &mut expvalo, + ho: &mut expval_ho, + n: (4*sys.nfij + sys.nvij + sys.ngi) as i32, + nsamp: match vmcparams.compute_energy_method { + EnergyComputationMethod::ExactSum => 1.0, + EnergyComputationMethod::MonteCarlo => sys.nmcsample as f64, + }, + nsamp_int: sys.mcsample_interval, + mu: -1, + visited: &mut visited, + pfaff_off: sys.ngi + sys.nvij, + jas_off: sys.ngi, + epsilon: vmcparams.epsilon, + }; + //let mut work_derivative = DerivativeOperator { + // o_tilde: &mut work_otilde, + // expval_o: &mut work_expvalo, + // ho: &mut work_expval_ho, + // n: (sys.nfij + sys.nvij + sys.ngi) as i32, + // nsamp: match vmcparams.compute_energy_method { + // EnergyComputationMethod::ExactSum => 1.0, + // EnergyComputationMethod::MonteCarlo => sys.nmcsample as f64, + // }, + // nsamp_int: sys.mcsample_interval, + // mu: -1, + // visited: &mut work_visited, + // pfaff_off: sys.ngi + sys.nvij, + // jas_off: sys.ngi, + // epsilon: vmcparams.epsilon, + //}; + for opt_iter in 0..vmcparams.noptiter { + let (mean_energy, _accumulated_states, deltae, correlation_time) = { + match vmcparams.compute_energy_method { + EnergyComputationMethod::MonteCarlo => compute_mean_energy(rng, initial_state, params, sys, &mut der), + EnergyComputationMethod::ExactSum => { + (compute_mean_energy_exact(params, sys, &mut der), Vec::with_capacity(0), 0.0, 0.0) + }, + } + }; + // Save energy, error and correlation_time. + output_energy_array[opt_iter * 3] = mean_energy; + output_energy_array[opt_iter * 3 + 1] = deltae; + output_energy_array[opt_iter * 3 + 2] = correlation_time; + + //// Copy out the relevant terms. + //work_derivative.mu = derivative.mu; + //let mut i = 0; + //for elem in derivative.visited.iter() { + // work_derivative.visited[i] = *elem; + // i += 1; + //} + //mapto_pairwf(&derivative, &mut work_derivative, sys); + + //let mut x0 = vec![0.0; sys.nfij + sys.nvij + sys.ngi]; + //x0[(sys.ngi + sys.nvij)..(sys.ngi + sys.nvij + sys.nfij)].copy_from_slice(params.fij); + //x0[sys.ngi..(sys.ngi + sys.nvij)].copy_from_slice(params.vij); + //x0[0..sys.ngi].copy_from_slice(params.gi); + + // 68 misawa + let mut b: Vec = vec![0.0; der.n as usize]; + unsafe { + let incx = 1; + let incy = 1; + daxpy(der.n, -mean_energy, der.expval_o, incx, der.ho, incy); + dcopy(der.n, der.ho, incx, &mut b, incy); + } + + let mut _flag: bool = true; + let ignored_columns = match vmcparams.optimise_energy_method { + EnergyOptimisationMethod::ExactInverse => { + exact_overlap_inverse(&der, &mut b, vmcparams.epsilon, vmcparams.nparams as i32, vmcparams.threshold) + }, + EnergyOptimisationMethod::ConjugateGradiant => { + panic!("Not implemented."); + //conjugate_gradiant(&der, &mut b, &mut x0, vmcparams.epsilon, vmcparams.kmax, vmcparams.nparams as i32, vmcparams.threshold, vmcparams.epsilon_cg) + }, + }; + let mut delta_alpha = vec![0.0; vmcparams.nparams]; + let mut j: usize = 0; + for i in 0..vmcparams.nparams { + if ignored_columns[i] { + continue; + } + delta_alpha[i] = b[j]; + j += 1; + if !::is_finite(delta_alpha[i]) { + _flag = false; + } + } + if vmcparams.optimise { + unsafe { + let incx = 1; + let incy = 1; + let alpha = - vmcparams.dt * ::exp(- (opt_iter as f64) * vmcparams.optimisation_decay); + if vmcparams.optimise_gutzwiller { + daxpy(sys.ngi as i32, alpha, &delta_alpha, incx, &mut params.gi, incy); + } + if vmcparams.optimise_jastrow { + daxpy(sys.nvij as i32, alpha, &delta_alpha[sys.ngi..vmcparams.nparams], incx, &mut params.vij, incy); + } + if vmcparams.optimise_orbital { + daxpy(sys.nfij as i32, alpha, &delta_alpha[sys.ngi + sys.nvij..vmcparams.nparams], incx, &mut params.fij, incy); + } + } + info!("Correctly finished optimisation iteration {}", opt_iter); + //info!("Rescaling the params."); + //let scale: f64 = unsafe { + // let incx = 1; + // let incy = 1; + // ddot(derivative.n, derivative.expval_o, incx, params.gi, incy) + //}; + //info!("Scale by : {}", scale); + //let ratio = 1.0 / (scale + 1.0); + //unsafe { + // let incx = 1; + // dscal(vmcparams.nparams as i32, ratio, params.gi, incx) + //} + //info!("Scaled params by ratio = {}", ratio); + + // JastrowGutzwiller Shifting + let mut shift = 0.0; + for i in 0..sys.ngi { + shift += params.gi[i]; + } + for i in 0..sys.nvij { + shift += params.vij[i]; + } + shift = shift / (sys.ngi + sys.nvij) as f64; + for i in 0..sys.ngi { + params.gi[i] -= shift; + } + for i in 0..sys.nvij { + params.vij[i] -= shift; + } + } + // HARD CODE vij = vji + // Slater Rescaling + unsafe { + let incx = 1; + let max = params.fij[idamax(sys.nfij as i32, ¶ms.fij, incx) - 1]; + if ::abs(max) < 1e-16 { + error!("Pfaffian params are all close to 0.0. Rescaling might bring noise."); + panic!("Undefined behavior."); + } + info!("Max was: {}", max); + dscal(sys.nfij as i32, 1.0 / max, &mut params.fij, incx); + } + //unsafe { + // dcopy( + // sys.nfij as i32, + // ¶ms.fij, + // 1, + // &mut params.fij[sys.nfij..2*sys.nfij], + // 1 + // ); + //} + zero_out_derivatives(&mut der, sys); + //print_delta_alpha(&delta_alpha, sys.ngi, sys.nvij, sys.nfij); + //let opt_delta = unsafe { + // let incx = 1; + // dnrm2(der.n, &delta_alpha, incx) + //}; + //error!("Changed params by norm {}", opt_delta); + //opt_progress_bar.inc(1); + //opt_progress_bar.set_message(format!("Changed params by norm: {:+>.05e} Current energy: {:+>.05e}", opt_delta, mean_energy)); + } + //opt_progress_bar.finish() + output_energy_array +} diff --git a/src/fock_state.rs b/src/fock_state.rs index 4200ff9..3602c1d 100644 --- a/src/fock_state.rs +++ b/src/fock_state.rs @@ -121,7 +121,8 @@ impl BitOps for I impl From for u8 { fn from(other: SpinState) -> Self { - if ARRAY_SIZE > 1 {error!("Casting SpinState of size {} to discards overflowing bits.", SIZE)} + //Dumbo + //if ARRAY_SIZE > 1 {error!("Casting SpinState of size {} to discards overflowing bits.", SIZE)} other.state[0] } } diff --git a/src/lib.rs b/src/lib.rs index 949c4dc..a3620cf 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -3,13 +3,6 @@ use pyo3::prelude::*; use blas::{dcopy, dgemm}; -/// Input file parsing util. -/// # Subfiles -/// * __`orbitale.csv`__ - Variationnal parameters for the orbital. In csv -/// format, 3 column, for $f_{ij}$. First column is `i`, second is `j`. The third -/// column is for the parameters identifier. -pub mod parse; - mod strings; // Have the FockState struct at the root. @@ -327,6 +320,10 @@ pub mod hamiltonian; /// TODOC pub mod optimisation; +/// Computation of the ground-state +/// TODOC +pub mod dvmc; + #[cfg(feature = "python-interface")] #[pymodule] fn impurity(m: &Bound<'_, PyModule>) -> PyResult<()> { diff --git a/src/monte_carlo.rs b/src/monte_carlo.rs index c2b1260..4ffcc83 100644 --- a/src/monte_carlo.rs +++ b/src/monte_carlo.rs @@ -217,7 +217,7 @@ fn accumulate_expvals(energy: &mut f64, state_energy: f64, der: &mut DerivativeO #[inline(always)] fn normalize(energy: &mut f64, energy_bootstraped: &mut f64, expval_o: &mut [f64], ho: &mut [f64], nsample: f64, nbootstrap: f64, nparams: usize) { - *energy *= 1.0 / nsample; + *energy *= 1.0 / (nsample - 1.0); *energy_bootstraped *= 1.0 / nsample; *energy_bootstraped *= 1.0 / nbootstrap; for i in 0..nparams { @@ -323,6 +323,14 @@ where Standard: Distribution derivatives.mu = 0; // Compute the derivative for the first element in the markov chain compute_derivative_operator(state, &pstate, derivatives, sys); + // Accumulate the first state into the markov chain + accumulated_states.push(state); + derivatives.visited[derivatives.mu as usize] += 1; + let state_energy = compute_hamiltonian(state, &pstate, proj, params, sys); + + accumulate_expvals(&mut energy, state_energy, derivatives, 1.0); + energy_error_estimation(state_energy, &mut previous_energies, &mut energy_sums, &mut + energy_quad_sums, &mut n_values, 0, error_estimation_level); for mc_it in 0..(sys.nmcsample * sys.mcsample_interval) { let mut proj_copy = proj; trace!("Before proposition: ~O_[0, {}] = {}", derivatives.mu + 1, derivatives.o_tilde[(derivatives.n * (derivatives.mu + 1)) as usize]); From 0a0cc66c964ed595f109558c5404e00ff243e688 Mon Sep 17 00:00:00 2001 From: Dimitri Bonanni-Surprenant Date: Thu, 27 Mar 2025 10:05:46 -0400 Subject: [PATCH 2/4] (fix) Conjugate gradient problem Conjugate gradient implementation was broken because - The output params step was in x0 -> now in b - Epsilon on the diagonal was poorly managed -> now added in the Ap computation --- src/bin/dvmc_pairwf.rs | 11 +++++++++-- src/optimisation.rs | 21 ++++++++++++++------- 2 files changed, 23 insertions(+), 9 deletions(-) diff --git a/src/bin/dvmc_pairwf.rs b/src/bin/dvmc_pairwf.rs index 5e2371c..efba5d8 100644 --- a/src/bin/dvmc_pairwf.rs +++ b/src/bin/dvmc_pairwf.rs @@ -16,7 +16,7 @@ const NVIJ: usize = SIZE*(SIZE - 1) / 2; const NGI: usize = SIZE; const NPARAMS: usize = NFIJ + NGI + NVIJ; const NELEC: usize = SIZE; -const NMCSAMP: usize = 1000; +const NMCSAMP: usize = 100; const NBOOTSTRAP: usize = 1; const NMCWARMUP: usize = 1000; const MCSAMPLE_INTERVAL: usize = 1; @@ -40,7 +40,7 @@ const OPTIMISE_JAST: bool = true; const OPTIMISE_ORB: bool = true; const SET_EXPVALO_ZERO: bool = false; const COMPUTE_ENERGY_METHOD: EnergyComputationMethod = EnergyComputationMethod::MonteCarlo; -const OPTIMISE_ENERGY_METHOD: EnergyOptimisationMethod = EnergyOptimisationMethod::ExactInverse; +const OPTIMISE_ENERGY_METHOD: EnergyOptimisationMethod = EnergyOptimisationMethod::ConjugateGradiant; pub enum EnergyOptimisationMethod { ExactInverse, @@ -357,6 +357,13 @@ fn main() { //spread_eigenvalues(&mut derivative); //println!("x0 in = {:?}", x0); let mut _flag: bool = true; + let mut bcopy = b.clone(); + //let ignored_columns = exact_overlap_inverse(&work_derivative, &mut bcopy, EPSILON_SHIFT, NPARAMS as i32, PARAM_THRESHOLD); + //conjugate_gradiant(&work_derivative, &mut b, &mut x0, EPSILON_SHIFT, KMAX, NPARAMS as i32, PARAM_THRESHOLD, EPSILON_CG); + //println!("Exact b: {:?}", bcopy); + //println!("CG b: {:?}", b); + //panic!("Stop"); + // Tmp to look at difference in CG and exact let ignored_columns = match OPTIMISE_ENERGY_METHOD { EnergyOptimisationMethod::ExactInverse => { exact_overlap_inverse(&work_derivative, &mut b, EPSILON_SHIFT, NPARAMS as i32, PARAM_THRESHOLD) diff --git a/src/optimisation.rs b/src/optimisation.rs index ffc773c..0520a80 100644 --- a/src/optimisation.rs +++ b/src/optimisation.rs @@ -5,13 +5,13 @@ use colored::Colorize; use crate::DerivativeOperator; -fn gradient(x: &[f64], otilde: &[f64], visited: &[usize], expval_o: &[f64], b: &mut [f64], dim: i32, mu: i32, nsamp: f64) { +fn gradient(x: &[f64], otilde: &[f64], visited: &[usize], expval_o: &[f64], b: &mut [f64], epsilon: f64, dim: i32, mu: i32, nsamp: f64) { let alpha = -1.0; let incx = 1; let incy = 1; let mut work: Vec = vec![0.0; dim as usize]; // Compute Ax - compute_w(&mut work, otilde, visited, expval_o, x, dim, mu, nsamp); + compute_w(&mut work, otilde, visited, expval_o, x, epsilon, dim, mu, nsamp); unsafe { // Compute b - Ax daxpy(dim, alpha, &work, incx, b, incy); @@ -26,7 +26,7 @@ fn update_x(x: &mut [f64], pk: &[f64], alpha: f64, dim: i32) { }; } -fn compute_w(w: &mut [f64], otilde: &[f64], visited: &[usize], expval_o: &[f64], p: &[f64], dim: i32, mu: i32, nsamp: f64) { +fn compute_w(w: &mut [f64], otilde: &[f64], visited: &[usize], expval_o: &[f64], p: &[f64], epsilon: f64, dim: i32, mu: i32, nsamp: f64) { // Computes Ap let incx = 1; let incy = 1; @@ -38,6 +38,8 @@ fn compute_w(w: &mut [f64], otilde: &[f64], visited: &[usize], expval_o: &[f64], unsafe { trace!("x_[n] = {:?}", p); trace!("mu = {}, n = {}", mu, dim); + // Copy p in w, so we can add epsilon*p to the result + dcopy(dim, p, incx, w, incy); // 80 misawa dgemv(b"T"[0], dim, mu, alpha, otilde, dim, p, incx, beta, &mut work, incy); for i in 0..mu as usize { @@ -47,7 +49,7 @@ fn compute_w(w: &mut [f64], otilde: &[f64], visited: &[usize], expval_o: &[f64], trace!("Len(work) = {}, a.n = {}, a.mu = {}", work.len(), dim, mu); trace!("~O_[0, mu] = {:?}", otilde.iter().step_by(mu as usize).copied().collect::>()); // Sometimes segfaults - dgemv(b"N"[0], dim, mu, 1.0 / nsamp, otilde, dim, &work, incx, beta, w, incy); + dgemv(b"N"[0], dim, mu, 1.0 / nsamp, otilde, dim, &work, incx, epsilon, w, incy); trace!("O_[m, mu] O^[T]_[mu, n] x_[n] = {:?}", w); let alpha = ddot(dim, expval_o, incx, p, incy); // 81 misawa @@ -360,7 +362,9 @@ pub fn exact_overlap_inverse(a: &DerivativeOperator, b: &mut [f64], epsilon: f64 /// Computes the solution of $A\mathbf{x}-\mathbf{b}=0$ /// TODOC -pub fn conjugate_gradiant(a: &DerivativeOperator, b: &mut [f64], x0: &mut [f64], _epsilon: f64, kmax: usize, dim: i32, thresh: f64, epsilon_convergence: f64) -> Vec +/// Output +/// b is the optimised parameter step. +pub fn conjugate_gradiant(a: &DerivativeOperator, b: &mut [f64], x0: &mut [f64], epsilon: f64, kmax: usize, dim: i32, thresh: f64, epsilon_convergence: f64) -> Vec { // PRE FILTER let mut ignore = vec![false; dim as usize]; @@ -408,7 +412,7 @@ pub fn conjugate_gradiant(a: &DerivativeOperator, b: &mut [f64], x0: &mut [f64], e *= epsilon_convergence; trace!("Error threshold e = {}", e); //println!("Error threshold e = {}", e); - gradient(x0, õ, a.visited, &expvalo, b, new_dim as i32, a.mu, a.nsamp); + gradient(x0, õ, a.visited, &expvalo, b, epsilon, new_dim as i32, a.mu, a.nsamp); let mut p = vec![0.0; new_dim]; let mut j: usize = 0; for i in 0..dim as usize { @@ -428,7 +432,7 @@ pub fn conjugate_gradiant(a: &DerivativeOperator, b: &mut [f64], x0: &mut [f64], trace!("p_{} : {:?}", k, p); //println!("r_{} : {:?}", k, b); //println!("p_{} : {:?}", k, p); - compute_w(&mut w, õ, a.visited, &expvalo, &p, new_dim as i32, a.mu, a.nsamp); + compute_w(&mut w, õ, a.visited, &expvalo, &p, epsilon, new_dim as i32, a.mu, a.nsamp); //println!("w_{} : {:?}", k, w); let nrm2rk = alpha_k(b, &p, &w, &mut alpha, new_dim as i32); trace!("alpha_{} : {}", k, alpha); @@ -450,5 +454,8 @@ pub fn conjugate_gradiant(a: &DerivativeOperator, b: &mut [f64], x0: &mut [f64], trace!("beta_{} : {}", k, beta); update_p(b, &mut p, beta, new_dim as i32); } + unsafe { + dcopy(new_dim as i32, x0, 1, b, 1); + } ignore } From 83e66e243680caea48499ffb512e1a3866230d76 Mon Sep 17 00:00:00 2001 From: Dimitri Bonanni-Surprenant Date: Mon, 7 Apr 2025 09:53:37 -0400 Subject: [PATCH 3/4] (feat) Pair wave-function final. --- src/bin/dvmc_pairwf.rs | 4 +-- src/bin/dvmc_pairwf2.rs | 67 ++++++++++++++++++++++++----------------- src/constants.rs | 2 +- src/dvmc.rs | 11 +++---- src/hoppings.rs | 2 +- 5 files changed, 48 insertions(+), 38 deletions(-) diff --git a/src/bin/dvmc_pairwf.rs b/src/bin/dvmc_pairwf.rs index efba5d8..55df29f 100644 --- a/src/bin/dvmc_pairwf.rs +++ b/src/bin/dvmc_pairwf.rs @@ -16,7 +16,7 @@ const NVIJ: usize = SIZE*(SIZE - 1) / 2; const NGI: usize = SIZE; const NPARAMS: usize = NFIJ + NGI + NVIJ; const NELEC: usize = SIZE; -const NMCSAMP: usize = 100; +const NMCSAMP: usize = 500; const NBOOTSTRAP: usize = 1; const NMCWARMUP: usize = 1000; const MCSAMPLE_INTERVAL: usize = 1; @@ -32,7 +32,7 @@ const OPTIMISATION_TIME_STEP: f64 = 1e-2; const OPTIMISATION_DECAY: f64 = 0.0; const NOPTITER: usize = 1000; const KMAX: usize = NPARAMS; -const PARAM_THRESHOLD: f64 = ::MIN; +const PARAM_THRESHOLD: f64 = 1e-5; //const PARAM_THRESHOLD: f64 = 0.0; const OPTIMISE: bool = true; const OPTIMISE_GUTZ: bool = true; diff --git a/src/bin/dvmc_pairwf2.rs b/src/bin/dvmc_pairwf2.rs index d3e3196..426c69e 100644 --- a/src/bin/dvmc_pairwf2.rs +++ b/src/bin/dvmc_pairwf2.rs @@ -10,7 +10,7 @@ use impurity::{VarParams, SysParams, generate_bitmask, FockState, RandomStateGen use impurity::dvmc::{variationnal_monte_carlo, EnergyOptimisationMethod, EnergyComputationMethod, VMCParams}; const SEED: u64 = 14; -const SIZE: usize = 16; +const SIZE: usize = 8; const NFIJ: usize = 4*SIZE*SIZE; const NVIJ: usize = SIZE*(SIZE - 1) / 2; const NGI: usize = SIZE; @@ -19,7 +19,7 @@ const NELEC: usize = SIZE; const NMCSAMP: usize = 1_000; const NBOOTSTRAP: usize = 1; const NMCWARMUP: usize = 1000; -const MCSAMPLE_INTERVAL: usize = 4; +const MCSAMPLE_INTERVAL: usize = 1; const _NTHREADS: usize = 1; const CLEAN_UPDATE_FREQUENCY: usize = 32; const TOLERENCE_SHERMAN_MORRISSON: f64 = 1e-12; @@ -43,7 +43,7 @@ const OPTIMISE_JAST: bool = true; const OPTIMISE_ORB: bool = true; const SET_EXPVALO_ZERO: bool = false; const COMPUTE_ENERGY_METHOD: EnergyComputationMethod = EnergyComputationMethod::MonteCarlo; -const OPTIMISE_ENERGY_METHOD: EnergyOptimisationMethod = EnergyOptimisationMethod::ExactInverse; +const OPTIMISE_ENERGY_METHOD: EnergyOptimisationMethod = EnergyOptimisationMethod::ConjugateGradiant; const ENERGY_CONV_AVERAGE_SAMPLE: usize = 20; // 4 et 2 Sites @@ -55,38 +55,49 @@ const ENERGY_CONV_AVERAGE_SAMPLE: usize = 20; // 0.0, 1.0, 1.0, 0.0 //]; -// 8 Sites +// 7 Sites //pub const HOPPINGS: [f64; SIZE*SIZE] = [ -// 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, -// 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, -// 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, -// 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, -// 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, -// 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, -// 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, -// 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, +// 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, +// 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, +// 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, +// 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, +// 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, +// 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, +// 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, //]; // 8 Sites pub const HOPPINGS: [f64; SIZE*SIZE] = [ - 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, - 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, - 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, - 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, - 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, - 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, - 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, + 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, + 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, + 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, + 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, + 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, + 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, ]; +// 16 Sites +//pub const HOPPINGS: [f64; SIZE*SIZE] = [ +// 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, +// 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, +// 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, +// 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, +// 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, +// 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, +// 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, +// 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, +// 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, +//]; + fn sq(a: f64) -> f64 { ::abs(a) * ::abs(a) diff --git a/src/constants.rs b/src/constants.rs index a3eb0cc..97fece5 100644 --- a/src/constants.rs +++ b/src/constants.rs @@ -1,4 +1,4 @@ -const SIZE: usize = 9; +const SIZE: usize = 17; const ARRAY_SIZE: usize = (SIZE + 7) / 8; /// System parameters /// TODOC diff --git a/src/dvmc.rs b/src/dvmc.rs index c3bbf90..f28d43d 100644 --- a/src/dvmc.rs +++ b/src/dvmc.rs @@ -125,10 +125,10 @@ where T: BitOps + From + Display + Debug, Standard: Distribution //} //mapto_pairwf(&derivative, &mut work_derivative, sys); - //let mut x0 = vec![0.0; sys.nfij + sys.nvij + sys.ngi]; - //x0[(sys.ngi + sys.nvij)..(sys.ngi + sys.nvij + sys.nfij)].copy_from_slice(params.fij); - //x0[sys.ngi..(sys.ngi + sys.nvij)].copy_from_slice(params.vij); - //x0[0..sys.ngi].copy_from_slice(params.gi); + let mut x0 = vec![0.0; sys.nfij + sys.nvij + sys.ngi]; + x0[(sys.ngi + sys.nvij)..(sys.ngi + sys.nvij + sys.nfij)].copy_from_slice(params.fij); + x0[sys.ngi..(sys.ngi + sys.nvij)].copy_from_slice(params.vij); + x0[0..sys.ngi].copy_from_slice(params.gi); // 68 misawa let mut b: Vec = vec![0.0; der.n as usize]; @@ -145,8 +145,7 @@ where T: BitOps + From + Display + Debug, Standard: Distribution exact_overlap_inverse(&der, &mut b, vmcparams.epsilon, vmcparams.nparams as i32, vmcparams.threshold) }, EnergyOptimisationMethod::ConjugateGradiant => { - panic!("Not implemented."); - //conjugate_gradiant(&der, &mut b, &mut x0, vmcparams.epsilon, vmcparams.kmax, vmcparams.nparams as i32, vmcparams.threshold, vmcparams.epsilon_cg) + conjugate_gradiant(&der, &mut b, &mut x0, vmcparams.epsilon, vmcparams.kmax, vmcparams.nparams as i32, vmcparams.threshold, vmcparams.epsilon_cg) }, }; let mut delta_alpha = vec![0.0; vmcparams.nparams]; diff --git a/src/hoppings.rs b/src/hoppings.rs index 78d8012..656b27a 100644 --- a/src/hoppings.rs +++ b/src/hoppings.rs @@ -10,7 +10,7 @@ pub fn generate_bitmask(transfer_matrix: &[f64], size: usize) -> Vec let one: u8 = 1; let mut j: usize = 0; while j < (size - 1 - i) { - if transfer_matrix[j + i + 1 + SIZE * j] != 0.0 { + if transfer_matrix[j + i + 1 + size * j] != 0.0 { if i == 0 { mask.state[(j + 2) / WORD_SIZE] ^= one << (WORD_SIZE - (j + 2) % WORD_SIZE); } else { From 935e739a289467e200e3f69c0f08b292785d4221 Mon Sep 17 00:00:00 2001 From: Dimitri Bonanni-Surprenant Date: Mon, 7 Apr 2025 09:54:53 -0400 Subject: [PATCH 4/4] (fix) Cagro changes --- Cargo.toml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Cargo.toml b/Cargo.toml index b5ca202..bfb75ff 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -51,6 +51,9 @@ harness = false [package.metadata.docs.rs] rustdoc-args = [ "--html-in-header", "katex-header.html" ] +[profile.release] +debug = true + [features] default = [] python-interface = ["dep:pyo3"]