nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。

立体角積分 - Lebedev quadrature

数値計算で立体角積分をしたいとき、安直に二重積分するよりも、Lebedev求積法を使った方が効率が良い。
Lebedev quadrature - Wikipedia
http://people.maths.ox.ac.uk/beentjes/Essays/QuadratureSphere.pdf

効率が良いという意味は、メッシュ数が少なくても精度が出るという意味で、言い換えれば同じ精度を出すのに二重積分よりもずっと速く計算が出来るということである。

求積法は、つまるところ、積分を「上手く精度が出るようにメッシュ間隔等を調整しかつ各点に重み付けをした和」に変換したものである。
Lebedev法の場合には、

 \displaystyle
\int sin \theta \, d\theta \, \int d\phi \, f(\theta,\phi) \approx 4 \pi \sum^N_{i=1} \omega_i \, f(\theta_i, \phi_i )

というように 4 \pi は後で掛けることを想定している。

Lebedev求積法のFortranのフリーのソースコードとしては、以下のものが知られている。
http://ccl.net/cca/software/SOURCES/FORTRAN/Lebedev-Laikov-Grids/Lebedev-Laikov.F

subroutineの構成としては、

get_oh( ... )
LD0006( X, Y, Z, W, N )
LD0014( X, Y, Z, W, N )
...
LD5810( X, Y, Z, W, N )

という感じで、使うのはLDXXXX(...)の方である。get_ohは外から使うことは無い。

使い方のポイントを箇条書きすると、

  • 引数の"X, Y, Z, W, N"は全てoutputであるため、こちらから何か値を用意する必要は無い。
  • ただし、"X, Y, Z, W"は配列サイズを事前に指定する必要がある
  • 例えばLD0006(...)の"0006"はメッシュ数を表していて、これを使いたければ少なくとも6より大きい配列サイズの"X, Y, Z, W"を用意しておく。
  • "N"はメッシュ数であるためLD0006(...)ではN=6が返って来る。
  • "X Y Z"がLebedevメッシュであり、 ( \theta, \phi )に直したければ、"acos(z)"や"atan2(y,x)"で自分で変換する。
  • "W"が重み付けで、"sum( W )"はほぼ1になるようになっている。
  • どのLDXXXX(...)を使うべきかは、計算時間と計算精度との相談になる。

gfortranでコンパイルしてWの和が1になるかチェックする簡単なテストプログラム。

program check_Lebedev
implicit none
double precision, allocatable :: x(:), y(:), z(:), w(:)
! double precision, allocatable :: theta(:), phi(:)
integer :: n, n_lbdv

write( *, * ) 'number of mesh:'
read( *, * ) n
allocate( x(n), y(n), z(n), w(n) )
select case( n )
    case(    6 )
        call LD0006( x, y, z, w, n_lbdv )
    ...
    case( 5810 )
        call LD5810( x, y, z, w, n_lbdv )
    case default
        write( *, * ) ' n /= n_lbdv'
        stop
end select

write( *, * ) sum( w )

! allocate( theta( n_lbdv ), phi( n_lbdv ) )
! theta = acos( z )
! phi = atan2( y, x )

end program check_Lebedev

include "Lebedev-Laikov_CommentOutBy!.F"

リンク先のLebedev-Laikov.Fはコメントアウトが行頭の"c"で行われているため、sedか何かを使って"!"でコメントアウトするようにする必要がある(適当なコンパイルオプションを使えば、そんなことしなくても良いのかもしれないが)。

これで球面調和関数とかを掛けた計算とかが比較的簡単に出来るようになると思われる。

Møller wave operator における 加藤の連鎖律(?)

N粒子系における二体ポテンシャル Vを、粒子ペアに番号を付けて、その間の相互作用の和( M=N(N-1)/2)として定義する。

 \displaystyle
\begin{align}
V &\equiv \sum^{M}_{k=1} v_k = V_i + \bar{ V }_i = V_M = \bar{ V_0 }\\
V_i &\equiv \sum^{ i }_{k=1} v_k \\
\bar{ V }_i &\equiv \sum^{M}_{k= i + 1 } v_k
\end{align}

二体ポテンシャルとしたが、結局は相互作用に何かしら番号が付けられて、全体がその和として与えらえるのであれば、いろんな応用が考えられる(と思う)。
このような相互作用に対するN粒子系のGreen関数を以下のように定義。

 \displaystyle
\begin{align}
G &\equiv ( z - H_0 - V )^{-1} = G_M \\
G_i &\equiv ( z - H_0 - V_i )^{-1}
\end{align}

