From c66fbe1807f2da1cc0affe0af44582514a1c7c0f Mon Sep 17 00:00:00 2001 From: gonzalezma Date: Mon, 17 Nov 2025 09:32:17 +0100 Subject: [PATCH 1/3] Model for a mixture of two homopolymers (case 0 in rpa model) --- sasmodels/models/binary_blend.py | 144 ++++++++++++++++++ .../models/rpa_binary_mixture_homopolymers.c | 59 +++++++ .../models/rpa_binary_mixture_homopolymers.py | 122 +++++++++++++++ 3 files changed, 325 insertions(+) create mode 100644 sasmodels/models/binary_blend.py create mode 100644 sasmodels/models/rpa_binary_mixture_homopolymers.c create mode 100644 sasmodels/models/rpa_binary_mixture_homopolymers.py diff --git a/sasmodels/models/binary_blend.py b/sasmodels/models/binary_blend.py new file mode 100644 index 00000000..2ecddffc --- /dev/null +++ b/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]], + ] diff --git a/sasmodels/models/rpa_binary_mixture_homopolymers.c b/sasmodels/models/rpa_binary_mixture_homopolymers.c new file mode 100644 index 00000000..de325863 --- /dev/null +++ b/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; +} diff --git a/sasmodels/models/rpa_binary_mixture_homopolymers.py b/sasmodels/models/rpa_binary_mixture_homopolymers.py new file mode 100644 index 00000000..91fe7b33 --- /dev/null +++ b/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.594 + (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", 5.5, [-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 From c0d2447a1259c6f58610d2f766e5d1c12b6aa5f4 Mon Sep 17 00:00:00 2001 From: gonzalezma Date: Mon, 17 Nov 2025 10:50:10 +0100 Subject: [PATCH 2/3] Fix indentation error in documentation --- sasmodels/models/rpa_binary_mixture_homopolymers.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sasmodels/models/rpa_binary_mixture_homopolymers.py b/sasmodels/models/rpa_binary_mixture_homopolymers.py index 91fe7b33..f6e9f2c8 100644 --- a/sasmodels/models/rpa_binary_mixture_homopolymers.py +++ b/sasmodels/models/rpa_binary_mixture_homopolymers.py @@ -73,8 +73,8 @@ .. 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. + a problem with the code or to the problems with the interface to select + and pass the correct parameters. References ---------- From 0ce5481e7e2ccf39e553c3d730c496d24ccc5fbe Mon Sep 17 00:00:00 2001 From: gonzalezma Date: Mon, 17 Nov 2025 13:33:29 +0100 Subject: [PATCH 3/3] Correct scattering length of PVME so that default parameters are ok for PSD/PVME mixture described by Hammouda --- sasmodels/models/rpa_binary_mixture_homopolymers.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sasmodels/models/rpa_binary_mixture_homopolymers.py b/sasmodels/models/rpa_binary_mixture_homopolymers.py index f6e9f2c8..f8e68a92 100644 --- a/sasmodels/models/rpa_binary_mixture_homopolymers.py +++ b/sasmodels/models/rpa_binary_mixture_homopolymers.py @@ -68,7 +68,7 @@ * 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.594 + * 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 @@ -111,7 +111,7 @@ ["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", 5.5, [-inf, inf], "", "scattering length 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"], ]