Converging on DRISM+KH+MDIIS now

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

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))

src/closure.rs

@@ -69,10 +69,10 @@ pub fn hyper_netted_chain(problem: &DataRs) -> Array3<f64> {
 pub fn kovalenko_hirata(problem: &DataRs) -> Array3<f64> {
     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

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;

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<f64>,
-        prev: &Array3<f64>,
-        res: &Array3<f64>,
-    ) -> Array3<f64> {
+    fn step_picard(&mut self, curr: &Array3<f64>, prev: &Array3<f64>) -> Array3<f64> {
         // calculate difference between current and previous solutions from RISM equation
         let diff = curr.clone() - prev.clone();
 
@@ -54,18 +51,16 @@ 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
     }
 
     fn step_mdiis(
         &mut self,
         curr: &Array3<f64>,
         prev: &Array3<f64>,
-        res: &Array3<f64>,
         gr: &Array3<f64>,
     ) -> Array1<f64> {
         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()