implementing diff forms of isotherm comp

2AUK · Jul 4, 2023 · 2ff163f · 2ff163f
1 parent c23c684
commit 2ff163f
Show file tree

Hide file tree

Showing 3 changed files with 27 additions and 9 deletions.
diff --git a/pyrism/data/cSPCE_XRISM.toml b/pyrism/data/cSPCE_XRISM.toml
@@ -16,7 +16,7 @@ depth = 12
 picard_damping = 0.5
 mdiis_damping = 0.2
 itermax = 10000
-tol = 1E-5
+tol = 1E-12
 diel = 78.497
 adbcor = 1.5
 

diff --git a/pyrism/rism_ctrl.py b/pyrism/rism_ctrl.py
@@ -90,7 +90,6 @@ def do_rism(self, verbose=False):
             self.solve(self.vv, dat2=self.uv, verbose=verbose)
         else:
             self.solve(self.vv, dat2=None, verbose=verbose)
-            print(self.__isothermal_compressibility(self.vv))
 
     def read_input(self):
         """ Reads .toml input file, populates vv and uv dataclasses
@@ -583,20 +582,30 @@ def epilogue(self, dat1, dat2=None):
             dat2.h = dat2.t + dat2.c
             self.SFED_calc(dat2, vv=dat1)
 
-    def __isothermal_compressibility(self, dat):
+    def isothermal_compressibility(self, dat):
         # WIP
+        ck = np.zeros((dat.npts, dat.ns1, dat.ns2), dtype=np.float64)
+
+        for i, j in np.ndindex(dat.ns1, dat.ns2):
+            ck[..., i, j] = dat.grid.dht(dat.c[..., i, j])
+
         ck0 = 0.0
         for i in range(dat.grid.npts):
-            ck0r = 0
+            ck0r = 0.0
             for j, k in np.ndindex(dat.ns1, dat.ns2):
                 if j == k:
-                    msym = 1
+                    msym = 1.0
                 else:
-                    msym = 2
-                ck0r += msym*np.diag(dat.p)[j]*np.diag(dat.p)[k]*dat.c[i, j, k]
-            ck0 += ck0r * dat.grid.ri[i] ** 2
+                    msym = 2.0
+                ck0r += np.diag(dat.p)[j]*np.diag(dat.p)[k]*ck[0, j, k]
+            ck0 += ck0r * dat.grid.ri[i] ** 2.0
         ck0 *= 4.0 * np.pi * dat.grid.d_r
-        return dat.B / (np.sum(dat.p) - ck0)
+
+        pck = np.sum(dat.p @ ck)
+        p = np.sum(dat.p)
+
+
+        return (dat.B / (np.sum(dat.p) - ck0), dat.B / (p * (1.0 - pck)))
 
         """
         total_dens = np.sum(dat.p)

diff --git a/pyrism/tests/test_thermo.py b/pyrism/tests/test_thermo.py
@@ -0,0 +1,9 @@
+from pyrism.rism_ctrl import *
+
+mol = RismController("../data/cSPCE_XRISM.toml")
+
+mol.initialise_controller()
+
+mol.do_rism(verbose=True)
+
+print(mol.isothermal_compressibility(mol.vv))