2020-03-25

Gaussianの半値幅。

数学

規格化されたGaussianは次の形で与えられる。
$\displaystyle g( x ) = \frac{ 1 }{ \sigma \sqrt{ 2 \pi } } e^{ - \frac{ ( x - \mu )^2 }{ 2 \sigma^2 } }$

規格化されている（積分値が1になっている）ことは、次のようにして確かめられる。
$\displaystyle \int^\infty_{-\infty} dx \, e^{ - a x^2 } = \sqrt{ \frac{ \pi }{ a } } \\ \displaystyle \therefore \int^\infty_{-\infty} dx \, g( x ) = \frac{ 1 }{ \sigma \sqrt{ 2 \pi } } \sqrt{ \frac{ \pi }{ \frac{ 1 }{ 2 \sigma^2 } } } = 1$

$\sigma$ を用いた表記は、統計学的な処理を行う際に便利である（と想像される）が、直感的にどれくらいの幅を持つのかが分かりにくい。
そこで、半値半幅 $\Gamma$ を用いた表記を考えることにする。
$\Gamma$ の定義は、次で与えられる。
$\displaystyle \frac{ g( \Gamma + \mu ) }{ \max{ g } } = \frac{ g( \Gamma + \mu ) }{ g( \mu ) } = e^{ - \frac{ \Gamma^2 }{ 2 \sigma^2 } } = \frac{1}{2}$

したがって、 $\sigma \rightarrow \Gamma$ の変換は、
$\displaystyle - \frac{ \Gamma^2 }{ 2 \sigma^2 } = \ln{ \frac{1}{2} } = - \ln{2} \\ \displaystyle \Gamma^2 = 2 \sigma^2 \ln{2} \\ \displaystyle \therefore \Gamma = \sigma \sqrt{ 2 \ln{2} } \Leftrightarrow \sigma = \frac{ \Gamma }{ \sqrt{ 2 \ln{2} } }$
数値的に見積もると、 $\Gamma \approx 1.2 \sigma$ となる。半値半幅 $\Gamma$ は $\sigma$ よりも気持ちちょっと大きい程度である。

元のGaussianの式を $\Gamma$ で表すと、
$\displaystyle g( x ) = \frac{1}{ \Gamma } \sqrt{ \frac{ \ln{2} }{ \pi } } e^{ - \frac{ \ln{2} }{ \Gamma^2 } ( x - \mu )^2 }$

（おそらくLorentz関数でもそうだと思うが、）Gaussian型のピーク形状の場合、半値全幅 $2\Gamma$ （文献によっては、 $\Gamma$ を半値全幅の表記に用いる場合もあるため、注意が必要）を分解能として扱うため、Gaussian broadening等で理論値に実験誤差を入れ込む場合には、半値幅の用いた表記の方が便利である。

2020-02-28

NiOの反強磁性結晶構造における対称性をspglibを使って求める。

通常の結晶構造であれば、データベースか何かから引っ張って来れるが、それをベースにした超周期構造を作ろうと思うと、対称性（空間群）が変わるため、手でやるには歯が立たない。
そこで、入力した原子座標に対して対称性を見つけてくれる spglib を使って、反強磁性でのNiOの対称性を決定することを試みる。

spglibはC言語、FORTRAN, pythonに対応しているが、いずれにおいても「実行ファイル（今風に言えばアプリケーション）」ではなく、（~libとあるように）ライブラリーとして提供されている。そのためGUIのようなものに結晶構造を入れるようなものではないため、多少プログラミングに慣れている必要がある。
Spglib — spglib 1.12.0 documentation

個人的には、pythonが簡便でオススメである。なので、ここではpythonでの例を紹介する。
Spglib for Python — spglib 1.12.0 documentation
pythonバージョンのインストールは、pipですぐ入る。

pip install spglib

これを使って、まずは通常のNiOの結晶構造で試してみる。

import numpy as np
import spglib

# lattice constant
a = 4.19

# structural information ( [lattice constants], [coordination], [atomic numbers] )
# Caution: it should be a tuple
nio = \
    ( [ [ a, 0, 0 ], [ 0, a, 0 ], [ 0, 0, a ] ],
      [ [ 0,   0,   0   ],  # Ni
        [ 0,   0.5, 0.5 ],  # Ni
        [ 0.5, 0,   0.5 ],  # Ni
        [ 0.5, 0.5, 0   ],  # Ni
        [ 0.5, 0,   0   ],  # O
        [ 0,   0.5, 0   ],  # O
        [ 0,   0,   0.5 ],  # O
        [ 0.5, 0.5, 0.5 ]   # O
      ],
      [28,] * 4 + [8,] * 4
    )

# just space group 
print( spglib.get_spacegroup( nio ) )
#Fm-3m (225)

# primitive cell
print( spglib.find_primitive( nio ) )
#(array([[0.   , 2.095, 2.095],
#       [2.095, 0.   , 2.095],
#       [2.095, 2.095, 0.   ]]), array([[0. , 0. , 0. ],
#       [0.5, 0.5, 0.5]]), array([28,  8], dtype=int32))

$2 \times 2 \times 2$ のsuper cell でもちゃんと元の構造の対称性を返してくれる。

