時間依存しないシュレーディンガー方程式と変分原理 1/2 （計算手法の説明）

時間依存しないシュレーディンガー方程式を数値的に解くための、変分原理について述べていきます。
変分原理はエネルギー固有値の上限を与えます(厳密な解よりも変分原理で求めた解は必ず大きくなる)。
基底関数を適切に選べば高精度の結果が与えられます。ですが、基底関数をどのように選べばよいか？の指標は無く、経験則に頼らざるを得ない状況が出てきます。

ここでは、変分原理によってなぜエネルギー固有値が求まるのかを述べて行きます。

もしも変分原理より正しい解を得たければ、shooting methodと呼ばれる方法も推奨します。
恐らくこちらの方が色々なポテンシャルに対して有効ですが、
・計算速度が遅い
・発散の可能性がある
・変分原理よりも経験値が必要
です。

まとめ

区間[a,b]で定義された時間依存しないシュレーディンガー方程式
$
\displaystyle \left[-\frac{1}{2}\frac{d^2}{dx^2}+V(x)\right]\psi(x)=E\psi(x)
$
の解である固有関数を
$
\displaystyle \psi(x)=\sum_n c_n \varphi_n(x) \ \ …(6)
$
の形で固定端条件($\psi(a)=\psi(b)=0$)の下で探します。
ここで、$\varphi_n(x)$は区間$[a,b]$で正規直交化されたsin関数系
$
\displaystyle \varphi_n(x)=\sqrt{\frac{2}{b-a}}\sin\left(\frac{n\pi}{b-a}(x-a)\right),\;\;\;(n=1,2,\cdots,N)
$
です。この時、行列要素$H_{ij}=\langle\varphi_i|-\frac{1}{2}\frac{d^2}{dx^2}+V(x)|\varphi_j\rangle$は
$
\begin{align}
H_{ij}&=K_{ij}+V_{ij} \\
&\hspace{1em}K_{ij}=\frac{1}{2}\left(\frac{j\pi}{b-a}\right)^2\delta_{ij} \\
&\hspace{1em}V_{ij}=\frac{1}{j}\sum_{m=2}^{M-1}\omega_m^{(GL)} \cos\left(\frac{\pi}{2}x_m^{(GL)}\right)
\sum_{n=0}^{j-1}(-1)^n \sin\left(\frac{i\pi}{2j}(x_m^{GL}+2n+1)\right)V(X_n^{(m)}) \\
&\hspace{2em}X_{n}^{(m)}=\frac{h}{2}(x_m^{(GL)}+2n+1)+a,\;\;h=\frac{b-a}{j}
\end{align}
$
と求められます($j\ge i$を推奨)。ここで、$\omega_m^{(GL)},\ x_m^{(GL)}$はM点のGauss-Lobatto求積法の区間[-1,1]における重みと求積点で、
$
\displaystyle \omega_m=\frac{2}{M(M-1)[P_{M-1}(x_m)]^2},\ \ x_m\ne \pm 1
$
であり、求積点は関数
$
\displaystyle \frac{P_{M-1}(x)+P_{M}(x)}{x-1}
$
のゼロ点です。$P_{n}(x)$は$n$次Legendre多項式を表します。

解説
時間依存から時間依存しないシュレーディンガー方程式へ
間接法
数値計算プログラム
正規直交系を用いない場合
研究分野での数値計算法について
参考文献

計算の具体例は
時間依存しないシュレーディンガー方程式と変分原理 2/2 （計算例）をご覧ください。

[adsense1]

解説

対象とする問題を導出します。何はともあれ、時間依存シュレーディンガー方程式から出発します。

時間依存から時間依存しないシュレーディンガー方程式へ

時間依存シュレーディンガー方程式は
$
\displaystyle i\frac{\partial \Psi(x,t)}{\partial t}=\hat{H}(x)\Psi(x,t)
$
と記述されます。ここで、簡単ために原子単位系($\hbar=1, m=1$)を考えます。
変数分離の形で解を仮定します。
$
\displaystyle \Psi(x,t)=\psi(x)U(t)
$
代入して移行すると以下の式を得ます。

