Skip to content

Feat pair wavefunction #3

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 4 commits into from
Apr 7, 2025
Merged

Feat pair wavefunction #3

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
6134638
(fix) MonteCarlo systematic error.
Mar 27, 2025
0a0cc66
(fix) Conjugate gradient problem
Mar 27, 2025
83e66e2
(feat) Pair wave-function final.
Apr 7, 2025
935e739
(fix) Cagro changes
Apr 7, 2025
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
3 changes: 3 additions & 0 deletions Cargo.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,9 @@ harness = false
[package.metadata.docs.rs]
rustdoc-args = [ "--html-in-header", "katex-header.html" ]

[profile.release]
debug = true

[features]
default = []
python-interface = ["dep:pyo3"]
Expand Down
13 changes: 10 additions & 3 deletions src/bin/dvmc_pairwf.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ const NVIJ: usize = SIZE*(SIZE - 1) / 2;
const NGI: usize = SIZE;
const NPARAMS: usize = NFIJ + NGI + NVIJ;
const NELEC: usize = SIZE;
const NMCSAMP: usize = 1_000;
const NMCSAMP: usize = 500;
const NBOOTSTRAP: usize = 1;
const NMCWARMUP: usize = 1000;
const MCSAMPLE_INTERVAL: usize = 1;
Expand All @@ -32,15 +32,15 @@ const OPTIMISATION_TIME_STEP: f64 = 1e-2;
const OPTIMISATION_DECAY: f64 = 0.0;
const NOPTITER: usize = 1000;
const KMAX: usize = NPARAMS;
const PARAM_THRESHOLD: f64 = <f64>::MIN;
const PARAM_THRESHOLD: f64 = 1e-5;
//const PARAM_THRESHOLD: f64 = 0.0;
const OPTIMISE: bool = true;
const OPTIMISE_GUTZ: bool = true;
const OPTIMISE_JAST: bool = true;
const OPTIMISE_ORB: bool = true;
const SET_EXPVALO_ZERO: bool = false;
const COMPUTE_ENERGY_METHOD: EnergyComputationMethod = EnergyComputationMethod::MonteCarlo;
const OPTIMISE_ENERGY_METHOD: EnergyOptimisationMethod = EnergyOptimisationMethod::ExactInverse;
const OPTIMISE_ENERGY_METHOD: EnergyOptimisationMethod = EnergyOptimisationMethod::ConjugateGradiant;

