diff --git a/pyrism/data/1.toml b/pyrism/data/1.toml index 9516ecda..bce2834f 100644 --- a/pyrism/data/1.toml +++ b/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] diff --git a/pyrism/data/Methanol_wat.toml b/pyrism/data/Methanol_wat.toml index 523e8e47..1aba856b 100644 --- a/pyrism/data/Methanol_wat.toml +++ b/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 diff --git a/pyrism/data/RISM-MOL_methane.toml b/pyrism/data/RISM-MOL_methane.toml new file mode 100644 index 00000000..00fa9603 --- /dev/null +++ b/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,],] \ No newline at end of file diff --git a/pyrism/rism_ctrl.py b/pyrism/rism_ctrl.py index 8d97abc7..759cca2a 100644 --- a/pyrism/rism_ctrl.py +++ b/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): diff --git a/pyrism/tests/test_thermo.py b/pyrism/tests/test_thermo.py index 4def72fd..aab145ed 100644 --- a/pyrism/tests/test_thermo.py +++ b/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()