11. scipyの基本と応用¶

scipyは科学技術計算を行うライブラリで、numpyで作成したデータに基づいて様々な数値解析を高速かつ容易に行うことができます。なお、numpy.linalgで利用できる関数は、同じものがscipy.linalgでも利用できます。scipy.linalg で呼び出す関数の方が、numpy.linalgより高速になっている場合があります。 scipyをimportした場合の線形代数計算は、scipy.linalgを使う方が良いです。

numpy同様、scipyは巨大なライブラリです。信号解析からデータ圧縮まで非常に多くの機能を有しますが、化学と関係のありそうなものだけピックアップして紹介します。

https://docs.scipy.org/doc/scipy/reference/index.html

今回の学習目標は以下のとおりです。

scipyを使ったパラメータの最適化法を理解する。
scipyで数値積分を行う。
scipyで電子軌道と関係する特殊関数を呼び出してみる。

最適化とパラメータフィッテング¶

https://docs.scipy.org/doc/scipy/reference/optimize.html

fsolve¶

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html

x軸との交点を求める。

\[f(x) = 2x^2 - 20x - 15\]

example-01¶

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

x = np.linspace(-10, 10, 1000)
def func(x):
    f = 3 * x**2 - 4 * x - 80
    return f

v1 = fsolve(func, 0)
v2 = fsolve(func, 2.5)
fig, ax = plt.subplots()
ax.plot(x, func(x))
ax.plot(x, np.zeros(x.shape))
ax.plot(v1, func(v1), "o")
ax.plot(v2, func(v2), "o")
fig.savefig("fsolve.png")
plt.show()

2つの関数の交点を求める。

\[ \begin{align}\begin{aligned}&x^2 + y^2 = 2\\&y = x^2 - 4x + 1\end{aligned}\end{align} \]

example-02¶

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

x = np.linspace(-10, 10, 1000)
def func(p):
    x, y = p[0], p[1]
    f = np.zeros(2)
    f[0] = x**2 + y**2 - 2
    f[1] = x**2 - 4*x + 1 - y
    return f

o1 = fsolve(func, [1, 1])
o2 = fsolve(func, [-1, -1])
# x**2 + y**2 = 2
theta = np.linspace(0, 2*np.pi, 1000)
r = np.sqrt(2)
x = np.sqrt(2) * np.cos(theta)
y = np.sqrt(2) * np.sin(theta)
# x**2 - 4*x2 + 1
x2 = np.linspace(-2, 2, 1000)
y2 = x2**2 - 4*x2 + 1
fig, ax = plt.subplots()
ax.plot(x, y)
ax.plot(x2, y2)
ax.plot(o1[0], o1[1], "o")
ax.plot(o2[0], o2[1], "o")
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
fig.savefig("fsolve2.png")
plt.show()

curve_fit (leastsq)¶

最小自乗法でパラメータ \(a, b, c\) の最適値を決める。

\[f(x) = a\exp(bx) + cx^2\]

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html

example-03¶

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

x = np.linspace(0, 5, 20)
y = 3*np.exp(-2*x) + 0.11*x**2
y = y + (np.random.random(20) - 1)*0.5
fig, ax = plt.subplots()

func = lambda x, a, b, c: a*np.exp(b*x) + c*x**2

p0 = [0, 0, 0] # 初期値
p, cov = curve_fit(func, x, y, p0)
print(p)   # [ 2.81148394 -2.62913599  0.09389399]
print(cov) # [[  3.15666320e-02  -2.55900954e-02   7.83101093e-06] 
           #  [ -2.55900954e-02   1.03468885e-01  -4.74614449e-05] 
           #  [  7.83101093e-06  -4.74614449e-05   1.25025906e-05]]

ax.plot(x, y, "ro", label="exp")
ax.plot(x, func(x, p[0], p[1], p[2]), "r-", label="fit")
fig.savefig("curve_fit.png")
plt.show()

