dumping thermodynamical data

pyrism/rism_ctrl.py

@@ -317,29 +317,29 @@ def read_input(self):
 
         solver_config = SolverConfig(solver, settings)
 
-        rism_job = RISMDriver(
-            name, data_config, operator_config, potential_config, solver_config
-        )
-        solutions = rism_job.do_rism("very verbose", False)
-
-        grid = Core.Grid(data_config.npts, data_config.radius)
-
-        self.write_solv(
-            name,
-            data_config.solvent_atoms,
-            data_config.solvent_atoms,
-            solutions.vv,
-            grid.ri,
-        )        # literature route
-
-        if solutions.uv:
-            self.write_solu(
-                name,
-                data_config.solute_atoms,
-                data_config.solvent_atoms,
-                solutions.uv,
-                grid.ri,
-            )
+        # rism_job = RISMDriver(
+        #     name, data_config, operator_config, potential_config, solver_config
+        # )
+        # solutions = rism_job.do_rism("very verbose", False)
+        #
+        # grid = Core.Grid(data_config.npts, data_config.radius)
+        #
+        # self.write_solv(
+        #     name,
+        #     data_config.solvent_atoms,
+        #     data_config.solvent_atoms,
+        #     solutions.vv,
+        #     grid.ri,
+        # )        # literature route
+        #
+        # if solutions.uv:
+        #     self.write_solu(
+        #         name,
+        #         data_config.solute_atoms,
+        #         data_config.solvent_atoms,
+        #         solutions.uv,
+        #         grid.ri,
+        #     )
 
         self.name = os.path.basename(self.fname).split(sep=".")[0]
         if "solvent" not in inp:

pyrism/tests/test_thermo.py

@@ -2,9 +2,7 @@
 from pyrism import rust_helpers
 import matplotlib.pyplot as plt
 
-print(rust_helpers.sum_as_string(5, 10))
-
-mol = RismController("../data/cSPCE_XRISM_methane.toml")
+mol = RismController("cSPCE_XRISM_methane.toml")
 
 mol.initialise_controller()
 

src/main.rs

@@ -1,5 +1,6 @@
 use librism::{
     driver::{RISMDriver, Verbosity},
+    thermo::TDDriver,
     writer::RISMWriter,
 };
 use std::path::PathBuf;
@@ -39,16 +40,37 @@ fn main() -> Result<(), lexopt::Error> {
     let args = parse_args()?;
     let mut driver = RISMDriver::from_toml(&args.input_file);
     let solutions = driver.execute(args.verbosity, args.compress);
-    let writer = RISMWriter::new(
-        &args
-            .input_file
-            .file_stem()
-            .unwrap()
-            .to_str()
-            .unwrap()
-            .to_string(),
-        solutions,
+    // let writer = RISMWriter::new(
+    //     &args
+    //         .input_file
+    //         .file_stem()
+    //         .unwrap()
+    //         .to_str()
+    //         .unwrap()
+    //         .to_string(),
+    //     solutions,
+    // );
+    let td = TDDriver::new(solutions);
+    println!(
+        "Isothermal Compressibility: {}",
+        td.isothermal_compressibility()
     );
-    writer.write();
+    println!(
+        "Molecular KB theory PMV: {} (A^3)",
+        td.kb_partial_molar_volume()
+    );
+    println!(
+        "Molecular KB theory PMV: {} (cm^3 / mol)",
+        td.kb_partial_molar_volume() / 1e24 * 6.022e23
+    );
+    println!(
+        "RISM KB theory PMV: {} (A^3)",
+        td.rism_kb_partial_molar_volume()
+    );
+    println!(
+        "RISM KB theory PMV: {} (cm^3 / mol)",
+        td.rism_kb_partial_molar_volume() / 1e24 * 6.022e23
+    );
+
     Ok(())
 }

src/thermo.rs

@@ -1,6 +1,6 @@
-use ndarray::Array;
-
-use crate::solution::Solutions;
+use crate::{data::Grid, solution::Solutions, transforms::fourier_bessel_transform_fftw};
+use ndarray::{s, Array, Axis, Zip};
+use std::f64::consts::PI;
 
 pub struct TDDriver {
     pub solutions: Solutions,
@@ -10,10 +10,78 @@ impl TDDriver {
     pub fn new(solutions: Solutions) -> Self {
         TDDriver { solutions }
     }
+
     pub fn isothermal_compressibility(&self) -> f64 {
         let vv = &self.solutions.vv;
-        let ck = Array::zeros(vv.correlations.cr.raw_dim());
+        let grid = Grid::new(vv.data_config.npts, vv.data_config.radius);
+        let rtok = 2.0 * PI * grid.dr;
+        let mut ck = Array::zeros(vv.correlations.cr.raw_dim());
+
+        Zip::from(vv.correlations.cr.lanes(Axis(0)))
+            .and(ck.lanes_mut(Axis(0)))
+            .par_for_each(|cr_lane, mut ck_lane| {
+                ck_lane.assign(&fourier_bessel_transform_fftw(
+                    rtok,
+                    &grid.rgrid.view(),
+                    &grid.kgrid.view(),
+                    &cr_lane.to_owned(),
+                ));
+            });
+
+        let c_sum = ck.slice(s![0, .., ..]).sum();
+        let p_sum = vv
+            .data_config
+            .solvent_species
+            .iter()
+            .fold(0.0, |acc, x| acc + x.dens);
+
+        1.0 / (p_sum * (1.0 - p_sum * c_sum))
+    }
+
+    pub fn kb_partial_molar_volume(&self) -> f64 {
+        let vv = &self.solutions.vv;
+        let uv = &self.solutions.uv.as_ref().unwrap();
+
+        let p_sum = vv
+            .data_config
+            .solvent_species
+            .iter()
+            .fold(0.0, |acc, x| acc + x.dens);
+
+        let n = uv.data_config.nsu.unwrap() as f64;
+
+        let hkvv_sum = vv.correlations.hk.slice(s![0, .., ..]).sum();
+        let hkuv_sum = uv.correlations.hk.slice(s![0, .., ..]).sum();
+
+        (1.0 / p_sum) + (hkvv_sum - hkuv_sum) / n
+    }
+
+    pub fn rism_kb_partial_molar_volume(&self) -> f64 {
+        let vv = &self.solutions.vv;
+        let uv = &self.solutions.uv.as_ref().unwrap();
+        let grid = Grid::new(vv.data_config.npts, vv.data_config.radius);
+        let rtok = 2.0 * PI * grid.dr;
+        let inv_beta = vv.data_config.temp * vv.data_config.kt;
+        let mut ck = Array::zeros(uv.correlations.cr.raw_dim());
+
+        Zip::from(uv.correlations.cr.lanes(Axis(0)))
+            .and(ck.lanes_mut(Axis(0)))
+            .par_for_each(|cr_lane, mut ck_lane| {
+                ck_lane.assign(&fourier_bessel_transform_fftw(
+                    rtok,
+                    &grid.rgrid.view(),
+                    &grid.kgrid.view(),
+                    &cr_lane.to_owned(),
+                ));
+            });
 
-        todo!()
+        let c_sum = ck.slice(s![0, .., ..]).sum();
+        let p_sum = vv
+            .data_config
+            .solvent_species
+            .iter()
+            .fold(0.0, |acc, x| acc + x.dens);
+        let compressibility = self.isothermal_compressibility();
+        inv_beta * compressibility * (1.0 - p_sum * c_sum)
     }
 }