source: sasmodels/sasmodels/models/stickyhardsphere.py @ 3e428ec

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 3e428ec was 3e428ec, checked in by Doucet, Mathieu <doucetm@…>, 9 years ago

Cleaning up

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