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

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since eb69cce was eb69cce, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

make model docs more consistent; build pdf docs

  • Property mode set to 100644
File size: 6.0 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    %\begin{align*} % isn't working with pdflatex
21    \begin{array}{rl}
22    \tau     &= \frac{1}{12\epsilon} \exp(u_o / kT) \\
23    \epsilon &= \Delta / (\sigma + \Delta) \\
24    \end{array}
25
26where the interaction potential is
27
28.. math::
29
30    U(r) = \begin{cases}
31        \infty & r < \sigma \\
32        -U_o   & \sigma \leq r \leq \sigma + \Delta \\
33        0      & r > \sigma + \Delta
34        \end{cases}
35
36The Percus-Yevick (PY) closure was used for this calculation, and is an
37adequate closure for an attractive interparticle potential. This solution
38has been compared to Monte Carlo simulations for a square well fluid, with
39good agreement.
40
41The true particle volume fraction, $\phi$, is not equal to $h$, which appears
42in most of the reference. The two are related in equation (24) of the
43reference. The reference also describes the relationship between this
44perturbation solution and the original sticky hard sphere (or adhesive
45sphere) model by Baxter.
46
47**NB**: The calculation can go haywire for certain combinations of the input
48parameters, producing unphysical solutions - in this case errors are
49reported to the command window and the $S(q)$ is set to -1 (so it will
50disappear on a log-log plot). Use tight bounds to keep the parameters to
51values that you know are physical (test them) and keep nudging them until
52the optimization does not hit the constraints.
53
54In sasview the effective radius will be calculated from the parameters
55used in the form factor $P(q)$ that this $S(q)$ is combined with.
56
57For 2D data the scattering intensity is calculated in the same way
58as 1D, where the $q$ vector is defined as
59
60.. math::
61
62    q = \sqrt{q_x^2 + q_y^2}
63
64.. figure:: img/stickyhardsphere_1d.jpg
65
66    1D plot using the default values (in linear scale).
67
68References
69----------
70
71S V G Menon, C Manohar, and K S Rao, *J. Chem. Phys.*, 95(12) (1991) 9186-9190
72"""
73
74# TODO: refactor so that we pull in the old sansmodels.c_extensions
75
76from numpy import inf
77
78name = "stickyhardsphere"
79title = "Sticky hard sphere structure factor, with Percus-Yevick closure"
80description = """\
81    [Sticky hard sphere structure factor, with Percus-Yevick closure]
82        Interparticle structure factor S(Q)for a hard sphere fluid with
83        a narrow attractive well. Fits are prone to deliver non-physical
84        parameters, use with care and read the references in the full manual.
85        In sasview the effective radius will be calculated from the
86        parameters used in P(Q).
87"""
88category = "structure-factor"
89
90#             ["name", "units", default, [lower, upper], "type","description"],
91parameters = [
92    #   [ "name", "units", default, [lower, upper], "type",
93    #     "description" ],
94    ["effect_radius", "Ang", 50.0, [0, inf], "volume",
95     "effective radius of hard sphere"],
96    ["volfraction", "", 0.2, [0, 0.74], "",
97     "volume fraction of hard spheres"],
98    ["perturb", "", 0.05, [0.01, 0.1], "",
99     "perturbation parameter, epsilon"],
100    ["stickiness", "", 0.20, [-inf, inf], "",
101     "stickiness, tau"],
102    ]
103
104# No volume normalization despite having a volume parameter
105# This should perhaps be volume normalized?
106form_volume = """
107    return 1.0;
108    """
109
110Iq = """
111    double onemineps,eta;
112    double sig,aa,etam1,etam1sq,qa,qb,qc,radic;
113    double lam,lam2,test,mu,alpha,beta;
114    double kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq;
115
116    onemineps = 1.0-perturb;
117    eta = volfraction/onemineps/onemineps/onemineps;
118
119    sig = 2.0 * effect_radius;
120    aa = sig/onemineps;
121    etam1 = 1.0 - eta;
122    etam1sq=etam1*etam1;
123    //C
124    //C  SOLVE QUADRATIC FOR LAMBDA
125    //C
126    qa = eta/12.0;
127    qb = -1.0*(stickiness + eta/etam1);
128    qc = (1.0 + eta/2.0)/etam1sq;
129    radic = qb*qb - 4.0*qa*qc;
130    if(radic<0) {
131        //if(x>0.01 && x<0.015)
132        //    Print "Lambda unphysical - both roots imaginary"
133        //endif
134        return(-1.0);
135    }
136    //C   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL
137    lam = (-1.0*qb-sqrt(radic))/(2.0*qa);
138    lam2 = (-1.0*qb+sqrt(radic))/(2.0*qa);
139    if(lam2<lam) {
140        lam = lam2;
141    }
142    test = 1.0 + 2.0*eta;
143    mu = lam*eta*etam1;
144    if(mu>test) {
145        //if(x>0.01 && x<0.015)
146        // Print "Lambda unphysical mu>test"
147        //endif
148        return(-1.0);
149    }
150    alpha = (1.0 + 2.0*eta - mu)/etam1sq;
151    beta = (mu - 3.0*eta)/(2.0*etam1sq);
152    //C
153    //C   CALCULATE THE STRUCTURE FACTOR
154    //C
155    kk = q*aa;
156    k2 = kk*kk;
157    k3 = kk*k2;
158    SINCOS(kk,ds,dc);
159    //ds = sin(kk);
160    //dc = cos(kk);
161    aq1 = ((ds - kk*dc)*alpha)/k3;
162    aq2 = (beta*(1.0-dc))/k2;
163    aq3 = (lam*ds)/(12.0*kk);
164    aq = 1.0 + 12.0*eta*(aq1+aq2-aq3);
165    //
166    bq1 = alpha*(0.5/kk - ds/k2 + (1.0 - dc)/k3);
167    bq2 = beta*(1.0/kk - ds/k2);
168    bq3 = (lam/12.0)*((1.0 - dc)/kk);
169    bq = 12.0*eta*(bq1+bq2-bq3);
170    //
171    sq = 1.0/(aq*aq +bq*bq);
172
173    return(sq);
174"""
175
176Iqxy = """
177    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
178    """
179
180# ER defaults to 0.0
181# VR defaults to 1.0
182
183oldname = 'StickyHSStructure'
184oldpars = dict()
185demo = dict(effect_radius=200, volfraction=0.2, perturb=0.05,
186            stickiness=0.2, effect_radius_pd=0.1, effect_radius_pd_n=40)
187
Note: See TracBrowser for help on using the repository browser.