source: sasmodels/sasmodels/models/be_polyelectrolyte.py @ c1e44e5

Last change on this file since c1e44e5 was c1e44e5, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

Add local link to source files. Refs #1263.

  • Property mode set to 100644
File size: 9.0 KB
Line 
1r"""
2.. note:: Please read the Validation section below.
3
4Definition
5----------
6This model calculates the structure factor of a polyelectrolyte solution with
7the RPA expression derived by Borue and Erukhimovich\ [#Borue]_.  Note however
8that the fitting procedure here does not follow the notation in that reference
9as 's' and 't' are **not** decoupled. Instead the scattering intensity $I(q)$
10is calculated as
11
12.. math::
13
14    I(q) = K\frac{q^2+k^2}{4\pi L_b\alpha ^2}
15    \frac{1}{1+r_{0}^4(q^2+k^2)(q^2-12hC_a/b^2)} + background
16
17    k^2 = 4\pi L_b(2C_s + \alpha C_a)
18
19    r_{0}^2 = \frac{b}{\alpha \sqrt{C_a 48\pi L_b}}
20
21where
22
23$K$ is the contrast factor for the polymer which is defined differently than in
24other models and is given in barns where 1 $barn = 10^{-24}$ $cm^2$.  $K$ is
25defined as:
26
27.. math::
28
29    K = a^2
30
31    a = b_p - (v_p/v_s) b_s
32
33where:
34
35- $b_p$ and $b_s$ are **sum of the scattering lengths of the atoms**
36  constituting the polymer monomer and the solvent molecules, respectively.
37
38- $v_p$ and $v_s$ are the partial molar volume of the polymer and the
39  solvent, respectively.
40
41- $L_b$ is the Bjerrum length (|Ang|) - **Note:** This parameter needs to be
42  kept constant for a given solvent and temperature!
43
44- $h$ is the virial parameter (|Ang^3|) - **Note:** See [#Borue]_ for the
45  correct interpretation of this parameter.  It incorporates second and third
46  virial coefficients and can be *negative*.
47
48- $b$ is the monomer length (|Ang|).
49
50- $C_s$ is the concentration of monovalent salt(1/|Ang^3| - internally converted from mol/L).
51
52- $\alpha$ is the degree of ionization (the ratio of charged monomers to the total
53  number of monomers)
54
55- $C_a$ is the polymer molar concentration (1/|Ang^3| - internally converted from mol/L)
56
57- $background$ is the incoherent background.
58
59For 2D data the scattering intensity is calculated in the same way as 1D,
60where the $\vec q$ vector is defined as
61
62.. math::
63
64    q = \sqrt{q_x^2 + q_y^2}
65
66Validation
67----------
68
69As of the last revision, this code is believed to be correct.  However it
70needs further validation and should be used with caution at this time.  The
71history of this code goes back to a 1998 implementation. It was recently noted
72that in that implementation, while both the polymer concentration and salt
73concentration were converted from experimental units of mol/L to more
74dimensionally useful units of 1/|Ang^3|, only the converted version of the
75polymer concentration was actually being used in the calculation while the
76unconverted salt concentration (still in apparent units of mol/L) was being
77used.  This was carried through to Sasmodels as used for SasView 4.1 (though
78the line of code converting the salt concentration to the new units was removed
79somewhere along the line). Simple dimensional analysis of the calculation shows
80that the converted salt concentration should be used, which the original code
81suggests was the intention, so this has now been corrected (for SasView 4.2).
82Once better validation has been performed this note will be removed.
83
84References
85----------
86
87.. [#Borue] V Y Borue, I Y Erukhimovich, *Macromolecules*, 21 (1988) 3240
88.. [#] J F Joanny, L Leibler, *Journal de Physique*, 51 (1990) 545
89.. [#] A Moussaid, F Schosseler, J P Munch, S Candau, *J. Journal de Physique
90   II France*, 3 (1993) 573
91.. [#] E Raphael, J F Joanny, *Europhysics Letters*, 11 (1990) 179
92
93Authorship and Verification
94----------------------------
95
96* **Author:** NIST IGOR/DANSE **Date:** pre 2010
97* **Last Modified by:** Paul Butler **Date:** September 25, 2018
98* **Last Reviewed by:** Paul Butler **Date:** September 25, 2018
99"""
100
101import numpy as np
102from numpy import inf, pi, sqrt
103
104name = "be_polyelectrolyte"
105title = "Polyelectrolyte with the RPA expression derived by Borue and Erukhimovich"
106description = """
107            Evaluate
108            F(x) = K 1/(4 pi Lb (alpha)^(2)) (q^(2)+k2)/(1+(r02)^(2))
109                 (q^(2)+k2) (q^(2)-(12 h C/b^(2)))
110
111            has 3 internal parameters :
112                   The inverse Debye Length: K2 = 4 pi Lb (2 Cs+alpha C)
113                   r02 =1/alpha/Ca^(0.5) (B/(48 pi Lb)^(0.5))
114                   Ca = 6.022136e-4 C
115            """
116category = "shape-independent"
117
118# pylint: disable=bad-whitespace, line-too-long
119#   ["name", "units", default, [lower, upper], "type", "description"],
120parameters = [
121    ["contrast_factor",       "barns",   10.0,  [-inf, inf], "", "Contrast factor of the polymer"],
122    ["bjerrum_length",        "Ang",      7.1,  [0, inf],    "", "Bjerrum length"],
123    ["virial_param",          "Ang^3", 12.0,  [-inf, inf], "", "Virial parameter"],
124    ["monomer_length",        "Ang",     10.0,  [0, inf],    "", "Monomer length"],
125    ["salt_concentration",    "mol/L",    0.0,  [-inf, inf], "", "Concentration of monovalent salt"],
126    ["ionization_degree",     "",         0.05, [0, inf],    "", "Degree of ionization"],
127    ["polymer_concentration", "mol/L",    0.7,  [0, inf],    "", "Polymer molar concentration"],
128    ]
129# pylint: enable=bad-whitespace, line-too-long
130
131
132def Iq(q,
133       contrast_factor,
134       bjerrum_length,
135       virial_param,
136       monomer_length,
137       salt_concentration,
138       ionization_degree,
139       polymer_concentration):
140    """
141    :params: see parameter table
142    :return: 1-D form factor for polyelectrolytes in low salt
143
144    parameter names, units, default values, and behavior (volume, sld etc) are
145    defined in the parameter table.  The concentrations are converted from
146    experimental mol/L to dimensionaly useful 1/A3 in first two lines
147    """
148
149    concentration_pol = polymer_concentration * 6.022136e-4
150    concentration_salt = salt_concentration * 6.022136e-4
151
152    k_square = 4.0 * pi * bjerrum_length * (2*concentration_salt +
153                                            ionization_degree * concentration_pol)
154
155    r0_square = 1.0/ionization_degree/sqrt(concentration_pol) * \
156                (monomer_length/sqrt((48.0*pi*bjerrum_length)))
157
158    term1 = contrast_factor/(4.0 * pi * bjerrum_length *
159                             ionization_degree**2) * (q**2 + k_square)
160
161    term2 = 1.0 + r0_square**2 * (q**2 + k_square) * \
162        (q**2 - (12.0 * virial_param * concentration_pol/(monomer_length**2)))
163
164    return term1/term2
165
166Iq.vectorized = True  # Iq accepts an array of q values
167
168def random():
169    """Return a random parameter set for the model."""
170    # TODO: review random be_polyelectrolyte model generation
171    pars = dict(
172        scale=10000, #background=0,
173        #polymer_concentration=0.7,
174        polymer_concentration=np.random.beta(5, 3), # around 70%
175        #salt_concentration=0.0,
176        # keep salt concentration extremely low
177        # and use explicit molar to match polymer concentration
178        salt_concentration=np.random.beta(1, 100)*6.022136e-4,
179        #contrast_factor=10.0,
180        contrast_fact=np.random.uniform(1, 100),
181        #bjerrum_length=7.1,
182        bjerrum_length=np.random.uniform(1, 10),
183        #virial_param=12.0,
184        virial_param=np.random.uniform(-1000, 30),
185        #monomer_length=10.0,
186        monomer_length=10.0**(4*np.random.beta(1.5, 3)),
187        #ionization_degree=0.05,
188        ionization_degree=np.random.beta(1.5, 4),
189    )
190    return pars
191
192demo = dict(scale=1, background=0.1,
193            contrast_factor=10.0,
194            bjerrum_length=7.1,
195            virial_param=12.0,
196            monomer_length=10.0,
197            salt_concentration=0.0,
198            ionization_degree=0.05,
199            polymer_concentration=0.7)
200
201tests = [
202
203    # Accuracy tests based on content in test/utest_other_models.py
204    # Note that these should some day be validated beyond this self validation
205    # (circular reasoning). -- i.e. the "good value," at least for those with
206    # non zero salt concentrations, were obtained by running the current
207    # model in SasView and copying the appropriate result here.
208    #    PDB -- Sep 26, 2018
209    [{'contrast_factor':       10.0,
210      'bjerrum_length':         7.1,
211      'virial_param':          12.0,
212      'monomer_length':        10.0,
213      'salt_concentration':     0.0,
214      'ionization_degree':      0.05,
215      'polymer_concentration':  0.7,
216      'background':             0.001,
217     }, 0.001, 0.0948379],
218
219    [{'contrast_factor':       10.0,
220      'bjerrum_length':       100.0,
221      'virial_param':           3.0,
222      'monomer_length':         5.0,
223      'salt_concentration':     1.0,
224      'ionization_degree':      0.1,
225      'polymer_concentration':  1.0,
226      'background':             0.0,
227     }, 0.1, 0.253469484],
228
229    [{'contrast_factor':       10.0,
230      'bjerrum_length':       100.0,
231      'virial_param':           3.0,
232      'monomer_length':         5.0,
233      'salt_concentration':     1.0,
234      'ionization_degree':      0.1,
235      'polymer_concentration':  1.0,
236      'background':             1.0,
237     }, 0.05, 1.738358122],
238
239    [{'contrast_factor':     100.0,
240      'bjerrum_length':       10.0,
241      'virial_param':         12.0,
242      'monomer_length':       10.0,
243      'salt_concentration':    0.1,
244      'ionization_degree':     0.5,
245      'polymer_concentration': 0.1,
246      'background':           0.01,
247     }, 0.5, 0.012881893],
248    ]
Note: See TracBrowser for help on using the repository browser.