# some coordinates are just added 1 to discribe them in the other unit cell
# therefore, lattice constans and coordinates should be multipled and divided by 2, respectively
nio_sc \
    = [ [ ( a, 0, 0 ), ( 0, a, 0 ), ( 0, 0, a ) ],
        [ [ 0,   0,   0   ],  # Ni 0, 0, 0
          [ 0,   0.5, 0.5 ],
          [ 0.5, 0,   0.5 ],
          [ 0.5, 0.5, 0   ],
          [ 1,   0,   0   ],  # Ni 1, 0, 0
          [ 1,   0.5, 0.5 ],
          [ 1.5, 0,   0.5 ],
          [ 1.5, 0.5, 0   ],
          [ 0,   1,   0   ],  # Ni 0, 1, 0
          [ 0,   1.5, 0.5 ],
          [ 0.5, 1,   0.5 ],
          [ 0.5, 1.5, 0   ],
          [ 0,   0,   1   ],  # Ni 0, 0, 1
          [ 0,   0.5, 1.5 ],
          [ 0.5, 0,   1.5 ],
          [ 0.5, 0.5, 1   ],
          [ 1,   1,   0   ],  # Ni 1, 1, 0
          [ 1,   1.5, 0.5 ],
          [ 1.5, 1,   0.5 ],
          [ 1.5, 1.5, 0   ],
          [ 1,   0,   1   ],  # Ni 1, 0, 1
          [ 1,   0.5, 1.5 ],
          [ 1.5, 0,   1.5 ],
          [ 1.5, 0.5, 1   ],
          [ 0,   1,   1   ],  # Ni 0, 1, 1
          [ 0,   1.5, 1.5 ],
          [ 0.5, 1,   1.5 ],
          [ 0.5, 1.5, 1   ],
          [ 1,   1,   1   ],  # Ni 1, 1, 1
          [ 1,   1.5, 1.5 ],
          [ 1.5, 1,   1.5 ],
          [ 1.5, 1.5, 1   ],
          [ 0.5, 0,   0   ],  # O  0, 0, 0
          [ 0,   0.5, 0   ],
          [ 0,   0,   0.5 ],
          [ 0.5, 0.5, 0.5 ],
          [ 1.5, 0,   0   ],  # O  1, 0, 0
          [ 1,   0.5, 0   ],
          [ 1,   0,   0.5 ],
          [ 1.5, 0.5, 0.5 ],
          [ 0.5, 1,   0   ],  # O  0, 1, 0
          [ 0,   1.5, 0   ],
          [ 0,   1,   0.5 ],
          [ 0.5, 1.5, 0.5 ],
          [ 0.5, 0,   1   ],  # O  0, 0, 1
          [ 0,   0.5, 1   ],
          [ 0,   0,   1.5 ],
          [ 0.5, 0.5, 1.5 ],
          [ 1.5, 1,   0   ],  # O  1, 1, 0
          [ 1,   1.5, 0   ],
          [ 1,   1,   0.5 ],
          [ 1.5, 1.5, 0.5 ],
          [ 1.5, 0,   1   ],  # O  1, 0, 1
          [ 1,   0.5, 1   ],
          [ 1,   0,   1.5 ],
          [ 1.5, 0.5, 1.5 ],
          [ 0.5, 1,   1   ],  # O  0, 1, 1
          [ 0,   1.5, 1   ],
          [ 0,   1,   1.5 ],
          [ 0.5, 1.5, 1.5 ],
          [ 1.5, 1,   1   ],  # O  1, 1, 1
          [ 1,   1.5, 1   ],
          [ 1,   1,   1.5 ],
          [ 1.5, 1.5, 1.5 ] ],
        [28,] * 32 + [8,] * 32
      ]

# check whether there is dupulication
for i, p1 in enumerate( nio_sc[1] ):
    for p2 in nio_sc[1][ i+1 : ]:
        if p1 == p2:
            print( p1 )

# lattice constants should be twice in 2*2*2 super cell
nio_sc[0] = [ list( a_ ) for a_ in np.array( nio_sc[0] ) * 2 ]

# coordinate should be divided by 2 
nio_sc[1] = [ list( p_ ) for p_ in np.array( nio_sc[1] ) / 2 ]

nio_sc = tuple( nio_sc )

print( spglib.get_spacegroup( nio_sc ) )
print( spglib.find_primitive( nio_sc ) )

座標の重複があると、'None'が返って来るため、チェックする必要がある。目で見て一生懸命探すのには限界がある。

最後に反強磁性だが、格子定数と座標はsuper cellのままでよく、原子番号で反強磁性秩序の情報を入れる。仮に磁化軸方向と平行な原子をNi(28)として、反平行なものを仮にPd(46)とする。酸素はそのままで良い。
NiOの反強磁性磁化軸は[111]方向なので、それを考慮して原子番号を対応させていく。

