物性物理におけるモンテカルロ法とその大規模化

藤堂眞治
東京大学物性研究所 計算物質科学研究センター
wistaria@issp.u-tokyo.ac.jp

powered by ALPS
http://alps.comp-phys.org/

1

Acknowledgment• 松尾春彦 (現 RIST), 設楽秀之 (富士通)• 諏訪秀麿 (現 Boston University), 本山裕一, 安田真也 (東大院工)• 川島直輝, 富田祐介, 渡辺宙志, Jie Lou (現 復旦大学), 正木晶子 (東大物性研),

鈴木隆史 (兵庫県立大), 原田健自 (京大), 坂倉耕太 (NEC)• Matthias Troyer (ETH), 五十嵐亮 (物性研), 他 (ALPS Collaboration)• 次世代ナノ統合シミュレーションソフトウェアの研究開発• HPCI戦略プログラム(SPIRE) 分野2 <新物質・エネルギー創成> (CMSI)• 理研計算科学研究機構(AICS) 「京」試験利用

2

CMSIは３つのコミュニティーの融合組織
分野２「新物質・エネルギー創成」

計算物質科学イニシアティブ
Computational Material Science Initiative (CMSI)
戦略機関 （代表 物性研）

分子科学研究所

物性研究所

金属材料研究所

協力機関（１１機関）

産官学連携（２機関）

産 物
総 材
研 機
構

分子科学

教育（９機関）豊 総 名 神 京
橋 研 古 戸 都
技 大 屋 大 大
科
大 学 学
大
学
物性科学

金
沢
大
学

大
阪
大
学

東
京
大
学

東
北
大
学

材料科学
7
3