$
\begin{align}
i\psi(x)\frac{\partial U(t)}{\partial t}=U(t)\hat{H}\psi(x) \\
i\frac{1}{U(t)}\frac{\partial U(t)}{\partial t}=\frac{1}{\psi(x)}\hat{H}\psi(x) = E
\end{align}
$

2番目の式の右辺$E$は何らかの定数です。
左辺は時間のみを含む式、真ん中の式は位置のみを含む式なので、いつどんな時でもこの関係式が成り立っていないといけないという条件からこれらはなんらかの定数でなければなりません。その定数を$E$と置きました。

すなわち、２つの式
$
\begin{align}
i\frac{\partial U(t)}{\partial t}=E U(t)\;\;\; \cdots(1)\\
\hat{H}\psi(x) = E\psi(x)\;\;\; \cdots(2)
\end{align}
$
が導けます。

$t$に関する式は$\hat{H}(x)$に依らず解けてしまい、
$U(t)=e^{-i E t}$(1)
です。位置に関する式は
$
\hat{H}\psi(x) = E\psi(x)
$(2)
を満たさなければなりません。この(2)が時間依存しないシュレーディンガー方程式と呼ばれているものです。
これを数値計算で解いていくことを考えます。

左から$\psi^*$を掛けて全空間で積分しますと、
$
\begin{align}
\int \psi^* \hat{H} \psi dx – E \int \psi^* \psi dx = 0\;\;\; \cdots(3)
\end{align}
$
という式を得ます。この問題は
条件$\int \psi^* \psi dx=1$の下で、$\int \psi^* \hat{H} \psi dx$の極値を求めよ、
という問いになります。
上記問題を考えた時、(3)と同一の式になります。

すなわち、原理的にはどんな$\psi(x)$を与えても$E(=E[\psi])$の値を計算できるわけです。
その無数に考えられる$\psi$の中で、特に(3)の解が停留する$\psi$、すなわち
$
\displaystyle \delta \int \psi^* (\hat{H}-E) \psi dx = 0
$
となる$\psi$を探せ、という問題になります。

↑がもしも停留する$\psi$が見付かれば、それはシュレーディンガー方程式と同一になるということは以下の通りに示せます。

$
\displaystyle \delta \int \psi^* (\hat{H}-E) \psi dx =0 \;\;\;\cdots (4)
$
厳密な$\psi$の変分に対して停留、すなわち$\psi\to \psi+\delta \psi$、もしくは複素平面で$\frac{\pi}{2}$ずれた方向$\psi\to \psi+i\delta \psi$に対して0でなければなりません。
(4)に対して２つの方向に対して変分を取ると、
$
\begin{align}
\int (\delta \psi)^* (\hat{H}-E)\psi dx &+ \int \psi^* (\hat{H}-E)(\delta\psi) dx =0 \\
-i\int (\delta \psi)^* (\hat{H}-E)\psi dx &+ i\int \psi^* (\hat{H}-E)(\delta\psi) dx =0
\end{align}
$
と2式が導けます。2式は変分をとった方向にかかわらず同時に満たされていなければなりません。同時に成立するためには
$
(\hat{H}-E)\psi=0,\;\; \text{or}\;\; \psi^*(\hat{H}^{\dagger}-E^*)=0
$
を満たしている必要があります。これは時間依存しないシュレーディンガー方程式と一致します。

以上をまとめますと、

厳密な$\psi$の変分はシュレーディンガー方程式に一致する
↓
近似的に$\psi$を変分するとシュレーディンガー方程式の近似的な解に帰着されるだろう
↓
近似的に$\psi$を構成し、変分問題を解く

代表的な方法として、
手法1)直接法…$\psi$の具体的な形を変分パラメータと共に与える
手法2)間接法…変分パラメータとしてある正規直交基底関数系を採用する
の２つがありますが、数値計算で使われる間接法を説明します。

間接法

間接法による変分原理は、

固有値の上限を与える（変分原理によるエネルギー固有値は必ず厳密なエネルギー固有値よりも大きくなる）
基底関数の数を増やすほど精度は上がる(完全系の場合)

という特性があります。

