debugging drism iterations

2AUK · Oct 10, 2023 · dd97c4a · dd97c4a
1 parent 1d97308
commit dd97c4a
Show file tree

Hide file tree

Showing 5 changed files with 14 additions and 19 deletions.
diff --git 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 = 16384
+npts = 1
 radius = 409.6 
 lam = 1
 
@@ -12,10 +12,10 @@ potential = "LJ"
 closure = "KH"
 IE = "DRISM"
 solver = "MDIIS"
-depth = 16
+depth = 2
 picard_damping = 0.5
 mdiis_damping = 0.6
-itermax = 10000
+itermax = 10
 tol = 1E-8
 diel = 78.497
 adbcor = 1.5

diff --git a/pyrism/IntegralEquations/DRISM.py b/pyrism/IntegralEquations/DRISM.py
@@ -41,6 +41,7 @@ def compute_vv(self):
             self.data_vv.p,
             self.chi,
         )
+        print(self.data_vv.h)
         self.data_vv.t = (
             self.data_vv.grid.idht(self.data_vv.h - ck)
             - self.data_vv.B * self.data_vv.ur_lr
@@ -121,10 +122,9 @@ def D_matrix(self):
         d1z = np.zeros((self.data_vv.ns1), dtype=np.float64)
         for ki, k in enumerate(self.data_vv.grid.ki):
             hck = self.h_c0 * np.exp(-np.power((self.adbcor * k / 2.0), 2.0))
-            i = -1
+            i = 0 
             for isp in self.data_vv.species:
                 for iat in isp.atom_sites:
-                    i += 1
                     k_coord = k * iat.coords
                     if k_coord[0] == 0.0:
                         d0x[i] = 1.0
@@ -138,6 +138,7 @@ def D_matrix(self):
                         d1z[i] = 0.0
                     else:
                         d1z[i] = Util.j1(k_coord[2])
+                    i+=1
             for i, j in np.ndindex((self.data_vv.ns1, self.data_vv.ns2)):
                 self.chi[ki, i, j] = (
                     d0x[i] * d0y[i] * d1z[i] * d0x[j] * d0y[j] * d1z[j] * hck

diff --git a/pyrism/Solvers/MDIIS.py b/pyrism/Solvers/MDIIS.py
@@ -64,22 +64,15 @@ def solve(self, RISM, Closure, lam, verbose=False):
                 print("diff: {diff}".format(diff=(c_A - c_prev).min()))
                 raise e
             if len(self.fr) < self.m:
-                print("Picard Step")
                 c_next = self.step_Picard(c_A, c_prev)
                 RMS = np.sqrt(
                     1
                     / self.data_vv.ns1
                     / self.data_vv.grid.npts
                     * np.power((c_A - c_prev).sum(), 2)
                 )
-                print(
-                    "Iteration: {i}\nRMS: {RMS}\nDiff: {diff}".format(
-                        i=i, RMS=RMS, diff=(c_A - c_prev).min()
-                    )
-                )
                 self.RMS_res.append(RMS)
             else:
-                print("MDIIS Step")
                 c_next = self.step_MDIIS(c_A, c_prev, self.data_vv.t + c_A)
                 c_next = np.reshape(c_next, c_prev.shape)
                 RMS = np.sqrt(
@@ -88,11 +81,6 @@ def solve(self, RISM, Closure, lam, verbose=False):
                     / self.data_vv.grid.npts
                     * np.power((c_A - c_prev).sum(), 2)
                 )
-                print(
-                    "Iteration: {i}\nRMS: {RMS}\nDiff: {diff}".format(
-                        i=i, RMS=RMS, diff=(c_A - c_prev).min()
-                    )
-                )
                 if RMS > 10 * min(self.RMS_res):
                     print("Restarting MDIIS")
                     min_index = self.RMS_res.index(min(self.RMS_res))
@@ -102,6 +90,12 @@ def solve(self, RISM, Closure, lam, verbose=False):
                     self.RMS_res.clear()
                 self.RMS_res.append(RMS)
                 self.RMS_res.pop(0)
+            self.converged(c_next, c_prev) 
+            print(
+                    "Iteration: {i} Convergence RMS: {RMS}\nDiff: {diff}".format(
+                        i=i, RMS=self.rms, diff=(c_next - c_prev).min()
+                    )
+                )
 
             # print(c_next)
             self.data_vv.c = c_next

diff --git a/src/driver.rs b/src/driver.rs
@@ -236,7 +236,6 @@ impl RISMDriver {
                                 i += 1;
                             }
                         }
-
                         for i in 0..self.data.nsv {
                             for j in 0..self.data.nsv {
                                 chi[[ki, i, j]] =

diff --git a/src/integralequation.rs b/src/integralequation.rs
@@ -134,7 +134,7 @@ fn rism_vv_equation_impl(
     ck = ck - b * uk_lr.to_owned();
 
     // Perform integral equation calculation in k-space
-    // H = (I - W * C * P)^-1 * (W * C * W)
+    // H = (I - W * C * P * Χ)^-1 * (W * C * W) + Χ
     Zip::from(hk.outer_iter_mut())
         .and(wk.outer_iter())
         .and(ck.outer_iter())
@@ -146,6 +146,7 @@ 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;