%(  
 
%4=
+^8

+`8


 9 5FK:D
SWMRIXZ\PRY\+_8!   NUYO\0I 'GHL#6/I.

TVRQ\YI [; -I7G
/&$A/ 

+]8



 
2, )* G<BJ!CE"a@@@@9 [9 /I1[ [32

>?

4

応用
（奔流）

デバイス

エネルギー変換

分子機能

重点；ナノ構造電子機能予測

重点；燃料電池物質

重点；全原子Simによるウイルス分子科学

高性能リチウムイオン電池
電子的～機械的性質のマルチ・シミュレーション
スピントロニクス、
マルチフェロイックス材料

重点；水素・メタンハイドレート
生成、安定

ナノ・生体系反応制御・
化学反応ダイナミクス

太陽電池高効率長寿命化

生体分子高次構造と機能

バイオマス酸素反応解析

光機能分子と非線形外場応答
分子の光物性

光誘起電子ダイナミクスと
光・電子機能量子デバイス

新材料探索

高効率エネルギー
変換構造用材料

ナノ構造体材料の高効率非平衡エネルギー変換

重点；相関の強い量子系新量子相探求とダイナミックス

電子相間の強い物質の新機構解明強相関電子系の励起ダイナミックス

新しい量子相・量子臨界現象

基礎
12
（源流）

ポリモルフ生起分子集団機能

重点；電子状態・動力学・熱揺らぎ融和と物質理論新凝縮分子科学系揺らぎとダイナミックス
分子における電子の動的過程と多体量子動力学

分子の微細量子構造予測

新量子相・新物質5Why quantum lattice models?
• QLM: quantum spin models, bosonic Hubbard models, etc• Effects of strong correlations in multi-degree-of-freedom systems• various types of long-range order• quantum disordered phase!
J xy ! + −− +z
H=(Si Sj + Si Sj ) + J
Siz Sjz
2!i,j"(quantum liquids, spin gap phases)• phase transitions and critical phenomena• 1st order transitions, continuous phasetransitions, quantum critical point• Universality• depends on few parameters: dimensionality,
symmetry of order parameter, etc• Powerful unbiased simulation algorithms• quantum Monte Carlo methods• exact diagonalization, DMRG, etc

!i,j"

6

Deconfined critical phenomena• Possibility of continuous phase transition between two symmetry broken

phasesNovel critical phenomena due to quantum interference
Neel$state$

VBS$state$

Emergence of U(1) symmetry

Senthil, et al (2004)

broken&SU(2)&rota+onal&&&&&&symmetry&not&broken&

not&broken&

la2ce&rota+onal&&&&&symmetry&broken&• SU(2) symmetric NCCP1 model (Kukulov et al 2008)

weak 1st order?
• SU(N) J-Q model (Lou et al 2009)

consistent with DCP scenario7

Spin ladder material Na2Fe2(C2O4)3(H2O)2
J. Phys. Soc. Jpn., Vol. 78, No. 12

(b)

(a)

• Fe2+ ions in octahedral crystal fieldOeffective S=1 spins at low TbC• Fitting experimental data by QMC results

for several theoretical models (ladders,dimers, etc)

J. Phys. Soc. Jpn., Vol. 78, No. 12

a
c

(c)

Fe

a
b

H. YAMAGUCHI et al.

H. YAMAGUCHI et al.
cSusceptibility (emu/mol)

Susceptibility (emu/mol)/ion)

QMC and the DMRG methods, respectively. For
Fig. 1. (Color online) (a) Crystal structure of Na2 Fe2 (C2 O4 )3 (H2 O)2 .
Fig. 2.
(a) z
(a)
H
=
100
Oe
2þ
these
calculations,
we
use
the
followingparameter
values,ions areHydrogen and sodium atoms are omitted for clarity. Fe
0.20
(C2 O
1288.1 GHz
connected by oxalate molecules and form two-leg ladders along the
gz ¼ 4:7, gx ¼ 2:0, gy ¼ 2:5, !a ¼ 0:30, and R ¼ 0:1{ 0:5 for
0.20
H//a FeO6 octahedra
a-axis, as shown by red broken lines. (b) Distorted"
977.2 GHzH k a (" ¼ 24" , ’ ¼ 90" ) and H ? a (" ¼ 102" , ’ ¼ H//
45
).
a
y
bridged by oxalate ion. (c) Arrangement of the ladders in the bc-plane,
0.15
847.0 GHzThe principal axis directions and their g-values are deterand red broken lines corresponding to rungs of the ladders.
0.15
xmined from the saturation magnetizations and the slopes of
730.5 GHzthe ESR resonance0.10 modes in the following isolated dimer
2.2
Experimental methods
model analysis. In these field directions, the two sites of the
0.10
H⊥ a
The magnetic susceptibilities of a single crystal of SIO
H//a ladder are magnetically equivalent. In the calculation, we
0.05
were measured at 100 Oe with a superconducting
quantumH⊥aaddthetemperature-independent
Van Vleck paramagnetism
T=1.3 K
experimentmagnetometer (Quantum
interference
device (SQUID)
0.05
H = 100 results
Oe
do notexperiment
reproduce
to both #k and #? . The calculated
calculation magnetizations in pulsed
Design MPMS-XL7). High-fieldtheexperimental
ones
well,
especiallytheupturn
of
#
k
0
50
0.00
magnetic fields up to about 53 T were measured utilizing a
15 Van20 Vleck25 para30
at approximately 5 K. In5 Fig. 106(b), the5
10
15
20 the 25
30
non-destructive
pulse
magnet,andmagnetizations
were
Temperature (K)manetic contributions are subtracted from both Mk and M? .
Temperature(K) using a standard pick-up
detected by an induction
method)3 (H2 O)2 at 1.3 K for(b)
Wesignals,
find a small anomaly at half of the saturation magnetng and weak
coil system.
z
H// aYamaguchi,
Kimura,Honda,
Okunishi, Todo, Kindo, Hagiwara(2009)(b)ization of R ¼ 0:14 calculated for H k a, which can be
ESRmeasurements on a single crystalH//aof SIO in pulsed
4
considered to originate
fromy a 1/2 magnetization plateau.
magnetic fields up to about 53 T were carried out using our Fig. 3.
24°
andThe calculated results
also poorly explain the experimental
8 H
pulse field ESR apparatus equipped with a non-destructive
3
B/ion)

the1623.6 GHzcb

YAMAMOTO,
ated by SSE. System sizes
are from LTODO,= 8 to AND
L = 16.MIYASHITA
denote linear extrapolations of data points.

0.3

1.4

PHYS

B.1.2
Mean-fi
MOTO, TODO, AND MIYASHITAPHYSICAL REVIEW B 79, 094503 $1
We also analyze
theTherefore we confirm thethe su- Bose-Hubbard Hamiltonian
We existence
analyze theofextended
on a
0.2
SUCCESSIVE PHASE TRANSITIONSAT FINITE… 22 In orderapproximation.
!
"
1
/
2.
se in the region of
cubic
latticeII. MODEL
B. Mean-field approximation 0.8
sublattice structure which
0.15 0.25
Sπ(L=12)by theSSE
We nalso
theNLordering
processes
derTTSofdensity.0.6by
Herethew
a j +a aia†j % + V= − t $a†i on
analyze the extended Bose-Hubbard H
Hamiltonianπ
in j analyzebySSEρ
• Interactings
TRANSITIONSAT FINITETEMPERATURES
soft-corebosons "ij#π(L=10)"ij#
0.1 220.2SIn
order to orderstudy
the superfluid
solid state,orderweapproximation.
andattice0.4
ρs(L=12)sublattice structure
which is characterizeda staggere
tonian for thisbyMF
is givenρ(L=10)
0.15
1
study the ordered† states †at finite temperatures.
sSF− !of then0.05
$1%NSwe adopt mean
fields0.2
for theH = − t $ai a j + aia j % + V nin j + U ni$ni − 1%der
i , density. HereII. MODEL

0.25

Sπ

!

ρs

Supersolid in extended Bose-Hubbard model
!

ρs

Sπ

T/U

!ithe
!i
!
!
ansition between
thenormalliquid
phase
and
2
=
H
0.1
MF H
"ij#
"ij#
order and superfluidorderatsublattices
A
and
B.
The
0
0
SS
s expected to belong to the universality
class
of
†
tonianfor thisoperators
MF is0.1
given
by
0.05
0
0.2
0.3
0.4
1
where
a
and
a
are
the
creationandannihilationofiiodel•because
the
phase
transition
occursinthe†+ U ! ni$n−boson1% − !$&anai†,' = "of% and n$1%
Supersolid=ai co-existence
!
†
(b)
,
=
a
T/U+ HA +i jij with i=−
zt$ai ai. The parameter t 0denotes
2 i of the orderi parameter
symmetry breaking=
HC,H
A + a A% & B +0.04 MF
0.045 A 0.05 B 0.055
0.06
the
hoppingmatrixelement, Uandand V are the on-site andlong-rangeorder(=solid)On the otherhand,
the phasetransitiony† of Z2. diagonal
t/U
!
is
the
nearest-neighbor
repulsions,
respectively,
and
0.3
1.4
aioff-diagonal
are thetocreationandannihilation
operatorsoflong-rangeorder(=superfluid)dai isandexpected
belongtothe
XY
universality
†
FIG.
6.sum
The
t-T phase
diagram for V / UU
= 1 / z and ! / U = 0.7† overchemicalpotential.
The
notation
"ij#
means
the
the
a†j ' = "parameter
%
and
n
=
a
a
.
The
parametertdenotesne the$&ai ,order
ij
i has
i
=
−
zt$a
+
a
%
&
+
zVn
m
+
1% T−S !nH
i theS(L=12)
A 0.25B
A$n
A −state
1.2
symmetryopticalof U!1". lattice
A TheAtransition
π A B
obtainedbytemperatures
of n
the
solid
3 SSE.
•
Experimental
realization:
"
2
,
where
L
is
Theandsystem size areis N=Lpping matrix element, nearest-neighbor
U and V are thepairs.
on-siteS(L=10)
plotted by solid circles and thoseT#s areπ of the superfluid staten, we study the temperaturedependenceof theseThe order parameter†ρ(L=12)
1thelengthofthesystem.
of
the
solid
= −datazt$a&A +HB the-neighbor repulsions, respectively, and ! is theplotted by open circles. The liness connect
points
for
B%the
B+a
0.2
eters.
ρs(L=10)
guide to the eyes.
state
al potential. The notation"ij#ismeans the sum over the
0.8
3-neighbor pairs. The system size is N = L , where L is
0.15
U solid for0.6!t / U
1
successive transitionsof superfluid and
†
iQ·$r
−r
%
j
k
gth A.
of Stochastic
the system.
The expansion
order parameterS#of= the2 !
solidseries+ aBAs%&zVnBinmthe
− !n
e
"nH
#,= − zt$a$2%
A + seen
A +SSE n
B$nB − 1%here
jn=kB0.045,!/UB= 0.7".
wassimulation,
2
N jk
0.1
C =su2zt&
we find again the suppression of the solid order by 0.4
theshow the results obtained by SSE. The simulaperfluid fraction. Namely, S" has a cusp at the superfluid
0.05
#
,
#
,
#
%
is
the
wave
vectorthat
represents
thedepictwhere
where
Q
=
$
z$=6%
is andthenumbeerformed in the1grandiQ·$r
canonical
ensemble using
0.2comtransition
point. We also
the phase
diagram−r%j k "n n #,e3.staggered$2% parameter
As for the ordersuperfluid!
# = N2=
j order.
kpare ofthattheto thatof SSE !Fig. 8".
a qualitatively
mATheyandshowmB arethe mean fi
12
es N = 103Sand
N jk
0
0!tc ,numbe&Ais&tation
zVm
Ce.g.,= 2zt
T c"
good fraction
agreement,
there
the
we adopt
Bose-Einstein
B − tetracritical
AmB ,pointvaluesofthethat of condensation
the order parameters state,of solid,S#, and
0 temperatures
0.1
0.2
0.4
and the critical
TS" and T0.3
#s are suppressed byspectively,, # , #% is ofthetemperature
wave vector
represents
= $#
, as
a function
forthatvariousvaluesz$=6% isof thenumberofnearest-neighborsites.
sQ
appearance
the
other
order
as
before.
(c)1 the † where
T/Umentioned†
http://www.uibk.ac.at/th-physik/qo$k=0superfluid= 2 ! "a j akm+Aakand
a j#Let$3% fields
edorder.
forthethe transition
order parameterof the
5!a",
we As
show
from normal
liquidstudythethe meancompetition
betweenthe solid and tosupermBus arecorresponding
the
e
=
"n
#
m
N
A
A
ordersofbythe
analyzingtheGinzburg
LandauYamamoto,
Todo,
Miyashita(2009)
wecondensation
fraction con- jk
tationnumberoperatorsfor
A !GL"
and freeB sitelidadopt
for t /Bose-Einstein
U = 0.02 in which
only S# appears
FIG. fluid
5.values
Temperaturedependenceof orderparameters
obt
energy. Since the order parameters
m
and
$
take
small
value
Similarly,&normalA and &
B correthe same way, the transition
from normal
liquid whilebyinspectively,
SSE!V= Uvicinity/simulationz". !a"of !tthe/ Utetracritical= 0.02, $/ U = 0.7":liquid-nointhemean-field
analysis,
the
SSE
we
in
the
point,
theGL
freeenergyis 91

多体量子系のためのアルゴリズム• 最近20年の進展• 厳密対角化、Lanczos法• 補助場量子モンテカルロ法 Scalapino, Hirsh, Sorella, Imada• 密度行列繰り込み群法 White• 動的平均場近似とその拡張 Kuramoto, Metzner, Vollhardt, Kotliar• 変分モンテカルロ法 Gros, Yokoyama, Shiba, Sorella, Ogata• 経路積分繰り込み群法 Kashima, Imada• ガウス基底モンテカルロ法 Gorney, Drummond, Assaad, Troyer, Imada• 第一原理計算との融合 Kotliar, Georges, Aryasetiawan, ...
• クラスターアルゴリズム量子モンテカルロ法 Evertz, Kawashima, Todo, ...

10World-line QMC for lattice models• Progress since 1990s• Direct simulation in imaginary time continuum• Alternative representation based on high-temperature expansion (SSE)• Simulation in grand-canonical ensemble• Solving critical slowing down by cluster updates (loop algorithm, etc)
• Worm (directed-loop) updates for systems with less symmetries• Large spin size S ≥ 1, bosons• Off-diagonal correlation functions• Solving negative sign problem (only for very special models)• Extended ensemble Monte Carlo method (quantum Wang-Landau)• Order-N method for long-range interaction• Non-trivial parallelization of loops using MPI and OpenMp• Using balance-condition instead of DBC

11

「量子モンテカルロ」とは何か?
tr[A exp( H)]A =tr[exp( H)]• 次の式を評価したい• ただし、H も A も超巨大な行列 (次元 〜 10106 )• テイラー展開の第 n 次の項を考える

tr H =
n

i1

i2

···

in

Hi1 ,i2 Hi2 ,i3 · · · Hin ,i1

106 n 個の項を確率的にサンプリング• この中に現われる (10)• それぞれの項は、ある種の「経路」(path)を表す
5

4

3

6

1

1

4

3
5
6
12

連続虚時間ループアルゴリズム量子モンテカルロ法(Continuous imaginary-time loop algorithm quantum Monte Carlo method)• 虚時間経路積分により d+1 次元古典系に焼き直す (世界線表示)• クラスター(ループ)アルゴリズムによる非局所更新 - 緩和は非常に速い
o=`

o=0

• 実際にプログラムの内部で行っていること• 指数分布に従って、一定密度で「横木」を生成• Union-Find アルゴリズムを使って、全てのクラスター(ループ)を認識
13

One Monte Carlo step in loop algorithm

cluster
identification

labeling

world-line configuration

graphloop flip

loops

next configuration

parallelization of loop algorithm is nontrivial!
• insertion of operators/cuts with constant density• cluster identification

Poisson process

non-local operations on trees/lists• cf) conventional local flip is easy to parallelize, but is suffered fromcritical slowing down ~(system size)2
14