さて、今解きたい問題は
条件
$
\displaystyle \int \psi^* \psi dx =1
$
の下で
$
\displaystyle \int \psi^* \hat{H} \psi dx
$
の極値を与える$\psi$を探すこと、です。

すなわち
$
\displaystyle \int \psi^* \hat{H} \psi dx – E \int \psi^* \psi dx =0 \ \ …(5)
$
を考えて変分します。求めたい波動関数は
$
\displaystyle \psi(x)=\sum_n c_n \varphi_n(x) \ \ …(6)
$
と展開できる、と考えます。ここで、$c_n$は変分パラメータ、$\varphi_n(x)$は区間[a,b]で定義されている正規直交関数系です。正規直交関数系$\{\varphi_n(x)\}$は以下の性質を持ちます。
$
\displaystyle \int_a^b \varphi_n(x) \varphi_m(x)dx =\delta_{nm}
$
この制約が無くても変分は可能ですが、以降の式が簡略になるので、正規直交関数だけを考えます。

(6)を(5)に代入すると、
$
\begin{align}
\sum_{n}\sum_{m} c_n^* c_m \left[H_{nm}-E \delta_{nm}\right]=0 \\
\end{align}\ \ …(7)
$
なります。ここで、
$
\begin{align}
H_{nm}&=\int_a^b \varphi_n(x) \hat{H} \varphi_m(x)dx
\end{align}
$
とおきました。(7)式を満たす係数$\{c_n\}$の内、自明な解($c_n=0\ (for\ all\ n)$)はすぐ見つかります。しかし、この解は$\psi(x)=0$であり、条件$
\displaystyle \int \psi^* \psi dx =1$を満たしません。なので解としては不適切です。

では非自明な解を見つけましょう。(7)式の解を探すため、極値を探します。今、係数の組を変分パラメータとみなしているので、(変分パラメータで微分したもの)=0を計算します。すなわち
$
\displaystyle \frac{\partial}{\partial c_n^*} (for\ all\ n)
$
を(7)式に作用させます。すると、
$
\begin{align}
\sum_{n}\sum_{m} \left[H_{nm}-E\delta_{nm}\right]c_m=0
\end{align}
$
となります。$n,m$を行列の行、列ととらえれば、

$
\begin{align}
\left(
\begin{array}{cccc}
H_{11} & H_{12} & \cdots & H_{1N}\\
H_{21} & H_{22} & \cdots & H_{2N}\\
\ddots & \ddots & \ldots & \ddots \\
H_{N1} & H_{N2} & \cdots & H_{NN}
\end{array}
\right)
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
=
E
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
\end{align}
$

となるわけです。ということは、左辺の正方行列を対角化すれば右辺のE,すなわち固有値が求まる、という事になります。

行列要素の計算

行列を対角化するためには行列要素$H_{nm}$が分かっていないとなりません。
これからこの行列要素を計算していきます。

ハミルトニアンが
$
\displaystyle \hat{H}=-\frac{1}{2}\frac{d^2 }{d x^2} + V(x)
$
として表現されている場合、行列要素は
$
\begin{align}
H_{nm}&=K_{nm}+V_{nm} \\
&\hspace{3em} K_{nm}=-\frac{1}{2}\int_a^b \varphi_n(x)\frac{d^2}{dx^2} \varphi_m(x)dx \\
&\hspace{3em} V_{nm}=\int_a^b \varphi_n(x)V(x) \varphi_m(x)dx
\end{align}
$
と書き表され、これらを数値計算で求めます。

$K_{nm}$の計算

