Model for a mixture of two homopolymers (case 0 in rpa model)

sasmodels/models/binary_blend.py

@@ -0,0 +1,144 @@
+r"""
+
+Two-polymer RPA model with a flat background.
+
+Definition
+----------
+This model calculates the scattering from a two polymer
+blend using the Random Phase Approximation (RPA).
+
+This is a revision of the *binary_blend* model posted to the
+Marketplace in May 2020 following User feedback.
+
+.. note:: The two polymers are assumed to be monodisperse and incompressible.
+
+.. note:: There is an equivalent model in sasmodels called *rpa_binary_mixture_homopolymers*.
+          While they do the same calculation (although *rpa_binary_mixture_homopolymers* is
+          coded in C), they use different input parameters (see the validation section in
+          the documentation of *rpa_binary_mixture_homopolymers* for details), so you can
+          choose one or another depending on your preferences.
+
+The scattered intensity $I(q)$ is calculated as [1,2]
+
+.. math::
+
+     \frac{(\rho_A - \rho_B)^2}{N_{Av} I(q)} = \text{scale} \cdot
+     \left[\frac{1}{\phi_A n_A v_A P_A(q)} + \frac{1}{\phi_B n_B v_B P_B(q)} -
+     \frac {2 \chi_{AB}}{v_0}\right] + \text{background}
+
+where $\rho_i$ is the scattering length density, $\phi_i$ is the volume
+fraction (such that $\phi_A$ + $\phi_B$ = 1), $n_i$ is the degree of
+polymerization (number of repeat units), $v_i$ is the molar specific volume
+of one *monomer*, $P_i(q)$ is Debye's Gaussian coil form factor
+for polymer $i$, $N_{Av}$ is Avogadro's number, $\chi_{AB}$ is the Flory-Huggins
+interaction parameter and $v_0$ is the reference volume,
+
+with
+
+.. math::
+
+     v_i = m_i / \delta_i , \;
+     v_0 = \sqrt{v_A \cdot v_B} , \;
+     Z = (q \cdot Rg_i)^2 , \;
+     P_i(q) = 2 \cdot [exp(-Z) + Z - 1] / Z^2 ,
+
+where $m_i$ is the molecular weight of a *repeat unit* (not of the polymer),
+$\delta_i$ is the mass density, and $Rg_i$ is the radius of gyration of polymer $i$.
+
+.. note:: This model works best when as few parameters as possible are allowed
+          to optimize. Indeed, most parameters should be known *a priori* anyhow!
+          The calculation should also be exact, meaning the *scale* parameter
+          should be left at 1.
+
+.. note:: Tip: Try alternately optimizing $n_A$ & $Rg_A$ and $n_B$ & $Rg_B$ and
+          only then optimizing $\chi_{AB}$.
+
+Acknowledgments
+----------------
+The author would like to thank James Cresswell and Richard Thompson for
+highlighting some issues with the model as it was originally coded.
+The molar specific volumes were being computed from the polymer molecular
+weight, not the weight of the repeat unit, meaning in most instances the
+values were grossly over-estimated, whilst the reference volume was fixed at
+a value which in most instances would have been too small. Both issues are
+corrected in this version of the model.
+
+References
+----------
+
+#.  Lapp, Picot & Benoit, *Macromolecules*, (1985), 18, 2437-2441 (Appendix) 
+#.  Hammouda, *The SANS Toolbox*, Chapters 28, 29 & 34 and Section H
+
+Authorship and Verification
+---------------------------
+
+* **Author:** Steve King **Date:** 07/05/2020
+* **Last Modified by:** Steve King **Date:** 12/09/2024
+* **Last Reviewed by:** Miguel A. Gonzalez **Date:** 16/11/2025
+
+"""
+import numpy as np
+from numpy import inf
+
+name = "binary_blend"
+title = "Two-polymer RPA model"
+description = """
+    Evaluates a two-
+    polymer RPA model
+"""
+category = "Polymers"
+structure_factor = False
+single = False
+
+#   ["name", "units", default, [lower, upper], "type","description"],
+#   lower limit for m_A/m_B set for CH repeat unit
+#   lower limit for n_A/n_B set for PS (see Macromolecules, 37(1), 2004, 161-166);
+#   could be reduced for more flexible polymers and vice versa
+#   lower limit for rg_A/rg_B set for C-C bond!
+parameters = [
+    ["chi_AB", "cm^3/mol", 0.001, [0, 100], "", "Flory interaction parameter"],
+    ["density_A", "g/cm^3", 1.05, [0.1, 3], "", "Mass density of A"],
+    ["density_B", "g/cm^3", 0.90, [0.1, 3], "", "Mass density of B"],
+    ["m_A", "g/mol", 112, [13, inf], "", "Mol weight of REPEAT A"],
+    ["m_B", "g/mol", 104, [13, inf], "", "Mol weight of REPEAT B"],
+    ["n_A", "", 465, [50, inf], "", "No. repeats of A"],
+    ["n_B", "", 501, [50, inf], "", "No. repeats of B"],
+    ["rg_A", "Ang", 59.3, [1.5, inf], "", "Radius of gyration of A"],
+    ["rg_B", "Ang", 59.3, [1.5, inf], "", "Radius of gyration of B"],
+    ["sld_A", "1e-6/Ang^2", 6.55, [-1, 7], "sld", "SLD of A"],
+    ["sld_B", "1e-6/Ang^2", 1.44, [-1, 7], "sld", "SLD of B"],
+    ["volfrac_A", "", 0.48, [0, 1], "", "Volume fraction of A"],
+]
+
+def Iq(q, chi_AB, density_A, density_B, m_A, m_B,
+    n_A, n_B, rg_A, rg_B, sld_A, sld_B, volfrac_A):
+
+    N_Av = 6.023E+23             # Avogadro number (mol^-1)
+    v_A = m_A / density_A        # Molar specific volume of A (cm^3 mol^-1)
+    v_B = m_B / density_B        # Molar specific volume of B (cm^3 mol^-1)
+    v_0 = np.sqrt(v_A * v_B)     # Reference volume (cm^3 mol^-1)
+
+    Z_A = (q * rg_A) * (q * rg_A)
+    Z_B = (q * rg_B) * (q * rg_B)
+
+    contrast = (sld_A - sld_B)   # Contrast (Ang^-2)
+
+    Pq_A = 2.0 * (np.exp(-1.0 * Z_A) - 1.0 + Z_A) / (Z_A * Z_A)
+    Pq_B = 2.0 * (np.exp(-1.0 * Z_B) - 1.0 + Z_B) / (Z_B * Z_B)
+
+    U_A = volfrac_A * n_A * v_A * Pq_A
+    U_B = (1.0 - volfrac_A) * n_B * v_B * Pq_B
+    chiterm = (2.0 * chi_AB) / v_0
+    prefactor = 1.0E+20 * (contrast * contrast) / N_Av
+    Inverse_Iq = (1.0 / prefactor) * ((1.0 / U_A) + (1.0 / U_B) - chiterm)
+    result = 1.0 / Inverse_Iq
+    return result
+
+Iq.vectorized = True  # Iq accepts an array of q values
+
+tests = [
+   [{'scale': 1.0, 'background' : 0.08, 'chi_AB': 0.0005, 'density_A': 1.05,
+     'density_B': 0.9, 'm_A': 112, 'm_B': 104, 'n_A' : 400, 'n_B' : 351,
+     'rg_A': 60, 'rg_B': 48, 'sld_A': 6.55, 'sld_B': 1.44, 'volfrac_A': 0.5},
+     [0.002009078240301904, 0.24943453429220586], [49.594719935513766, 0.5713619727360871]],
+   ]

