HPC-Phys 2023.2.8

大規模固有値問題の解法に関わる諸問題

吉田 大輔, 肥山 詠美子
理化学研究所
仁科加速器科学研究センター
ストレンジネス核物理研究室

量子多体問題の第一原理計算で目指すこと

2

定常状態のSchrödinger方程式“量子多体”問題

̂ = T ̂ + V̂ Ψ = EΨ
HΨ
()

時間発展

束縛状態(実固有E状態)
共鳴散乱(複素固有E状態)

反応等

物質のスケールによっていろいろな力の場(有効場)がある
Nuclei

分子

Electron■ Coulomb力 (長距離)■ スピン軌道(LS)力

QuarksNucleons

原子

■ 中心力

ハドロン原子核■ 短距離の強い斥力■ テンソル力(非中心)■ クォーク模型によるカラー閉じ込め力■ スピン軌道(LS)力変分法(非経験的手法)によって、最良のエネルギー固有値 E と固有関数 Ψ を得る！
さまざまな物理量の予測が可能： ⟨A⟩ = ∫ Ψ*(x) Â Ψ(x)dx

変分法

3定常状態のSchrödinger方程式

̂ = T ̂ + V̂ Ψ = EΨ
HΨ
()

r3
r2

「完全系による展開」を再現するようなr1相対運動を記述するためには基底系を選べたなら…

永年方程式を解けば良い

N = 粒子数-1 (空間積分は3N)

H11

H12 ⋯ H1N

H21

H22

⋮

⋮

H2N
⋮

−E

HN1 HN2 ⋯ HNN
̂ J (r1, r2, ⋯)dr
HIJ = ψ*(r1, r2, ⋯)Hψ
I
∫

S11

S12 ⋯ S1N

ψ1(r1, r2, ⋯)S21

S22

S2N

ψ2(r1, r2, ⋯)SN1

SN2 ⋯ SNNψN (r1, r2, ⋯)⋮⋮⋮⋮

SIJ = ψ*(r1, r2, ⋯)ψJ (r1, r2, ⋯)drI∫「完全系による展開」に、いかにして近づけるか
それを目指したとき、どのような問題に直面するか=0

