- Categories
- Sphere
- poly_mic_exv_disph
- poly_mic_exv_disph.py
poly_mic_exv_disph - poly_mic_exv_disph.py
r"""
Definition
----------
For block copolymer micelles, the small angle scattering tends to show features from the particles and the polymer chains on different length scales, via different regions in q. A fundamental framework for describing both the polymer and particle scattering was developed by Jan Skov Pedersen and colleagues. The model framework can be adapted for differences in the shape of the nanoparticle (e.g. spherical, ellipsoidal, cylindrical), the concentration profile of the corona (or shell), and the types of polymer physics describing the chain interactions. In this model, we describe a spherical micelle created from a diblock copolymer; the model specifically includes a defined, though malleable, density profile for the corona and can model both idealized chains (e.g. $\theta$ condition) and polymer chains with excluded volume interactions (inter-chain and intra-chain).
A schematic and other details of the theory can be found in the reference by Rehmann et al, which will be linked below after publication.
The scattering equation is:
\[I(q) = \frac{\phi_{poly}}{N_{agg} \cdot \ v_{poly,total}}\lbrack N_{agg}^{2} \cdot \beta_{core}^{2} \cdot \left( {\ A}_{core}\left( q,\ r_{core}\ \right) \cdot e^{\frac{{- q^{2} \cdot \left( \frac{t_{itf}}{4} \right)}^{2}}{2}} \right)^{2}\]
\[+ \ 2 \cdot N_{agg}^{2}{\cdot \beta}_{core}{{\cdot A}_{core}\left( q,\ r_{core} \right) \cdot \beta_{corona} \cdot A}_{corona}\left( q,\rho_{corona\ profile} \right)e^{- q^{2} \cdot \left( \frac{t_{itf}}{4} \right)^{2}}\]
\[+ \ N_{agg} \cdot \left( N_{agg} - P_{chain,eff}(q = 0) \right) \cdot \beta_{corona}^{2}\left( A_{corona}\left( q,\rho_{corona\ profile}\ \right) \cdot e^{\frac{{- q^{2} \cdot \left( \frac{t_{itf}}{4} \right)}^{2}}{2}} \right)^{2}\]
\[+ \ N_{agg} \cdot \beta_{corona}^{2} \cdot P_{chain,\ eff}\left( q,\nu,\ \xi_{eff},\ v_{cc} \right)\rbrack\]
Each of the terms is described below.
$\phi_{poly}$ is the volume fraction of the total polymer, usually described by scale in the SasView model if using absolute units, or incorporated into the scale if using arbitrary units.
$N_{agg}$ is the aggregation number or the number of polymer chains per one micelle. Note that this is not a fit parameter, but is calculated from the core radius ($r_{core}$), the volume fraction of solvent in the core ($f_{core}$), and volume of the core polymer ($v_{poly,core}$) using:
\[N_{agg} = \frac{4}{3}\pi\cdot\frac{r^3_{core}}{v_{poly,core}}\cdot (1 - f_{core})\]
$v_{poly,total}$ is the volume of a single polymer chain, generally calculated from the volume of each block (e.g. = $v_{poly,core} + v_{poly,\ corona}$). An estimate of this volume can be calculated using $v_{poly,i} = M_{i} \cdot 10^{24}/(\rho_{m,i} \cdot N_{Avo})$, in which $M_i$ is the molar mass and $\rho_{m,i}$ is the mass density of the polymer block. These volumes are set parameters in the SasView code and should not be fit, like the SLDs.
$\beta_i$ is the excess scattering length where $i$ denotes either the core or corona and is defined as $\beta_i =v_i\cdot(\rho_{SLD,i} -\rho_{SLD,0})$. $\rho_{SLD,0}$ is the solvent scattering length density.
$A_{core}$ is the amplitude of the "shape function" $F(q)$ for the core. Since we are assuming this is a spherical particle, it is:
\[A_{core} = \ 3 \cdot \frac{j_{1}\left( q \cdot r_{core} \right)}{q \cdot r_{core}}\]
$t_{itf}$ is the thickness of the interfacial region between the core and corona; the exponential function allows for a smooth decay between the scattering length density profiles of the core and the corona. For most systems, this term is close to 0 and a reasonable starting case is to set the parameter to 0.
$A_{corona}$ is the amplitude of the "shape function" ($F(q)$) for the corona. It describes the scattering from the particle aspects of the corona. It is defined as:
\[A_{corona}\left( q,\rho_{corona\ profile}\ \right) = \frac{4 \cdot \pi\int_{}^{}{\rho_{corona,\ profile}(r) \cdot j_{0}(q \cdot r) \cdot r^{2}dr\ \ }}{4 \cdot \pi\int_{}^{}{\rho_{corona,\ profile}(r) \cdot r^{2}dr}}\]
The corona profile can be whatever density profile one wants, but in this specific code, we use the combination of 2 cubic b-splines below:
\[\rho_{corona,\ \ profile} = (1 - \alpha) \cdot \rho_{1}\left( r,\ r_{core},\ s \right) + \alpha \cdot \rho_{2}(r,\ r_{core},\ s)\]
for \(r_{core} \leq r < r_{core} + s\):
\[\rho_{1}(r) = \frac{4\left( r - \left( r_{core} + s \right) \right)^{3}\ - \left( r\ - \left( r_{core} + 2s \right) \right)^{3}}{4s^{3}}\]
\[\rho_{2}(r) = \frac{- \left( r - \left( r_{core} + s \right) \right)^{3}\ }{4s^{3}}\]
For \(r_{core} + s \leq r < r_{core} + 2s\):
\[\rho_{1}(r) = \frac{- \left( r\ - \left( r_{core} + 2s \right) \right)^{3}}{4s^{3}}\]
Elsewhere:
\[\rho_{1}(r) = 0\]
\[\rho_{2}(r) = 0\]
Note that the profile ends at $r_{core} + 2s$, the thickness of the corona is $2s$, and $\alpha$ is associated with steepness of the decay of the volume fraction profile.
$P_{chain,eff}$ is the form factor of the effective polymer chain. This describes the scattering of the polymer chains at high q. It is calculated with:
\[P_{chain,eff} = \dfrac{P_{exv}}{1 + v_{cc}\cdot P_{exv}} \]
$v_{cc}$ describes the effective concentration of chains in the corona. It allows the chain scattering to be modulated as if the chain is in a semi-dilute solution; as concentration increases, $v_{cc}$ increases, and the chain scattering contributes less to the polymer micelle scattering. To model dilute chain scattering, simply set $v_{cc}= 0$.
The excluded volume polymer chain scattering has two parameters which describe it: \(\nu\) is the excluded volume coefficient and \(\xi_{eff}\) is the effective size of the polymer. This scattering function is similar to the $excluded\_volume\_polymer$ model in SasView, although the effective polymer size (a diameter) is used here instead of the radius of gyration of the polymer chain (more details in the reference).
\[P_{exv}(q) = \frac{1}{\nu \cdot U^{\ \frac{1}{2} \cdot \nu}\ \ }\left\lbrack \gamma\left( \frac{1}{2\nu},\ U \right) - \frac{1}{U^{\ \frac{1}{2} \cdot \nu}} \cdot \gamma\left( \frac{1}{\nu},\ U \right) \right\rbrack\]
\[\gamma(x,\ U) = \int_{0}^{U}{e^{- t} \cdot t^{\ x - 1}\ dt\ }\ \]
\[U = q^{2} \cdot \left( \frac{\xi_{eff}}{2} \right)^{2} \cdot \frac{(2\nu + 1)(2\nu + 2)}{6} = q^{2} \cdot \frac{b_{k}^{2}n_{k}^{2 \cdot \nu}}{6}\]
in which \(\gamma(x,\ U)\) is the incomplete $\gamma$ function (defined in the interval above); U is a scattering variable defined as either a function of the effective polymer size (\(\xi_{eff}\)) and \(\nu\) or as a function of \(\nu\), the Kuhn length $b_k$, and the number of Kuhn segments $n_k$.
$P_{exv}(q)$ in this form can describe different solvent conditions through $\nu$: ideal solvents (e.g. $\theta$ condition) are $\nu = \frac{1}{2}$; poor solvents are $\nu = \frac{1}{3}$, and good solvents are described when $\nu = \frac{3}{5}$. This describes the scattering at high q as $I(q) \propto q^{-1/\nu}$. To model an ideal chain, simply set \(\nu = \frac{1}{2}\) to recapitulate the Deybe scattering equation for polymers with Gaussian chain statistics. This model allows one to fit the value for \(\nu\), but there tends to be issues with correlations between \(\nu,\ \xi_{eff}\), and the background. We recommend choosing $\nu$ and not fitting it.
Additional Practical Considerations
Due to the number of parameters and other considerations, we strongly suggest using this model with the advanced DREAM algorithm. In doing so, we suggest the following initialization parameters:
samples >= 1.5e6
this can be a little less if you use the burn feature, which we generally set to 1e4. The aim is to have > 1e6 samples after the convergence criteria is met.
initializer = lhs
i.e. latin hyper square, which uses the parameter range you defined for the fit to choose randomized starting parameters. Be sure to set the ranges for fit parameters to reasonable values.
One should also take a careful look at the parameter traces after the fit: there may be outlier chains in the result that require trimming after the fact or more than one set of parameters with similar likelihoods (see the reference for more information about bifurcation with this model). These features are easier to access with the sasmodels program than in the SasView application GUI, where one can re-set the fit, use burn after the fact, or select a subset of results.
References
----------
1. K. Rehmann and K. Weigandt *Generalized Form Factor Model for Block Copolymer Micelles* (2026) submitted.
2. J.S. Pedersen; C. Svaneborg; K. Almdal; I.W. Hamley; R.N. Young, *A Small-Angle Neutron and X-ray Contrast Variation Scattering Study of the Structure of Block Copolymer Micelles: Corona Shape and Excluded Interactions* Macromolecules 2003 36 (2), 416-433 DOI: 10.1021/ma0204913
3. B. Hammouda *SANS from homogeneous polymer mixtures: a unified overview.* chapter in: "Polymer Characteristics. Advances in Polymer Science," (1993) https://doi.org/10.1007/BFb0025862
4. C. Svaneborg and J.S. Pedersen *Block copolymer micelle coronas as quasi-two-dimensional dilute of semi-dilute polymer solutions* Physical Review E Vol 54 (2001)
5. J.S. Pedersen and M.C. Gerstenberg *The structure of P85 Pluronic block copolymer micelles determined by small-angle neutron scattering* Colloids and Surfaces A (2003)
Validation
----------
References
----------
1. K. Rehmann and K. Weigandt *Generalized Form Factor Model for Block Copolymer Micelles* (2026) in prep.
2. J.S. Pedersen; C. Svaneborg; K. Almdal; I.W. Hamley;
R.N. Young, *A Small-Angle Neutron and X-ray Contrast Variation Scattering
Study of the Structure of Block Copolymer Micelles: Corona Shape and Excluded
Interactions* Macromolecules 2003 36 (2), 416-433
DOI: 10.1021/ma0204913
3. B. Hammouda *SANS from homogeneous polymer mixtures: a unified overview.*
chapter in: "Polymer Characteristics. Advances in Polymer Science," (1993)
https://doi.org/10.1007/BFb0025862
4. C. Svaneborg and J.S. Pedersen *Block copolymer micelle coronas as quasi-two-dimensional dilute of
semi-dilute polymer solutions* Physical Review E Vol 54 (2001)
5. J.S. Pedersen and M.C. Gerstenberg *The structure of P85 Pluronic block copolymer micelles determined by small-angle neutron scattering* Colloids and Surfaces A (2003)
Authorship and Verification
----------------------------
* **Author:** Kelsi M.S. Rehmann
* **Last Modified by: Kelsi M.S. Rehmann
* **Last Reviewed by:** Reviewer Name Here **Date:** Date Here
"""
import numpy as np
from numpy import inf
# Added imports below:
import sasmodels
# from sasmodels.special import (sas_3j1x_x , sas_gammainc, sas_sinx_x,
# sas_gamma)
from sasmodels.special import (sas_3j1x_x , sas_sinx_x, )
import scipy
from scipy.special import (gamma as sas_gamma, gammainc as sas_gammainc)
name = "poly_mic_exv_disph"
title = ('An updated spherical, diblock copolymer micelle model that '
+'includes excluded volume interactions for the chains in the corona '
+ 'and specifies the corona density profile')
description = """\
An updated spherical, diblock copolymer micelle model with excluded volume interactions and a specific corona density
I(q) = (scale/N_agg * v_poly,tot) * {
[N_agg * beta_core * A_core(q, r_core)* exp(-q^2 * (t_itf/4)^2/2)]^2
+[ N_agg *beta_corona*P_chain,eff(q, nu, xi_eff)]
+ N_agg*(N_agg - P_chain,eff(q=0)) * [beta_cororona * A_cororona(q, r_core, s, alpha)* exp(-q^2 * ((t_itf/4)^2)/2)]^2
+ [2* N_agg^2 * beta_core* A_core* beta_corona *A_corona * exp(-(q^2*(t_itf/4))^2)]
} + background
scale: polymer volume fraction (if I(q) on absolute scale)
N_agg : number of polymer chains in 1 micelle
v_core: volume of 1 core chain
v_corona: volume of 1 corona chain
sld_core: neutron scattering length density of micelle core polymer
sld_corona: neutron scattering length density of micelle corona polymer
sld_solvent: neutron scattering length density of solvent
f_core: volume fraction of solvent in the core
t_itf: Width of the core-corona interface
radius_core: radius of the core
xi_eff: effective blob size of polymer chains in the corona
nu: excluded volume parameter (polymer-solvent interactions)
v_cc: effective chain-chain interaction parameter in the corona
s: width of the corona profile/2, fitting parameter for the radial density profile
alpha : fitting parameter for shape of the radial density profile
"""
category = "plugin"
# single= True # This is in the writing a model document; it defaults to True
structure_factor = False
# ["name", "units", default, [lower, upper], "type", "description"],
parameters = [["v_core", "Ang^3", 2.533e3, [1, inf], "volume",
"volume of 1 core chain"],
["v_corona", "Ang^3", 6.989e3, [1, inf], "volume",
"volume of 1 corona chain"],
# v_polymer can be estimated by molar mass of polymer[g/mol]
# *(10^24 Ang^3/cm^3)/
# (mass density of the polymer[g/cm^3]*N_Avogadro)
# or by using v_monomer* n (n=degree of polymerization)
["sld_core", "1e-6/Ang^2", 1.410, [-inf, inf], "sld",
"neutron scattering length density of micelle core polymer"],
["sld_corona", "1e-6/Ang^2", 0.671, [-inf, inf], "sld",
"neutron scattering length density of micelle corona polymer"],
["sld_solvent", "1e-6/Ang^2", 6.393, [-inf, inf], "sld",
"neutron scattering length density of solvent"],
["f_core", "None", 0, [0, 1], "volume",
"volume fraction of solvent in the core"],
["t_itf", "Ang", 5, [0, inf], "",
"Width of the core-corona interface"] ,
["radius_core", "Ang", 35, [1, inf], "volume",
"radius of the core"],
["xi_eff", "Ang", 100, [1, inf], "volume",
"effective blob size of polymer chains in the corona"],
["nu", "None", 0.5, [0.3, 1], "",
"excluded volume parameter (polymer-solvent interactions)"],
["v_cc", "None", 2.3, [0, inf], "",
"effective chain-chain interaction parameter in the corona"],
["s", "Ang", 100, [1, inf], "volume",
"width of the corona profile/2, " +
"fitting parameter for radial density profile"],
# E.g. R_micelle = radius_core + 2*s
["alpha", "None", 0.2, [0, 1], "",
"fitting parameter for shape of the radial density profile"]
]
# Note: you cannot put preliminary calculations outside of Iq definition
def Iq(q,
v_core, v_corona, sld_core, sld_corona, sld_solvent,
# The parameters above this line should not be fit, unless with a
# tight range and prior expressed as the standard deviation of an
# associated measurement
f_core, t_itf, radius_core, xi_eff, nu, v_cc, s, alpha
):
"""
Parameters:
q: input scattering vectors, units 1/Ang
v_core: volume of 1 core chain
v_corona: volume of 1 corona chain
sld_core: neutron scattering length density of micelle core polymer
sld_corona: neutron scattering length density of micelle corona polymer
sld_solvent: neutron scattering length density of solvent
f_core: volume fraction of solvent in the core
t_itf: Width of the core-corona interface
radius_core: radius of the core
xi_eff: effective blob size of polymer chains in the corona
nu: excluded volume parameter (polymer-solvent interactions)
v_cc: effective chain-chain interaction parameter in the corona
s: width of the corona profile/2, fitting parameter for
radial density profile
alpha : fitting parameter for shape of the radial density profile]
Returns:
I(q): 1D scattering intensity at q, units 1/cm
"""
# Section 1 : Set-up
# 1.1 Preliminary Calculations
# Based on volume of polymer in core = N_agg*v_core and = 4*Pi/3 *r_core^3/(1-f_core)
N_agg = ((4.0*np.pi/3.0)*(radius_core)**3)*((1-f_core)/v_core)
# Excess scattering length, beta
beta_core = v_core*(sld_core-sld_solvent) #1e-6 /Ang^2 unit conversion applied at the end
beta_corona = v_corona*(sld_corona-sld_solvent) #1e-6 /Ang^2 unit conversion applied at the end
# Term describing interface of core-corona, e.g. a "smoothly
# decaying density at the surface" with an interfacial width of
# t_itf. Because the term takes the form of a Gaussian in this code,
# the 'interfacial thickness' t_itf is considered to be
# +/- 2 * sigma of the associated Gaussian, which we call sigma_itf.
# The term is convoluted with both the core and corona so that
# there is not a sharp drop in scattering length density where the
# profile is = 0, e.g. at the interface.
sigma_itf = t_itf/4
int_dk = np.exp(-1*np.float_power((q*sigma_itf),2)/2)
# Section 2: Calculate each of the sub-terms for the total model scattering
# 2.1 A_corona is the form factor amplitude of the radial density profile
# of the corona chains assuming the associated scattering is centro-symmetric;
# specifically, it is "the normalized Fourier Transform of the
# configurationally average radial density
# distribution"
# (ref: Pedersen et al Macromolecules 2003)
# There are two approaches for calculating A_corona: a numeric
# calculation, in which one calls the profile function (rho_bspline
# below) and calculates the associated integrals numerically.
# The second approach uses the analytical equations for the Fourier
# Transform and normalization integrals from the literature.
# Note that the ref. from Pedersen and Gerstenberg and ref. from
# Pedersen ... Young have slightly different forms for the
# 2 cubic bspline model, and thus the FT equations are different.
# The analytical approach below uses the equations from ref. from
# Pedersen and Gerstenberg,
# using the normalization approach described in the
# ref. from Rehmann et al.
# This code uses the analytical approach,
# but can be adapted to use the numerical
# approach if preferred. However, a test of the numerical approach
# took much longer with DREAM (1-2 hr runs became ~30 hr runs).
# If you change the corona profile, we advise determining an
# appropriate analytical result and checking it with
# the numerical result, then changing the A_corona calculation below
# 2.1.1
# The code below can be moved to a Python User Interface / Jupyter
# Notebook to graph the associated radial density profile of
# the corona based on the parameters from this fitting procedure.
# def rho_bspline2_lambda(x, radius_core_ ,s_ ,alpha_ ):
# """
# rho_bspline2_lambda calculates the radial density profile of the corona
# based on the 2 cubic b spline model cited above, using a lambda
# function. This function can be copied into another Python User Interface
# to produce a graph of your profile based on the fit parameters.
# To run:
# Define a set of r's from r_core to r_core + 2s
# r_array= np.linspace(r_core, r_core+2*s_test, spacing)
# Call the profile and the profile*r^2
# # Define the volume of the one repeat unit (v_poly below)
# # and the number of repeat units (z_pzn below)
# Perform the integrals and normalize for volume fraction.
# rho_spline = lambda x: rho_bspline2_lambda(x, r_core, s_test, alpha_test)
# rho_spline_r2 = lambda r: rho_spline(r)*r*r
# C_spline = scipy.integrate.quad(rho_spline_r2 , r_core, r_core+2*s_test)
# data_rhospline_volfrac = (np.vectorize(rho_spline)(r_array) *
# (1/(4*np.pi*C_spline[0]))*(N_agg_test*v_poly*z_pzn))
# """
# if (x >= radius_core_) & (x < (radius_core_ + s_)):
# rho_1= ((4*np.power((x - radius_core_ - s_),3)
# - np.power((x - radius_core_ - 2*s_),3))
# /(4*np.power(s_,3)))
# elif ( (x >= (radius_core_ + s_)) & (x < (radius_core_ + 2*s_)) ):
# rho_1= (-1*np.power((x - radius_core_ - 2*s_),3)
# / (4*np.power(s_,3)))
# else:
# rho_1 =0
# if ( (x >= (radius_core_)) & (x < (radius_core_+s_)) ): # change in Gerstenberg; Bang; Kelley
# rho_2=(-1*np.power((x - radius_core_ - s_ ),3)
# /(4*np.power(s_,3)))
# else:
# rho_2 =0
# rho_chain_ans = (1 - alpha_)*rho_1 + alpha_*rho_2
# return rho_chain_ans
# ============================================================================
# Code below calculates A_corona numerically as part of the
# fitting code; this is much slower!
# def Acorona_num(profile, q_test, r_start, r_end):
# profile is a lambda function of q and r
# pro_j0_r2 = lambda r: profile(r)*(np.sin(q_test*r)/(q_test*r))*r*r
# pro_r2 = lambda r: profile(r)*r*r
# numerator = scipy.integrate.quad_vec(pro_j0_r2, r_start, r_end)[0]
# denominator = scipy.integrate.quad(pro_r2, r_start, r_end)[0]
# # Divided out 4pi/4pi
# Acorona = (numerator/denominator)
# return Acorona
# Rm = radius_core + 2*s
# profile_func = lambda x: rho_bspline2_lambda(x,r_core, s_test, alpha_test)
# A_corona= Acorona_num(profile_func, q, radius_core, Rm)
# =============================================================================
# Analytical approach to calculating the Fourier Transform (FT) of the radial density
# profile in the corona
# Normalization Constants for FT of rho_1 below
Cn1= (np.float_power(s,4)*(15*np.float_power(radius_core,2)
+ 14*radius_core*s
+ 5*np.float_power(s,2)))
V1 = (Cn1/5)/(4*np.float_power(s,3))
# Fourier Transform of rho_1 above
S1 = 24*np.cos(q*(radius_core+2*s))/np.float_power(q,6)
S1 += 6*(radius_core + 2*s)*np.sin(q*(radius_core+2*s))/np.float_power(q,5)
S1 -= 96*np.cos(q*(radius_core + s))/np.float_power(q,6)
S1 -= 24*(radius_core+s)*np.sin(q*(radius_core+s))/np.float_power(q,5)
S1 += (4*(np.float_power(q,4)*radius_core*np.float_power(s,3)
+ (3*np.float_power(q,2)*(radius_core*s))
+ 18)*np.cos(q*radius_core)/ np.float_power(q,6))
S1 -= (2*((2*np.float_power(q,2)*np.float_power(s,3)
- 9*(radius_core-2*s))*np.sin(q*radius_core))/np.float_power(q,5))
S1 = S1/(4*np.float_power(s,3))
# Normalization Constants for FT of rho_2 below
Cn2 = (np.float_power(s,4)*(15*np.float_power(radius_core,2)
+ 6*radius_core*s + np.float_power(s,2)))
V2 = (Cn2/15)/(16*np.float_power(s,3))
# Fourier Transform of rho_2 above
S2 = 96*np.cos(q*(radius_core + s))/(np.float_power(q,6))
S2 += 24*((radius_core+s)*np.sin(q*(radius_core+s)))/(np.float_power(q,5))
# Kelley et al. / Pedersen & Gertsenberg have an s in the above 2 equations
# instead of Bang et al's 2s (likely a typo)
S2 += (4*( (np.float_power(q,4)*radius_core*np.float_power(s,3)
- 6*np.float_power(q,2)*s*(radius_core - s) - 24)*np.cos(q*radius_core))
/(np.float_power(q,6)))
S2 += (4*((np.float_power(q,2)* np.float_power(s,2))*(3*radius_core- s)
-6*(radius_core- 3*s))*np.sin(q*radius_core)/np.float_power(q,5))
S2 = S2/(16*np.float_power(s,3))
# Together
A_corona= ((1-alpha)*S1 + alpha*S2)/((1-alpha)*V1 +alpha*V2)
# =============================================================================
# 3.2 : Self-correlation of Core Term
Term_core = np.float_power(N_agg*beta_core*sas_3j1x_x(q*radius_core)
*int_dk,2)
# 3.3 : Self-correlation of Chains Term
# Form factor of the corona chain(s) with excluded volume effects from
# ref. by Hammouda, similar to polymer_excl_volume model in SasView
# We use xi_eff as an effective blob size (a diameter), but in the modeling
# we use a radius as part of the scattering variable
# (usually Rg for a polymer in solution)
R_eff = xi_eff/2
# U is scattering dependant variable
U = ((q*R_eff)**2)*((2.0*nu+1.0)*(2.0*nu+2.0)/6)
# F_chain_exv = ((1/(nu*U**(1/(2*nu))))
# * sas_gammainc((1/(2*nu)), U) - 1/(nu*U**(1/nu))
# * (sas_gammainc((1/nu), U)))
# In the SasView documentation, it says that SciPy gamma inc is
# regularized, and needs to be re-normalized. So use the function
# below with _sc)
F_chain_exv_sc = ( (1/(nu*np.power(U,(1/(2*nu))))
* sas_gammainc((1/(2*nu)), U)
* sas_gamma((1/(2*nu))))
- ((1/(nu*np.power(U,(1/nu))))
*(sas_gammainc((1/nu), U)*(sas_gamma(1/nu)))) )
# Scale F_chain by v_cc expression, which represents
# increased concentration of chains, and associated chain-chain
# interactions, in the corona
# Note: There is an empirical calculation for v_cc in ref by
# Pedersen.. Young in the appendix and ref. by Svaneborg and Pedersen
# using an estimate of the reduced surface coverage. The authors note
# that the parameters for the power law change based on the methodology
# used, so this should be taken as a first order estimate only
F_chain_exv_sc = F_chain_exv_sc/(1 + v_cc*F_chain_exv_sc)
# Determine the limit of F_chain_exv as q-> 0
# To simplify, choose the first 10 q values, and fit a line to them
intv = 10
# log_qs, log_Fchain = np.log(q[0:intv]), np.log(F_chain_exv_sc[0:intv])
# Fit a line to the log-log scaling of Fchain vs q, then store the intercept
# F_0 = scipy.stats.linregress(log_qs,log_Fchain).intercept
# Revert to the true value from log value
# F_chain_exv_0 = np.exp(F_0)
# This will also work:
F_chain_exv_0 = scipy.stats.linregress(q[0:intv], F_chain_exv_sc[0:intv]).intercept
Term_Chain = N_agg*(beta_corona**2)*F_chain_exv_sc
# 3.4: Cross Term between core and chains
Term_CoreCorona = (2*N_agg**2*beta_core*beta_corona
* (sas_3j1x_x(q*radius_core) * int_dk)*(int_dk * A_corona))
# 3.5: Centro-Symmetric Corona term
Term_CoronaCorona = (N_agg*(N_agg - F_chain_exv_0)
*np.float_power((beta_corona * A_corona* int_dk),2))
# Section 4: Add the terms and normalize
F_mic = Term_core + Term_Chain + Term_CoreCorona + Term_CoronaCorona
# I_micelle needs to be converted from Ang^2 to cm^2, then appropriately
# scaled to get to 1/cm for absolute data. The units of I_micelle
# are "scattering length density ^2", so scattering length/volume of
# scattering object. SasView divides the I(q) below by form_volume
# (units = Ang^3) and assumed scale = volume fraction, such that
# number density = scale/form_volume. Because of this, the
# volume fraction and form_volume are linked; additionally, the
# form_volume is linked to your choice of scattering length density / beta.
# Because beta above assumes we are only using the dry polymer SLD,
# the form_volume must be for the dry polymer volume. As a consequence,
# the scale should be equal volume fraction of dry polymer,
# not the volume fraction of the full micelle (e.g with water
# in the core) as typically used by colloidal scientists.
I = F_mic*10**-4 # Reasonably confident this converts from 1e-12 Ang^2 to 1/cm^2,
# assuming form_volume is in Ang^3
# Section 5 : Adding in structure factor,
# The structure factor proposed by Pedersen et al has the same form as
# the static decoupling approximation (e.g. beta approximation), so an
# additional code can be written to use SasView's implementation of the
# s.d. approximation. Currently (5/2026) this requires the rest of the code
# to be written in C, which will be done at a later date and uploaded as separate code
# in the Model Marketplace. Setting structure_factor = True above will not
# make this code compatible with the structure factor models in SasView.
return I
# End of model
Iq.vectorized = True
"""
Form_volume takes as its input the list of parameters marked as "volume" parameters from
the parameters list above. It does this in the specific order you defined them. You can either run it as
def form_volume(*list_v) where list_v auto-generates said parameter list,
or list them out. If you list them out, it must be in the same order as defined for volume above.
e.e. if v_core is first and then v_corona, then below it must be v_core, v_corona. The list below
is not calling them by name, but by index.
For reference, in this mode, *list_fv is generated in this order:
list_fv = [v_core, v_corona, f_core, radius_core, xi_eff, s]
"""
def form_volume(*list_fv):
"""
Volume of the TOTAL particle shape including all "shells" and
hollow portions (e.g. a hollow sphere like a vesicle). This is
used to compute absolute scattering intensity, to weight
polydisperse parameter contributions and when multiplying by
a structure factor. Note that for hollow particles you will
ALSO need to define a shell volume.
KMSR: This should be dry polymer volume, based on the way we defined the SLDs.
"""
# If using f_core, use one of the equations below
r_core_v = list_fv[3]
N_agg_v = ((4.0*np.pi/3.0)*(r_core_v)**3)*((1-list_fv[2])/list_fv[0])
# N_agg_v =int(N_agg_v + 0.5) -> this results in strange histograms
v_poly_v = list_fv[0] +list_fv[1] #v_core + v_corona
volume = N_agg_v *v_poly_v
return volume
# Tests, assuming scale =1 and background = 0
tests = [
[{"v_core" : 2521.0, "v_corona" : 6960.0, "sld_core": 1.417, "sld_corona": 0.674, "sld_solvent": 6.388, "f_core": 0.0,
"t_itf": 5.0, "radius_core": 35, "xi_eff": 35.6*2, "nu": 0.6, "v_cc": 2, "s": 70, "alpha": 0.4},
[0.001025, 0.00107625, 0.00113006, 0.00118657, 0.00124589, 0.00130819, 0.0013736, 0.00144228, 0.00151439, 0.00159011,
0.002018, 0.0106, 0.5], #qs
[2.05102458e+03, 2.05057919e+03, 2.05008830e+03, 2.04954724e+03, 2.04895096e+03, 2.04829365e+03, 2.04756926e+03, 2.04677093e+03,
2.04589119e+03, 2.04492166e+03, 2.03856574e+03, 1.63941866e+03, 2.19539328e-01] #Iq
],
[{"v_core" : 2521.0, "v_corona" : 6960.0, "sld_core": 1.417, "sld_corona": 0.674, "sld_solvent": 6.388, "f_core": 0.2,
"t_itf": 0.0, "radius_core": 52, "xi_eff": 100, "nu": 0.6, "v_cc": 10, "s": 130, "alpha": 0.8},
[0.001025, 0.00107625, 0.00113006, 0.00118657, 0.00124589, 0.00130819, 0.0013736, 0.00144228, 0.00151439, 0.00159011,
0.002018, 0.0106, 0.5], #qs
[1.40625622e+04, 1.40538756e+04, 1.40443058e+04, 1.40337611e+04, 1.40221484e+04, 1.40093543e+04, 1.39952642e+04, 1.39797472e+04,
1.39626618e+04, 1.39438497e+04, 1.38209617e+04, 7.52708803e+03, 3.16136631e-01]] #Iq
],
Back to Model
Download