source: sasmodels/sasmodels/models/hayter_msa.py @ 9947865

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 9947865 was 9947865, checked in by smk78, 6 months ago

Changed low limit for charge in hayter_msa from 0 to 1e-6

  • Property mode set to 100644
File size: 7.2 KB
Line 
1# Note: model title and parameter table are inserted automatically
2r"""
3Calculates the interparticle structure factor for a system of charged,
4spheroidal, objects in a dielectric medium [1,2]. When combined with an
5appropriate form factor $P(q)$, this allows for inclusion of the
6interparticle interference effects due to screened Coulombic
7repulsion between the charged particles.
8
9.. note::
10
11   This routine only works for charged particles! If the charge is set
12   to zero the routine may self-destruct! For uncharged particles use
13   the :ref:`hardsphere` $S(q)$ instead. The upper limit for the charge
14   is limited to 200e to avoid numerical instabilities.
15   
16.. note::
17
18   Earlier versions of SasView did not incorporate the so-called
19   $\beta(q)$ ("beta") correction [3] for polydispersity and non-sphericity.
20   This is only available in SasView versions 4.2.2 and higher.
21
22The salt concentration is used to compute the ionic strength of the solution
23which in turn is used to compute the Debye screening length. There is no
24provision for entering the ionic strength directly. **At present the
25counterions are assumed to be monovalent**, though it should be possible
26to simulate the effect of multivalent counterions by increasing the salt
27concentration.
28
29Over the range 0 - 100 C the dielectric constant $\kappa$ of water may be
30approximated with a maximum deviation of 0.01 units by the empirical
31formula [4]
32
33.. math::
34
35    \kappa = 87.740 - 0.40008 T + 9.398x10^{-4} T^2 - 1.410x10^{-6} T^3
36   
37where $T$ is the temperature in celsius.
38
39In SasView the effective radius may be calculated from the parameters
40used in the form factor $P(q)$ that this $S(q)$ is combined with.
41
42The computation uses a Taylor series expansion at very small rescaled $qR$, to
43avoid some serious rounding error issues, this may result in a minor artefact
44in the transition region under some circumstances.
45
46For 2D data, the scattering intensity is calculated in the same way as 1D,
47where the $q$ vector is defined as
48
49.. math::
50
51    q = \sqrt{q_x^2 + q_y^2}
52
53
54References
55----------
56
57.. [#] J B Hayter and J Penfold, *Molecular Physics*, 42 (1981) 109-118
58
59.. [#] J P Hansen and J B Hayter, *Molecular Physics*, 46 (1982) 651-656
60
61.. [#] M Kotlarchyk and S-H Chen, *J. Chem. Phys.*, 79 (1983) 2461-2469
62
63.. [#] C G Malmberg and A A Maryott, *J. Res. Nat. Bureau Standards*, 56 (1956) 2641
64
65Source
66------
67
68`hayter_msa.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/hayter_msa.py>`_
69
70`hayter_msa.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/hayter_msa.c>`_
71
72Authorship and Verification
73----------------------------
74
75* **Author:**
76* **Last Modified by:**
77* **Last Reviewed by:** Steve King **Date:** March 28, 2019
78* **Source added by :** Steve King **Date:** March 25, 2019
79"""
80
81import numpy as np
82from numpy import inf
83
84category = "structure-factor"
85structure_factor = True
86single = False  # double precision only!
87
88#  dp[0] = 2.0*radius_effective();
89#  dp[1] = fabs(charge());
90#  dp[2] = volfraction();
91#  dp[3] = temperature();
92#  dp[4] = concentration_salt();
93#  dp[5] = dielectconst();
94
95
96
97
98name = "hayter_msa"
99title = "Hayter-Penfold Rescaled Mean Spherical Approximation (RMSA) structure factor for charged spheres"
100description = """\
101    [Hayter-Penfold RMSA charged sphere interparticle S(Q) structure factor]
102        Interparticle structure factor S(Q) for charged hard spheres.
103    This routine only works for charged particles! For uncharged particles
104    use the hardsphere S(q) instead. The "beta(q)" correction is available
105    in versions 4.2.2 and higher.
106"""
107
108
109# pylint: disable=bad-whitespace, line-too-long
110#             [ "name", "units", default, [lower, upper], "type", "description" ],
111#
112# NOTE: SMK, 28Mar19 The upper limit for charge is set to 200 to avoid instabilities noted by PK in
113#       Ticket #1152. Also see the thread in Ticket 859. The docs above also note that charge=0 will
114#       cause problems, yet the default parameters allowed it! After discussions with PK I have
115#       changed it to (an arbitarily) small but non-zero value.  But I haven't changed the low limit
116#       in function random() below.
117#
118parameters = [
119    ["radius_effective", "Ang", 20.75,   [0, inf],    "volume", "effective radius of charged sphere"],
120    ["volfraction",   "None",     0.0192, [0, 0.74],   "", "volume fraction of spheres"],
121    ["charge",        "e",   19.0,    [0.000001, 200],    "", "charge on sphere (in electrons)"],
122    ["temperature",   "K",  318.16,   [0, 450],    "", "temperature, in Kelvin, for Debye length calculation"],
123    ["concentration_salt",      "M",    0.0,    [0, inf], "", "conc of salt, moles/litre, 1:1 electolyte, for Debye length"],
124    ["dielectconst",  "None",    71.08,   [-inf, inf], "", "dielectric constant (relative permittivity) of solvent, kappa, default water, for Debye length"]
125    ]
126# pylint: enable=bad-whitespace, line-too-long
127
128source = ["hayter_msa.c"]
129# No volume normalization despite having a volume parameter
130# This should perhaps be volume normalized?
131form_volume = """
132    return 1.0;
133    """
134
135def random():
136    """Return a random parameter set for the model."""
137    # TODO: too many failures for random hayter_msa parameters
138    pars = dict(
139        scale=1, background=0,
140        radius_effective=10**np.random.uniform(1, 4.7),
141        volfraction=10**np.random.uniform(-2, 0),  # high volume fraction
142        charge=min(int(10**np.random.uniform(0, 1.3)+0.5), 200),
143        temperature=10**np.random.uniform(0, np.log10(450)), # max T = 450
144        #concentration_salt=10**np.random.uniform(-3, 1),
145        dialectconst=10**np.random.uniform(0, 6),
146        #charge=10,
147        #temperature=318.16,
148        concentration_salt=0.0,
149        #dielectconst=71.08,
150    )
151    return pars
152
153# default parameter set,  use  compare.sh -midQ -linear
154# note the calculation varies in different limiting cases so a wide range of
155# parameters will be required for a thorough test!
156# odd that the default st has concentration_salt zero
157demo = dict(radius_effective=20.75,
158            charge=19.0,
159            volfraction=0.0192,
160            temperature=318.16,
161            concentration_salt=0.05,
162            dielectconst=71.08,
163            radius_effective_pd=0.1,
164            radius_effective_pd_n=40)
165#
166# attempt to use same values as old sasview unit test at Q=.001 was 0.0712928,
167# then add lots new ones assuming values from new model are OK, need some
168# low Q values to test the small Q Taylor expansion
169tests = [
170    [{'scale': 1.0,
171      'background': 0.0,
172      'radius_effective': 20.75,
173      'charge': 19.0,
174      'volfraction': 0.0192,
175      'temperature': 298.0,
176      'concentration_salt': 0,
177      'dielectconst': 78.0,
178      'radius_effective_pd': 0},
179     [0.00001, 0.0010, 0.01, 0.075], [0.0711646, 0.0712928, 0.0847006, 1.07150]],
180    [{'scale': 1.0,
181      'background': 0.0,
182      'radius_effective': 20.75,
183      'charge': 19.0,
184      'volfraction': 0.0192,
185      'temperature': 298.0,
186      'concentration_salt': 0.05,
187      'dielectconst': 78.0,
188      'radius_effective_pd': 0.1,
189      'radius_effective_pd_n': 40},
190     [0.00001, 0.0010, 0.01, 0.075], [0.450272, 0.450420, 0.465116, 1.039625]]
191    ]
192# ADDED by:  RKH  ON: 16Mar2016 converted from sasview, new Taylor expansion at smallest rescaled Q
Note: See TracBrowser for help on using the repository browser.