- Categories
- Sphere
- poly_mic_exv_disph_scl
- poly_mic_exv_disph_scl.py
poly_mic_exv_disph_scl - poly_mic_exv_disph_scl.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. In this specific code, we use the generalized scaling profile of Halperin and multiplied it by an exponential cut off function such that the density $\rho(r)=0$ at $r= r_{core}+ t_{corona}$
\[\rho(r) = r^{\frac{-1(3\cdot\nu -1)}{\nu}}\cdot(1- e^{\frac{ r - (r_{core} + t_{corona})}{\sigma_m \cdot (r_{core} + t_{corona})}} ) \]
when $\nu= \frac{3}{5}$, the exponent is $-\frac{4}{3}$. $\sigma_m$ is effectively a tuning parameter that shifts the contribution of the profile from the scaling behavior to the exponential term. 0.1 gives an appropriate trade-off between the scaling behavior near the core and the drop off near the end of the micelle .
$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 correlations between $v_{cc}= 0$ and $\xi_{eff}$, we recommend either setting $v_{cc}= 0$ and fitting $\xi_{eff}$, or setting $\xi_{eff}$= the radius of gyration of the equivalent homopolymer in solution and fitting $v_{cc}$.
Similarly, 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.
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 *New Implementation of a 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
----------
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_scl"
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 using the generalized scaling'
+ 'relationship of Halperin, modified by KMSR')
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 correlation length of polymer chains in the corona
nu: excluded volume parameter (polymer-solvent interactions)
v_cc: effective chain-chain interaction parameter in the corona
t_corona: width of the corona profile, fitting parameter for the radial density profile
sigma_m : tuning parameter for exponential cut off of scaling profile; do not fit
"""
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", 0, [0, inf], "",
"Width of the core-corona interface"] ,
["radius_core", "Ang", 35, [1, inf], "volume",
"radius of the core"],
["xi_eff", "Ang", 50, [1, inf], "volume",
"effective correlation length 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"],
["t_corona", "Ang", 100, [1, inf], "volume",
"width of the corona profile, " +
"fitting parameter for radial density profile"],
["sigma_m", "None", 0.1, [0.01, 1], "",
"tuning parameter for exponential cut off of scaling profile; do not fit"]
]
# 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, t_corona, sigma_m
):
"""
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 correlation length of polymer chains in the corona
nu: excluded volume parameter (polymer-solvent interactions)
v_cc: effective chain-chain interaction parameter in the corona
t_corona: width of the corona profile, fitting parameter for
radial density profile
sigma_m: 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 (gen_Halperin
# below) and calculates the associated integrals numerically.
# The numerical calculation usually is more computationally time consuming than
# the analytical calculation; in the case of the generalized Halperin profile,
# the analytical calculation includes incomplete Gamma functions
# with imaginary numbers, which are not available in scipy.
# They are available with mpmath, which does not use vectorization, thus slowing
# down the calculation considerably. Since the functions used for
# the numerical calculation use scipy and numpy with vectorization,
# that approach is faster. There may be a way to approximate the calculation
# to increase the speed.
# 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 gen_Halperin(x, r_core, t_corona1, nu1, sigma_m1):
exp_pow = -1*(3*nu1 -1)/nu1
r_total = r_core + t_corona1
rho_corona = x**(exp_pow)*(1-(np.exp((x-r_total)/(sigma_m1*r_total))))
return rho_corona
# Code below calculates A_corona numerically as part of the
# fitting code
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 + t_corona
profile_func = lambda x: gen_Halperin(x, radius_core, t_corona, nu, sigma_m)
A_corona= Acorona_num(profile_func, q, radius_core, Rm)
# =============================================================================
# 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 correlation length, considered a radius here
# because it can be on the order of Rg for a polymer in solution.
# U is scattering dependant variable
U = ((q*xi_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 that equation 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
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, t_corona]
"""
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
Back to Model
Download