ここで次の量を考える。

 \displaystyle
\begin{align}
1 + G \bar{ V }_i &= G G^{-1} \left( 1 + G \bar{ V }_i \right)
    = G \left( G^{-1} + \bar{ V }_i \right) = G \left(  z - H_0 - V + \bar{ V }_i \right) \\
    &= G \left(  z - H_0 - V_i \right) = G \left( G_i \right)^{-1} \\
1 + G_i v_i &= G_i \left( G_i \right)^{-1} \left( 1 + G_i v_i \right)
    = G_i \left( \left( G_i \right)^{-1} + v_i \right) = G_i \left(  z - H_0 - V_i + v_i \right) \\
    &= G_i \left(  z - H_0 - V_{i-1} \right) = G_i \left( G_{i-1} \right)^{-1}
\end{align}

したがって、 G = G_Mに注意すると、

 \displaystyle
\begin{align}
G_M \left( G_i \right)^{-1} 
    &= G_M \left( G_{M-1} \right)^{-1} G_{M-1} \left( G_i \right)^{-1} \\
    &= G_M \left( G_{M-1} \right)^{-1} G_{M-1} \cdots \left( G_{i+1} \right)^{-1} G_{i+1} \left( G_i \right)^{-1}
\end{align}

すなわち、

 \displaystyle
\begin{align}
1 + G_M \bar{ V }_i
    &= \left( 1 + G_M v_M \right) \left( 1 + G_{M-1} v_{M-1} \right) \cdots \left( 1 + G_{i+1} v_{i+1} \right)
\end{align}

これを i = 0の場合に適用すれば、

 \displaystyle
\begin{align}
1 + G V &= 1 + G \sum^{M}_{k=1} v_k \\
    &= \left( 1 + G_M v_M \right) \left( 1 + G_{M-1} v_{M-1} \right) \cdots \left( 1 + G_1 v_1 \right)
\end{align}

という様に、和を積で表現し直すことが出来る。

 1 + GVはMøller wave operator  \Omegaとして知られており、 \Omega_i = 1 + G_i v_i と定義すれば、

 \displaystyle
\begin{align}
\Omega = \prod^{M}_{k=1} \Omega_k
\end{align}

と表せる。これを加藤の連鎖律と呼ぶらしい。

相互作用のペアを大きめに取っておいて、上手い具合に分割したら上手いこと相互作用を繰り込んだりして遊べそうな気がする。

球対称関数を別の位置で球平均する

 rを原点から測った距離とし、原点から見て球対称な関数を F( r )とする。
原点から見て位置 {\bf R}_iにあるサイト iがあるとする。
サイト iを原点に取り直した任意の位置ベクトルを {\bf r}_iと定義する。
やりたいことは、 F( r )をサイト iから見ると球対称ではないので、サイト iから見たときの球対称成分 f( r_i )を作りたいということである。
それはすなわちサイト i上で球平均することに対応する。
 {\bf r}_i = {\bf r} - {\bf R}_i、もしくは {\bf r} = {\bf r}_i + {\bf R}_i(こっちの方が直観的かも)であるから、

 \displaystyle
\begin{align}
f( r_i ) = \left( \int d \hat{\bf r}_i \, F( | {\bf r}_i + {\bf R}_i | ) \right) / \left( \int d \hat{\bf r}_i \right)
    = \frac{ 1 }{ 4 \pi } \int d \hat{\bf r}_i \, F( | {\bf r}_i + {\bf R}_i | )
\end{align}

球対称であるから、 {\bf R}_i z軸上に取っても結果は変わらない。
この時、 F( | {\bf r}_i + {\bf R}_i | )は方位角に依存しないから、

 \displaystyle
\begin{align}
f( r_i ) = \frac{ 1 }{ 2 } \int^{\pi}_{0} {\rm sin} \theta \, d \theta \, F( | {\bf r}_i + {\bf R}_i | )
\end{align}

ここで、 \theta {\bf r}_i - {\bf R}_iの成す角である。
マイナスが気になるかも知れないが、最終的にはその動径成分の R_iのみを使うので、符号はあまり重要でない。
 | {\bf r}_i + {\bf R}_i | = \sqrt{ r^2_i + R^2_i - 2 r_i R_i {\rm cos} \theta }であることに注意すると、

 \displaystyle