nio_af \
    = [ [ ( a, 0, 0 ), ( 0, a, 0 ), ( 0, 0, a ) ],
        [ [ 0,   0,   0   ],  # Ni 0, 0, 0
          [ 0,   0.5, 0.5 ],
          [ 0.5, 0,   0.5 ],
          [ 0.5, 0.5, 0   ],
          [ 1,   0,   0   ],  # Ni 1, 0, 0
          [ 1,   0.5, 0.5 ],
          [ 1.5, 0,   0.5 ],
          [ 1.5, 0.5, 0   ],
          [ 0,   1,   0   ],  # Ni 0, 1, 0
          [ 0,   1.5, 0.5 ],
          [ 0.5, 1,   0.5 ],
          [ 0.5, 1.5, 0   ],
          [ 0,   0,   1   ],  # Ni 0, 0, 1
          [ 0,   0.5, 1.5 ],
          [ 0.5, 0,   1.5 ],
          [ 0.5, 0.5, 1   ],
          [ 1,   1,   0   ],  # Ni 1, 1, 0
          [ 1,   1.5, 0.5 ],
          [ 1.5, 1,   0.5 ],
          [ 1.5, 1.5, 0   ],
          [ 1,   0,   1   ],  # Ni 1, 0, 1
          [ 1,   0.5, 1.5 ],
          [ 1.5, 0,   1.5 ],
          [ 1.5, 0.5, 1   ],
          [ 0,   1,   1   ],  # Ni 0, 1, 1
          [ 0,   1.5, 1.5 ],
          [ 0.5, 1,   1.5 ],
          [ 0.5, 1.5, 1   ],
          [ 1,   1,   1   ],  # Ni 1, 1, 1
          [ 1,   1.5, 1.5 ],
          [ 1.5, 1,   1.5 ],
          [ 1.5, 1.5, 1   ],
          [ 0.5, 0,   0   ],  # O  0, 0, 0
          [ 0,   0.5, 0   ],
          [ 0,   0,   0.5 ],
          [ 0.5, 0.5, 0.5 ],
          [ 1.5, 0,   0   ],  # O  1, 0, 0
          [ 1,   0.5, 0   ],
          [ 1,   0,   0.5 ],
          [ 1.5, 0.5, 0.5 ],
          [ 0.5, 1,   0   ],  # O  0, 1, 0
          [ 0,   1.5, 0   ],
          [ 0,   1,   0.5 ],
          [ 0.5, 1.5, 0.5 ],
          [ 0.5, 0,   1   ],  # O  0, 0, 1
          [ 0,   0.5, 1   ],
          [ 0,   0,   1.5 ],
          [ 0.5, 0.5, 1.5 ],
          [ 1.5, 1,   0   ],  # O  1, 1, 0
          [ 1,   1.5, 0   ],
          [ 1,   1,   0.5 ],
          [ 1.5, 1.5, 0.5 ],
          [ 1.5, 0,   1   ],  # O  1, 0, 1
          [ 1,   0.5, 1   ],
          [ 1,   0,   1.5 ],
          [ 1.5, 0.5, 1.5 ],
          [ 0.5, 1,   1   ],  # O  0, 1, 1
          [ 0,   1.5, 1   ],
          [ 0,   1,   1.5 ],
          [ 0.5, 1.5, 1.5 ],
          [ 1.5, 1,   1   ],  # O  1, 1, 1
          [ 1,   1.5, 1   ],
          [ 1,   1,   1.5 ],
          [ 1.5, 1.5, 1.5 ] ],
        [ 28, 46, 46, 46 ]
      + [ 46, 28, 28, 28 ] * 3
      + [ 28, 46, 46, 46 ] * 3
      + [ 46, 28, 28, 28 ]
      + [ 8, ] * 32
      ]

# check whether there is dupulication
for i, p1 in enumerate( nio_af[1] ):
    for p2 in nio_af[1][ i+1 : ]:
        if p1 == p2:
            print( p1 )

# lattice constants should be twice in 2*2*2 super cell
nio_af[0] = [ list( a_ ) for a_ in np.array( nio_af[0] ) * 2 ]

# coordinate should be divided by 2 
nio_af[1] = [ list( p_ ) for p_ in np.array( nio_af[1] ) / 2 ]

nio_af = tuple( nio_af )

print( spglib.get_spacegroup( nio_af ) )
#R-3m (166)

print( spglib.find_primitive( nio_af ) )
#(array([[ 1.48138871,  0.85528017,  4.83819526],
#       [-1.48138871,  0.85528017,  4.83819526],
#       [ 0.        , -1.71056034,  4.83819526]]), array([[-3.70074342e-17,  0.00000000e+00,  0.00000000e+00],
#       [ 5.00000000e-01,  5.00000000e-01,  5.00000000e-01],
#       [ 7.50000000e-01,  7.50000000e-01,  7.50000000e-01],
#       [ 2.50000000e-01,  2.50000000e-01,  2.50000000e-01]]), array([28, 46,  8,  8], dtype=int32))

対称性がFm-3mからR-3mに落ちているのが分かる。

ここで、結晶軸がベクトル表記のままだと結晶構造作成プログラムでは不便な場合があるので、格子定数の形式に変換すると、

prim_nio_af = spglib.find_primitive( nio_af )
axes = prim_nio_af[0]

abc = [ np.linalg.norm( axis ) for axis in axes ]
print( abc )
#[5.131681011130759, 5.131681011130759, 5.131681011130759]

# alpha beta gamma
abg =[]
abg.append( np.degrees( np.arccos( np.dot( axes[1], axes[2] ) / np.linalg.norm(axes[1]) / np.linalg.norm(axes[2]) ) ) )
abg.append( np.degrees( np.arccos( np.dot( axes[2], axes[0] ) / np.linalg.norm(axes[2]) / np.linalg.norm(axes[0]) ) ) )
abg.append( np.degrees( np.arccos( np.dot( axes[0], axes[1] ) / np.linalg.norm(axes[0]) / np.linalg.norm(axes[1]) ) ) )
print( abg )
#[33.55730976192072, 33.55730976192072, 33.55730976192072]

これで、第一原理計算プログラムで反強磁性が計算できるようになった。

2020-02-09

gfortran での変数の次元の制約。

プログラミング

何も考えずにプログラムを作っていたら、配列の次元が7個を超えるとコンパイルできないらしい。

program test_array7

real( kind( 0d0 ) ) :: r( 1, 2, 3, 4, 5, 6, 7, 8 )

end program test_array7

これを gfortran でコンパイルすると

test_array7.f90:3:46:

 real( kind( 0d0 ) ) :: r( 1, 2, 3, 4, 5, 6, 7, 8 )
                                              1
Error: Array specification at (1) has more than 7 dimensions

添字が11個ある変数の計算しようとしていたので、これはちょっと困った。。。
構造体とかを使って、何とか逃げるしかないか。。。
ちなみに構造体なら7個以上の要素があっても大丈夫であった。

program test_type7                                                              
 
implicit none
 
type :: test
    real( kind( 0d0 ) ) :: r1
    real( kind( 0d0 ) ) :: r2
    real( kind( 0d0 ) ) :: r3
    real( kind( 0d0 ) ) :: r4
    real( kind( 0d0 ) ) :: r5
    real( kind( 0d0 ) ) :: r6
    real( kind( 0d0 ) ) :: r7
    real( kind( 0d0 ) ) :: r8