sasmodels/models/rpa_binary_mixture_homopolymers.c

@@ -0,0 +1,59 @@
+double Iq(double q,
+    double phiA, double nA, double vA, double lA, double bA,
+    double nB, double vB, double lB, double bB,
+    double chiAB
+    );
+
+static inline double debye(double X) {
+  if (X < 1e-3) {
+    // 2*(e^-X - 1 + X)/X^2 ≈ 1 - X/3 + X^2/12 - X^3/60 + ...
+    const double X2 = X*X;
+    return 1.0 - X/3.0 + X2/12.0 - X2*X/60.0;
+  } else {
+    // use expm1 for better accuracy when X is small-to-moderate
+    return 2.0*(expm1(-X) + X)/ (X*X);
+  }
+}
+
+double Iq(double q,
+    double phiA,  // VOL FRACTION
+    double nA,    // DEGREE OF POLYMERIZATION
+    double vA,    // SPECIFIC VOLUME
+    double lA,    // SCATT. LENGTH
+    double bA,    // SEGMENT LENGTH
+    double nB,    // DEGREE OF POLYMERIZATION
+    double vB,    // SPECIFIC VOLUME
+    double lB,    // SCATT. LENGTH
+    double bB,    // SEGMENT LENGTH
+    double chiAB  // CHI PARAM
+    )
+{
+  const double q2 = q*q;
+
+  // Volume fraction of polymer B
+  const double phiB = 1.0 - phiA;
+
+  // Calculate Q^2 * Rg^2 for each homopolymer assuming random walk
+  const double XA = q2 * bA * bA * nA / 6.0;
+  const double XB = q2 * bB * bB * nB / 6.0;
+
+  //calculate all partial structure factors Pij and normalize n^2
+  const double PAA = debye(XA);
+  const double PBB = debye(XB);
+  const double S0AA = nA * phiA * vA * PAA;
+  const double S0BB = nB * phiB * vB *PBB;
+
+  const double v0 = sqrt(vA*vB); // Reference volume
+
+  const double v11 = (1.0/S0BB) - (2.0*chiAB/v0);
+
+  const double Nav = 6.022141e+23;
+  const double contrast = (lA/vA -lB/vB) * (lA/vA -lB/vB) * Nav;
+
+  double intensity = contrast * S0AA / (1.0 + v11*S0AA);
+
+  //rescale for units of lA^2 and lB^2 (fm^2 to cm^2)
+  intensity *= 1.0e-26;
+
+  return intensity;
+}