\begin{align}
&\frac{ d | {\bf r}_i + {\bf R}_i | }{ d \theta }
    = \frac{ r_i R_i }{ \sqrt{ r^2_i + R^2_i - 2 r_i R_i {\rm cos} \theta } } {\rm sin} \theta
    = \frac{ r_i R_i }{ | {\bf r}_i + {\bf R}_i | } {\rm sin} \theta
\\
&\therefore \; {\rm sin} \theta \, d \theta
    = \frac{ | {\bf r}_i + {\bf R}_i | }{ r_i R_i  } \, d | {\bf r}_i + {\bf R}_i | 
\end{align}

よって、

 \displaystyle
\begin{align}
f( r_i ) &= \frac{ 1 }{ 2 r_i R_i } \int^{ R_i + r_i }_{ R_i - r_i } \, d | {\bf r}_i + {\bf R}_i | \, | {\bf r}_i + {\bf R}_i | \, F( | {\bf r}_i + {\bf R}_i | ) \\
    &= \frac{ 1 }{ 2 r_i R_i } \int^{ R_i + r_i }_{ R_i - r_i } \, dr \, r \, F( r )
\end{align}

 r_iの関数として f( r_i )を求めたいので、ここでの扱いは言わばパラメータのように外から与えられるものであり、積分変数ではないため積分の外に出せる。

例として、クーロンポテンシャル F( r ) = Z / r \equiv V( r )を考えると、

 \displaystyle
\begin{align}
V( r_i ) &= \frac{ Z }{ 2 r_i R_i } \left( R_i + r_i - ( R_i - r_i ) \right) = \frac{ Z }{ R_i } \equiv \bar{ V }
\end{align}

となり、違うサイトから見ると定数ポテンシャルとして与えられることになる。
したがって、点電荷が一次元に正負(絶対値は同じ大きさ Z)で交互に等間隔 Rで無限に並んでいて、ある点電荷に働くトータルの静電ポテンシャルは、球平均すれば

 \displaystyle
\begin{align}
V = - \frac{ 2 Z }{ R } \left( 1 - \frac{ 1 }{ 2 } + \frac{ 1 }{ 3 } - \cdots \right) = - \frac{ 2 Z }{ R } {\rm ln}2
\end{align}

と求めることが出来る。
ただ、これは ln( 1 + x )マクローリン展開、つまり x \approx 0に対して速く収束する展開に、 x = 1を代入して得られたものを使用しているわけだから、本当に点電荷がほぼ無限に並んでいないと(ほぼ無限に展開しないと)値の一致は悪いと思われる。クーロンポテンシャルの収束が遅いのがよく分かる。

デルタ関数のエンタングルメントを完備関係式で断つ?

例えば球面調和関数

 \displaystyle