まず、微分を含んだ
$
\displaystyle K_{nm}=-\frac{1}{2}\int_a^b \varphi_n(x)\frac{d^2}{dx^2} \varphi_m(x)dx
$
の項を考えます。微分は数値計算で行いません。
というのも正規直交関数系はこちらで用意する知っている関数系なので、求めることが出来てしまうからです。
多くの正規直交関数系には漸化式と微分方程式が存在します。一番簡単な基底としてsin関数系を考えてみましょう。
区間$[a,b]$で正規直交化されたsin関数系は
$
\displaystyle \varphi_n(x)=\sqrt{\frac{2}{b-a}}\sin\left(\frac{n\pi}{b-a}(x-a)\right),\;\;\;(n=1,2,\cdots,N)
$
と書き表せます。sin関数なので、2階微分はもちろん
$
\displaystyle \frac{d^2}{dx^2}\varphi_n(x)=-\left(\frac{n\pi}{b-a}\right)^2\varphi_n(x)
$
を満たします。これを積分の中に代入して、
$
\begin{align}
K_{nm}&=-\frac{1}{2}\int_a^b \varphi_n(x)\frac{d^2}{dx^2} \varphi_m(x)dx \\
&=\frac{1}{2}\left(\frac{n\pi}{b-a}\right)^2\int_a^b \varphi_n(x) \varphi_m(x)dx \\
&=\frac{1}{2}\left(\frac{n\pi}{b-a}\right)^2\delta_nm
\end{align}
$
となります。積分を計算することなく求める事が出来ました。

より一般には一回微分の存在しない形、もしくは漸化式によって手計算で微分を含む積分を片付けてしまったり、部分積分を用いて何とか計算します。
正規直交系であれば、うまく整理することによって誤差が全く入らず、かつ高速に計算できます。

$V_{nm}$の計算

さて、残る問題は
$
\displaystyle V_{nm}=\int_a^b \varphi_n(x)V(x) \varphi_m(x)dx
$
の計算です。基底関数の次数が低い場合、この計算は何も問題がありません。
調和ポテンシャル$V(x)=x^2$の時、$n=1,m=2$の時の被積分関数は

と非常になだらかですが、$n=5,m=99$の時の被積分関数は

となります。なので細かい分点をとることになり、非常に積分が難しいです。

離散フーリエ変換によって求められそうな気もしますが、周期境界条件が満たされていないことは明らかなのであきらめましょう。
有限区間で早く振動する関数の積分は難しいです。二重指数関数型数値積分も振動する関数は苦手です。

こういった場合の計算方法として、
・ゼロ点の位置が分かっている場合、そこで区間を分割して積分する
という直観的な方法があります[2]。
ゼロ点の位置が分かっていれば、その間の被積分関数の振る舞いはシンプルな形になるので滑らかな関数として考えられます。

また、端点での値がゼロだと分かっている場合、Gauss-Lobatto求積法が非常に有用です。

Gauss-Lobatto求積法による積分の計算

Gauss-Lobatto求積法はGauss-Legendre求積法と同じく被積分関数を多項式で近似します[1]。
Gauss-Legendre求積法と異なる点は端点で分点を持つということです。端点を明示的に書けば、積分は
$
\displaystyle \int_{-1}^1 f(x)dx\sim \frac{2}{N(N-1)}(f(-1)+f(1))+\sum_{i=2}^{N-1}\omega_i f(x_i)
$
と書けます。ここで、$\omega_i,x_i$はGauss-Lobatto求積法における重み、求積点です。
Gauss-Lobatto求積法の重みは
$
\displaystyle \omega_i=\frac{2}{N(N-1)[P_{N-1}(x_i)]^2},\ \ x_i\ne \pm 1
$
であり、求積点は
$
\displaystyle \frac{P_{N-1}(x)+P_{N}(x)}{x-1}
$
のゼロ点です。$P_{n}(x)$は$n$次Legendre多項式を表します。

Gauss-Lobatto求積法は
$f(\pm 1)$が既知の値(特に0)を持つ場合
に便利なことが分かります。
関数$f(\pm 1)$がゼロであれば、その点を計算する必要性がなくなりますし、見かけ上の特異点であってもこちらで値を指定すれば参照されることはなく、計算できないという状況に陥りません。

ゼロ点に挟まれる区間の積分

