different parameters for PMV

pyrism/data/1.toml

@@ -4,7 +4,7 @@ kT = 1.0
 kU = 0.00198720414667
 charge_coeff = 167101.0
 npts = 16384
-radius = 10.24
+radius = 40.96
 lam = 1
 
 [params]

pyrism/data/Methanol_wat.toml

@@ -1,11 +1,11 @@
 #Taken from Lee and Maggiora 1993
 
 [system]
-temp = 298.15
+temp = 300
 kT = 1.0
 kU = 0.00198720414667
 charge_coeff = 167101.0
-npts = 2048
+npts = 16384
 radius = 20.48
 lam = 1
 
@@ -16,6 +16,7 @@ IE = "XRISM"
 solver = "MDIIS"
 depth = 16
 picard_damping = 0.5
+mdiis_damping = 0.5
 itermax = 10000
 tol = 1E-9
 diel = 78.497
@@ -25,24 +26,23 @@ adbcor = 1.475
 nsv = 3
 nspv = 1
 
-#Pettitt and Rossky TIP3P water model
 [solvent.water]
 dens = 0.03334
 ns = 3
- "O_w" = [
-     [76.48937, 3.15, -0.834],
-     [0.0, 0.0, 0.0]
- ]
+"O" = [
+    [78.2526594071, 3.166, -0.8476], 
+    [0.00000000e+00, 0.00000000e+00, 0.00000000e+00]
+]
 
- "H1_w" = [
-     [23.14810, 0.4, 0.417],
-     [9.57200000e-01, 0.0, 0.0]
- ]
+"H1" = [
+    [23.1635928747, 0.8, 0.4238],
+    [1.00000000e+00, 0.00000000e+00, 0.00000000e+00]
+] 
 
- "H2_w" = [
-     [23.14810, 0.4, 0.417],
-     [-2.39988000e-01, 9.26627000e-01, 0.00000000e+00]
- ]
+"H2" = [
+    [23.1635928747, 0.8, 0.4238],
+    [-3.33314000e-01, 9.42816000e-01, 0.00000000e+00]
+]
 
 [solute]
 nsu = 6

pyrism/data/RISM-MOL_methane.toml

@@ -0,0 +1,56 @@
+[system]
+temp = 300
+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 = 15
+picard_damping = 0.1
+mdiis_damping = 0.6
+itermax = 10000
+tol = 1E-8
+diel = 78.497
+adbcor = 1.5
+
+[solvent]
+nsv = 3
+nspv = 1
+
+[solute]
+nsu = 5
+nspu = 1
+
+[solvent.water]
+dens = 0.03334
+ns = 3
+"O" = [
+    [78.2526594071, 3.166, -0.8476], 
+    [0.00000000e+00, 0.00000000e+00, 0.00000000e+00]
+]
+
+"H1" = [
+    [23.1635928747, 0.8, 0.4238],
+    [1.00000000e+00, 0.00000000e+00, 0.00000000e+00]
+] 
+
+"H2" = [
+    [23.1635928747, 0.8, 0.4238],
+    [-3.33314000e-01, 9.42816000e-01, 0.00000000e+00]
+]
+
+[solute.methane]
+dens = 0.0
+ns = 5
+C1-0 = [ [ 33.2347202115, 3.5, -0.24,], [ -0.4458, -0.0120, 0.0,],]
+H1-1 = [ [ 15.1066910052, 2.5, 0.06,], [ -0.0823, -1.0402, 0.0,],]
+H2-2 = [ [ 15.1066910052, 2.5, 0.06,], [ -0.0823, 0.5020, 0.8904,],]
+H3-3 = [ [ 15.1066910052, 2.5, 0.06,], [ -0.0823, 0.5020, -0.8904,],]
+H4-4 = [ [ 15.1066910052, 2.5, 0.06,], [ -1.5363, -0.0120, 0.0,],]

pyrism/rism_ctrl.py

@@ -625,30 +625,36 @@ def isothermal_compressibility(self, dat):
     def kb_partial_molar_volume(self):
 
         uv = self.uv
+        vv = self.vv
 
         ck = np.zeros((uv.npts, uv.ns1, uv.ns2), dtype=np.float64)
-        #hk = np.zeros((uv.npts, uv.ns1, uv.ns2), dtype=np.float64)
+        hk_vv = np.zeros((vv.npts, vv.ns1, vv.ns2), dtype=np.float64)
+        hk_uv = np.zeros((uv.npts, uv.ns1, uv.ns2), dtype=np.float64)
+
+        for i, j in np.ndindex(vv.ns1, vv.ns2):
+            hk_vv[..., i, j] = vv.grid.dht((vv.t + vv.c)[..., i, j])
 
         for i, j in np.ndindex(uv.ns1, uv.ns2):
-            ck[..., i, j] = uv.grid.dht(uv.c[..., i, j])
-            #hk[..., i, j] = uv.grid.dht((uv.t + uv.c)[..., i, j])
-        hk = self.uv.h_k
+            ck[..., i, j] = uv.grid.dht(uv.c[..., i, j])   
+            hk_uv[..., i, j] = uv.grid.dht((uv.t + uv.c)[..., i, j])
+        #hk_vv = self.vv.h_k
+        #hk_uv = self.uv.h_k
 
         compres = self.isothermal_compressibility(self.vv)
 
         r = self.uv.grid.ri[:, np.newaxis, np.newaxis]
         ck0 = self.integrate(self.uv.c * r * r, 4.0 * np.pi * self.uv.grid.d_r)
-        rhvv = self.integrate(self.vv.h[:, 0, 0] * r * r, 4.0 * np.pi * self.uv.grid.d_r)
-        rhuv = self.integrate(self.uv.h[:, 0, :] * r * r, 4.0 * np.pi * self.uv.grid.d_r)
-        khvv = np.sum(hk[0,0,0])
-        khuv = np.sum(hk[0,:,0])
+        rhvv = self.integrate(self.vv.h * r * r, 4.0 * np.pi * self.uv.grid.d_r)
+        rhuv = self.integrate(self.uv.h * r * r, 4.0 * np.pi * self.uv.grid.d_r)
+        khvv = np.sum(hk_vv[0, ...])
+        khuv = np.sum(hk_uv[0, ...])
         pv = self.vv.p[0][0]
         pvec = np.diag(self.vv.p)
 
         inv_B = self.uv.kT * self.uv.T
         ck0_direct = np.sum(ck[0, ...] @ self.vv.p)
 
-        return (1.0 / pv) + (rhvv - rhuv) / self.uv.ns1
+        return (1.0 / pv) + (khvv - khuv) / self.uv.ns1
 
     def rism_kb_partial_molar_volume(self):
 

pyrism/tests/test_thermo.py

@@ -7,9 +7,6 @@
 
 mol.do_rism(verbose=True)
 
-#plt.plot(mol.vv.grid.ri, mol.uv.h[:, 0, 0])
-#plt.show()
-
 print("isothermal compressibility (1/A^3):", (mol.isothermal_compressibility(mol.vv)))
 
 pressure, pressure_plus = mol.pressure()