end type
 
type( test ) :: t
 
end program test_type7

2020-01-26

Fourier変換の符号や係数の気持ち。

数学量子力学

自分で毎回忘れるので、まとめておく。

Fourier変換の符号や係数は、変換と逆変換とで対になっていれば良いため（語弊があるが）不定性があり、分野によって異なる（気がしている）。

一応、自分ルールとして、

位置 $x$ および時間 $t$ を基準とし、波数 $k$ もしくは（角）振動数 $\omega$ に変換にする方向を「順」方向。
したがって、 $k \rightarrow x$ または $\omega \rightarrow t$ は「逆」方向。

と考える。

量子力学的な観点でここでは考えていく。
運動量およびエネルギーの固有関数は、それぞれ次のような平面波で表される。
$\displaystyle \phi( x; p ) = \langle x | p \rangle = N e^{ i k x } \quad ( p = \hbar k ) \\ \displaystyle \left( \hat p( x ) \phi( x; p ) = -i\hbar \frac{ d }{ dx } \phi( x; p ) = p \phi( x; p ) \right)$
$\displaystyle \psi( x, t; E ) = \psi( x; E ) e^{ - i \omega t } \quad ( E = \hbar \omega ) \\ \displaystyle \left( i\hbar \frac{ d }{ dt } \psi( x, t; E ) = E \psi( x, t; E) \right)$

したがって、任意の $x$ や $t$ の関数をこれらの固有関数の重なりで表現すると思うと、
$\displaystyle f( x ) = c_x \int dk \, f(k) e^{ i k x } \\ \displaystyle f( t ) = c_t \int d\omega \, f(\omega) e^{ - i \omega t }$
という位相を採用するのが適当であろう。
注意するべきは、これらは「逆変換」の式である。何となく「逆」とついていると後から決まるような感じがするが、ここではそうではない。

これによって、「順」方向側の変換の位相が定まる。
$\displaystyle f( k ) = c_k \int dx \, f(x) e^{ - i k x } \\ \displaystyle f( \omega ) = c_\omega \int dt \, f(t) e^{ i \omega t }$
「基準となる関数から、知りたい平面波の係数を抜き出す」というイメージになるだろうか。

残るは係数であるが、
$\displaystyle c_x c_k = \frac{1}{2\pi} \\ \displaystyle c_t c_\omega = \frac{1}{2\pi}$
という関係であるため、一意には決まらない。

ここで、Fourier変換のモチベーションに戻ると、基準である $x, t$ の関数はわかっているから、知りたいのは $f(k), f(\omega)$ である。
そのため、 $f(k) = \cdots$ , $f(\omega) = \cdots$ という式をパッと使いたくなるし、頻出することが予想される。それゆえの「順」方向である。
したがって、順方向側の係数を簡単にして、逆方向側にメンドクサイのを押し付ける方が便利（なように感じるの）である。
よって、
$\displaystyle c_x = \frac{1}{2\pi}, \quad c_k = 1 \\ \displaystyle c_t = \frac{1}{2\pi}, \quad c_\omega = 1$

最終形は、
$\displaystyle f( k ) = \int dx \, f(x) e^{ - i k x } \Leftrightarrow f( x ) = \frac{1}{2\pi} \int dk \, f(k) e^{ i k x } \\ \displaystyle f( \omega ) = \int dt \, f(t) e^{ i \omega t } \Leftrightarrow f( t ) = \frac{1}{2\pi} \int d\omega \, f(\omega) e^{ - i \omega t }$

2020-01-20

実空間多重散乱からKKR法への接続。

固体物理量子力学数学

あるセル内の散乱によるセルT行列を $\mathcal T$ と置くと、全てのサイトによる散乱Green関数 $G$ は次のように書ける。
$\displaystyle G^{\alpha_i \beta_j}_{LL'} = \left[ \left( g_0^{-1} - \mathcal T \right)^{-1} \right]^{\alpha_i \beta_j}_{LL'} = g^{\alpha_i \beta_j}_{0,LL'} + \sum_{m} \sum_{\gamma_1, \gamma_2} \sum_{L_1, L_2} g^{\alpha_i \gamma_{1m}}_{0,L L_1} \mathcal T^{(m) \gamma_1 \gamma_2}_{L_1L_2} G^{\gamma_{2m} \beta_j}_{L_2 L'}$
ここで、 $g_0$ は自由散乱状態におけるGreen関数である。また、セル $i$ に属するサイト $\alpha$ を $\alpha_i$ という風に表した。

