source: sasmodels/test/rpa.py @ e66c9f9

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

rpa model

  • Property mode set to 100644
File size: 5.7 KB
Line 
1from collections import namedtuple
2import numpy as np
3from numpy import sqrt, exp, expm1
4
5AVOGADRO = 6.022e23
6
7Polymer = namedtuple("Polymer", "n phi v a b".split())
8def E(q, poly):
9    qrsq = (q*Rg(poly))**2
10    retval = exp(-qrsq)
11    return retval
12
13def F(q, poly):
14    qrsq = (q*Rg(poly))**2
15    retval = -expm1(-qrsq)/qrsq
16    return retval
17
18def P_ii(q, poly):
19    qrsq = (q*Rg(poly))**2
20    retval = 2 * (expm1(-qrsq) + qrsq)/qrsq**2
21    return retval
22
23def P_ij(q, poly_list):
24    i, j = poly_list[0], poly_list[-1]
25    retval = F(q, i) * np.prod([E(q,p) for p in poly_list[1:-1]]) * F(q, j)
26    return retval
27
28def Rg(poly):
29    return sqrt(poly.n/6.)*poly.a
30
31def S0_ii(q, poly):
32    retval = poly.n*poly.phi*poly.v*P_ii(q, poly)
33    return retval
34
35def S0_ij(q, poly_list):
36    i,j = poly_list[0], poly_list[-1]
37    retval = sqrt(i.n*i.phi*i.v*j.n*j.phi*j.v) * P_ij(q, poly_list)
38    return retval
39
40def drho(poly, base):
41    return (poly.b/poly.v -  base.b/base.v)*1e-13*sqrt(AVOGADRO)
42
43def binary_blend(q, C, D, Kcd):
44    """
45    de Gennes, 1979
46    """
47    S0cc = S0_ii(q, C)
48    S0dd = S0_ii(q, D)
49    vcc = 1/S0dd - 2*Kcd #/v0
50    Scc = S0cc/(1 + vcc*S0cc)
51    rhocd = drho(C,D)
52    Iq = rhocd**2 * Scc
53    return Iq
54
55def ternary_blend(q, B, C, D, Kbc, Kbd, Kcd):
56    S0bb = S0_ii(q, B)
57    S0cc = S0_ii(q, C)
58    S0dd = S0_ii(q, D)
59    vbb = 1/S0dd - 2*Kbd
60    vcc = 1/S0dd - 2*Kcd
61    vbc = 1/S0dd + Kbd - Kbc - Kcd
62    rhobd = drho(B,D)
63    rhocd = drho(C,D)
64    det = (1+vbb*S0bb)*(1+vcc*S0cc) - vbc**2*S0bb*S0cc
65    Sbb = S0bb*(1+vcc*S0cc)/det
66    Scc = S0cc*(1+vbb*S0bb)/det
67    Sbc = -S0bb*vbc*S0cc/det
68    Iq = rhobd**2*Sbb + rhocd**2*Scc + 2*rhobd*rhocd*Sbc
69    return Iq
70
71def diblock_copolymer(q, C, D, Kcd):
72    """
73    Leibler, 1980
74    """
75    S0cc = S0_ii(q, C)
76    S0dd = S0_ii(q, D)
77    S0cd = S0_ij(q, [C, D])
78    Scc = (S0cc*S0dd - S0cd**2)/((S0cc+S0dd + 2*S0cd)-2*Kcd*(S0cc+S0dd-2*S0cd))
79    rhocd = drho(C,D)
80    Iq = rhocd**2 * Scc
81    return Iq
82
83def matrix_form(q, case_num, polys,
84                Kab=0., Kac=0., Kad=0., Kbc=0., Kbd=0., Kcd=0.):
85    if case_num < 2:
86        C, D = polys
87        A = B = D
88    elif case_num < 5:
89        B, C, D = polys
90        A = D
91    else:
92        A, B, C, D = polys
93
94    rho = np.matrix([[drho(p, D)] for p in (A,B,C)])
95    S0aa = S0_ii(q, A)
96    S0bb = S0_ii(q, B)
97    S0cc = S0_ii(q, C)
98    S0ab = S0_ij(q, [A, B])
99    S0ac = S0_ij(q, [A, B, C])
100    S0bc = S0_ij(q, [B, C])
101    if case_num == 4: # No a-c interaction in triblock copolymer
102        S0ac *= 0.0
103    elif case_num == 9: # No a-c or a-d interaction in tetrablock copolymer
104        S0ac *= 0.0
105
106    S0 = np.array([[S0aa, S0ab, S0ac], [S0ab, S0bb, S0bc], [S0ac, S0bc, S0cc]])
107    S0dd = S0_ii(q, D)
108    vaa = 1./S0dd - 2*Kad
109    vbb = 1./S0dd - 2*Kbd
110    vcc = 1./S0dd - 2*Kcd
111    vab = 1./S0dd + Kab - Kad - Kbd
112    vac = 1./S0dd + Kac - Kad - Kcd
113    vbc = 1./S0dd + Kbc - Kbd - Kcd
114    v = np.array([[vaa, vab, vac], [vab, vbb, vbc], [vac, vbc, vcc]])
115
116    Iq = np.empty_like(q)
117    for k, qk in enumerate(q):
118        S0_k = S0[:,:,k].M
119        v_k = v[:,:,k].M
120        S_k = np.linalg.inv(1 + S0_k*v_k)*S0_k
121        Iq[k] = rho.T * S_k * rho
122
123def build_pars(case_num, polys, **interactions):
124    def set_poly(x, poly):
125        pars["N"+x] = poly.n
126        pars["Phi"+x] = poly.phi
127        pars["v"+x] = poly.v
128        pars["b"+x] = poly.a
129        pars["L"+x] = poly.b
130    pars = interactions.copy()
131    pars["case_num"] = case_num
132    polys = list(reversed(polys))
133    if len(polys) >= 4: set_poly("a",polys[3])
134    if len(polys) >= 3: set_poly("b",polys[2])
135    if len(polys) >= 2: set_poly("c",polys[1])
136    if len(polys) >= 1: set_poly("d",polys[0])
137    return pars
138
139def sasmodels_rpa(q, pars):
140    from sasmodels.models import rpa
141    from sasmodels.core import load_model
142    from sasmodels.direct_model import DirectModel
143    from sasmodels.data import empty_data1D
144    data = empty_data1D(q, resolution=0.0)
145    model = load_model(rpa, dtype="double", platform="dll")
146    #model = load_model(rpa, dtype="single", platform="ocl")
147    M = DirectModel(data, model)
148    return M(**pars)
149
150def sasview_rpa(q, pars):
151    from sasmodels.models import rpa
152    from sasmodels.compare import eval_sasview
153    from sasmodels.data import empty_data1D
154    data = empty_data1D(q, resolution=0.0)
155    M = eval_sasview(rpa, data)
156    return M(**pars)
157
158def demo():
159    import sys
160    case_num = 0 if len(sys.argv) < 2 else int(sys.argv[1])
161
162    B = Polymer(n=525,phi=0.57,v=97.5,b=-4.99,a=8)
163    C = Polymer(n=525,phi=0.57,v=97.5,b=-4.99,a=8)
164    D = Polymer(n=1105,phi=0.43,v=81.9,b=53.1,a=2)
165    q = np.logspace(-4,-1,400)
166    #q = np.array([0.1])
167    Kab=Kac=Kad=0.0
168    Kcd = 0.0106*0.0035 - 1.84e-5
169    Kbd = Kcd + 2e-4
170    Kbc = (Kcd + 1e-4)*0.5
171    K = dict(Kab=Kab,Kac=Kac,Kad=Kad,Kbc=Kbc,Kbd=Kbd,Kcd=Kcd)
172    if case_num == 0:
173        Iq = binary_blend(q, C, D, Kcd)
174    elif case_num == 1:
175        Iq = diblock_copolymer(q, C, D, Kcd)
176    elif case_num == 2:
177        Iq = ternary_blend(q, B, C, D, Kbc, Kbd, Kcd)
178    else:
179        raise ValueError("Case %d not implmented"%case_num)
180
181    pars = build_pars(case_num, [B, C, D], **K)
182    print "eval sasmodels"
183    Iq_sasmodels = sasmodels_rpa(q, pars)
184    print "eval sasview"
185    Iq_sasview = sasview_rpa(q, pars)
186    print 1./Iq[0], 1./Iq_sasmodels[0], 1./Iq_sasview[0]
187
188    #return
189    import pylab
190    pylab.subplot(121)
191    pylab.loglog(q, Iq, label='direct')
192    pylab.loglog(q, Iq_sasmodels, label='sasmodels')
193    pylab.loglog(q, Iq_sasview, label='sasview')
194    pylab.legend()
195    pylab.subplot(122)
196    pylab.loglog(q, abs(Iq_sasview - Iq_sasmodels)/Iq_sasmodels, label='sasview-sasmodels')
197    pylab.loglog(q, abs(Iq_sasmodels - Iq)/Iq, label='sasmodels-direct')
198    pylab.legend()
199    pylab.show()
200
201if __name__ == "__main__":
202    demo()
Note: See TracBrowser for help on using the repository browser.