Setting up to merge to main

Moving the original Python code to a reference folder.
Need to flesh out the python library interface.

reference/pyrism/Closures/__init__.py

@@ -0,0 +1 @@
+from .closure_dispatcher import Closure

reference/pyrism/Closures/closure_dispatcher.py

@@ -0,0 +1,22 @@
+import numpy as np
+from .closure_routines import *
+
+
+class Closure(object):
+    closure_dispatcher = {
+        "HNC": HyperNetted_Chain,
+        "KH": KovalenkoHirata,
+        "PSE-1": PSE_1,
+        "PSE-2": PSE_2,
+        "PSE-3": PSE_3,
+        "PY": PercusYevick,
+        "rHNC": renormalized_HyperNetted_Chain,
+        "rPY": renormalized_PercusYevick,
+        "rKH": renormalized_KovalenkoHirata,
+    }
+
+    def __init__(self, clos):
+        self.closure = clos
+
+    def get_closure(self):
+        return self.closure_dispatcher[self.closure]

reference/pyrism/Closures/closure_object.py

@@ -0,0 +1,39 @@
+#!/usr/bin/env python3
+
+from pyrism.Core import RISM_Obj
+import numpy as np
+from dataclasses import dataclass, field
+
+
+@dataclass
+class ClosureObject:
+    dat_vv: RISM_Obj
+    dat_uv: RISM_Obj = None
+
+    def compute_GF(self, p, energy_unit):
+        if self.dat_uv is not None:
+            return self.__GF_impl(
+                self.dat_uv.t,
+                self.dat_uv.h,
+                self.dat_uv.c,
+                self.dat_uv.B,
+                self.dat_uv.T,
+                self.dat_uv.grid.ri,
+                self.dat_uv.grid.d_r,
+                p,
+                energy_unit,
+            )
+        else:
+            raise RuntimeError(
+                "No free energy functional for solvent-solvent interactions implemented."
+            )
+
+    def __GF_impl(self, t, h, c, B, T, r, dr, p, energy_unit):
+        mu = (
+            4.0
+            * np.pi
+            * p
+            * dr
+            * np.sum(np.power(r, 2)[:, None, None] * ((0.5 * c * h) - c))
+        )
+        return mu / B * energy_unit

reference/pyrism/Closures/closure_routines.py

@@ -0,0 +1,69 @@
+import numpy as np
+from pyrism.Core import RISM_Obj
+
+
+def renormalized_HyperNetted_Chain(data):
+    return (
+        np.exp(-(data.B * data.u_sr) + data.tau + data.Q_r) - 1.0 - data.tau - data.Q_r
+    )
+
+
+def renormalized_PercusYevick(data):
+    return (
+        np.exp(-(data.B * data.u_sr + data.Q_r)) * (1.0 + data.tau)
+        - data.tau
+        - 1.0
+        - data.Q_r
+    )
+
+
+def renormalized_KovalenkoHirata(data):
+    return np.where(
+        (-(data.B * data.u_sr) + data.tau + data.Q_r) <= 0,
+        np.exp(-(data.B * data.u_sr) + data.tau + data.Q_r) - 1.0 - data.tau - data.Q_r,
+        -(data.B * data.u_sr + data.tau + data.Q_r) - data.tau - data.Q_r,
+    )
+
+
+def HyperNetted_Chain(data):
+    return np.exp(-(data.B * data.u_sr) + data.t) - 1.0 - data.t
+
+
+def KovalenkoHirata(data):
+    return np.where(
+        (-(data.B * data.u_sr) + data.t) <= 0,
+        np.exp(-(data.B * data.u_sr) + data.t) - 1.0 - data.t,
+        -(data.B * data.u_sr),
+    )
+
+
+def KobrynGusarovKovalenko(data):
+    zeros = np.zeros_like(data.t)
+    return np.maximum(zeros, -(data.B * data.u_sr))
+
+
+def PercusYevick(data):
+    return np.exp(-(data.B * data.u_sr)) * (1.0 + data.t) - data.t - 1.0
+
+
+def PSE_n(data, n):
+    t_fac = 0
+    for i in range(n):
+        t_fac += np.power((-(data.B * data.u_sr) + data.t), i) / np.math.factorial(i)
+    return np.where(
+        (-(data.B * data.u_sr) + data.t) <= 0,
+        np.exp(-(data.B * data.u_sr) + data.t) - 1.0 - data.t,
+        t_fac - 1.0 - data.t,
+    )
+
+
+def PSE_1(data):
+    return PSE_n(data, 1)
+
+
+def PSE_2(data):
+    return PSE_n(data, 2)
+
+
+def PSE_3(data):
+    return PSE_n(data, 3)

reference/pyrism/Core/Data.py

