2019-12-09

fortran のオプショナル引数

練習として、オプショナル引数を使ったサンプルプログラムを作った。

program ex

call ex_sub( 1, z = 3 )

contains

subroutine ex_sub( x, y, z )
implicit none
!integer, intent( in ) :: x, y, z
integer, optional :: x, y, z

if ( .not. present( x ) ) x = 0
if ( .not. present( y ) ) y = 0
if ( .not. present( z ) ) z = 0

print *, x, y, z

return
end subroutine ex_sub

end program ex

省略された引数の値を初期化したいので、引数を変更することを禁止するintentとは相性が悪いと思われる。

2019-12-09

fortran における複素数の入力

プログラミング

以下のファイルの中身を複素数の変数に読み込ませたいとする。

# test input file: test.inp
1.00000  2.00000

結構探したが、fortranで外部から複素数を読み込ませるときの上手い方法が見つからなかった。
結局、実数変数に一時的に格納し、後で複素数に変換するという形しか無さそう。

program test

implicit none
integer, parameter :: INPUT = 10
integer, parameter :: DP = selected_real_kind( 8 )
complex( DP ), parameter :: I = cmplx( 0, 1 )

real( DP ) :: temp( 2 )
complex( DP ) :: c
character( 99 ) :: filename

filename = 'test.inp'
open( INPUT, file = filename, status = 'old' )
read( INPUT, * ) temp( 1 : 2 )
c = cmplx( temp( 1 ), temp( 2 ) )
! or c = temp( 1 ) + I * temp( 2 )

end program test

2019-12-09

fortran で End of file を検出する。

プログラミング

fortranでファイルを限界まで読み続けると、End of file でエラーが返って来る。この辺がなんとも原始的な香りを感じさせる。
この End of file を検出/処理するためには、以下の二通りが考えられる。

END で行番号へ飛ぶ。

Fortran 入門: ファイル操作

implicit none
character( 99 ) :: input

open( UNIT, FILE, STATUS = 'old' )
do while ( .true. )
    read( UNIT, FORMAT,  END = 999 ) input
    print *, input
end do
999  print *, 'end of file'

End of file 時に、指定した行番号へ飛び、処理を続ける。
しかし、GoTo文がダサいという認識の下では、行番号に飛ぶのは度し難く感じる。

IOSTAT で検出する。

【fortran】ファイル読み込みのコツ open, read | 理系夫婦の方程式

implicit none
integer :: ios
character( 99 ) :: input

open( UNIT, FILE, STATUS = 'old' )
do while ( .true. )
    read( UNIT, FORMAT,  IOSTAT = ios ) input
    if ( ios == -1 ) then
        print *, 'end of file'
        exit
    end if
    print *, input
end do

End of file 時に、IOSTATに -1 が代入され、処理が続けられる。こちらの方がスマートであろう。
ちなみに、End of file 時には変数(input)への代入は行われないため、最後の行の情報がそのまま変数inputの中に入っていることになる。

2019-12-07

LAPACK・LAPACK95の使い方。

プログラミング

LAPACKの使い方でいつも困るので自分用も兼ねてまとめる。
LAPACKのインストールの仕方については沢山サイトがあるので割愛。

サブルーチンのリスト（マニュアル）

LAPACK95 -- Fortran95 interface to LAPACK
Index of LAPACK95 Routines
LAPACK95のサブルーチンのリスト。
スタンダードな感じの無味乾燥としたページ。
LAPACK95についてはこれで事足りる気がする。

LAPACKサンプルプログラム目次
nagによるサブルーチン（LAPACK）のリスト。
サンプルプログラムのソースも載っているが、nag独自のmoduleもガンガン使っているので、そのままコピペしても動かない。
入力データと出力結果の例が見れるのは良いと思うが、ページの見方に慣れが必要。
結局はリンク先のマニュアルページを見ることになる。

