Skip to content

(fix) Filtering issue with diagonal truncation #5

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Duumbo merged 1 commit into from
Jul 7, 2025
Merged

(fix) Filtering issue with diagonal truncation #5

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions src/bin/dvmc_pairwf.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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);
Expand All @@ -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];
Expand Down
120 changes: 76 additions & 44 deletions src/bin/dvmc_pairwf2.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,21 +15,21 @@ 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;
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;
Expand All @@ -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 = <f64>::EPSILON;
const PARAM_THRESHOLD: f64 = 1e-3;
//const PARAM_THRESHOLD: f64 = 0.0;
//const PARAM_THRESHOLD: f64 = -<f64>::INFINITY;
const OPTIMISE: bool = true;
const OPTIMISE_GUTZ: bool = true;
const OPTIMISE_JAST: bool = true;
Expand All @@ -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
Expand Down Expand Up @@ -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 {
<f64>::abs(a) * <f64>::abs(a)
}
Expand Down Expand Up @@ -242,6 +245,7 @@ fn main() {
let mut all_params: Vec<f64> = 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);
Expand All @@ -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;
Expand All @@ -270,6 +275,32 @@ fn main() {
parameters.fij[i +3*SIZE*SIZE] = 0.0;
}
}
let mut max = <f64>::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 {
Expand All @@ -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<BitSize> = {
Expand Down
38 changes: 35 additions & 3 deletions src/density.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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};

Expand Down Expand Up @@ -77,7 +77,39 @@ where T: BitOps + std::fmt::Debug + std::fmt::Display + From<SpinState> + 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'|pf>/<x|pf>: {}, |x'> = {}, |x> = {}", pfaffian_ratio, new_state, previous_state);
(pfaffian_ratio, b_vec, col)
}

pub fn fast_internal_product_spin_change<T>(
previous_state: &FockState<T>,
new_state: &FockState<T>,
previous_pstate: &PfaffianState,
hopping: &(usize, usize, Spin, Spin),
previous_proj: &mut f64,
params: &VarParams,
) -> (f64, Vec<f64>, usize)
where T: BitOps + std::fmt::Debug + std::fmt::Display + From<SpinState> + std::ops::Shl<usize, Output = T> + 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 = &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);
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.
Expand Down Expand Up @@ -116,7 +148,7 @@ where T: BitOps + std::fmt::Debug + std::fmt::Display + From<SpinState> + 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.
Expand Down
Loading
Loading