@@ -0,0 +1,97 @@
+from dataclasses import dataclass, field
+import numpy as np
+from .Grid import Grid
+from .Site import Site
+from .Species import Species
+
+
+@dataclass
+class RISM_Obj(object):
+    # Initial parameters required to instantiate the other attributes
+    T: float
+    kT: float
+    kU: float
+    amph: float
+    ns1: int
+    ns2: int
+    nsp1: int
+    nsp2: int
+    npts: int
+    radius: float
+    nlam: int
+
+    # Set of attributes that are iterated during the RISM calculation
+    c: np.ndarray = field(init=False)
+    c_prev: np.ndarray = field(init=False)
+    t: np.ndarray = field(init=False)
+    h: np.ndarray = field(init=False)
+    g: np.ndarray = field(init=False)
+    h_k: np.ndarray = field(init=False)
+
+    # Set of attributes that remain constant during the RISM calculation
+    B: float = field(init=False)
+    u: np.ndarray = field(init=False)
+    u_sr: np.ndarray = field(init=False)
+    ur_lr: np.ndarray = field(init=False)
+    uk_lr: np.ndarray = field(init=False)
+    w: np.ndarray = field(init=False)
+    p: np.ndarray = field(init=False)
+    grid: Grid = field(init=False)
+    species: list = field(init=False, default_factory=list)
+    atoms: list = field(init=False, default_factory=list)
+    Q_r: np.ndarray = field(init=False)  # XRISM-DB
+    Q_k: np.ndarray = field(init=False)  # XRISM-DB
+    tau: np.ndarray = field(init=False)  # XRISM-DB
+
+    def calculate_beta(self):
+        self.B = 1 / self.T / self.kT
+
+    def __post_init__(self):
+        self.calculate_beta()
+        self.c_prev = np.zeros((self.npts, self.ns1, self.ns2), dtype=np.float64)
+        self.t = np.zeros((self.npts, self.ns1, self.ns2), dtype=np.float64)
+        self.h = np.zeros((self.npts, self.ns1, self.ns2), dtype=np.float64)
+        self.h_k = np.zeros_like(self.h)
+        self.Q_k = np.zeros_like(self.h)
+        self.Q_r = np.zeros_like(self.h)
+        self.tau = np.zeros_like(self.h)
+        self.u = np.zeros((self.npts, self.ns1, self.ns2), dtype=np.float64)
+        self.u_sr = np.zeros_like(self.u)
+        self.ur_lr = np.zeros_like(self.u)
+        self.uk_lr = np.zeros_like(self.u)
+        self.c = np.zeros((self.npts, self.ns1, self.ns2), dtype=np.float64)
+        self.w = np.zeros((self.npts, self.ns1, self.ns1), dtype=np.float64)
+        self.g = np.zeros((self.npts, self.ns1, self.ns2), dtype=np.float64)
+        self.p = np.zeros((self.ns1, self.ns2), dtype=np.float64)
+        self.grid = Grid(self.npts, self.radius)
+
+    def __iter__(self):
+        return iter(
+            (
+                self.T,
+                self.kT,
+                self.kU,
+                self.amph,
+                self.ns1,
+                self.ns2,
+                self.nsp1,
+                self.nsp2,
+                self.npts,
+                self.radius,
+                self.nlam,
+                self.c,
+                self.t,
+                self.h,
+                self.h_k,
+                self.g,
+                self.B,
+                self.u,
+                self.u_sr,
+                self.ur_lr,
+                self.uk_lr,
+                self.w,
+                self.p,
+                self.grid.ri,
+                self.grid.ki,
+            )
+        )

reference/pyrism/Core/Grid.py

@@ -0,0 +1,83 @@
+#!/usr/bin/env python3
+"""
+grid.py
+Defines Grid object to handle generation and transforms
+
+"""
+from dataclasses import dataclass, field
+import numpy as np
+from scipy.fftpack import dst, idst
+from .Transforms import discrete_hankel_transform, inverse_discrete_hankel_transform
+
+
+@dataclass(init=True)
+class Grid:
+    npts: int
+    radius: float
+    ri: np.ndarray = field(init=False)
+    ki: np.ndarray = field(init=False)
+    d_r: float = field(init=False)
+    d_k: float = field(init=False)
+
+    def __post_init__(self):
+        self.ri = np.zeros(self.npts, dtype=float)
+        self.ki = np.zeros(self.npts, dtype=float)
+        self.d_r = self.radius / float(self.npts)
+        self.d_k = 2 * np.pi / (2 * float(self.npts) * self.d_r)
+        self.generate_grid()
+
+    def generate_grid(self):
+        """
+        Generates r-space and k-space grids to compute functions over
+
+        Parameters
+        ----------
+
+        Returns
+        -------
+
+        ri: ndarray
+           r-space grid
+        ki: ndarry
+           k-space grid
+        """
+
+        for i in np.arange(0, int(self.npts)):
+            self.ri[i] = (i + 0.5) * self.d_r
+            self.ki[i] = (i + 0.5) * self.d_k
+
+    def dht(self, fr: np.ndarray) -> np.ndarray:
+        """
+        Discrete Hankel Transform
+
+        Parameters
+        ----------
+
+        fr: ndarray
+           Function to be transformed from r-space to k-space
+
+        Returns
+        -------
+
+        fk: ndarray
+           Transformed function from r-space to k-space
+        """
+        return discrete_hankel_transform(self.ri, self.ki, fr, self.d_r)
+
+    def idht(self, fk: np.ndarray) -> np.ndarray:
+        """
+        Inverse Discrete Hankel Transform
+
+        Parameters
+        ----------
+
+        fk: ndarray
+           Function to be transformed from k-space to r-space
+
+        Returns
+        -------
+
+        fr: ndarray
+           Transformed function from k-space to r-space
+        """
+        return inverse_discrete_hankel_transform(self.ri, self.ki, fk, self.d_k)

reference/pyrism/Core/Site.py

@@ -0,0 +1,9 @@
+import numpy as np
+from dataclasses import dataclass, field
+
+
+@dataclass
+class Site(object):
+    atom_type: str
+    params: list
+    coords: np.ndarray

reference/pyrism/Core/Species.py

@@ -0,0 +1,19 @@
+import numpy as np
+from dataclasses import dataclass, field
+
+
+@dataclass
+class Species(object):
+    species_name: str
+    dens: float = field(init=False)
+    ns: int = field(init=False)
+    atom_sites: list = field(default_factory=list)
+
+    def add_site(self, atom_site):
+        self.atom_sites.append(atom_site)
+
+    def set_density(self, density):
+        self.dens = density
+
+    def set_numsites(self, ns):
+        self.ns = ns