1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
|
from __future__ import print_function
import numpy as np
from numpy import exp, sin, degrees, radians, pi, sqrt
from sasmodels.weights import Dispersion as BaseDispersion
class Dispersion(BaseDispersion):
r"""
Maier-Saupe dispersion on orientation.
.. math:
w(\theta) = e^{a{\cos^2 \theta}}
This provides a close match to the gaussian distribution for
low angles, but the tails are limited to $\pm 90^\circ$. For $a \ll 1$
the distribution is approximately uniform. The usual polar coordinate
projection applies, with $\theta$ weights scaled by $\cos \theta$
and $\phi$ weights unscaled.
This is equivalent to a cyclic gaussian distribution
$w(\theta) = e^{-sin^2(\theta)/(2\a^2)}.
Note that the incorrect distribution is used for $a=0$. With the
distribution parameter labelled as *width*, the value *width=0* is
is assumed to be completely oriented and no polydispersity distribution
is generated. However a value of $a=0$ should be completely unoriented.
Any small value (e.g., 0.01 or lower) should suffice.
The order parameter $P_2$ is defined as
.. math:
P(a, \beta) = \frac{e^{a \cos^2 \beta}}{
4\pi \int_0^{\pi/2} e^{a \cos^2 \beta} \sin \beta\,d\beta}
P_2 = 4\pi \int_0^{\pi/2} \frac{3 \cos^2 \beta - 1)}{2}
P(a, \beta) \sin \beta\,d\beta
where $a$ is the distribution width $\sigma$ handed to the weights function.
There is a somewhat complicated closed form solution
.. math:
P_2 = \frac{3e^a}{2\sqrt{\pi a} E_a} - \frac{3}{4a} - \frac{1}{2}
where $E_a = \mathrm{erfi}(\sqrt{a})$ is the imaginary error function,
which can be coded in python as::
from numpy import pi, sqrt, exp
from scipy.special import erfi
def P_2(a):
# protect against overflow with a Taylor series at infinity
if a <= 700:
r = exp(a)/sqrt(pi*a)/erfi(sqrt(a))
else:
r = 1/((((6.5525/a + 1.875)/a + 0.75)/a + 0.5)/a + 1)
return 1.5*r - 0.75/a - 0.5
Given an order parameter $S = P_2(a)$, one can also solve for the
equivalent $a$:
from scipy.optimize import fsolve
def P_2_inv(S):
return fsolve(lambda x: P_2(x) - S, 1.0)[0]
References
----------
[1] Hardouin, et al., 1995. SANS study of a semiflexible main chain
liquid crystalline polyether. *Macromolecules* 28, 5427-5433.
"""
type = "maier_saupe"
default = dict(npts=35, width=1, nsigmas=None)
# Note: center is always zero for orientation distributions
def _weights(self, center, sigma, lb, ub):
# use the width parameter as the value for Maier-Saupe "a"
# and find the equivalent width sigma
a = sigma
sigma = 1./sqrt(2.*a)
# Limit width to +/- 90 degrees
width = min(self.nsigmas*sigma, pi/2)
x = np.linspace(-width, width, self.npts)
# Truncate the distribution in case the parameter value is limited
x[(x >= radians(lb)) & (x <= radians(ub))]
# Return orientation in degrees with Maier-Saupe weights
# Note: weights are normalized to sum to 1, so we can scale
# by an arbitrary scale factor c = exp(m) to get:
# w = exp(m*cos(x)**2)/c = exp(-m*sin(x)**2)
return degrees(x), exp(-a*sin(x)**2)
def P_2(a):
from numpy import pi, sqrt, exp
from scipy.special import erfi
# Double precision e^x overflows so do a Taylor expansion around infinity
# for large values of a. With five terms it is accurate to 12 digits
# at a = 700, and 7 digits at a = 75.
if a <= 700:
r = exp(a)/sqrt(pi*a)/erfi(sqrt(a))
else:
r = 1/((((6.5525/a + 1.875)/a + 0.75)/a + 0.5)/a + 1)
return 1.5*r - 0.75/a - 0.5
def P_2_inv(S):
from scipy.optimize import fsolve
return fsolve(lambda x: P_2(x) - S, 1.0)[0]
def P_2_numerical(a):
from scipy.integrate import romberg
from numpy import cos, sin, pi, exp
def num(beta):
return (1.5 * cos(beta)**2 - 0.5) * exp(a * cos(beta)**2) * sin(beta)
def denom(beta):
return exp(a * cos(beta)**2) * sin(beta)
return romberg(num, 0, pi/2) / romberg(denom, 0, pi/2)
if __name__ == "__main__":
import sys
a = float(sys.argv[1])
sigma = 1/(2*radians(a)**2)
#print("P_2", P_2(a), "difference from integral", P_2(a) - P_2_numerical(a))
print("a=%g, sigma=%g, P_2=%g, P_2_inv(P_2(a))-a=%g"
% (a, sigma, P_2(a), P_2_inv(P_2(a))-a))
|