sasmodels/models/rpa_binary_mixture_homopolymers.py

@@ -0,0 +1,122 @@
+r"""
+
+RPA model for a binary mixture of homopolymers
+
+Definition
+----------
+
+Calculates the macroscopic scattering intensity for a
+two polymer blend using the Random Phase Approximation (RPA).
+
+.. note:: There is an equivalent model in sasmodels called *binary_blend*.
+          While they do the same calculation (although *binary_blend* is
+          coded purely in python), they use different input parameters (see
+          the validation section below for details), so you can choose
+          one or another depending on your preferences.
+
+The model is based on the papers by Akcasu *et al.* [1] and by
+Hammouda [2], assuming the polymer follows Gaussian statistics such
+that $R_g^2 = n b^2/6$, where $b$ is the statistical segment length and $n$ is
+the number of statistical segment lengths. A nice tutorial on how these are
+constructed and implemented can be found in chapters 28, 31 and 34, and Part H,
+of Hammouda's 'SANS Toolbox' [3].
+
+The scattered intensity $I(q)$ is calculated as[2,3]
+
+.. math::
+
+     \frac{\left(l_A/v_A - l_B/v_B\right)^2 N_{Av}}{I(q)} = \text{scale} \cdot
+     \left[\frac{1}{\phi_A n_A v_A P_A(q)} + \frac{1}{\phi_B n_B v_B P_B(q)} -
+     \frac{2 \chi_{AB}}{v_0}\right] + \text{background}
+
+where $l_i$, $v_i$, and $b_i$ are the scattering length,
+molar specific volume and segment length of a single *monomer*,
+$\phi_i$ is the volume fraction (noting that $\phi_A$ + $\phi_B$ = 1),
+$n_i$ is the degree of polymerization (number of repeat units),
+$N_{Av}$ is Avogadro's number, $\chi_{AB}$ is the Flory-Huggins
+interaction parameter, $v_0 = \sqrt{v_A \cdot v_B}$ is the
+reference volume, and $P_i(q)$ is Debye's Gaussian coil form
+factor for an isolated chain of polymer $i$:
+
+.. math::
+
+     P_i(q) = \frac{2 \cdot \left[exp(-q^2 \cdot Rg_i^2) + q^2 \cdot Rg_i^2 - 1\right]}
+     {\left(q \cdot Rg_i\right)^4}
+
+.. note::
+    * In general, the degrees of polymerization, the volume
+      fractions, the molar volumes and the scattering lengths for each
+      component are obtained from other methods and held fixed, while the *scale*
+      parameter should be held equal to unity.
+    * The variables to fit are normally the segment lengths $b_a$ and $b_b$,
+      and the Flory-Huggins interaction parameter $\chi_{AB}$.
+
+Validation
+----------
+
+The default parameters correspond to the example case of a binary blend
+of deuterated polysterene (PSD) and poly (vinyl methyl ether) (PVME) at 60 C,
+reported in section 8.2 of [2]. Thus, the generated scattering curve
+can be compared with the one shown in Fig. 11 of [2].
+
+The same result is also obtained with the *binary_blend* model, which
+uses alternative input entries:
+
+  * Mass density and molecular weight instead of molar volume (use 1.12 $g/cm^3$
+    and 112 g/mol for PSD and 1.05 $g/cm^3$ and 58 g/mol for PVME).
+
+  * Radius of gyration instead of the segment length of the monomer (use 136.3 |Ang| for
+    PSD and 128.2 |Ang| for PVME).
+
+  * Scattering length densities instead of scattering lengths (use 6.42 and 0.359
+    (in units of $10^{-6}/A^2$) for PSD and PVME, respectively).
+
+.. warning:: The original *rpa* model applied to case 0 did not produce the same
+             result. However, at present is not completely clear if this is due to
+             a problem with the code or to the problems with the interface to select
+             and pass the correct parameters.
+
+References
+----------
+
+#. A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993) 4136
+#. B. Hammouda, *Advances in Polymer Science* 106 (1993) 87
+#. B. Hammouda, *SANS Toolbox*
+   https://www.nist.gov/system/files/documents/2023/04/14/the_sans_toolbox.pdf
+   (accessed 15 November 2025).
+
+Authorship and Verification
+----------------------------
+
+* **Author:** Boualem Hammouda - NIST IGOR/DANSE **Date:** pre 2010
+* **Converted to sasmodels by:** Paul Kienzle **Date:** July 18, 2016
+* **Single case rewritten by:** Miguel A. Gonzalez **Date:** November 16, 2025
+"""
+
+from numpy import inf
+
+name = "rpa_binary_mixture_homopolymers"
+title = "Random Phase Approximation for a binary mixture of homopolymers"
+description = """
+Intensity of a binary mixture of homopolymers calculated within the RPA
+"""
+category = "Polymers"
+
+#   ["name", "units", default, [lower, upper], "type","description"],
+parameters = [
+    ["phiA", "", 0.484, [0, 1], "", "volume fraction of polymer A"],
+    ["nA", "", 1741.0, [1, inf], "", "Degree of polymerization of polymer A"],
+    ["vA", "cm^3/mol", 100.0, [0, inf], "", "molar volume of polymer A"],
+    ["lA", "fm", 106.5, [-inf, inf], "", "scattering length of polymer A"],
+    ["bA", "Ang", 8.0, [0, inf], "", "segment length of polymer A"],
+    ["nB", "", 2741.0, [1, inf], "", "Degree of polymerization of polymer B"],
+    ["vB", "cm^3/mol", 55.4, [0, inf], "", "molar volume of polymer B"],
+    ["lB", "fm", 3.3, [-inf, inf], "", "scattering length of polymer B"],
+    ["bB", "Ang", 6.0, [0, inf], "", "segment length of polymer B"],
+    ["chiAB", "cm^3/mol", -0.021, [-inf, inf], "", "A:B interaction parameter"],
+]
+
+source = ["rpa_binary_mixture_homopolymers.c"]
+single = False
+
+# TODO: no random parameters generated for RPA