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

Last change on this file since be0942c 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: 6.8 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
80Authorship and Verification
81----------------------------
82
83* **Author:**
84* **Last Modified by:**
85* **Last Reviewed by:** Steve King **Date:** March 27, 2019
86"""
87
88# TODO: refactor so that we pull in the old sansmodels.c_extensions
89
90import numpy as np
91from numpy import inf
92
93name = "stickyhardsphere"
94title = "'Sticky' hard sphere structure factor with Percus-Yevick closure"
95description = """\
96    [Sticky hard sphere structure factor, with Percus-Yevick closure]
97        Interparticle structure factor S(Q) for a hard sphere fluid
98    with a narrow attractive well. Fits are prone to deliver non-
99    physical parameters; use with care and read the references in
100    the model documentation.The "beta(q)" correction is available
101    in versions 4.2.2 and higher.
102"""
103category = "structure-factor"
104structure_factor = True
105
106single = False
107#             ["name", "units", default, [lower, upper], "type","description"],
108parameters = [
109    #   [ "name", "units", default, [lower, upper], "type",
110    #     "description" ],
111    ["radius_effective", "Ang", 50.0, [0, inf], "volume",
112     "effective radius of hard sphere"],
113    ["volfraction", "", 0.2, [0, 0.74], "",
114     "volume fraction of hard spheres"],
115    ["perturb", "", 0.05, [0.01, 0.1], "",
116     "perturbation parameter, tau"],
117    ["stickiness", "", 0.20, [-inf, inf], "",
118     "stickiness, epsilon"],
119    ]
120
121def random():
122    """Return a random parameter set for the model."""
123    pars = dict(
124        scale=1, background=0,
125        radius_effective=10**np.random.uniform(1, 4.7),
126        volfraction=np.random.uniform(0.00001, 0.74),
127        perturb=10**np.random.uniform(-2, -1),
128        stickiness=np.random.uniform(0, 1),
129    )
130    return pars
131
132# No volume normalization despite having a volume parameter
133# This should perhaps be volume normalized?
134form_volume = """
135    return 1.0;
136    """
137
138Iq = """
139    double onemineps,eta;
140    double sig,aa,etam1,etam1sq,qa,qb,qc,radic;
141    double lam,lam2,test,mu,alpha,beta;
142    double kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq;
143
144    onemineps = 1.0-perturb;
145    eta = volfraction/onemineps/onemineps/onemineps;
146
147    sig = 2.0 * radius_effective;
148    aa = sig/onemineps;
149    etam1 = 1.0 - eta;
150    etam1sq=etam1*etam1;
151    //C
152    //C  SOLVE QUADRATIC FOR LAMBDA
153    //C
154    qa = eta/6.0;
155    qb = stickiness + eta/etam1;
156    qc = (1.0 + eta/2.0)/etam1sq;
157    radic = qb*qb - 2.0*qa*qc;
158    if(radic<0) {
159        //if(x>0.01 && x<0.015)
160        //    Print "Lambda unphysical - both roots imaginary"
161        //endif
162        return(-1.0);
163    }
164    //C   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL
165    radic = sqrt(radic);
166    lam = (qb-radic)/qa;
167    lam2 = (qb+radic)/qa;
168    if(lam2<lam) {
169        lam = lam2;
170    }
171    test = 1.0 + 2.0*eta;
172    mu = lam*eta*etam1;
173    if(mu>test) {
174        //if(x>0.01 && x<0.015)
175        // Print "Lambda unphysical mu>test"
176        //endif
177        return(-1.0);
178    }
179    alpha = (1.0 + 2.0*eta - mu)/etam1sq;
180    beta = (mu - 3.0*eta)/(2.0*etam1sq);
181    //C
182    //C   CALCULATE THE STRUCTURE FACTOR
183    //C
184    kk = q*aa;
185    k2 = kk*kk;
186    k3 = kk*k2;
187    SINCOS(kk,ds,dc);
188    //ds = sin(kk);
189    //dc = cos(kk);
190    aq1 = ((ds - kk*dc)*alpha)/k3;
191    aq2 = (beta*(1.0-dc))/k2;
192    aq3 = (lam*ds)/(12.0*kk);
193    aq = 1.0 + 12.0*eta*(aq1+aq2-aq3);
194    //
195    bq1 = alpha*(0.5/kk - ds/k2 + (1.0 - dc)/k3);
196    bq2 = beta*(1.0/kk - ds/k2);
197    bq3 = (lam/12.0)*((1.0 - dc)/kk);
198    bq = 12.0*eta*(bq1+bq2-bq3);
199    //
200    sq = 1.0/(aq*aq +bq*bq);
201
202    return(sq);
203"""
204
205demo = dict(radius_effective=200, volfraction=0.2, perturb=0.05,
206            stickiness=0.2, radius_effective_pd=0.1, radius_effective_pd_n=40)
207#
208tests = [
209    [{'scale': 1.0, 'background': 0.0, 'radius_effective': 50.0,
210      'perturb': 0.05, 'stickiness': 0.2, 'volfraction': 0.1,
211      'radius_effective_pd': 0},
212     [0.001, 0.003], [1.09718, 1.087830]],
213    ]
Note: See TracBrowser for help on using the repository browser.