From bc67c181d8bc000dfc513889f7c8713c7c4558cb Mon Sep 17 00:00:00 2001 From: Dimitri Bonanni-Surprenant Date: Mon, 7 Apr 2025 13:36:37 -0400 Subject: [PATCH 1/7] WIP --- src/dvmc.rs | 58 +++++---- src/lib.rs | 32 +++-- src/optimisation.rs | 299 +++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 356 insertions(+), 33 deletions(-) diff --git a/src/dvmc.rs b/src/dvmc.rs index f28d43d..d422f3b 100644 --- a/src/dvmc.rs +++ b/src/dvmc.rs @@ -62,30 +62,42 @@ pub fn variationnal_monte_carlo( where T: BitOps + From + Display + Debug, Standard: Distribution { let mut output_energy_array = vec![0.0; vmcparams.noptiter * 3]; - let mut otilde: Vec = vec![0.0; (4*sys.nfij + sys.nvij + sys.ngi) * (sys.nmcsample + 1)]; - //let mut work_otilde: Vec = vec![0.0; (sys.nfij + sys.nvij + sys.ngi) * (sys.nmcsample + 1)]; - let mut expvalo: Vec = vec![0.0; 4*sys.nfij + sys.nvij + sys.ngi]; - //let mut work_expvalo: Vec = vec![0.0; sys.nfij + sys.nvij + sys.ngi]; - let mut expval_ho: Vec = vec![0.0; 4*sys.nfij + sys.nvij + sys.ngi]; - //let mut work_expval_ho: Vec = vec![0.0; sys.nfij + sys.nvij + sys.ngi]; - let mut visited: Vec = vec![0; sys.nmcsample + 1]; - //let mut work_visited: Vec = vec![0; sys.nmcsample + 1]; - let mut der = DerivativeOperator { - o_tilde: &mut otilde, - expval_o: &mut expvalo, - ho: &mut expval_ho, - n: (4*sys.nfij + sys.nvij + sys.ngi) as i32, - nsamp: match vmcparams.compute_energy_method { + //let mut otilde: Vec = vec![0.0; (4*sys.nfij + sys.nvij + sys.ngi) * (sys.nmcsample + 1)]; + ////let mut work_otilde: Vec = vec![0.0; (sys.nfij + sys.nvij + sys.ngi) * (sys.nmcsample + 1)]; + //let mut expvalo: Vec = vec![0.0; 4*sys.nfij + sys.nvij + sys.ngi]; + ////let mut work_expvalo: Vec = vec![0.0; sys.nfij + sys.nvij + sys.ngi]; + //let mut expval_ho: Vec = vec![0.0; 4*sys.nfij + sys.nvij + sys.ngi]; + ////let mut work_expval_ho: Vec = vec![0.0; sys.nfij + sys.nvij + sys.ngi]; + //let mut visited: Vec = vec![0; sys.nmcsample + 1]; + ////let mut work_visited: Vec = vec![0; sys.nmcsample + 1]; + let mut der = DerivativeOperator::new( + (sys.ngi + sys.nvij + sys.nfij) as i32, + -1, + match vmcparams.compute_energy_method { EnergyComputationMethod::ExactSum => 1.0, EnergyComputationMethod::MonteCarlo => sys.nmcsample as f64, }, - nsamp_int: sys.mcsample_interval, - mu: -1, - visited: &mut visited, - pfaff_off: sys.ngi + sys.nvij, - jas_off: sys.ngi, - epsilon: vmcparams.epsilon, - }; + sys.mcsample_interval, + sys.ngi + sys.nvij, + sys.ngi, + vmcparams.epsilon + ); + //let mut der = DerivativeOperator { + // o_tilde: &mut otilde, + // expval_o: &mut expvalo, + // ho: &mut expval_ho, + // n: (4*sys.nfij + sys.nvij + sys.ngi) as i32, + // nsamp: match vmcparams.compute_energy_method { + // EnergyComputationMethod::ExactSum => 1.0, + // EnergyComputationMethod::MonteCarlo => sys.nmcsample as f64, + // }, + // nsamp_int: sys.mcsample_interval, + // mu: -1, + // visited: &mut visited, + // pfaff_off: sys.ngi + sys.nvij, + // jas_off: sys.ngi, + // epsilon: vmcparams.epsilon, + //}; //let mut work_derivative = DerivativeOperator { // o_tilde: &mut work_otilde, // expval_o: &mut work_expvalo, @@ -135,8 +147,8 @@ where T: BitOps + From + Display + Debug, Standard: Distribution unsafe { let incx = 1; let incy = 1; - daxpy(der.n, -mean_energy, der.expval_o, incx, der.ho, incy); - dcopy(der.n, der.ho, incx, &mut b, incy); + daxpy(der.n, -mean_energy, &der.expval_o, incx, &mut der.ho, incy); + dcopy(der.n, &der.ho, incx, &mut b, incy); } let mut _flag: bool = true; diff --git a/src/lib.rs b/src/lib.rs index a3620cf..78f15a3 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -51,23 +51,41 @@ impl<'a> std::fmt::Display for VarParams<'a> { /// 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], +pub struct DerivativeOperator { /// Number of variationnal parameters pub n: i32, /// Number of monte-carlo sampling pub mu: i32, pub nsamp: f64, pub nsamp_int: usize, - pub visited: &'a mut [usize], pub pfaff_off: usize, pub jas_off: usize, pub epsilon: f64, + pub o_tilde: Box<[f64]>, + pub expval_o: Box<[f64]>, + pub ho: Box<[f64]>, + pub visited: Box<[usize]>, } -impl<'a> std::fmt::Display for DerivativeOperator<'a> { +impl DerivativeOperator { + fn new(n: i32, mu: i32, nsamp: f64, nsamp_int: usize, pfaff_off: usize, jas_off: usize, epsilon: f64) -> Self { + DerivativeOperator { + o_tilde: vec![0.0; (n * mu) as usize].into_boxed_slice(), + expval_o: vec![0.0; n as usize].into_boxed_slice(), + ho: vec![0.0; n as usize].into_boxed_slice(), + visited: vec![0; nsamp as usize].into_boxed_slice(), + mu, + n, + nsamp, + nsamp_int, + pfaff_off, + jas_off, + epsilon + } + } +} + +impl<'a> std::fmt::Display for DerivativeOperator { fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::result::Result<(), std::fmt::Error> { let width = 10; let mut expval = " = ".to_owned(); @@ -99,7 +117,7 @@ fn _save_otilde(der: &DerivativeOperator) { let mut c = vec![0.0; (der.n * der.n) as usize]; println!("dim = {}", der.n * der.n); unsafe { - dgemm(b"N"[0], b"T"[0], der.n, der.n, der.mu, 1.0, der.o_tilde, der.n, der.o_tilde, der.n, 0.0, &mut c, der.n); + dgemm(b"N"[0], b"T"[0], der.n, der.n, der.mu, 1.0, &der.o_tilde, der.n, &der.o_tilde, der.n, 0.0, &mut c, der.n); } let mut outstr = "".to_owned(); outstr.push_str(&format!(" = ")); diff --git a/src/optimisation.rs b/src/optimisation.rs index 0520a80..d589875 100644 --- a/src/optimisation.rs +++ b/src/optimisation.rs @@ -5,6 +5,299 @@ use colored::Colorize; use crate::DerivativeOperator; +/// Encoding of a map to a reduced representation +/// # Map encoding +/// This struct encodes a Map to and from the general representation with +/// all the parameters to and from a reduced representation where some parameters +/// are multiples of other. +/// This struct is usefull to reduce the dimension of the optimisation and also +/// reduces the amount of linearly dependant columns, which helps with conditionning. +pub struct GenParameterMap { + /// Number of independant parameters. + pub dim: i32, + /// Dimension of general representation + pub gendim: i32, + /// Number of general parameters + pub n_genparams: i32, + /// Number of independant gutzwiller parameters. + pub n_independant_gutzwiller: usize, + /// Number of independant jastrow parameters. + pub n_independant_jastrow: usize, + /// The projector from the general representation to the reduced one. + /// Should behave: + /// Ax => \tilde{x} where every exquivalent component of x are sumed in + /// the lowest index and zero elsewhere. + /// Example: + /// A = ((1, 1), (0, 0)) + /// x = (a, b) + /// \tilde{x} = (a + b, 0) + /// A will be triangular, so only the upper triangle needs to be stored. + /// If you transpose this matrix, it will reverse this operation somewhat + /// A^T = ((1, 0), (1, 0) + /// A\tilde{x} = (a + b, a + b) + /// This indeed copies the parameter x_0 to x_1. If x_1 is set to half x_0, + /// then + /// A = ((1, 2), (0, 0)) + /// x = (a, b), Ax = (a + 2b) + /// Then the coefficient wise-inverse of the transpose is usefull + /// A^T^{-1} = ((1, 0), (1/2, 0)) + /// x = (a, 0) + /// \tilde{x} = (a, 1/2 a) + /// It is then usefull to store only the lower triangle of the "inverse" of + /// the transpose, as coefficient don't matter when projecting, but do + /// matter when recovering. + /// + /// Indexing + /// n = dim (dim - 1) / 2 + i - (dim - j) * ( dim - j - 1) / 2 + /// Can only have one non zero value per row + /// + /// PERFORMANCE: + /// If performance is an issue, an idea would be to write this matrix inside + /// an integer. This is not done here as it removes to possibility to make + /// a parameter say half another one. + pub projector: Box<[f64]>, +} + +pub trait ReducibleGeneralRepresentation { + /// Overrides b and x0 + /// Safety: + /// b and x0 must have lenght self.gendim + fn mapto_general_representation(&self, der: &DerivativeOperator, b: &mut [f64], x0: &mut [f64],) -> DerivativeOperator; + /// Overrides b and x0 + /// Safety: + /// b and x0 must have lenght self.gendim + /// On exit, the range [self.dim..self.gendim] contains garbage. + fn mapto_reduced_representation(&self, der: &DerivativeOperator, b: &mut [f64], x0: &mut [f64],) -> DerivativeOperator; + /// Overrides b and x0 + /// Safety: + /// b and x0 must have lenght self.gendim + fn update_general_representation(&self, der: &DerivativeOperator, out_der: &mut DerivativeOperator, b: &mut [f64], x0: &mut [f64],); + /// Overrides b and x0 + /// Safety: + /// b and x0 must have lenght self.gendim + /// On exit, the range [self.dim..self.gendim] contains garbage. + fn update_reduced_representation(&self, der: &DerivativeOperator, out_der: &mut DerivativeOperator, b: &mut [f64], x0: &mut [f64],); +} + +impl<'a> ReducibleGeneralRepresentation for GenParameterMap { + fn mapto_general_representation(&self, der: &DerivativeOperator, b: &mut [f64], x0: &mut [f64]) -> DerivativeOperator + { + if self.dim == 0 { + error!("Tried to map from an empty parameter representation, this will + break conjugate gradient or diagonalisation. Dim = {}", self.dim); + panic!("Undefined behavior."); + } + + let mut out_b = vec![0.0; self.gendim as usize]; + let mut out_x = vec![0.0; self.gendim as usize]; + let n = self.gendim as usize; + let mut out_der = DerivativeOperator::new(self.gendim, der.mu, der.nsamp, der.nsamp_int, + self.n_independant_gutzwiller + self.n_independant_jastrow, + self.n_independant_gutzwiller, der.epsilon); + + // Copy the visited vector + for i in 0..der.mu as usize { + out_der.visited[i] = der.visited[i]; + } + + // Copy all other data + for j in 0..n { + for i in 0..j { + let pij = self.projector[ + (n * (n - 1) / 2) + i - ((n - j) * (n - j - 1) /2) + ]; + if pij != 0.0 { + out_der.expval_o[j] = pij * der.expval_o[i]; + out_der.ho[j] = pij * der.ho[i]; + out_b[j] = pij * b[i]; + out_x[j] = pij * x0[i]; + unsafe { + let incx = self.dim; + let incy = self.gendim; + dcopy( + der.mu, + &der.o_tilde[j..(self.dim*der.mu) as usize], + incx, + &mut out_der.o_tilde[i..(self.gendim * der.mu) as usize], + incy + ); + } + // Rest of the row is zero as stated in doc + break; + } + } + } + unsafe { + let incx = 1; + let incy = 1; + dcopy(self.gendim, &out_b, incx, b, incy); + dcopy(self.gendim, &out_x, incx, x0, incy); + } + + // Should return b and x0 as well + out_der + } + + fn mapto_reduced_representation(&self, der: &DerivativeOperator, b: &mut [f64], x0: &mut [f64]) -> DerivativeOperator + { + if self.dim == 0 { + error!("Tried to map from an empty parameter representation, this will + break conjugate gradient or diagonalisation. Dim = {}", self.dim); + panic!("Undefined behavior."); + } + let n = self.gendim as usize; + let mut out_der = DerivativeOperator::new(self.dim, der.mu, der.nsamp, der.nsamp_int, + self.n_independant_gutzwiller + self.n_independant_jastrow, + self.n_independant_gutzwiller, der.epsilon); + + // Copy the visited vector + for i in 0..der.mu as usize { + out_der.visited[i] = der.visited[i]; + } + + // Copy all other data + // iter over transpose's diagonal + let mut i: usize = 0; + for j in 0..n { + let pjj = self.projector[ + (n * (n - 1) / 2) + j - ((n - j) * (n - j - 1) /2) + ]; + if pjj != 0.0 { // If 0, parameter does not map to itself. + if pjj != 1.0 { + error!("Parameters must map to themselfs with multiple 1."); + panic!("Undefined behavior."); + } + out_der.expval_o[i] = der.expval_o[j]; + out_der.ho[i] = der.ho[j]; + // Modify the vectors inplace, we won't iterate on past values + b[i] = b[j]; + x0[i] = x0[j]; + unsafe { + let incx = self.gendim; + let incy = self.dim; + dcopy( + der.mu, + &der.o_tilde[j..n*der.mu as usize], + incx, + &mut out_der.o_tilde[i..(self.dim * der.mu) as usize], + incy + ); + } + // Filling the independant parameter index + i += 1; + } + } + if i != self.dim as usize { + error!("Parameter map did not yield a dimension of {}, instead got {}", self.dim, i); + panic!("Undefined behavior."); + } + + out_der + } + + fn update_general_representation(&self, der: &DerivativeOperator, out_der: &mut DerivativeOperator, b: &mut [f64], x0: &mut [f64]) { + if self.dim == 0 { + error!("Tried to map from an empty parameter representation, this will + break conjugate gradient or diagonalisation. Dim = {}", self.dim); + panic!("Undefined behavior."); + } + + let mut out_b = vec![0.0; self.gendim as usize]; + let mut out_x = vec![0.0; self.gendim as usize]; + let n = self.gendim as usize; + + // Copy the visited vector + for i in 0..der.mu as usize { + out_der.visited[i] = der.visited[i]; + } + + // Copy all other data + for j in 0..n { + for i in 0..j { + let pij = self.projector[ + (n * (n - 1) / 2) + i - ((n - j) * (n - j - 1) /2) + ]; + if pij != 0.0 { + out_der.expval_o[j] = pij * der.expval_o[i]; + out_der.ho[j] = pij * der.ho[i]; + out_b[j] = pij * b[i]; + out_x[j] = pij * x0[i]; + unsafe { + let incx = self.dim; + let incy = self.gendim; + dcopy( + der.mu, + &der.o_tilde[j..(self.dim*der.mu) as usize], + incx, + &mut out_der.o_tilde[i..(self.gendim * der.mu) as usize], + incy + ); + } + // Rest of the row is zero as stated in doc + break; + } + } + } + unsafe { + let incx = 1; + let incy = 1; + dcopy(self.gendim, &out_b, incx, b, incy); + dcopy(self.gendim, &out_x, incx, x0, incy); + } + } + + fn update_reduced_representation(&self, der: &DerivativeOperator, out_der: &mut DerivativeOperator, b: &mut [f64], x0: &mut [f64]) { + if self.dim == 0 { + error!("Tried to map from an empty parameter representation, this will + break conjugate gradient or diagonalisation. Dim = {}", self.dim); + panic!("Undefined behavior."); + } + let n = self.gendim as usize; + + // Copy the visited vector + for i in 0..der.mu as usize { + out_der.visited[i] = der.visited[i]; + } + + // Copy all other data + // iter over transpose's diagonal + let mut i: usize = 0; + for j in 0..n { + let pjj = self.projector[ + (n * (n - 1) / 2) + j - ((n - j) * (n - j - 1) /2) + ]; + if pjj != 0.0 { // If 0, parameter does not map to itself. + if pjj != 1.0 { + error!("Parameters must map to themselfs with multiple 1."); + panic!("Undefined behavior."); + } + out_der.expval_o[i] = der.expval_o[j]; + out_der.ho[i] = der.ho[j]; + // Modify the vectors inplace, we won't iterate on past values + b[i] = b[j]; + x0[i] = x0[j]; + unsafe { + let incx = self.gendim; + let incy = self.dim; + dcopy( + der.mu, + &der.o_tilde[j..n*der.mu as usize], + incx, + &mut out_der.o_tilde[i..(self.dim * der.mu) as usize], + incy + ); + } + // Filling the independant parameter index + i += 1; + } + } + if i != self.dim as usize { + error!("Parameter map did not yield a dimension of {}, instead got {}", self.dim, i); + panic!("Undefined behavior."); + } + } +} + fn gradient(x: &[f64], otilde: &[f64], visited: &[usize], expval_o: &[f64], b: &mut [f64], epsilon: f64, dim: i32, mu: i32, nsamp: f64) { let alpha = -1.0; let incx = 1; @@ -308,7 +601,7 @@ pub fn exact_overlap_inverse(a: &DerivativeOperator, b: &mut [f64], epsilon: f64 // 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); + 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); @@ -412,7 +705,7 @@ pub fn conjugate_gradiant(a: &DerivativeOperator, b: &mut [f64], x0: &mut [f64], e *= epsilon_convergence; trace!("Error threshold e = {}", e); //println!("Error threshold e = {}", e); - gradient(x0, õ, a.visited, &expvalo, b, epsilon, new_dim as i32, a.mu, a.nsamp); + gradient(x0, õ, &a.visited, &expvalo, b, epsilon, new_dim as i32, a.mu, a.nsamp); let mut p = vec![0.0; new_dim]; let mut j: usize = 0; for i in 0..dim as usize { @@ -432,7 +725,7 @@ pub fn conjugate_gradiant(a: &DerivativeOperator, b: &mut [f64], x0: &mut [f64], trace!("p_{} : {:?}", k, p); //println!("r_{} : {:?}", k, b); //println!("p_{} : {:?}", k, p); - compute_w(&mut w, õ, a.visited, &expvalo, &p, epsilon, new_dim as i32, a.mu, a.nsamp); + compute_w(&mut w, õ, &a.visited, &expvalo, &p, epsilon, new_dim as i32, a.mu, a.nsamp); //println!("w_{} : {:?}", k, w); let nrm2rk = alpha_k(b, &p, &w, &mut alpha, new_dim as i32); trace!("alpha_{} : {}", k, alpha); From bc7b7718625af468130d0e0c7c78d81e8300b70a Mon Sep 17 00:00:00 2001 From: Dimitri Bonanni-Surprenant Date: Tue, 8 Apr 2025 09:20:40 -0400 Subject: [PATCH 2/7] (feat) General parameter map (untested) --- src/bin/dvmc_pairwf.rs | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/bin/dvmc_pairwf.rs b/src/bin/dvmc_pairwf.rs index 55df29f..ca7984b 100644 --- a/src/bin/dvmc_pairwf.rs +++ b/src/bin/dvmc_pairwf.rs @@ -87,7 +87,7 @@ fn _save_otilde(fp: &mut File, der: &DerivativeOperator) { let mut c = vec![0.0; (der.n * der.n) as usize]; println!("dim = {}", der.n * der.n); unsafe { - dgemm(b"N"[0], b"T"[0], der.n, der.n, der.mu, 1.0, der.o_tilde, der.n, der.o_tilde, der.n, 0.0, &mut c, der.n); + dgemm(b"N"[0], b"T"[0], der.n, der.n, der.mu, 1.0, &der.o_tilde, der.n, &der.o_tilde, der.n, 0.0, &mut c, der.n); } let mut outstr = "".to_owned(); outstr.push_str(&format!(" = ")); @@ -262,9 +262,9 @@ fn main() { let mut visited: Vec = vec![0; NMCSAMP + 1]; let mut work_visited: Vec = vec![0; NMCSAMP + 1]; let mut derivative = DerivativeOperator { - o_tilde: &mut otilde, - expval_o: &mut expvalo, - ho: &mut expval_ho, + o_tilde: otilde.into_boxed_slice(), + expval_o: expvalo.into_boxed_slice(), + ho: expval_ho.into_boxed_slice(), n: (4*NFIJ + NVIJ + NGI) as i32, nsamp: match COMPUTE_ENERGY_METHOD { EnergyComputationMethod::ExactSum => 1.0, @@ -272,15 +272,15 @@ fn main() { }, nsamp_int: MCSAMPLE_INTERVAL, mu: -1, - visited: &mut visited, + visited: visited.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: EPSILON_SHIFT, }; let mut work_derivative = DerivativeOperator { - o_tilde: &mut work_otilde, - expval_o: &mut work_expvalo, - ho: &mut work_expval_ho, + o_tilde: work_otilde.into_boxed_slice(), + expval_o: work_expvalo.into_boxed_slice(), + ho: work_expval_ho.into_boxed_slice(), n: (NFIJ + NVIJ + NGI) as i32, nsamp: match COMPUTE_ENERGY_METHOD { EnergyComputationMethod::ExactSum => 1.0, @@ -288,7 +288,7 @@ fn main() { }, nsamp_int: MCSAMPLE_INTERVAL, mu: -1, - visited: &mut work_visited, + visited: work_visited.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: EPSILON_SHIFT, @@ -347,8 +347,8 @@ fn main() { unsafe { let incx = 1; let incy = 1; - daxpy(work_derivative.n, -mean_energy, work_derivative.expval_o, incx, work_derivative.ho, incy); - dcopy(work_derivative.n, work_derivative.ho, incx, &mut b, incy); + daxpy(work_derivative.n, -mean_energy, &work_derivative.expval_o, incx, &mut work_derivative.ho, incy); + dcopy(work_derivative.n, &work_derivative.ho, incx, &mut b, incy); } //save_otilde(&mut der_fp, &derivative); //save_otilde(&mut wder_fp, &work_derivative); From e6eda192676b25da4ebe30939ba183d2f3606cc2 Mon Sep 17 00:00:00 2001 From: Dimitri Bonanni-Surprenant Date: Tue, 8 Apr 2025 10:08:08 -0400 Subject: [PATCH 3/7] (fix) Compilations errors in tests. --- src/bin/boboche2000.rs | 18 ++++--- src/bin/dvmc_fastupdate.rs | 23 +++++---- src/bin/dvmc_pairwf.rs | 22 +++++---- src/bin/dvmc_pairwf2.rs | 7 ++- src/constants.rs | 49 +++++++++++++++++++ src/dvmc.rs | 3 +- src/monte_carlo.rs | 8 +-- src/pfaffian.rs | 10 ++++ tests/orbitales_def/test_csv_parse_error.csv | 2 - tests/orbitales_def/test_int_parse_error.csv | 2 - tests/orbitales_def/test_parse_2x2.csv | 2 - tests/orbitales_def/test_parse_3x3.csv | 4 -- tests/parse_orbitale_test.rs | 51 -------------------- tests/vmc_2sites.rs | 22 +++++---- tests/vmc_2sites_analytic_conv.rs | 42 ++++++++-------- tests/vmc_2sites_exact_sum.rs | 38 ++++++++------- tests/vmc_clean_updates.rs | 18 ++++--- 17 files changed, 171 insertions(+), 150 deletions(-) delete mode 100644 tests/orbitales_def/test_csv_parse_error.csv delete mode 100644 tests/orbitales_def/test_int_parse_error.csv delete mode 100644 tests/orbitales_def/test_parse_2x2.csv delete mode 100644 tests/orbitales_def/test_parse_3x3.csv delete mode 100644 tests/parse_orbitale_test.rs diff --git a/src/bin/boboche2000.rs b/src/bin/boboche2000.rs index 14d1995..853e513 100644 --- a/src/bin/boboche2000.rs +++ b/src/bin/boboche2000.rs @@ -84,26 +84,28 @@ fn main() { hopping_bitmask: &bitmask, clean_update_frequency: CLEAN_UPDATE_FREQUENCY, nmcwarmup: NMCWARMUP, + nwarmupchains: 1, nmcsample: NMCSAMP, tolerance_sherman_morrison: TOLERENCE_SHERMAN_MORRISSON, tolerance_singularity: TOLERANCE_SINGULARITY, pair_wavefunction: false, + _opt_iter: 0, }; log_system_parameters(&sys); - 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 otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let visited: Vec = vec![0; NMCSAMP]; let mut der = DerivativeOperator { - o_tilde: &mut otilde, - expval_o: &mut expvalo, - ho: &mut expval_ho, + o_tilde: otilde.into_boxed_slice(), + expval_o: expvalo.into_boxed_slice(), + ho: expval_ho.into_boxed_slice(), n: (NFIJ + NVIJ + NGI) as i32, nsamp: NMCSAMP as f64, nsamp_int: 1, mu: -1, - visited: &mut visited, + visited: visited.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: 0.0, diff --git a/src/bin/dvmc_fastupdate.rs b/src/bin/dvmc_fastupdate.rs index f979ab5..5c29566 100644 --- a/src/bin/dvmc_fastupdate.rs +++ b/src/bin/dvmc_fastupdate.rs @@ -19,6 +19,7 @@ const NELEC: usize = SIZE; const NMCSAMP: usize = 10_000; const NBOOTSTRAP: usize = 1; const NMCWARMUP: usize = 1000; +const NWARMUPCHAINS: usize = 1; const MCSAMPLE_INTERVAL: usize = 1; const _NTHREADS: usize = 6; const CLEAN_UPDATE_FREQUENCY: usize = 2; @@ -160,10 +161,12 @@ fn main() { nmcsample: NMCSAMP, nbootstrap: NBOOTSTRAP, nmcwarmup: NMCWARMUP, + nwarmupchains: NWARMUPCHAINS, mcsample_interval: MCSAMPLE_INTERVAL, tolerance_sherman_morrison: TOLERENCE_SHERMAN_MORRISSON, tolerance_singularity: TOLERENCE_SINGULARITY, pair_wavefunction: false, + _opt_iter: 0, }; let mut rng = Mt64::new(SEED); @@ -270,19 +273,19 @@ 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 otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let visited: Vec = vec![0; NMCSAMP]; let mut derivative = DerivativeOperator { - o_tilde: &mut otilde, - expval_o: &mut expvalo, - ho: &mut expval_ho, + o_tilde: otilde.into_boxed_slice(), + expval_o: expvalo.into_boxed_slice(), + ho: expval_ho.into_boxed_slice(), n: (NFIJ + NVIJ + NGI) as i32, nsamp: NMCSAMP as f64, nsamp_int: MCSAMPLE_INTERVAL, mu: -1, - visited: &mut visited, + visited: visited.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: EPSILON_SPREAD, @@ -320,8 +323,8 @@ fn main() { unsafe { let incx = 1; let incy = 1; - daxpy(derivative.n, -mean_energy, derivative.expval_o, incx, derivative.ho, incy); - dcopy(derivative.n, derivative.ho, incx, &mut b, incy); + daxpy(derivative.n, -mean_energy, &derivative.expval_o, incx, &mut derivative.ho, incy); + dcopy(derivative.n, &derivative.ho, incx, &mut b, incy); } spread_eigenvalues(&mut derivative); exact_overlap_inverse(&derivative, &mut b, EPSILON_CG, NPARAMS as i32, PARAMTHRESHOLD); diff --git a/src/bin/dvmc_pairwf.rs b/src/bin/dvmc_pairwf.rs index ca7984b..89b07e2 100644 --- a/src/bin/dvmc_pairwf.rs +++ b/src/bin/dvmc_pairwf.rs @@ -19,6 +19,7 @@ const NELEC: usize = SIZE; const NMCSAMP: usize = 500; const NBOOTSTRAP: usize = 1; const NMCWARMUP: usize = 1000; +const NWARMUPCHAINS: usize = 1; const MCSAMPLE_INTERVAL: usize = 1; const _NTHREADS: usize = 1; const CLEAN_UPDATE_FREQUENCY: usize = 32; @@ -192,7 +193,7 @@ fn main() { // Initialize logger env_logger::init(); let bitmask = generate_bitmask(&HOPPINGS, SIZE); - let system_params = SysParams { + let mut system_params = SysParams { size: SIZE, nelec: NELEC, array_size: (SIZE + 7) / 8, @@ -207,10 +208,12 @@ fn main() { nmcsample: NMCSAMP, nbootstrap: NBOOTSTRAP, nmcwarmup: NMCWARMUP, + nwarmupchains: NWARMUPCHAINS, mcsample_interval: MCSAMPLE_INTERVAL, tolerance_sherman_morrison: TOLERENCE_SHERMAN_MORRISSON, tolerance_singularity: TOLERENCE_SINGULARITY, pair_wavefunction: true, + _opt_iter: 0, }; let mut rng = Mt64::new(SEED); @@ -253,14 +256,14 @@ fn main() { tmp }; - let mut otilde: Vec = vec![0.0; (4*NFIJ + NVIJ + NGI) * (NMCSAMP + 1)]; - let mut work_otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * (NMCSAMP + 1)]; - let mut expvalo: Vec = vec![0.0; 4*NFIJ + NVIJ + NGI]; - let mut work_expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut expval_ho: Vec = vec![0.0; 4*NFIJ + NVIJ + NGI]; - let mut work_expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut visited: Vec = vec![0; NMCSAMP + 1]; - let mut work_visited: Vec = vec![0; NMCSAMP + 1]; + let otilde: Vec = vec![0.0; (4*NFIJ + NVIJ + NGI) * (NMCSAMP + 1)]; + let work_otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * (NMCSAMP + 1)]; + let expvalo: Vec = vec![0.0; 4*NFIJ + NVIJ + NGI]; + let work_expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let expval_ho: Vec = vec![0.0; 4*NFIJ + NVIJ + NGI]; + let work_expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let visited: Vec = vec![0; NMCSAMP + 1]; + let work_visited: Vec = vec![0; NMCSAMP + 1]; let mut derivative = DerivativeOperator { o_tilde: otilde.into_boxed_slice(), expval_o: expvalo.into_boxed_slice(), @@ -305,6 +308,7 @@ fn main() { .progress_chars("##-")); for opt_iter in 0..NOPTITER { + system_params._opt_iter = opt_iter; let (mean_energy, accumulated_states, deltae, correlation_time) = { match COMPUTE_ENERGY_METHOD { EnergyComputationMethod::MonteCarlo => compute_mean_energy(&mut rng, state, ¶meters, &system_params, &mut derivative), diff --git a/src/bin/dvmc_pairwf2.rs b/src/bin/dvmc_pairwf2.rs index 426c69e..5827de9 100644 --- a/src/bin/dvmc_pairwf2.rs +++ b/src/bin/dvmc_pairwf2.rs @@ -19,6 +19,7 @@ const NELEC: usize = SIZE; const NMCSAMP: usize = 1_000; const NBOOTSTRAP: usize = 1; const NMCWARMUP: usize = 1000; +const NWARMUPCHAINS: usize = 1; const MCSAMPLE_INTERVAL: usize = 1; const _NTHREADS: usize = 1; const CLEAN_UPDATE_FREQUENCY: usize = 32; @@ -137,7 +138,7 @@ fn main() { let mut rng = Mt64::new(SEED); for nsweep_iter in 0..NRATIO_POINTS { - let system_params = SysParams { + let mut system_params = SysParams { size: SIZE, nelec: NELEC, array_size: (SIZE + 7) / 8, @@ -152,10 +153,12 @@ fn main() { nmcsample: NMCSAMP, nbootstrap: NBOOTSTRAP, nmcwarmup: NMCWARMUP, + nwarmupchains: NWARMUPCHAINS, mcsample_interval: MCSAMPLE_INTERVAL, tolerance_sherman_morrison: TOLERENCE_SHERMAN_MORRISSON, tolerance_singularity: TOLERENCE_SINGULARITY, pair_wavefunction: true, + _opt_iter: 0, }; println!("U = {}, T = {}", system_params.cons_u, system_params.cons_t); @@ -197,7 +200,7 @@ fn main() { tmp }; - let e_array = variationnal_monte_carlo(&mut rng, state, &mut parameters, &system_params, &vmcparams); + let e_array = variationnal_monte_carlo(&mut rng, state, &mut parameters, &mut system_params, &vmcparams); write_energy(&mut fp, &e_array); log_energy_convs(&e_array, &mut paramsfp); diff --git a/src/constants.rs b/src/constants.rs index 97fece5..5784503 100644 --- a/src/constants.rs +++ b/src/constants.rs @@ -16,12 +16,61 @@ pub struct SysParams<'a> { pub hopping_bitmask: &'a [SpinState], pub nmcwarmup: usize, pub nmcsample: usize, + pub nwarmupchains: usize, pub nbootstrap: usize, pub mcsample_interval: usize, pub clean_update_frequency: usize, pub tolerance_singularity: f64, pub tolerance_sherman_morrison: f64, pub pair_wavefunction: bool, + pub _opt_iter: usize, +} + +impl<'a> SysParams<'a> { + pub fn new( + size: usize, + nelec: usize, + array_size: usize, + cons_u: f64, + cons_t: f64, + nfij: usize, + nvij: usize, + ngi: usize, + transfert_matrix: &'a [f64], + hopping_bitmask: &'a [SpinState], + nmcwarmup: usize, + nmcsample: usize, + nwarmupchains: usize, + nbootstrap: usize, + mcsample_interval: usize, + clean_update_frequency: usize, + tolerance_singularity: f64, + tolerance_sherman_morrison: f64, + pair_wavefunction: bool, + ) -> Self { + SysParams { + size, + nelec, + array_size, + cons_t, + cons_u, + nfij, + nvij, + ngi, + transfert_matrix, + hopping_bitmask, + nmcwarmup, + nmcsample, + nwarmupchains, + nbootstrap, + mcsample_interval, + clean_update_frequency, + tolerance_sherman_morrison, + tolerance_singularity, + pair_wavefunction, + _opt_iter: 0, + } + } } pub type BitStruct = u8; diff --git a/src/dvmc.rs b/src/dvmc.rs index d422f3b..8f6cdc9 100644 --- a/src/dvmc.rs +++ b/src/dvmc.rs @@ -56,7 +56,7 @@ pub fn variationnal_monte_carlo( rng: &mut R, initial_state: FockState, params: &mut VarParams, - sys: &SysParams, + sys: &mut SysParams, vmcparams: &VMCParams ) -> Vec where T: BitOps + From + Display + Debug, Standard: Distribution @@ -115,6 +115,7 @@ where T: BitOps + From + Display + Debug, Standard: Distribution // epsilon: vmcparams.epsilon, //}; for opt_iter in 0..vmcparams.noptiter { + sys._opt_iter = opt_iter; let (mean_energy, _accumulated_states, deltae, correlation_time) = { match vmcparams.compute_energy_method { EnergyComputationMethod::MonteCarlo => compute_mean_energy(rng, initial_state, params, sys, &mut der), diff --git a/src/monte_carlo.rs b/src/monte_carlo.rs index 4ffcc83..5d02fa7 100644 --- a/src/monte_carlo.rs +++ b/src/monte_carlo.rs @@ -313,9 +313,11 @@ where Standard: Distribution let mut n_values = vec![0; error_estimation_level]; let mut energy_bootstraped = 0.0; - 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); + 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 diff --git a/src/pfaffian.rs b/src/pfaffian.rs index ed769a6..a4be29a 100644 --- a/src/pfaffian.rs +++ b/src/pfaffian.rs @@ -658,9 +658,11 @@ mod tests { clean_update_frequency: 0, nmcwarmup: 0, nmcsample: 0, + nwarmupchains: 0, tolerance_sherman_morrison: 0.0, tolerance_singularity: 0.0, pair_wavefunction: false, + _opt_iter: 0, }; let mut rng = SmallRng::seed_from_u64(42); // Size of the system @@ -797,9 +799,11 @@ mod tests { clean_update_frequency: 0, nmcwarmup: 0, nmcsample: 0, + nwarmupchains: 0, tolerance_sherman_morrison: 0.0, tolerance_singularity: 0.0, pair_wavefunction: false, + _opt_iter: 0, }; let mut params = vec![0.0; 4 * SIZE * SIZE]; // params[i+8*j] = f_ij @@ -840,9 +844,11 @@ mod tests { clean_update_frequency: 0, nmcwarmup: 0, nmcsample: 0, + nwarmupchains: 0, tolerance_sherman_morrison: 0.0, tolerance_singularity: 0.0, pair_wavefunction: false, + _opt_iter: 0, }; let mut params = vec![0.0; 4 * SIZE * SIZE]; // params[i+8*j] = f_ij @@ -902,9 +908,11 @@ mod tests { clean_update_frequency: 0, nmcwarmup: 0, nmcsample: 0, + nwarmupchains: 0, tolerance_sherman_morrison: 0.0, tolerance_singularity: 0.0, pair_wavefunction: false, + _opt_iter: 0, }; let mut params = vec![0.0; 4 * SIZE * SIZE]; // params[i+8*j] = f_ij @@ -982,9 +990,11 @@ mod tests { clean_update_frequency: 0, nmcwarmup: 0, nmcsample: 0, + nwarmupchains: 0, tolerance_sherman_morrison: 0.0, tolerance_singularity: 0.0, pair_wavefunction: false, + _opt_iter: 0, }; let mut params = vec![0.0; 4 * SIZE * SIZE]; // params[i+8*j] = f_ij diff --git a/tests/orbitales_def/test_csv_parse_error.csv b/tests/orbitales_def/test_csv_parse_error.csv deleted file mode 100644 index 96e3e80..0000000 --- a/tests/orbitales_def/test_csv_parse_error.csv +++ /dev/null @@ -1,2 +0,0 @@ -i,j,n -0,1,1,1 diff --git a/tests/orbitales_def/test_int_parse_error.csv b/tests/orbitales_def/test_int_parse_error.csv deleted file mode 100644 index f24926c..0000000 --- a/tests/orbitales_def/test_int_parse_error.csv +++ /dev/null @@ -1,2 +0,0 @@ -i,j,n -0,1,k diff --git a/tests/orbitales_def/test_parse_2x2.csv b/tests/orbitales_def/test_parse_2x2.csv deleted file mode 100644 index eb2b5e3..0000000 --- a/tests/orbitales_def/test_parse_2x2.csv +++ /dev/null @@ -1,2 +0,0 @@ -i,j,n -0,1,1 diff --git a/tests/orbitales_def/test_parse_3x3.csv b/tests/orbitales_def/test_parse_3x3.csv deleted file mode 100644 index ca3e18a..0000000 --- a/tests/orbitales_def/test_parse_3x3.csv +++ /dev/null @@ -1,4 +0,0 @@ -i,j,n -0,1,1 -0,2,2 -1,2,3 diff --git a/tests/parse_orbitale_test.rs b/tests/parse_orbitale_test.rs deleted file mode 100644 index 9f77418..0000000 --- a/tests/parse_orbitale_test.rs +++ /dev/null @@ -1,51 +0,0 @@ -use std::path::Path; - -use impurity::parse::orbitale::parse_orbitale_def; - -#[test] -#[ignore] -fn test_parse_2x2() { - let size = 2; - let fp = Path::new("tests/orbitales_def/test_parse_2x2.csv"); - let fij = parse_orbitale_def(&fp.to_path_buf(), size).unwrap(); - assert_eq!(fij.len(), 4); - assert_eq!(fij[1], -fij[2]); - assert!(fij[1] <= 0.0); -} - -#[test] -#[ignore] -fn test_parse_3x3() { - let size = 3; - let fp = Path::new("tests/orbitales_def/test_parse_3x3.csv"); - let fij = parse_orbitale_def(&fp.to_path_buf(), size).unwrap(); - assert_eq!(fij.len(), 9); - assert_eq!(fij[1], -fij[3]); - assert_eq!(fij[6], -fij[2]); - assert_eq!(fij[5], -fij[7]); - assert!(fij[1] <= 0.0); - assert!(fij[2] <= 0.0); - assert!(fij[5] <= 0.0); -} - -#[test] -fn test_int_parse_error() { - let size = 2; - let fp = Path::new("tests/orbitales_def/test_int_parse_error.csv"); - let fij = parse_orbitale_def(&fp.to_path_buf(), size); - match fij { - Ok(_) => panic!("Should have errored."), - Err(_) => println!("Error as expected."), - } -} - -#[test] -fn test_csv_parse_error() { - let size = 2; - let fp = Path::new("tests/orbitales_def/test_csv_parse_error.csv"); - let fij = parse_orbitale_def(&fp.to_path_buf(), size); - match fij { - Ok(_) => panic!("Should have errored."), - Err(_) => println!("Error as expected."), - } -} diff --git a/tests/vmc_2sites.rs b/tests/vmc_2sites.rs index 66e6be8..0eb4024 100644 --- a/tests/vmc_2sites.rs +++ b/tests/vmc_2sites.rs @@ -357,9 +357,11 @@ fn comupte_energy_from_all_states() { clean_update_frequency: CLEAN_UPDATE_FREQUENCY, nmcwarmup: NMCWARMUP, nmcsample: NMCSAMP, + nwarmupchains: 1, tolerance_sherman_morrison: TOL_SHERMAN, tolerance_singularity: TOL_SINGULARITY, pair_wavefunction: false, + _opt_iter: 0, }; let mut fij = [ 0.0, 0.0, 0.0, 0.0, @@ -423,19 +425,19 @@ fn comupte_energy_from_all_states() { energy_individual_state(&states_names[i], ¶meters), 1e-12); } close(mean_energy, analytic(¶meters), 1e-12); - 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 otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let visited: Vec = vec![0; NMCSAMP]; let mut der = DerivativeOperator { - o_tilde: &mut otilde, - expval_o: &mut expvalo, - ho: &mut expval_ho, + o_tilde: otilde.into_boxed_slice(), + expval_o: expvalo.into_boxed_slice(), + ho: expval_ho.into_boxed_slice(), n: (NFIJ + NVIJ + NGI) as i32, nsamp: NMCSAMP as f64, nsamp_int: 1, mu: -1, - visited: &mut visited, + visited: visited.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: 0.0, @@ -466,7 +468,7 @@ fn comupte_energy_from_all_states() { // Test derivatives let exp_val = analytic_derivatives_expval(¶meters); - print_der(der.expval_o, &exp_val, sys.ngi+sys.nvij+sys.nfij); + print_der(&der.expval_o, &exp_val, sys.ngi+sys.nvij+sys.nfij); let psi = norm(¶meters); println!("Norm: {:10.4e}", psi); for i in 0..sys.ngi+sys.nvij+sys.nfij { @@ -474,7 +476,7 @@ fn comupte_energy_from_all_states() { } let exp_val_ho = analytic_ho_expval(¶meters); - print_der(der.ho, &exp_val_ho, sys.ngi+sys.nvij+sys.nfij); + print_der(&der.ho, &exp_val_ho, sys.ngi+sys.nvij+sys.nfij); let psi = norm(¶meters); println!("Norm: {:10.4e}", psi); for i in 0..sys.ngi+sys.nvij+sys.nfij { diff --git a/tests/vmc_2sites_analytic_conv.rs b/tests/vmc_2sites_analytic_conv.rs index 7ea254f..460f10d 100644 --- a/tests/vmc_2sites_analytic_conv.rs +++ b/tests/vmc_2sites_analytic_conv.rs @@ -326,9 +326,11 @@ fn comupte_energy_from_all_states() { clean_update_frequency: CLEAN_UPDATE_FREQUENCY, nmcwarmup: NMCWARMUP, nmcsample: NMCSAMP, + nwarmupchains: 1, tolerance_sherman_morrison: TOL_SHERMAN, tolerance_singularity: TOL_SINGULARITY, pair_wavefunction: false, + _opt_iter: 0, }; let mut fij = [ 0.0, 0.0, 0.0, 0.0, @@ -393,36 +395,36 @@ fn comupte_energy_from_all_states() { } close(mean_energy, analytic(¶meters), 1e-12); mean_energy = mean_energy / norm(¶meters); - 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 otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let visited: Vec = vec![0; NMCSAMP]; let mut der = DerivativeOperator { - o_tilde: &mut otilde, - expval_o: &mut expvalo, - ho: &mut expval_ho, + o_tilde: otilde.into_boxed_slice(), + expval_o: expvalo.into_boxed_slice(), + ho: expval_ho.into_boxed_slice(), n: (NFIJ + NVIJ + NGI) as i32, nsamp: NMCSAMP as f64, nsamp_int: 1, mu: -1, - visited: &mut visited, + visited: visited.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: 0.0, }; - let mut otilde_pair: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; - let mut expvalo_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut expval_ho_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut visited_pair: Vec = vec![0; NMCSAMP]; + let otilde_pair: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let expvalo_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let expval_ho_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let visited_pair: Vec = vec![0; NMCSAMP]; let mut der_pair = DerivativeOperator { - o_tilde: &mut otilde_pair, - expval_o: &mut expvalo_pair, - ho: &mut expval_ho_pair, + o_tilde: otilde_pair.into_boxed_slice(), + expval_o: expvalo_pair.into_boxed_slice(), + ho: expval_ho_pair.into_boxed_slice(), n: (sys.size*sys.size + NVIJ + NGI) as i32, nsamp: NMCSAMP as f64, nsamp_int: 1, mu: -1, - visited: &mut visited_pair, + visited: visited_pair.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: 0.0, @@ -448,7 +450,7 @@ fn comupte_energy_from_all_states() { // Test derivatives let exp_val = analytic_derivatives_expval(¶meters); println!("Checking "); - print_der(der_pair.expval_o, &exp_val, sys.ngi+sys.nvij+sys.nfij); + print_der(&der_pair.expval_o, &exp_val, sys.ngi+sys.nvij+sys.nfij); let psi = norm(¶meters); println!("Norm: {:10.4e}", psi); for i in 0..sys.ngi+sys.nvij+sys.nfij { @@ -458,7 +460,7 @@ fn comupte_energy_from_all_states() { let exp_val_ho = analytic_ho_expval(¶meters); println!("Checking "); - print_der(der_pair.ho, &exp_val_ho, sys.ngi+sys.nvij+sys.nfij); + print_der(&der_pair.ho, &exp_val_ho, sys.ngi+sys.nvij+sys.nfij); let psi = norm(¶meters); println!("Norm: {:10.4e}", psi); for i in 0..sys.ngi+sys.nvij+sys.nfij { @@ -478,8 +480,8 @@ fn comupte_energy_from_all_states() { unsafe { let incx = 1; let incy = 1; - daxpy(der_pair.n, -mean_energy, der_pair.expval_o, incx, der_pair.ho, incy); - dcopy(der_pair.n, der_pair.ho, incx, &mut b, incy); + daxpy(der_pair.n, -mean_energy, &der_pair.expval_o, incx, &mut der_pair.ho, incy); + dcopy(der_pair.n, &der_pair.ho, incx, &mut b, incy); } println!("x0 = {:?}", x0); } diff --git a/tests/vmc_2sites_exact_sum.rs b/tests/vmc_2sites_exact_sum.rs index 048d67a..29eb7b6 100644 --- a/tests/vmc_2sites_exact_sum.rs +++ b/tests/vmc_2sites_exact_sum.rs @@ -322,9 +322,11 @@ fn comupte_energy_from_all_states() { clean_update_frequency: CLEAN_UPDATE_FREQUENCY, nmcwarmup: NMCWARMUP, nmcsample: NMCSAMP, + nwarmupchains: 1, tolerance_sherman_morrison: TOL_SHERMAN, tolerance_singularity: TOL_SINGULARITY, pair_wavefunction: false, + _opt_iter: 0, }; let mut fij = [ 0.0, 0.0, 0.0, 0.0, @@ -389,36 +391,36 @@ fn comupte_energy_from_all_states() { } close(mean_energy, analytic(¶meters), 1e-12); mean_energy = mean_energy / norm(¶meters); - 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 otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let visited: Vec = vec![0; NMCSAMP]; let mut der = DerivativeOperator { - o_tilde: &mut otilde, - expval_o: &mut expvalo, - ho: &mut expval_ho, + o_tilde: otilde.into_boxed_slice(), + expval_o: expvalo.into_boxed_slice(), + ho: expval_ho.into_boxed_slice(), n: (NFIJ + NVIJ + NGI) as i32, nsamp: NMCSAMP as f64, nsamp_int: 1, mu: -1, - visited: &mut visited, + visited: visited.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: 0.0, }; - let mut otilde_pair: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; - let mut expvalo_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut expval_ho_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut visited_pair: Vec = vec![0; NMCSAMP]; + let otilde_pair: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let expvalo_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let expval_ho_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let visited_pair: Vec = vec![0; NMCSAMP]; let mut der_pair = DerivativeOperator { - o_tilde: &mut otilde_pair, - expval_o: &mut expvalo_pair, - ho: &mut expval_ho_pair, + o_tilde: otilde_pair.into_boxed_slice(), + expval_o: expvalo_pair.into_boxed_slice(), + ho: expval_ho_pair.into_boxed_slice(), n: (sys.size*sys.size + NVIJ + NGI) as i32, nsamp: NMCSAMP as f64, nsamp_int: 1, mu: -1, - visited: &mut visited_pair, + visited: visited_pair.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: 0.0, @@ -444,7 +446,7 @@ fn comupte_energy_from_all_states() { // Test derivatives let exp_val = analytic_derivatives_expval(¶meters); println!("Checking "); - print_der(der_pair.expval_o, &exp_val, sys.ngi+sys.nvij+sys.nfij); + print_der(&der_pair.expval_o, &exp_val, sys.ngi+sys.nvij+sys.nfij); let psi = norm(¶meters); println!("Norm: {:10.4e}", psi); for i in 0..sys.ngi+sys.nvij+sys.nfij { @@ -454,7 +456,7 @@ fn comupte_energy_from_all_states() { let exp_val_ho = analytic_ho_expval(¶meters); println!("Checking "); - print_der(der_pair.ho, &exp_val_ho, sys.ngi+sys.nvij+sys.nfij); + print_der(&der_pair.ho, &exp_val_ho, sys.ngi+sys.nvij+sys.nfij); let psi = norm(¶meters); println!("Norm: {:10.4e}", psi); for i in 0..sys.ngi+sys.nvij+sys.nfij { diff --git a/tests/vmc_clean_updates.rs b/tests/vmc_clean_updates.rs index 855c7b2..4e67ab4 100644 --- a/tests/vmc_clean_updates.rs +++ b/tests/vmc_clean_updates.rs @@ -86,25 +86,27 @@ fn monte_carlo_first_iteration() { clean_update_frequency: CLEAN_UPDATE_FREQUENCY, nmcwarmup: NMCWARMUP, nmcsample: NMCSAMP, + nwarmupchains: 1, tolerance_sherman_morrison: TOLERENCE_SHERMAN_MORRISSON, tolerance_singularity: TOLERANCE_SINGULARITY, pair_wavefunction: false, + _opt_iter: 0, }; log_system_parameters(&sys); - 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 otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let visited: Vec = vec![0; NMCSAMP]; let mut der = DerivativeOperator { - o_tilde: &mut otilde, - expval_o: &mut expvalo, - ho: &mut expval_ho, + o_tilde: otilde.into_boxed_slice(), + expval_o: expvalo.into_boxed_slice(), + ho: expval_ho.into_boxed_slice(), n: (NFIJ + NVIJ + NGI) as i32, nsamp: NMCSAMP as f64, nsamp_int: 1, mu: -1, - visited: &mut visited, + visited: visited.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: 0.0, From b95d600ef85fceb92a157f40095286d9a450504a Mon Sep 17 00:00:00 2001 From: Dimitri Bonanni-Surprenant Date: Sun, 20 Apr 2025 11:13:58 -0400 Subject: [PATCH 4/7] (feat) General parameter map (tested for N=2), works until nparams = 3 --- src/bin/dvmc_pairwf.rs | 2 +- src/bin/dvmc_pairwf2.rs | 117 +++++++++++++++++++++++++++++-------- src/dvmc.rs | 49 +++++++++------- src/lib.rs | 2 +- src/optimisation.rs | 125 ++++++++++++++++++++++++++++------------ src/parse/mod.rs | 4 -- src/parse/orbitale.rs | 116 ------------------------------------- 7 files changed, 213 insertions(+), 202 deletions(-) delete mode 100644 src/parse/mod.rs delete mode 100644 src/parse/orbitale.rs diff --git a/src/bin/dvmc_pairwf.rs b/src/bin/dvmc_pairwf.rs index 89b07e2..0f5bf6b 100644 --- a/src/bin/dvmc_pairwf.rs +++ b/src/bin/dvmc_pairwf.rs @@ -467,5 +467,5 @@ fn main() { opt_progress_bar.inc(1); opt_progress_bar.set_message(format!("Changed parameters by norm: {:+>.05e} Current energy: {:+>.05e}", opt_delta, mean_energy)); } - opt_progress_bar.finish() + opt_progress_bar.finish(); } diff --git a/src/bin/dvmc_pairwf2.rs b/src/bin/dvmc_pairwf2.rs index 5827de9..afc7a9b 100644 --- a/src/bin/dvmc_pairwf2.rs +++ b/src/bin/dvmc_pairwf2.rs @@ -1,3 +1,4 @@ +use impurity::optimisation::GenParameterMap; //use blas::dcopy; //use log::{debug, info}; use rand_mt::Mt64; @@ -10,7 +11,7 @@ use impurity::{VarParams, SysParams, generate_bitmask, FockState, RandomStateGen use impurity::dvmc::{variationnal_monte_carlo, EnergyOptimisationMethod, EnergyComputationMethod, VMCParams}; const SEED: u64 = 14; -const SIZE: usize = 8; +const SIZE: usize = 2; const NFIJ: usize = 4*SIZE*SIZE; const NVIJ: usize = SIZE*(SIZE - 1) / 2; const NGI: usize = SIZE; @@ -27,7 +28,7 @@ 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 INITIAL_RATIO_UT: f64 = 8.0; const FINAL_RATIO_UT: f64 = 32.0; const NRATIO_POINTS: usize = 1; const EPSILON_CG: f64 = 1e-16; @@ -36,7 +37,7 @@ const OPTIMISATION_TIME_STEP: f64 = 1e-2; const OPTIMISATION_DECAY: f64 = 0.0; const NOPTITER: usize = 1000; const KMAX: usize = NPARAMS; -const PARAM_THRESHOLD: f64 = ::MIN; +const PARAM_THRESHOLD: f64 = 1e-5; //const PARAM_THRESHOLD: f64 = 0.0; const OPTIMISE: bool = true; const OPTIMISE_GUTZ: bool = true; @@ -46,15 +47,19 @@ 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; +const N_INDEP_PARAMS: usize = 3; +const N_GUTZ: usize = NGI; +const N_JAST: usize = NVIJ; +const PAIRWF: bool = false; // 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 -//]; +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] = [ @@ -68,16 +73,16 @@ const ENERGY_CONV_AVERAGE_SAMPLE: usize = 20; //]; // 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, -]; +//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] = [ @@ -98,6 +103,53 @@ pub const HOPPINGS: [f64; SIZE*SIZE] = [ // 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, //]; +// + +// Pairwf +const PARAMS_PROJECTOR: [f64; (NFIJ + NVIJ + NGI) * (NFIJ + NVIJ + NGI - 1) / 2 + NFIJ + NVIJ + NGI] = [ + /* g0 */ 1.0, // 0, + /* g1 */ 1.0, 0.0, // + /* v01 */ 0.0, 0.0, 0.0, + /* f00uu */ 0.0, 0.0, 0.0, 0.0, + /* f01uu */ 0.0, 0.0, 0.0, 0.0, 0.0, + /* f10uu */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* f11uu */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* f00ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, + /* f01ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, + /* f10ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, + /* f11ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, + /* f00du */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + /* f01du */ 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, + /* f10du */ 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, + /* f11du */ 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, + /* f00dd */ 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, + /* f01dd */ 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, + /* f10dd */ 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, + /* f11dd */ 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, +]; + +// General rep +//const PARAMS_PROJECTOR: [f64; (NFIJ + NVIJ + NGI) * (NFIJ + NVIJ + NGI - 1) / 2 + NFIJ + NVIJ + NGI] = [ +// /* g0 */ 1.0, // 0, +// /* g1 */ 0.0, 1.0, // +// /* v01 */ 0.0, 0.0, 1.0, +// /* f00uu */ 0.0, 0.0, 0.0, 1.0, +// /* f01uu */ 0.0, 0.0, 0.0, 0.0, 1.0, +// /* f10uu */ 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, +// /* f11uu */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, +// /* f00ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, +// /* f01ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, +// /* f10ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, +// /* f11ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, +// /* f00du */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, +// /* f01du */ 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, +// /* f10du */ 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, +// /* f11du */ 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, +// /* f00dd */ 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, +// /* f01dd */ 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, +// /* f10dd */ 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, 1.0, +// /* f11dd */ 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, 1.0, +//]; fn sq(a: f64) -> f64 { @@ -157,7 +209,7 @@ fn main() { mcsample_interval: MCSAMPLE_INTERVAL, tolerance_sherman_morrison: TOLERENCE_SHERMAN_MORRISSON, tolerance_singularity: TOLERENCE_SINGULARITY, - pair_wavefunction: true, + pair_wavefunction: PAIRWF, _opt_iter: 0, }; println!("U = {}, T = {}", system_params.cons_u, system_params.cons_t); @@ -174,6 +226,15 @@ fn main() { vij, size: SIZE }; + for i in 0..SIZE*SIZE { + parameters.fij[i] = 0.0; + parameters.fij[i +SIZE*SIZE] = 0.5; + parameters.fij[i +2*SIZE*SIZE] = 0.0; + parameters.fij[i +3*SIZE*SIZE] = 0.0; + } + parameters.gi[1] = 0.0; + parameters.gi[0] = parameters.gi[1]; + println!("{:?}", parameters.fij); let vmcparams = VMCParams { dt: OPTIMISATION_TIME_STEP, @@ -183,7 +244,7 @@ fn main() { epsilon: EPSILON_SHIFT, epsilon_cg: EPSILON_CG, noptiter: NOPTITER, - nparams: NGI + NVIJ + NFIJ, + nparams: N_INDEP_PARAMS, optimise: OPTIMISE, optimise_gutzwiller: OPTIMISE_GUTZ, optimise_jastrow: OPTIMISE_JAST, @@ -199,9 +260,19 @@ fn main() { } tmp }; + let param_map = GenParameterMap { + dim: N_INDEP_PARAMS as i32, + gendim: (NFIJ + NGI + NVIJ) as i32, + n_genparams: (NFIJ + NGI + NVIJ) as i32, + n_independant_gutzwiller: N_GUTZ, + n_independant_jastrow: N_JAST, + projector: Box::new(PARAMS_PROJECTOR), + }; - let e_array = variationnal_monte_carlo(&mut rng, state, &mut parameters, &mut system_params, &vmcparams); + let e_array = variationnal_monte_carlo(&mut rng, state, &mut parameters, &mut system_params, &vmcparams, ¶m_map); write_energy(&mut fp, &e_array); + println!("gi = {:?}", parameters.gi); + println!("fij = {:?}", parameters.fij); log_energy_convs(&e_array, &mut paramsfp); diff --git a/src/dvmc.rs b/src/dvmc.rs index 8f6cdc9..8aac065 100644 --- a/src/dvmc.rs +++ b/src/dvmc.rs @@ -6,7 +6,7 @@ use blas::{idamax, daxpy, dcopy, dscal}; use crate::{VarParams, DerivativeOperator, SysParams, FockState, BitOps}; use crate::monte_carlo::{compute_mean_energy, compute_mean_energy_exact}; -use crate::optimisation::{conjugate_gradiant, exact_overlap_inverse}; +use crate::optimisation::{conjugate_gradiant, exact_overlap_inverse, GenParameterMap, ReducibleGeneralRepresentation}; #[derive(Debug)] pub enum EnergyOptimisationMethod { @@ -39,10 +39,10 @@ pub struct VMCParams { } fn zero_out_derivatives(der: &mut DerivativeOperator, sys: &SysParams) { - for i in 0.. (sys.nfij + sys.nvij + sys.ngi) * sys.nmcsample { + for i in 0.. (der.n as usize) * sys.nmcsample { der.o_tilde[i] = 0.0; } - for i in 0..sys.nfij + sys.nvij + sys.ngi { + for i in 0..der.n as usize { der.expval_o[i] = 0.0; der.ho[i] = 0.0; } @@ -57,31 +57,28 @@ pub fn variationnal_monte_carlo( initial_state: FockState, params: &mut VarParams, sys: &mut SysParams, - vmcparams: &VMCParams + vmcparams: &VMCParams, + param_map: &GenParameterMap, ) -> Vec where T: BitOps + From + Display + Debug, Standard: Distribution { let mut output_energy_array = vec![0.0; vmcparams.noptiter * 3]; - //let mut otilde: Vec = vec![0.0; (4*sys.nfij + sys.nvij + sys.ngi) * (sys.nmcsample + 1)]; - ////let mut work_otilde: Vec = vec![0.0; (sys.nfij + sys.nvij + sys.ngi) * (sys.nmcsample + 1)]; - //let mut expvalo: Vec = vec![0.0; 4*sys.nfij + sys.nvij + sys.ngi]; - ////let mut work_expvalo: Vec = vec![0.0; sys.nfij + sys.nvij + sys.ngi]; - //let mut expval_ho: Vec = vec![0.0; 4*sys.nfij + sys.nvij + sys.ngi]; - ////let mut work_expval_ho: Vec = vec![0.0; sys.nfij + sys.nvij + sys.ngi]; - //let mut visited: Vec = vec![0; sys.nmcsample + 1]; - ////let mut work_visited: Vec = vec![0; sys.nmcsample + 1]; + + let mut x0 = vec![0.0; sys.ngi + sys.nvij + sys.nfij].into_boxed_slice(); + let mut b = vec![0.0; sys.ngi + sys.nvij + sys.nfij].into_boxed_slice(); let mut der = DerivativeOperator::new( - (sys.ngi + sys.nvij + sys.nfij) as i32, + vmcparams.nparams as i32, -1, match vmcparams.compute_energy_method { EnergyComputationMethod::ExactSum => 1.0, EnergyComputationMethod::MonteCarlo => sys.nmcsample as f64, }, sys.mcsample_interval, - sys.ngi + sys.nvij, - sys.ngi, + param_map.n_independant_jastrow + param_map.n_independant_gutzwiller, + param_map.n_independant_gutzwiller, vmcparams.epsilon ); + let mut work_der = param_map.mapto_general_representation(&der, &mut x0); //let mut der = DerivativeOperator { // o_tilde: &mut otilde, // expval_o: &mut expvalo, @@ -115,12 +112,15 @@ where T: BitOps + From + Display + Debug, Standard: Distribution // epsilon: vmcparams.epsilon, //}; for opt_iter in 0..vmcparams.noptiter { + //println!("{:?}", params.fij); + //println!("x0 = {:?}", x0); + //println!("b = {:?}", b); sys._opt_iter = opt_iter; let (mean_energy, _accumulated_states, deltae, correlation_time) = { match vmcparams.compute_energy_method { - EnergyComputationMethod::MonteCarlo => compute_mean_energy(rng, initial_state, params, sys, &mut der), + EnergyComputationMethod::MonteCarlo => compute_mean_energy(rng, initial_state, params, sys, &mut work_der), EnergyComputationMethod::ExactSum => { - (compute_mean_energy_exact(params, sys, &mut der), Vec::with_capacity(0), 0.0, 0.0) + (compute_mean_energy_exact(params, sys, &mut work_der), Vec::with_capacity(0), 0.0, 0.0) }, } }; @@ -138,13 +138,16 @@ where T: BitOps + From + Display + Debug, Standard: Distribution //} //mapto_pairwf(&derivative, &mut work_derivative, sys); - let mut x0 = vec![0.0; sys.nfij + sys.nvij + sys.ngi]; x0[(sys.ngi + sys.nvij)..(sys.ngi + sys.nvij + sys.nfij)].copy_from_slice(params.fij); x0[sys.ngi..(sys.ngi + sys.nvij)].copy_from_slice(params.vij); x0[0..sys.ngi].copy_from_slice(params.gi); + param_map.update_reduced_representation(&work_der, &mut der, &mut x0); + //println!("Expval Full = {:?}", work_der.ho); + //println!("Expval Reduced = {:?}", der.ho); + //panic!("Stop."); // 68 misawa - let mut b: Vec = vec![0.0; der.n as usize]; + //let mut b: Vec = vec![0.0; der.n as usize]; unsafe { let incx = 1; let incy = 1; @@ -161,7 +164,7 @@ where T: BitOps + From + Display + Debug, Standard: Distribution conjugate_gradiant(&der, &mut b, &mut x0, vmcparams.epsilon, vmcparams.kmax, vmcparams.nparams as i32, vmcparams.threshold, vmcparams.epsilon_cg) }, }; - let mut delta_alpha = vec![0.0; vmcparams.nparams]; + let mut delta_alpha = vec![0.0; param_map.gendim as usize]; let mut j: usize = 0; for i in 0..vmcparams.nparams { if ignored_columns[i] { @@ -173,6 +176,11 @@ where T: BitOps + From + Display + Debug, Standard: Distribution _flag = false; } } + //panic!("Stop!"); + param_map.update_delta_alpha_reduced_to_gen(&mut delta_alpha); + //println!("b = {:?}", b); + //println!("delta_alpha = {:?}", delta_alpha); + //panic!("Stop."); if vmcparams.optimise { unsafe { let incx = 1; @@ -241,6 +249,7 @@ where T: BitOps + From + Display + Debug, Standard: Distribution // ); //} zero_out_derivatives(&mut der, sys); + zero_out_derivatives(&mut work_der, sys); //print_delta_alpha(&delta_alpha, sys.ngi, sys.nvij, sys.nfij); //let opt_delta = unsafe { // let incx = 1; diff --git a/src/lib.rs b/src/lib.rs index 78f15a3..faa762d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -70,7 +70,7 @@ pub struct DerivativeOperator { impl DerivativeOperator { fn new(n: i32, mu: i32, nsamp: f64, nsamp_int: usize, pfaff_off: usize, jas_off: usize, epsilon: f64) -> Self { DerivativeOperator { - o_tilde: vec![0.0; (n * mu) as usize].into_boxed_slice(), + o_tilde: vec![0.0; n as usize * nsamp as usize].into_boxed_slice(), expval_o: vec![0.0; n as usize].into_boxed_slice(), ho: vec![0.0; n as usize].into_boxed_slice(), visited: vec![0; nsamp as usize].into_boxed_slice(), diff --git a/src/optimisation.rs b/src/optimisation.rs index d589875..0084d7d 100644 --- a/src/optimisation.rs +++ b/src/optimisation.rs @@ -62,25 +62,26 @@ pub trait ReducibleGeneralRepresentation { /// Overrides b and x0 /// Safety: /// b and x0 must have lenght self.gendim - fn mapto_general_representation(&self, der: &DerivativeOperator, b: &mut [f64], x0: &mut [f64],) -> DerivativeOperator; + fn mapto_general_representation(&self, der: &DerivativeOperator, x0: &mut [f64],) -> DerivativeOperator; /// Overrides b and x0 /// Safety: /// b and x0 must have lenght self.gendim /// On exit, the range [self.dim..self.gendim] contains garbage. - fn mapto_reduced_representation(&self, der: &DerivativeOperator, b: &mut [f64], x0: &mut [f64],) -> DerivativeOperator; + fn mapto_reduced_representation(&self, der: &DerivativeOperator, x0: &mut [f64],) -> DerivativeOperator; /// Overrides b and x0 /// Safety: /// b and x0 must have lenght self.gendim - fn update_general_representation(&self, der: &DerivativeOperator, out_der: &mut DerivativeOperator, b: &mut [f64], x0: &mut [f64],); + fn update_general_representation(&self, der: &DerivativeOperator, out_der: &mut DerivativeOperator, x0: &mut [f64],); /// Overrides b and x0 /// Safety: /// b and x0 must have lenght self.gendim /// On exit, the range [self.dim..self.gendim] contains garbage. - fn update_reduced_representation(&self, der: &DerivativeOperator, out_der: &mut DerivativeOperator, b: &mut [f64], x0: &mut [f64],); + fn update_reduced_representation(&self, der: &DerivativeOperator, out_der: &mut DerivativeOperator, x0: &mut [f64],); + fn update_delta_alpha_reduced_to_gen(&self, x0: &mut [f64]); } impl<'a> ReducibleGeneralRepresentation for GenParameterMap { - fn mapto_general_representation(&self, der: &DerivativeOperator, b: &mut [f64], x0: &mut [f64]) -> DerivativeOperator + fn mapto_general_representation(&self, der: &DerivativeOperator, x0: &mut [f64]) -> DerivativeOperator { if self.dim == 0 { error!("Tried to map from an empty parameter representation, this will @@ -88,12 +89,14 @@ impl<'a> ReducibleGeneralRepresentation for GenParameterMap { panic!("Undefined behavior."); } - let mut out_b = vec![0.0; self.gendim as usize]; let mut out_x = vec![0.0; self.gendim as usize]; let n = self.gendim as usize; let mut out_der = DerivativeOperator::new(self.gendim, der.mu, der.nsamp, der.nsamp_int, self.n_independant_gutzwiller + self.n_independant_jastrow, self.n_independant_gutzwiller, der.epsilon); + if der.mu < 0 { + return out_der; + } // Copy the visited vector for i in 0..der.mu as usize { @@ -101,36 +104,40 @@ impl<'a> ReducibleGeneralRepresentation for GenParameterMap { } // Copy all other data + let mut k = 0; + let mut flag = false; for j in 0..n { - for i in 0..j { + for i in n..n { let pij = self.projector[ - (n * (n - 1) / 2) + i - ((n - j) * (n - j - 1) /2) + j + (i * (i+1) / 2) ]; if pij != 0.0 { - out_der.expval_o[j] = pij * der.expval_o[i]; - out_der.ho[j] = pij * der.ho[i]; - out_b[j] = pij * b[i]; - out_x[j] = pij * x0[i]; + out_der.expval_o[i] = pij * der.expval_o[k]; + out_der.ho[i] = pij * der.ho[k]; + out_x[i] = pij * x0[k]; unsafe { let incx = self.dim; let incy = self.gendim; dcopy( der.mu, - &der.o_tilde[j..(self.dim*der.mu) as usize], + &der.o_tilde[i..(self.dim*der.mu) as usize], incx, - &mut out_der.o_tilde[i..(self.gendim * der.mu) as usize], + &mut out_der.o_tilde[k..(self.gendim * der.mu) as usize], incy ); } - // Rest of the row is zero as stated in doc - break; + flag = true; + //break; } } + if flag{ + k += 1; + flag = false; + } } unsafe { let incx = 1; let incy = 1; - dcopy(self.gendim, &out_b, incx, b, incy); dcopy(self.gendim, &out_x, incx, x0, incy); } @@ -138,7 +145,7 @@ impl<'a> ReducibleGeneralRepresentation for GenParameterMap { out_der } - fn mapto_reduced_representation(&self, der: &DerivativeOperator, b: &mut [f64], x0: &mut [f64]) -> DerivativeOperator + fn mapto_reduced_representation(&self, der: &DerivativeOperator, x0: &mut [f64]) -> DerivativeOperator { if self.dim == 0 { error!("Tried to map from an empty parameter representation, this will @@ -149,6 +156,9 @@ impl<'a> ReducibleGeneralRepresentation for GenParameterMap { let mut out_der = DerivativeOperator::new(self.dim, der.mu, der.nsamp, der.nsamp_int, self.n_independant_gutzwiller + self.n_independant_jastrow, self.n_independant_gutzwiller, der.epsilon); + if der.mu < 0 { + return out_der; + } // Copy the visited vector for i in 0..der.mu as usize { @@ -160,7 +170,7 @@ impl<'a> ReducibleGeneralRepresentation for GenParameterMap { let mut i: usize = 0; for j in 0..n { let pjj = self.projector[ - (n * (n - 1) / 2) + j - ((n - j) * (n - j - 1) /2) + j + (j * (j+1) / 2) ]; if pjj != 0.0 { // If 0, parameter does not map to itself. if pjj != 1.0 { @@ -170,7 +180,6 @@ impl<'a> ReducibleGeneralRepresentation for GenParameterMap { out_der.expval_o[i] = der.expval_o[j]; out_der.ho[i] = der.ho[j]; // Modify the vectors inplace, we won't iterate on past values - b[i] = b[j]; x0[i] = x0[j]; unsafe { let incx = self.gendim; @@ -195,16 +204,16 @@ impl<'a> ReducibleGeneralRepresentation for GenParameterMap { out_der } - fn update_general_representation(&self, der: &DerivativeOperator, out_der: &mut DerivativeOperator, b: &mut [f64], x0: &mut [f64]) { + fn update_general_representation(&self, der: &DerivativeOperator, out_der: &mut DerivativeOperator, x0: &mut [f64]) { if self.dim == 0 { error!("Tried to map from an empty parameter representation, this will break conjugate gradient or diagonalisation. Dim = {}", self.dim); panic!("Undefined behavior."); } - let mut out_b = vec![0.0; self.gendim as usize]; let mut out_x = vec![0.0; self.gendim as usize]; let n = self.gendim as usize; + out_der.mu = der.mu; // Copy the visited vector for i in 0..der.mu as usize { @@ -212,47 +221,90 @@ impl<'a> ReducibleGeneralRepresentation for GenParameterMap { } // Copy all other data + let mut k = 0; + let mut flag = false; for j in 0..n { - for i in 0..j { + for i in j..n { let pij = self.projector[ - (n * (n - 1) / 2) + i - ((n - j) * (n - j - 1) /2) + j + (i * (i+1) / 2) ]; if pij != 0.0 { - out_der.expval_o[j] = pij * der.expval_o[i]; - out_der.ho[j] = pij * der.ho[i]; - out_b[j] = pij * b[i]; - out_x[j] = pij * x0[i]; + out_der.expval_o[i] = pij * der.expval_o[k]; + out_der.ho[i] = pij * der.ho[k]; + out_x[i] = pij * x0[k]; unsafe { let incx = self.dim; let incy = self.gendim; dcopy( der.mu, - &der.o_tilde[j..(self.dim*der.mu) as usize], + &der.o_tilde[i..(self.dim*der.mu) as usize], incx, - &mut out_der.o_tilde[i..(self.gendim * der.mu) as usize], + &mut out_der.o_tilde[k..(self.gendim * der.mu) as usize], incy ); } - // Rest of the row is zero as stated in doc - break; + flag = true; + //break; } } + if flag{ + k += 1; + flag = false; + } } unsafe { let incx = 1; let incy = 1; - dcopy(self.gendim, &out_b, incx, b, incy); dcopy(self.gendim, &out_x, incx, x0, incy); } } - fn update_reduced_representation(&self, der: &DerivativeOperator, out_der: &mut DerivativeOperator, b: &mut [f64], x0: &mut [f64]) { + fn update_delta_alpha_reduced_to_gen(&self, x0: &mut [f64]) { if self.dim == 0 { error!("Tried to map from an empty parameter representation, this will break conjugate gradient or diagonalisation. Dim = {}", self.dim); panic!("Undefined behavior."); } let n = self.gendim as usize; + let mut out_x0 = vec![0.0; n]; + println!("da = {:?}", x0); + + // Copy all other data + // iter over transpose's diagonal + // Copy all other data + let mut k = 0; + let mut flag = false; + for j in 0..n { + for i in j..n { + let pij = self.projector[ + j + (i * (i+1) / 2) + ]; + if pij != 0.0 { + out_x0[i] = pij * x0[k]; + flag = true; + //break; + } + } + if flag{ + k += 1; + flag = false; + } + } + unsafe { + let incx = 1; + let incy = 1; + dcopy(self.gendim, &out_x0, incx, x0, incy); + } + } + + fn update_reduced_representation(&self, der: &DerivativeOperator, out_der: &mut DerivativeOperator, x0: &mut [f64]) { + if self.dim == 0 { + error!("Tried to map from an empty parameter representation, this will + break conjugate gradient or diagonalisation. Dim = {}", self.dim); + panic!("Undefined behavior."); + } + let n = self.gendim as usize; + out_der.mu = der.mu; // Copy the visited vector for i in 0..der.mu as usize { @@ -261,10 +313,10 @@ impl<'a> ReducibleGeneralRepresentation for GenParameterMap { // Copy all other data // iter over transpose's diagonal - let mut i: usize = 0; + let mut i = 0; for j in 0..n { let pjj = self.projector[ - (n * (n - 1) / 2) + j - ((n - j) * (n - j - 1) /2) + j + (j * (j+1) / 2) ]; if pjj != 0.0 { // If 0, parameter does not map to itself. if pjj != 1.0 { @@ -274,7 +326,6 @@ impl<'a> ReducibleGeneralRepresentation for GenParameterMap { out_der.expval_o[i] = der.expval_o[j]; out_der.ho[i] = der.ho[j]; // Modify the vectors inplace, we won't iterate on past values - b[i] = b[j]; x0[i] = x0[j]; unsafe { let incx = self.gendim; diff --git a/src/parse/mod.rs b/src/parse/mod.rs deleted file mode 100644 index 229ff7e..0000000 --- a/src/parse/mod.rs +++ /dev/null @@ -1,4 +0,0 @@ -pub use orbitale::parse_orbitale_def as orbitale_def; - -/// For orbitale.csv -pub mod orbitale; diff --git a/src/parse/orbitale.rs b/src/parse/orbitale.rs deleted file mode 100644 index 3185813..0000000 --- a/src/parse/orbitale.rs +++ /dev/null @@ -1,116 +0,0 @@ -use csv::{Reader, StringRecord}; -use derive_more::{Constructor, Error}; -use rand::Rng; -use std::fmt; -use std::fs::File; -use std::path::PathBuf; - -/// Parse the orbital variationnal definition -/// # Arguments -/// * __`fp`__ - File path to the definition file, in csv format. -/// * __`s`__ - Size of the system. -pub fn parse_orbitale_def(fp: &PathBuf, s: usize) -> Result> { - // Parameter reference vector. column packed. - let mut fij: Vec = vec![0.0; 4 * s * s]; - // Parameter vector. Ordered by "n". - let mut orbitale_params: Vec = Vec::new(); - let mut n_orbitale_params: usize = 0; - // RNG thread - let mut rng = rand::thread_rng(); - - // Parse input file. - let file = File::open(fp)?; - let mut reader = Reader::from_reader(file); - for (k, result) in reader.records().enumerate() { - let rec = result?; - // Should have 3 column. - if rec.len() != 3 { - let d = format!( - "Error at line {} in orbitale.csv, invalid number of elements.", - k - ); - let mut e = OrbitaleParseError::new("Invalid number of argument on a line.".to_owned()); - e.details = d; - return Err(e); - } - - // Parse element coordinates. - let i = parse_single_elem(&rec, 0, k)?; - let j = parse_single_elem(&rec, 1, k)?; - - let n = rec.get(2).unwrap().parse::(); - match n { - Ok(param) => { - if n_orbitale_params < param { - n_orbitale_params += 1; - orbitale_params.push(rng.gen::() * 2.0 - 1.0); - } - let var_param: f64 = orbitale_params[param - 1]; - fij[i + j * s] = var_param; - fij[j + i * s] = -var_param; - fij[i + j * s + s*s] = var_param; - fij[j + i * s + s*s] = -var_param; - fij[i + j * s + 2*s*s] = var_param; - fij[j + i * s + 2*s*s] = -var_param; - fij[i + j * s + 3*s*s] = var_param; - fij[j + i * s + 3*s*s] = -var_param; - } - Err(error) => { - // Add error message. - let d = format!("Error at line {}, invalid parameter identifier.", k); - let mut e = OrbitaleParseError::from(error); - e.details = d; - return Err(e); - } - }; - } - Ok(fij) -} - -fn parse_single_elem(line: &StringRecord, col: usize, l: usize) -> Result { - match line.get(col).unwrap().parse::() { - Ok(v) => Ok(v), - Err(error) => { - let details = format!( - "Expected valid coordinates in orbitale.csv at line {}, col {}", - l, col - ); - let mut e = OrbitaleParseError::from(error); - e.details = details; - Err(e) - } - } -} - -type Result = std::result::Result; - -/// Error in the orbital params definition. -#[derive(Debug, Clone, Error, Constructor)] -pub struct OrbitaleParseError { - pub details: String, -} - -impl fmt::Display for OrbitaleParseError { - fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - write!(f, "Parsing error encountered in orbitale.csv.") - } -} - -impl From for OrbitaleParseError { - fn from(err: std::num::ParseIntError) -> Self { - println!("{}", err); - OrbitaleParseError::new("Expected to parse an integer.".to_owned()) - } -} - -impl From for OrbitaleParseError { - fn from(err: csv::Error) -> Self { - OrbitaleParseError::new(err.to_string()) - } -} - -impl From for OrbitaleParseError { - fn from(err: std::io::Error) -> Self { - OrbitaleParseError::new(err.to_string()) - } -} From 1e55bffb11f264550a1331f8d3509b99f4465c1b Mon Sep 17 00:00:00 2001 From: Dimitri Bonanni-Surprenant Date: Tue, 22 Apr 2025 13:00:46 -0400 Subject: [PATCH 5/7] (fix) Any number of sites possible --- src/bin/dvmc_pairwf2.rs | 211 ++++++++++++++++++++++++---------------- src/constants.rs | 2 +- src/dvmc.rs | 20 ++-- src/fock_state.rs | 40 ++++++-- src/gutzwiller.rs | 4 +- src/hamiltonian.rs | 2 +- src/hoppings.rs | 2 + src/jastrow.rs | 6 +- src/optimisation.rs | 1 - 9 files changed, 178 insertions(+), 110 deletions(-) diff --git a/src/bin/dvmc_pairwf2.rs b/src/bin/dvmc_pairwf2.rs index afc7a9b..86b8941 100644 --- a/src/bin/dvmc_pairwf2.rs +++ b/src/bin/dvmc_pairwf2.rs @@ -10,14 +10,16 @@ use std::io::Write; use impurity::{VarParams, SysParams, generate_bitmask, FockState, RandomStateGeneration}; use impurity::dvmc::{variationnal_monte_carlo, EnergyOptimisationMethod, EnergyComputationMethod, VMCParams}; +type BitSize = u32; + const SEED: u64 = 14; -const SIZE: usize = 2; +const SIZE: usize = 16; const NFIJ: usize = 4*SIZE*SIZE; const NVIJ: usize = SIZE*(SIZE - 1) / 2; const NGI: usize = SIZE; const NPARAMS: usize = NFIJ + NGI + NVIJ; const NELEC: usize = SIZE; -const NMCSAMP: usize = 1_000; +const NMCSAMP: usize = 1000; const NBOOTSTRAP: usize = 1; const NMCWARMUP: usize = 1000; const NWARMUPCHAINS: usize = 1; @@ -32,12 +34,12 @@ 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-2; +const EPSILON_SHIFT: f64 = 1e-3; 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 = 1e-5; +const PARAM_THRESHOLD: f64 = 1e-3; //const PARAM_THRESHOLD: f64 = 0.0; const OPTIMISE: bool = true; const OPTIMISE_GUTZ: bool = true; @@ -47,19 +49,20 @@ 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; -const N_INDEP_PARAMS: usize = 3; +const N_INDEP_PARAMS: usize = SIZE*SIZE + NGI + NVIJ; const N_GUTZ: usize = NGI; const N_JAST: usize = NVIJ; const PAIRWF: bool = false; +const CONV_PARAM_THRESHOLD: f64 = 1e-10; // 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 -]; +//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] = [ @@ -84,73 +87,118 @@ pub const HOPPINGS: [f64; SIZE*SIZE] = [ // 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, //]; -// 16 Sites +// 9 Sites //pub const HOPPINGS: [f64; SIZE*SIZE] = [ -// 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, -// 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, -// 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, -// 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, -// 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, -// 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, -// 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, -// 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, -// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, -// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, -// 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, -// 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, -// 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, -// 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, +// 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, +// 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.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, 0.0, 0.0, 1.0, 0.0, +// 1.0, 0.0, 0.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, +// 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.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, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //]; -// -// Pairwf -const PARAMS_PROJECTOR: [f64; (NFIJ + NVIJ + NGI) * (NFIJ + NVIJ + NGI - 1) / 2 + NFIJ + NVIJ + NGI] = [ - /* g0 */ 1.0, // 0, - /* g1 */ 1.0, 0.0, // - /* v01 */ 0.0, 0.0, 0.0, - /* f00uu */ 0.0, 0.0, 0.0, 0.0, - /* f01uu */ 0.0, 0.0, 0.0, 0.0, 0.0, - /* f10uu */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* f11uu */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* f00ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, - /* f01ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, - /* f10ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, - /* f11ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, - /* f00du */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - /* f01du */ 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, - /* f10du */ 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, - /* f11du */ 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, - /* f00dd */ 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, - /* f01dd */ 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, - /* f10dd */ 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, - /* f11dd */ 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, +// 10 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, +// 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.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, 0.0, 0.0, 1.0, 0.0, 0.0, +// 1.0, 0.0, 0.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, +// 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.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, 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, +//]; + +// 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, ]; +// -// General rep +// Pairwf //const PARAMS_PROJECTOR: [f64; (NFIJ + NVIJ + NGI) * (NFIJ + NVIJ + NGI - 1) / 2 + NFIJ + NVIJ + NGI] = [ // /* g0 */ 1.0, // 0, -// /* g1 */ 0.0, 1.0, // -// /* v01 */ 0.0, 0.0, 1.0, -// /* f00uu */ 0.0, 0.0, 0.0, 1.0, -// /* f01uu */ 0.0, 0.0, 0.0, 0.0, 1.0, -// /* f10uu */ 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, -// /* f11uu */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, +// /* g1 */ 1.0, 0.0, // +// /* v01 */ 0.0, 0.0, 0.0, +// /* f00uu */ 0.0, 0.0, 0.0, 0.0, +// /* f01uu */ 0.0, 0.0, 0.0, 0.0, 0.0, +// /* f10uu */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* f11uu */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // /* f00ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, // /* f01ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, -// /* f10ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, -// /* f11ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, -// /* f00du */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, -// /* f01du */ 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, -// /* f10du */ 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, -// /* f11du */ 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, -// /* f00dd */ 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, -// /* f01dd */ 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, -// /* f10dd */ 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, 1.0, -// /* f11dd */ 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, 1.0, +// /* f10ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, +// /* f11ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, +// /* f00du */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* f01du */ 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, +// /* f10du */ 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, +// /* f11du */ 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, +// /* f00dd */ 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, +// /* f01dd */ 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, +// /* f10dd */ 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, +// /* f11dd */ 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, //]; +// 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; + 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) @@ -212,7 +260,7 @@ fn main() { pair_wavefunction: PAIRWF, _opt_iter: 0, }; - println!("U = {}, T = {}", system_params.cons_u, system_params.cons_t); + //println!("U = {}, T = {}", system_params.cons_u, system_params.cons_t); let mut all_params: Vec = Vec::with_capacity(NGI + NVIJ + NFIJ); for _ in 0..(NGI + NVIJ + NFIJ) { @@ -226,15 +274,13 @@ fn main() { vij, size: SIZE }; - for i in 0..SIZE*SIZE { - parameters.fij[i] = 0.0; - parameters.fij[i +SIZE*SIZE] = 0.5; - parameters.fij[i +2*SIZE*SIZE] = 0.0; - parameters.fij[i +3*SIZE*SIZE] = 0.0; - } - parameters.gi[1] = 0.0; - parameters.gi[0] = parameters.gi[1]; - println!("{:?}", parameters.fij); + //for i in 0..SIZE*SIZE { + // parameters.fij[i] = 0.0; + // parameters.fij[i +SIZE*SIZE] = 0.5; + // parameters.fij[i +2*SIZE*SIZE] = 0.0; + // parameters.fij[i +3*SIZE*SIZE] = 0.0; + //} + //println!("{:?}", parameters.fij); let vmcparams = VMCParams { dt: OPTIMISATION_TIME_STEP, @@ -251,10 +297,11 @@ fn main() { optimise_orbital: OPTIMISE_ORB, compute_energy_method: COMPUTE_ENERGY_METHOD, optimise_energy_method: OPTIMISE_ENERGY_METHOD, + conv_param_threshold: CONV_PARAM_THRESHOLD, }; - let state: FockState = { - let mut tmp: FockState = FockState::generate_from_nelec(&mut rng, NELEC, SIZE); + 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); } @@ -270,9 +317,7 @@ fn main() { }; let e_array = variationnal_monte_carlo(&mut rng, state, &mut parameters, &mut system_params, &vmcparams, ¶m_map); - write_energy(&mut fp, &e_array); - println!("gi = {:?}", parameters.gi); - println!("fij = {:?}", parameters.fij); + //write_energy(&mut fp, &e_array); log_energy_convs(&e_array, &mut paramsfp); diff --git a/src/constants.rs b/src/constants.rs index 5784503..6acf1d2 100644 --- a/src/constants.rs +++ b/src/constants.rs @@ -1,4 +1,4 @@ -const SIZE: usize = 17; +const SIZE: usize = 32; const ARRAY_SIZE: usize = (SIZE + 7) / 8; /// System parameters /// TODOC diff --git a/src/dvmc.rs b/src/dvmc.rs index 8aac065..8c686f3 100644 --- a/src/dvmc.rs +++ b/src/dvmc.rs @@ -2,7 +2,7 @@ use log::{info, error}; use rand::distributions::{Distribution, Standard}; use rand::Rng; use std::fmt::{Display, Debug}; -use blas::{idamax, daxpy, dcopy, dscal}; +use blas::{idamax, daxpy, dcopy, dscal, dnrm2}; use crate::{VarParams, DerivativeOperator, SysParams, FockState, BitOps}; use crate::monte_carlo::{compute_mean_energy, compute_mean_energy_exact}; @@ -36,6 +36,7 @@ pub struct VMCParams { pub optimise_gutzwiller: bool, pub optimise_jastrow: bool, pub optimise_orbital: bool, + pub conv_param_threshold: f64, } fn zero_out_derivatives(der: &mut DerivativeOperator, sys: &SysParams) { @@ -112,9 +113,6 @@ where T: BitOps + From + Display + Debug, Standard: Distribution // epsilon: vmcparams.epsilon, //}; for opt_iter in 0..vmcparams.noptiter { - //println!("{:?}", params.fij); - //println!("x0 = {:?}", x0); - //println!("b = {:?}", b); sys._opt_iter = opt_iter; let (mean_energy, _accumulated_states, deltae, correlation_time) = { match vmcparams.compute_energy_method { @@ -142,9 +140,6 @@ where T: BitOps + From + Display + Debug, Standard: Distribution x0[sys.ngi..(sys.ngi + sys.nvij)].copy_from_slice(params.vij); x0[0..sys.ngi].copy_from_slice(params.gi); param_map.update_reduced_representation(&work_der, &mut der, &mut x0); - //println!("Expval Full = {:?}", work_der.ho); - //println!("Expval Reduced = {:?}", der.ho); - //panic!("Stop."); // 68 misawa //let mut b: Vec = vec![0.0; der.n as usize]; @@ -178,9 +173,6 @@ where T: BitOps + From + Display + Debug, Standard: Distribution } //panic!("Stop!"); param_map.update_delta_alpha_reduced_to_gen(&mut delta_alpha); - //println!("b = {:?}", b); - //println!("delta_alpha = {:?}", delta_alpha); - //panic!("Stop."); if vmcparams.optimise { unsafe { let incx = 1; @@ -229,6 +221,14 @@ where T: BitOps + From + Display + Debug, Standard: Distribution } // HARD CODE vij = vji // Slater Rescaling + let opt_delta = unsafe { + let incx = 1; + dnrm2(der.n, &delta_alpha, incx) + }; + if opt_delta <= vmcparams.conv_param_threshold { + println!("Exit early, achieved convergence within {} iteration, update now under supplied threshold.", opt_iter+1); + return output_energy_array; + } unsafe { let incx = 1; let max = params.fij[idamax(sys.nfij as i32, ¶ms.fij, incx) - 1]; diff --git a/src/fock_state.rs b/src/fock_state.rs index 3602c1d..ad58e4e 100644 --- a/src/fock_state.rs +++ b/src/fock_state.rs @@ -127,6 +127,18 @@ impl From for u8 { } } +impl From for u16 { + fn from(other: SpinState) -> Self { + const NBYTES: usize = 2; + if ARRAY_SIZE > NBYTES {error!("Casting SpinState of size {} to discards overflowing bits.", SIZE)} + let mut tmp: [u8; NBYTES] = [0x00; NBYTES]; + for i in 0..ARRAY_SIZE { + tmp[i] = other.state[i]; + } + ::from_be_bytes(tmp) + } +} + impl From for u32 { fn from(other: SpinState) -> Self { const NBYTES: usize = 4; @@ -135,7 +147,7 @@ impl From for u32 { for i in 0..ARRAY_SIZE { tmp[i] = other.state[i]; } - ::from_ne_bytes(tmp) + ::from_be_bytes(tmp) } } @@ -147,7 +159,7 @@ impl From for u64 { for i in 0..ARRAY_SIZE { tmp[i] = other.state[i]; } - ::from_ne_bytes(tmp) + ::from_be_bytes(tmp) } } @@ -159,7 +171,7 @@ impl From for u128 { for i in 0..ARRAY_SIZE { tmp[i] = other.state[i]; } - ::from_ne_bytes(tmp) + ::from_be_bytes(tmp) } } @@ -375,6 +387,15 @@ impl std::ops::Shl for SpinState { } } +impl fmt::Display for SpinState { + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + for i in 0..ARRAY_SIZE { + write!(f, "{:0b} ", self.state[i])?; + } + write!(f, "") + } +} + /// The Fock state structure. Encodes the spins positions. /// # Definition /// This structure has two different fields, the spin up and spin down component. @@ -392,14 +413,14 @@ impl std::ops::Shl for SpinState { /// assert_eq!(state_5_both << 2, state_20_both); /// ``` #[derive(Debug, Eq, PartialEq, Clone, Copy)] -pub struct FockState where T: From +pub struct FockState where T: From + std::fmt::Display { pub spin_up: T, pub spin_down: T, pub n_sites: usize, } -impl + From> std::ops::Shl for FockState { +impl + From + std::fmt::Display> std::ops::Shl for FockState { type Output = Self; fn shl(self, u: usize) -> Self::Output { @@ -411,7 +432,7 @@ impl + From> std::ops::Shl } } -impl + From> std::ops::Shr for FockState { +impl + From + std::fmt::Display> std::ops::Shr for FockState { type Output = Self; fn shr(self, u: usize) -> Self::Output { @@ -423,7 +444,7 @@ impl + From> std::ops::Shr } } -impl fmt::Display for FockState { +impl fmt::Display for FockState { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { write!(f, "|")?; for i in 0..self.n_sites { @@ -449,7 +470,7 @@ pub trait Hopper { fn generate_all_exchange(self: &Self) -> Vec<(usize, usize)>; } -impl> Hopper for FockState { +impl + std::fmt::Display> Hopper for FockState { fn generate_all_hoppings(self: &FockState, bitmask: &[SpinState]) -> Vec<(usize, usize, Spin)> { let sup = self.spin_up; let sdo = self.spin_down; @@ -460,6 +481,7 @@ impl> Hopper for FockState { for j in 1..self.n_sites /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(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; @@ -562,7 +584,7 @@ pub trait RandomStateGeneration { } impl RandomStateGeneration for FockState -where T: BitOps, +where T: BitOps + std::fmt::Display, Standard: Distribution { fn generate_from_nelec(rng: &mut R, nelec: usize, max_size: usize) -> FockState { diff --git a/src/gutzwiller.rs b/src/gutzwiller.rs index 85693fb..1ebd3b4 100644 --- a/src/gutzwiller.rs +++ b/src/gutzwiller.rs @@ -34,7 +34,7 @@ pub fn compute_gutzwiller_exp( n_sites: usize, ) -> f64 where - T: BitOps + From + std::ops::Shl + Debug, + T: BitOps + From + std::ops::Shl + Debug + std::fmt::Display, { // sum_i g_i n_i up n_i down let mut gutzwiller_sites = fock_state.spin_up & fock_state.spin_down; @@ -56,7 +56,7 @@ pub fn compute_gutzwiller_der( der: &mut DerivativeOperator, ) where - T: BitOps + From + std::ops::Shl + Debug, + T: BitOps + From + std::ops::Shl + Debug + std::fmt::Display, { // sum_i g_i n_i up n_i down let mut gutzwiller_sites = fock_state.spin_up & fock_state.spin_down; diff --git a/src/hamiltonian.rs b/src/hamiltonian.rs index dff737b..5f7b49e 100644 --- a/src/hamiltonian.rs +++ b/src/hamiltonian.rs @@ -19,7 +19,7 @@ use crate::{VarParams, Spin, SysParams}; /// $$ pub fn potential(state: FockState, proj: f64, pstate: &PfaffianState, sys: &SysParams) -> f64 where - T: BitOps, + T: BitOps + std::fmt::Display, { let pot = ((state.spin_up & state.spin_down).count_ones() as f64) * sys.cons_u; trace!("Output potential = {:.2} for state |x> = {}", pot, state); diff --git a/src/hoppings.rs b/src/hoppings.rs index 656b27a..eee2c0d 100644 --- a/src/hoppings.rs +++ b/src/hoppings.rs @@ -1,3 +1,5 @@ +// Problem when compiling with N = ARRAY_SIZE, one array is out of bound, +// Should look into this, seems like a waste to use u32 for 16 Sites for example pub fn generate_bitmask(transfer_matrix: &[f64], size: usize) -> Vec { const WORD_SIZE: usize = 8; let all_zeros = SpinState{state: [0x00; ARRAY_SIZE]}; diff --git a/src/jastrow.rs b/src/jastrow.rs index 87d966e..eda0c0b 100644 --- a/src/jastrow.rs +++ b/src/jastrow.rs @@ -113,7 +113,7 @@ fn jastrow_undo_update( index_skip: usize, n_sites: usize, ) where - T: BitOps, + T: BitOps + std::fmt::Display, { if spin_mask.check(index_j) { //*spin_mask = *spin_mask & (::ones() >> (index_j + 1)); @@ -152,7 +152,7 @@ fn jastrow_do_update( index_skip: usize, n_sites: usize, ) where - T: BitOps, + T: BitOps + std::fmt::Display, { if spin_mask.check(index_j) { //*spin_mask = *spin_mask & (::ones() >> (index_j + 1)); @@ -192,7 +192,7 @@ fn jastrow_single_update( index_i: usize, sign: bool, ) where - T: BitOps, + T: BitOps + std::fmt::Display, { if spin_mask.check(index_j) & spin_mask.check(index_i) { let (n1, n2) = (fock_state.spin_up, fock_state.spin_down); diff --git a/src/optimisation.rs b/src/optimisation.rs index 0084d7d..00b6d49 100644 --- a/src/optimisation.rs +++ b/src/optimisation.rs @@ -267,7 +267,6 @@ impl<'a> ReducibleGeneralRepresentation for GenParameterMap { } let n = self.gendim as usize; let mut out_x0 = vec![0.0; n]; - println!("da = {:?}", x0); // Copy all other data // iter over transpose's diagonal From 1d870617a3f9167c212db7d7032892e0119c774e Mon Sep 17 00:00:00 2001 From: Dimitri Bonanni-Surprenant Date: Wed, 7 May 2025 20:17:48 -0400 Subject: [PATCH 6/7] (feat) Parallel initial impl. Works for 1, 4, 8 threads only, somehow --- Cargo.lock | 7 ++ Cargo.toml | 3 + src/bin/dvmc_pairwf2.rs | 84 ++++++++++++++---------- src/density.rs | 8 +-- src/dvmc.rs | 141 +++++++++++++++++++++++++++++++++++++--- src/fock_state.rs | 20 +++--- src/gutzwiller.rs | 4 +- src/hamiltonian.rs | 4 +- src/jastrow.rs | 12 ++-- src/monte_carlo.rs | 13 ++-- src/pfaffian.rs | 2 +- 11 files changed, 223 insertions(+), 75 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index a7bad1b..3afd84f 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -252,6 +252,12 @@ version = "0.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "2da6da31387c7e4ef160ffab6d5e7f00c42626fe39aea70a7b0f1773f7dd6c1b" +[[package]] +name = "closure" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d6173fd61b610d15a7566dd7b7620775627441c4ab9dac8906e17cb93a24b782" + [[package]] name = "colorchoice" version = "1.0.1" @@ -875,6 +881,7 @@ dependencies = [ "assert", "blas", "blas-sys", + "closure", "colored", "criterion", "csv", diff --git a/Cargo.toml b/Cargo.toml index bfb75ff..ade86d5 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -28,6 +28,7 @@ intel-mkl-src = {version = "0.8.1", features = ["mkl-dynamic-lp64-iomp"]} #openblas-src = "0.10.8" # For python interface pyo3 = {version = "0.21.2", features = ["extension-module"], optional = true} +# Parallelism rayon = "1.10.0" # Mersenne Twister rand_mt = {version = "4.2.2", features = ["rand_core"]} @@ -35,6 +36,8 @@ rand_mt = {version = "4.2.2", features = ["rand_core"]} indicatif = "0.17.8" # Colored terminal text colored = "3.0.0" +# Selectively move into closure +closure = "0.3.0" [dependencies.rand] version = "0.8.5" diff --git a/src/bin/dvmc_pairwf2.rs b/src/bin/dvmc_pairwf2.rs index 86b8941..f35114a 100644 --- a/src/bin/dvmc_pairwf2.rs +++ b/src/bin/dvmc_pairwf2.rs @@ -1,8 +1,10 @@ use impurity::optimisation::GenParameterMap; //use blas::dcopy; //use log::{debug, info}; -use rand_mt::Mt64; +use rand_mt::{Mt19937GenRand64, Mt64}; use rand::Rng; + +use rayon::ThreadPoolBuilder; //use indicatif::{ProgressBar, ProgressStyle}; use std::fs::File; use std::io::Write; @@ -13,18 +15,18 @@ use impurity::dvmc::{variationnal_monte_carlo, EnergyOptimisationMethod, EnergyC type BitSize = u32; const SEED: u64 = 14; -const SIZE: usize = 16; +const SIZE: usize = 2; 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 = 250; const NBOOTSTRAP: usize = 1; -const NMCWARMUP: usize = 1000; +const NMCWARMUP: usize = 100; const NWARMUPCHAINS: usize = 1; const MCSAMPLE_INTERVAL: usize = 1; -const _NTHREADS: usize = 1; +const NTHREADS: usize = 8; const CLEAN_UPDATE_FREQUENCY: usize = 32; const TOLERENCE_SHERMAN_MORRISSON: f64 = 1e-12; const TOLERENCE_SINGULARITY: f64 = 1e-12; @@ -56,13 +58,13 @@ const PAIRWF: bool = false; const CONV_PARAM_THRESHOLD: f64 = 1e-10; // 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 -//]; +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] = [ @@ -115,24 +117,24 @@ const CONV_PARAM_THRESHOLD: f64 = 1e-10; //]; // 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, -]; +//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, +//]; // // Pairwf @@ -228,6 +230,7 @@ fn log_energy_convs(e: &[f64], fp: &mut File) { fn main() { let mut fp = File::create("u_t_sweep").unwrap(); + let pool = ThreadPoolBuilder::new().num_threads(NTHREADS).build_global().unwrap(); writeln!(fp, "{}", format!("# {} {} {}", SIZE, NMCSAMP, NOPTITER)).unwrap(); let mut paramsfp = File::create("params").unwrap(); writeln!(paramsfp, "{}", format!("# {} {} {}", SIZE, NMCSAMP, NOPTITER)).unwrap(); @@ -235,7 +238,13 @@ fn main() { // Initialize logger env_logger::init(); let bitmask = generate_bitmask(&HOPPINGS, SIZE); - let mut rng = Mt64::new(SEED); + let mut rng = Vec::new(); + let mut rngs = Vec::new(); + for i in 0..NTHREADS { + rng.push(Mt64::new(SEED + i as u64)); + let rng_ptr: *mut Mt19937GenRand64 = &mut rng[i] as *mut _; + rngs.push(unsafe{&mut *rng_ptr}); + } for nsweep_iter in 0..NRATIO_POINTS { let mut system_params = SysParams { @@ -264,7 +273,7 @@ fn main() { let mut all_params: Vec = Vec::with_capacity(NGI + NVIJ + NFIJ); for _ in 0..(NGI + NVIJ + NFIJ) { - all_params.push(rng.gen()); + all_params.push(rngs[0].gen()); } let (gi, params) = all_params.split_at_mut(NGI); let (vij, fij) = params.split_at_mut(NVIJ); @@ -298,15 +307,20 @@ fn main() { compute_energy_method: COMPUTE_ENERGY_METHOD, optimise_energy_method: OPTIMISE_ENERGY_METHOD, conv_param_threshold: CONV_PARAM_THRESHOLD, + nthreads: NTHREADS }; let state: FockState = { - let mut tmp: FockState = FockState::generate_from_nelec(&mut rng, NELEC, SIZE); + let mut tmp: FockState = FockState::generate_from_nelec(&mut rngs[0], NELEC, SIZE); while tmp.spin_up.count_ones() != tmp.spin_down.count_ones() { - tmp = FockState::generate_from_nelec(&mut rng, NELEC, SIZE); + tmp = FockState::generate_from_nelec(&mut rngs[0], NELEC, SIZE); } tmp }; + let mut states_vec = Vec::with_capacity(NTHREADS); + for i in 0..NTHREADS { + states_vec.push(state); + } let param_map = GenParameterMap { dim: N_INDEP_PARAMS as i32, gendim: (NFIJ + NGI + NVIJ) as i32, @@ -316,7 +330,7 @@ fn main() { projector: Box::new(PARAMS_PROJECTOR), }; - let e_array = variationnal_monte_carlo(&mut rng, state, &mut parameters, &mut system_params, &vmcparams, ¶m_map); + let e_array = variationnal_monte_carlo(&mut rngs, &mut states_vec, &mut parameters, &mut system_params, &vmcparams, ¶m_map); //write_energy(&mut fp, &e_array); log_energy_convs(&e_array, &mut paramsfp); diff --git a/src/density.rs b/src/density.rs index 3f5587e..6b3e911 100644 --- a/src/density.rs +++ b/src/density.rs @@ -12,7 +12,7 @@ pub fn compute_internal_product( params: &VarParams, sys: &SysParams, ) -> f64 -where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl +where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl + Send { let (mut pfaffian_state, jastrow_exp, gutz_exp) = { ( @@ -34,7 +34,7 @@ pub fn compute_internal_product_parts( params: &VarParams, sys: &SysParams, ) -> (PfaffianState, f64) -where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl +where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl + Send { let (mut pfaffian_state, jastrow_exp, gutz_exp) = { ( @@ -55,7 +55,7 @@ pub fn fast_internal_product_no_otilde( previous_proj: &mut f64, params: &VarParams, ) -> (f64, Vec, usize) -where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl +where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl + Send { // Rename things. let previous_i = hopping.0; @@ -94,7 +94,7 @@ pub fn fast_internal_product( previous_proj: &mut f64, params: &VarParams, ) -> (f64, Vec, usize) -where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl +where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl + Send { // Rename things. let previous_i = hopping.0; diff --git a/src/dvmc.rs b/src/dvmc.rs index 8c686f3..424c0e8 100644 --- a/src/dvmc.rs +++ b/src/dvmc.rs @@ -1,8 +1,11 @@ use log::{info, error}; use rand::distributions::{Distribution, Standard}; use rand::Rng; +use rayon::iter::PanicFuse; use std::fmt::{Display, Debug}; use blas::{idamax, daxpy, dcopy, dscal, dnrm2}; +use rayon::ThreadPool; +use std::thread; use crate::{VarParams, DerivativeOperator, SysParams, FockState, BitOps}; use crate::monte_carlo::{compute_mean_energy, compute_mean_energy_exact}; @@ -31,6 +34,7 @@ pub struct VMCParams { pub dt: f64, pub optimisation_decay: f64, pub nparams: usize, + pub nthreads: usize, pub threshold: f64, pub optimise: bool, pub optimise_gutzwiller: bool, @@ -39,6 +43,42 @@ pub struct VMCParams { pub conv_param_threshold: f64, } +fn merge_derivatives(der_vec: &[DerivativeOperator], nmc: usize, nthreads: usize) -> DerivativeOperator { + let mut out_der = DerivativeOperator::new( + der_vec[0].n, + 0, + nmc as f64, + der_vec[0].nsamp_int, + der_vec[0].pfaff_off, + der_vec[0].jas_off, + der_vec[0].epsilon + ); + let mut k = 0; + for i in 0..nthreads { + out_der.mu += der_vec[i].mu; + unsafe { + let incx = 1; + let incy = 1; + dcopy( + der_vec[i].n * der_vec[i].mu, + &der_vec[i].o_tilde, + incx, + &mut out_der.o_tilde[(der_vec[i].n * der_vec[i].mu) as usize * i..der_vec[i].n as usize * nmc], + incy + ); + daxpy(der_vec[i].n, 1.0 / nthreads as f64, &der_vec[i].expval_o, incx, &mut out_der.expval_o, incy); + daxpy(der_vec[i].n, 1.0 / nthreads as f64, &der_vec[i].ho, incx, &mut out_der.ho, incy); + } + for j in 0..der_vec[i].mu as usize { + out_der.visited[k] = der_vec[i].visited[j]; + k += 1; + + } + } + + out_der +} + fn zero_out_derivatives(der: &mut DerivativeOperator, sys: &SysParams) { for i in 0.. (der.n as usize) * sys.nmcsample { der.o_tilde[i] = 0.0; @@ -53,15 +93,86 @@ fn zero_out_derivatives(der: &mut DerivativeOperator, sys: &SysParams) { der.mu = -1; } -pub fn variationnal_monte_carlo( - rng: &mut R, - initial_state: FockState, +fn update_initial_state + Display + Debug + Send + Sync>( + states: &mut [FockState], + acc_states: &[FockState], + nmcsample: usize, + nthreads: usize +) { + for i in 0..nthreads { + states[i] = acc_states[i * nmcsample]; + } +} + +fn parallel_monte_carlo( + rngs: &mut [&mut R], + initial_state: &[FockState], + work_der_vec: &mut [DerivativeOperator], + vmcparams: &VMCParams, + params: & VarParams, + sys: & SysParams +) -> (f64, Vec>, f64, f64) +where + T: BitOps + From + Display + Debug + Send + Sync, + R: Rng + ?Sized + Send + Sync, + Standard: Distribution + Send +{ + let mut res_vec = Vec::new(); + // Scoped thread because threads need to accept not static stack parameters + thread::scope(|scope| { + let threads: Vec<_> = (0..vmcparams.nthreads) + .map(|idx| { + // Wrestle the borrow checker + let wder_ptr = &mut work_der_vec[idx] as *mut _; + let rng_ptr = &mut *rngs[idx] as *mut R; + let state_ptr = &initial_state[idx] as *const _; + let rng = unsafe { &mut *rng_ptr}; + let wder = unsafe { &mut *wder_ptr}; + let state = unsafe { & *state_ptr}; + scope.spawn( + || { + compute_mean_energy( + rng, + *state, + params, + sys, + wder + ) + } + ) + }) + .collect(); + + for handle in threads { + res_vec.push(handle.join().unwrap()); + } + }); + let mut out_energy = 0.0; + let mut out_states_vec = Vec::new(); + let mut out_de = 0.0; + let mut out_corrtime = 0.0; + for thread in 0..vmcparams.nthreads { + //println!("Energy of thread {} = {}", thread, res_vec[thread].0); + out_energy += res_vec[thread].0; + out_de += res_vec[thread].2; + out_corrtime += res_vec[thread].3; + out_states_vec.append(&mut res_vec[thread].1); + } + + //panic!("Stop!"); + //(res_vec[0].0, res_vec[0].1.clone(), res_vec[0].2, res_vec[0].3) + (out_energy / vmcparams.nthreads as f64, out_states_vec, out_de / vmcparams.nthreads as f64, out_corrtime / vmcparams.nthreads as f64) +} + +pub fn variationnal_monte_carlo( + rng: &mut [&mut R], + initial_state: &mut [FockState], params: &mut VarParams, sys: &mut SysParams, vmcparams: &VMCParams, param_map: &GenParameterMap, ) -> Vec -where T: BitOps + From + Display + Debug, Standard: Distribution +where T: BitOps + From + Display + Debug + Send + Sync, Standard: Distribution + std::marker::Send { let mut output_energy_array = vec![0.0; vmcparams.noptiter * 3]; @@ -72,14 +183,18 @@ where T: BitOps + From + Display + Debug, Standard: Distribution -1, match vmcparams.compute_energy_method { EnergyComputationMethod::ExactSum => 1.0, - EnergyComputationMethod::MonteCarlo => sys.nmcsample as f64, + EnergyComputationMethod::MonteCarlo => (sys.nmcsample * vmcparams.nthreads) as f64, }, sys.mcsample_interval, param_map.n_independant_jastrow + param_map.n_independant_gutzwiller, param_map.n_independant_gutzwiller, vmcparams.epsilon ); - let mut work_der = param_map.mapto_general_representation(&der, &mut x0); + let mut work_der_vec = Vec::new(); + for i in 0..vmcparams.nthreads { + let mut work_der = param_map.mapto_general_representation(&der, &mut x0); + work_der_vec.push(work_der); + } //let mut der = DerivativeOperator { // o_tilde: &mut otilde, // expval_o: &mut expvalo, @@ -116,12 +231,15 @@ where T: BitOps + From + Display + Debug, Standard: Distribution sys._opt_iter = opt_iter; let (mean_energy, _accumulated_states, deltae, correlation_time) = { match vmcparams.compute_energy_method { - EnergyComputationMethod::MonteCarlo => compute_mean_energy(rng, initial_state, params, sys, &mut work_der), + EnergyComputationMethod::MonteCarlo => { + parallel_monte_carlo(rng, initial_state, &mut work_der_vec, vmcparams, params, sys) + }, EnergyComputationMethod::ExactSum => { - (compute_mean_energy_exact(params, sys, &mut work_der), Vec::with_capacity(0), 0.0, 0.0) + (compute_mean_energy_exact(params, sys, &mut work_der_vec[0]), Vec::with_capacity(0), 0.0, 0.0) }, } }; + 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; @@ -139,7 +257,8 @@ where T: BitOps + From + Display + Debug, Standard: Distribution x0[(sys.ngi + sys.nvij)..(sys.ngi + sys.nvij + sys.nfij)].copy_from_slice(params.fij); x0[sys.ngi..(sys.ngi + sys.nvij)].copy_from_slice(params.vij); x0[0..sys.ngi].copy_from_slice(params.gi); - param_map.update_reduced_representation(&work_der, &mut der, &mut x0); + //let work_der = merge_derivatives(&work_der_vec, sys.nmcsample, vmcparams.nthreads); + param_map.update_reduced_representation(&work_der_vec[0], &mut der, &mut x0); // 68 misawa //let mut b: Vec = vec![0.0; der.n as usize]; @@ -249,7 +368,9 @@ where T: BitOps + From + Display + Debug, Standard: Distribution // ); //} zero_out_derivatives(&mut der, sys); - zero_out_derivatives(&mut work_der, sys); + for i in 0..vmcparams.nthreads { + zero_out_derivatives(&mut work_der_vec[i], sys); + } //print_delta_alpha(&delta_alpha, sys.ngi, sys.nvij, sys.nfij); //let opt_delta = unsafe { // let incx = 1; diff --git a/src/fock_state.rs b/src/fock_state.rs index ad58e4e..c2df100 100644 --- a/src/fock_state.rs +++ b/src/fock_state.rs @@ -5,6 +5,7 @@ use rand::Rng; use rand::distributions::{Distribution, Standard}; use std::fmt; use log::error; +use std::marker::Send; use strings::{UPARROW, DOWNARROW}; @@ -40,7 +41,8 @@ pub trait BitOps: std::ops::Shr + std::ops::Shl + std::ops::BitOr + - From + From + + Sync { /// Provides the number of leading zeros in the bitstring. This gives the /// position of the first set bit in the string. This method is consistent @@ -72,7 +74,7 @@ pub trait BitOps: /// built-in methods. [BitOps::set] and [BitOps::check] are implemented by /// shifting a bitmask. impl BitOps for I - where I: PrimInt + std::ops::BitXorAssign + std::ops::BitAndAssign + From + From + where I: PrimInt + std::ops::BitXorAssign + std::ops::BitAndAssign + From + From + Sync { #[inline(always)] fn leading_zeros(self) -> u32 { @@ -413,14 +415,14 @@ impl fmt::Display for SpinState { /// assert_eq!(state_5_both << 2, state_20_both); /// ``` #[derive(Debug, Eq, PartialEq, Clone, Copy)] -pub struct FockState where T: From + std::fmt::Display +pub struct FockState where T: From + std::fmt::Display + Send { pub spin_up: T, pub spin_down: T, pub n_sites: usize, } -impl + From + std::fmt::Display> std::ops::Shl for FockState { +impl + From + std::fmt::Display + Send> std::ops::Shl for FockState { type Output = Self; fn shl(self, u: usize) -> Self::Output { @@ -432,7 +434,7 @@ impl + From + std::fmt::Display> } } -impl + From + std::fmt::Display> std::ops::Shr for FockState { +impl + From + std::fmt::Display + Send> std::ops::Shr for FockState { type Output = Self; fn shr(self, u: usize) -> Self::Output { @@ -444,7 +446,7 @@ impl + From + std::fmt::Display> } } -impl fmt::Display for FockState { +impl fmt::Display for FockState { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { write!(f, "|")?; for i in 0..self.n_sites { @@ -470,7 +472,7 @@ pub trait Hopper { fn generate_all_exchange(self: &Self) -> Vec<(usize, usize)>; } -impl + std::fmt::Display> Hopper for FockState { +impl + std::fmt::Display + Send> Hopper for FockState { fn generate_all_hoppings(self: &FockState, bitmask: &[SpinState]) -> Vec<(usize, usize, Spin)> { let sup = self.spin_up; let sdo = self.spin_down; @@ -565,7 +567,7 @@ impl + std::fmt::Display> Hopper for FockState { } // Interface for random state generation -impl Distribution> for Standard where Standard: Distribution { +impl Distribution> for Standard where Standard: Distribution { fn sample(&self, rng: &mut R) -> FockState { let su = rng.gen::(); let sd = rng.gen::(); @@ -584,7 +586,7 @@ pub trait RandomStateGeneration { } impl RandomStateGeneration for FockState -where T: BitOps + std::fmt::Display, +where T: BitOps + std::fmt::Display + Send, Standard: Distribution { fn generate_from_nelec(rng: &mut R, nelec: usize, max_size: usize) -> FockState { diff --git a/src/gutzwiller.rs b/src/gutzwiller.rs index 1ebd3b4..a79664b 100644 --- a/src/gutzwiller.rs +++ b/src/gutzwiller.rs @@ -34,7 +34,7 @@ pub fn compute_gutzwiller_exp( n_sites: usize, ) -> f64 where - T: BitOps + From + std::ops::Shl + Debug + std::fmt::Display, + T: BitOps + From + std::ops::Shl + Debug + std::fmt::Display + Send, { // sum_i g_i n_i up n_i down let mut gutzwiller_sites = fock_state.spin_up & fock_state.spin_down; @@ -56,7 +56,7 @@ pub fn compute_gutzwiller_der( der: &mut DerivativeOperator, ) where - T: BitOps + From + std::ops::Shl + Debug + std::fmt::Display, + T: BitOps + From + std::ops::Shl + Debug + std::fmt::Display + Send, { // sum_i g_i n_i up n_i down let mut gutzwiller_sites = fock_state.spin_up & fock_state.spin_down; diff --git a/src/hamiltonian.rs b/src/hamiltonian.rs index 5f7b49e..38ca019 100644 --- a/src/hamiltonian.rs +++ b/src/hamiltonian.rs @@ -19,7 +19,7 @@ use crate::{VarParams, Spin, SysParams}; /// $$ pub fn potential(state: FockState, proj: f64, pstate: &PfaffianState, sys: &SysParams) -> f64 where - T: BitOps + std::fmt::Display, + T: BitOps + std::fmt::Display + Send, { let pot = ((state.spin_up & state.spin_down).count_ones() as f64) * sys.cons_u; trace!("Output potential = {:.2} for state |x> = {}", pot, state); @@ -40,7 +40,7 @@ where /// $$ 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 + T: BitOps + From + std::fmt::Debug + std::fmt::Display + Send { let hops = state.generate_all_hoppings(&sys.hopping_bitmask); diff --git a/src/jastrow.rs b/src/jastrow.rs index eda0c0b..ba90bcf 100644 --- a/src/jastrow.rs +++ b/src/jastrow.rs @@ -34,7 +34,7 @@ pub fn compute_jastrow_exp( n_sites: usize, ) -> f64 where - T: BitOps + std::fmt::Display, + T: BitOps + std::fmt::Display + Send, { let mut jastrow_out = 0.0; let mut regular_nor = !(fock_state.spin_up ^ fock_state.spin_down); @@ -73,7 +73,7 @@ pub fn compute_jastrow_der( n_sites: usize, ) where - T: BitOps + std::fmt::Display, + T: BitOps + std::fmt::Display + Send, { let mut regular_nor = !(fock_state.spin_up ^ fock_state.spin_down); regular_nor.mask_bits(n_sites); @@ -113,7 +113,7 @@ fn jastrow_undo_update( index_skip: usize, n_sites: usize, ) where - T: BitOps + std::fmt::Display, + T: BitOps + std::fmt::Display + Send, { if spin_mask.check(index_j) { //*spin_mask = *spin_mask & (::ones() >> (index_j + 1)); @@ -152,7 +152,7 @@ fn jastrow_do_update( index_skip: usize, n_sites: usize, ) where - T: BitOps + std::fmt::Display, + T: BitOps + std::fmt::Display + Send, { if spin_mask.check(index_j) { //*spin_mask = *spin_mask & (::ones() >> (index_j + 1)); @@ -192,7 +192,7 @@ fn jastrow_single_update( index_i: usize, sign: bool, ) where - T: BitOps + std::fmt::Display, + T: BitOps + std::fmt::Display + Send, { if spin_mask.check(index_j) & spin_mask.check(index_i) { let (n1, n2) = (fock_state.spin_up, fock_state.spin_down); @@ -272,7 +272,7 @@ pub fn fast_update_jastrow( previous_j: usize, new_j: usize, ) where - T: BitOps + std::fmt::Display, + T: BitOps + std::fmt::Display + Send, { // Undo previous fock state // 1/2 sum {i\neq j} v_{ij}(n_i - 1)(n_j - 1) diff --git a/src/monte_carlo.rs b/src/monte_carlo.rs index 5d02fa7..198724c 100644 --- a/src/monte_carlo.rs +++ b/src/monte_carlo.rs @@ -2,6 +2,7 @@ use blas::daxpy; use log::{error, info, trace, warn}; use rand::distributions::{Distribution, Standard}; use rand::Rng; +use std::sync::Arc; use crate::gutzwiller::compute_gutzwiller_der; use crate::jastrow::compute_jastrow_der; @@ -34,7 +35,7 @@ use crate::hamiltonian::{kinetic, potential}; #[inline(always)] fn propose_hopping > +T: BitOps + std::fmt::Display + std::fmt::Debug + From + Send> ( state: &FockState, pfaff_state: &PfaffianState, @@ -54,7 +55,7 @@ T: BitOps + std::fmt::Display + std::fmt::Debug + From> } #[inline(always)] -fn compute_hamiltonian(state: FockState, pstate: &PfaffianState, proj: f64, params: &VarParams, sys: &SysParams) -> f64 { +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); @@ -62,7 +63,7 @@ fn compute_hamiltonian(state: F } #[inline(always)] -fn compute_derivative_operator> +fn compute_derivative_operator + Send> (state: FockState, pstate: &PfaffianState, der: &mut DerivativeOperator, sys: &SysParams) { compute_gutzwiller_der(state, sys.size, der); @@ -71,7 +72,7 @@ fn compute_derivative_operator>( +fn make_update + Send>( n: &mut usize, proj: &mut f64, proj_copy: &mut f64, @@ -126,7 +127,7 @@ fn warmup( params: &VarParams, sys: &SysParams ) -where T: BitOps + From + std::fmt::Debug + std::fmt::Display, +where T: BitOps + From + std::fmt::Debug + std::fmt::Display + Send, R: Rng + ?Sized, Standard: Distribution { @@ -284,7 +285,7 @@ pub fn compute_mean_energy_exact(params: &VarParams, sys: &SysParams, der: &mut pub fn compute_mean_energy > +T: BitOps + std::fmt::Debug + std::fmt::Display + From + Send> (rng: &mut R, initial_state: FockState, params: &VarParams, sys: &SysParams, derivatives: &mut DerivativeOperator) -> (f64, Vec>, f64, f64) where Standard: Distribution { diff --git a/src/pfaffian.rs b/src/pfaffian.rs index a4be29a..684332b 100644 --- a/src/pfaffian.rs +++ b/src/pfaffian.rs @@ -127,7 +127,7 @@ fn transpose(a: &Vec, n: usize) -> Vec{ /// * __`state`__ - The state of the system. pub fn construct_matrix_a_from_state(fij: &[f64], state: FockState, sys: &SysParams) -> PfaffianState where - T: BitOps + std::fmt::Display, + T: BitOps + std::fmt::Display + Send, { // Fij upup, updown, downup, downdown let n = state.spin_up.count_ones() as usize + state.spin_down.count_ones() as usize; From e0656dc26f77341b18b4663a0f371d8d55954ec9 Mon Sep 17 00:00:00 2001 From: Dimitri Bonanni-Surprenant Date: Fri, 9 May 2025 14:14:01 -0400 Subject: [PATCH 7/7] (fix) Use after free fixed --- src/bin/dvmc_pairwf2.rs | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/bin/dvmc_pairwf2.rs b/src/bin/dvmc_pairwf2.rs index f35114a..955be9f 100644 --- a/src/bin/dvmc_pairwf2.rs +++ b/src/bin/dvmc_pairwf2.rs @@ -3,6 +3,7 @@ use impurity::optimisation::GenParameterMap; //use log::{debug, info}; use rand_mt::{Mt19937GenRand64, Mt64}; use rand::Rng; +use std::mem; use rayon::ThreadPoolBuilder; //use indicatif::{ProgressBar, ProgressStyle}; @@ -21,12 +22,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 = 250; +const NMCSAMP: usize = 200; const NBOOTSTRAP: usize = 1; const NMCWARMUP: usize = 100; const NWARMUPCHAINS: usize = 1; const MCSAMPLE_INTERVAL: usize = 1; -const NTHREADS: usize = 8; +const NTHREADS: usize = 12; const CLEAN_UPDATE_FREQUENCY: usize = 32; const TOLERENCE_SHERMAN_MORRISSON: f64 = 1e-12; const TOLERENCE_SINGULARITY: f64 = 1e-12; @@ -239,12 +240,10 @@ fn main() { env_logger::init(); let bitmask = generate_bitmask(&HOPPINGS, SIZE); let mut rng = Vec::new(); - let mut rngs = Vec::new(); for i in 0..NTHREADS { rng.push(Mt64::new(SEED + i as u64)); - let rng_ptr: *mut Mt19937GenRand64 = &mut rng[i] as *mut _; - rngs.push(unsafe{&mut *rng_ptr}); } + let mut rngs: Vec<&mut Mt64> = rng.iter_mut().collect(); for nsweep_iter in 0..NRATIO_POINTS { let mut system_params = SysParams { @@ -336,4 +335,5 @@ fn main() { log_energy_convs(&e_array, &mut paramsfp); } + mem::drop(rng); }