Parallelization of loop algorithm• distributing world-line configuration to several processors• (d+1)-dimensional lattice is split into Np slices with equal

imaginary time thickness• labeling process can be performed independently on eachprocess15

Parallelization of loop algorithm• binary-tree algorithm for cluster identification

N ⇥• overhead of parallelization is (N log Np )/= Np log Np /Np(negligible at very low temperatures)

16

Further Optimization(a)• Introduction of butterfly-type communication• eliminates overhead in distribute process• majority vote trick to determine loop

node 7slave

node 6masternode 5slave

node 4

masternode 3slave

node 2masternode 7slave

node 6

node 5master

slave

slave

node 0masternode 4

node 3slave

node 2

directions with minimum communication• efficient for high-dimensional torus• Reduction of number of stages• ‘bucket brigade’ adjacent communication• effective for low-dimensional torus• Hybrid parallelization with MPI and OpenMP• multi-thread parallelization with Open MP• inter-node parallelization with MPI• reduces memory and communication overhead in each node
• fine-grained parallelization of cluster identification by introducing
asynchronous wait-free union-find algorithm
node 1

masterunify

updateunify

stage 1

node 1

masterbroadcast

updateunify

stage 2

node 0broadcastrank 7

rank 6

rank 6

rank 5rank 5rank 4rank 4rank 3rank 3rank 2