パラメータの境界制限を指定してフィッテングを行う。

example-04¶

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

x = np.linspace(0, 50, 100)
y = 10*np.exp(-0.2*x) - 1*np.exp(-0.05*x) 
y = y + (np.random.random(x.size) - 0.5)*0.1

f = lambda x, a, b: a*np.exp(b*x)
func = lambda x, a, b, c, d: f(x, a, b) + f(x, c, d)

# 境界制限なし
p0 = [0, 0, 0, 0]
p, cov = curve_fit(func, x, y, p0)
print(p)
fig, ax = plt.subplots(ncols=1, nrows=2)
ax[0].plot(x, y, "ro", label="exp")
ax[0].plot(x, func(x, *p), "r-", label="fit") # *p -> p[0], p[1], p[2], p[3]
ax[0].plot(x, f(x, p[0], p[1]), "-", label="exp1")
ax[0].plot(x, f(x, p[2], p[3]), "-", label="exp2")
ax[0].legend()
ax[0].set_title("non constraint")

# 境界制限付き
b = [(-np.inf, -np.inf, -np.inf, -np.inf), (np.inf, np.inf, 100, np.inf)]
p, cov = curve_fit(func, x, y, p0, bounds=b)
print(p)
ax[1].set_title("constraint")
ax[1].plot(x, y, "ro", label="exp")
ax[1].plot(x, func(x, *p), "r-", label="fit") # *p -> p[0], p[1], p[2], p[3]
ax[1].plot(x, f(x, p[0], p[1]), "-", label="exp1")
ax[1].plot(x, f(x, p[2], p[3]), "-", label="exp2")
ax[1].legend()
fig.savefig("curve_fit_const.png")
plt.show()

minimize¶

パラメータ同士の関係またはパラメータの範囲に拘束条件を課して最適値を求める。

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html

\[\begin{split}&-3x^3 - x^2y + \frac{1}{y} = 3z - 20\\ &x + y \gt z\\ &x^2 + y^2 = -20\\ &-10 \le x \le 30\\ &-10 \le y \le 30\\ &-10 \le z \le 30\\\end{split}\]

example-05¶

import numpy as np
from scipy.optimize import minimize


def func(p):
    x, y, z = p
    res = -3 * x**3 - x**2 * y  - 1 / (x*y) - 3 * z + 20
    return res


def cond_ineq(p):
    x, y, z = p
    return x + y - z


def cond_eq(p):
    x, y, z = p
    return x**2 + y**2 + 20

bounds = ((-10, 30), (-10, 30), (-10, 30))
conds = ({'type': 'ineq', 'fun': cond_ineq},
         {'type': 'eq', 'fun': cond_eq})

p0 = [5, 5, 5]  # 初期値
o = minimize(func, p0, constraints=conds, bounds=bounds)
print("func", func(o.x))
print("cond_ineq", cond_ineq(o.x))
print("cond_eq", cond_eq(o.x))

数値積分¶

以下の関数について、数値解をscipy.integrate.quad()を使って定積分を行います。

https://docs.scipy.org/doc/scipy/reference/integrate.html

\[f(x) = \int_0^2 (2x^2 + 3x - 1)dx\]

example-06¶

from scipy.integrate import quad

def func(x):
    return 2*x**2 + 3*x - 1

v, err = quad(func, 0, 1)
print(v)

積分範囲に \(\infty\) を含む定積分も可能です。関数い追加の引数を与えたい場合は、quad()にargsで渡します。gauss積分をscipyで計算して解析解と比較してみます。

\[\int_{-\infty}^\infty \exp(-a(x-b)^2)dx\]

example-07¶

from scipy.integrate import quad
import numpy as np


def func(x, a, b):
    return np.exp(-a * (x-b)**2)


a, b = 3, 2
v, err = quad(func, -np.infty, np.infty, args=(a, b))
print("scipy integral: {:.6f}".format(v))
print("gauss integral: {:.6f}".format(2 * np.sqrt(np.pi / a) / 2))