LAPACK
NECによるサブルーチン（LAPACK）のリスト。
一ページにサブルーチンの名前と簡単な説明が全部載っているので、何があるかザっと見るには良いページだと思う。
サブルーチン名をクリックした先にある詳細な説明も見易くて良い。

LAPACKとLAPACK95について（１）

LAPACKでは、サブルーチンの名前は「"接頭辞" + "名前"」になっており、変数の型によって接頭辞を変える必要がある。
とは言っても、数値計算では倍精度の実数か複素数を使うことになるので、慣れれば気にならない。
一方、LAPACK95では、「"LA_" + "名前"」がサブルーチンの名前になっていて、変数の型を気にせずに呼び出すことができる。

LAPACKとLAPACK95について（２）

LAPACK・LAPACK95を利用するには、対応したモジュールを使用することをプログラム内で宣言する必要がある。

use f77_lapack  ! LAPACK
use f95_lapack  ! LAPACK95

両方を宣言する必要はない。LAPACKのサブルーチンを使いたいのか、LAPACK95のサブルーチンを使いたいのか、どちらかだけで良い。
この「f77_lapack, f95_lapack」というのは、moduleファイル（*.mod）ファイルの名前でもある。コンパイルする際には「f77_lapack.mod, f95_lapack.mod」が存在するPATHをincludeする必要がある。上手くコンパイルできないときは、一度確認することをお勧めする。筆者の場合には、"include"ディレクトリの中の"lapack95_modules"にmoduleファイルが入っていたため、"include"ディレクトリを指定しても上手く行かなかった。

LAPACKとLAPACK95について（３）

LAPACK と LAPACK95 で引数が異なることがある。
大体の場合、LAPACK95の方が楽ちんだと（勝手に）思っている。
例えば、複素エルミート行列の固有値と固有ベクトルを求める「ZHEEV or LA_HEEV」を見てみると、

! LAPACK
SUBROUTINE ZHEEV (JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO)

! LAPACK95
SUBROUTINE LA_HEEV( A, W, JOBZ=jobz, UPLO=uplo, INFO=info )

LA_HEEVの場合には、二個引数を入れるだけで動いてくれる（後ろの３つはoptionalだから指定しなくても可）が、ZHEEVは毎回全ての引数を用意しないといけない。
また、両者で引数の順番が違うことにも注意。LAPACKの方を知っているからといって、LAPACK95をマニュアル無しで使おうとすると失敗する。

試しにLA_HEEVを使うと、

program check_lapack95

use f95_lapack, only : LA_HEEV

implicit none

integer, parameter :: dp = selected_real_kind( 8 )

real( dp )    :: W( 2 )
complex( dp ) :: A( 2, 2 )

! input/output matrix
A = reshape( [ cmplx( 1, 0 ), cmplx( 0, 2 ), cmplx( 0, -2 ), cmplx( 1, 0 )  ], shape( A ), order = [ 2, 1 ] )

! LA_HEEV( A, W, (JOBS, UPLO, INFO) )
call LA_HEEV( A, W, 'V' )  ! 'V': eigenvector will be stored in A

! A( 1, : ) = i/sqrt{2}, -1/sqrt{2}
! A( 2, : ) = i/sqrt{2}, 1/sqrt{2}
! W = -1, 3

end program check_lapack95

コンパイルする際は、includeファイルとlibraryに関するコンパイルオプションをお忘れなく。

2019-11-06

数演算子の波数表示。

量子力学

数演算子
$\displaystyle n_i = a^\dagger_i a_i$
の波数表示を求める。

まず、素直に $n_i$ をフーリエ級数展開で波数表示にしてみる。
$\displaystyle n_i = \sum_q n_q e^{ i q \cdot r_i } \qquad \left( n_q = \frac{ 1 }{ N } \sum_i n_i e^{ - i q \cdot r_i } \right) \\ \displaystyle a_i = \frac{1}{ \sqrt{N} } \sum_k a_k e^{ i k \cdot r_i } \qquad \left( a_k = \frac{1}{ \sqrt{N} } \sum_i a_i e^{ - i k \cdot r_i } \right)$

