source: sasmodels/sasmodels/models/stickyhardsphere.py @ db1d9d5

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since db1d9d5 was db1d9d5, checked in by Paul Kienzle <pkienzle@…>, 8 months ago

merge with master

  • Property mode set to 100644
File size: 7.0 KB
Line 
1# Note: model title and parameter table are inserted automatically
2r"""
3Calculates the interparticle structure factor for a hard sphere fluid
4with a narrow, attractive, potential well. Unlike the :ref:`squarewell`
5model, here a perturbative solution of the Percus-Yevick closure
6relationship is used. The strength of the attractive well is described
7in terms of "stickiness" as defined below.
8
9The perturbation parameter (perturb), $\tau$, should be fixed between 0.01
10and 0.1 and the "stickiness", $\epsilon$, allowed to vary to adjust the
11interaction strength. The "stickiness" is defined in the equation below and is
12a function of both the perturbation parameter and the interaction strength.
13$\epsilon$ and $\tau$ are defined in terms of the hard sphere diameter $(\sigma = 2 R)$,
14the width of the square well, $\Delta$ (having the same units as $R$\ ),
15and the depth of the well, $U_o$, in units of $kT$. From the definition, it
16is clear that smaller $\epsilon$ means a stronger attraction.
17
18.. math::
19
20    \epsilon     &= \frac{1}{12\tau} \exp(u_o / kT) \\
21    \tau &= \Delta / (\sigma + \Delta)
22
23where the interaction potential is
24
25.. math::
26
27    U(r) = \begin{cases}
28        \infty & r < \sigma \\
29        -U_o   & \sigma \leq r \leq \sigma + \Delta \\
30        0      & r > \sigma + \Delta
31        \end{cases}
32
33The Percus-Yevick (PY) closure is used for this calculation, and is an
34adequate closure for an attractive interparticle potential. The solution
35has been compared to Monte Carlo simulations for a square well fluid, with
36good agreement.
37
38The true particle volume fraction, $\phi$, is not equal to $h$ which appears
39in most of reference [1]. The two are related in equation (24). Reference
40[1] also describes the relationship between this perturbative solution and
41the original sticky hard sphere (or "adhesive sphere") model of Baxter [2].
42
43.. note::
44
45   The calculation can go haywire for certain combinations of the input
46   parameters, producing unphysical solutions. In this case errors are
47   reported to the command window and $S(q)$ is set to -1 (so it will
48   disappear on a log-log plot!).
49
50   Use tight bounds to keep the parameters to values that you know are
51   physical (test them), and keep nudging them until the optimization
52   does not hit the constraints.
53
54.. note::
55
56   Earlier versions of SasView did not incorporate the so-called
57   $\beta(q)$ ("beta") correction [3] for polydispersity and non-sphericity.
58   This is only available in SasView versions 4.2.2 and higher.
59
60In SasView the effective radius may be calculated from the parameters
61used in the form factor $P(q)$ that this $S(q)$ is combined with.
62
63For 2D data the scattering intensity is calculated in the same way
64as 1D, where the $q$ vector is defined as
65
66.. math::
67
68    q = \sqrt{q_x^2 + q_y^2}
69
70
71References
72----------
73
74.. [#] S V G Menon, C Manohar, and K S Rao, *J. Chem. Phys.*, 95(12) (1991) 9186-9190
75
76.. [#] R J Baxter, *J. Chem. Phys.*, 49 (1968), 2770-2774
77
78.. [#] M Kotlarchyk and S-H Chen, *J. Chem. Phys.*, 79 (1983) 2461-2469
79
80Source
81------
82
83`stickyhardsphere.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/stickyhardsphere.py>`_
84
85Authorship and Verification
86----------------------------
87
88* **Author:**
89* **Last Modified by:**
90* **Last Reviewed by:** Steve King **Date:** March 27, 2019
91* **Source added by :** Steve King **Date:** March 25, 2019
92"""
93
94# TODO: refactor so that we pull in the old sansmodels.c_extensions
95
96import numpy as np
97from numpy import inf
98
99name = "stickyhardsphere"
100title = "'Sticky' hard sphere structure factor with Percus-Yevick closure"
101description = """\
102    [Sticky hard sphere structure factor, with Percus-Yevick closure]
103        Interparticle structure factor S(Q) for a hard sphere fluid
104    with a narrow attractive well. Fits are prone to deliver non-
105    physical parameters; use with care and read the references in
106    the model documentation.The "beta(q)" correction is available
107    in versions 4.2.2 and higher.
108"""
109category = "structure-factor"
110structure_factor = True
111
112single = False
113#             ["name", "units", default, [lower, upper], "type","description"],
114parameters = [
115    #   [ "name", "units", default, [lower, upper], "type",
116    #     "description" ],
117    ["radius_effective", "Ang", 50.0, [0, inf], "volume",
118     "effective radius of hard sphere"],
119    ["volfraction", "", 0.2, [0, 0.74], "",
120     "volume fraction of hard spheres"],
121    ["perturb", "", 0.05, [0.01, 0.1], "",
122     "perturbation parameter, tau"],
123    ["stickiness", "", 0.20, [-inf, inf], "",
124     "stickiness, epsilon"],
125    ]
126
127def random():
128    """Return a random parameter set for the model."""
129    pars = dict(
130        scale=1, background=0,
131        radius_effective=10**np.random.uniform(1, 4.7),
132        volfraction=np.random.uniform(0.00001, 0.74),
133        perturb=10**np.random.uniform(-2, -1),
134        stickiness=np.random.uniform(0, 1),
135    )
136    return pars
137
138# No volume normalization despite having a volume parameter
139# This should perhaps be volume normalized?
140form_volume = """
141    return 1.0;
142    """
143
144Iq = """
145    double onemineps,eta;
146    double sig,aa,etam1,etam1sq,qa,qb,qc,radic;
147    double lam,lam2,test,mu,alpha,beta;
148    double kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq;
149
150    onemineps = 1.0-perturb;
151    eta = volfraction/onemineps/onemineps/onemineps;
152
153    sig = 2.0 * radius_effective;
154    aa = sig/onemineps;
155    etam1 = 1.0 - eta;
156    etam1sq=etam1*etam1;
157    //C
158    //C  SOLVE QUADRATIC FOR LAMBDA
159    //C
160    qa = eta/6.0;
161    qb = stickiness + eta/etam1;
162    qc = (1.0 + eta/2.0)/etam1sq;
163    radic = qb*qb - 2.0*qa*qc;
164    if(radic<0) {
165        //if(x>0.01 && x<0.015)
166        //    Print "Lambda unphysical - both roots imaginary"
167        //endif
168        return(-1.0);
169    }
170    //C   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL
171    radic = sqrt(radic);
172    lam = (qb-radic)/qa;
173    lam2 = (qb+radic)/qa;
174    if(lam2<lam) {
175        lam = lam2;
176    }
177    test = 1.0 + 2.0*eta;
178    mu = lam*eta*etam1;
179    if(mu>test) {
180        //if(x>0.01 && x<0.015)
181        // Print "Lambda unphysical mu>test"
182        //endif
183        return(-1.0);
184    }
185    alpha = (1.0 + 2.0*eta - mu)/etam1sq;
186    beta = (mu - 3.0*eta)/(2.0*etam1sq);
187    //C
188    //C   CALCULATE THE STRUCTURE FACTOR
189    //C
190    kk = q*aa;
191    k2 = kk*kk;
192    k3 = kk*k2;
193    SINCOS(kk,ds,dc);
194    //ds = sin(kk);
195    //dc = cos(kk);
196    aq1 = ((ds - kk*dc)*alpha)/k3;
197    aq2 = (beta*(1.0-dc))/k2;
198    aq3 = (lam*ds)/(12.0*kk);
199    aq = 1.0 + 12.0*eta*(aq1+aq2-aq3);
200    //
201    bq1 = alpha*(0.5/kk - ds/k2 + (1.0 - dc)/k3);
202    bq2 = beta*(1.0/kk - ds/k2);
203    bq3 = (lam/12.0)*((1.0 - dc)/kk);
204    bq = 12.0*eta*(bq1+bq2-bq3);
205    //
206    sq = 1.0/(aq*aq +bq*bq);
207
208    return(sq);
209"""
210
211demo = dict(radius_effective=200, volfraction=0.2, perturb=0.05,
212            stickiness=0.2, radius_effective_pd=0.1, radius_effective_pd_n=40)
213#
214tests = [
215    [{'scale': 1.0, 'background': 0.0, 'radius_effective': 50.0,
216      'perturb': 0.05, 'stickiness': 0.2, 'volfraction': 0.1,
217      'radius_effective_pd': 0},
218     [0.001, 0.003], [1.09718, 1.087830]],
219    ]
Note: See TracBrowser for help on using the repository browser.