注釈

記号を定義して代数計算(因数分解や微分積分、微分方程式を解くことまで!) 行う方法は、Sympy(scipy.Symbolとは別パッケージ)でまとめて扱います。 Sympyが使えれば大学1年くらいまでの数学演習は全部pythonで解けるかも?

特殊関数¶

面倒な特殊関数もscipyを使えば簡単に計算してくれます。

https://docs.scipy.org/doc/scipy/reference/special.html

特殊関数で解析解が得られる水素原子の波動関数を計算するために、ルジャンドル(Legendre)陪多項式、ラゲール(Laguerre)多項式、球面調和関数 (Spherical harmonics)を計算してみます。

https://ja.wikipedia.org/wiki/%E6%B0%B4%E7%B4%A0%E5%8E%9F%E5%AD%90

\[\begin{split}&\Psi(r, \theta, \phi) = \sqrt{ \left(\frac{2}{na_0}\right)^3 \frac{(n-l-1)!}{2n(n+l)!} } \exp\left(-\frac{\rho}{2}\right) \rho^l L^{2l+1}_{n-l-1}(\rho) Y_{lm}(\theta, \phi)\\ &\rho = \frac{2r}{na_0}\\ &Y_{lm}(\theta, \phi) = \left\{ \begin{array}{ll} \frac{i}{\sqrt{2}}\left(Y_l^m - (-1)^mY_l^{-m} \right) & ({\rm if}\ m<0) \\ Y_l^m & ({\rm if}\ m=0)\\ \frac{1}{\sqrt{2}}\left(Y_l^{-m} - (-1)^mY_l^{m} \right) & ({\rm if}\ m>0) \\ \end{array} \right.\end{split}\]

ルジャンドル(Legendre)陪多項式¶

球面調和関数 \(Y^l_m (m=0)\) の結果に等しいルジャンドルの多項式です。

https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.legendre.html

直行多項式であるため量子力学でたまに出てきます。以下の微分方程式を満たす \(P_n(x)\) です。

\[ \begin{align}\begin{aligned}P_v^m = (-1)^m(1-x^2)^{m/2}\frac{d^m}{dx^m}P_v(x)\\P_v = \sum_{k=0}^{\infty}\frac{(-v)_k(v+1)_k}{(k!)^2}\left(\frac{1-x}{2}\right)^k\end{aligned}\end{align} \]

\(x=\cos(\theta)\) で \(P_n(x)\) を計算してみます。

example-08¶

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import lpmv

# scipy.special.lpmv(m, v, x)
# order m(磁気量子数; m), degree v(方位量子数; l)
theta = np.linspace(0, 2*np.pi, 1000)
x = np.cos(theta)
fig, ax = plt.subplots()
ax.plot(theta * 180/np.pi, lpmv(0, 1, x), label="$l=1, m=0$")
ax.plot(theta * 180/np.pi, lpmv(1, 1, x), label="$l=1, m=1$")
ax.plot(theta * 180/np.pi, lpmv(0, 2, x), label="$l=2, m=0$")
ax.plot(theta * 180/np.pi, lpmv(1, 2, x), label="$l=2, m=1$")
ax.plot(theta * 180/np.pi, lpmv(2, 2, x), label="$l=2, m=2$")
ax.legend()
ax.set_xlabel("angle/deg.")
ax.set_ylabel("$P^m_l$")
fig.savefig("legendre.png")
plt.show()

水素の波動関数の関数系である \(x = \cos(\theta)\) として、lとmのパターンをいくつか計算してプロットしてみます。

ラゲール(Laguerre)多項式¶

水素原子の波動関数の動径方向の減衰を決める関数です。

https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.laguerre.html

球面調和関数(Spherical harmonics)¶

球の表面形状を記述するための直行多項式です。軌道の形を決めために使われます。degree(方位量指数lに相当)とorder(磁気量子数mに相当)、極座標の \(\theta\) と \(\phi\) を与えると、複素数として \(Y^m_n\) が得られます。 \(P^m_n(\cos(\theta)\) はルジャンドルの多項式です。

\[Y^m_n(\theta, \phi) = \sqrt{ \frac{2n+1}{4\pi} \frac{(n-m)!}{(n+m)!} } \exp(i\theta) P_n^m(\cos(\phi))\]

https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.sph_harm.html

scipy.special.sph_harm(degree, order, theta, phi)を呼び出せば計算してくれます。引き数は以下のとおり。

order: 方位量子数m。 \(Y^m_n\) の \(m\)
degree: 磁気量子数l。\(Y^m_n\) の \(n\)
theta: 極座標の偏角。一般に記号として \(\phi [0, 2\pi]\) を使うので注意
phi: 極座標の偏角。一般に記号として \(\theta [0, \pi]\) を使うので注意

m, n, theta, phiを与えると球の表面までの距離が得られます。matplotlibで3次元プロットしてみます。実空間では実数のみ扱うので、realパートだけプロットします。軌道の形(\(\theta\) と \(\phi\) 依存性)は、球面調和関数で表現されます。

直行座標で可視化するためには、 \(|Y^m_n|\) の先端がなぞる面を表示します。 \(Y^m_n\) の正負は、色を変えて表現することとします。

example-09¶

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import sph_harm
from mpl_toolkits.mplot3d import Axes3D

def PlotData(Ylm, theta, phi, lm, outfile=None):
    """
    データのプロット
    """
    l, m = lm[0], lm[1]
    # |Ylm|がなぞるベクトルの先端を計算
    r = np.abs(Ylm.real)
    x = r * np.sin(theta) * np.cos(phi)
    y = r * np.sin(theta) * np.sin(phi)
    z = r * np.cos(theta)
    color = np.ones(r.shape)  # 正負の色つけのため   
    color[Ylm < 0] = -1       # 正負の色つけのため
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.plot_surface(x, y, z, facecolors=plt.cm.jet(color),
                    rstride=1, cstride=1)
    plt.title("$Y^{{m={}}}_{{l={}}}$".format(m, l), size=20)
    length = 0.5
    ax.set_xlim(-length, length)
    ax.set_ylim(-length, length)
    ax.set_zlim(-length, length)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    plt.setp(ax.get_xticklabels(), visible=False)
    plt.setp(ax.get_yticklabels(), visible=False)
    plt.setp(ax.get_zticklabels(), visible=False)
    if outfile is not None:
        fig.savefig(outfile)
        print(outfile, "was created.")

def CalcYlm(lm):
    l, m = lm[0], lm[1]
    # order=m degree=lで球面調和関数を計算
    t = np.linspace(0, np.pi, 50)
    p = np.linspace(0, 2 * np.pi, 100)
    theta, phi = np.meshgrid(t, p)
    if m == 0:
        Y = sph_harm(m, l, phi, theta)
    if m > 0:
        Ya = sph_harm(-m, l, phi, theta)
        Yb = sph_harm(m, l, phi, theta)
        Y = 1/np.sqrt(2) * (Ya + (-1)**m*Yb)
    if m < 0:
        Ya = sph_harm(m, l, phi, theta)
        Yb = sph_harm(-m, l, phi, theta)
        Y = 1j/np.sqrt(2) * (Ya - (-1)**m*Yb)
    return Y, theta, phi

# (l, m)の組み合わせを好きなだけ
lmset = [(1, -1), (1, 0), (1, 1),
         (2, -2), (2, -1), (2, 0), (2, 1), (2, 2),
         (3, -3), (3, -2), (3, -1), (3, 0), (3, 1), (3, 2), (3, 3)]

for lm in lmset:
    Y, theta, phi = CalcYlm(lm)
    # lとmから出力ファイル名を決める
    outfile = "l{}_m{:+d}.png".format(lm[0], lm[1])
    PlotData(Y, theta, phi, lm, outfile)
# plt.show()

プログラムを実行するとp, d, f軌道に相当する球面調和関数の図を作成します。画像をならべた結果は以下のとおりです。

水素の波動関数¶

\(n, l, m\) を与えて \(|\psi^2|\) を計算します。 matplotlibは、\(|\psi^2(x, y, z)|\) のような3次元の密度データをプロットできないので、可視化はgaussian cubeファイルを作成してvestaで行うこととします。cubeファイルは、小さな体積に空間を分割(voxel)して各voxelに密度データを割り当てるグリッドデータです。空間(voxel)を定義するためのベクトルと分割数を必要とします。

あまり教科書にでてこないg軌道(n=6, l=5, m=0,1,2,3,4,5)のcubeファイルを作成してみます。m=-1,-2,-3,-4,-5は位相が変わるだけで軌道の形は同じです。

example-10¶

import numpy as np
from scipy.special import sph_harm
from scipy.special import factorial
from scipy.special import assoc_laguerre

a0 = 0.529177  # global, Distance: Bohr unit a0 = 0.529177 Ang.

def OutputCube(nlm, L, psi2, outfile):
    """
    gaussian cube fileで書き出し
    a, b, c vectorはBohr単位で指定する
    """
    d = (L[1] - L[0]) / L[2] / a0 # Bohr単位
    p = (L[1] - L[0])/a0 * 0.5
    o = open(outfile, "w")
    o.write("hydrogen orbital cube file\n")
    o.write("(n, l, m) = {} {} {}".format(nlm[0], nlm[2], nlm[2]))
    o.write("\n")
    o.write("{:6d} {:8f} {:8f} {:8f}\n".format(1, 0, 0, 0)) # 原子数、origin
    o.write("{:6d} {:8f} {:8f} {:8f}\n".format(L[2], d, 0, 0)) # a vector bohr
    o.write("{:6d} {:8f} {:8f} {:8f}\n".format(L[2], 0, d, 0)) # b vector bohr
    o.write("{:6d} {:8f} {:8f} {:8f}\n".format(L[2], 0, 0, d)) # c vector bohr
    o.write("{:6d} {:8f} {:8f} {:8f} {:8f}\n".format(1, 0, p, p, p)) # Hの座標
    psi2 = psi2.reshape(-1) * 100
    for i, rho in enumerate(psi2):
        o.write("{:12g} ".format(rho))
        if i % 6 == 5:
            o.write("\n")
    print(outfile, "was created.")

def CalcPsi(nlm, L):
    n, l, m = nlm[0], nlm[1], nlm[2]
    # voxcelデータの生成
    xi = np.linspace(L[0], L[1], L[2])
    yi = np.linspace(L[0], L[1], L[2])
    zi = np.linspace(L[0], L[1], L[2])
    x, y, z = np.meshgrid(xi, yi, zi)
    # 直行座標 -> 極座標
    r = np.sqrt(x**2 + y**2 + z**2)
    theta = np.arccos(z / r)
    phi = np.arctan2(y, x)
    # 係数
    rho = 2*r/a0/n
    coef1 = np.sqrt((2 / (n * a0))**3 * factorial(n - l - 1) /
                 (2 * n * factorial(n + l)))
    coef2 = np.exp(-rho/2) * rho**l
    # ラゲール(Laguerre)
    L = assoc_laguerre(rho, n - l - 1, 2 * l + 1)
    # 球面調和関数(spherical harmonic function)
    if m == 0:
        Y = sph_harm(m, l, phi, theta)
    if m > 0:
        Ya = sph_harm(-m, l, phi, theta)
        Yb = sph_harm(m, l, phi, theta)
        Y = 1/np.sqrt(2) * (Ya + (-1)**m*Yb)
    if m < 0:
        Ya = sph_harm(m, l, phi, theta)
        Yb = sph_harm(-m, l, phi, theta)
        Y = 1j/np.sqrt(2) * (Ya - (-1)**m*Yb)
    psi =coef1 * coef2 * L * Y
    psi2 = psi.real * np.conj(psi.real)
    # 色つけのための符号
    s = np.ones(psi.shape)
    s[Y < 0] = -1
    psi2 = psi2 * s
    return psi2

# 主量子数, 方位量子数, 磁気量子数
nlmset = [(6, 5, 0), (6, 5, 1), (6, 5, 2), (6, 5, 3),  (6, 5, 4), (6, 5, 5)]
nlmset = [(6, 5, -5), (6, 5, -4), (6, 5, -3), (6, 5, -2),  (6, 5, -1)]
# xyz軸の min(Ang.), max(Ang.), voxcelの分割数
L = (-30, 30, 100)
for nlm in nlmset:
    outfile = "nml{}{}{}.cube".format(nlm[0], nlm[1], nlm[2])
    psi2 = CalcPsi(nlm, L)
    OutputCube(nlm, L, psi2, outfile)

プログラムを実行して生成されるcubefileです。vestaで開いて確認してみてください。

vestaで可視化して作成した図をならべた結果です。

cubeファイルは少なくとも原子座標を1つ含まなければvestaで可視化できないようです。そのためH原子を1つだけ中心に置いています。上の図ではShow modelを offにしてH原子を表示していません。

vesta上でboundaryを指定すれば、任意の面でスライスすることもできます。(n,l,m)=(6,5,5)のデータをzmin=0.5からzmax=0.5でスライスした場合です。

クイズ¶

Q1¶

実験データ(xy_fit.dat)を2つのガウス関数でピーク分離することを考える。

\[y = p_0\exp\left(-\frac{x-p_1}{p_2}\right) + p_3\exp\left(-\frac{x-p_4}{p_5}\right)\]

ただし2つのガウス関数の面積は等しいとして以下の拘束条件を満たすように \(p\) を求めます。ガウス関数の面積はquad()で求めます。

\[\int_{-\infty}^{\infty} p_0\exp\left(-\frac{x-p_1}{p_2}\right) = \int_{-\infty}^{\infty} p_3\exp\left(-\frac{x-p_4}{p_5}\right)\]

minimizeは初期値が最適値から遠いとうまく最適値を求めてくれません。そこで、最初にleast_squareで拘束条件なしでパラメータを求め、求めたパラメータを初期値としてminimizeで拘束条件付きでフィッテングします。

結果は以下のとおりです。

least_squares fitting:
p[0] = 2.01
p[1] = -2.00
p[2] = 1.00
p[3] = 1.00
p[4] = 0.99
p[5] = 3.00
area ratio gauss1/gauss2: 1.16
minimize fitting:
p[0] = 1.97
p[1] = -2.03
p[2] = 0.92
p[3] = 0.99
p[4] = 0.90
p[5] = 3.63
area ratio gauss1/gauss2: 1.00

注釈

quadで積分するのはコストがかかります。面積で規格化された下式のガウス関数を使えば、わざわざ数値積分する必要がないので、効率的にフィッテングできるでしょう。

\[\begin{split}&f(x) = p_0\sqrt{\frac{\log(2)}{\pi}}\frac{1}{p_2}\exp\left[ -\log(2)\left\{\frac{x - p_1}{p_2} \right\}^2 \right]\\ &\int_{-\infty}^{\infty} f(x) dx = p_0\end{split}\]

Q2¶

量子数が \((n, l, m) = (3, 2, 0)\) の場合、原子核から3.0 Ang.以内の距離に電子が存在する確率を計算する。微小体積の確率を求めて足し合わせで求めてください。3.0 Ang.以内の領域を分割する分割数 \(N\) を10, 50, 100, 200と増やして2桁の精度が確認できれば十分です。

N=10 7.77%
N=50 8.85%
N=100 8.88%
N=200 8.89%

答え 8.9%