rank 2

rank 1

rank 1

rank 0rank 0

updatestage 0

unify

updatestage 1

unify

updatestage 3rank 7unify(b)updatestage 217

Asynchronous wait-free union-find algorithm
1.find root of each cluster/tree
2.unify two clusters
3.compress path to the new root

• locking whole clusters (e.g. by using“omp critical”) is too expensive• finding root and path compression are“thread-safe”• wait-free unification by using CAS(compare-and-swap) atomic operation8
7

speed up

6
5
4
3
2
1
1

2

3

4

5

6

7

8

number of threads
18

並列化チューニング状況• 京におけるチューニング状況

ノード数
24576
通信時間
13.7%
並列化オーバーヘッド 54.4%

100000

10000

82944
15%
58%

9

1000

8
7
wall time [sec/MCS]

speed up

京フルシステムでの並列性能推定

100

10

6
5
4
3
2

1

1

1

10

100
1000 10000 100000
number of nodes

0

1

10

100
1000 10000 100000
number of nodes• ウィークスケーリング：24576ノード:12288ノード = 98.66%• MIPS値：0.164 PIPS (理論ピーク性能の10.4%)• FLOPS値：7.63 TFLOPS (理論ピーク性能の 0.243%)
19

ALPS/looper による「京」全系計算• マイルストーン：S=5 反強磁性ハイゼンベルグ鎖の量子モンテカルロ計算→ 相関長 〜 2×105、エネルギーギャップ 〜 5 × 10-5
→ 膨大な数の頂点からなるランダムグラフ中の全ての連結成分の識別問題• 頂点・辺の数 〜 1013、クラスター数 〜 1012
• 線形演算ではなく、並列グラフ処理により量子統計物理の重要問題を解く• Graph500 の “Large” 〜 “Huge” 問題サイズ(メモリ 0.14〜1.1PB)に匹敵• 2011/11の Graph500 では最大でも“Medium” サイズ(メモリ 17TB)• TOP500の一位の「京」フルノードでしかできない非自明並列化• 単なるベンチマークだけでなく、フルノードを使ったプロダクション