$n_q$ に $n_i$ を代入することで $n_q$ を求める。
$\displaystyle n_q = \frac{ 1 }{ N } \sum_i n_i e^{ - i q \cdot r_i } = \frac{ 1 }{ N } \sum_i a^\dagger_i a_i e^{ - i q \cdot r_i } \\ \displaystyle \qquad = \frac{ 1 }{ N^2 } \sum_i \sum_{ k_1, k_2 } a^\dagger_{k_1} a_{k_2} e^{ i ( - k_1 + k_2 ) \cdot r_i } e^{ - i q \cdot r_i } \\ \displaystyle \qquad = \frac{ 1 }{ N } \sum_{ k_1, k_2 } a^\dagger_{k_1} a_{k_2} \delta_{ k_2, k_1 + q } \\ \displaystyle \qquad = \frac{ 1 }{ N } \sum_{ k_1 } a^\dagger_{k_1} a_{k_1 + q} \\ \displaystyle \qquad = \frac{ 1 }{ N } \sum_{ k } a^\dagger_{k } a_{k + q}$

元に戻るかを確認すると、
$\displaystyle \sum_q n_q e^{ i q \cdot r_i } = \frac{ 1 }{ N } \sum_q \sum_{ k } a^\dagger_{k} a_{k + q} e^{ i q \cdot r_i } \\ \displaystyle \qquad = \frac{ 1 }{ N^2 } \sum_q \sum_{ k } \sum_{ l, m } a^\dagger_{l} a_{m} e^{ i k \cdot r_l } e^{ - i ( k + q ) \cdot r_m } e^{ i q \cdot r_i } \\ \displaystyle \qquad = \frac{ 1 }{ N^2} \sum_q \sum_{ k } \sum_{ l, m } a^\dagger_{l} a_{m} e^{ i k \cdot ( r_l - r_m ) } e^{ i q \cdot ( r_i - r_m ) } \\ \displaystyle \qquad = \sum_{ l, m } a^\dagger_{l} a_{m} \delta_{ l,m } \delta_{ i, m } \\ \displaystyle \qquad = a^\dagger_{i} a_{i} = n_i$

$n_q$ の対称性として、 $n_{-q} = n_q^\dagger$ が使える。
$\displaystyle n_{-q} = \sum_i n_i e^{ i q \cdot r _i } \\ \displaystyle \qquad = \sum_i \left( n_i \right)^\dagger e^{ i q \cdot r _i } \\ \displaystyle \qquad = \left( \sum_i n_i e^{ - i q \cdot r _i } \right)^\dagger \\ \displaystyle \qquad = n_{q}^\dagger$