pub enum EnergyOptimisationMethod {
ExactInverse,
Expand Down Expand Up @@ -357,6 +357,13 @@ fn main() {
//spread_eigenvalues(&mut derivative);
//println!("x0 in = {:?}", x0);
let mut _flag: bool = true;
let mut bcopy = b.clone();
//let ignored_columns = exact_overlap_inverse(&work_derivative, &mut bcopy, EPSILON_SHIFT, NPARAMS as i32, PARAM_THRESHOLD);
//conjugate_gradiant(&work_derivative, &mut b, &mut x0, EPSILON_SHIFT, KMAX, NPARAMS as i32, PARAM_THRESHOLD, EPSILON_CG);
//println!("Exact b: {:?}", bcopy);
//println!("CG b: {:?}", b);
//panic!("Stop");
// Tmp to look at difference in CG and exact
let ignored_columns = match OPTIMISE_ENERGY_METHOD {
EnergyOptimisationMethod::ExactInverse => {
exact_overlap_inverse(&work_derivative, &mut b, EPSILON_SHIFT, NPARAMS as i32, PARAM_THRESHOLD)
Expand Down
206 changes: 206 additions & 0 deletions src/bin/dvmc_pairwf2.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,206 @@
//use blas::dcopy;
//use log::{debug, info};
use rand_mt::Mt64;
use rand::Rng;
//use indicatif::{ProgressBar, ProgressStyle};
use std::fs::File;
use std::io::Write;

use impurity::{VarParams, SysParams, generate_bitmask, FockState, RandomStateGeneration};
use impurity::dvmc::{variationnal_monte_carlo, EnergyOptimisationMethod, EnergyComputationMethod, VMCParams};

const SEED: u64 = 14;
const SIZE: usize = 8;
const NFIJ: usize = 4*SIZE*SIZE;
const NVIJ: usize = SIZE*(SIZE - 1) / 2;
const NGI: usize = SIZE;
const NPARAMS: usize = NFIJ + NGI + NVIJ;
const NELEC: usize = SIZE;
const NMCSAMP: usize = 1_000;
const NBOOTSTRAP: usize = 1;
const NMCWARMUP: usize = 1000;
const MCSAMPLE_INTERVAL: usize = 1;
const _NTHREADS: usize = 1;
const CLEAN_UPDATE_FREQUENCY: usize = 32;
const TOLERENCE_SHERMAN_MORRISSON: f64 = 1e-12;
const TOLERENCE_SINGULARITY: f64 = 1e-12;
const CONS_U: f64 = 1.0;
const CONS_T: f64 = 1.0;
const INITIAL_RATIO_UT: f64 = 1.0;
const FINAL_RATIO_UT: f64 = 32.0;
const NRATIO_POINTS: usize = 1;
const EPSILON_CG: f64 = 1e-16;
const EPSILON_SHIFT: f64 = 1e-2;
const OPTIMISATION_TIME_STEP: f64 = 1e-2;
const OPTIMISATION_DECAY: f64 = 0.0;
const NOPTITER: usize = 1000;
const KMAX: usize = NPARAMS;
const PARAM_THRESHOLD: f64 = <f64>::MIN;
//const PARAM_THRESHOLD: f64 = 0.0;
const OPTIMISE: bool = true;
const OPTIMISE_GUTZ: bool = true;
const OPTIMISE_JAST: bool = true;
const OPTIMISE_ORB: bool = true;
const SET_EXPVALO_ZERO: bool = false;
const COMPUTE_ENERGY_METHOD: EnergyComputationMethod = EnergyComputationMethod::MonteCarlo;
const OPTIMISE_ENERGY_METHOD: EnergyOptimisationMethod = EnergyOptimisationMethod::ConjugateGradiant;
const ENERGY_CONV_AVERAGE_SAMPLE: usize = 20;

// 4 et 2 Sites
//pub const HOPPINGS: [f64; SIZE*SIZE] = [
// //0.0, 1.0, 1.0, 0.0,
// 0.0, 1.0, 1.0, 0.0,
// 1.0, 0.0, 0.0, 1.0,
// 1.0, 0.0, 0.0, 1.0,
// 0.0, 1.0, 1.0, 0.0
//];

// 7 Sites
//pub const HOPPINGS: [f64; SIZE*SIZE] = [
// 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0,
// 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0,
// 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0,
// 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0,
// 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
// 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0,
// 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0,
//];

// 8 Sites
pub const HOPPINGS: [f64; SIZE*SIZE] = [
0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0,
1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0,
0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0,
1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0,
1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0,
0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0,
0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0,
0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0,
];

// 16 Sites
//pub const HOPPINGS: [f64; SIZE*SIZE] = [
// 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0,
// 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
// 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
// 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,
// 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
// 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
// 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
// 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0,
// 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0,
// 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0,
// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0,
// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0,
// 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0,
// 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0,
// 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0,
// 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0,
//];


fn sq(a: f64) -> f64 {
<f64>::abs(a) * <f64>::abs(a)
}

fn write_energy(fp: &mut File, e: &[f64]) {
let mut mean_energy = 0.0;
let mut mean_error = 0.0;
let mut mean_corr = 0.0;
for i in 0..ENERGY_CONV_AVERAGE_SAMPLE {
mean_energy += e[NOPTITER*3 - 3*i - 3];
mean_error += e[NOPTITER*3 - 3*i - 2];
mean_corr += e[NOPTITER*3 - 3*i - 1];
}
mean_energy *= 1.0 / ENERGY_CONV_AVERAGE_SAMPLE as f64;
mean_error *= 1.0 / ENERGY_CONV_AVERAGE_SAMPLE as f64;
mean_corr *= 1.0 / ENERGY_CONV_AVERAGE_SAMPLE as f64;

writeln!(fp, "{}", format!("{} {} {}", mean_energy, mean_error, mean_corr)).unwrap();
}

fn log_energy_convs(e: &[f64], fp: &mut File) {
for i in 0..NOPTITER {
writeln!(fp, "{}", format!("{} {} {}", e[3*i], e[3*i +1], e[3*i+2])).unwrap();
}
}

fn main() {
let mut fp = File::create("u_t_sweep").unwrap();
writeln!(fp, "{}", format!("# {} {} {}", SIZE, NMCSAMP, NOPTITER)).unwrap();
let mut paramsfp = File::create("params").unwrap();
writeln!(paramsfp, "{}", format!("# {} {} {}", SIZE, NMCSAMP, NOPTITER)).unwrap();
let mut _save: bool = true;
// Initialize logger
env_logger::init();
let bitmask = generate_bitmask(&HOPPINGS, SIZE);
let mut rng = Mt64::new(SEED);

for nsweep_iter in 0..NRATIO_POINTS {
let system_params = SysParams {
size: SIZE,
nelec: NELEC,
array_size: (SIZE + 7) / 8,
cons_t: -CONS_T,
cons_u: CONS_T * INITIAL_RATIO_UT + (nsweep_iter as f64 / NRATIO_POINTS as f64) * CONS_T * (FINAL_RATIO_UT - INITIAL_RATIO_UT),
nfij: NFIJ,
nvij: NVIJ,
ngi: NGI,
transfert_matrix: &HOPPINGS,
hopping_bitmask: &bitmask,
clean_update_frequency: CLEAN_UPDATE_FREQUENCY,
nmcsample: NMCSAMP,
nbootstrap: NBOOTSTRAP,
nmcwarmup: NMCWARMUP,
mcsample_interval: MCSAMPLE_INTERVAL,
tolerance_sherman_morrison: TOLERENCE_SHERMAN_MORRISSON,
tolerance_singularity: TOLERENCE_SINGULARITY,
pair_wavefunction: true,
};
println!("U = {}, T = {}", system_params.cons_u, system_params.cons_t);

let mut all_params: Vec<f64> = Vec::with_capacity(NGI + NVIJ + NFIJ);
for _ in 0..(NGI + NVIJ + NFIJ) {
all_params.push(rng.gen());
}
let (gi, params) = all_params.split_at_mut(NGI);
let (vij, fij) = params.split_at_mut(NVIJ);
let mut parameters = VarParams {
fij,
gi,
vij,
size: SIZE
};

let vmcparams = VMCParams {
dt: OPTIMISATION_TIME_STEP,
optimisation_decay: OPTIMISATION_DECAY,
threshold: PARAM_THRESHOLD,
kmax: KMAX,
epsilon: EPSILON_SHIFT,
epsilon_cg: EPSILON_CG,
noptiter: NOPTITER,
nparams: NGI + NVIJ + NFIJ,
optimise: OPTIMISE,
optimise_gutzwiller: OPTIMISE_GUTZ,
optimise_jastrow: OPTIMISE_JAST,
optimise_orbital: OPTIMISE_ORB,
compute_energy_method: COMPUTE_ENERGY_METHOD,
optimise_energy_method: OPTIMISE_ENERGY_METHOD,
};

let state: FockState<u8> = {
let mut tmp: FockState<u8> = FockState::generate_from_nelec(&mut rng, NELEC, SIZE);
while tmp.spin_up.count_ones() != tmp.spin_down.count_ones() {
tmp = FockState::generate_from_nelec(&mut rng, NELEC, SIZE);
}
tmp
};

let e_array = variationnal_monte_carlo(&mut rng, state, &mut parameters, &system_params, &vmcparams);
write_energy(&mut fp, &e_array);

log_energy_convs(&e_array, &mut paramsfp);

}
}
2 changes: 1 addition & 1 deletion src/constants.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
const SIZE: usize = 4;
const SIZE: usize = 17;
const ARRAY_SIZE: usize = (SIZE + 7) / 8;
/// System parameters
/// TODOC
Expand Down
Loading
Loading