ラン(平衡状態のシミュレーション)• ハイブリッド並列化による高速化 〜 105• ループアルゴリズムによる高速化 〜L2〜 1014

実質的には 1019 の高速化
20

The ALPS project
ALPS = Algorithms and Libraries for Physics Simulations• International collaboration for developing open-source softwares for

simulation of quantum lattice models, such as quantum spin systems,
electron systems, etc• ALPS Libraries = collection of genericParameter XML Files

Lattice XML File

Model XML FileC++ libraries• ALPS Applications = collection of

Quantum Lattice Modelapplication packages using modern
algorithms such as QMC, DMRG, ED, etc

Quantum Monte Carlo

Exact Diagonalization

DMRG• ALPS Framework = environment for

executing large-scale parallel simulations
including XML schemas, tools, scheduler, etc

Outputs in XML Format

http://alps.comp-phys.org/
21

ALPSプロジェクトの目標• 計算物性物理の現状• 研究グループ毎にそれぞれ異ったコード• シミュレーションを行なう模型毎に異ったコード• アルゴリズムはますます複雑に• 互換性のない入出力形式• ALPSプロジェクトの目標• 最新のアルゴリズムを用いた ”community code” の開発• 大規模並列計算などのためのC++ライブラリ・フレームワーク開発• 統一入出力形式の提案とそれにもとづくデータ解析ツールの作成• 計算物理の専門家だけでなく、理論家・実験家にも使えるシミュレー

ションソフトウェア

22

A Simulation using ALPS
Lattice XML Filesquare lattice

Model XML File
H = J

