source: sasmodels/sasmodels/models/stickyhardsphere.py @ 0507e09

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 0507e09 was 0507e09, checked in by smk78, 9 months ago

Added link to source code to each model. Closes #883

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