source: sasmodels/sasmodels/models/be_polyelectrolyte.py @ 284bdd4

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 284bdd4 was 284bdd4, checked in by butler, 13 months ago

Fix docs and code for be_polyelectrolyte.py

Docs were fixed as explained in ticket #1155 as was the code. A
validation note was added to docs to indicate change and the fact that
full validation has not really been done. Also all unit tests that had
non-zero salt concentration values have been commented out till we
decide how to handle the change in terms of validating unit tests.

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