ゼロ点に挟まれる区間の積分は、結局以下の式になります[2]。
$
\displaystyle \varphi_j(x)=\sqrt{\frac{2}{b-a}}\sin\left(\frac{j\pi}{b-a}\right)
$
とすると、積分は
$
\begin{align}
&\int_a^b \varphi_i(x) V(x) \varphi_j(x) dx \\
&=\sum_{n=0}^{j-1} \int_{x_n}^{x_{n+1}} \varphi_i(x) V(x) \varphi_j(x)dx \\
&\approx \sum_{n=0}^{j-1}\sum_{m=2}^{M-1}\frac{h}{2} \omega_m^{(GL)}\varphi_i(X_n^{(m)}) V(X_n^{(m)}) \varphi_j(X_n^{(m)}) \\
&=\frac{h}{2}\sum_{m=2}^{M-1}\omega_m^{(GL)}\varphi_j(X_0^{(m)})S_j^{(m)} \\
\end{align}
$
ここで、
$
\begin{align}
&x_n=nh+a,\ h=\frac{b-a}{j},\ (n=0,1,\cdots,j) \\
&X_n^{(m)}=\frac{h}{2}x_m^{(GL)}+\frac{x_{n+1}+x_n}{2} \\
&X_{n+1}^{(m)}=X_n^{(m)}+h \\
&S_j^{(m)}=\sum_{n=0}^{j-1}(-1)^n \varphi_i(X_n^{(m)})V(X_n^{(m)})
\end{align}
$
と置きました。
$x_m^{(GL)},\ \omega_m^{(GL)}$は区間$[-1,1]$でのM点Gauss-Lobatto求積法の求積点と重みです。
計算を可能な限り減らすため、$\varphi_j(x)$の周期性を用いています。

数値計算プログラム

行列要素を組み立てた後は適当なルーチンを用いて対角化します。
今回の場合、基底関数をsin関数に選んだため、強制的に固定端条件が課されます。
なので、この基底関数系を用いるといかなる問題も端点$a,b$に無限に高いポテンシャルの壁がある、と読み替えることが出来ます。
よって行列は対称行列、エルミートです。故に固有値は必ず実数です。

この問題を解く場合、実対称行列を対角化するルーチン、lapackのdsyevdを用いるのが良いでしょう。
対角化のプログラムに関する詳細は対角化へどうぞ。

以下のプログラムは行列要素を上記手法で作り、対角化し、得られた波動関数を規格化し、出力するものです。

▼ここクリックでこの場に展開

program main
! sikinote
! date : 2017/04/09 (yyyy/mm/dd)
! author : sikino
implicit none
integer::N,Nout
double precision::a,b,s
double precision,allocatable::H(:,:),E(:)

a=-30d0 ! Left side range
b= 30d0 ! Right side range (a < b)
N=200 ! Number of sin basis function
! with boundary condition psi(a)=psi(b)=0
Nout=10 ! Number of output wave functions

!$ call omp_set_num_threads(1)

allocate(H(1:N,1:N),E(1:N)); H=0d0; E=0d0

! Construct matrix element of Hamiltonian
! Basis function : orthnormalized sin function [a,b]
! Integral method : Gauss-lobatto quadrature
! using the zero points of integrand function.
! (12 point for each zeros interval)
call Hamiltonian(N,H,a,b)

! Diagonalize hermite matrix using lapack "dsyevd" for real hermitian matrix
call Diagonalize(N,H,E)
! Normalize eigenvector
call Normalize(N,H)

! Output eigenenergy to a terminal & files
call Eigenenergy_output(N,E)
! Output wavefunction to files
call Eigenfunction_output(N,H,E,Nout,a,b)
! Output potential to a file
call Potential_output(400,1.2*a,1.2*b)

stop
end program main

function V(x)
implicit none
double precision::V,x

V=0.5d0*x*x
!V=0d0

return
end function V

!==================

subroutine Eigenfunction_output(N,H,E,Nout,a,b)
implicit none
integer,intent(in)::Nout,N
double precision,intent(in)::a,b,H(1:N,1:N),E(1:N)
! a,b : output range, (a < b)

double precision::x,dx,wfx,evec(1:N)
integer::i,j,Np
character(48)::fname

Np=1000
dx=(b-a)/dble(Np)
do j=1,Nout
call intchar(j,fname)
open(11,file=trim(fname))
evec(1:N)=H(1:N,j)
do i=0,Np
x=a+i*dx
call wavefunction(x,N,evec,a,b,wfx)
write(11,'(2e20.10e2,e14.7e2)')x,wfx,E(j)
enddo
close(11)
enddo