rc = {rci } is the relative coordinatesin c-th Jacobi coordinate system. The parameter↵ rshows...
and={rcci } is the relative coordinates in c-th Jacobi coordinate sows ...
↵ shows
The...many-body wave①function can be expanded in terms
変分法で「厳密」な解を得るためには
4 of
he many-body wave function
canJacobi
be expanded
in terms
of channel functions c foraThegivencoordinate
systemmany-bodywave function
can be expanded in terms of chen Jacobi coordinate systema given Jacobi coordinate systemNc↵「完全系による展開」に、いかにして近づけるかXX↵
NcNormalizationf↵ AN, rc2 , ..., rcX
XLM (r1 , r2 , ..., rN ) =↵c↵c (rc1factorsXXf ↵ A↵
LM (r1 , r2 , ..., rN ) =c (rc1 , rc2 , ..., rcN )↵c=1 (r , r , ..., r )(anti-)symmetrization(r,r,...,r)=fA
LM
12N↵ ↵c c1c2cN
試行波動関数↵c=1

c=1N ↵cX
X=
f↵ AN↵c↵XX
c (rc1 , rc2 , ..., rcNX
X=f ↵ A↵⇥ ⇥LM (rc=1(1.1)
c (rc1 , rc2 , ..., rcN ) ↵c ),=f ↵ A↵c (rc1 , rc2 , ..., rcN )
角度部分
動径部分
Nc↵

r3

↵

factors

↵

c=1

変数分離

↵
c=1
where A↵ is theoperator
todeterminecombination
coefficients toR(r1, r2, r3, ⋯)Θ(r1̂ , r2̂ , r3̂ , ⋯)e A↵ is the operator to determine
combination coefficients
imposethe symmetrizationAor↵ anti-symmetrization,
for atogiven
particleconfiguration
reim
whereis the operator
todeterminecombination
coefficients
tor2
どのような動径関数で展開するかor anti-symmetrization, foragivenparticleconfiguration
representedby ↵,and {f↵}areornormalization
factors forconfigurations,whichshouldbetionanti-symmetrization,for ↵a givenparticle configuration
repre
ormalization factors for ↵theconfigurations,which shouldbe the number
factors tonormalize
total
wavefunction.
Thechannel
c=implicitl
are normalization
factors for↵ configurations,whichc(↵)
shouldbe th■ 調和振動(harmonic
oscillator,HO)関数
otal wave function.
The channel
numberc = c(↵) implicitly
contains
information ofr1
theparticleconfiguration↵.
HO2
the total wave
channelψ (rfunction.
) = hk(ri; rcThe, α) exp−α(rinumber− rc) /2) c = c(↵) implicitly c(iarticle configuration ↵.
the particle
configuration↵. c , the angular part ⇥LM (rc ) is usually
In a channel
function
簡便であるが大きな展開が必要
n a channel function c , thepart ⇥LMharmonics(rc ) is usuallya vector-coupled
productofInangularsolidspherical
suitablefortherelative
coordinata channel function c , the angular part ⇥LM (rc ) is usuallyav
lid spherical harmonics suitable
relativeoftheJacobimomenta,system,explicitly-correlatedGaussian,ECG)関数
is,forin thetermsof thecoordinates
coupling
thetheangularof whichsolid■ 顕に相関するガウス(sphericalharmonics
suitableoffor
relative
coordinates
ECGangular momenta,
2
h is, in terms of the coupling of theψ(r)=Aexp(−αriαi)
which is, in terms ofh the couplingof the angular momenta, ih⇥LM (rc ) =h HOより少ない展開で厳密解に近づける[Y`1 ,m1 (r̂c1i)Y`2 ,m2 (r̂c2 )]L1,2 ,M1,2 Y`3 ,m3 (r̂c3i ) L12,3 ,M[Y`1 ,m1 (r̂c1 )Y`2 ,m2 (r̂c2 )]L1,2 ,M1,2 Y`3 ,m3 (r̂c3 )... Y`N ,mN (r̂cN )M (rc ) =⇥LM (rc ) =[Y`1 ,m1 (r̂c1 )Y`2L,mr̂c2)]Y`3 ,m3 (r̂c3LM)
2 (,M
12,312,3 L1,2 ,M1,2L12,3 ,M12
「厳密」=「良い基底」で展開する= C`c1 ,mc1 ,`c2 ,mc2 ;L12 ,M12 CL12 ,M12 ,`c3 ,mc3 ;L12,3 ,M12,3 ...

rc = {rci } is the relative coordinatesin c-th Jacobi coordinate system. The parameter↵ rshows...
and={rcci } is the relative coordinates in c-th Jacobi coordinate sows ...
Nc↵XX↵shows...
The many-body wave②function can be expanded in termsof
変分法で「厳密」な解を得るためには
5=f↵ A↵ in termsrcN ) ⇥ ⇥(1.1)
c (rc1 , rc2LM (rc ),he many-body wave function
canJacobi
be expandedof, ...,
channelfunctionsc for
aThegiven
coordinatesystem
many-bodywave↵
c=1 function can be expanded in terms of chen Jacobi coordinate systema given Jacobi coordinate system
Nc↵「完全系による展開」に、いかにして近づけるかwhere Athe operator
to determine
coefficientsto imposeX
the symmetrizaX↵ is secular
Theequation
is givenNcombination
by↵cNormalization(rrN ) = representedf↵ AN(rand, r{fX
Xtion or anti-symmetrization,
for
a
given
particleconfiguration} rc
LM
1 , r2 , ...,↵c↵ by c↵,c1factors
c2 ,↵...,
02
3
2
3
1
23X
X
(r1 , r2 , H..., rNfactors) H= forf↵configurations,A, rSc2 , ...,rcN )beLM normalization↵c (rc1which↵c=1 (rare↵shouldthefactors
to normalize(anti-)symmetrization
···H
S···
S
(r,r,...,r)=fA
11121N
11121Nccc
LM12N↵ ↵c1 , rc21, ..., rcN )factors7↵ channelc=1
equation
istotal
given
Bwave
6 by
7
6C 6 2 7 ofthe試行波動関数
function.
The
numberc=c(↵)
implicitly
contains
information
H
H···
H
S
S···S
↵
↵
c=1
21
22
2N
21222Ncc
B6767Nc 7 C 6X
XE
the particle
configuration↵.
767C
6B6 3..
..2.. 7 = 02313↵
変数分離Nc5@44 .
4c1 ,. rc25, ..., rcN.
=f↵ AN↵c↵ 5AX
Xc (rH12 · · · In H
SS···S
1N
11121N
1ccX
X
a channel
function,theangularpart⇥(r)isusuallyavector-coupled
productcLM
c=f
A(r,r,...,r)⇥⇥(r),(1.1)
767C
67H
H···
H
S
S···S
↵
↵cc1
c2
cN
LMc↵
c=1
N
1N2NN
N
1N
2N
N
N,c ..., r
H22 ·of· · solid
H2NSS···Scccccccc222N
2=fA(r,rc 7c6 21
7C67↵↵cc1c2cN )spherical
harmonics
suitablefortherelative
coordinatesoftheJacobisystem,
角度部分
動径部分
c=1
7 E 6 .. ↵7C 6 .. 7 ↵= 0 (1.8)
c=1which is, in 5
terms of4 wherethetheoperator
angularmomenta,
5A4determine
5. coupling.
A↵ is ofthe
tocombination
coefficients toR(r1, r2, r3, ⋯)Θ(r1̂ , r2̂ , r3̂ , ⋯)h determineeAcombination
toimposethe symmetrizaH· operator
HNc Nc toS· · SNccoefficientsNc↵
2 is· ·theNc 1 orSanti-symmetrization,
NNc
Nc ai givention
forparticleconfigurationc2whereA↵ isthe· operator
to determine
combination
coefficients toreim⇥LM(rc ) =Infinitesimally-shifted[Yfor(r̂↵c3represented)configurations,... byY`N↵,(r̂expansio)should
1.2
lobeor anti-symmetrization,
given
particle`1are,mac1 )Y`2 ,m,mcN {f↵}1 (r̂2 (r̂c2 )]Lconfiguration
3GaussianNandfactorsforwhichbe,M1,2 Y`3 ,m1,2
tion
ornormalization
anti-symmetrization,fora givenparticleconfiguration
repre
L12,3 ,M12,3LM
ormalization factors for ↵theconfigurations,whichshouldbethefactors
tonormalize
total wave function.
The
channel numberwhichc = c(↵)
implicitl
are normalization
factors for↵ configurations,shouldbe th
角運動量 M(L)の固有状態として許されるような球面調和関数のカップリングで展開= The
C`c1 ,mCL12Gaussian
otal wave Anfunction.
channelnumberc,M=12 ,`c(↵)
information
ofis,M12;L↵.
alternative
expression
thebasis
especiallyfor`>
0c1 ,`c2 ,mc2 ;L12of
c3 ,mc3implicitly
12,3 ,Mfunction
12,3 ... contains
theparticle
configuration
the total
wave function.
The channel
number c = c(↵) implicitly cnitesimally-shifted
Gaussian
lobeexpansion
特殊関数積分のアプローチは、多体への拡張が困難(ボトルネック)
article configuration ↵.
2the`c2 )...Y⌫n↵.
rcN⇥Y`particle()r̂channel)Y(r̂(r̂cN),
In(rfunction
theangularpart ⇥LM (rc ) is (1.2)
usuallyc1 configuration`,m`c , cNi
c1 ,mc1ac2
c2= Nn` ri eY,mn`m i`m (r̂i )nexpressiona channelof
functionangularpart especially⇥LMharmonics(rc ) isforusuallya vector-coupled
product
c , theGaussianlobe関数を用いる展開法により、多体の積分が可能になった．
the Gaussianfunction`>
0is
ofbasissolid
spherical
suitablefortherelative
coordinatInachannelfunction,theangularpart⇥(r)isusually
avc
LMc
where
C
is
the
Clebsch-Gordancoefficient.
`i ,mi ,`j ,mj ;L,Msuitable for the relative coordinatesK
lid spherical harmonics
oftheJacobimomenta,system,X
whichis,intermsofthecouplingof
theangular2
12of
solid
sphericalharmonics
suitablefortherelative⌫n (ri "i D`m,k )coordinates`⌫n r iNn` lim
C`m,k. corre(rterms= ofNn`
e inYthe(r̂channel
radialfunction
are represented
by eexpansions by theh is,theriparts
coupling
theangularmomenta,n`mini ) The`mofi) =`which is, in terms ofh"!0the(⌫couplingof the angular momenta, in "i )k=1hi2)Y`2 ,m2 (r̂c2 )]L ,M Y`3 ,m3 (r̂c3 )K(rc ) =⇥[Y(r̂LM`,mc1X
11hi L12,3 ,M1,21,21⌫n (ri 1"i D`m,k )= = We[YN`1n`,m(r̂c1 )Y)]`m,kY`for(r̂c3), r. GTFs.
..., rYr̂cN )Ce(1.9)M (rc )任意の粒子数でハミルトニアン行列`2 ,m`2 (r̂c2the,m3s-type`, N,mN (が生成可能̂
1lim
3=should
calculateintegralL,MHψ*(r,⋯)Hψ(r⋯)dr
1,21,2(rc ) =[Y`IJ1 ,m1∫(r̂c1r̂c2)]Y`3 ,m3 (r̂c3LM)
1 `2L2,mJ L1 ,M2"!0 (⌫⇥I )Y2 (,M12,312,3n "LMi)
1,21,2L12,3 ,M12
k=1= C`c1 ,mc1 ,`c2 ,mc2 ;L12 ,M12 CL12 ,M12 ,`c3 ,mc3 ;L12,3 ,M12,3 ...

「厳密な」解に辿り着くために必要な基底関数系の次元
動径波動関数

6■ 粒子数-1 個の相対座標■ 座標系の取り方 (生成・解離過程を適切に記述する)
N = 4 体問題の例 (3個の相対座標)

N2N3 = 20

N2N1= 20

N3N1N1 = 20⋯

N3N1N2N3N3N2N2N1⋯+ 同種粒子の交換

Nchannel ~ 80

粒子数、粒子種に依る

Nbasis ∼ Nchannel × N1 ⋅ N2 ⋅ N3= 640,000

Nbasis ∼ Nchannel ×N−1

∏
i=1Ni

N (-body)# channels# basisTotal35~ 30~ 1302020

104-105
106-1078~ 1000

20> 1010!

Wagrnin「厳密な」解に辿り着くための計算要件

7

対角化演算が計算時間の大部分を占めるようになる

105

CPU時間(秒)104

我々が目指す多体問題
で必要になる領域

行列生成

103
102
101

対角化

100
10-1

粗い解
-2

10

102

厳密解
103

104

105

基底関数の数106

107* スーパーコンピュータシステムITO(九州大学 情報基盤研究開発センター),
1 node, 36 core, 168 GB.

3体問題でのベンチマークデータ高速な大規模行列(1010次元以上)の対角化が必要

8定常状態のSchrödinger方程式“量子多体”問題

̂ = T ̂ + V̂ Ψ = EΨHΨ()
H11H12 ⋯ H1NH21

H22H2N⋮⋮⋮HN1 HN2 ⋯ HNN−E

時間発展

束縛状態(実固有E状態)共鳴散乱(複素固有E状態)

S11

S12 ⋯ S1N

ψ1(r1, r2, ⋯)S21S22S2N

ψ2(r1, r2, ⋯)⋮SN1⋮⋮SN2 ⋯ SNN

⋮

反応等=0

ψN (r1, r2, ⋯)計算要件✔ 対称(定常状態)、複素固有値問題(共鳴状態)→ 非対称✔ 基本的には密(dense)行列✔ 8-byte integer 長の大きな配列 + 演算領域が必要✔ 最低固有値から数〜20個程度(?)、固有値・固有ベクトルを得られればよい・並列版対角化アルゴリズムが必要
・“条件数”問題はクリアできる？

An alternative expressionof the Gaussianfunction especially for ` > 0 is
itesimally-shifted
Gaussian
lobebasis
expansion