Y [ S iz S jz + ( S i+ S j- + S i- S j+ )/ 2 ] - h Yi S iz

< i,j >

Parameter File
LATTICE = "square lattice"
MODEL= "spin"
L
= 16
Jxy
= 1
Jz
= 2
SWEEPS = 10000
THERMALIZATION = 1000
{ T = 0.1 }{ T = 0.2 }{ T = 0.5 }{ T = 1.0 }parameter2xml
tool<LATTICES><LATTICE name="square lattice" dimension="2"><PARAMETER name="a" default="1"/><BASIS><VECTOR>a 0</VECTOR><VECTOR>0 a</VECTOR></BASIS></LATTICE><UNITCELL name="simple2d" dimension="2"><VERTEX/><EDGE><SOURCE vertex="1" offset="0 0"/><TARGET vertex="1" offset="0 1"/></EDGE><EDGE><SOURCE vertex="1" offset="0 0"/><TARGET vertex="1" offset="1 0"/></EDGE></UNITCELL><LATTICEGRAPH name="square lattice"><FINITELATTICE><LATTICE ref="square lattice"/>(0,1)<EXTENT dimension="1" size="L"/><EXTENT dimension="2" size="L"/><BOUNDARY type="periodic"/></FINITELATTICE><UNITCELL ref="simple2d"/>
(0,0)
(1,0)
</LATTICEGRAPH></LATTICES>

Parameter XML Fileapplication programs<MODELS><BASIS name="spin"><SITEBASIS name="spin"><PARAMETER name="local_S" default="1/2"/><QUANTUMNUMBER name="S" min="local_S" max="local_S"/><QUANTUMNUMBER name="Sz" min="-S" max="S"/><OPERATOR name="Sz" matrixelement="Sz"/><OPERATOR name="Splus" matrixelement="sqrt(S*(S+1)-Sz*(Sz+1))"><CHANGE quantumnumber="Sz" change="1"/></OPERATOR><OPERATOR name="Sminus" matrixelement="sqrt(S*(S+1)-Sz*(Sz-1))"><CHANGE quantumnumber="Sz" change="-1"/></OPERATOR></SITEBASIS></BASIS><HAMILTONIAN name="spin"><PARAMETER name="J" default="1"/><PARAMETER name="h" default="0"/><BASIS ref="spin"/><SITETERM>-h * Sz(i)</SITETERM><BONDTERM source="i" target="j">J * (Sz(i)*Sz(j) + (Splus(i)*Sminus(j)+Sminus(i)*Splus(j)))/2</BONDTERM></HAMILTONIAN></MODELS>

Quantum Lattice Model

Quantum Monte Carlo

Exact Diagonalization

DMRG

DMFT

Python based evaluation tools
Plots

Outputs in HDF5 & XML

23

ALPS Lattice XML
periodic chain with length L<LATTICE name="chain lattice" dimension="1"><BASIS><VECTOR> 1 </VECTOR></BASIS></LATTICE><UNITCELL name="simple1d" dimension="1" vertices=”1”><EDGE><SOURCE vertex="1" offset="0"/><TARGET vertex="1" offset="1"/></EDGE></UNITCELL><LATTICEGRAPH name = "chain lattice"><FINITELATTICE><LATTICE ref="chain lattice"/><PARAMETER name=”L”/><EXTENT size ="L"/><BOUNDARY type="periodic"/></FINITELATTICE><UNITCELL ref="simple1d"/></LATTICEGRAPH>

24

Model XML for describing Hamiltonian
XXZ spin model with two types of bonds (e.g. nearest and next nearest neighbor interactions)<HAMILTONIAN name=”spin”><PARAMETER name=”Jz” default=”J”/><PARAMETER name=”Jxy” default=”J”/><PARAMETER name=”J” default=”1”/><PARAMETER name=”Jz'” default=”J'”/><PARAMETER name=”Jxy'” default=”J'”/><PARAMETER name=”J'” default=”0”/><PARAMETER name=”h” default=”0”/><BASIS ref=”spin”/><SITETERM site=”i”> -h * Sz(i) </SITETERM><BONDTERM type=”0” source=”i” target=”j”>
Jz * Sz(i) * Sz(j) + Jxy * (Splus(x)*Sminus(y)+Sminus(x)*Splus(y)) / 2</BONDTERM><BONDTERM type=”1” source=”i” target=”j”>
Jz' * Sz(i) * Sz(j) + Jxy' * (Splus(x)*Sminus(y)+Sminus(x)*Splus(y)) / 2</BONDTERM></HAMILTONIAN>

H=

!

!i,j"