return
end subroutine Eigenfunction_output

subroutine Eigenenergy_output(N,E)
implicit none
integer,intent(in)::N
double precision,intent(in)::E(1:N)

integer::i,Nout

Nout=10
if(N.lt.10)Nout=N

do i=1,Nout
write(6,*)i,E(i)
enddo
open(10,file="eigenenergy.d")
do i=1,N
write(10,*)i,E(i)
enddo
close(10)

return
end subroutine Eigenenergy_output

subroutine Potential_output(N,a,b)
implicit none
integer,intent(in)::N
double precision,intent(in)::a,b

integer::i
double precision::h,x,V
external::V

open(20,file="potential.d")
h=(b-a)/dble(N)
do i=0,N
x=a+i*h
write(20,*)x,V(x)
enddo
close(20)

return
end subroutine Potential_output

subroutine wavefunction(x,N,evec,a,b,wf)
implicit none
integer,intent(in)::N
double precision,intent(in)::a,b,x,evec(1:N)
double precision,intent(out)::wf

integer::i
double precision::sinbasis
external::sinbasis

wf=0d0
do i=1,N
wf=wf+evec(i)*sinbasis(i,x,a,b)
enddo

return
end subroutine wavefunction

subroutine Normalize(N,H)
implicit none
integer,intent(in)::N
double precision,intent(inout)::H(1:N,1:N)

integer::i,k
double precision::s

do k=1,N
s=0d0
do i=1,N
s=s+H(i,k)*H(i,k)
enddo
H(1:N,k)=H(1:N,k)/sqrt(s)
enddo

return
end subroutine Normalize

subroutine Hamiltonian(N,H,a,b)
implicit none
integer,intent(in)::N
double precision,intent(in)::a,b
double precision,intent(out)::H(1:N,1:N)

integer::i,j
double precision,parameter::pi=3.14159265358979323846264338327950288d0
double precision::Vij,c,kin
external::Vij
!$ integer::omp_get_thread_num
!$ external::omp_get_thread_num

kin=pi/(b-a)
!$omp parallel do schedule(static,1)
do i=1,N
c=i*kin
H(i,i)=0.5d0*c*c+Vij(i,i,a,b)
enddo
!$omp end parallel do

!$omp parallel do schedule(static,1)
do i=1,N
do j=i+1,N
H(i,j)=Vij(i,j,a,b)
enddo
enddo
!$omp end parallel do

return
end subroutine Hamiltonian

function Vij(i,j,a,b)
implicit none
integer,intent(in)::i,j
double precision,intent(in)::a,b
double precision::Vij

! i,j>=1
!
! Vij=\int_a^b \varphi_i(x) V(x) \varphi_j(x) dx
!
! Using zero points of \varphi_j(x) with Gauss-Lobatto quadrature.
! (j>=i)
!
! \varphi_n(x) : Orthnormalized sin basis for [a,b]
! = \sqrt{\frac{2}{b-a}}\sin{(n\pi\frac{x-a}{b-a})}
!
integer::n,m,k,l,sgn
integer,parameter::Ngl=12
double precision,parameter::pi=3.14159265358979323846264338327950288d0
double precision::xgl(1:Ngl),wgl(1:Ngl),p2x(1:Ngl),xz(1:Ngl)
double precision::h,Smj,s,ijp2
double precision::V
external::V

k=i; l=j
if(i.gt.j)then
k=j; l=i
endif

h=(b-a)/l
call gauss_lobatto_xw(Ngl,xgl,wgl)
ijp2=0.5d0*pi*k/dble(l)
p2x(1:Ngl)=0.5d0*pi*xgl(1:Ngl)
xz(1:Ngl)=a+0.5d0*h*(xgl(1:Ngl)+1d0)

s=0d0
do m=2,Ngl-1
Smj=0d0; sgn=1
do n=0,l-1
Smj=Smj+sgn*sin(ijp2*(xgl(m)+2d0*n+1d0))*V(xz(m)+n*h)
sgn=-sgn
enddo
s=s+wgl(m)*cos(p2x(m))*Smj
enddo
Vij=s/dble(l)

