source: sasmodels/sasmodels/models/squarewell.py @ b6422c7

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since b6422c7 was b297ba9, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

lint

  • Property mode set to 100644
File size: 4.7 KB
Line 
1# Note: model title and parameter table are inserted automatically
2r"""
3This calculates the interparticle structure factor for a square well fluid
4spherical particles. The mean spherical approximation (MSA) closure was
5used for this calculation, and is not the most appropriate closure for
6an attractive interparticle potential. This solution has been compared
7to Monte Carlo simulations for a square well fluid, showing this calculation
8to be limited in applicability to well depths $\epsilon < 1.5$ kT and
9volume fractions $\phi < 0.08$.
10
11Positive well depths correspond to an attractive potential well. Negative
12well depths correspond to a potential "shoulder", which may or may not be
13physically reasonable. The stickyhardsphere model may be a better choice in
14some circumstances. Computed values may behave badly at extremely small $qR$.
15
16The well width $(\lambda)$ is defined as multiples of the particle diameter
17$(2 R)$.
18
19The interaction potential is:
20
21  .. image:: img/squarewell.png
22
23.. math::
24
25    U(r) = \begin{cases}
26    \infty & r < 2R \\
27    -\epsilon & 2R \leq r < 2R\lambda \\
28    0 & r \geq 2R\lambda
29    \end{cases}
30
31where $r$ is the distance from the center of the sphere of a radius $R$.
32
33In sasview the effective radius may be calculated from the parameters
34used in the form factor $P(q)$ that this $S(q)$ is combined with.
35
36For 2D data: The 2D scattering intensity is calculated in the same way as 1D,
37where the $q$ vector is defined as
38
39.. math::
40
41    q = \sqrt{q_x^2 + q_y^2}
42
43References
44----------
45
46R V Sharma, K C Sharma, *Physica*, 89A (1977) 213.
47"""
48
49import numpy as np
50from numpy import inf
51
52name = "squarewell"
53title = "Square well structure factor, with MSA closure"
54description = """\
55    [Square well structure factor, with MSA closure]
56        Interparticle structure factor S(Q)for a hard sphere fluid with
57        a narrow attractive well. Fits are prone to deliver non-physical
58        parameters, use with care and read the references in the full manual.
59        In sasview the effective radius will be calculated from the
60        parameters used in P(Q).
61"""
62category = "structure-factor"
63structure_factor = True
64single = False
65
66#single = False
67#             ["name", "units", default, [lower, upper], "type","description"],
68parameters = [
69    #   [ "name", "units", default, [lower, upper], "type",
70    #     "description" ],
71    ["radius_effective", "Ang", 50.0, [0, inf], "volume",
72     "effective radius of hard sphere"],
73    ["volfraction", "", 0.04, [0, 0.08], "",
74     "volume fraction of spheres"],
75    ["welldepth", "kT", 1.5, [0.0, 1.5], "",
76     "depth of well, epsilon"],
77    ["wellwidth", "diameters", 1.2, [1.0, inf], "",
78     "width of well in diameters (=2R) units, must be > 1"],
79    ]
80
81# No volume normalization despite having a volume parameter
82# This should perhaps be volume normalized?
83form_volume = """
84    return 1.0;
85    """
86
87Iq = """
88// single precision is very poor at extreme small Q, would need a Taylor series
89        double req,phis,edibkb,lambda,struc;
90        double sigma,eta,eta2,eta3,eta4,etam1,etam14,alpha,beta,gamm;
91        double x,sk,sk2,sk3,sk4,t1,t2,t3,t4,ck;
92        double S,C,SL,CL;
93        x= q;
94
95        req = radius_effective;
96        phis = volfraction;
97        edibkb = welldepth;
98        lambda = wellwidth;
99
100        sigma = req*2.;
101        eta = phis;
102        eta2 = eta*eta;
103        eta3 = eta*eta2;
104        eta4 = eta*eta3;
105        etam1 = 1. - eta;
106        etam14 = etam1*etam1*etam1*etam1;
107        // temp borrow sk for an intermediate calc
108        sk = 1.0 +2.0*eta;
109        sk *= sk;
110        alpha = (  sk + eta3*( eta-4.0 )  )/etam14;
111        beta = -(eta/3.0) * ( 18. + 20.*eta - 12.*eta2 + eta4 )/etam14;
112        gamm = 0.5*eta*( sk + eta3*(eta-4.) )/etam14;
113
114        //  calculate the structure factor
115
116        sk = x*sigma;
117        sk2 = sk*sk;
118        sk3 = sk*sk2;
119        sk4 = sk3*sk;
120        SINCOS(sk,S,C);
121        SINCOS(lambda*sk,SL,CL);
122        t1 = alpha * sk3 * ( S - sk * C );
123        t2 = beta * sk2 * 2.0*( sk*S - (0.5*sk2 - 1.)*C - 1.0 );
124        t3 = gamm*( ( 4.0*sk3 - 24.*sk ) * S - ( sk4 - 12.0*sk2 + 24.0 )*C + 24.0 );
125        t4 = -edibkb*sk3*(SL +sk*(C - lambda*CL) - S );
126        ck =  -24.0*eta*( t1 + t2 + t3 + t4 )/sk3/sk3;
127        struc  = 1./(1.-ck);
128
129        return(struc);
130"""
131
132def random():
133    """Return a random parameter set for the model."""
134    pars = dict(
135        scale=1, background=0,
136        radius_effective=10**np.random.uniform(1, 4.7),
137        volfraction=np.random.uniform(0.00001, 0.08),
138        welldepth=np.random.uniform(0, 1.5),
139        wellwidth=np.random.uniform(1, 1.2),
140    )
141    return pars
142
143demo = dict(radius_effective=50, volfraction=0.04, welldepth=1.5,
144            wellwidth=1.2, radius_effective_pd=0, radius_effective_pd_n=0)
145#
146tests = [
147    [{'scale': 1.0, 'background': 0.0, 'radius_effective': 50.0,
148      'volfraction': 0.04, 'welldepth': 1.5, 'wellwidth': 1.2,
149      'radius_effective_pd': 0}, [0.001], [0.97665742]],
150    ]
151# ADDED by: converting from sasview RKH  ON: 16Mar2016
Note: See TracBrowser for help on using the repository browser.