!
J
xy
(Si+ Sj− + Si− Sj+ )] −[Jz Siz Sjz +
hSiz
2
i
25

ALPS framework and “K Computer”
ALPS/looper library

tools
applications
domain-specificlibraries
numerics•

parallelization within MCS

XML	 manipulation	
 	
 	
 	
 	
 GUI	
 	
 	
 	
 	
 database	
 	
 
MC	
 	
 	
 	
 	
 QMC	
 	
 	
 	
 	
 	
 ED	
 	
 	
 	
 	
 	
 DMRG
lattice	 	 	
 model	 	 	
 observables	 	 	
 scheduler
random	 	 	 ublasiterative	
 eigenvalue	
 solver

Boost	
 library
generic	
 C++
C	
 /	
 Fortrangraph

ALPS/parapack library
multiple-level
parallelization• adaptors for extended
ensembles•

serialization	 	 	
 XML/XSLT	 BLAS	
 	 	 	 	
 LAPACK	
 	
 	
 	
 	
 MPI

26

ALPSの展開• 上方展開 (大規模化・高性能化・並列化)• 中核アプリ「大規模並列量子モンテカルロ法」(ALPS/looper) の高度化
• 高並列スケジューラ(ALPS/parapack)のハイブリッド多重並列化• Fortran, C プログラムのための API 作成• 下方展開• 実験家・理論家による幅広い利用を視野に• Windows/Macintosh 用バイナリインストーラの開発• ワークフロー・履歴管理システムとの統合• GUI (グラフィカルユーザインターフェース)の開発• マニュアル・チュートリアルの日本語化
27

Windows, Mac OS X 用バイナリインストーラ

for Mac OS X

for Windows28

ワークフロー・履歴管理システムとの統合

http://www.vistrails.org/

29

ALPS version 2.1 & ALPS/looper version 4.0• JSTAT, P05001 (2011)
The ALPS project release 2.0:
open source software for strongly correlated systems
B. Bauer1 L. D. Carr2 H.G. Evertz3 A. Feiguin4 J. Freire5
S. Fuchs6 L. Gamper1 J. Gukelberger1 E. Gull7 S. Guertler8
A. Hehn1 R. Igarashi9,10 S. Isakov1 D. Koop5 P.N. Ma1
P. Mates1,5 H. Matsuo17 O. Parcollet12 G. Pawlowski13
J.D. Picon14 L. Pollet11,1 T. Pruschke6 E. Santos5
V.W. Scarola15 U. Schollwöck16 C. Silva5 B. Surer1 S. Todo17,10
S. Trebst18 M. Troyer1 ‡ M. L. Wall2 P. Werner1 S. Wessel19,20

Theoretische Physik, ETH Zurich, 8093 Zurich, Switzerland
Department of Physics, Colorado School of Mines, Golden, CO 80401, USA
3
Institut für Theoretische Physik, Technische Universität Graz, A-8010 Graz, Austria
4
Department of Physics and Astronomy, University of Wyoming, Laramie, Wyoming
82071, USA
5
Scientific Computing and Imaging Institute, University of Utah, Salt Lake City,
Utah 84112, USA
6
Institut für Theoretische Physik, Georg-August-Universität Göttingen, Göttingen,Germany
7
Columbia University, New York, NY 10027, USA
8
Bethe Center for Theoretical Phyics, Universität Bonn, Bonn, Germany
9Center for Computational Science & e-Systems, Japan Atomic Energy Agency,
110-0015 Tokyo, Japan
10
Core Research for Evolutional Science and Technology, Japan Science andTechnology Agency, 332-0012 Kawaguchi, Japan
11
Physics Department, Harvard University, Cambridge 02138, Massachusetts, USA
12
Institut de Physique Theorique, CEA/DSM/IPhT-CNRS/URA 2306, CEA-Saclay,
F-91191 Gif-sur-Yvette, France
13
Faculty of Physics, A. Mickiewicz University, ul. Umultowska 85, 61-614 Poznan,
Poland14
Institute of Theoretical Physics, EPF Lausanne, CH-1015 Lausanne, Switzerland
15
Department of Physics, Virginia Tech, Blacksburg, Virginia 24061, USA
16
Arnold Sommerfeld Center for Theoretical Physics and Center for NanoScience,
University of Munich, Theresienstrasse 37, 80333 Munich, Germany
17
Department of Applied Physics, University of Tokyo, 113-8656 Tokyo, Japan18
Microsoft Research, Station Q, University of California, Santa Barbara, CA 93106,
USA
19
Institute for Solid State Theory, RWTH Aachen University, 52056 Aachen,
Germany
20
Institut für Theoretische Physik III, Universität Stuttgart, Pfaffenwaldring 57,
70550 Stuttgart, Germany
12•

http://alps.comp-phys.org/•

http://wistaria.comp-phys.org/alps-looper/‡ Corresponding author

30

応用
（奔流）

デバイス

エネルギー変換

分子機能

重点；ナノ構造電子機能予測

重点；燃料電池物質

重点；全原子Simによるウイルス分子科学

高性能リチウムイオン電池
電子的～機械的性質のマルチ・シミュレーション
スピントロニクス、
マルチフェロイックス材料

重点；水素・メタンハイドレート
生成、安定

ナノ・生体系反応制御・
化学反応ダイナミクス

太陽電池高効率長寿命化

生体分子高次構造と機能

バイオマス酸素反応解析

光機能分子と非線形外場応答
分子の光物性

光誘起電子ダイナミクスと
光・電子機能量子デバイス

新材料探索

高効率エネルギー
変換構造用材料

ナノ構造体材料の高効率非平衡エネルギー変換

重点；相関の強い量子系新量子相探求とダイナミックス

電子相間の強い物質の新機構解明強相関電子系の励起ダイナミックス

新しい量子相・量子臨界現象

基礎
12
（源流）

ポリモルフ生起分子集団機能

重点；電子状態・動力学・熱揺らぎ融和と物質理論新凝縮分子科学系揺らぎとダイナミックス
分子における電子の動的過程と多体量子動力学

分子の微細量子構造予測

新量子相・新物質
31

CMSIで開発中のアプリケーション (41本)• 第1部会• MACE, DDMRG, 2D-Extended-DMRG, ALPS/looper, GELLAN, MC-

MOZ, GAMESS拡張• 第2部会• RSDFT, CONQUEST, RS-CPMD, RSPACE, Multi-probe electron transport

simulator, hybridQMCL, PIQUANDY, CPVO, HiLAPW, OpenMX, TC++,
CASINO, QMAS, xTAPP, FEMTECK• 第3部会• modylas, GAMESS-FMO, REM, PIMD, ermod, DC, Direct SAC-CI• 第4部会• STATE, TOMBO, feram, CP2K, MoIDS, CPMD, BigDFT, 3D-RISM,

Amber9拡張, Machikaneyama2002, HiRUNE• 他• hnpack, NANIWA

赤: ナノGC中核アプリ(6), 太字: 京で試験利用中(29)32物質科学分野の代表的なアルゴリズム• 密度汎関数法 (平面波基底、実空間基底)、第一原理分子動力学法、時間依存密度汎関数法、 実空間・実時間電子電磁場ダイナミクス法、分子動力学 (短距離相互作用、長距離相互作用)、モンテカルロ法 (古典、量子)、量子化学計算、厳密対角化• 物質科学のアプリケーションの特徴• 相関の強い系の平衡状態・定常状態・アンサンブルに興味がある (1トラジェ

クトリの計算だけではだめ)• シミュレーションする系の次元は3次元に限らない。量子的な相関など、空間的に局所的であるとも限らない• 遠くの相関を取り込む/平衡状態にたどりつくための工夫として、非局所的な操

作/演算が多くの場合必要• メモリバンド幅を要求するアプリケーションが多い• 大域通信・バタフライ型通信など高次元のネットワークを要求するものも多い
33

「サイエンス・モデル」と開発環境• 計算物質科学では、少数のコミュニティーコードではなく、むしろ研究グ

ループ毎に多数のコードが存在 → 新しい手法・アルゴリズムがどんどん試されている• 「スモールサイエンス」的な要素も強い• 現象から「有効模型」の抽出 〜 様々なアイデアにもとづいたアルゴリ

ズムの開発 〜 プログラミング 〜 シミュレーションの実行 〜 結果の解
析〜「有効模型」へのフィードバックをごく少人数で• 「ギリギリのチューニングをしてxx%」ではなく「3ヶ月でyy%の実効効率」が出せるような開発環境が重要• PC 〜 クラスタ 〜 スパコン 〜 京 〜 次世代へのシームレスな環境が必要• 「京」のような大規模なスパコンを使いこなすためには、ソフトウェア開発においても個人を越えて組織的に進めることが今後不可欠に34

Summary & Outlook• HPCI戦略プログラム分野2 <新物質・エネルギー創成> CMSI
• 強相関量子多体系研究のフロンティア• 新量子相と量子臨界相の探求・解明• 量子モンテカルロの概要・並列化・大規模化• 整数演算・グラフアルゴリズムの重要性• 計算機の性能向上 + アルゴリズムのブレークスルーによる相乗効果• ALPSプロジェクト• 物質科学分野のアプリケーション

35