return
end function Vij

function sinbasis(n,x,a,b)
implicit none
integer,intent(in)::n
double precision,intent(in)::x,a,b
! Orthnormalized sin basis for [a,b]

double precision,parameter::pi=3.14159265358979323846264338327950288d0
double precision::sinbasis

sinbasis=sqrt(2d0/(b-a))*sin(n*pi*(x-a)/(b-a))

return
end function sinbasis

subroutine Diagonalize(N,A,ev)
implicit none
integer,intent(in)::N
double precision,intent(inout)::A(1:N,1:N)
double precision,intent(out)::ev(1:N)
!for Hermite matrix
!---sikinote

integer::lda,lwork,liwork,info
double precision,allocatable::work(:)
double precision::qw(1:3)
integer,allocatable::iwork(:)
integer::qiw(1:3)
character(1)::job,uplo

job="V"; uplo="U"; lda=N

!size query of dsyevd or zheevd parameter.
call dsyevd(job,uplo,N,0,lda,0,qw,-1,qiw,-1,info)
if(info.ne.0)then
write(6,'(A)')" something wrong at dsyevd(size query) in Diagonalize"
write(6,'(A)')" program stop"
write(6,'(A)')" info --> ",info
stop
endif

!just to be safe, +1
lwork=idint(qw(1))+1; liwork=qiw(1)+1
allocate(work(1:lwork),iwork(1:liwork))
work(1:lwork)=0d0; iwork(1:liwork)=0

!diagonalize.
call dsyevd(job,uplo,N,A,lda,ev,work,lwork,iwork,liwork,info)
if(info.ne.0)then
write(6,'(A)')" something wrong at subroutine dsyevd in Diagonalize"
write(6,'(A)')" program stop"
write(6,'(A)')" info --> ",info
stop
endif

deallocate(work,iwork)

return
end subroutine Diagonalize

subroutine gauss_lobatto_xw(N,x,w)
implicit none
integer,intent(in)::N
double precision,intent(out)::x(1:*),w(1:*)

! Avarable for
! N= 7 12 20

integer::i,M
if(N.eq.7)then
x(4)=0d0
x(5)=0.46884879347071421380377188190876632941d0
x(6)=0.83022389627856692987203221396746513959d0
x(7)=1d0

w(4)=0.48761904761904761904761904761904761905d0
w(5)=0.4317453812098626234178710222813622779d0
w(6)=0.2768260473615659480107004062900662935d0
w(7)=0.047619047619047619047619047619047619048d0
elseif(N.eq.12)then
x( 7)=0.1365529328549275548640618557396938969d0
x( 8)=0.39953094096534893226434979156696690053d0
x( 9)=0.63287615303186067766240485444365585824d0
x(10)=0.81927932164400667834864158171690266069d0
x(11)=0.94489927222288222340758013830321871361d0
x(12)=1d0

w( 7)=0.2714052409106961770002883384996028404d0
w( 8)=0.2512756031992012802932444121475961833d0
w( 9)=0.2125084177610211453583020773668662702d0
w(10)=0.15797470556437011516467106270034026825d0
w(11)=0.09168451741319613066834259413407928631d0
w(12)=0.015151515151515151515151515151515151515d0
elseif(N.eq.20)then
x(11)=0.08054593723882183797594451815955446302d0
x(12)=0.23955170592298649518240135692708807194d0
x(13)=0.39235318371390929938647470381582436667d0
x(14)=0.53499286403188626164813596182898398301d0
x(15)=0.66377640229031128984640332297115885248d0
x(16)=0.77536826095205587041431752759469134337d0
x(17)=0.86687797808995014130984721461628521396d0
x(18)=0.93593449881266543571618158493062692992d0
x(19)=0.98074370489391417192544643858423091523d0
x(20)=1d0