\begin{align}
\delta( \hat{\bf r} - \hat{\bf r}' )
    &= \delta( \hat{\bf r}' - \hat{\bf r} ) \\
    &= \sum_L Y_{L} ( \hat{\bf r} ) Y^*_L ( \hat{\bf r}' ) 
    = \sum_L Y^*_{L} ( \hat{\bf r} ) Y_L ( \hat{\bf r}' ) \\
    &= \sum_L Y_{L} ( \hat{\bf r}' ) Y^*_L ( \hat{\bf r} ) 
    = \sum_L Y^*_{L} ( \hat{\bf r}' ) Y_L ( \hat{\bf r} ) 
\end{align}

これはある意味 \hat{\bf r},\, \hat{\bf r} エンタングルメントを断ったと言えるか?
添字が一個なので、ベクトルに拡張すれば、

 \displaystyle
\begin{align}
\delta( \hat{\bf r} - \hat{\bf r}' )
 = \left(
        Y_{L_1} ( \hat{\bf r} ), Y_{L_2} ( \hat{\bf r} ), \cdots
    \right)
    \left(
        \begin{array}{c}
            Y^*_{L_1}( \hat{\bf r}' ) \\ Y^*_{L_2}( \hat{\bf r}' ) \\ \vdots
        \end{array}
    \right)
 = \left(
        Y^*_{L_1} ( \hat{\bf r} ), Y^*_{L_2} ( \hat{\bf r} ), \cdots
    \right)
    \left(
        \begin{array}{c}
            Y_{L_1}( \hat{\bf r}' ) \\ Y_{L_2}( \hat{\bf r}' ) \\ \vdots
        \end{array}
    \right) \\
 = \left(
        Y_{L_1} ( \hat{\bf r}' ), Y_{L_2} ( \hat{\bf r}' ), \cdots
    \right)
    \left(
        \begin{array}{c}
            Y^*_{L_1}( \hat{\bf r} ) \\ Y^*_{L_2}( \hat{\bf r} ) \\ \vdots
        \end{array}
    \right)
 = \left(
        Y^*_{L_1} ( \hat{\bf r}' ), Y^*_{L_2} ( \hat{\bf r}' ), \cdots
    \right)
    \left(
        \begin{array}{c}
            Y_{L_1}( \hat{\bf r} ) \\ Y_{L_2}( \hat{\bf r} ) \\ \vdots
        \end{array}
    \right)
\end{align}

よって、ある種の行列化によって変数分離が可能である。ポイントは、もともと行列ではなくただの数だったので、「左に横ベクトル・右に縦ベクトル」の順番は変更出来ないという点である。
あまりエンタングルメントについて詳しくない(というか無知に等しい)が、エンタングルメントを断つ(解す)とは多分こういうことなのかと。
上では非局所相関を局所相関の積で表すというようなことをしたわけだが、この結果を用いると、いわゆる多重極展開と呼ばれる、動径成分と立体角成分に分ける展開が可能になる。
(これもある種のエンタングルメント解しか?)

 \displaystyle
\begin{align}
F( r, \hat{\bf r} )
    &= \int F( r, \hat{\bf r}' ) \delta( \hat{\bf r} - \hat{\bf r}' ) d \hat{\bf r}'
    = \int F( r, \hat{\bf r}' ) \sum_L Y_{L} ( \hat{\bf r} ) Y^*_L ( \hat{\bf r}' ) d \hat{\bf r}' \\
    &= \sum_L \left( \int F( r, \hat{\bf r}' ) Y^*_L ( \hat{\bf r}' ) d \hat{\bf r}' \right) Y_{L}( \hat{\bf r} ) 
    \equiv \sum_L f_L( r ) Y_{L}( \hat{\bf r} ) \\
    &= \left(
        f_{L_1}( r ), f_{L_2}( r ), \cdots
    \right)
    \left(
        \begin{array}{c}
            Y_{L_1}( \hat{\bf r} ) \\ Y_{L_2}( \hat{\bf r} ) \\ \vdots
        \end{array}
    \right)
\end{align}

これは「連続量の相関を離散量の話に置き換える」と言えるだろうか?
ちなみに逆もおそらく可能(離散量の相関を連続量の話に置き換える)。

 \displaystyle
\begin{align}
\delta_{ L, L' }
    &= \int d \hat{\bf r} Y^*_L( \hat{\bf r} ) Y_{L'}( \hat{\bf r} ) \\
g_{N L}
    &= \sum_{L'} g_{NL'} \delta_{L' L}
    = \sum_{L'} g_{N L'} \int d \hat{\bf r} Y^*_{L'}( \hat{\bf r} ) Y_{L}( \hat{\bf r} ) \\
    &= \int d \hat{\bf r} \left( \sum_{L'} g_{N L'} Y^*_{L'}( \hat{\bf r} ) \right) Y_{L}( \hat{\bf r} )
    \equiv \int d \hat{\bf r} G_{N}( \hat{\bf r} ) Y_{L}( \hat{\bf r} )
\end{align}

連続と離散の絡み合いは楽しいので、もっと深めて行きたい。

Gaunt積分のまとめ

Gaunt積分の自分的メモ。
Clebsch–Gordan coefficients - Wikipedia
3-j symbol - Wikipedia

Gaunt積分(多分一般的な呼び名ではない。が、自分の分野ではよく出てくる。)
 \displaystyle
\begin{align}
 & \int d\hat{\bf r} \, Y^*_{ l_3 m_3 }(\hat{\bf r}) Y_{ l_1 m_1 }(\hat{\bf r}) Y_{ l_2 m_2 }(\hat{\bf r}) \\
 & \quad = \sqrt{ \frac{ (2 l_1 + 1 ) ( 2 l_2 + 1 ) }{ 4 \pi ( 2 l_3 + 1 ) } }
        < l_1,0;l_2,0|l_3,0 > < l_1,m_1;l_2,m_2|l_3,m_3 >
\\
 & \quad = (-1)^{ m_3 } \sqrt{ \frac{ (2 l_1 + 1 ) ( 2 l_2 + 1 ) ( 2 l_3 + 1 ) }{ 4 \pi } }
        \left( \begin{array}{ccc}  l_1 & l_2 & l_3 \\ 0 & 0 & 0 \end{array} \right)
        \left( \begin{array}{ccc}  l_1 & l_2 & l_3 \\ m_1 & m_2 & - m_3 \end{array} \right)
\end{align}

Gaunt係数(Wikipediaとかに載っているのはこっち。)
 \displaystyle
\begin{align}
 & \int d\hat{\bf r} \, Y_{ l_1 m_1 }(\hat{\bf r}) Y_{ l_2 m_2 }(\hat{\bf r}) Y_{ l_3 m_3 }(\hat{\bf r}) \\
 & \quad = \sqrt{ \frac{ (2 l_1 + 1 ) ( 2 l_2 + 1 ) ( 2 l_3 + 1 ) }{ 4 \pi } }
        \left( \begin{array}{ccc}  l_1 & l_2 & l_3 \\ 0 & 0 & 0 \end{array} \right)
        \left( \begin{array}{ccc}  l_1 & l_2 & l_3 \\ m_1 & m_2 & m_3 \end{array} \right)
\end{align}

Clebsh-Gordan係数と3j記号の関係
 \displaystyle
\begin{align}
 < l_1,m_1;l_2,m_2|l_3,m_3 > =
        (-1)^{ l_1 - l_2 } (-1)^{ m_3 } \sqrt{ 2 l_3 + 1 }
        \left( \begin{array}{ccc}  l_1 & l_2 & l_3 \\ m_1 & m_2 & - m_3 \end{array} \right)
\end{align}

3j記号の列の交換
 \displaystyle
\begin{align}
\left( \begin{array}{ccc}  l_1 & l_2 & l_3 \\ m_1 & m_2 & m_3 \end{array} \right)
    &= (-1)^{ l_1 + l_2 + l_3 } \left( \begin{array}{ccc}  l_2 & l_1 & l_3 \\ m_2 & m_1 & m_3 \end{array} \right)
\\
    &= (-1)^{ l_1 + l_2 + l_3 } \left( \begin{array}{ccc}  l_1 & l_3 & l_2 \\ m_1 & m_3 & m_2 \end{array} \right)
\end{align}

行列要素を視覚的にチェックする方法

行列要素をなんとか可視化したくて、三次元プロット散布図で表現してみた。
Pythonでデータを可視化するmatplotlib初級(3D散布図編) - しぷぜん

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

N = 20
row = range( N )
column = range( N )

X, Y = np.meshgrid( row, column )
X = X.reshape( N * N )
Y = Y.reshape( N * N )

fig = plt.figure()
ax = fig.add_subplot( 111, projection='3d' )
for i in range(3):
    matrix_reshape_1d = np.random.rand( N * N )
    # default marker size = 20
    ax.scatter( X, Y, [i] * ( N * N ), s = matrix_reshape_1d * 20. )
plt.show()

f:id:koideforest:20180501190738p:plain

例えば行列全体が時間に依存していて、それの時間発展を見たい時には良さそう。
複素数の場合には、色で位相を表現することになると思うが、色は配列で渡せないので、一点一点をいちいちプロットしないと行けない気がする。行列サイズが大きいと、おそらく途中で限界が来るだろう。

Legendre多項式の直交性を python (scipy) で確かめる

いつもnumpyの多項式の使い方を忘れるので、Legendre多項式の直交関係式で練習する。
ルジャンドル多項式 - Wikipedia


\displaystyle
\int^{1}_{-1} P_{ l' }( x ) \, P_{ l }( x ) \, dx = \frac{ 2 }{ 2 l + 1 } \delta_{ l' l }

scipy.special.legendreでLegendre多項式を作ると、poly1dというタイプの関数が帰ってくるので、そのままxを代入できる(下記リンク or サンプルスクリプト参照)。
scipy.special.legendre — SciPy v1.0.0 Reference Guide
numpy.poly1d — NumPy v1.14 Manual

積分は一次元積分のquadを使用。
Python SciPy : SciPy の積分関数の基本的使い方 | org-技術

from scipy.special import legendre
from scipy.integrate import quad

def Plx( l, x ):
    # legendre( l ) >>> poly1d([ ... ])
    # poly1d([ ... ])( x ) >>> float
    return legendre( l )( x )

def orthogonal_relation( lp, l ):
    # quad >>> ( result, error )
    return quad( lambda x: Plx( lp, x ) * Plx( l, x ), -1., 1. ) [0]

if __name__ == '__main__'
    for lp in range( 5 ):
        for l in range( 5 ):
            print( 'lp, l = ', lp, l )
            print ( 'numerical', orthogonal_relation( lp, l ) )
            if lp == l:
                print( 'analytical', 2. / ( 2 * l + 1. ) )
            else:
                print( 'analytical', 0. )