diff --git a/src/bin/dvmc_pairwf.rs b/src/bin/dvmc_pairwf.rs index 0f5bf6b..bc566e8 100644 --- a/src/bin/dvmc_pairwf.rs +++ b/src/bin/dvmc_pairwf.rs @@ -35,6 +35,7 @@ const NOPTITER: usize = 1000; const KMAX: usize = NPARAMS; const PARAM_THRESHOLD: f64 = 1e-5; //const PARAM_THRESHOLD: f64 = 0.0; +const FILTER_BEFORE_CUT: bool = true; const OPTIMISE: bool = true; const OPTIMISE_GUTZ: bool = true; const OPTIMISE_JAST: bool = true; @@ -361,7 +362,7 @@ fn main() { //spread_eigenvalues(&mut derivative); //println!("x0 in = {:?}", x0); let mut _flag: bool = true; - let mut bcopy = b.clone(); + //let 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); @@ -373,7 +374,7 @@ fn main() { exact_overlap_inverse(&work_derivative, &mut b, EPSILON_SHIFT, NPARAMS as i32, PARAM_THRESHOLD) }, EnergyOptimisationMethod::ConjugateGradiant => { - conjugate_gradiant(&work_derivative, &mut b, &mut x0, EPSILON_SHIFT, KMAX, NPARAMS as i32, PARAM_THRESHOLD, EPSILON_CG) + conjugate_gradiant(&work_derivative, &mut b, &mut x0, EPSILON_SHIFT, KMAX, NPARAMS as i32, PARAM_THRESHOLD, EPSILON_CG, FILTER_BEFORE_CUT) }, }; let mut delta_alpha = vec![0.0; NPARAMS]; diff --git a/src/bin/dvmc_pairwf2.rs b/src/bin/dvmc_pairwf2.rs index 7f4636a..0b95082 100644 --- a/src/bin/dvmc_pairwf2.rs +++ b/src/bin/dvmc_pairwf2.rs @@ -15,8 +15,8 @@ use impurity::dvmc::{variationnal_monte_carlo, EnergyOptimisationMethod, EnergyC type BitSize = u128; const SEED: u64 = 142; -const SIZE_N: usize = 6; -const SIZE_M: usize = 6; +const SIZE_N: usize = 4; +const SIZE_M: usize = 4; // SIZE = SIZE_N x SIZE_M const SIZE: usize = SIZE_N*SIZE_M; const NFIJ: usize = 4*SIZE*SIZE; @@ -24,12 +24,12 @@ 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 = 10000; const NBOOTSTRAP: usize = 1; -const NMCWARMUP: usize = 500; -const NWARMUPCHAINS: usize = 1; +const NMCWARMUP: usize = 1000; +const NWARMUPCHAINS: usize = 10; const MCSAMPLE_INTERVAL: usize = 1; -const NTHREADS: usize = 12; +const NTHREADS: usize = 1; const CLEAN_UPDATE_FREQUENCY: usize = 32; const TOLERENCE_SHERMAN_MORRISSON: f64 = 1e-12; const TOLERENCE_SINGULARITY: f64 = 1e-12; @@ -39,13 +39,16 @@ const INITIAL_RATIO_UT: f64 = 8.0; const FINAL_RATIO_UT: f64 = 32.0; const NRATIO_POINTS: usize = 1; const EPSILON_CG: f64 = 1e-16; -const EPSILON_SHIFT: f64 = 1e-3; +const EPSILON_SHIFT: f64 = 1e-4; const OPTIMISATION_TIME_STEP: f64 = 1e-2; const OPTIMISATION_DECAY: f64 = 0.0; -const NOPTITER: usize = 1000; +const NOPTITER: usize = 10000; const KMAX: usize = NPARAMS; +const FILTER_BEFORE_SHIFT: bool = true; // Better false (16 sites) +//const PARAM_THRESHOLD: f64 = ::EPSILON; const PARAM_THRESHOLD: f64 = 1e-3; //const PARAM_THRESHOLD: f64 = 0.0; +//const PARAM_THRESHOLD: f64 = -::INFINITY; const OPTIMISE: bool = true; const OPTIMISE_GUTZ: bool = true; const OPTIMISE_JAST: bool = true; @@ -58,12 +61,12 @@ const N_JAST: usize = NVIJ; const PAIRWF: bool = false; const CONV_PARAM_THRESHOLD: f64 = 1e-100; -const N_INDEP_PARAMS: usize = NFIJ + NGI + NVIJ; -//const N_INDEP_PARAMS: usize = SIZE*SIZE + NGI + NVIJ; +//const N_INDEP_PARAMS: usize = NFIJ + NGI + NVIJ; +const N_INDEP_PARAMS: usize = SIZE*SIZE + NGI + NVIJ; //const N_INDEP_PARAMS: usize = 3; -const SET_VIJ_ZERO: bool = false; -const SET_GI_EQUAL: bool = false; -const SET_PAIR_PFAFFIAN: bool = false; +const SET_VIJ_ZERO: bool = true; +const SET_GI_EQUAL: bool = true; +const SET_PAIR_PFAFFIAN: bool = true; pub const HOPPINGS: [f64; SIZE*SIZE] = { // Constructs hopping matrix for SITES_N*SITES_M @@ -118,47 +121,47 @@ pub const HOPPINGS: [f64; SIZE*SIZE] = { //]; // General rep -const PARAMS_PROJECTOR: [f64; (NFIJ + NVIJ + NGI) * (NFIJ + NVIJ + NGI - 1) / 2 + NFIJ + NVIJ + NGI] = { - let mut param = [0.0; (NFIJ + NVIJ + NGI) * (NFIJ + NVIJ + NGI - 1) / 2 + NFIJ + NVIJ + NGI]; - let mut j = 0; - let mut n = 0; - while j < NFIJ + NVIJ + NGI { - param[j + (j * (j+1) / 2)] = 1.0; - j += 1; - n += 1; - } - if n != N_INDEP_PARAMS { - panic!("Number of set independant params is not correct."); - } - param -}; - - -// General pairwf rep //const PARAMS_PROJECTOR: [f64; (NFIJ + NVIJ + NGI) * (NFIJ + NVIJ + NGI - 1) / 2 + NFIJ + NVIJ + NGI] = { // let mut param = [0.0; (NFIJ + NVIJ + NGI) * (NFIJ + NVIJ + NGI - 1) / 2 + NFIJ + NVIJ + NGI]; // let mut j = 0; // let mut n = 0; -// while j < NVIJ + NGI { -// param[j + (j * (j+1) / 2)] = 1.0; -// j += 1; -// n += 1; -// } -// let mut j = NVIJ + NGI; -// while j < SIZE * SIZE + NVIJ + NGI{ -// j += SIZE * SIZE; +// while j < NFIJ + NVIJ + NGI { // param[j + (j * (j+1) / 2)] = 1.0; -// j -= SIZE * SIZE; // j += 1; // n += 1; // } -// if n != N_INDEP_PARAMS { -// panic!("Number of set independant params is not correct."); -// } +// if n != N_INDEP_PARAMS { +// panic!("Number of set independant params is not correct."); +// } // param //}; +// General pairwf rep +const PARAMS_PROJECTOR: [f64; (NFIJ + NVIJ + NGI) * (NFIJ + NVIJ + NGI - 1) / 2 + NFIJ + NVIJ + NGI] = { + let mut param = [0.0; (NFIJ + NVIJ + NGI) * (NFIJ + NVIJ + NGI - 1) / 2 + NFIJ + NVIJ + NGI]; + let mut j = 0; + let mut n = 0; + while j < NVIJ + NGI { + param[j + (j * (j+1) / 2)] = 1.0; + j += 1; + n += 1; + } + let mut j = NVIJ + NGI; + while j < SIZE * SIZE + NVIJ + NGI{ + j += SIZE * SIZE; + param[j + (j * (j+1) / 2)] = 1.0; + j -= SIZE * SIZE; + j += 1; + n += 1; + } + if n != N_INDEP_PARAMS { + panic!("Number of set independant params is not correct."); + } + param +}; + + fn _sq(a: f64) -> f64 { ::abs(a) * ::abs(a) } @@ -242,6 +245,7 @@ fn main() { let mut all_params: Vec = Vec::with_capacity(NGI + NVIJ + NFIJ); for _ in 0..(NGI + NVIJ + NFIJ) { all_params.push(rngs[0].gen()); + //all_params.push(1.0); } let (gi, params) = all_params.split_at_mut(NGI); let (vij, fij) = params.split_at_mut(NVIJ); @@ -251,7 +255,8 @@ fn main() { vij, size: SIZE }; - let g = parameters.gi[0]; + //let g = parameters.gi[0]; + let g = 0.0; if SET_GI_EQUAL { for i in 0..NGI { parameters.gi[i] = g; @@ -270,6 +275,32 @@ fn main() { parameters.fij[i +3*SIZE*SIZE] = 0.0; } } + let mut max = ::MIN; + for i in 0..NGI { + if max < parameters.gi[i] { + max = parameters.gi[i]; + } + } + for i in 0..NVIJ { + if max < parameters.vij[i] { + max = parameters.vij[i]; + } + } + for i in 0..NFIJ { + if max < parameters.fij[i] { + max = parameters.fij[i]; + } + } + // Scale parameters + for i in 0..NGI { + parameters.gi[i] /= max; + } + for i in 0..NVIJ { + parameters.vij[i] /= max; + } + for i in 0..NFIJ { + parameters.fij[i] /= max; + } //println!("{:?}", parameters.fij); let vmcparams = VMCParams { @@ -288,7 +319,8 @@ fn main() { compute_energy_method: COMPUTE_ENERGY_METHOD, optimise_energy_method: OPTIMISE_ENERGY_METHOD, conv_param_threshold: CONV_PARAM_THRESHOLD, - nthreads: NTHREADS + nthreads: NTHREADS, + filter_before_shift: FILTER_BEFORE_SHIFT }; let state: FockState = { diff --git a/src/density.rs b/src/density.rs index 6b3e911..ed144ac 100644 --- a/src/density.rs +++ b/src/density.rs @@ -3,7 +3,7 @@ use log::trace; use pyo3::{pyfunction, PyResult}; use crate::jastrow::{compute_jastrow_exp, fast_update_jastrow}; -use crate::gutzwiller::{compute_gutzwiller_exp, fast_update_gutzwiller}; +use crate::gutzwiller::{compute_gutzwiller_exp, fast_update_gutzwiller, fast_update_gutzwiller_spin_change}; use crate::pfaffian::{construct_matrix_a_from_state, get_pfaffian_ratio, PfaffianState}; use crate::{BitOps, FockState, Spin, SpinState, VarParams, SysParams}; @@ -77,7 +77,39 @@ where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::o } } - get_pfaffian_ratio(previous_pstate, previous_i, new_i, spin, fij) + get_pfaffian_ratio(previous_pstate, previous_i, new_i, spin, 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); + (pfaffian_ratio, b_vec, col) +} + +pub fn fast_internal_product_spin_change( + previous_state: &FockState, + new_state: &FockState, + previous_pstate: &PfaffianState, + hopping: &(usize, usize, Spin, Spin), + previous_proj: &mut f64, + params: &VarParams, +) -> (f64, Vec, usize) +where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl + Send +{ + // Rename things. + let previous_i = hopping.0; + let new_i = hopping.1; + let previous_spin = hopping.2; + let new_spin = hopping.3; + let n_sites = new_state.n_sites; + + let (pfaffian_ratio, b_vec, col) = { + let fij = ¶ms.fij; + let vij = ¶ms.vij; + let gi = ¶ms.gi; + fast_update_jastrow(previous_proj, vij, previous_state, new_state, n_sites, previous_i, new_i); + fast_update_gutzwiller_spin_change(previous_proj, gi, previous_state, previous_i, new_i, previous_spin, new_spin); + get_pfaffian_ratio(previous_pstate, previous_i, new_i, previous_spin, new_spin, fij) }; // Combine to get the internal product. @@ -116,7 +148,7 @@ where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::o } } - get_pfaffian_ratio(previous_pstate, previous_i, new_i, spin, fij) + get_pfaffian_ratio(previous_pstate, previous_i, new_i, spin, spin, fij) }; // Combine to get the internal product. diff --git a/src/dvmc.rs b/src/dvmc.rs index 7046711..5b82c35 100644 --- a/src/dvmc.rs +++ b/src/dvmc.rs @@ -39,6 +39,7 @@ pub struct VMCParams { pub optimise_jastrow: bool, pub optimise_orbital: bool, pub conv_param_threshold: f64, + pub filter_before_shift: bool, } fn merge_derivatives(der: &mut DerivativeOperator, der_vec: &mut [DerivativeOperator], param_map: &GenParameterMap, nmc: usize, nthreads: usize, x0: &mut [f64]){ @@ -214,7 +215,7 @@ where T: BitOps + From + Display + Debug + Send + Sync, Standard: Distributi vmcparams.nparams as i32, -1, match vmcparams.compute_energy_method { - EnergyComputationMethod::ExactSum => 1.0, + EnergyComputationMethod::ExactSum => sys.nmcsample as f64, EnergyComputationMethod::MonteCarlo => (sys.nmcsample * vmcparams.nthreads) as f64, }, sys.mcsample_interval, @@ -227,20 +228,34 @@ where T: BitOps + From + Display + Debug + Send + Sync, Standard: Distributi let work_der = param_map.mapto_general_representation(&der, &mut x0); work_der_vec.push(work_der); } + //match vmcparams.compute_energy_method { + // EnergyComputationMethod::MonteCarlo => { + // }, + // EnergyComputationMethod::ExactSum => { + // der.nsamp = 1.0; + // }, + //} for opt_iter in 0..vmcparams.noptiter { + //println!("gi = {:?}", params.gi); + //println!("vij = {:?}", params.vij); + //println!("fij = {:?}", params.fij); sys._opt_iter = opt_iter; let (mean_energy, _accumulated_states, deltae, correlation_time) = { match vmcparams.compute_energy_method { EnergyComputationMethod::MonteCarlo => { - parallel_monte_carlo(rng, initial_state, &mut work_der_vec, ¶m_map, &mut der.o_tilde, vmcparams, params, sys) + let (mean_energy, _accumulated_states, deltae, correlation_time) = + parallel_monte_carlo(rng, initial_state, &mut work_der_vec, ¶m_map, &mut der.o_tilde, vmcparams, params, sys); + update_initial_state(initial_state, &_accumulated_states, sys.nmcsample, vmcparams.nthreads); + (mean_energy, _accumulated_states, deltae, correlation_time) }, EnergyComputationMethod::ExactSum => { - (compute_mean_energy_exact(params, sys, &mut work_der_vec[0]), Vec::with_capacity(0), 0.0, 0.0) + let (mean_energy, _accumulated_states, deltae, correlation_time) = + (compute_mean_energy_exact(params, sys, &mut work_der_vec[0]), Vec::with_capacity(0), 0.0, 0.0); + param_map.parallel_reduced_representation_otilde(&work_der_vec[0], &mut der.o_tilde); + (mean_energy, _accumulated_states, deltae, correlation_time) }, } }; - println!("Finished computing an energy."); - update_initial_state(initial_state, &_accumulated_states, sys.nmcsample, vmcparams.nthreads); // Save energy, error and correlation_time. output_energy_array[opt_iter * 3] = mean_energy; output_energy_array[opt_iter * 3 + 1] = deltae; @@ -259,6 +274,8 @@ where T: BitOps + From + Display + Debug + Send + Sync, Standard: Distributi x0[sys.ngi..(sys.ngi + sys.nvij)].copy_from_slice(params.vij); x0[0..sys.ngi].copy_from_slice(params.gi); merge_derivatives(&mut der, &mut work_der_vec, param_map, sys.nmcsample, vmcparams.nthreads, &mut x0); + //println!("{:?}", work_der_vec[0].o_tilde); + //println!("{:?}", der.o_tilde); let mut mu_vec: Vec = vec![0; vmcparams.nthreads]; for i in 0..vmcparams.nthreads { mu_vec[i] = work_der_vec[i].mu; @@ -277,6 +294,7 @@ where T: BitOps + From + Display + Debug + Send + Sync, Standard: Distributi daxpy(der.n, -mean_energy, &der.expval_o, incx, &mut der.ho, incy); dcopy(der.n, &der.ho, incx, &mut b, incy); } + //println!("gm = {:?}", b); let b_nrm = unsafe { let incx = 1; dnrm2(param_map.dim, &b, incx) @@ -289,39 +307,62 @@ where T: BitOps + From + Display + Debug + Send + Sync, Standard: Distributi 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) + exact_overlap_inverse(&der, &mut b, vmcparams.epsilon, vmcparams.nparams as i32, + vmcparams.threshold) }, EnergyOptimisationMethod::ConjugateGradiant => { - 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, + vmcparams.filter_before_shift) }, }; let mut delta_alpha = vec![0.0; param_map.gendim as usize]; let mut j: usize = 0; + let mut norm2 = 0.0; for i in 0..vmcparams.nparams { if ignored_columns[i] { continue; } delta_alpha[i] = b[j]; + norm2 += b[j] * b[j]; j += 1; if !::is_finite(delta_alpha[i]) { _flag = false; + error!("Parameter update contains NaN or Inf at iter {}.", opt_iter); + panic!("Undefined behavior"); } } - //panic!("Stop!"); + println!("norm2 = {}", norm2); + //if norm2 > 10000.0 { + // unsafe { + // dscal(vmcparams.nparams as i32, ::sqrt(1e-2 / norm2), &mut delta_alpha, 1); + // } + //} + //unsafe { + // dscal(vmcparams.nparams as i32, ::sqrt(1e0 / norm2), &mut delta_alpha, 1); + //} + //if opt_iter == 74 { + // println!("da = {:?}", delta_alpha); + // println!("norm2 = {}", norm2); + //} + param_map.update_delta_alpha_reduced_to_gen(&mut delta_alpha); if vmcparams.optimise { unsafe { let incx = 1; let incy = 1; let alpha = - vmcparams.dt * ::exp(- (opt_iter as f64) * vmcparams.optimisation_decay); + //let alpha = 1.0; 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); + 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); + 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); diff --git a/src/fock_state.rs b/src/fock_state.rs index c2df100..2f20b95 100644 --- a/src/fock_state.rs +++ b/src/fock_state.rs @@ -201,8 +201,10 @@ impl BitOps for SpinState { } fn mask_bits(&mut self, by: usize) { - let meta_by = by / 8; - let granular_by = by % 8; + let mask_by = SIZE - by; + let meta_by = mask_by / 8; + let granular_by = mask_by % 8; + println!("{}", granular_by); for i in 0..meta_by { self.state[ARRAY_SIZE - 1 - i] &= 0; } diff --git a/src/green.rs b/src/green.rs new file mode 100644 index 0000000..5c8f43e --- /dev/null +++ b/src/green.rs @@ -0,0 +1,478 @@ +use log::{error, info, trace, warn}; +use rand::distributions::{Distribution, Standard}; +use rand::Rng; + +use crate::{BitOps, FockState, RandomStateGeneration, Spin, SysParams, VarParams}; +use crate::density::{compute_internal_product_parts, fast_internal_product, fast_internal_product_spin_change}; +use crate::pfaffian::{update_pstate, PfaffianState}; +use crate::hamiltonian::{kinetic, potential}; + +pub enum Projector { + Identity, + ElectronExcitation { + i: usize, + sigmai: Spin, + }, + PairExcitation { + i: usize, + sigmai: Spin, + j: usize, + sigmaj: Spin + }, +} + +//fn propose_exchange +//> +//( +// state: &FockState, +// pfaff_state: &PfaffianState, +// previous_proj: &mut f64, +// exchange: &mut (usize, usize), +// rng: &mut R, +// params: &VarParams, +// sys: &SysParams, +//) -> (f64, FockState, Vec, usize) +// where Standard: Distribution +//{ +// let state2 = state.generate_exchange(rng, exchange); +// 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) +//} + +#[inline(always)] +fn propose_hopping + + Send> +( + 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); + 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) +} + +#[inline(always)] +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, proj, pstate, sys); + trace!("Hamiltonian application = {} for state: |x> = {}", e, state); + e / (pstate.pfaff * ::exp(proj)) +} + +#[inline(always)] +fn make_update + Send>( + n: &mut usize, + proj: &mut f64, + proj_copy: &mut f64, + proj_copy_persistent: &mut f64, + ratio: &f64, + ratio_prod: &mut f64, + state: &mut FockState, + state2: &FockState, + hop: &(usize, usize, Spin), + col: Vec, + colidx: usize, + pstate: &mut PfaffianState, + params: &VarParams, + sys: &SysParams +) { + // Clean update once in a while + // written in this way for branch prediction. + if *n < sys.clean_update_frequency { + *state = *state2; + *proj = *proj_copy; + update_pstate(pstate, *hop, col, colidx); + *ratio_prod *= ratio; + } else { + let tmp_pfaff = pstate.pfaff; + (*pstate, *proj) = compute_internal_product_parts(*state2, params, sys); + let inverse_proj = ::exp(*proj_copy_persistent - *proj); + let err = ::abs(::abs(tmp_pfaff * *ratio * *ratio_prod * 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}. Expected {}, got {}", err, sys.tolerance_sherman_morrison, tmp_pfaff * *ratio * *ratio_prod * inverse_proj, pstate.pfaff); + } + *n = 0; + *state = *state2; + *proj_copy_persistent = *proj; + *ratio_prod = 1.0; + } +} + +#[inline(always)] +fn warmup( + rng: &mut R, + state: &mut FockState, + hop: &mut (usize, usize, Spin), + proj: &mut f64, + proj_copy_persistent: &mut f64, + ratio_prod: &mut f64, + pstate: &mut PfaffianState, + n_accepted_updates: &mut usize, + params: &VarParams, + sys: &SysParams +) +where T: BitOps + From + std::fmt::Debug + std::fmt::Display + Send, + R: Rng + ?Sized, + Standard: Distribution +{ + // Warmup + for _ in 0..sys.nmcwarmup { + let mut proj_copy = *proj; + let (mut ratio, state2, col, colidx) = propose_hopping(&state, &pstate, &mut proj_copy, hop, rng, params, sys); + trace!("Current state: {}", state); + trace!("Proposed state: {}", state2); + trace!("Ratio: {}", ratio); + ratio *= ::exp(proj_copy - *proj); + let w = rng.gen::(); + if ::abs(ratio) * ::abs(ratio) >= w { + // We ACCEPT + trace!("Accept."); + *n_accepted_updates += 1; + make_update( + n_accepted_updates, + proj, + &mut proj_copy, + proj_copy_persistent, + &ratio, + ratio_prod, + state, + &state2, + &hop, + col, + colidx, + pstate, + params, + sys + ); + + } + } + +} + +#[inline(always)] +fn accumulate_expvals(energy: &mut f64, state_energy: f64, rho: f64) { + // Accumulate Energy + *energy += state_energy; +} + +#[inline(always)] +fn normalize(energy: &mut f64, energy_bootstraped: &mut f64, nsample: f64, nbootstrap: f64) { + *energy *= 1.0 / (nsample - 1.0); + *energy_bootstraped *= 1.0 / nsample; + *energy_bootstraped *= 1.0 / nbootstrap; +} + +#[inline(always)] +fn sq(x: f64) -> f64 +{ + x * x +} + +fn general_hopping + Send> +(state: &FockState, hop: &(usize, usize, Spin, Spin)) -> Option> +{ + let i = hop.0; + let j = hop.1; + let mut state2 = state.clone(); + match hop.3 { + Spin::Up => { + if state.spin_up.check(j) { + return None; + } + state2.spin_up.set(j); + }, + Spin::Down => { + if state.spin_down.check(j) { + return None; + } + state2.spin_down.set(j); + }, + } + match hop.2 { + Spin::Up => { + if !state.spin_up.check(i) { + return None; + } + state2.spin_up.set(i); + }, + Spin::Down => { + if !state.spin_down.check(i) { + return None; + } + state2.spin_down.set(i); + }, + } + Some(state2) +} + +// C^t_{i,sigmai} C_{j,sigmaj} C^t_{k,sigmak} C_{l,sigmal} +#[inline(always)] +fn opperate_by_correlator_twice + + Send> +( + state: &FockState, + pfaff_state: &PfaffianState, + previous_proj: &mut f64, + hop: &(usize, usize, Spin, Spin), + params: &VarParams, +) -> Option<(f64, FockState, Vec, usize)> +{ + let state2 = general_hopping(state, hop); + match state2 { + None => { + None + }, + Some(x) => { + let (ratio_ip, updated_column, col_idx) = { + fast_internal_product_spin_change(state, &x, pfaff_state, &hop, previous_proj, params) + }; + Some((ratio_ip, x, updated_column, col_idx)) + }, + } +} + +// C^t_{i,sigmai} C_{j,sigmaj} +#[inline(always)] +fn opperate_by_correlator + + Send> +( + state: &FockState, + pfaff_state: &PfaffianState, + previous_proj: &mut f64, + hop: &(usize, usize, Spin, Spin), + params: &VarParams, +) -> Option<(f64, FockState, Vec, usize)> +{ + let state2 = general_hopping(state, hop); + match state2 { + None => { + None + }, + Some(x) => { + let (ratio_ip, updated_column, col_idx) = { + fast_internal_product_spin_change(state, &x, pfaff_state, &hop, previous_proj, params) + }; + Some((ratio_ip, x, updated_column, col_idx)) + }, + } +} + +#[inline(always)] +fn accumulate_correlators + Send> +(state: FockState, pstate: &PfaffianState, proj: f64, params: &VarParams, n_sites: usize, expval_corr: &mut [f64]) +{ + // Loop over all correlators + // These loops are unrolled to allow the compiler to optimise hard coded spins => less branches + // since everything is inlined + // up up + for i in 0..n_sites { + for j in 0..n_sites { + let corr = (i, j, Spin::Up, Spin::Up); + let mut proj_copy = proj; + let (ratio, _, _, _) = match opperate_by_correlator(&state, pstate, &mut proj_copy, &corr, params) { + None => { + // Accumulate nothing if hopping is impossible. + continue; + }, + Some(r) => { + r + } + }; + let total_ratio = ratio * ::exp(proj_copy - proj); + expval_corr[j + i * n_sites] += total_ratio; + } + } + // up down + for i in 0..n_sites { + for j in 0..n_sites { + let corr = (i, j, Spin::Up, Spin::Down); + let mut proj_copy = proj; + let (ratio, _, _, _) = match opperate_by_correlator(&state, pstate, &mut proj_copy, &corr, params) { + None => { + // Accumulate nothing if hopping is impossible. + continue; + }, + Some(r) => { + r + } + }; + let total_ratio = ratio * ::exp(proj_copy - proj); + expval_corr[j + i * n_sites + n_sites*n_sites] += total_ratio; + } + } + // down up + for i in 0..n_sites { + for j in 0..n_sites { + let corr = (i, j, Spin::Down, Spin::Up); + let mut proj_copy = proj; + let (ratio, _, _, _) = match opperate_by_correlator(&state, pstate, &mut proj_copy, &corr, params) { + None => { + // Accumulate nothing if hopping is impossible. + continue; + }, + Some(r) => { + r + } + }; + let total_ratio = ratio * ::exp(proj_copy - proj); + expval_corr[j + i * n_sites + 2*n_sites*n_sites] += total_ratio; + } + } + // down down + for i in 0..n_sites { + for j in 0..n_sites { + let corr = (i, j, Spin::Down, Spin::Down); + let mut proj_copy = proj; + let (ratio, _, _, _) = match opperate_by_correlator(&state, pstate, &mut proj_copy, &corr, params) { + None => { + // Accumulate nothing if hopping is impossible. + continue; + }, + Some(r) => { + r + } + }; + let total_ratio = ratio * ::exp(proj_copy - proj); + expval_corr[j + i * n_sites + 3*n_sites*n_sites] += total_ratio; + } + } +} + +// C^t C +pub fn compute_mean_correlator + + Send> +(rng: &mut R, initial_state: FockState, projection: Projector, params: &VarParams, sys: &SysParams) -> (f64, Vec>) +where Standard: Distribution +{ + let mut expval_correlators = vec![0.0; 4*sys.size*sys.size]; + let mut state = initial_state; + let mut accumulated_states: Vec> = Vec::new(); + let (mut pstate, mut proj) = compute_internal_product_parts(state, params, sys); + 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 _energy_sq: f64 = 0.0; + let mut proj_copy_persistent = proj; + let mut ratio_prod = 1.0; + let mut sample_counter: usize = 0; + if ::log2(sys.nmcsample as f64) <= 6.0 { + error!("Not enough monte-carlo sample for an accurate error estimation. NMCSAMPLE = {}", sys.nmcsample); + panic!("NMCSAMPLE is less than 64."); + } + let mut energy_bootstraped = 0.0; + + if sys.nwarmupchains > sys._opt_iter { + info!("Starting the warmup phase."); + warmup(rng, &mut state, &mut hop, &mut proj, &mut proj_copy_persistent, &mut ratio_prod, &mut + pstate, &mut n_accepted_updates, params, sys); + } + + info!("Starting the sampling phase."); + // MC Sampling + // Accumulate the first state into the markov chain + accumulated_states.push(state); + let state_energy = compute_hamiltonian(state, &pstate, proj, params, sys); + + accumulate_expvals(&mut energy, state_energy, 1.0); + for mc_it in 0..(sys.nmcsample * sys.mcsample_interval) { + let mut proj_copy = proj; + let (mut ratio, state2, col, colidx) = propose_hopping(&state, &pstate, &mut proj_copy, &mut hop, rng, params, sys); + ratio *= ::exp(proj_copy - proj); + let w = rng.gen::(); + if ::abs(ratio) * ::abs(ratio) >= w { + // We ACCEPT + trace!("Accept."); + n_accepted_updates += 1; + make_update( + &mut n_accepted_updates, + &mut proj, + &mut proj_copy, + &mut proj_copy_persistent, + &ratio, + &mut ratio_prod, + &mut state, + &state2, + &hop, + col, + colidx, + &mut pstate, + params, + sys + ); + } + if sample_counter >= sys.mcsample_interval { + accumulated_states.push(state); + match projection { + Projector::Identity => { + accumulate_correlators(state, &pstate, proj, ¶ms, sys.size, &mut expval_correlators); + }, + Projector::ElectronExcitation { i, sigmai } => { + let excitation = (i, i, sigmai, sigmai); + let x = general_hopping(&state, &excitation); + match x { + None => { + // If None, then n_isigmai |x> = 0. Accumulate nothing + }, + Some(_) => { + accumulate_correlators(state, &pstate, proj, ¶ms, sys.size, &mut expval_correlators); + }, + } + }, + Projector::PairExcitation { i, sigmai, j, sigmaj } => { + let excitation1 = (i, i, sigmai, sigmai); + let excitation2 = (j, j, sigmaj, sigmaj); + let res1 = general_hopping(&state, &excitation1); + let res2 = general_hopping(&state, &excitation2); + match res1 { + None => { + // If None, then n_isigmai |x> = 0. Accumulate nothing + }, + Some(_) => { + match res2 { + None => { + // If None, then n_jsigmaj |x> = 0. Accumulate nothing + }, + Some(_) => { + accumulate_correlators(state, &pstate, proj, ¶ms, sys.size, &mut expval_correlators); + }, + } + }, + } + }, + } + let state_energy = compute_hamiltonian(state, &pstate, proj, params, sys); + + accumulate_expvals(&mut energy, state_energy, 1.0); + sample_counter = 0; + if mc_it >= (sys.nmcsample - sys.nbootstrap) * sys.mcsample_interval { + energy_bootstraped += energy; + } + } + sample_counter += 1; + + } + info!("Final Energy: {:.2}", energy); + normalize(&mut energy, &mut energy_bootstraped, sys.nmcsample as f64, sys.nbootstrap as f64); + info!("Final Energy normalized: {:.2}", energy); + (energy_bootstraped, accumulated_states) +} diff --git a/src/gutzwiller.rs b/src/gutzwiller.rs index a79664b..f7d0011 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, DerivativeOperator, FockState}; +use crate::{BitOps, DerivativeOperator, FockState, Spin}; /// Computes the gutzwiller factor for a single fock state. /// # Arguments @@ -117,6 +117,44 @@ pub fn fast_update_gutzwiller( } } +#[inline(always)] +pub fn fast_update_gutzwiller_spin_change( + previous_gutz: &mut f64, + gutzwiller_params: &[f64], + previous_fock: &FockState, + previous_index: usize, + new_index: usize, + previous_spin: Spin, + new_spin: Spin, +) where + T: BitOps + std::fmt::Display + std::marker::Send, +{ + match previous_spin { + Spin::Up => { + if previous_fock.spin_down.check(previous_index) { + *previous_gutz -= gutzwiller_params[previous_index]; + } + }, + Spin::Down => { + if previous_fock.spin_up.check(previous_index) { + *previous_gutz -= gutzwiller_params[previous_index]; + } + }, + } + match new_spin { + Spin::Up => { + if previous_fock.spin_down.check(previous_index) { + *previous_gutz += gutzwiller_params[new_index]; + } + }, + Spin::Down => { + if previous_fock.spin_up.check(previous_index) { + *previous_gutz += gutzwiller_params[new_index]; + } + }, + } +} + #[cfg(feature = "python-interface")] #[pyfunction] pub fn gutzwiller_exponent( @@ -382,6 +420,7 @@ mod test { let mut rng = SmallRng::seed_from_u64(42); // This is a random test, run it five times. for test_iter in 0..10 { + println!("Test iter = {}", test_iter); let mut e_up = [0; ARRAY_SIZE]; let mut e_down = [0; ARRAY_SIZE]; for i in 0..ARRAY_SIZE { diff --git a/src/jastrow.rs b/src/jastrow.rs index ba90bcf..0ebac90 100644 --- a/src/jastrow.rs +++ b/src/jastrow.rs @@ -87,15 +87,19 @@ where let k = indices[nk]; if n1.check(i) ^ n2.check(k) { if k > i { - der.o_tilde[der.jas_off + k*(k-1)/2 + i + (der.n * der.mu) as usize] = -0.5; + //der.o_tilde[der.jas_off + k*(k-1)/2 + i + (der.n * der.mu) as usize] = -0.5; + der.o_tilde[der.jas_off + k*(k-1)/2 + i + (der.n * der.mu) as usize] = -1.0; } else { - der.o_tilde[der.jas_off + i*(i-1)/2 + k + (der.n * der.mu) as usize] = -0.5; + //der.o_tilde[der.jas_off + i*(i-1)/2 + k + (der.n * der.mu) as usize] = -0.5; + der.o_tilde[der.jas_off + i*(i-1)/2 + k + (der.n * der.mu) as usize] = -1.0; } } else { if k > i { - der.o_tilde[der.jas_off + k*(k-1)/2 + i + (der.n * der.mu) as usize] = 0.5; + //der.o_tilde[der.jas_off + k*(k-1)/2 + i + (der.n * der.mu) as usize] = 0.5; + der.o_tilde[der.jas_off + k*(k-1)/2 + i + (der.n * der.mu) as usize] = 1.0; } else { - der.o_tilde[der.jas_off + i*(i-1)/2 + k + (der.n * der.mu) as usize] = 0.5; + //der.o_tilde[der.jas_off + i*(i-1)/2 + k + (der.n * der.mu) as usize] = 0.5; + der.o_tilde[der.jas_off + i*(i-1)/2 + k + (der.n * der.mu) as usize] = 1.0; } } } @@ -445,7 +449,7 @@ mod test { n_sites: usize ) -> f64 where - T: BitOps + From, + T: BitOps + From + Send + Sync + std::fmt::Display, { let mut zeta: Vec = vec![]; for i in 0..n_sites { diff --git a/src/lib.rs b/src/lib.rs index faa762d..1267d7d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -342,6 +342,8 @@ pub mod optimisation; /// TODOC pub mod dvmc; +pub mod green; + #[cfg(feature = "python-interface")] #[pymodule] fn impurity(m: &Bound<'_, PyModule>) -> PyResult<()> { diff --git a/src/monte_carlo.rs b/src/monte_carlo.rs index 76c6ccd..0775787 100644 --- a/src/monte_carlo.rs +++ b/src/monte_carlo.rs @@ -238,12 +238,40 @@ pub fn compute_mean_energy_exact(params: &VarParams, sys: &SysParams, der: &mut error!("The derivative operator current row was mu = {} on entry, is it reinitialized?", der.mu); } //let all_states = initial_state.generate_all_hoppings(sys.hopping_bitmask); - let all_states = vec![ - FockState{spin_up: 128u8, spin_down: 128u8, n_sites: 2}, - FockState{spin_up: 64u8, spin_down: 128u8, n_sites: 2}, - FockState{spin_up: 64u8, spin_down: 64u8, n_sites: 2}, - FockState{spin_up: 128u8, spin_down: 64u8, n_sites: 2}, - ]; + // Compute dimension of the available states. + let n_factorial = (1..sys.size + 1) + .try_fold(1usize, |acc, x| acc.checked_mul(x)) + .expect("Size too big for exact Sum. Overflowed in number of states to compute."); + let ne_factorial = (1..sys.nelec/2 + 1) + .try_fold(1usize, |acc, x| acc.checked_mul(x)) + .expect("Size too big for exact Sum. Overflowed in number of states to compute."); + let n_ne_factorial = (1..sys.size - sys.nelec/2 + 1) + .try_fold(1usize, |acc, x| acc.checked_mul(x)) + .expect("Size too big for exact Sum. Overflowed in number of states to compute."); + let n_choose_ne = n_factorial / ( + ne_factorial.checked_mul(n_ne_factorial) + .expect("Size too big for exact Sum. Overflowed in number of states to compute.") + ); + let dim = n_choose_ne * n_choose_ne; + + println!("Number of states to compute: {}", dim); + + // Generate all states. + let mut states: Vec = vec![]; + for mut s in (0..=::MAX).step_by(1<<(8-sys.size)) { + s.mask_bits(sys.size); + if s.count_ones()as usize == sys.nelec/2 { + states.push(s); + } + } + let mut all_states: Vec> = vec![]; + for i in 0..n_choose_ne { + for j in 0..n_choose_ne { + all_states.push( + FockState { spin_up: states[i], spin_down: states[j], n_sites: sys.size } + ); + } + } let mut energy = 0.0; let mut norm = 0.0; @@ -254,34 +282,48 @@ pub fn compute_mean_energy_exact(params: &VarParams, sys: &SysParams, der: &mut let (pstate, proj) = compute_internal_product_parts(*state2, params, sys); norm += sq(::abs(pstate.pfaff) * ::exp(proj)); let rho = ::abs(pstate.pfaff) * ::exp(proj); - println!("rho = {}", sq(rho)); compute_derivative_operator(*state2, &pstate, der, sys); state_energy = compute_hamiltonian(*state2, &pstate, proj, params, sys) * sq(rho); accumulate_expvals(&mut energy, state_energy, der, sq(rho)); der.visited[der.mu as usize] = 1; - let mut outstr = "".to_owned(); - outstr.push_str(&format!("H_mu O_k mu =")); for i in 0..der.n as usize { - outstr.push_str(&format!("{} ", der.o_tilde[i + (der.n * der.mu) as usize] * state_energy)); der.o_tilde[i + (der.n * der.mu) as usize] *= rho; } - println!("{}", outstr); der.mu += 1; } - println!("energy = {}", energy); for i in 0.. der.n as usize { der.expval_o[i] *= 1.0 / norm; der.ho[i] *= 1.0 / norm; } + //println!("OTILDE = \n{}", _save_otilde(&der.o_tilde, der.mu as usize, der.n as usize)); for i in 0..der.mu as usize { for j in 0..der.n as usize { der.o_tilde[j + i * der.n as usize] *= 1.0 / ::sqrt(norm); } } + //println!(" = {:?}", der.expval_o); + //println!(" = {:?}", der.ho); + //println!(" = {:?}", energy / norm); energy / norm } +fn _save_otilde(der: &[f64], mu: usize, n: usize) -> String { + let width = 11; + let mut o_tilde = "".to_owned(); + for m in 0..mu { + for i in 0..n { + if i == m { + o_tilde.push_str(&format!("{:>width$.04e}", der[i + m * n])); + } else { + o_tilde.push_str(&format!("{:>width$.04e}", der[i + m * n])); + } + } + o_tilde.push_str("\n"); + } + o_tilde +} + pub fn compute_mean_energy + Send> diff --git a/src/optimisation.rs b/src/optimisation.rs index 9501b79..2189969 100644 --- a/src/optimisation.rs +++ b/src/optimisation.rs @@ -385,13 +385,14 @@ impl<'a> ReducibleGeneralRepresentation for GenParameterMap { } } -fn gradient(x: &[f64], otilde: &[f64], visited: &[usize], expval_o: &[f64], b: &mut [f64], epsilon: f64, dim: i32, mu: i32, nsamp: f64) { +fn gradient(x: &[f64], otilde: &[f64], visited: &[usize], expval_o: &[f64], b: &mut [f64], diag_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, epsilon, dim, mu, nsamp); + compute_w(&mut work, otilde, visited, expval_o, x, diag_epsilon, dim, mu, nsamp); + //println!("Ax = {:?}", work); unsafe { // Compute b - Ax daxpy(dim, alpha, &work, incx, b, incy); @@ -406,8 +407,12 @@ 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], epsilon: f64, dim: i32, mu: i32, nsamp: f64) { +fn compute_w(w: &mut [f64], otilde: &[f64], visited: &[usize], expval_o: &[f64], p: &[f64], diag_epsilon: &[f64], dim: i32, mu: i32, nsamp: f64) { // Computes Ap + if dim == 0 { + error!("Cannot compute the matrix product of dimension 0. Something happened with the cutting of dimensions that is not accounted for"); + panic!("Undefined behavior."); + } let incx = 1; let incy = 1; let alpha = 1.0; @@ -429,12 +434,16 @@ 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, epsilon, w, incy); + dgemv(b"N"[0], dim, mu, 1.0 / nsamp, otilde, dim, &work, incx, beta, w, incy); + //dgemv(b"N"[0], dim, mu, 1.0, otilde, dim, &work, incx, beta, w, incy); trace!("O_[m, mu] O^[T]_[mu, n] x_[n] = {:?}", w); let alpha = ddot(dim, expval_o, incx, p, incy); // 81 misawa daxpy(dim, - alpha, expval_o, incx, w, incy); } + for i in 0..dim as usize { + w[i] += p[i] * diag_epsilon[i]; + } } fn update_r(rk: &mut [f64], w: &[f64], alphak: f64, dim: i32) { @@ -497,44 +506,69 @@ pub fn spread_eigenvalues(a: &mut DerivativeOperator) { } } -fn prefilter_overlap_matrix(a: &DerivativeOperator, _ignore_idx: &mut [bool], dim: i32, _diag_threshold: f64) -> usize { +fn prefilter_overlap_matrix(a: &DerivativeOperator, _ignore_idx: &mut [bool], dim: i32, _diag_threshold: f64, epsilon: f64, filter_before_shift: bool) -> (usize, Box<[f64]>) { // Loop over diagonal elements of S_km // Reminder, S_{kk} = 1/N_{\rm MC.} \sum_\mu \tilde{O}^*_{k\mu}\tilde{O}^T_{\mu k} - // \Re{\expval{O_k}}^2 - let skip_param_count: usize = 0; - let mut diag_elem = vec![0.0; dim as usize]; + let mut skip_param_count: usize = 0; + let mut diag_elem = vec![0.0; dim as usize].into_boxed_slice(); + let mut diag_eps = vec![0.0; dim as usize].into_boxed_slice(); for k in 0..dim as usize { // Start with 1/N_{\rm MC.} \sum_\mu \tilde{O}^*_{k\mu}\tilde{O}^T_{\mu k} - let z1: f64 = unsafe { - ddot( - a.mu, - &a.o_tilde[k..a.mu as usize * (k + 1)], - a.n, - &a.o_tilde[k..a.mu as usize * (k + 1)], - a.n - ) + let z1: f64 = { + let mut tmp = 0.0; + for j in 0..a.mu as usize { + let elem = a.o_tilde[k + j*(dim as usize)]; + tmp += elem * elem * a.visited[j] as f64; + } + tmp + //tmp / a.mu as f64 }; // Now \Re{\expval{O_k}}^2 let z2: f64 = a.expval_o[k] * a.expval_o[k]; - diag_elem[k] = z1 - z2; + if filter_before_shift { + diag_elem[k] = z1 - z2; + } else { + diag_elem[k] = (z1 - z2) * (1.0 + epsilon); + } + } let mut max_elem = ::MIN; + let mut min_elem = ::MAX; for k in 0..dim as usize { if diag_elem[k] > max_elem { max_elem = diag_elem[k]; } + if diag_elem[k] < min_elem { + min_elem = diag_elem[k]; + } + } + if min_elem < _diag_threshold * max_elem { + for k in 0..dim as usize { + if diag_elem[k] < _diag_threshold * max_elem { + skip_param_count += 1; + _ignore_idx[k] = true; + } + } + } + let mut n = 0; + for i in 0..dim as usize { + if _ignore_idx[i] { + continue; + } + if filter_before_shift { + diag_eps[n] = diag_elem[i] * epsilon; + } else { + diag_eps[n] = diag_elem[i] * epsilon / (1.0 + epsilon); + } + n += 1; } - //for k in 0..dim as usize { - // if diag_elem[k] < threshold { - // skip_param_count += 1; - // ignore_idx[k] = true; - // } - //} - dim as usize - skip_param_count + + (dim as usize - skip_param_count, diag_eps) } @@ -545,15 +579,19 @@ fn cpy_segmented_matrix_to_dense(a: &DerivativeOperator, output_otilde: &mut [f6 panic!("Will panic during this call during the copy."); } for k in 0..dim as usize { + if j > nparams_opt { + error!("Supplied dimension does not match allocated memory."); + panic!("Undefined behavior"); + } if ignore_idx[k] { continue; } unsafe { dcopy( a.mu, - &a.o_tilde[k..a.mu as usize * (k + 1)], + &a.o_tilde[k..(a.mu * a.n) as usize + k - a.n as usize], a.n, - &mut output_otilde[j..a.mu as usize * (j+1)], + &mut output_otilde[j..nparams_opt * a.mu as usize - nparams_opt + j], nparams_opt as i32 ); } @@ -567,6 +605,8 @@ fn compute_s_explicit(otilde: &[f64], expval_o: &[f64], visited: &[usize], dim: // Computes Ap let incx = 1; let alpha = 1.0/nsamp; + // Duumbo + //let alpha = 1.0; let beta = -(1.0); let mut otilde_w_visited = vec![0.0; (dim * mu) as usize]; for i in 0..mu as usize { @@ -574,17 +614,22 @@ fn compute_s_explicit(otilde: &[f64], expval_o: &[f64], visited: &[usize], dim: otilde_w_visited[j + i*dim as usize] = otilde[j + i*dim as usize] * visited[i] as f64; } } + //println!("{}", _save_otilde(õ, mu as usize, dim as usize)); + //println!("{:?}", visited); // Temp work vector let mut work = vec![0.0; (dim * dim) as usize]; unsafe { // 80 misawa dger(dim, dim, 1.0, &expval_o, incx, &expval_o, incx, &mut work, dim); - dgemm(b"N"[0], b"T"[0], dim, dim, mu, alpha, õ_w_visited, dim, õ_w_visited, dim, beta, &mut work, dim); + dgemm(b"N"[0], b"T"[0], dim, dim, mu, alpha, õ, dim, õ_w_visited, dim, beta, &mut work, dim); } // Shift smallest eigenvalues + //println!("Before diag filter"); + //println!("{}", _save_otilde(&work, dim as usize, dim as usize)); for i in 0..dim as usize { - work[i + dim as usize * i] += epsilon; + work[i + dim as usize * i] *= 1.0 + epsilon; + //work[i + dim as usize * i] += epsilon; } work } @@ -668,7 +713,7 @@ fn _compute_matrix_product(s: &mut [f64], eigenvectors: &[f64], eigenvalues: &[f } for i in 0..dim as usize{ for j in 0..dim as usize { - if eigenvalues[i] < 1e-1 { + if eigenvalues[i] < 1e-6 { work[j + dim as usize * i] *= 0.0; work_direct[j + dim as usize * i] *= 0.0; continue; @@ -691,17 +736,53 @@ fn _compute_matrix_product(s: &mut [f64], eigenvectors: &[f64], eigenvalues: &[f pub fn exact_overlap_inverse(a: &DerivativeOperator, b: &mut [f64], epsilon: f64, dim: i32, thresh: f64) -> Vec{ // PRE FILTER let mut ignore = vec![false; dim as usize]; - //println!("{}", save_otilde(a.o_tilde, a.mu as usize, a.n as usize)); - let mut unfiltered_s = compute_s_explicit(&a.o_tilde, &a.expval_o, &a.visited, dim, a.mu, a.nsamp, epsilon); //println!("dim = {}, Unfiltered S = ", dim); //println!("{}", save_otilde(&unfiltered_s, dim as usize, dim as usize)); - let new_dim = prefilter_overlap_matrix(a, &mut ignore, dim, thresh); - //println!("ignore : {:?}", ignore); - let mut otilde = vec![0.0; new_dim * a.mu as usize]; - let mut expvalo = vec![0.0; new_dim]; + let (new_dim, _) = prefilter_overlap_matrix(a, &mut ignore, dim, thresh, epsilon, true); + let mut otilde = vec![0.0; new_dim * a.mu as usize].into_boxed_slice(); + let mut expvalo = vec![0.0; new_dim].into_boxed_slice(); cpy_segmented_matrix_to_dense(a, &mut otilde, &mut expvalo, &ignore, dim, new_dim); + //println!("{}", _save_otilde(&a.o_tilde, a.mu as usize, a.n as usize)); + let mut filtered_s = compute_s_explicit(õ, &expvalo, &a.visited, new_dim as i32, a.mu, a.nsamp, epsilon); + let mut work = vec![0.0; new_dim]; + unsafe { + dgemv(b"N"[0], new_dim as i32, new_dim as i32, 1.0, &filtered_s, new_dim as i32, &b, 1, 0.0, &mut work, 1); + } + println!("Ax = {:?}", work); + println!("{}", _save_otilde(&filtered_s, new_dim as usize, new_dim as usize)); + println!("a.mu = {}, a.n = {}", a.mu, a.n); + println!("dim = {}, new_dim = {}", dim, new_dim); + let mut _unfiltered_s = compute_s_explicit(&a.o_tilde, &a.expval_o, &a.visited, dim as i32, a.mu, a.nsamp, epsilon); + //for i in 0..dim as usize { + // if ignore[i] { + // println!("{}", _save_otilde(&filtered_s, new_dim, new_dim)); + // println!("ignore : {:?}", ignore); + // println!("{}", _save_otilde(&unfiltered_s, dim as usize, dim as usize)); + // panic!("stop"); + // break; + + // } + //} + // Test matrix are equal + //for i in 0..dim as usize { + // for j in 0.. dim as usize { + // if ! (filtered_s[j + i * dim as usize] == unfiltered_s[j + i * dim as usize]) { + // panic!("Assertion"); + // } + // } + //} + //println!("dim = {}, Unfiltered S = ", dim); + //println!("{}", _save_otilde(&unfiltered_s, dim as usize, dim as usize)); + //let new_dim = prefilter_overlap_matrix(a, &mut ignore, dim, thresh); + ////println!("ignore : {:?}", ignore); + //let mut otilde = vec![0.0; new_dim * a.mu as usize]; + //let mut expvalo = vec![0.0; new_dim]; + //cpy_segmented_matrix_to_dense(a, &mut otilde, &mut expvalo, &ignore, dim, new_dim); //let filtered_s = compute_s_explicit(õ, &expvalo, a.visited, new_dim as i32, a.mu, a.nsamp); - let mut eigenvalues = diagonalize_dense_matrix(&mut unfiltered_s, dim); + let eigenvalues = diagonalize_dense_matrix(&mut filtered_s, new_dim as i32); + let _eigenvalues = diagonalize_dense_matrix(&mut _unfiltered_s, dim as i32); + //_compute_matrix_product(&mut s_copy, &unfiltered_s, &eigenvalues, dim); + //println!("S^[-1] = \n{}", _save_otilde(&s_copy, dim as usize, dim as usize)); //let eigenvectors = &unfiltered_s; //println!("S = \n{}", save_otilde(&s_copy, dim as usize, dim as usize)); //compute_matrix_product(&mut s_copy, &eigenvectors, &eigenvalues, dim); @@ -715,32 +796,37 @@ pub fn exact_overlap_inverse(a: &DerivativeOperator, b: &mut [f64], epsilon: f64 //} //println!("UD^[-1]U^[T] = \n{}", save_otilde(&s_copy, dim as usize, dim as usize)); //println!("x0 = {:?}", x0_raw); - //println!("eigenvalues: {:?}", eigenvalues); + println!("filted eigenvalues: {:?}", eigenvalues); + println!("eigenvalues: {:?}", _eigenvalues); // Remove problematic eigenvalue - let mut max = ::MIN; - for e in eigenvalues.iter() { - if *e > max { - max = *e; - } - } - let threshold = thresh * max; - let mut i = 0; - for e in eigenvalues.iter_mut() { - if *e < threshold { - *e = 0.0; - ignore[i] = true; - } - i+=1; - } - //Invert matrix - for e in eigenvalues.iter_mut() { - if *e == 0.0 { - continue; - } - *e = 1.0 / *e; - } - compute_delta_from_eigenvalues(b, &unfiltered_s, &eigenvalues, dim); + //let mut max = ::MIN; + ////println!("{:?}", eigenvalues); + //for e in eigenvalues.iter() { + // if *e > max { + // max = *e; + // } + //} + //let threshold = thresh * max; + //let mut i = 0; + //for e in eigenvalues.iter_mut() { + // if *e < threshold { + // *e = 0.0; + // //ignore[i] = true; + // } + // i+=1; + //} + ////println!("{:?}", ignore); + ////Invert matrix + //for e in eigenvalues.iter_mut() { + // if *e == 0.0 { + // continue; + // } + // *e = 1.0 / *e; + //} + compute_delta_from_eigenvalues(b, &filtered_s, &eigenvalues, new_dim as i32); + //println!("b = {:?}", b); + return ignore; } @@ -748,18 +834,22 @@ pub fn exact_overlap_inverse(a: &DerivativeOperator, b: &mut [f64], epsilon: f64 /// TODOC /// 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 +pub fn conjugate_gradiant(a: &DerivativeOperator, b: &mut [f64], x0: &mut [f64], epsilon: f64, kmax: usize, dim: i32, thresh: f64, epsilon_convergence: f64, filter_before_shift: bool) -> Vec { // PRE FILTER let mut ignore = vec![false; dim as usize]; //println!("{}", save_otilde(a.o_tilde, a.mu as usize, a.n as usize)); //println!("dim = {}, Unfiltered S = ", dim); //println!("{}", save_otilde(&unfiltered_s, dim as usize, dim as usize)); - let new_dim = prefilter_overlap_matrix(a, &mut ignore, dim, thresh); + let (new_dim, diag_epsilon) = prefilter_overlap_matrix(a, &mut ignore, dim, thresh, epsilon, filter_before_shift); + //println!("diag : {:?}", diag_epsilon); //println!("ignore : {:?}", ignore); - let mut otilde = vec![0.0; new_dim * a.mu as usize]; - let mut expvalo = vec![0.0; new_dim]; + let mut otilde = vec![0.0; new_dim * a.mu as usize].into_boxed_slice(); + let mut expvalo = vec![0.0; new_dim].into_boxed_slice(); cpy_segmented_matrix_to_dense(a, &mut otilde, &mut expvalo, &ignore, dim, new_dim); + //let mut filtered_s = compute_s_explicit(õ, &expvalo, &a.visited, new_dim as i32, a.mu, a.nsamp, epsilon); + //println!("{}", _save_otilde(&filtered_s, new_dim as usize, new_dim as usize)); + //let diag_epsilon = compute_diag_epsilon(õ, &expvalo, epsilon, new_dim, a.mu); //let filtered_s = compute_s_explicit(õ, &expvalo, a.visited, new_dim as i32, a.mu, a.nsamp); //let eigenvectors = &unfiltered_s; //println!("S = \n{}", save_otilde(&s_copy, dim as usize, dim as usize)); @@ -782,9 +872,16 @@ pub fn conjugate_gradiant(a: &DerivativeOperator, b: &mut [f64], x0: &mut [f64], //println!("{}", save_otilde(&filtered_s, new_dim as usize, new_dim as usize)); //println!("{}", save_otilde(õ, a.mu as usize, new_dim as usize)); + //unsafe { + // dcopy(new_dim as i32, &b, 1, x0, 1); + //} + unsafe { + dscal(new_dim as i32, 0.0, x0, 1); + } + trace!("Initial guess x_0: {:?}", x0); trace!("Initial equation b: {:?}", b); - let mut w: Vec = vec![0.0; new_dim]; + let mut w = vec![0.0; new_dim].into_boxed_slice(); // Error threshold let mut e = 0.0; for i in 0..dim as usize { @@ -796,8 +893,8 @@ 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, epsilon, new_dim as i32, a.mu, a.nsamp); - let mut p = vec![0.0; new_dim]; + gradient(x0, õ, &a.visited, &expvalo, b, &diag_epsilon, new_dim as i32, a.mu, a.nsamp); + let mut p = vec![0.0; new_dim].into_boxed_slice(); let mut j: usize = 0; for i in 0..dim as usize { if ignore[i] { @@ -816,15 +913,15 @@ 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, epsilon, new_dim as i32, a.mu, a.nsamp); + compute_w(&mut w, õ, &a.visited, &expvalo, &p, &diag_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); //println!("alpha_{} : {}", k, alpha); //if alpha < 0.0 { - // //error!("Input overlap matrix S was not positive-definite."); - // break; - // //panic!("p^T S p < 0.0"); + // error!("Input overlap matrix S was not positive-definite."); + // //break; + // panic!("p^T S p < 0.0"); //} update_x(x0, &p, alpha, new_dim as i32); trace!("x_{} : {:?}", k+1, x0); @@ -841,5 +938,6 @@ pub fn conjugate_gradiant(a: &DerivativeOperator, b: &mut [f64], x0: &mut [f64], unsafe { dcopy(new_dim as i32, x0, 1, b, 1); } + //println!("b = {:?}", b); ignore } diff --git a/src/pfaffian.rs b/src/pfaffian.rs index 684332b..b8e21f3 100644 --- a/src/pfaffian.rs +++ b/src/pfaffian.rs @@ -229,7 +229,9 @@ pub fn compute_pfaffian_derivative(pstate: &PfaffianState, der: &mut DerivativeO 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) as usize] = -a[ii * n + (jj + off)]; + //der.o_tilde[der.pfaff_off + indices2[jj] + size * indices[ii] + size*size + (der.n * der.mu) as usize] = a[ii * n + (jj + off)]; trace!("~O_[{}, {}] = {}", der.pfaff_off + indices2[jj] + size * indices[ii] + size * size, der.mu, -a[ii * n + (jj + off)]); + //der.o_tilde[der.pfaff_off + indices[ii] + size * indices2[jj] + 2*size*size + (der.n * der.mu) as usize] = a[(jj + off) * n + ii]; der.o_tilde[der.pfaff_off + indices[ii] + size * indices2[jj] + 2*size*size + (der.n * der.mu) as usize] = -a[(jj + off) * n + ii]; trace!("~O_[{}, {}] = {}", der.pfaff_off + indices[ii] + size * indices2[jj] + 2 * size * size, der.mu, -a[(jj + off) * n + ii]); } @@ -256,6 +258,190 @@ pub fn compute_pfaffian_derivative(pstate: &PfaffianState, der: &mut DerivativeO } } +/// Gets the ratio of pfaffian after an exchange. +/// # Fields +/// * __`previous_pstate`__ - The pfaffian state to update. +/// * __`previous_i`__ - The initial index of the jumping electron. +/// * __`new_i`__ - The index of the electron after the jump. +/// +/// # Returns +/// * __`pfaff_up`__ - The new pfaffian ratio after the update. +/// * __`new_b`__ - The new column and row in the matrix $A$. +/// * __`col`__ - The index of the column and row that changed in the matrix $A$. +pub fn get_pfaffian_ratio_exchange( + previous_pstate: &PfaffianState, + previous_i: usize, + new_i: usize, + previous_spin: Spin, + new_spin: Spin, + fij: &[f64], +) -> (f64, Vec, Vec, usize, usize) { + // Rename + let indx_up = &previous_pstate.indices.0; + trace!("Up : {:?}", indx_up); + let indx_down = &previous_pstate.indices.1; + trace!("Down : {:?}", indx_down); + let n_sites = previous_pstate.n_sites; + let n_elec = previous_pstate.n_elec; + + // Gen new vector b + // X_{ij}=F_{ij}^{\sigma_i,\sigma_j} - F_{ji}^{\sigma_j,\sigma_i}, tahara2008 + // +0 -> upup, +SIZE^2 -> updown, +2*SIZE^2 -> downup, +3*SIZE^2 -> down down + let mut new_b1: Vec = Vec::with_capacity(n_elec); + for iup in indx_up.iter() { + match previous_spin { + Spin::Up => { + if *iup == previous_i { + trace!("Pushed 0.0"); + new_b1.push(0.0); + continue; + } + trace!("Pushed X_[{}, {}], sector up up", new_i, iup); + new_b1.push( + fij[new_i + n_sites * iup] + -fij[iup + n_sites * new_i]); + } + Spin::Down => { + trace!("Pushed X_[{}, {}], sector up down", new_i, iup); + new_b1.push( + fij[new_i + n_sites * iup + n_sites*n_sites] + -fij[iup + n_sites * new_i + 2*n_sites*n_sites]); + } + }; + } + for idown in indx_down.iter() { + match previous_spin { + Spin::Up => { + trace!("Pushed X_[{}, {}], sector down up", new_i, idown); + new_b1.push( + fij[new_i + n_sites * idown + 2*n_sites*n_sites] + -fij[idown + n_sites * new_i + n_sites*n_sites]); + } + Spin::Down => { + if *idown == previous_i { + new_b1.push(0.0); + trace!("Pushed 0.0"); + continue; + } + trace!("Pushed X_[{}, {}], sector down up", new_i, idown); + new_b1.push( + fij[new_i + n_sites * idown + 3*n_sites*n_sites] + -fij[idown + n_sites * new_i + 3*n_sites*n_sites]); + } + }; + } + let mut new_b2: Vec = Vec::with_capacity(n_elec); + for iup in indx_up.iter() { + match new_spin { + Spin::Up => { + if *iup == previous_i { + trace!("Pushed 0.0"); + new_b2.push(0.0); + continue; + } + trace!("Pushed X_[{}, {}], sector up up", new_i, iup); + new_b2.push( + fij[new_i + n_sites * iup] + -fij[iup + n_sites * new_i]); + } + Spin::Down => { + trace!("Pushed X_[{}, {}], sector up down", new_i, iup); + new_b2.push( + fij[new_i + n_sites * iup + n_sites*n_sites] + -fij[iup + n_sites * new_i + 2*n_sites*n_sites]); + } + }; + } + for idown in indx_down.iter() { + match new_spin { + Spin::Up => { + trace!("Pushed X_[{}, {}], sector down up", new_i, idown); + new_b2.push( + fij[new_i + n_sites * idown + 2*n_sites*n_sites] + -fij[idown + n_sites * new_i + n_sites*n_sites]); + } + Spin::Down => { + if *idown == previous_i { + new_b2.push(0.0); + trace!("Pushed 0.0"); + continue; + } + trace!("Pushed X_[{}, {}], sector down up", new_i, idown); + new_b2.push( + fij[new_i + n_sites * idown + 3*n_sites*n_sites] + -fij[idown + n_sites * new_i + 3*n_sites*n_sites]); + } + }; + } + + // Get the column to replace. + trace!("Making hopping ({}, {}, {}, {})", previous_i, new_i, previous_spin, new_spin); + trace!("Index: up {:?}, down {:?}", indx_up, indx_down); + let col1 = match previous_spin { + Spin::Up => indx_up.iter().position(|&r| r == previous_i).unwrap(), + Spin::Down => indx_down.iter().position(|&r| r == previous_i).unwrap() + indx_up.len(), + }; + let col2 = match new_spin { + Spin::Up => indx_up.iter().position(|&r| r == new_i).unwrap(), + Spin::Down => indx_down.iter().position(|&r| r == new_i).unwrap() + indx_up.len(), + }; + + // Compute the updated pfaffian. + trace!("Need to update col {}", col1); + trace!("X^[-1] = {}", previous_pstate); + trace!("X^[-1]_[col, i] = {:?}", + &previous_pstate.inv_matrix[n_elec * col1..n_elec + n_elec * col1], +); + let mut c_matrix = vec![0.0; 4]; + unsafe { + c_matrix[0] = ddot( + n_elec as i32, + &new_b1, + 1, + &previous_pstate.inv_matrix[n_elec * col1..n_elec + n_elec * col1], + 1, + ); + c_matrix[1] = ddot( + n_elec as i32, + &new_b2, + 1, + &previous_pstate.inv_matrix[n_elec * col1..n_elec + n_elec * col1], + 1, + ); + c_matrix[2] = ddot( + n_elec as i32, + &new_b1, + 1, + &previous_pstate.inv_matrix[n_elec * col2..n_elec + n_elec * col2], + 1, + ); + c_matrix[3] = ddot( + n_elec as i32, + &new_b2, + 1, + &previous_pstate.inv_matrix[n_elec * col2..n_elec + n_elec * col2], + 1, + ); + } + // Determinant of update + let det = c_matrix[0] * c_matrix[3] - c_matrix[1] * c_matrix[2]; + // Correction + let mut y = vec![0.0; n_elec]; + let correction = unsafe { + let trans = b"N"[0]; + let m = n_elec as i32; + let incx = 0; + let incy = 0; + let alpha = previous_pstate.inv_matrix[col1 + n_elec * col2]; + let beta = 0.0; + dgemv(trans, m, m, alpha, &previous_pstate.inv_matrix, m, &new_b1, incx, beta, &mut y, incy); + ddot(m, &new_b2, incx, &y, incy) + }; + let pfaff_up = det + previous_pstate.inv_matrix[col1 + n_elec * col2]*new_b1[col2] + correction; + trace!("pfaffu_up = {}", pfaff_up); + (pfaff_up, new_b1, new_b2, col1, col2) +} + /// Gets the ratio of pfaffian after an update. /// # Fields /// * __`previous_pstate`__ - The pfaffian state to update. @@ -271,7 +457,8 @@ pub fn get_pfaffian_ratio( previous_pstate: &PfaffianState, previous_i: usize, new_i: usize, - spin: Spin, + previous_spin: Spin, + new_spin: Spin, fij: &[f64], ) -> (f64, Vec, usize) { // Rename @@ -287,7 +474,7 @@ pub fn get_pfaffian_ratio( // +0 -> upup, +SIZE^2 -> updown, +2*SIZE^2 -> downup, +3*SIZE^2 -> down down let mut new_b: Vec = Vec::with_capacity(n_elec); for iup in indx_up.iter() { - match spin { + match new_spin { Spin::Up => { if *iup == previous_i { trace!("Pushed 0.0"); @@ -308,7 +495,7 @@ pub fn get_pfaffian_ratio( }; } for idown in indx_down.iter() { - match spin { + match new_spin { Spin::Up => { trace!("Pushed X_[{}, {}], sector down up", new_i, idown); new_b.push( @@ -330,9 +517,9 @@ pub fn get_pfaffian_ratio( } // Get the column to replace. - trace!("Making hopping ({}, {}, {})", previous_i, new_i, spin); + trace!("Making hopping ({}, {}, {}, {})", previous_i, new_i, previous_spin, new_spin); trace!("Index: up {:?}, down {:?}", indx_up, indx_down); - let col = match spin { + let col = match previous_spin { Spin::Up => indx_up.iter().position(|&r| r == previous_i).unwrap(), Spin::Down => indx_down.iter().position(|&r| r == previous_i).unwrap() + indx_up.len(), }; @@ -755,10 +942,10 @@ mod tests { println!("Spin is up: {}", is_spin_up); let tmp = if is_spin_up { - get_pfaffian_ratio(&pfstate, initial_index, final_index, Spin::Up, ¶ms) + get_pfaffian_ratio(&pfstate, initial_index, final_index, Spin::Up, Spin::Up, ¶ms) } else { - get_pfaffian_ratio(&pfstate, initial_index, final_index, Spin::Down, ¶ms) + get_pfaffian_ratio(&pfstate, initial_index, final_index, Spin::Down, Spin::Down, ¶ms) }; println!("Ratio: {}", tmp.0); println!("B col: {:?}", tmp.1); @@ -885,7 +1072,7 @@ mod tests { }; let pfstate2 = construct_matrix_a_from_state(¶ms, state2, &sys); println!("Inverse Matrix: {}", pfstate2); - let pfaff_ratio = get_pfaffian_ratio(&pfstate, 6, 5, Spin::Up, ¶ms).0; + let pfaff_ratio = get_pfaffian_ratio(&pfstate, 6, 5, Spin::Up, Spin::Up, ¶ms).0; close(pfstate.pfaff * pfaff_ratio, pfstate2.pfaff, 1e-12); } @@ -965,7 +1152,7 @@ mod tests { }; let pfstate2 = construct_matrix_a_from_state(¶ms, state2, &sys); println!("Inverse Matrix: {}", pfstate2); - let tmp = get_pfaffian_ratio(&pfstate, 6, 5, Spin::Down, ¶ms); + let tmp = get_pfaffian_ratio(&pfstate, 6, 5, Spin::Down, Spin::Down, ¶ms); println!("B: {:?}", tmp.1); println!("Ratio: {}", tmp.0); close(pfstate.pfaff * tmp.0, pfstate2.pfaff, 1e-12); @@ -1051,7 +1238,7 @@ mod tests { let pfstate2 = construct_matrix_a_from_state(¶ms, state2, &sys); println!("Inverse Matrix: {}", pfstate2); println!("------------- Proposed Update ------------------"); - let tmp = get_pfaffian_ratio(&pfstate, 6, 5, Spin::Down, ¶ms); + let tmp = get_pfaffian_ratio(&pfstate, 6, 5, Spin::Down, Spin::Down, ¶ms); println!("Ratio: {}", tmp.0); println!("B col: {:?}", tmp.1); close(pfstate.pfaff * tmp.0, pfstate2.pfaff, 1e-12); diff --git a/tests/vmc_2sites.rs b/tests/vmc_2sites.rs index 0eb4024..907dc7f 100644 --- a/tests/vmc_2sites.rs +++ b/tests/vmc_2sites.rs @@ -248,7 +248,7 @@ fn analytic_derivatives_expval(par: &VarParams) -> Vec { let mut out_der = vec![0.0; SIZE + 1 + 4*SIZE*SIZE]; out_der[0] = sq(individual_state(&State::F10, par)); out_der[1] = sq(individual_state(&State::F5, par)); - out_der[2] = 0.5 * (- out_der[0] - out_der[1]); + out_der[2] = - out_der[0] - out_der[1]; //out_der[3] = //out_der[4] = //out_der[5] = @@ -287,7 +287,7 @@ fn analytic_ho_expval(par: &VarParams) -> Vec { * energy_individual_state(&State::F10, par); out_der[1] = individual_state(&State::F5, par) * energy_individual_state(&State::F5, par); - out_der[2] = 0.5 * (- out_der[0] - out_der[1]); + out_der[2] = - out_der[0] - out_der[1]; //out_der[3] = //out_der[4] = //out_der[5] = @@ -471,6 +471,7 @@ fn comupte_energy_from_all_states() { print_der(&der.expval_o, &exp_val, sys.ngi+sys.nvij+sys.nfij); let psi = norm(¶meters); println!("Norm: {:10.4e}", psi); + println!("Comparing , tol: {}", 1e-2); for i in 0..sys.ngi+sys.nvij+sys.nfij { close(der.expval_o[i] * psi, exp_val[i], 1e-2); } @@ -479,6 +480,7 @@ fn comupte_energy_from_all_states() { print_der(&der.ho, &exp_val_ho, sys.ngi+sys.nvij+sys.nfij); let psi = norm(¶meters); println!("Norm: {:10.4e}", psi); + println!("Comparing , tol: {}", 1e-2); for i in 0..sys.ngi+sys.nvij+sys.nfij { close(der.ho[i] * psi, exp_val_ho[i], 2e-2); } diff --git a/tests/vmc_2sites_analytic_conv.rs b/tests/vmc_2sites_analytic_conv.rs index 460f10d..cca4777 100644 --- a/tests/vmc_2sites_analytic_conv.rs +++ b/tests/vmc_2sites_analytic_conv.rs @@ -254,7 +254,7 @@ fn analytic_ho_expval(par: &VarParams) -> Vec { * energy_individual_state(&State::F10, par); out_der[1] = individual_state(&State::F5, par) * energy_individual_state(&State::F5, par); - out_der[2] = 0.5 * (- out_der[0] - out_der[1]); + out_der[2] = - out_der[0] - out_der[1]; //out_der[3] = //out_der[4] = //out_der[5] = @@ -279,7 +279,7 @@ fn analytic_derivatives_expval(par: &VarParams) -> Vec { let mut out_der = vec![0.0; SIZE + 1 + 4*SIZE*SIZE]; out_der[0] = sq(individual_state(&State::F10, par)); out_der[1] = sq(individual_state(&State::F5, par)); - out_der[2] = 0.5 * (- out_der[0] - out_der[1]); + out_der[2] = - out_der[0] - out_der[1]; // fij out_der[3] = 1.0 * { ::exp(par.gi[0] - par.vij[0]) / individual_state(&State::F10, par) diff --git a/tests/vmc_2sites_exact_sum.rs b/tests/vmc_2sites_exact_sum.rs index 29eb7b6..2b1afaa 100644 --- a/tests/vmc_2sites_exact_sum.rs +++ b/tests/vmc_2sites_exact_sum.rs @@ -249,7 +249,7 @@ fn analytic_ho_expval(par: &VarParams) -> Vec { * energy_individual_state(&State::F10, par); out_der[1] = individual_state(&State::F5, par) * energy_individual_state(&State::F5, par); - out_der[2] = 0.5 * (- out_der[0] - out_der[1]); + out_der[2] = - out_der[0] - out_der[1]; //out_der[3] = //out_der[4] = //out_der[5] = @@ -274,7 +274,7 @@ fn analytic_derivatives_expval(par: &VarParams) -> Vec { let mut out_der = vec![0.0; SIZE + 1 + 4*SIZE*SIZE]; out_der[0] = sq(individual_state(&State::F10, par)); out_der[1] = sq(individual_state(&State::F5, par)); - out_der[2] = 0.5 * (- out_der[0] - out_der[1]); + out_der[2] = - out_der[0] - out_der[1]; // fij out_der[3] = 1.0 * { ::exp(par.gi[0] - par.vij[0]) / individual_state(&State::F10, par)