w(11)=0.1607432863878457490077267264490839339d0
w(12)=0.1565801026474754871581698967936433723d0
w(13)=0.14836155407091682581471301373396660493d0
w(14)=0.13630048235872418448978079298903198174d0
w(15)=0.12070922762867472509942970500239335987d0
w(16)=0.10199149969945081568378120573288652087d0
w(17)=0.0806317639961196031447768461137205717d0
w(18)=0.05718180212756682600475362717324278251d0
w(19)=0.0322371231884889414916050281172940301d0
w(20)=0.0052631578947368421052631578947368421053d0
end if

if(mod(N,2).eq.0)then
M=N/2
else
M=(N-1)/2
endif
do i=1,M
x(i)=-x(N-i+1)
w(i)= w(N-i+1)
enddo

return
end subroutine gauss_lobatto_xw

subroutine intchar(i,fname)
implicit none
integer,intent(in)::i
character(*),intent(out)::fname

if(i.le.9)then
write(fname,'(A,i0,A)')"wf_00",i,".d"
elseif(i.le.99)then
write(fname,'(A,i0,A)')"wf_0",i,".d"
elseif(i.le.999)then
write(fname,'(A,i0,A)')"wf_",i,".d"
else
write(6,*)"Not supported, at subroutine filename."
stop
endif

return
end subroutine intchar

[adsense2]

コンパイルと実行は

$ ifort -mkl main.f90
$ ./a.out

などとしてください。lapackを用いますので、それと一緒にコンパイルしてください。

ラグランジュの未定乗数について

ラグランジュの未定乗数は特に物理的な意味を持たない魔法の変数、として置きますが、今回の場合は未定乗数がエネルギー固有値というちゃんとした意味を持つのです。

正規直交系を用いない場合

$
\begin{align}
\sum_{n}\sum_{m} c_n^* c_m \left[H_{nm}-E I_{nm}\right]=0 \\
\sum_{n}\sum_{m} c_n^* c_m \left[J_{nm} H_{nm}-E \right]=0
\end{align}\ \ …(7)
$
となります。ここで、
$
\begin{align}
H_{nm}&=\int_a^b \varphi_n(x) \hat{H} \varphi_m(x)dx\\
J_{nm}&=I_{nm}^{-1}\\
I_{nm}&=\int_a^b \varphi_n(x) \varphi_m(x)dx
\end{align}
$
と置きました。行列の形では、
$
\begin{align}
\left(
\begin{array}{cccc}
J_{11}H_{11} & J_{12}H_{12} & \cdots & J_{1N}H_{1N}\\
J_{21}H_{21} & J_{22}H_{22} & \cdots & J_{2N}H_{2N}\\
\ddots & \ddots & \ldots & \ddots \\
J_{N1}H_{N1} & J_{N2}H_{N2} & \cdots & J_{NN}H_{NN}
\end{array}
\right)
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
=
E
\left(
\begin{array}{cccc}
1 & 0 & \cdots & 0\\
0 & 1 & \cdots & 0\\
\ddots & \ddots & \ldots & \ddots \\
0 & 0 & \cdots & 1
\end{array}
\right)
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
\end{align}
$
と書けます。

$
\begin{align}
\sum_{n}\sum_{m} c_n^* c_m \left[H_{nm}-E I_{nm}\right]=0 \\
\sum_{n}\sum_{m} c_n^* c_m \left[I_{nm}^{-1} H_{nm}-E \right]=0
\end{align}
$

研究分野での数値計算法について

実際の研究の場では対角化はしますが、上記方法によって行列要素を作ることはしません。遅いからです。
特別な基底関数を用いると行列要素$H_{nm}$が1~2項、しかも非常に簡単な式に置き換えることが出来ます。精度も速度も上記方法以上です。
この計算方法に関してはこれ以上は述べることはしませんが、こういった方法も存在する、とだけ述べておきます。

参考文献

[1]P. J. Davis, P. Rabinowitz著, 森正武訳, 「計算機による数値積分法」, Gauss-Lobatto求積法 p.88
[2]P. J. Davis, P. Rabinowitz著, 森正武訳, 「計算機による数値積分法」, ゼロ点に挟まれる区間の積分 p.131
[3]ランダウ＝リフシッツ,「量子力学１」第3刷, p81-84

（更新停止）

時間依存しないシュレーディンガー方程式と変分原理 1/2 （計算手法の説明）

まとめ