diff --git a/cSPCE_DRISM.toml b/cSPCE_DRISM.toml index acccbd8f..013e5e77 100644 --- a/cSPCE_DRISM.toml +++ b/cSPCE_DRISM.toml @@ -3,7 +3,7 @@ temp = 298 kT = 1.0 kU = 0.0019872158728366637 charge_coeff = 167101.0 -npts = 1 +npts = 16384 radius = 409.6 lam = 1 @@ -13,12 +13,12 @@ closure = "KH" IE = "DRISM" solver = "MDIIS" depth = 12 -picard_damping = 0.1 -mdiis_damping = 0.6 -itermax = 300 +picard_damping = 0.5 +mdiis_damping = 0.5 +itermax = 10000 tol = 1E-8 diel = 78.497 -adbcor = 0.1 +adbcor = 1.5 [solvent] nsv = 3 diff --git a/pyrism/Solvers/MDIIS.py b/pyrism/Solvers/MDIIS.py index 72746b03..7b430332 100644 --- a/pyrism/Solvers/MDIIS.py +++ b/pyrism/Solvers/MDIIS.py @@ -29,7 +29,6 @@ class MDIIS(SolverObject): def step_Picard(self, curr, prev): self.fr.append(curr.flatten()) self.res.append((curr - prev).flatten()) - print(self.damp_picard) return prev + self.damp_picard * (curr - prev) def step_MDIIS(self, curr, prev, gr): @@ -73,7 +72,6 @@ def solve(self, RISM, Closure, lam, verbose=False): * np.power((c_A - c_prev).sum(), 2) ) self.RMS_res.append(RMS) - print(self.RMS_res) else: c_next = self.step_MDIIS(c_A, c_prev, self.data_vv.t + c_A) c_next = np.reshape(c_next, c_prev.shape) @@ -83,7 +81,6 @@ def solve(self, RISM, Closure, lam, verbose=False): / self.data_vv.grid.npts * np.power((c_A - c_prev).sum(), 2) ) - print(min(self.RMS_res)) if RMS > 10 * min(self.RMS_res): print("Restarting MDIIS") min_index = self.RMS_res.index(min(self.RMS_res)) diff --git a/src/closure.rs b/src/closure.rs index aeccd081..3ad11132 100644 --- a/src/closure.rs +++ b/src/closure.rs @@ -69,10 +69,10 @@ pub fn hyper_netted_chain(problem: &DataRs) -> Array3 { pub fn kovalenko_hirata(problem: &DataRs) -> Array3 { let mut out = Array::zeros(problem.correlations.tr.raw_dim()); par_azip!((a in &mut out, &b in &(-problem.system.beta * &problem.interactions.u_sr + &problem.correlations.tr), &c in &problem.correlations.tr) { - if b < 0.0 { - *a = b.exp() + if b <= 0.0 { + *a = b.exp() - 1.0 - c } else { - *a = b + c + *a = b - c } }); out diff --git a/src/integralequation.rs b/src/integralequation.rs index c11baa64..c23213cd 100644 --- a/src/integralequation.rs +++ b/src/integralequation.rs @@ -146,7 +146,6 @@ fn rism_vv_equation_impl( let wcw = w_bar.dot(&ck_matrix.dot(&w_bar)); hk_matrix.assign(&(inverted_iwcp.dot(&wcw) + chi_matrix)); }); - println!("{:#?}", hk); // Compute t(k) = h(k) - c(k) let tk = &hk - ck; diff --git a/src/mdiis.rs b/src/mdiis.rs index 7f9c5d66..2b1db7b5 100644 --- a/src/mdiis.rs +++ b/src/mdiis.rs @@ -1,7 +1,9 @@ use crate::data::DataRs; use crate::operator::Operator; use crate::solver::{Solver, SolverError, SolverSettings}; -use log::{debug, warn}; +use fftw::plan::*; +use fftw::types::*; +use log::warn; use ndarray_linalg::Solve; use numpy::ndarray::{Array, Array1, Array2, Array3}; use std::collections::VecDeque; @@ -39,12 +41,7 @@ impl MDIIS { } } - fn step_picard( - &mut self, - curr: &Array3, - prev: &Array3, - res: &Array3, - ) -> Array3 { + fn step_picard(&mut self, curr: &Array3, prev: &Array3) -> Array3 { // calculate difference between current and previous solutions from RISM equation let diff = curr.clone() - prev.clone(); @@ -54,10 +51,9 @@ impl MDIIS { // push flattened difference into residual array self.res - .push_back(Array::from_iter(res.clone().into_iter())); + .push_back(Array::from_iter(diff.clone().into_iter())); // return Picard iteration step - println!("{}", self.picard_damping); prev + self.picard_damping * diff } @@ -65,7 +61,6 @@ impl MDIIS { &mut self, curr: &Array3, prev: &Array3, - res: &Array3, gr: &Array3, ) -> Array1 { let mut a = Array2::zeros((self.m + 1, self.m + 1)); @@ -109,7 +104,7 @@ impl MDIIS { // push flattened difference into residual array self.res - .push_back(Array::from_iter(res.clone().into_iter())); + .push_back(Array::from_iter(diff.clone().into_iter())); self.fr.pop_front(); self.res.pop_front(); @@ -133,27 +128,24 @@ impl Solver for MDIIS { let c_prev = problem.correlations.cr.clone(); (operator.eq)(problem); let c_a = (operator.closure)(&problem); - let g_a = &c_a + 1.0 + &problem.correlations.tr; let mut c_next; if self.fr.len() < self.m { //println!("Picard Step"); - c_next = self.step_picard(&c_a, &c_prev, &(&g_a - 1.0 - &problem.correlations.hr)); + c_next = self.step_picard(&c_a, &c_prev); let rmse = compute_rmse(ns1, ns2, npts, &c_a, &c_prev); //println!("\tMDIIS RMSE: {}", rmse); - self.rms_res.push_back(rmse); - debug!("RMS Vector:\n{:#?}", self.rms_res); + self.rms_res.push_back(rmse) } else { + let gr = &c_a + &problem.correlations.tr; //println!("MDIIS Step"); - let gr = &problem.correlations.tr + &c_a + 1.0; c_next = self - .step_mdiis(&c_a, &c_prev, &(&g_a - 1.0 - &problem.correlations.hr), &gr) + .step_mdiis(&c_a, &c_prev, &gr) .into_shape(shape) .expect("could not reshape array into original shape"); let rmse = compute_rmse(ns1, ns2, npts, &c_a, &c_prev); //println!("\tMDIIS RMSE: {}", rmse); let rmse_min = self.rms_res.iter().fold(f64::INFINITY, |a, &b| a.min(b)); - println!("rmse-min: {}", rmse_min); let min_index = self .rms_res .iter()