working on pressure impl

pyrism/data/cSPCE_XRISM_ethene.toml

@@ -0,0 +1,57 @@
+[system]
+temp = 298.15
+kT = 1.0
+kU = 0.0019872158728366637
+charge_coeff = 167101.0
+npts = 16384
+radius = 409.6
+lam = 1
+
+[params]
+potential = "LJ"
+closure = "KH"
+IE = "XRISM"
+solver = "MDIIS"
+depth = 12
+picard_damping = 0.3
+mdiis_damping = 0.3
+itermax = 10000
+tol = 1E-5
+diel = 78.497
+adbcor = 1.5
+
+[solvent]
+nsv = 3
+nspv = 1
+
+[solute]
+nsu = 6
+nspu = 1
+
+[solvent.water]
+dens = 0.03334
+ns = 3
+"O" = [
+    [78.15, 3.1657, -0.8476000010975563], 
+    [0.00000000e+00, 0.00000000e+00, 0.00000000e+00]
+]
+
+"H1" = [
+    [7.815, 1.1657, 0.4238],
+    [1.00000000e+00, 0.00000000e+00, 0.00000000e+00]
+] 
+
+"H2" = [
+    [7.815, 1.1657, 0.4238],
+    [-3.33314000e-01, 9.42816000e-01, 0.00000000e+00]
+]
+
+[solute.ethene]
+dens = 0.0
+ns = 6
+C = [ [ 43.276879937664646, 1.6998347539724155, -0.21799999999999997,], [ -0.667, 0.0, 0.0,],]
+C1 = [ [ 43.276879937664646, 1.6998347539724155, -0.21799999999999997,], [ 0.667, 0.0, 0.0,],]
+H = [ [ 7.548293033434298, 1.2998212293675426, 0.10899999999999999,], [ -1.221, -0.929, 0.071,],]
+H1 = [ [ 7.548293033434298, 1.2998212293675426, 0.10899999999999999,], [ -1.221, 0.929, -0.071,],]
+H2 = [ [ 7.548293033434298, 1.2998212293675426, 0.10899999999999999,], [ 1.221, 0.929, -0.071,],]
+H3 = [ [ 7.548293033434298, 1.2998212293675426, 0.10899999999999999,], [ 1.221, -0.929, 0.071,],]

pyrism/rism_ctrl.py

@@ -658,6 +658,26 @@ def __virial_pressure(self):
 
         return (self.vv.kT * self.vv.T * np.sum(self.vv.p)) - self.vv.grid.d_r * 2.0 / 3.0 * np.pi * np.power(np.sum(self.vv.p), 2) * np.sum(pressure)
 
+
+    def pressure(self):
+        uv = self.vv
+        nu = uv.ns1
+        nv = uv.ns2
+
+        p0 = uv.p[0][0]
+
+        B = 1.0 / uv.T / 1.380658E-23
+
+        ck = np.zeros((uv.npts, uv.ns1, uv.ns2), dtype=np.float64)
+
+        for i, j in np.ndindex(uv.ns1, uv.ns2):
+            ck[..., i, j] = uv.grid.dht(uv.c[..., i, j])
+
+        zk = ((nu - 1) / p0) + ck[0, ...]
+
+        pressure = nu * p0 / B - (np.sum(zk) / B / 2.0)
+
+        return pressure * 1e24
 @jit
 def build_Ur_impl(npts, ns1, ns2, sr_pot, mix, cou, atoms1, atoms2, r, charge_coeff, lam=1):
     """Tabulates full short-range and Coulombic potential

pyrism/tests/test_thermo.py

@@ -1,7 +1,7 @@
 from pyrism.rism_ctrl import *
 import matplotlib.pyplot as plt
 
-mol = RismController("../data/cSPCE_XRISM.toml")
+mol = RismController("../data/cSPCE_XRISM_ethene.toml")
 
 mol.initialise_controller()
 
@@ -11,3 +11,5 @@
 #plt.show()
 
 print("{:.5E}".format(mol.isothermal_compressibility(mol.vv)))
+
+print(mol.pressure())