nano_exit

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

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

定数ポテンシャルを与える球対称電子密度について。

数学 固体物理 電磁気

Poisson方程式の解は、
\displaystyle V( \vec r ) = \int d\vec r' \, \frac{ \rho ( \vec r' ) }{ | \vec r - \vec r' | }
と書ける（原子単位系）。

多重極子展開および角運動量展開を利用すると、
\displaystyle V( \vec r ) = \sum_L \frac{ 4\pi }{ 2 l + 1 } \int d\vec r' \, \frac{ r_<^l }{ r_>^{l+1} } Y_L( \hat r ) Y_L^*( \hat r' ) \rho ( \vec r' ) \\ \displaystyle \qquad = \sum_L \frac{ 4\pi }{ 2 l + 1 } Y_L( \hat r ) \int d\vec r' \, \frac{ r_<^l }{ r_>^{l+1} } Y_L^*( \hat r' ) \sum_{L'} \rho_{l'} ( r' ) Y_{L'}( \hat r' ) \\ \displaystyle \qquad = \sum_L \frac{ 4\pi }{ 2 l + 1 } Y_L( \hat r ) \int^\infty_0 dr' \, r'^2 \frac{ r_<^l }{ r_>^{l+1} } \, \rho_{l} ( r' )

ここで、元々 \rho( \vec r ) = \rho( r ) の時、
\displaystyle \rho( r ) = \rho_0( r ) Y_{00} = \frac{ \rho_0( r ) }{ \sqrt{ 4 \pi } } \\ \displaystyle \therefore \rho_0( r ) = \sqrt{ 4 \pi } \rho( r )

\displaystyle V( r ) = 4\pi \int^\infty_0 dr' \, r'^2 \frac{ 1 }{ r_> } \, \rho( r' ) = 4\pi \left( \frac{ 1 }{ r } \int^r_0 dr' \, r'^2 \, \rho( r' ) + \int^\infty_r dr' \, r' \, \rho( r' ) \right)

\rho( r ) が半径 Rまで値を持つ場合において、その内側で V( r < R ) が定数になる \rho の形を求めることを試みる。
定数であるから、微分するとポテンシャルはゼロになる条件を使って、
\displaystyle \frac{1}{4\pi} \frac{ d V }{ d r } = 0 = - \frac{ 1 }{ r^2 } \int^r_0 dr' \, r'^2 \, \rho( r' ) + \frac{1}{r} r^2 \, \rho( r ) + R \, \rho( R ) - r \, \rho( r ) \\ \displaystyle \therefore \int^r_0 dr' \, r'^2 \, \rho( r' ) = r^2 R \, \rho( R )

もう一度微分すると、
\displaystyle r^2 \, \rho( r ) = 2 r R \, \rho( R ) \\ \displaystyle \therefore \rho( r ) = 2 \rho( R ) \frac{ R }{ r } \\ \displaystyle \therefore \rho( R ) = 0 \rightarrow \rho( r ) = 0

したがって、定数ポテンシャルを与えるような電子密度は存在しないことが証明される。

仮に、一様な電荷 \rho( r ) = \rhoを与えたとしても、
\displaystyle V( r < R ) = 4\pi \rho \left( \frac{ 1 }{ r } \int^r_0 dr' \, r'^2 + \int^R_r dr' \, r' \right) = \rho \left( \frac{1}{ r } \frac{ 4\pi r^3 }{ 3 } + \frac{ 4\pi R^2 - 4\pi r^2 }{ 2 } \right) \\ \displaystyle V( r \ge R ) = 4\pi \rho \frac{ 1 }{ r } \int^R_0 dr' \, r'^2 = \rho \frac{1}{ r } \frac{ 4\pi R^3 }{ 3 } = \frac{ Q }{ r } \quad ( \rho \equiv Q \frac{ 3 }{ 4 \pi R^3 } )
という様に、位置に依存した関数となってしまう。

化学ポテンシャルの温度依存性。

固体物理 数学 量子力学

前回、Sommerfeld展開についてまとめた。
koideforest.hatenadiary.com
\displaystyle I = \int^{\infty}_0 d\varepsilon \, g( \varepsilon ) f( \varepsilon ) \approx \int^{\mu}_0 d\varepsilon \, g( \varepsilon ) + \frac{ \pi^2 }{ 6 } \frac{ 1 }{ \beta^2 } g'( \mu )

今回は、 g = \rho（状態密度）として適用し、化学ポテンシャルの温度依存性について調べる。
\displaystyle \frac{N}{2} = \int^{\varepsilon_F}_0 d\varepsilon \, \rho( \varepsilon ) \\ \displaystyle \frac{N}{2} \approx \int^{\mu}_0 d\varepsilon \, \rho( \varepsilon ) + \frac{ \pi^2 }{ 6 } \frac{ 1 }{ \beta^2 } \rho'( \mu ) \\ \displaystyle \therefore \int^{\varepsilon_F}_0 d\varepsilon \, \rho( \varepsilon ) - \int^{\mu}_0 d\varepsilon \, \rho( \varepsilon ) \approx \frac{ \pi^2 }{ 6 } \frac{ 1 }{ \beta^2 } \rho'( \mu )

ここで、 \varepsilon_F - \mu \ll 1 とすると、Tayler展開を用いて次のように書ける。
\displaystyle G( \mu ) \equiv \int^{\mu}_0 d\varepsilon \, \rho( \varepsilon ) \approx G( \varepsilon_F ) + G'( \varepsilon_F )( \mu - \varepsilon_F ) \\ \displaystyle \qquad = G( \varepsilon_F ) + \rho'( \varepsilon_F )( \mu - \varepsilon_F ) \\ \displaystyle \therefore \int^{\varepsilon_F}_0 d\varepsilon \, \rho( \varepsilon ) - \int^{\mu}_0 d\varepsilon \, \rho( \varepsilon ) \approx - \rho( \varepsilon_F )( \mu - \varepsilon_F ) \\ \displaystyle - ( \mu - \varepsilon_F ) \approx \frac{ \pi^2 }{ 6 } \frac{ 1 }{ \beta^2 } \frac{ \rho'( \mu ) }{ \rho( \varepsilon_F ) } \\ \displaystyle \mu \approx \varepsilon_F \left( 1 - \frac{ \pi^2 }{ 6 } \frac{ \varepsilon_F }{ ( \varepsilon_F \beta )^2 } \frac{ \rho'( \mu ) }{ \rho( \varepsilon_F ) } \right)

ここで、自由電子の状態密度が
\displaystyle \rho( \varepsilon ) = \frac{3}{4} N \varepsilon_F^{-3/2} \sqrt{\varepsilon}
であるため、 \rho'( \varepsilon_F ) \propto \varepsilon_F^{-1} となることに注意すると、 \beta \varepsilon \gg 1の低温条件（低温と言えど普通室温も含まれる）で、
\displaystyle \mu \approx \varepsilon_F \left( 1 - \frac{ \pi^2 }{ 6 } \frac{ \varepsilon_F }{ ( \varepsilon_F \beta )^2 } \frac{ \rho'( \varepsilon_F ) }{ \rho( \varepsilon_F ) } \right)
とできる。これは第二項が既に O( ( \beta \varepsilon_F )^{-2})のオーダーであるため、 \rho'( \mu ) のTayler展開は一次で十分となる。

式から読み取れることをまとめると、

  1. \rho'(\varepsilon_F) = 0のときには、化学ポテンシャルは温度で変化しない。
  2. 化学ポテンシャルが増えるか減るかは、 \rho'(\varepsilon_F) の符号で決まる。
  3. ただし、低温領域( \beta \mu \gg 1 )での近似式なので、状態密度が \varepsilon_Fを中心に左右非対称であれば \rho'(\varepsilon_F) = 0でも高温条件で温度変化が起きる。これはFermi分布関数の形から言えることである。
  4. 逆に言えば、状態密度が \varepsilon_Fを中心に左右対称の場合、化学ポテンシャルは変化しない。
  5. 自由電子の場合には、常に \rho'(\varepsilon_F) > 0 であるため、化学ポテンシャルは減少する。

化学ポテンシャルは一般に計算しにくい量なので、何となくの振る舞いがわかるだけでも、イメージし易くなって良いと思います。

Sommerfeld展開。

固体物理 量子力学 数学

Sommerfeld展開は、Fermi分布関数 f( \varepsilon ) が掛かった関数の積分値 Iに関する展開式で、以下のように表される。
\displaystyle f( \varepsilon ) = \frac{ 1 }{ e^{ \beta ( \varepsilon - \mu ) } + 1 } \quad \left( \beta = \frac{ 1 }{ k_B T } \right) \\ \displaystyle I = \int^\infty_0 d\varepsilon \, g( \epsilon ) f( \varepsilon ) = \int^\mu_0 d\varepsilon \, g( \varepsilon ) + \frac{ \pi^2 }{ 6 } \frac{ 1 }{ \beta^2 } g'( \mu ) + \cdots
室温程度であれば、( g(\varepsilon)が十分滑らかな関数である限り) 第二項までの打ち切りが良い近似として成り立つことが知られている。
この展開を導出してみることにする。

先に数学的な項目をまとめ、その後に導出を議論する。

  • 公式等々

この展開の導出には、以前の記事で紹介した公式（？）が役に立つ。
koideforest.hatenadiary.com
\displaystyle \int^\infty_0 \frac{ x^{p-1} }{ e^x + 1 } dx = \left( 1 - \frac{ 1 }{ 2^{p-1} }\right) \zeta( p ) \Gamma( p )

次の式変形に慣れておくと、この先で少し楽になる。
\displaystyle - \frac{ d }{ dx } \frac{ 1 }{ e^x + 1 } = \frac{ e^x }{ ( e^x + 1 )^2 } \\ \displaystyle \qquad = \frac{ 1 }{ e^x + 2 + e^{-x} } = \frac{ 1 }{ ( e^x + 1 )( e^{ -x } + 1 ) } = \frac{1}{2} \frac{ 1 }{ 1 + \cosh( x ) }
したがって、 - \frac{ d }{ dx } \frac{ 1 }{ e^x + 1 } は偶関数である。

考えたいのは、次の積分である。
\displaystyle \int^\infty_{-\infty} x^{p} \left( - \frac{ d }{ dx } \frac{ 1 }{ e^x + 1 } \right) dx
この積分が求まれば、Sommerfeld展開は7割ぐらい終わったと言っても良いと思う。

p が奇数だと、被積分関数が奇関数なのでゼロになる。
一方、 p が偶数 p = p_e だとすると、
\displaystyle \int^\infty_{-\infty} x^{p_e} \left( - \frac{ d }{ dx } \frac{ 1 }{ e^x + 1 } \right) dx = 2 \int^\infty_0 x^{p_e} \left( - \frac{ d }{ dx } \frac{ 1 }{ e^x + 1 } \right) dx \\ \displaystyle \qquad = 2 \left[ - \left. \frac{ x^{p_e} }{ e^x + 1 } \right|^\infty_0 + p_e \int^\infty_0 \frac{ x^{p_e - 1} }{ e^x + 1 } dx \right]
第一項はゼロになるため、第二項のみが残る。
第二項では積分公式が使えるので、
\displaystyle \int^\infty_{-\infty} x^{p_e} \left( - \frac{ d }{ dx } \frac{ 1 }{ e^x + 1 } \right) dx = 2 p_e \left( 1 - \frac{ 1 }{ 2^{p_e-1} }\right) \zeta( p_e ) \Gamma( p_e ) \\ \displaystyle \qquad \left( = \left( 2 - \frac{ 1 }{ 2^{p_e - 2} }\right) \zeta( p_e ) \Gamma( p_e + 1 ) \right)
最後の式変形までは要らないかもしれないが、一応 \Gamma( p_e + 1 ) = p_e \Gamma( p_e ) を使ってまとめた。

例えば、 p_e = 2 のとき、 \zeta( 2 ) = \pi^2 / 6 , \Gamma( n ) = ( n - 1 )! であるため、
\displaystyle \int^\infty_{-\infty} x^2 \left( - \frac{ d }{ dx } \frac{ 1 }{ e^x + 1 } \right) dx = \left( 2 - 1 \right) \zeta( 2 ) \Gamma( 3 ) = \frac{ \pi^2 }{ 6 } \cdot 2 = \frac{ \pi^2 }{ 3 }
と求まる。

ただし、 p = 0 は例外扱いしなければならない（不定になる）。
\displaystyle \int^\infty_{-\infty} x^{0} \left( - \frac{ d }{ dx } \frac{ 1 }{ e^x + 1 } \right) dx = \int^\infty_{-\infty} \left( - \frac{ d }{ dx } \frac{ 1 }{ e^x + 1 } \right) dx \\ \displaystyle = 2 \int^\infty_0 \left( - \frac{ d }{ dx } \frac{ 1 }{ e^x + 1 } \right) dx = - 2 \left. \frac{ 1 }{ e^x + 1 } \right|^\infty_0 = 1

  • Sommerfeld展開

まず、 g( \varepsilon)の積分値を定義しておく。
\displaystyle G( \varepsilon ) \equiv \int^\varepsilon_0 d\varepsilon' \, g( \varepsilon' )
G( 0 ) = 0 であることがポイントである。

これによって、積分 Iを Gで表すことが出来る。
\displaystyle I = \left. G( \varepsilon ) f( \varepsilon ) \right|^\infty_0 - \int^\infty_0 d\varepsilon \, G( \varepsilon ) \frac{ d f( \varepsilon )}{ d \varepsilon }
第一項はゼロになるため、第二項が重要であることがわかる。。
ここで、 -\frac{ d f( \varepsilon )}{ d \varepsilon } は \varepsilon = \mu 近傍でのみ値を持つ（ T = 0 でdelta関数に帰着する）ので、積分の中の G(\varepsilon) についても \varepsilon = \mu 近傍だけ考えれば良い。

したがって、 G(\varepsilon) を \varepsilon = \mu 近傍でTayler展開すると、
\displaystyle G(\varepsilon) = G( \mu + ( \varepsilon - \mu ) ) = G( \mu ) + G'( \mu ) ( \varepsilon - \mu ) + \frac{ 1 }{ 2! } G''( \mu )( \varepsilon - \mu )^2 + \cdots \\ \displaystyle \qquad \equiv \sum_{n=0} \frac{ G^{(n)}( \mu ) }{ n! } ( \varepsilon - \mu )^n \\ \displaystyle \therefore I = \sum_{n=0} \frac{ G^{(n)}( \mu ) }{ n! } \int^\infty_0 d\varepsilon \, ( \varepsilon - \mu )^n \left( - \frac{ d f( \varepsilon )}{ d \varepsilon } \right) \\ \displaystyle \qquad = \sum_{n=0} \frac{ G^{(n)}( \mu ) }{ n! } \int^\infty_{-\mu} d\varepsilon' \, ( \varepsilon' )^n \left( - \frac{ d }{ d \varepsilon' } \frac{ 1 }{ e^{ \beta \varepsilon'} + 1 } \right)

ここで、上で導出した公式を利用するために、 \beta \mu \gg 1 であると仮定する。
\varepsilon' ( = \varepsilon - \mu ) = 0 近傍が重要であることを思い出すと、積分範囲の下限を無限大に飛ばしても積分値への影響は少ないはずなので、良い近似として成り立つと期待できる。
これによって、公式が利用できる。
\displaystyle I \approx \sum_{n=0} \frac{ G^{(n)}( \mu ) }{ n! } \frac{ 1 }{ \beta^n }\int^\infty_{-\infty} d\varepsilon' \, ( \beta \varepsilon' )^n \left( - \beta \frac{ d }{ d \beta \varepsilon' } \frac{ 1 }{ e^{ \beta \varepsilon'} + 1 } \right) \\ \displaystyle \qquad = \sum_{n=0} \frac{ G^{(n)}( \mu ) }{ n! } \frac{ 1 }{ \beta^n }\int^\infty_{-\infty} dx \, x^n \left( - \frac{ d }{ dx } \frac{ 1 }{ e^{ x } + 1 } \right) \\ \displaystyle \qquad = G( \mu ) + \sum_{k = 1} \frac{ G^{(2k)}( \mu ) }{ (2k)! } \frac{ 1 }{ \beta^{2k} } \left( 2 - \frac{ 1 }{ 2^{2k - 2} }\right) \zeta( 2k ) \Gamma( 2k + 1 ) \\ \displaystyle \qquad = G( \mu ) + \sum_{k =1} \frac{ g^{(2k - 1)}( \mu ) }{ (2k)! } \frac{ 1 }{ \beta^{2k} } \left( 2 - \frac{ 1 }{ 2^{2k - 2} }\right) \zeta( 2k ) \Gamma( 2k + 1 )
これで、Sommerfeld展開の導出が完了した。

OPW（直交化平面波法）

固体物理 量子力学

そもそもの大前提として、系が周期性を持ち、ハミルトニアンが（周期性を保つ）並進操作に対して不変であるとする。
\displaystyle H | \Psi_k \rangle = E(k) | \Psi_k \rangle \\ \displaystyle [ T_{R}, H ] = 0 \\ \displaystyle T_{R} | \Psi_k \rangle = e^{ i k \cdot R} | \Psi_k \rangle

平面波で基底を展開すると、強く局在したcore状態との直交性を記述しようとするために高い波数が必要になってしまい、実際に計算するにはとても不便である。
そのため、予めcore状態と直交させた平面波（OPW）で固有状態を展開することを考える。
\displaystyle \phi_k( x ) = \frac{ 1 }{ \sqrt{ \Omega } } e^{ikx} \\ \displaystyle \psi_k( x ) = \phi_k( x ) - \sum_R \sum_{c(R)} \langle \phi_{c(R)} | \phi_k \rangle \phi_{c(R)}( x )

\psi_kが実際にcore状態と直交していることを確認する。
\displaystyle \langle \psi_k | \phi_{c(R)} \rangle = \langle \phi_k | \phi_{c(R)} \rangle - \sum_{R'} \sum_{c(R')} \langle \phi_k | \phi_{c(R')} \rangle \langle \phi_{c(R')} | \phi_{c(R)} \rangle \\ \displaystyle \qquad = \langle \phi_k | \phi_{c(R)} \rangle - \sum_{R'} \sum_{c(R')} \langle \phi_k | \phi_{c(R')} \rangle \delta_{c(R'), c(R)} \\ \displaystyle \qquad = \langle \phi_k | \phi_{c(R)} \rangle - \langle \phi_k | \phi_{c(R)} \rangle = 0 \\ \displaystyle \because \langle \phi_{c(R')} | \phi_{c(R)} \rangle = \delta_{c(R'), c(R)}
core状態と言っているが、OPWではcore状態が自然軌道の性質を持つことが要求されている（他のサイトの軌道と直交する）点が重要である。
そのため、semi-core状態があると、必要な波数がやはり多くなることが予想される。

ただし、OPW同士は直交していないことに注意。
\displaystyle S_{k,k'} \equiv \langle \psi_k | \psi_{k'} \rangle \\ \displaystyle \qquad = \langle \phi_{k} | \phi_{k'} \rangle - \sum_{R_1} \sum_{c(R_1)} \langle \phi_k | \phi_{c(R_1)} \rangle \langle \phi_{c(R_1)} | \phi_{k'} \rangle \\ \displaystyle \qquad \qquad - \sum_{R_2} \sum_{c(R_2)} \langle \phi_{c(R_2)} | \phi_{k'} \rangle \langle \phi_{k} | \phi_{c(R_2)} \rangle \\ \displaystyle \qquad \qquad + \sum_{R_1, R_2} \sum_{c(R_1), c(R_2)} \langle \phi_{c(R_1)} | \phi_k \rangle \langle \phi_{k'} | \phi_{c(R_2)} \rangle \langle \phi_{c(R_1)} | \phi_{c(R_2)} \rangle \\ \displaystyle \qquad = \delta_{k,k'} - 2 \sum_{R} \sum_{c(R)} \langle \phi_k | \phi_{c(R)} \rangle \langle \phi_{c(R)} | \phi_{k'} \rangle \\ \displaystyle \qquad \qquad + \sum_{R_1, R_2} \sum_{c(R_1), c(R_2)} \langle \phi_{c(R_1)} | \phi_k \rangle \langle \phi_{k'} | \phi_{c(R_2)} \rangle \delta_{R_1,R_2} \delta_{c(R_1),c(R_2) } \\ \displaystyle \qquad = \delta_{k,k'} - \sum_{R} \sum_{c(R)} \langle \phi_k | \phi_{c(R)} \rangle \langle \phi_{c(R)} | \phi_{k'} \rangle

もちろん、OPWと平面波も直交していない。
\displaystyle \langle \phi_k | \psi_{k'} \rangle \\ \displaystyle \qquad = \langle \phi_{k} | \phi_{k'} \rangle - \sum_{R} \sum_{c(R)} \langle \phi_{c(R)} | \phi_{k'} \rangle \langle \phi_{k} | \phi_{c(R)} \rangle \\ \displaystyle \qquad = \delta_{k,k'} - \sum_{R} \sum_{c(R)} \langle \phi_k | \phi_{c(R)} \rangle \langle \phi_{c(R)} | \phi_{k'} \rangle \\ \displaystyle \qquad = S_{k,k'}

（保証があるわけではないが） \psi_kで真の固有状態 \Psi_kを展開出来るとすると、
\displaystyle \Psi_k ( x ) \approx \sum_K a_{k,K} \psi_{k+K}( x ) = \Phi_k( x ) - \sum_R \sum_{c(R)} \langle \phi_{c(R)} | \Phi_k \rangle \phi_{c(R)}( x ) \\ \displaystyle \Phi_k ( x ) = \sum_K a_{k,K} \phi_{k+K}( x )

ここで、core軌道がハミルトニアンの固有状態になっているとすると（周りの影響によってcore軌道のエネルギーが影響を受けないとすると）、
\displaystyle H | \phi_{c(R)} \rangle \approx E_{c(R)} | \phi_{c(R)} \rangle

\displaystyle H \Psi_k( x ) = E( k ) \Psi_k( x ) \\ \displaystyle \qquad = H \Phi_k( x ) - \sum_R \sum_{c(R)} \langle \phi_{c(R)} | \Phi_k \rangle E_{c(R)} \phi_{c(R)}( x ) \\ \displaystyle \qquad = \sum_K \left( H \phi_{k+K}( x ) - \sum_R \sum_{c(R)} \langle \phi_{c(R)} | \phi_{k+K} \rangle E_{c(R)} \phi_{c(R)}( x ) \right) a_{ k, K } \\ \displaystyle \langle \phi_{k+K'} | H | \Psi_k \rangle = E( k ) \sum_K S_{ k+K', k + K } a_{k, K} \\ \displaystyle \qquad = \sum_K \left( H_{k+K',k+K} - \sum_R \sum_{c(R)} E_{c(R)} \langle \phi_{k+K'} | \phi_{c(R)} \rangle \langle \phi_{c(R)} | \phi_{k+K} \rangle \right) a_{ k, K } \\ \displaystyle \qquad \equiv \sum_K \tilde{H}_{k+K',k+K} a_{ k, K } \\ \displaystyle \therefore S^{-1} \tilde H \, a = E \, a

aは平面波の展開係数でもあるから、平面波が感じる有効ハミルトニアンは \tilde Hであると解釈出来る。
Hに比べて、 \tilde Hは電子間相互作用が弱くなっているため、結晶中で自由電子に近い振る舞いが発現する理由になっている。

力学的相似（ビリアル定理）。

力学

（系の）ビリアルとは以下の量のことを呼ぶ（らしい）。
\displaystyle B \equiv \left\langle \sum_a \vec r_a \cdot \frac{ \partial U }{ \partial \vec r_a } \right\rangle

古典（質点）力学においては、変数は時間であるため、これは時間平均量として解釈出来る。
\displaystyle \left\langle f \right\rangle = \lim_{\tau \rightarrow \infty} \frac{ 1 }{ \tau } \int^\tau_0 dt\, f( t )

ビリアル定理とは、

  1. ポテンシャルエネルギー Uが座標 \vec rの同次関数
  2. 系の運動が有限領域内に限られている

場合に成り立つ、運動エネルギー Tとポテンシャルエネルギー Uの関係であり、次の形で与えられる。
\displaystyle 2 \left\langle T \right\rangle = k \left\langle U \right\rangle \quad ( U( \alpha \vec r ) = \alpha^k U( \vec r ) )

例えば、クーロンポテンシャルでは k = -1であるため、
\displaystyle 2 \left\langle T \right\rangle = - \left\langle U \right\rangle \\ \displaystyle \left\langle T \right\rangle = - E \quad ( \left\langle T \right\rangle + \left\langle U \right\rangle = \left\langle E \right\rangle = E )
という関係があることが分かる。
逆にこれを利用することで、例えばオリジナルのポテンシャルでは解きにくいから近似を入れるような時に、「ビリアル定理をなるべく破らないように近似を入れる」という指針を立てることが出来る。

以下、ビリアル定理を証明する。

運動エネルギー Tは、速度 vに対して二次の同次関数であるため、同次関数に対するオイラーの定理が成り立つ。
\displaystyle \sum_a \vec v_a \cdot \frac{ \partial T }{ \partial \vec v_a } = 2 T

確認しておくと、
\displaystyle \sum_a \vec v_a \cdot \frac{ \partial T }{ \partial \vec v_a } = \sum_a \vec v_a \cdot \frac{ \partial }{ \partial \vec v_a } \left( \frac{1}{2} m v_a^2 \right) = \sum_a m v_a^2 = 2 T \\ \displaystyle \left( \frac{ \partial }{ \partial \vec v } \equiv \left( \frac{ \partial }{ \partial v_x }, \frac{ \partial }{ \partial v_y }, \frac{ \partial }{ \partial v_z} \right) \right)

同様に、ポテンシャル Uが位置 rに対して k次の同次関数である条件の下では、
\displaystyle \sum_a \vec r_a \cdot \frac{ \partial U }{ \partial \vec r_a } = k \, U \\ \therefore B = k \left\langle U \right\rangle
と表される。

ここからは、上手い式変形を使って、 Tと Uが繋がるところを見る。
運動量 \vec p = m \vec vを用いれば
\displaystyle 2T = \sum_a \vec p_a \cdot \vec v_a = \sum_a \vec p_a \cdot \frac{ d \vec r_a }{ d t } = \frac{ d }{ d t } \left( \sum_a \vec p_a \cdot \vec r_a \right) - \sum_a \left( \frac{ d \vec p_a }{ d t } \right) \cdot \vec r_a \\ \displaystyle \qquad = \frac{ d }{ d t } \left( \sum_a \vec p_a \cdot \vec r_a \right) + \sum_a \left( \frac{ \partial U }{ \partial \vec r_a } \right) \cdot \vec r_a \\ \displaystyle \because \frac{ d \vec p }{ d t } = - \frac{ \partial U }{ \partial \vec r }
最後の式変形には、運動方程式を用いている。

一応、もう一度まとめておくと、
\displaystyle 2T = \frac{ d }{ d t } \left( \sum_a \vec p_a \cdot \vec r_a \right) + \sum_a \left( \frac{ \partial U }{ \partial \vec r_a } \right) \cdot \vec r_a
ここまでは厳密に成り立つ。 Uが同次関数であることも、有限領域の運動であることも使っていない。

Uが同次関数であれば、
\displaystyle 2T = \frac{ d }{ d t } \left( \sum_a \vec p_a \cdot \vec r_a \right) + k U
と簡単になる。

更に、運動が有限領域に限られる場合、（時間）平均を取ると、
\displaystyle 2\left\langle T \right\rangle = \lim_{\tau \rightarrow \infty } \frac{ 1 }{ \tau } \sum_a \left( \vec p_a( \tau ) \cdot \vec r_a( \tau ) - \vec p_a( 0 ) \cdot \vec r_a( 0 ) \right) + k \langle U \rangle \\ \displaystyle \qquad = k \langle U \rangle = B \\ \displaystyle \because \vec p_a( \tau ) \cdot \vec r_a( \tau ) - \vec p_a( 0 ) \cdot \vec r_a( 0 ) \ll \infty
よってビリアル定理が得られた。

参考文献：ランダウ=リフシッツ（力学）

力学的相似。

力学

ランダウ=リフシッツ（力学）に載っていた問題。

一般座標 qと時間 tと質量 mをスケール変換することを考える。
\displaystyle \vec r \rightarrow \alpha \, \vec r \\ \displaystyle t \rightarrow \beta \, t \\ \displaystyle m \rightarrow \gamma \, m

このスケール変換に対して、ラグランジュ方程式が不変に保たれる条件を考える。

速度等の普遍的な量は以下のスケール変換を受ける。
\displaystyle \vec v = \frac{ d \vec r }{d t } \rightarrow \frac{\alpha}{ \beta } \, \vec v \\ \displaystyle \vec p = m \vec v \rightarrow \frac{ \alpha \gamma }{ \beta } \, \vec v \\ \displaystyle T = \sum_a \frac{p_a^2}{2m} \rightarrow \left( \frac{ \alpha }{ \beta } \right)^2 \gamma \, T

ポテンシャル Uに関しても、以下のようなスケール変換を受けるとする。
\displaystyle U \rightarrow c( \alpha, \beta, \gamma ) \, U \equiv c \, U

ラグランジュアンおよびラグランジュ方程式は以下で与えられる。
\displaystyle L = T - U \\ \displaystyle \frac{ d }{ d t } \left( \frac{ \partial L }{ \partial v } \right) - \frac{ \partial L }{ \partial r } = 0

ラグランジュ方程式は、全体を定数倍しても不変に保たれるので、そのような関係を満たすような kを考える。
ラグランジュアン L が定数倍であるためには、
\displaystyle L \rightarrow \left( \frac{ \alpha }{ \beta } \right)^2 \gamma \, T - c \, U \\ \displaystyle \therefore \left( \frac{ \alpha }{ \beta } \right)^2 \gamma = c
という条件が課せられる。

このとき、以下のようにラグランジュ方程式が単に定数倍化されて不変に保たれることが分かる。
\displaystyle \beta^{-1} \frac{ d }{ d t } \left( c \left( \frac{ \alpha }{ \beta } \right)^{-1} \frac{ \partial L }{ \partial v } \right) - c \alpha^{-1} \frac{ \partial L }{ \partial r } = \frac{ c }{ \alpha } \left( \frac{ d }{ d t } \left( \frac{ \partial L }{ \partial v } \right) - \frac{ \partial L }{ \partial r } \right) = 0

したがって、ラグランジュ方程式を不変にするスケール変換の条件は、
\displaystyle \left( \frac{ \alpha }{ \beta } \right)^2 \gamma = c \\ \displaystyle \frac{ \alpha }{ \beta } \sqrt{ \gamma } = \sqrt{ c } \\ \displaystyle \beta = \alpha \sqrt{ \frac{ \gamma }{ c } }

スケール後の変数をプライム付きで表すと、
\displaystyle \frac{ t' }{ t } = \left( \frac{ r' }{ r } \right) \sqrt{ \frac{ m' }{ m } \frac{ U }{ U' } }

このスケール変換による特性を用いて、以下の問題に答えることが出来る。

  • ポテンシャルが同一で、同じ軌跡を描く運動において、質量を変化させたときの運動時間の変化。

\displaystyle \frac{ t' }{ t } = \sqrt{ \frac{ m' }{ m } }
質量の平方根だけ余計に時間が掛かる（遅くなる）ことがわかる。

  • ポテンシャルをスケール変化させても運動が同じ軌跡を描く場合の運動時間変化（質量変化無し）。

\displaystyle \frac{ t' }{ t } = \sqrt{ \frac{ U }{ U' } }
ポテンシャルに対しては、質量の平方根の逆数だけ早くなることがわかる。

例えば、長さを変えたときに体積がどれくらいのオーダーで大きくなるかが分かっているだけでもだいぶ便利なように、各変数の詳細な量が分からなくても、変化するオーダーが一瞬で分かってしまうのは、驚嘆に値する。