From a18dd3bb992aba2703fe4fec9bdfc00b05690cb3 Mon Sep 17 00:00:00 2001 From: Dimitri Bonanni-Surprenant Date: Fri, 23 Aug 2024 14:55:03 -0400 Subject: [PATCH] feat: parameter optimisation --- src/bin/dvmc.rs | 149 -------------------- src/bin/dvmc_fastupdate.rs | 277 ++++++++++++++++--------------------- src/constants.rs | 25 ++++ src/density.rs | 79 ++++++++--- src/fock_state.rs | 12 +- src/gutzwiller.rs | 22 ++- src/hamiltonian.rs | 20 +-- src/hoppings.rs | 31 +++++ src/jastrow.rs | 30 +++- src/lib.rs | 34 ++++- src/monte_carlo.rs | 207 +++++++++++++++++++++++++++ src/optimisation.rs | 130 +++++++++++++++++ src/params.rs | 106 -------------- src/pfaffian.rs | 40 +++++- 14 files changed, 713 insertions(+), 449 deletions(-) delete mode 100644 src/bin/dvmc.rs create mode 100644 src/constants.rs create mode 100644 src/monte_carlo.rs create mode 100644 src/optimisation.rs delete mode 100644 src/params.rs diff --git a/src/bin/dvmc.rs b/src/bin/dvmc.rs deleted file mode 100644 index b6098ae..0000000 --- a/src/bin/dvmc.rs +++ /dev/null @@ -1,149 +0,0 @@ -use rand::Rng; -use std::ptr::addr_of; -use log::{debug, info, trace}; - -use impurity::{FockState, RandomStateGeneration, VarParams, Spin}; -use impurity::{FIJ, GI, VIJ, SIZE, CONS_U, CONS_T}; -use impurity::density::compute_internal_product; -use impurity::hamiltonian::potential; - -const NELEC: usize = 4; -const NMCSAMP: usize = 10_000; -const NMCWARMUP: usize = 100; - -fn propose_hopping(state: &FockState, rng: &mut R, params: &VarParams) -> (f64, FockState) { - let n0: usize = 0; - let n1: usize = 0; - let n2: Spin = Spin::Up; - let state2 = state.generate_hopping(rng, SIZE as u32, &mut (n0, n1, n2)); - let ip2 = compute_internal_product(state2, params); - (ip2, state2) -} - -//fn generate_random_params(rng: &mut R) -> VarParams { -// let mut fij: Vec = Vec::with_capacity(4*SIZE*SIZE); -// let mut vij: Vec = Vec::with_capacity(SIZE * SIZE); -// for _ in 0..SIZE * SIZE { -// vij.push(rng.gen()) -// } -// for _ in 0..4*SIZE * SIZE { -// fij.push(rng.gen()) -// } -// let mut gi: Vec = Vec::with_capacity(SIZE); -// for _ in 0..SIZE { -// gi.push(rng.gen()) -// } -// VarParams{fij, vij, gi} -//} - -fn compute_hamiltonian(state: FockState, _ip: f64, _params: &VarParams) -> f64 { - //let kin = kinetic(state, params); - //let e = (kin / ::exp(ip)) + potential(state); - let e = potential(state); - trace!("Inside compute hamiltonian Energy: {} for state: {}", e, state); - e -} - -fn log_system_parameters() { - info!("System parameter SIZE = {}", SIZE); - info!("System parameter NELEC = {}", NELEC); - info!("System parameter NMCSAMP = {}", NMCSAMP); - info!("System parameter NMCWARMUP = {}", NMCWARMUP); - info!("System parameter CONS_U = {}", CONS_U); - info!("System parameter CONS_T = {}", CONS_T); - for i in 0..4*SIZE*SIZE { - unsafe { - if FIJ[i] == 0.0 {continue;} - if i < SIZE*SIZE { - debug!("F_[{},{}]^[up, up]={}", i/SIZE, i%SIZE, FIJ[i]); - } - else if i < 2*SIZE*SIZE { - debug!("F_[{},{}]^[up, down]={}", i/SIZE - SIZE, i%SIZE, FIJ[i]); - } - else if i < 3*SIZE*SIZE { - debug!("F_[{},{}]^[down, up]={}", i/SIZE - 2*SIZE, i%SIZE, FIJ[i]); - } - else { - debug!("F_[{},{}]^[down, down]={}", i/SIZE - 3*SIZE, i%SIZE, FIJ[i]); - } - } - } - for i in 0..SIZE*SIZE { - unsafe { - if VIJ[i] == 0.0 {continue;} - debug!("V_[{},{}]={}", i/SIZE, i%SIZE, VIJ[i]); - } - } - for i in 0..SIZE { - unsafe { - if GI[i] == 0.0 {continue;} - debug!("G_[{}]={}", i, GI[i]); - } - } -} - -fn main() { - // Initialize logger - env_logger::init(); - log_system_parameters(); - - let mut rng = rand::thread_rng(); - //let parameters = generate_random_params(&mut rng); - let parameters = unsafe { VarParams { - fij: addr_of!(FIJ) as *const f64, - gi: addr_of!(GI) as *const f64, - vij: addr_of!(VIJ) as *const f64 - }}; - - let mut 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 - }; - - info!("Initial State: {}", state); - info!("Initial Nelec: {}, {}", state.spin_down.count_ones(), state.spin_up.count_ones()); - info!("Nsites: {}", state.n_sites); - - let mut lip = compute_internal_product(state, ¶meters); - let mut energy: f64 = 0.0; - - info!("Starting the warmup phase."); - // Warmup - for _ in 0..NMCWARMUP { - let (lip2, state2) = propose_hopping(&state, &mut rng, ¶meters); - trace!("Current state: {}", state); - trace!("Proposed state: {}", state2); - let ratio = ::exp(lip2 - lip); - trace!("Ratio: {}", ratio); - let w = rng.gen::(); - if ratio >= w { - // We ACCEPT - trace!("Accept."); - state = state2; - lip = lip2; - } - } - - info!("Starting the sampling phase."); - // MC Sampling - for _ in 0..NMCSAMP { - let (lip2, state2) = propose_hopping(&state, &mut rng, ¶meters); - trace!("Current state: {}", state); - trace!("Proposed state: {}", state2); - let ratio = ::exp(lip2 - lip); - trace!("Ratio: {}", ratio); - let w = rng.gen::(); - if ratio >= w { - // We ACCEPT - trace!("Accept."); - state = state2; - lip = lip2; - } - energy += compute_hamiltonian(state, lip, ¶meters); - } - energy = energy / NMCSAMP as f64; - info!("Final Energy: {}", energy); -} diff --git a/src/bin/dvmc_fastupdate.rs b/src/bin/dvmc_fastupdate.rs index 418c8e8..950dc9e 100644 --- a/src/bin/dvmc_fastupdate.rs +++ b/src/bin/dvmc_fastupdate.rs @@ -1,49 +1,83 @@ -use impurity::pfaffian::{update_pstate, PfaffianState}; -use rand::Rng; use std::ptr::addr_of; -use log::{debug, info, trace, warn}; +use std::slice::from_raw_parts as slice; +use blas::{daxpy, dcopy, dscal}; +use impurity::optimisation::conjugate_gradiant; +use log::{debug, info}; -use impurity::{FockState, RandomStateGeneration, VarParams, Spin}; -use impurity::{FIJ, GI, VIJ, SIZE, CONS_U, CONS_T}; -use impurity::density::{compute_internal_product_parts, fast_internal_product}; -use impurity::hamiltonian::{potential, kinetic}; +use impurity::{generate_bitmask, DerivativeOperator, FockState, RandomStateGeneration, SysParams, VarParams}; +use impurity::monte_carlo::compute_mean_energy; +const SIZE: usize = 4; +const NFIJ: usize = 4*SIZE*SIZE; +const NVIJ: usize = SIZE*SIZE; +const NGI: usize = SIZE; +const NPARAMS: usize = NFIJ + NGI + NVIJ; const NELEC: usize = SIZE; -const NMCSAMP: usize = 100_000; -const NMCWARMUP: usize = 10_000; -const CLEAN_UPDATE_FREQUENCY: usize = 0; +const NMCSAMP: usize = 1_000; +const NMCWARMUP: usize = 1_000; +const CLEAN_UPDATE_FREQUENCY: usize = 16; 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 EPSILON_CG: f64 = 1e-10; -fn propose_hopping( - state: &FockState, - pfaff_state: &PfaffianState, - previous_proj: &mut f64, - hop: &mut (usize, usize, Spin), - rng: &mut R, - params: &VarParams -) -> (f64, FockState, Vec, usize) { - let state2 = state.generate_hopping(rng, SIZE as u32, hop); - let (ratio_ip, updated_column, col_idx) = { - fast_internal_product(state, &state2, pfaff_state, &hop, previous_proj, params) - }; - (ratio_ip, state2, updated_column, col_idx) -} +pub const HOPPINGS: [f64; SIZE*SIZE] = [ + 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 +]; -fn compute_hamiltonian(state: FockState, pstate: &PfaffianState, proj: f64, params: &VarParams) -> f64 { - let kin = kinetic(state, pstate, proj, params); - let e = kin + potential(state); - trace!("Hamiltonian application = {} for state: |x> = {}", e, state); - e -} +pub static mut FIJ: [f64; NFIJ] = [ + 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, +1.078313550425773 , +0.007172274681240365, +0.028714778076311877, +0.09168843535310542 , +0.04813118562079141 , +1.0625398526882723 , +0.08433353658389342 , +0.002722470871706029, +0.07270002762085896 , +0.026989164590497917, +0.007555596176108393, +0.046284058565227465, +0.011127921360085048, +0.07287939415825727 , +0.08138828369394709 , +0.012799567556772274, + 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, +]; -fn log_system_parameters() { +pub static mut VIJ: [f64; NVIJ] = [ + 0.0, 2.0, 2.0, 2.0, + 2.0, 0.0, 2.0, 2.0, + 2.0, 2.0, 0.0, 2.0, + 2.0, 2.0, 2.0, 0.0, +]; + +pub static mut GI: [f64; NGI] = [ + -1.0, -1.0, -1.0, -1.0 +]; + +fn log_system_parameters(sys: &SysParams) { info!("System parameter SIZE = {}", SIZE); info!("System parameter NELEC = {}", NELEC); info!("System parameter NMCSAMP = {}", NMCSAMP); info!("System parameter NMCWARMUP = {}", NMCWARMUP); - info!("System parameter CONS_U = {}", CONS_U); - info!("System parameter CONS_T = {}", CONS_T); + info!("System parameter CONS_U = {}", sys.cons_u); + info!("System parameter CONS_T = {}", sys.cons_t); debug!("System parameter CLEAN_UPDATE_FREQUENCY = {}", CLEAN_UPDATE_FREQUENCY); debug!("System parameter TOLERENCE_SHERMAN_MORRISSON = {}", TOLERENCE_SHERMAN_MORRISSON); for i in 0..4*SIZE*SIZE { @@ -77,42 +111,38 @@ fn log_system_parameters() { } } -fn get_sign(s1: &FockState, hop: &(usize, usize, Spin)) -> usize { - let (a, b) = if hop.0 < hop.1 { - (hop.1, hop.0) - } else { - (hop.0, hop.1) - }; - trace!("First mask: {:08b}, second mask: {:08b}", !(::MAX >> a), ::MAX >> (b + 1)); - let mask = { - !(::MAX >> a) & (::MAX >> (b + 1)) - }; - let n_ones = match hop.2 { - Spin::Up => { - (s1.spin_up & mask).count_ones() - }, - Spin::Down => { - (s1.spin_down & mask).count_ones() - } - }; - - n_ones as usize -} - fn main() { // Initialize logger env_logger::init(); - log_system_parameters(); + let bitmask = generate_bitmask(&HOPPINGS, SIZE); + let system_params = SysParams { + size: SIZE, + nelec: NELEC, + array_size: (SIZE + 7) / 8, + cons_t: CONS_T, + cons_u: CONS_U, + nfij: NFIJ, + nvij: NVIJ, + ngi: NGI, + transfert_matrix: &HOPPINGS, + hopping_bitmask: &bitmask, + clean_update_frequency: CLEAN_UPDATE_FREQUENCY, + nmcsample: NMCSAMP, + nmcwarmup: NMCWARMUP, + tolerance_sherman_morrison: TOLERENCE_SHERMAN_MORRISSON, + tolerance_singularity: TOLERENCE_SINGULARITY + }; + log_system_parameters(&system_params); let mut rng = rand::thread_rng(); //let parameters = generate_random_params(&mut rng); let parameters = unsafe { VarParams { - fij: addr_of!(FIJ) as *const f64, - gi: addr_of!(GI) as *const f64, - vij: addr_of!(VIJ) as *const f64 + fij: slice(addr_of!(FIJ) as *const f64, NFIJ), + gi: slice(addr_of!(GI) as *const f64, NGI), + vij: slice(addr_of!(VIJ) as *const f64, NVIJ) }}; - let mut state: FockState = { + 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); @@ -120,106 +150,43 @@ fn main() { tmp }; + let mut otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let mut expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let mut expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let mut visited: Vec = vec![0; NMCSAMP]; + let mut derivative = DerivativeOperator { + o_tilde: &mut otilde, + expval_o: &mut expvalo, + ho: &mut expval_ho, + n: (NFIJ + NVIJ + NGI) as i32, + nsamp: NMCSAMP as f64, + mu: -1, + visited: &mut visited, + pfaff_off: NGI + NVIJ, + jas_off: NGI, + }; + info!("Initial State: {}", state); info!("Initial Nelec: {}, {}", state.spin_down.count_ones(), state.spin_up.count_ones()); info!("Nsites: {}", state.n_sites); - let (mut pstate, mut proj) = compute_internal_product_parts(state, ¶meters); - let mut hop: (usize, usize, Spin) = (0, 0, Spin::Up); - let mut _lip = ::ln(::abs(::exp(proj) * pstate.pfaff)) * 2.0; - let mut n_accepted_updates: usize = 0; - let mut energy: f64 = 0.0; - let mut proj_copy_persistent = proj; - let mut sign: usize = 0; - - info!("Starting the warmup phase."); - // Warmup - for _ in 0..NMCWARMUP { - let mut proj_copy = proj; - let (ratio, state2, col, colidx) = propose_hopping(&state, &pstate, &mut proj_copy, &mut hop, &mut rng, ¶meters); - trace!("Current state: {}", state); - trace!("Proposed state: {}", state2); - trace!("Ratio: {}", ratio); - let w = rng.gen::(); - if ::abs(ratio) * ::abs(ratio) >= w { - // We ACCEPT - trace!("Accept."); - n_accepted_updates += 1; + let mean_energy = compute_mean_energy(&mut rng, state, ¶meters, &system_params, &mut derivative); - // Clean update once in a while - if n_accepted_updates < CLEAN_UPDATE_FREQUENCY { - state = state2; - proj = proj_copy; - update_pstate(&mut pstate, hop, col, colidx); - } else { - let tmp_pfaff = pstate.pfaff; - trace!("PfaffianState before clean update: {:?}", pstate); - sign = (sign + get_sign(&state, &hop)) % 2; - trace!("Number of electrons between hopping: {}", get_sign(&state, &hop)); - (pstate, proj) = compute_internal_product_parts(state2, ¶meters); - let inverse_proj = ::exp(proj_copy_persistent - proj); - if sign != 0 { - pstate.pfaff *= -1.0; - } - trace!("Number of electrons between hopping: {}", get_sign(&state, &hop)); - (pstate, proj) = compute_internal_product_parts(state2, ¶meters); - if pstate.pfaff*pstate.pfaff < TOLERENCE_SINGULARITY { - warn!("Updated matrix is probably singular, got pfaffian {:.2e} and Tolerence is : {:e}.", pstate.pfaff, TOLERENCE_SINGULARITY); - } - trace!("PfaffianState after clean update: {:?}", pstate); - let err = ::abs(tmp_pfaff * ratio * inverse_proj) - ::abs(pstate.pfaff); - if err >= TOLERENCE_SHERMAN_MORRISSON { - warn!("Sherman-Morrisson update error of {:.2e} on computed pfaffian. Tolerence is : {:e}. Ratio was {}, states were: {} -> {}", err, TOLERENCE_SHERMAN_MORRISSON, ratio, state, state2); - } - n_accepted_updates = 0; - state = state2; - proj_copy_persistent = proj; - } - } - } + let mut x0 = vec![0.0; NFIJ + NVIJ + NGI]; + x0[(NGI + NVIJ)..(NGI + NVIJ + NFIJ)].copy_from_slice(parameters.fij); + x0[NGI..(NGI + NVIJ)].copy_from_slice(parameters.vij); + x0[0..NGI].copy_from_slice(parameters.gi); - info!("Starting the sampling phase."); - // MC Sampling - for _ in 0..NMCSAMP { - let mut proj_copy = proj; - let (ratio, state2, col, colidx) = propose_hopping(&state, &pstate, &mut proj_copy, &mut hop, &mut rng, ¶meters); - trace!("Current state: {}", state); - trace!("Proposed state: {}", state2); - trace!("Ratio: {}", ratio); - let w = rng.gen::(); - if ::abs(ratio) * ::abs(ratio) >= w { - // We ACCEPT - trace!("Accept."); - n_accepted_updates += 1; - // Clean update once in a while - // written in this way for branch prediction. - if n_accepted_updates < CLEAN_UPDATE_FREQUENCY { - state = state2; - proj = proj_copy; - update_pstate(&mut pstate, hop, col, colidx); - } else { - let tmp_pfaff = pstate.pfaff; - sign = (sign + get_sign(&state, &hop)) % 2; - trace!("Number of electrons between hopping: {}", get_sign(&state, &hop)); - (pstate, proj) = compute_internal_product_parts(state2, ¶meters); - let inverse_proj = ::exp(proj_copy_persistent - proj); - if sign != 0 { - pstate.pfaff *= -1.0; - } - trace!("Number of electrons between hopping: {}", get_sign(&state, &hop)); - (pstate, proj) = compute_internal_product_parts(state2, ¶meters); - let err = ::abs(tmp_pfaff * ratio * inverse_proj) - ::abs(pstate.pfaff); - if err >= TOLERENCE_SHERMAN_MORRISSON { - warn!("Sherman-Morrisson update error of {:.2e} on computed pfaffian. Tolerence is : {:e}. Ratio was {}", err, TOLERENCE_SHERMAN_MORRISSON, ratio); - } - n_accepted_updates = 0; - state = state2; - proj_copy_persistent = proj; - } - } - energy += compute_hamiltonian(state, &pstate, proj, ¶meters); + let mut b: Vec = vec![0.0; derivative.n as usize]; + unsafe { + let incx = 1; + let incy = 1; + dscal(derivative.n, 1.0 / (NMCSAMP as f64), derivative.ho, incx); + dscal(derivative.n, 1.0 / (NMCSAMP as f64), derivative.expval_o, incx); + daxpy(derivative.n, -mean_energy, derivative.expval_o, incx, derivative.ho, incy); + dcopy(derivative.n, derivative.ho, incx, &mut b, incy); } - info!("Final Energy: {:.2}", energy); - energy = energy / NMCSAMP as f64; - info!("Final Energy normalized: {:.2}", energy); + conjugate_gradiant(&derivative, &mut b, &mut x0, EPSILON_CG, 4, NPARAMS as i32); + info!("Need to update parameters with: {:?}", x0); + info!("Visited: {:?}", derivative.visited); } diff --git a/src/constants.rs b/src/constants.rs new file mode 100644 index 0000000..4a116cc --- /dev/null +++ b/src/constants.rs @@ -0,0 +1,25 @@ +const SIZE: usize = 4; +const ARRAY_SIZE: usize = (SIZE + 7) / 8; +/// System parameters +/// TODOC +#[derive(Debug)] +pub struct SysParams<'a> { + pub size: usize, + pub nelec: usize, + pub array_size: usize, + pub cons_u: f64, + pub cons_t: f64, + pub nfij: usize, + pub nvij: usize, + pub ngi: usize, + pub transfert_matrix: &'a [f64], + pub hopping_bitmask: &'a [SpinState], + pub nmcwarmup: usize, + pub nmcsample: usize, + pub clean_update_frequency: usize, + pub tolerance_singularity: f64, + pub tolerance_sherman_morrison: f64, +} + +pub type BitStruct = u8; +pub const NBITS: usize = 8; diff --git a/src/density.rs b/src/density.rs index 433de40..06f4e63 100644 --- a/src/density.rs +++ b/src/density.rs @@ -2,25 +2,22 @@ use log::trace; #[cfg(feature = "python-interface")] use pyo3::{pyfunction, PyResult}; -use std::slice::from_raw_parts as slice; - use crate::jastrow::{compute_jastrow_exp, fast_update_jastrow}; use crate::gutzwiller::{compute_gutzwiller_exp, fast_update_gutzwiller}; use crate::pfaffian::{construct_matrix_a_from_state, get_pfaffian_ratio, PfaffianState}; -use crate::{FockState, VarParams, BitOps, Spin, SpinState}; -use crate::{NFIJ, NVIJ, NGI}; +use crate::{BitOps, FockState, Spin, SpinState, VarParams}; pub fn compute_internal_product( state: FockState, - params: &VarParams + params: &VarParams, ) -> f64 where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl { - let (mut pfaffian_state, jastrow_exp, gutz_exp) = unsafe { + let (mut pfaffian_state, jastrow_exp, gutz_exp) = { ( - construct_matrix_a_from_state(slice(params.fij, NFIJ), state), - compute_jastrow_exp(state, slice(params.vij, NVIJ), state.n_sites), - compute_gutzwiller_exp(state, slice(params.gi, NGI), state.n_sites) + construct_matrix_a_from_state(params.fij, state), + compute_jastrow_exp(state, params.vij, state.n_sites), + compute_gutzwiller_exp(state, params.gi, state.n_sites) ) }; let pfaffian = pfaffian_state.pfaff; @@ -33,28 +30,72 @@ where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl pub fn compute_internal_product_parts( state: FockState, - params: &VarParams + params: &VarParams, ) -> (PfaffianState, f64) where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl { - let (mut pfaffian_state, jastrow_exp, gutz_exp) = unsafe { + let (mut pfaffian_state, jastrow_exp, gutz_exp) = { ( - construct_matrix_a_from_state(slice(params.fij, NFIJ), state), - compute_jastrow_exp(state, slice(params.vij, NVIJ), state.n_sites), - compute_gutzwiller_exp(state, slice(params.gi, NGI), state.n_sites) + construct_matrix_a_from_state(params.fij, state), + compute_jastrow_exp(state, params.vij, state.n_sites), + compute_gutzwiller_exp(state, params.gi, state.n_sites) ) }; pfaffian_state.rebuild_matrix(); (pfaffian_state, jastrow_exp + gutz_exp) } +pub fn fast_internal_product_no_otilde( + previous_state: &FockState, + new_state: &FockState, + previous_pstate: &PfaffianState, + hopping: &(usize, usize, Spin), + previous_proj: &mut f64, + params: &VarParams, +) -> (f64, Vec, usize) +where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl +{ + // Rename things. + let previous_i = hopping.0; + let new_i = hopping.1; + let spin = hopping.2; + let n_sites = new_state.n_sites; + + // Take a copy of the projectors + let proj_cp = *previous_proj; + + let (pfaffian_ratio, b_vec, col) = { + let fij = params.fij; + let vij = params.vij; + let gi = params.gi; + fast_update_jastrow(previous_proj, vij, previous_state, new_state, n_sites, previous_i, new_i); + match spin { + Spin::Up => { + fast_update_gutzwiller(previous_proj, gi, &previous_state.spin_down, previous_i, new_i); + }, + Spin::Down => { + fast_update_gutzwiller(previous_proj, gi, &previous_state.spin_up, previous_i, new_i); + } + + } + get_pfaffian_ratio(previous_pstate, previous_i, new_i, spin, fij) + }; + + // Combine to get the internal product. + trace!("Fast Projector value: {}, for state: {}", previous_proj, new_state); + trace!("Fast Computed /: {}, |x'> = {}, |x> = {}", pfaffian_ratio, new_state, previous_state); + let true_ratio = pfaffian_ratio * ::exp(*previous_proj - proj_cp); + trace!("Fast Computed /: {}, |x'> = {}, |x> = {}", true_ratio, new_state, previous_state); + (true_ratio, b_vec, col) +} + pub fn fast_internal_product( previous_state: &FockState, new_state: &FockState, previous_pstate: &PfaffianState, hopping: &(usize, usize, Spin), previous_proj: &mut f64, - params: &VarParams + params: &VarParams, ) -> (f64, Vec, usize) where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl { @@ -67,10 +108,10 @@ where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::o // Take a copy of the projectors let proj_cp = *previous_proj; - let (pfaffian_ratio, b_vec, col) = unsafe { - let fij = slice(params.fij, NFIJ); - let vij = slice(params.vij, NVIJ); - let gi = slice(params.gi, NGI); + let (pfaffian_ratio, b_vec, col) = { + let fij = params.fij; + let vij = params.vij; + let gi = params.gi; fast_update_jastrow(previous_proj, vij, previous_state, new_state, n_sites, previous_i, new_i); match spin { Spin::Up => { diff --git a/src/fock_state.rs b/src/fock_state.rs index 071b8ea..7d444ef 100644 --- a/src/fock_state.rs +++ b/src/fock_state.rs @@ -432,11 +432,11 @@ fn rot_right(s: T, n: usize, size: usize) -> T { } pub trait Hopper { - fn generate_all_hoppings(self: &Self) -> Vec<(usize, usize, Spin)>; + fn generate_all_hoppings(self: &Self, bitmask: &[SpinState]) -> Vec<(usize, usize, Spin)>; } impl> Hopper for FockState { - fn generate_all_hoppings(self: &FockState) -> Vec<(usize, usize, Spin)> { + fn generate_all_hoppings(self: &FockState, bitmask: &[SpinState]) -> Vec<(usize, usize, Spin)> { // This is the linear periodic version of the function let sup = self.spin_up; let sdo = self.spin_down; @@ -446,7 +446,7 @@ impl> Hopper for FockState { for j in 1..SIZE/2 + 1 { // See note 2024-06-08, rotate right gives all the horizontal links. // We need to go to j up to SIZE / 2 - let bitmask = ::from(HOP_BITMASKS[j - 1]); + let bitmask = ::from(bitmask[j - 1]); let mut possible_hoppings_up = (rot_right(sup, j, self.n_sites) ^ sup) & bitmask; let mut possible_hoppings_do = (rot_right(sdo, j, self.n_sites) ^ sdo) & bitmask; @@ -502,7 +502,7 @@ impl Distribution> for Standard where Standard: Distribu pub trait RandomStateGeneration { fn generate_from_nelec(rng: &mut R, nelec: usize, max_size: usize) -> Self; - fn generate_hopping(self: &Self, rng: &mut R, max_size: u32, out_idx: &mut (usize, usize, Spin)) -> Self; + fn generate_hopping(self: &Self, rng: &mut R, max_size: u32, out_idx: &mut (usize, usize, Spin), sys: &SysParams) -> Self; } impl> RandomStateGeneration for FockState where Standard: Distribution { @@ -529,9 +529,9 @@ impl> RandomStateGeneration for FockState where Standard } - fn generate_hopping(self: &FockState, rng: &mut R, max_size: u32, out_idx: &mut (usize, usize, Spin)) -> FockState { + fn generate_hopping(self: &FockState, rng: &mut R, max_size: u32, out_idx: &mut (usize, usize, Spin), sys: &SysParams) -> FockState { // This is cheap, don't sweat it - let all_hops = self.generate_all_hoppings(); + let all_hops = self.generate_all_hoppings(&sys.hopping_bitmask); let rand_index = rng.gen_range(0..all_hops.len()); let hop = all_hops[rand_index]; diff --git a/src/gutzwiller.rs b/src/gutzwiller.rs index 53e8419..9015bf0 100644 --- a/src/gutzwiller.rs +++ b/src/gutzwiller.rs @@ -3,7 +3,7 @@ use log::trace; #[cfg(feature = "python-interface")] use pyo3::prelude::*; -use crate::{BitOps, FockState}; +use crate::{BitOps, DerivativeOperator, FockState}; /// Computes the gutzwiller factor for a single fock state. /// # Arguments @@ -50,6 +50,26 @@ where gutz_out } +pub fn compute_gutzwiller_der( + fock_state: FockState, + n_sites: usize, + der: &mut DerivativeOperator, +) +where + T: BitOps + From + std::ops::Shl + Debug, +{ + // sum_i g_i n_i up n_i down + let mut gutzwiller_sites = fock_state.spin_up & fock_state.spin_down; + gutzwiller_sites.mask_bits(n_sites); + let mut i = gutzwiller_sites.leading_zeros() as usize; + while i < n_sites { + der.o_tilde[i * (der.mu + 1) as usize] = 1.0; + trace!("Computed gutzwiller derivative O_[{}, {}] = {}", i, der.mu+1, 1); + gutzwiller_sites.set(i); + i = gutzwiller_sites.leading_zeros() as usize; + } +} + /// Computes the gutzwiller fast update. /// # Arguments /// * __`previous_gutz`__ - The previous Gutzwiller coefficient $\ln P_{\text{G}}$. diff --git a/src/hamiltonian.rs b/src/hamiltonian.rs index cdb6495..f07f726 100644 --- a/src/hamiltonian.rs +++ b/src/hamiltonian.rs @@ -1,9 +1,9 @@ -use log::{debug, trace}; +use log::trace; -use crate::density::fast_internal_product; +use crate::density::fast_internal_product_no_otilde; use crate::pfaffian::PfaffianState; -use crate::{BitOps, FockState, SpinState, CONS_U, Hopper, HOPPINGS, SIZE}; -use crate::{VarParams, CONS_T, Spin}; +use crate::{BitOps, FockState, Hopper, SpinState}; +use crate::{VarParams, Spin, SysParams}; /// Computes the potential term of the Hamiltonian. /// # Arguments @@ -17,11 +17,11 @@ use crate::{VarParams, CONS_T, Spin}; /// $$ /// H_U=U\sum_i n_{i\uparrow}n_{i\downarrow} /// $$ -pub fn potential(state: FockState) -> f64 +pub fn potential(state: FockState, sys: &SysParams) -> f64 where T: BitOps, { - let pot = ((state.spin_up & state.spin_down).count_ones() as f64) * CONS_U; + let pot = ((state.spin_up & state.spin_down).count_ones() as f64) * sys.cons_u; trace!("Output potential = {:.2} for state |x> = {}", pot, state); pot } @@ -38,11 +38,11 @@ where /// $$ /// H_T=-t\sum_{,\sigma}c^\dagger_{i\sigma}c_{j\sigma}+c^\dagger_{j\sigma}c_{i\sigma} /// $$ -pub fn kinetic(state: FockState, previous_pstate: &PfaffianState, previous_proj: f64, params: &VarParams) -> f64 +pub fn kinetic(state: FockState, previous_pstate: &PfaffianState, previous_proj: f64, params: &VarParams, sys: &SysParams) -> f64 where T: BitOps + From + std::fmt::Debug + std::fmt::Display { - let hops = state.generate_all_hoppings(); + let hops = state.generate_all_hoppings(&sys.hopping_bitmask); let mut kin = 0.0; for hop in hops.into_iter() { @@ -58,10 +58,10 @@ where } }; let mut proj = previous_proj; - let (ratio, _col, _colidx) = fast_internal_product(&state, &f_state, previous_pstate, &hop, &mut proj, params); + let (ratio, _col, _colidx) = fast_internal_product_no_otilde(&state, &f_state, previous_pstate, &hop, &mut proj, params); trace!("Projection state: |x'> = {}, z = {}", f_state, ratio); trace!("Adding kinetic term t_[i,j]: |x> = {}, |x'> = {}, hop = ({}, {}, {:?}) Computed / = {}", state, f_state, hop.0, hop.1, hop.2, ratio); - kin += ratio*CONS_T*HOPPINGS[hop.0 + hop.1*SIZE]; + kin += ratio*sys.cons_t*sys.transfert_matrix[hop.0 + hop.1*sys.size]; } trace!("Output kinetic = {:.2} for state |x> = {}", kin, state); diff --git a/src/hoppings.rs b/src/hoppings.rs index e69de29..54815b4 100644 --- a/src/hoppings.rs +++ b/src/hoppings.rs @@ -0,0 +1,31 @@ +pub fn generate_bitmask(transfer_matrix: &[f64], size: usize) -> Vec { + const WORD_SIZE: usize = 8; + let all_zeros = SpinState{state: [0x00; ARRAY_SIZE]}; + let mut hop_tmp: Vec = Vec::with_capacity(size / 2); + // Index for array + let mut i: usize = 0; + while i < size / 2 { + hop_tmp.push({ + let mut mask = all_zeros; + 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 i == 0 { + mask.state[(j + 2) / WORD_SIZE] ^= one << (WORD_SIZE - (j + 2) % WORD_SIZE); + } else { + mask.state[(j + 1) / WORD_SIZE] ^= one << (WORD_SIZE - (j + 1) % WORD_SIZE); + } + } + j += 1; + } + if i == 0 && transfer_matrix[SIZE - 1] != 0.0 { + mask.state[0] ^= one << (WORD_SIZE - 1); + } + // Index for HOPPINGS + mask + }); + i += 1; + } + hop_tmp +} diff --git a/src/jastrow.rs b/src/jastrow.rs index 11a87e3..94f1f7e 100644 --- a/src/jastrow.rs +++ b/src/jastrow.rs @@ -2,7 +2,7 @@ use log::trace; #[cfg(feature = "python-interface")] use pyo3::prelude::*; -use crate::{BitOps, FockState}; +use crate::{BitOps, DerivativeOperator, FockState}; /// Computes the Jastrow exponent for a single fock state. /// # Arguments @@ -59,6 +59,34 @@ where jastrow_out } +pub fn compute_jastrow_der( + fock_state: FockState, + der: &mut DerivativeOperator, + n_sites: usize, +) +where + T: BitOps + std::fmt::Display, +{ + let mut regular_nor = !(fock_state.spin_up ^ fock_state.spin_down); + regular_nor.mask_bits(n_sites); + let mut i: usize = regular_nor.leading_zeros() as usize; + let mut indices: Vec = Vec::with_capacity(n_sites); + while i < n_sites { + indices.push(i); + regular_nor.set(i); + for nk in 0..indices.len() - 1 { + let (n1, n2) = (fock_state.spin_up, fock_state.spin_down); + let k = indices[nk]; + if n1.check(i) ^ n2.check(k) { + der.o_tilde[der.jas_off + i + k * n_sites + (der.n * (der.mu + 1)) as usize] -= 0.5; + } else { + der.o_tilde[der.jas_off + i + k * n_sites + (der.n * (der.mu + 1)) as usize] += 0.5; + } + } + i = regular_nor.leading_zeros() as usize; + } +} + fn jastrow_undo_update( spin_mask: &mut T, diff --git a/src/lib.rs b/src/lib.rs index 862201c..ce33a48 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -12,7 +12,33 @@ pub mod parse; include!("fock_state.rs"); // Include the params definition -include!("params.rs"); +include!("constants.rs"); +include!("hoppings.rs"); + + +/// Collection of the variationnal parameters +/// TODOC +pub struct VarParams<'a> { + pub fij: &'a [f64], + pub vij: &'a [f64], + pub gi: &'a [f64], +} + +/// The operator $O_k$ +/// TODOC +pub struct DerivativeOperator<'a> { + pub o_tilde: &'a mut [f64], + pub expval_o: &'a mut [f64], + pub ho: &'a mut [f64], + /// Number of variationnal parameters + pub n: i32, + /// Number of monte-carlo sampling + pub mu: i32, + pub nsamp: f64, + pub visited: &'a mut [usize], + pub pfaff_off: usize, + pub jas_off: usize, +} /// Module to calculate pfaffian /// # Usage @@ -57,6 +83,8 @@ include!("params.rs"); /// $$ pub mod pfaffian; +pub mod monte_carlo; + pub mod density; /// Calculate Jastrow coefficients @@ -123,6 +151,10 @@ pub mod gutzwiller; /// $$ pub mod hamiltonian; +/// The optimisation of the ground state +/// TODOC +pub mod optimisation; + #[cfg(feature = "python-interface")] #[pymodule] fn impurity(m: &Bound<'_, PyModule>) -> PyResult<()> { diff --git a/src/monte_carlo.rs b/src/monte_carlo.rs new file mode 100644 index 0000000..d94cddd --- /dev/null +++ b/src/monte_carlo.rs @@ -0,0 +1,207 @@ +use blas::daxpy; +use log::{info, trace, warn}; +use rand::distributions::{Distribution, Standard}; +use rand::Rng; + +use crate::gutzwiller::compute_gutzwiller_der; +use crate::jastrow::compute_jastrow_der; +use crate::{BitOps, DerivativeOperator, FockState, RandomStateGeneration, Spin, SysParams, VarParams}; +use crate::density::{compute_internal_product_parts, fast_internal_product}; +use crate::pfaffian::{compute_pfaffian_derivative, update_pstate, PfaffianState}; +use crate::hamiltonian::{kinetic, potential}; + +fn propose_hopping +> +( + state: &FockState, + pfaff_state: &PfaffianState, + previous_proj: &mut f64, + hop: &mut (usize, usize, Spin), + rng: &mut R, + params: &VarParams, + sys: &SysParams, +) -> (f64, FockState, Vec, usize) + where Standard: Distribution +{ + let state2 = state.generate_hopping(rng, sys.size as u32, hop, sys); + let (ratio_ip, updated_column, col_idx) = { + fast_internal_product(state, &state2, pfaff_state, &hop, previous_proj, params) + }; + (ratio_ip, state2, updated_column, col_idx) +} + +fn compute_hamiltonian(state: FockState, pstate: &PfaffianState, proj: f64, params: &VarParams, sys: &SysParams) -> f64 { + let kin = kinetic(state, pstate, proj, params, sys); + let e = kin + potential(state, sys); + trace!("Hamiltonian application = {} for state: |x> = {}", e, state); + e +} + +fn get_sign(s1: &FockState, hop: &(usize, usize, Spin)) -> usize { + let (a, b) = if hop.0 < hop.1 { + (hop.1, hop.0) + } else { + (hop.0, hop.1) + }; + trace!("First mask: {:08b}, second mask: {:08b}", !(::MAX >> a), ::MAX >> (b + 1)); + let mask = { + !(::ones() >> a) & (::ones() >> (b + 1)) + }; + let n_ones = match hop.2 { + Spin::Up => { + (s1.spin_up & mask).count_ones() + }, + Spin::Down => { + (s1.spin_down & mask).count_ones() + } + }; + + n_ones as usize +} + +fn compute_derivative_operator> +(state: FockState, pstate: &PfaffianState, der: &mut DerivativeOperator, params: &VarParams, sys: &SysParams) +{ + compute_gutzwiller_der(state, sys.size, der); + compute_jastrow_der(state, der, sys.size); + compute_pfaffian_derivative(pstate, der, sys); +} + +pub fn compute_mean_energy +> +(rng: &mut R, initial_state: FockState, params: &VarParams, sys: &SysParams, derivatives: &mut DerivativeOperator) -> f64 +where Standard: Distribution +{ + let mut state = initial_state; + let (mut pstate, mut proj) = compute_internal_product_parts(state, params); + let mut hop: (usize, usize, Spin) = (0, 0, Spin::Up); + let mut _lip = ::ln(::abs(::exp(proj) * pstate.pfaff)) * 2.0; + let mut n_accepted_updates: usize = 0; + let mut energy: f64 = 0.0; + let mut proj_copy_persistent = proj; + let mut sign: usize = 0; + + info!("Starting the warmup phase."); + // Warmup + for _ in 0..sys.nmcwarmup { + let mut proj_copy = proj; + let (ratio, state2, col, colidx) = propose_hopping(&state, &pstate, &mut proj_copy, &mut hop, rng, params, sys); + trace!("Current state: {}", state); + trace!("Proposed state: {}", state2); + trace!("Ratio: {}", ratio); + let w = rng.gen::(); + if ::abs(ratio) * ::abs(ratio) >= w { + // We ACCEPT + trace!("Accept."); + n_accepted_updates += 1; + + // Clean update once in a while + if n_accepted_updates < sys.clean_update_frequency { + state = state2; + proj = proj_copy; + update_pstate(&mut pstate, hop, col, colidx); + } else { + let tmp_pfaff = pstate.pfaff; + trace!("PfaffianState before clean update: {:?}", pstate); + sign = (sign + get_sign(&state, &hop)) % 2; + trace!("Number of electrons between hopping: {}", get_sign(&state, &hop)); + (pstate, proj) = compute_internal_product_parts(state2, params); + let inverse_proj = ::exp(proj_copy_persistent - proj); + if sign != 0 { + pstate.pfaff *= -1.0; + } + if pstate.pfaff*pstate.pfaff < sys.tolerance_singularity { + warn!("Updated matrix is probably singular, got pfaffian {:.2e} and Tolerence is : {:e}.", pstate.pfaff, sys.tolerance_singularity); + } + trace!("PfaffianState after clean update: {:?}", pstate); + let err = ::abs(tmp_pfaff * ratio * inverse_proj) - ::abs(pstate.pfaff); + if err >= sys.tolerance_sherman_morrison { + warn!("Sherman-Morrisson update error of {:.2e} on computed pfaffian. Tolerence is : {:e}. Ratio was {}, states were: {} -> {}", err, sys.tolerance_sherman_morrison, ratio, state, state2); + } + n_accepted_updates = 0; + state = state2; + proj_copy_persistent = proj; + } + } + } + + info!("Starting the sampling phase."); + // MC Sampling + for i in 0..derivatives.n as usize { + derivatives.o_tilde[i] *= 0.0; + } + for mc_it in 0..sys.nmcsample { + let mut proj_copy = proj; + trace!("Before proposition: ~O_[0, {}] = {}", derivatives.mu + 1, derivatives.o_tilde[(derivatives.n * (derivatives.mu + 1)) as usize]); + let (ratio, state2, col, colidx) = propose_hopping(&state, &pstate, &mut proj_copy, &mut hop, rng, params, sys); + trace!("After proposition: ~O_[0, {}] = {}", derivatives.mu + 1, derivatives.o_tilde[(derivatives.n * (derivatives.mu + 1)) as usize]); + trace!("Current state: {}", state); + trace!("Proposed state: {}", state2); + trace!("Ratio: {}", ratio); + let w = rng.gen::(); + if ::abs(ratio) * ::abs(ratio) >= w { + // We ACCEPT + trace!("Accept."); + n_accepted_updates += 1; + // Clean update once in a while + // written in this way for branch prediction. + if n_accepted_updates < sys.clean_update_frequency { + derivatives.mu += 1; + state = state2; + proj = proj_copy; + update_pstate(&mut pstate, hop, col, colidx); + } else { + derivatives.mu += 1; + let tmp_pfaff = pstate.pfaff; + sign = (sign + get_sign(&state, &hop)) % 2; + trace!("Number of electrons between hopping: {}", get_sign(&state, &hop)); + (pstate, proj) = compute_internal_product_parts(state2, params); + let inverse_proj = ::exp(proj_copy_persistent - proj); + if sign != 0 { + pstate.pfaff *= -1.0; + } + let err = ::abs(tmp_pfaff * ratio * inverse_proj) - ::abs(pstate.pfaff); + if pstate.pfaff*pstate.pfaff < sys.tolerance_singularity { + warn!("Updated matrix is probably singular, got pfaffian {:.2e} and Tolerence is : {:e}.", pstate.pfaff, sys.tolerance_singularity); + } + trace!("PfaffianState after clean update: {:?}", pstate); + if err >= sys.tolerance_sherman_morrison { + warn!("Sherman-Morrisson update error of {:.2e} on computed pfaffian. Tolerence is : {:e}. Ratio was {}", err, sys.tolerance_sherman_morrison, ratio); + } + n_accepted_updates = 0; + state = state2; + proj_copy_persistent = proj; + } + // Compute the derivative operator + compute_derivative_operator(state, &pstate, derivatives, params, sys); + } + derivatives.visited[(derivatives.mu + 1) as usize] += 1; + // Accumulate energy + let state_energy = compute_hamiltonian(state, &pstate, proj, params, sys); + energy += state_energy; + // Accumulate + unsafe { + let incx = 1; + let incy = 1; + daxpy( + derivatives.n, + state_energy, + &derivatives.o_tilde[(derivatives.n * (derivatives.mu + 1)) as usize .. (derivatives.n * (derivatives.mu + 2)) as usize], + incx, + &mut derivatives.ho, + incy + ); + } + // Accumulate the derivative operators + for i in 0 .. derivatives.n as usize { + derivatives.expval_o[i] += derivatives.o_tilde[i + (derivatives.n * (derivatives.mu + 1)) as usize]; + } + + } + info!("Final Energy: {:.2}", energy); + energy = energy / sys.nmcsample as f64; + info!("Final Energy normalized: {:.2}", energy); + energy +} diff --git a/src/optimisation.rs b/src/optimisation.rs new file mode 100644 index 0000000..02d2b7a --- /dev/null +++ b/src/optimisation.rs @@ -0,0 +1,130 @@ +use blas::{daxpy, dcopy, ddot, dgemv, dscal}; +use log::trace; + +use crate::DerivativeOperator; + +fn gradient(x: &[f64], a: &DerivativeOperator, b: &mut [f64], dim: i32) { + 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, a, x); + unsafe { + // Compute b - Ax + daxpy(dim, alpha, &work, incx, b, incy); + } +} + +fn update_x(x: &mut [f64], pk: &[f64], alpha: f64, dim: i32) { + let incx = 1; + let incy = 1; + unsafe { + daxpy(dim, alpha, &pk, incx, x, incy); + }; +} + +fn compute_w(w: &mut [f64], a: &DerivativeOperator, p: &[f64]) { + // Computes Ap + let incx = 1; + let incy = 1; + let alpha = 1.0; + let beta = 0.0; + + // Temp work vector + let mut work: Vec = vec![0.0; a.mu as usize]; + unsafe { + trace!("x_[n] = {:?}", p); + trace!("mu = {}, n = {}", a.mu, a.n); + dgemv(b"T"[0], a.n, a.mu, alpha, a.o_tilde, a.n, p, incx, beta, &mut work, incy); + for i in 0..a.mu as usize { + work[i] *= a.visited[i] as f64; + } + trace!("O^[T]_[mu, n] x_[n] = {:?}", work); + trace!("Len(work) = {}, a.n = {}, a.mu = {}", work.len(), a.n, a.mu); + trace!("~O_[0, mu] = {:?}", a.o_tilde.iter().step_by(a.mu as usize).copied().collect::>()); + // Sometimes segfaults + dgemv(b"N"[0], a.n, a.mu, 1.0 / a.nsamp, a.o_tilde, a.n, &work, incx, beta, w, incy); + trace!("O_[m, mu] O^[T]_[mu, n] x_[n] = {:?}", w); + let alpha = ddot(a.n, a.expval_o, incx, p, incy); + daxpy(a.n, - alpha, a.expval_o, incx, w, incy); + } +} + +fn update_r(rk: &mut [f64], w: &[f64], alphak: f64, dim: i32) { + let incx = 1; + let incy = 1; + unsafe { + daxpy(dim, - alphak, w, incx, rk, incy) + } +} + +fn alpha_k(r_k: &[f64], p_k: &[f64], w: &[f64], alphak: &mut f64, dim: i32) -> f64 { + let incx = 1; + let incy = 1; + unsafe { + let rkrk = ddot(dim, r_k, incx, r_k, incy); + *alphak = rkrk / ddot(dim, p_k, incx, w, incy); + rkrk + } +} + +fn beta_k(r_k: &[f64], rkrk: f64, dim: i32) -> f64 { + let incx = 1; + let incy = 1; + unsafe { + ddot(dim, r_k, incx, r_k, incy) / rkrk + } +} + +fn update_p(r_k: &[f64], p_k: &mut [f64], beta: f64, dim: i32) { + let alpha = 1.0; + let incx = 1; + let incy = 1; + unsafe { + dscal(dim, beta, p_k, incx); + daxpy(dim, alpha, r_k, incx, p_k, incy); + } +} + +/// 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) { + trace!("Initial guess x_0: {:?}", x0); + trace!("Initial equation b: {:?}", b); + let mut w: Vec = vec![0.0; dim as usize]; + // Error threshold + let e = unsafe { + let incx = 1; + let incy = 1; + ddot(dim, b, incx, b, incy) * epsilon + }; + trace!("Error threshold e = {}", e); + gradient(x0, a, b, dim); + let mut p = Vec::with_capacity(dim as usize); + unsafe { + let incx = 1; + let incy = 1; + p.set_len(dim as usize); + dcopy(dim, b, incx, &mut p, incy) + } + let mut alpha = 0.0; + + for k in 0..kmax { + trace!("r_{} : {:?}", k, b); + trace!("p_{} : {:?}", k, p); + compute_w(&mut w, a, &p); + let nrm2rk = alpha_k(b, &p, &w, &mut alpha, dim); + trace!("alpha_{} : {}", k, alpha); + update_x(x0, &p, alpha, dim); + trace!("x_{} : {:?}", k+1, x0); + update_r(b, &w, alpha, dim); + let beta = beta_k(b, nrm2rk, dim); + if beta * nrm2rk < e { + trace!("Achieved convergence at {} iterations", k); + break; + } + trace!("beta_{} : {}", k, beta); + update_p(b, &mut p, beta, dim); + } +} diff --git a/src/params.rs b/src/params.rs deleted file mode 100644 index 5af33f4..0000000 --- a/src/params.rs +++ /dev/null @@ -1,106 +0,0 @@ -/// Size of the system. -pub const SIZE: usize = 4; -pub const ARRAY_SIZE: usize = (SIZE + 7) / 8; -pub type BitStruct = u8; -pub const NBITS: usize = 8; - -/// Hubbard's model $U$ parameter -pub static CONS_U: f64 = 1.0; -/// Hubbard's model $t$ parameter -pub static CONS_T: f64 = -1.0; - -/// Collection of the variationnal parameters -pub struct VarParams { - pub fij: *const f64, - pub vij: *const f64, - pub gi: *const f64, -} - -pub const NFIJ: usize = 4*SIZE*SIZE; -pub const NVIJ: usize = SIZE*SIZE; -pub const NGI: usize = SIZE; - -pub static mut FIJ: [f64; NFIJ] = [ - 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, - //1.0, 2.0, 3.0, 4.0, - //1.0, -0.5, 1.0, -0.5, - //10.0, -1.5, 2.0, -3.5, - //2.0, -3.5, 10.0, -4.5, -1.078313550425773 , -0.007172274681240365, -0.028714778076311877, -0.09168843535310542 , -0.04813118562079141 , -1.0625398526882723 , -0.08433353658389342 , -0.002722470871706029, -0.07270002762085896 , -0.026989164590497917, -0.007555596176108393, -0.046284058565227465, -0.011127921360085048, -0.07287939415825727 , -0.08138828369394709 , -0.012799567556772274, - 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, -]; - -pub static mut VIJ: [f64; NVIJ] = [ - 0.0, 2.0, 2.0, 2.0, - 2.0, 0.0, 2.0, 2.0, - 2.0, 2.0, 0.0, 2.0, - 2.0, 2.0, 2.0, 0.0, -]; - -pub static mut GI: [f64; NGI] = [ - -1.0, -1.0, -1.0, -1.0 -]; - -pub const HOPPINGS: [f64; SIZE*SIZE] = [ - 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 -]; - -pub const HOP_BITMASKS: [SpinState; SIZE / 2] = { - const WORD_SIZE: usize = 8; - let all_zeros = SpinState{state: [0x00; ARRAY_SIZE]}; - let mut hop_tmp: [SpinState; SIZE / 2] = [all_zeros; SIZE / 2]; - // Index for array - let mut i: usize = 0; - while i < SIZE / 2 { - hop_tmp[i] = { - let mut mask = all_zeros; - let one: u8 = 1; - let mut j: usize = 0; - while j < (SIZE - 1 - i) { - if HOPPINGS[j + i + 1 + SIZE * j] != 0.0 { - if i == 0 { - mask.state[(j + 2) / WORD_SIZE] ^= one << (WORD_SIZE - (j + 2) % WORD_SIZE); - } else { - mask.state[(j + 1) / WORD_SIZE] ^= one << (WORD_SIZE - (j + 1) % WORD_SIZE); - } - } - j += 1; - } - if i == 0 && HOPPINGS[SIZE - 1] != 0.0 { - mask.state[0] ^= one << (WORD_SIZE - 1); - } - // Index for HOPPINGS - mask - }; - i += 1; - } - hop_tmp -}; diff --git a/src/pfaffian.rs b/src/pfaffian.rs index 065d481..d3b032e 100644 --- a/src/pfaffian.rs +++ b/src/pfaffian.rs @@ -1,4 +1,4 @@ -use crate::{BitOps, FockState, Spin}; +use crate::{BitOps, DerivativeOperator, FockState, Spin, SysParams}; use blas::{daxpy, ddot, dgemv, dger, dgemm}; use lapack::{dgetrf, dgetri}; use log::{error, trace}; @@ -187,6 +187,7 @@ where // Invert matrix. let pfaffian_value = compute_pfaffian_wq(&mut a.clone(), n as i32); invert_matrix(&mut a, n as i32); + trace!("Computed log abs pfaffian {} for state {}", ::ln(::abs(pfaffian_value)), state); trace!("Computed pfaffian {} for state {}", pfaffian_value, state); @@ -199,6 +200,43 @@ where } } +pub fn compute_pfaffian_derivative(pstate: &PfaffianState, der: &mut DerivativeOperator, sys: &SysParams) +{ + // Temporary bindings + let indices = &pstate.indices.0; + let indices2 = &pstate.indices.1; + let size = sys.size; + let a = &pstate.inv_matrix; + let n = pstate.n_elec; + let off = indices.len(); + + trace!("Packing derivatives from inverse matrix."); + for jj in 0..indices2.len() { + for ii in 0..indices.len() { + der.o_tilde[der.pfaff_off + indices2[jj] + size * indices[ii] + size*size + (der.n * (der.mu + 1)) as usize] += a[ii * n + (jj + off)]; + der.o_tilde[der.pfaff_off + indices[ii] + size * indices2[jj] + 2*size*size + (der.n * (der.mu + 1)) as usize] += a[(jj + off) * n + ii]; + } + } + for jj in 0..indices.len() { + for ii in 0..indices.len() { + if indices[ii] == indices[jj] { + continue; + } + der.o_tilde[der.pfaff_off + indices[ii] + size * indices[jj] + (der.n * (der.mu + 1)) as usize] += a[ii + jj * n]; + der.o_tilde[der.pfaff_off + indices[jj] + size * indices[ii] + (der.n * (der.mu + 1)) as usize] += a[jj + ii * n]; + } + } + for jj in 0..indices2.len() { + for ii in 0..indices2.len() { + if indices2[ii] == indices2[jj] { + continue; + } + der.o_tilde[der.pfaff_off + indices2[ii] + size * indices2[jj] + 3*size*size + (der.n * (der.mu + 1)) as usize] += a[ii + off + (jj + off) * n]; + der.o_tilde[der.pfaff_off + indices2[jj] + size * indices2[ii] + 3*size*size + (der.n * (der.mu + 1)) as usize] += a[jj + off + (ii + off) * n]; + } + } +} + /// Gets the ratio of pfaffian after an update. /// # Fields /// * __`previous_pstate`__ - The pfaffian state to update.