もしくは以下のようにしても示せる。
$\displaystyle n_{-q} = \sum_k a^\dagger_k a_{k-q} \\ \displaystyle \qquad = \sum_{k'} a^\dagger_{ k' + q } a_{ k' } \\ \displaystyle \qquad = \sum_{k'} \left( a^\dagger_{ k'} a_{ k' + q } \right)^\dagger \\ \displaystyle \qquad = \left( \sum_{k'} a^\dagger_{ k'} a_{ k' + q } \right)^\dagger \\ \displaystyle \qquad = n_q^\dagger$

特別に $q = 0$ のとき、
$\displaystyle n_{q=0} = \frac{1}{N} \sum_k a_k^\dagger a_k = \frac{1}{N} \sum_k n_k = \frac{ \hat N }{ N }$
というように、電子数変化の割合を表す演算子になる。

2019-11-06

平均場近似の心。

量子力学固体物理

演算子、もしくは確率変数を $X, Y$ と置き、それらの期待値（平均値）をそれぞれ $\bar X, \bar Y$ すると、
$\displaystyle X Y = ( X - \bar X ) ( Y - \bar Y ) + X \bar Y + \bar X Y - \bar X \bar Y$
という式が成り立つ。
これはまだ厳密である。

ここで、「 $X, Y$ は平均値 $\bar X, \bar Y$ にそれぞれ近いから、「揺らぎ」である $X - \bar X \equiv \delta X$ , $Y - \bar Y \equiv \delta Y$ は凄く小さい」と考えるのが平均場近似である。
ポイントは、 $(\delta X) (\delta Y)$ というように微小量の二次になっているため、「二次微小量を落とす」という視点を持つことである。
そうでなければ、 $X \bar Y - \bar X \bar Y = (\delta X) \bar Y$ という項を残して良い理由が無くなってしまう。
これによって、
$\displaystyle X Y \approx X \bar Y + \bar X Y - \bar X \bar Y$
という近似形が求まる。

2019-11-06

オンサイトクーロンの波数表示。

量子力学固体物理

ハバードモデルにあるオンサイトクーロン項の波数表示
$\displaystyle U \sum_i n_{i,\uparrow} n_{i,\downarrow} = \frac{U}{N} \left( \prod^4_{j=1} \sum_{k_j} \right) a^\dagger_{k_1,\uparrow} a^\dagger_{k_2,\downarrow} a_{k_3,\downarrow} a_{k_4,\uparrow} \delta_{ k_1 + k_2, k_3 + k_4 }$
を証明する。

証明には、以下の関係式が必要である。
$\displaystyle a_i = \frac{ 1 }{ \sqrt{N} } \sum_k a_k e^{ i k \cdot R_i } \qquad \left( a_k = \frac{ 1 }{ \sqrt{N} } \sum_i a_i e^{ - i k \cdot R_i } \right) \\ \displaystyle \sum_i e^{ i k \cdot R_i } = N \delta_{ k, 0 }$

これらは以下の記事で扱ったことがある。
koideforest.hatenadiary.com
koideforest.hatenadiary.com

よって、
$\displaystyle \sum_i n_{i,\uparrow} n_{i,\downarrow} \\ \displaystyle \qquad = \sum_i a^\dagger_{i,\uparrow} a_{i,\uparrow} a^\dagger_{i,\downarrow} a_{i,\downarrow} \\ \displaystyle \qquad = \frac{1}{N^2} \sum_i \left( \prod^4_{j=1} \sum_{k_j} \right) a^\dagger_{k_1,\uparrow} a_{k_2,\uparrow} a^\dagger_{k_3,\downarrow} a_{k_4,\downarrow} e^{ i ( - k_1 + k_2 - k_3 + k_4 ) \cdot R_i } \\ \displaystyle \qquad = \frac{1}{N} \left( \prod^4_{j=1} \sum_{k_j} \right) a^\dagger_{k_1,\uparrow} a_{k_2,\uparrow} a^\dagger_{k_3,\downarrow} a_{k_4,\downarrow} \delta_{ k_1 + k_3, k_2 + k_4 } \\ \displaystyle \qquad = \frac{1}{N} \left( \prod^4_{j=1} \sum_{k_j} \right) a^\dagger_{k_1,\uparrow} a^\dagger_{k_3,\downarrow} a_{k_4,\downarrow} a_{k_2,\uparrow} \delta_{ k_1 + k_3, k_2 + k_4 } \\ \displaystyle \qquad = \frac{1}{N} \left( \prod^4_{j=1} \sum_{k_j} \right) a^\dagger_{k_1,\uparrow} a^\dagger_{k_2,\downarrow} a_{k_3,\downarrow} a_{k_4,\uparrow} \delta_{ k_1 + k_2, k_3 + k_4 }$

nano_exit

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

fortran のオプショナル引数

fortran における複素数の入力

fortran で End of file を検出する。

LAPACK・LAPACK95の使い方。

数演算子の波数表示。

平均場近似の心。

オンサイトクーロンの波数表示。