無限小変位ガウスローブ (Infinitesimally-Shifted
Gaussian Lobe) 法
2n`m (ri )Nn` ri` e ⌫n ri Y`m (r̂i )=xpression of the Gaussian basis function especially for ` > 0 isKX21⌫n (ri "i D`m,k )`⌫n ri2C`m,k e.
Y`m (r̂i ) = Nn` limn`m (ri ) = Nn` ri e`"!0 (⌫n "i )
k=1
K
2
1
11 X1
⌫n (ri "i D`m,k )
+εN
lim
Ce.
(1.9)
ε=0.3
ε ==1.0ε=0.6`m,kWe n`should
calculate theintegral for s-type GTFs.
"!0 (⌫ " )`pzn i9

k=1

lculate the integral for s-type GTFs.
01.3

0
Integral of0 the radial function

ISGL expansion

Integralin Jacobi coordinates
ral of1.3.1
the radialfunction-1−ε-1-1 a and b is written in terms oA
Jacobicoordinatetransformation
betweenchannels-4-2
024-4-2
024-4-2
02
4
ral inJacobi
coordinatesb-channel expression aspz 型関数の例
2 3 232 3 232 3a1channels
a1a1
ate transformation
betweenaandbiswritten
interms
ofthera1
Jb1 Jb2 Jb32 · · ·1 rb1ab,1 2 ab,1 ✓ab,1 · · · 2 rb1= z exp(−αzε) ) − exp (✓−α(z + ε)7)6a2a2a2) ∼ 7 6[exp7(−α(zsion as6rψa2ℓ=1,m=07 6(z)6 −ab,27]
JJJrrb2ab,2ab,2b2b2b3
6 7 6 b1
76 7 676 74πα(16ra3372= 63J a3 2J a3 J a3
7 6rb3 73=26 ab,37 6rb3 7 .
23✓ab,3ab,3b2b3
a1a1
a1 4
54 · ·5· 454 5Jb1Jb2
Jb3· ·.· 5 rb14 b1✓rab,1ab,1ab,1b1..
..
..
..
.
.
.
.
.
a2a2
6J a2 微小変位ローブ関数をそのまま積分すると桁落ちするが、
7 6rb2 7 . 6 ab,2 ab,2 ✓. ab,2 . 7 6rb2 7.
.
.
.
JJb2b36 b1
76 7 676 7(1.10)
6J a3 J a3 J a3
7 6rb3 7 = 6 ab,3 ab,3 ✓ab,37 6rb3 7 .
b2b3
4 b1 級数展開を応用することにより解析的な積分が可能．
54 5 454 5..
..
..
..
..
.. a.
.