セルT行列 $\mathcal T$ はサイトT行列 $t$ を用いて次のように書ける。
$\displaystyle g^{(i)\alpha \beta}_{LL'} = \left[ \left( g_0^{-1} - t \right)^{-1} \right]^{\alpha_i \beta_i}_{LL'} \\ \displaystyle \mathcal T^{\alpha_i \beta_j}_{LL'} = \mathcal T^{(i) \alpha \beta}_{LL'} \delta_{ij} \\ \displaystyle \mathcal T^{(i) \alpha \beta}_{LL'} = t^{\alpha_i}_L \delta_{\alpha_i \beta_i} \delta_{LL'} + t^{\alpha_i}_L g^{(i) \alpha \beta}_{LL'} t^{\beta_i}_{L'} = \left[ t \left( I - g_0 t \right)^{-1} \right]^{\alpha_i \beta_i}_{LL'}$

セルT行列に関する詳しい話は前回まとめた。
koideforest.hatenadiary.com

ここで、系が並進対称性を持ち、ユニットセルが周期的に繰り返しているとすると、セルT行列はセルに依存しなくなる（ $\mathcal T^{(i)} = \mathcal T$ , $t^{\alpha_i} = t^\alpha$ ）。

Fourier級数展開によって、セルの添字を結晶運動量に変換する。
$\displaystyle G^{\alpha_i \beta_j}_{LL'} = \frac{1}{N^2} \sum_{k,k'} e^{ i k \cdot R_i } G^{\alpha \beta}_{LL'}( k, k') e^{ - i k' \cdot R_j } \\ \displaystyle G^{\alpha \beta}_{LL'}( k, k') = \sum_{i, j} e^{ - i k \cdot R_i } G^{\alpha_i \beta_j}_{LL'} e^{ i k' \cdot R_j }$

$g_0$ はサイト間の差にしか依存しないため、
$\displaystyle g^{\alpha_i \beta_j}_{0,LL'} = g_{0, LL'}( ( R_\alpha + R_i ) - ( R_\beta + R_j ) ) \\ \displaystyle \qquad = g_{0, LL'}( ( R_\alpha - R_\beta ) + ( R_i - R_j ) ) = \frac{1}{N} \sum_{k} e^{ i k \cdot ( R_i - R_j ) } g^{\alpha \beta}_{0, LL'}( k ) \\ \displaystyle g^{\alpha \beta}_{0,LL'}( k ) = \sum_{i} e^{ - i k \cdot R_i } g^{\alpha_i \beta_0}_{0,LL'}$

したがって、
$\displaystyle \frac{1}{N^2} \sum_{k,k'} e^{ i k \cdot R_i } G^{\alpha \beta}_{LL'}( k, k') e^{ - i k' \cdot R_j } \\ \displaystyle \qquad = \frac{1}{N} \sum_{k} e^{ i k \cdot ( R_i - R_j ) } g^{\alpha \beta}_{0, LL'}( k ) \\ \displaystyle \qquad \qquad +\frac{1}{N^3} \sum_{ k, k', k''}\sum_{m} \sum_{\gamma_1, \gamma_2} \sum_{L_1, L_2} e^{ i k \cdot ( R_i - R_m ) } e^{ i k' \cdot R_m } e^{ - i k'' \cdot R_j } \\ \displaystyle \qquad \qquad \qquad \times g^{\alpha \gamma_{1}}_{0,L L_1}( k ) \mathcal T^{ \gamma_1 \gamma_2}_{L_1L_2} G^{\gamma_{2} \beta}_{L_2 L'}( k', k'') \\ \displaystyle \qquad = \frac{1}{N} \sum_{k} e^{ i k \cdot ( R_i - R_j ) } g^{\alpha \beta}_{0, LL'}( k ) \\ \displaystyle \qquad \qquad +\frac{1}{N^2} \sum_{ k, k', k''} \sum_{\gamma_1, \gamma_2} \sum_{L_1, L_2} e^{ i k \cdot R_i } e^{ - i k'' \cdot R_j } \delta_{ k, k' } \\ \displaystyle \qquad \qquad \qquad \times g^{\alpha \gamma_{1}}_{0,L L_1}( k ) \mathcal T^{ \gamma_1 \gamma_2}_{L_1L_2} G^{\gamma_{2} \beta}_{L_2 L'}( k', k'') \\ \displaystyle \qquad = \frac{1}{N} \sum_{k} e^{ i k \cdot ( R_i - R_j ) } g^{\alpha \beta}_{0, LL'}( k ) \\ \displaystyle \qquad \qquad +\frac{1}{N^2} \sum_{ k, k''} \sum_{\gamma_1, \gamma_2} \sum_{L_1, L_2} e^{ i k \cdot R_i } e^{ - i k'' \cdot R_j } g^{\alpha \gamma_{1}}_{0,L L_1}( k ) \mathcal T^{ \gamma_1 \gamma_2}_{L_1L_2} G^{\gamma_{2} \beta}_{L_2 L'}( k, k'') \\ \displaystyle \qquad = \frac{1}{N} \sum_{k,k'} e^{ i k \cdot R_i } g^{\alpha \beta}_{0, LL'}( k ) \delta_{kk'} e^{ - i k' \cdot R_j } \\ \displaystyle \qquad \qquad +\frac{1}{N^2} \sum_{ k, k'} \sum_{\gamma_1, \gamma_2} \sum_{L_1, L_2} e^{ i k \cdot R_i } g^{\alpha \gamma_{1}}_{0,L L_1}( k ) \mathcal T^{ \gamma_1 \gamma_2}_{L_1L_2} G^{\gamma_{2} \beta}_{L_2 L'}( k, k') e^{ - i k' \cdot R_j }$

よって
$\displaystyle G^{\alpha \beta}_{LL'}( k, k' ) = N g^{\alpha \beta}_{0, LL'}( k ) \delta_{kk'} + \sum_{\gamma_1, \gamma_2} \sum_{L_1, L_2} g^{\alpha \gamma_{1}}_{0,L L_1}( k ) \mathcal T^{ \gamma_1 \gamma_2}_{L_1L_2} G^{\gamma_{2} \beta}_{L_2 L'}( k, k') \\ \displaystyle \therefore G^{\alpha \beta}_{LL'}( k, k' ) = \left[ ( I - g_0( k ) \mathcal T )^{-1} N g( k ) \delta_{kk'} \right]^{\alpha \beta}_{LL'} = N \delta_{kk'} \left[ ( g_0^{-1}( k ) - \mathcal T )^{-1} \right]^{\alpha \beta}_{LL'}$

逆行列計算に計算コストが掛かるため、ユニットセル内の原子だけに抑えられるのはかなりのメリットだと思われる。

2020-01-20

T行列のサイトT行列展開とセルT行列展開。

固体物理数学

T行列は元々サイトポテンシャル $v$ の展開で以下の様に定義されている。
$\displaystyle T = \sum_\alpha v_\alpha + \sum_{\alpha, \beta} v_\alpha g_0 v_\beta + \sum_{\alpha, \beta, \gamma} v_\alpha g_0 v_\beta g_0 v_\gamma + \cdots$

しかし、「違うサイトへの伝搬」と「同じサイト内での多重散乱」が混在しており、実際の計算には非常に不便である。
そこで、同じサイトのポテンシャル散乱を繰り込んだサイトT行列 $t$ を定義すると、サイト内の多重散乱とサイト間の伝搬を分けることができる。これがいわゆる多重散乱展開である。
$\displaystyle t_\alpha = \sum_\alpha v_\alpha + \sum_{\alpha} v_\alpha g_0 v_\alpha + \sum_\alpha v_\alpha g_0 v_\alpha g_0 v_\alpha \cdots \\ \displaystyle T = \sum_\alpha t_\alpha + \sum_{\alpha} \sum_{ \beta (\neq \alpha ) } t_\alpha g_0 t_\beta + \sum_{\alpha} \sum_{ \beta (\neq \alpha ) } \sum_{ \gamma (\neq \beta ) } t_\alpha g_0 t_\beta g_0 t_\gamma + \cdots$
サイトT行列間の伝搬の際には、かならず異なるサイトへ行くため、散乱のイメージが付きやすい。
注意として、あるサイトから他サイトに伝搬した後、もう一度同じサイトに戻って来ることは可能である。

この考えを更に進めて、サイトを各集合に分類し、それをセルと呼ぶことにすると、セル内の多重散乱とセル間の散乱に分けることができる。セル内の多重散乱を繰り込んだものをセルT行列 $\mathcal T$ と定義すると、
$\displaystyle \mathcal T_i = \sum_{a_i} t_{a_i} + \sum_{a_i} \sum_{ b_i (\neq a_i ) } t_{a_i} g_0 t_{b_i} + \sum_{a_i} \sum_{ b_i (\neq a_i) } \sum_{ c_i (\neq b_i) } t_{a_i} g_0 t_{b_i} g_0 t_{c_i} + \cdots \\ \displaystyle T = \sum_i \mathcal T_i + \sum_i \sum_{j (\neq i)} \mathcal T_i g_0 \mathcal T_j + \sum_{i} \sum_{j (\neq i) } \sum_{ k (\neq j) } \mathcal T_i g_0 \mathcal T_j g_0 \mathcal T_k + \cdots$
ここで、セル $i$ に含まれてるサイト $a_i$ は、和を取ると全サイトの和に等しくなるように定義している。
$\displaystyle \sum_\alpha = \sum_i \sum_{a_i}$
注意として、セル間の散乱といえども、伝搬はあくまでサイト間である。「サイト同士が必ず異なるセルに属するような散乱」を「セル間の散乱」とここでは呼んでいる。

まとめると、
$\displaystyle T = \sum_\alpha v_\alpha + \sum_{\alpha, \beta} v_\alpha g_0 v_\beta + \sum_{\alpha, \beta, \gamma} v_\alpha g_0 v_\beta g_0 v_\gamma + \cdots \\ \displaystyle \qquad = \sum_\alpha t_\alpha + \sum_{\alpha} \sum_{ \beta (\neq \alpha ) } t_\alpha g_0 t_\beta + \sum_{\alpha} \sum_{ \beta (\neq \alpha ) } \sum_{ \gamma (\neq \beta ) } t_\alpha g_0 t_\beta g_0 t_\gamma + \cdots \\ \displaystyle \qquad = \sum_i \mathcal T_i + \sum_i \sum_{j (\neq i)} \mathcal T_i g_0 \mathcal T_j + \sum_{i} \sum_{j (\neq i) } \sum_{ k (\neq j) } \mathcal T_i g_0 \mathcal T_j g_0 \mathcal T_k + \cdots$

以下、もう少し細かく説明する。

サイトT行列展開

サイトの和から、あるサイト $a$ だけを抜き出して $t_a$ を作ることを考える。
$\displaystyle \sum_\alpha v_\alpha = v_a + \sum_{\alpha (\neq a)} v_\alpha \\ \displaystyle \sum_{\alpha, \beta} v_\alpha g_0 v_\beta = v_a g_0 v_a + \sum_{\alpha (\neq a)} ( v_a g_0 v_\alpha + v_\alpha g_0 v_a ) + \sum_{\alpha (\neq a)} \sum_{\beta (\neq a)} v_\alpha g_0 v_\beta \\ \displaystyle \sum_{\alpha, \beta, \gamma} v_\alpha g_0 v_\beta g_0 v_\gamma = v_a g_0 v_a g_0 v_a + \sum_{\alpha (\neq a)} ( v_a g_0 v_a g_0 v_\alpha + v_a g_0 v_\alpha g_0 v_a + v_\alpha g_0 v_a g_0 v_a ) \\ \displaystyle \qquad + \sum_{\alpha (\neq a)} \sum_{\beta (\neq a)} ( v_a g_0 v_\alpha g_0 v_\beta + v_\beta g_0 v_a g_0 v_\alpha + v_\alpha g_0 v_\beta g_0 v_a ) \\ \displaystyle \qquad + \sum_{\alpha (\neq a)} \sum_{\beta (\neq a)} \sum_{\gamma (\neq a)} v_\alpha g_0 v_\beta g_0 v_\gamma$

したがって、三次までで整理すると、
$\displaystyle T \approx \sum_\alpha v_\alpha + \sum_{\alpha, \beta} v_\alpha g_0 v_\beta + \sum_{\alpha, \beta, \gamma} v_\alpha g_0 v_\beta g_0 v_\gamma \\ \displaystyle \qquad = (v_a + v_a g_0 v_a + v_a g_0 v_a g_0 v_a ) \\ \displaystyle \qquad \qquad + \sum_{\alpha (\neq a)} ( (v_a + v_a g_0 v_a ) g_0 v_\alpha + v_\alpha g_0 (v_a + v_a g_0 v_a ) ) + \sum_{\alpha (\neq a)} v_a g_0 v_\alpha g_0 v_a \\ \displaystyle \qquad \qquad + \sum_{\alpha (\neq a)} \sum_{\beta (\neq a)} ( v_a g_0 v_\alpha g_0 v_\beta + v_\beta g_0 v_a g_0 v_\alpha + v_\alpha g_0 v_\beta g_0 v_a ) \\ \displaystyle \qquad \qquad + \sum_{\alpha (\neq a)} v_\alpha + \sum_{\alpha (\neq a)} \sum_{\beta (\neq a)} v_\alpha g_0 v_\beta + \sum_{\alpha (\neq a)} \sum_{\beta (\neq a)} \sum_{\gamma (\neq a)} v_\alpha g_0 v_\beta g_0 v_\gamma \\ \displaystyle \qquad \rightarrow t_a + \sum_{\alpha (\neq a)} ( t_a g_0 v_\alpha + v_\alpha g_0 t_a ) + \sum_{\alpha (\neq a)} t_a g_0 v_\alpha g_0 t_a \\ \displaystyle \qquad \qquad + \sum_{\alpha (\neq a)} \sum_{\beta (\neq a)} ( t_a g_0 v_\alpha g_0 v_\beta + v_\beta g_0 t_a g_0 v_\alpha + v_\alpha g_0 v_\beta g_0 t_a ) \\ \displaystyle \qquad \qquad + \sum_{\alpha (\neq a)} v_\alpha + \sum_{\alpha (\neq a)} \sum_{\beta (\neq a)} v_\alpha g_0 v_\beta + \sum_{\alpha (\neq a)} \sum_{\beta (\neq a)} \sum_{\gamma (\neq a)} v_\alpha g_0 v_\beta g_0 v_\gamma$
級数の中で $v_a$ がやがて $t_a$ になっていくことが何となく感じられたと思う。
ポイントは、 $a$ から移るときには必ず異なるサイトへ行くようになっていることである。
他のサイトに関しても同様に $t$ を作るように抜き出していくと、最終的にサイトT行列展開が得られる。

セル行列展開

サイトT行列展開の時と同様に、サイトの和からあるセル $i$ に含まれるサイト集合 $\{a_i\}$ を抜き出して $\mathcal T_i$ を作ることを考える。
$\displaystyle \sum_{\alpha} t_{\alpha} = \sum_{a_i} t_{a_i} + \sum_{ j (\neq i )} \sum_{a_j} t_{a_j} \\ \displaystyle \sum_{\alpha} \sum_{\beta (\neq \alpha) } t_\alpha g_0 t_\beta = \sum_{a_i} \sum_{ b_i (\neq a_i ) } t_{a_i} g_0 t_{b_i} \\ \displaystyle \qquad + \sum_{a_i} \left( \sum_{ j (\neq i) }\sum_{a_j} \right) ( t_{a_i} g_0 t_{a_j} + t_{a_j} g_0 t_{a_i} ) \\ \displaystyle \qquad + \left( \sum_{j (\neq i)} \sum_{a_j} \right) \left( \sum_{ k (\neq i) } \sum_{b_k (\neq a_j)} \right) t_{a_j} g_0 t_{b_k} \\ \displaystyle \sum_{\alpha} \sum_{\beta (\neq \alpha)} \sum_{\gamma (\neq \beta)} t_\alpha g_0 t_\beta g_0 t_\gamma = \sum_{a_i} \sum_{b_i (\neq a_i)} \sum_{c_i (\neq b_i)} t_{a_i} g_0 t_{b_i} g_0 t_{c_i} \\ \displaystyle \qquad + \sum_{a_i} \sum_{b_i (\neq a_i )} \left( \sum_{ j (\neq i)} \sum_{a_j} \right) ( t_{a_i} g_0 t_{b_i} g_0 t_{a_j} + t_{b_i} g_0 t_{a_j} g_0 t_{a_i} + t_{a_j} g_0 t_{a_i} g_0 t_{b_i} ) \\ \displaystyle \qquad + \sum_{a_i} \left( \sum_{j (\neq i)} \sum_{a_j} \right) \left( \sum_{k (\neq i)} \sum_{b_k (\neq a_j)} \right) ( t_{a_i} g_0 t_{a_j} g_0 t_{b_k} + t_{b_k} g_0 t_{a_i} g_0 t_{a_j} + t_{a_j} g_0 t_{b_k} g_0 t_{a_i} ) \\ \displaystyle \qquad + \left( \sum_{ j (\neq i) } \sum_{a_j} \right) \left( \sum_{ k (\neq i) } \sum_{ b_k (\neq a_j) } \right) \left( \sum_{ m (\neq i) } \sum_{c_m (\neq b_k)} \right) t_{a_j} g_0 t_{b_k} g_0 t_{c_m}$
異なるセルに属していれば、サイトが被ることは無い。

したがって、三次までで整理すると、
$\displaystyle T \approx \sum_\alpha t_\alpha + \sum_{\alpha} \sum_{\beta (\neq \alpha)} t_\alpha g_0 t_\beta + \sum_{\alpha} \sum_{\beta (\neq \alpha)} \sum_{\gamma (\neq \beta)} t_\alpha g_0 t_\beta g_0 t_\gamma \\ \displaystyle \qquad \rightarrow \mathcal T_i + \left( \sum_{ j (\neq i) } \sum_{a_j} \right) ( \mathcal T_i g_0 t_{a_j} + t_{a_j} g_0 \mathcal T_i ) + \left( \sum_{ j (\neq i ) } \sum_{a_j} \right) \mathcal T_i g_0 t_{a_j} g_0 \mathcal T_i$
$\displaystyle \qquad \qquad + \left( \sum_{ j (\neq i ) } \sum_{ a_j } \right) \left( \sum_{k (\neq i)} \sum_{ b_k ( \neq a_j ) } \right) ( \mathcal T_i g_0 t_{a_j} g_0 t_{b_k} + t_{b_k} g_0 \mathcal T_i g_0 t_{a_j} + t_{a_j} g_0 t_{b_k} g_0 \mathcal T_i ) \\ \displaystyle \qquad \qquad + \sum_{ j (\neq i )} \sum_{a_j} t_{a_j} + \left( \sum_{j (\neq i)} \sum_{a_j} \right) \left( \sum_{ k (\neq i) } \sum_{b_k (\neq a_j)} \right) t_{a_j} g_0 t_{b_k} \\ \displaystyle \qquad \qquad + \left( \sum_{ j (\neq i) } \sum_{a_j} \right) \left( \sum_{ k (\neq i) } \sum_{ b_k (\neq a_j) } \right) \left( \sum_{ m (\neq i) } \sum_{c_m (\neq b_k)} \right) t_{a_j} g_0 t_{b_k} g_0 t_{c_m}$
級数の中で $\sum_{a_i} t_{a_i}$ がやがて $\mathcal T_i$ になっていくことが何となく感じられたと思う。
セル $i$ のみを抜き出したため、 $j, k$ の和は $j=k$ となる場合があることに注意。そのために同じサイトへの伝搬が無い様に和が複雑になっている。
他のセルに関しても同様に $\mathcal T$ を作るように抜き出していくと、最終的にセルT行列展開が得られる。

2020-01-16

python で Lebedev quadrature を使う。

数学プログラミング

注意：情報を更新しました。
koideforest.hatenadiary.com

注意：以下の内容は古いので、テストコードはそのままでは通りません。

以前にLebedev求積法について紹介した。
koideforest.hatenadiary.com

最近ではpythonでもLebedev求積法が使えるmoduleがあるらしい。
quadpy · PyPI

とりあえずテストコードを書いた。

import numpy as np
import quadpy

def f_cart( xyz ):
    # np.shape( xyz ) -> ( 3, number of points )
    # n_points = np.shape( xyz )[1]
    return np.ones( n_points )
    # DO NOT WRITE "return 1"

def f_sphere( phi, theta ):
    # phi: azimuth
    # theta: polar
    # np.shape( phi ) = np.shape( theta ) = ( number of points )
    n_points = np.size( theta )
    return np.ones( n_points )
    # DO NOT WRITE "return 1"

scheme = quadpy.sphere.lebedev_019()

# \int^\pi_0 d\theta \sin\theta \int^{2\pi}_0 d\phi = 4 \pi

xyz_center = ( 0., 0., 0. )
radius = 1.
val_cart = scheme.integrate( f_cart , xyz_center, radius )

val_sphere  = scheme.integrate_spherical( f_sphere )

print( val_cart, val_sphere, 4*np.pi )

以前に紹介したFortranのコードとは異なり、こちらは内部で $4\pi$ を掛けてくれるため、被積分関数にだけ気を付ければ良い。

注意するべきは、被積分関数の仕様である。
quadpyの中身の一部を抜粋すると以下のようになる。

class SphereScheme:
    def __init__(self, name, weights, points, azimuthal_polar, degree, citation):
        self.name = name
        self.degree = degree
        self.citation = citation

        if weights.dtype == numpy.float64:
            self.weights = weights
        else:
            assert weights.dtype in [numpy.dtype("O"), numpy.int64]
            self.weights = weights.astype(numpy.float64)
            self.weights_symbolic = weights

        if points.dtype == numpy.float64:
            self.points = points
        else:
            assert points.dtype in [numpy.dtype("O"), numpy.int64]
            self.points = points.astype(numpy.float64)
            self.points_symbolic = points

        if azimuthal_polar.dtype == numpy.float64:
            self.azimuthal_polar = azimuthal_polar
        else:
            assert azimuthal_polar.dtype in [numpy.dtype("O"), numpy.int64]
            self.azimuthal_polar = azimuthal_polar.astype(numpy.float64)
            self.azimuthal_polar_symbolic = azimuthal_polar
        return

    def integrate(self, f, center, radius, dot=numpy.dot):
        """Quadrature where `f` is defined in Cartesian coordinates.
        """
        center = numpy.array(center)
        rr = numpy.multiply.outer(radius, self.points)
        rr = numpy.swapaxes(rr, 0, -2)
        ff = numpy.array(f((rr + center).T))
        return area(radius) * dot(ff, self.weights)

    def integrate_spherical(self, f, dot=numpy.dot):
        """Quadrature where `f` is a function of the spherical coordinates
        `azimuthal` and `polar` (in this order).
        """
        flt = numpy.vectorize(float)
        rr = numpy.swapaxes(flt(self.azimuthal_polar), 0, -2)
        ff = numpy.array(f(rr.T[0], rr.T[1]))
        return area(1.0) * dot(ff, flt(self.weights))

def area(radius):
    return 4 * numpy.pi * numpy.array(radius) ** 2

重要なのは "dot(ff, self.weights)" であり、和を内積で処理しているため、"ff" は "numpy.ndarray" の型でないといけない。
つまり、定義した被積分関数はスカラーでなく配列を返すように工夫する必要がある。
上で挙げたテストコードで、被積分関数が "1" だからといって単純に "1" のみを返すようにすると、"dot(1, self.weights) -> self.weights" となるため、積分値（スカラー）ではなく（各点の重みを表す）配列が返って来てしまう。
また、直交座標を用いる "integrate" では、"ff" に渡す座標の配列が "self.points" に対して転置を取っているため、配列の順番に気を付ける必要がある。

余談：
”numpy.swapaxes(rr, 0, -2)” で何をやっているのか不明だが、とりあえず気にしなくても動くっぽい。
swapaxes
numpy.swapaxes — NumPy v1.15 Manual

nano_exit

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