source: sasmodels/sasmodels/models/hardsphere.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@…>, 7 months ago

merge with master

  • Property mode set to 100644
File size: 7.4 KB
Line 
1# Note: model title and parameter table are inserted automatically
2r"""
3Calculates the interparticle structure factor for monodisperse
4spherical particles interacting through hard sphere (excluded volume)
5interactions. This $S(q)$ may also be a reasonable approximation for
6other particle shapes that freely rotate (but see the note below),
7and for moderately polydisperse systems.
8
9.. note::
10
11   This routine is intended for uncharged particles! For charged
12   particles try using the :ref:`hayter-msa` $S(q)$ instead.
13
14.. note::
15
16   Earlier versions of SasView did not incorporate the so-called
17   $\beta(q)$ ("beta") correction [1] for polydispersity and non-sphericity.
18   This is only available in SasView versions 4.2.2 and higher.
19
20radius_effective is the effective hard sphere radius.
21volfraction is the volume fraction occupied by the spheres.
22
23In SasView the effective radius may be calculated from the parameters
24used in the form factor $P(q)$ that this $S(q)$ is combined with.
25
26For numerical stability the computation uses a Taylor series expansion
27at very small $qR$, but there may be a very minor glitch at the
28transition point in some circumstances.
29
30This S(q) uses the Percus-Yevick closure relationship [2] where the
31interparticle potential $U(r)$ is
32
33.. math::
34
35    U(r) = \begin{cases}
36    \infty & r < 2R \\
37    0 & r \geq 2R
38    \end{cases}
39
40where $r$ is the distance from the center of a sphere of a radius $R$.
41
42For a 2D plot, the wave transfer is defined as
43
44.. math::
45
46    q = \sqrt{q_x^2 + q_y^2}
47
48
49References
50----------
51
52.. [#] M Kotlarchyk & S-H Chen, *J. Chem. Phys.*, 79 (1983) 2461-2469
53
54.. [#] J K Percus, J Yevick, *J. Phys. Rev.*, 110, (1958) 1
55
56Source
57------
58
59`hardsphere.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/hardsphere.py>`_
60
61Authorship and Verification
62----------------------------
63
64* **Author:**
65* **Last Modified by:**
66* **Last Reviewed by:**
67* **Source added by :** Steve King **Date:** March 25, 2019
68"""
69
70import numpy as np
71from numpy import inf
72
73name = "hardsphere"
74title = "Hard sphere structure factor, with Percus-Yevick closure"
75description = """\
76    [Hard sphere structure factor, with Percus-Yevick closure]
77        Interparticle S(Q) for random, non-interacting spheres.
78    May be a reasonable approximation for other particle shapes
79    that freely rotate, and for moderately polydisperse systems
80    . The "beta(q)" correction is available in versions 4.2.2
81    and higher.
82    radius_effective is the hard sphere radius
83    volfraction is the volume fraction occupied by the spheres.
84"""
85category = "structure-factor"
86structure_factor = True
87single = False # TODO: check
88
89#             ["name", "units", default, [lower, upper], "type","description"],
90parameters = [["radius_effective", "Ang", 50.0, [0, inf], "",
91               "effective radius of hard sphere"],
92              ["volfraction", "", 0.2, [0, 0.74], "",
93               "volume fraction of hard spheres"],
94             ]
95
96Iq = r"""
97      double D,A,B,G,X,X2,X4,S,C,FF,HARDSPH;
98      // these are c compiler instructions, can also put normal code inside the "if else" structure
99      #if FLOAT_SIZE > 4
100      // double precision
101      // orig had 0.2, don't call the variable cutoff as PAK already has one called that!
102      // Must use UPPERCASE name please.
103      // 0.05 better, 0.1 OK
104      #define CUTOFFHS 0.05
105      #else
106      // 0.1 bad, 0.2 OK, 0.3 good, 0.4 better, 0.8 no good
107      #define CUTOFFHS 0.4
108      #endif
109
110      if(fabs(radius_effective) < 1.E-12) {
111               HARDSPH=1.0;
112//printf("HS1 %g: %g\n",q,HARDSPH);
113               return(HARDSPH);
114      }
115      // removing use of pow(xxx,2) and rearranging the calcs
116      // of A, B & G cut ~40% off execution time ( 0.5 to 0.3 msec)
117      X = 1.0/( 1.0 -volfraction);
118      D= X*X;
119      A= (1.+2.*volfraction)*D;
120      A *=A;
121      X=fabs(q*radius_effective*2.0);
122
123      if(X < 5.E-06) {
124                 HARDSPH=1./A;
125//printf("HS2 %g: %g\n",q,HARDSPH);
126                 return(HARDSPH);
127      }
128      X2 =X*X;
129      B = (1.0 +0.5*volfraction)*D;
130      B *= B;
131      B *= -6.*volfraction;
132      G=0.5*volfraction*A;
133
134      if(X < CUTOFFHS) {
135      // RKH Feb 2016, use Taylor series expansion for small X
136      // else no obvious way to rearrange the equations to avoid
137      // needing a very high number of significant figures.
138      // Series expansion found using Mathematica software. Numerical test
139      // in .xls showed terms to X^2 are sufficient
140      // for 5 or 6 significant figures, but I put the X^4 one in anyway
141            //FF = 8*A +6*B + 4*G - (0.8*A +2.0*B/3.0 +0.5*G)*X2 +(A/35. +B/40. +G/50.)*X4;
142            // refactoring the polynomial makes it very slightly faster (0.5 not 0.6 msec)
143            //FF = 8*A +6*B + 4*G + ( -0.8*A -2.0*B/3.0 -0.5*G +(A/35. +B/40. +G/50.)*X2)*X2;
144
145            FF = 8.0*A +6.0*B + 4.0*G + ( -0.8*A -B/1.5 -0.5*G +(A/35. +0.0125*B +0.02*G)*X2)*X2;
146
147            // combining the terms makes things worse at smallest Q in single precision
148            //FF = (8-0.8*X2)*A +(3.0-X2/3.)*2*B + (4+0.5*X2)*G +(A/35. +B/40. +G/50.)*X4;
149            // note that G = -volfraction*A/2, combining this makes no further difference at smallest Q
150            //FF = (8 +2.*volfraction + ( volfraction/4. -0.8 +(volfraction/100. -1./35.)*X2 )*X2 )*A  + (3.0 -X2/3. +X4/40.)*2.*B;
151            HARDSPH= 1./(1. + volfraction*FF );
152//printf("HS3 %g: %g\n",q,HARDSPH);
153            return(HARDSPH);
154      }
155      X4=X2*X2;
156      SINCOS(X,S,C);
157
158// RKH Feb 2016, use version FISH code as is better than original sasview one
159// at small Q in single precision, and more than twice as fast in double.
160      //FF=A*(S-X*C)/X + B*(2.*X*S -(X2-2.)*C -2.)/X2 + G*( (4.*X2*X -24.*X)*S -(X4 -12.*X2 +24.)*C +24. )/X4;
161      // refactoring the polynomial here & above makes it slightly faster
162
163      FF=  (( G*( (4.*X2 -24.)*X*S -(X4 -12.*X2 +24.)*C +24. )/X2 + B*(2.*X*S -(X2-2.)*C -2.) )/X + A*(S-X*C))/X ;
164      HARDSPH= 1./(1. + 24.*volfraction*FF/X2 );
165
166      // changing /X and /X2 to *MX1 and *MX2, no significantg difference?
167      //MX=1.0/X;
168      //MX2=MX*MX;
169      //FF=  (( G*( (4.*X2 -24.)*X*S -(X4 -12.*X2 +24.)*C +24. )*MX2 + B*(2.*X*S -(X2-2.)*C -2.) )*MX + A*(S-X*C)) ;
170      //HARDSPH= 1./(1. + 24.*volfraction*FF*MX2*MX );
171
172// grouping the terms, was about same as sasmodels for single precision issues
173//     FF=A*(S/X-C) + B*(2.*S/X - C +2.0*(C-1.0)/X2) + G*( (4./X -24./X3)*S -(1.0 -12./X2 +24./X4)*C +24./X4 );
174//     HARDSPH= 1./(1. + 24.*volfraction*FF/X2 );
175// remove 1/X2 from final line, take more powers of X inside the brackets, stil bad
176//      FF=A*(S/X3-C/X2) + B*(2.*S/X3 - C/X2 +2.0*(C-1.0)/X4) + G*( (4./X -24./X3)*S -(1.0 -12./X2 +24./X4)*C +24./X4 )/X2;
177//      HARDSPH= 1./(1. + 24.*volfraction*FF );
178//printf("HS4 %g: %g\n",q,HARDSPH);
179      return(HARDSPH);
180   """
181
182def random():
183    """Return a random parameter set for the model."""
184    pars = dict(
185        scale=1, background=0,
186        radius_effective=10**np.random.uniform(1, 4),
187        volfraction=10**np.random.uniform(-2, 0),  # high volume fraction
188    )
189    return pars
190
191# Q=0.001 is in the Taylor series, low Q part, so add Q=0.1,
192# assuming double precision sasview is correct
193tests = [
194    [{'scale': 1.0, 'background' : 0.0, 'radius_effective' : 50.0,
195      'volfraction' : 0.2, 'radius_effective_pd' : 0},
196     [0.001, 0.1], [0.209128, 0.930587]],
197]
198# ADDED by: RKH  ON: 16Mar2016  using equations from FISH as better than
199# orig sasview, see notes above. Added Taylor expansions at small Q.
Note: See TracBrowser for help on using the repository browser.