Transcript ppt - KEK

理論化学の現状と将来の
可能性について
平尾 公彦
原子分子の基礎方程式
Schrödingerの波動方程式から80年！
N個の原子核とn個の電子から成る
系のSchrödinger方程式は
H  E
n
Z Z
2
i2
Z
1
H  
   
2  ,i ri i j rij  , R
 2M
i
N
  (r1, r2 ,, rn ; R1, R2 ,, RN )
原子核の運動－反応動力学
電子の運動－電子状態理論
原子核の運動－反応動力学
電子の運動－電子状態理論
原子分子の基礎方程式
原子核と電子の運動の分離
Born-Oppenheimer近似
n個の電子から成る系のSchrödinger方程式は
H  E
n
i2
Z
1
H      
2  ,i ri i  j rij
i
  (r1, r2 ,, rn ; R)
分子系ではハミルトニアンが厳密に定義される
Diracの予言
1929年、P. Dirac
「量子力学の出現で、化学の問題の解決に必要な基
礎的物理法則は全て分かった。困難は、ただ、これら
の法則を厳密に適用すると複雑すぎて解ける望みのな
い方程式に行きついてしまうことである」
原理的に化学が解けたにしても、化学の問題として
知りたいことはいくらでもある。
Quantum Chemistry has come a long way !
ある有限規模の計算の可能性のなかで、どうした
ら化学的精度(kcal/mol)を得るかというのが理論化
学者の問題であった。私たちは物理学者とは異な
るセンスで、しかし同じような細心さで内容を分
析し、理論を作り上げてきた。
理論化学の成功はDiracの予想をはるかにこえる
ものであったと思われる。Diracが考え及ばな
かった最大の要素はコンピュータの驚異的な発達
である。
Schrödinger (Dirac) 方程式を近似的に解く手段
変分法
摂動論
クラスター展開法
有効な数学的方法であるだけでなく、それぞれに
自然哲学的要素をもっている。
多電子系の近似理論
多電子系では、2階の微分方程式を解析的展開法
を利用して、行列式の固有値問題という代数の問
題に帰着させて近似的に解いている。波動関数/
電子密度、エネルギーは電子の軌道関数の汎関数
で表される
  i ,
E  E i 
また、軌道関数も原子軌道で展開される
   C
m
i
p
p
pi
近似
自然現象を説明する理論なるものは、所詮ひと
つの近似以上にでることはあり得ない。ある近
似の平面を選んでそこで自然を切り、その平面
への投影を見ているにすぎない。
願わくば、ものの本質が映し出される平面で自
然を切りたいものである。
分子シミュレーションの2つのアプローチ
 波動関数法（ab initio MO法）
体系的理論、近似を上げれば厳密解に収束
小さな系では極めて精度の高い結果を与える
計算時間の系の大きさ(N)依存性が大きい
（大きな系では計算時間がかかりすぎる）
 電子密度法（密度汎関数法）
簡潔で概念的、大きな系に適用可能
半定量的理論で、精度は汎関数に依存
2次の分子物性、van der Waals力や励起状態の記述が悪い
分子シミュレーションの2つのアプローチ
 波動関数法（ab initio MO法）
体系的理論、近似を上げれば厳密解に収束
小さな系では極めて精度の高い結果を与える
計算時間の系の大きさ(N)依存性が大きい
（大きな系では計算時間がかかりすぎる）
課題－効率的なアルゴリズムの開発
 電子密度法（密度汎関数法）
簡潔で概念的、大きな系に適用可能
半定量的理論、精度は汎関数に依存
2次の分子物性、van der Waals力や励起状態の記述が悪い
課題－汎関数の開発
私たちが開発してきた分子理論
密度汎関数法
波動関数理論
MRMP法, MCQDPT法
CASVB法, QCAS法, EBS法
SPS-SCF, PT
Linear Scaling 法
（GFC, Local MP2,平面波）
OP相関汎関数
Parameter-free交換汎関数
長距離交換補正(LC-GGA)
自己相互作用補正法
Van der Waals
Dual-level DFT
UTChem
相対論的分子理論
プログラム
RESC, DK3法
Dirac-HF法, DC-CASPT2法
Dirac-Kohn-Sham法
Ab initio Model Potential
相対論基底関数
ダイナミクス理論
Ab intio Dynamics
QM/MM分子動力学
量子ダイナミクス
Quantal cumulant dynamics
時間依存DFT
Energy lowering due to relativity
(a.u.)
Cu (29)
Ag (47)
14.2 au
115 au
Relativistic effects become significant
for second-row transition elements
and larger.
Au (79)
1149 au
金が黄色に輝くのも、水銀が常温で液体なのも、
Z
鉛がやわらかいのも、白金の触媒効果もみんな
Difference in atomic SCF energies as a function of
相対論効果
the nuclear charge
相対論的分子理論
4成分相対論的Hamilton演算子に基づく理論
計算負荷の高い相対論的方程式を効率的に解くための新しい計算スキーム
相対論分子積分の高速計算法
世界最速の積分アルゴリズム
相対論積分の計算時間
LLLL+LLSS+SSSS
LLLL+LLSS
LLLL
Present
1.37
0.77
0.21
MOLFDIR
76.35
21.16
1.63
Dirac-HF法, DC-CASPT2法
高速な4成分相対論分子理論
Dirac-Kohn-Sham密度汎関数法
世界最初のGauss型4成分密度汎関数理論
相対論的近似Hamilton演算子に基づく理論
重原子分子に対する高精度・大規模計算が可能
RESC法
本質を衝くモデル，simpleなアイデアにも拘わらず高精度
高次Douglas-Kroll(DK)法
2成分相対論的分子理論の決定版
60倍高速
相対論的分子理論
複雑なスペクトルの決定
磁気的性質
NMR化学シフト
光電子スペクトル
2.80
2.70
われわれの理論によってすべての元素を対象
分光学的定数の決定
化学反応の追跡
とした定量的理論計算が可能になった
触媒設計
O
C
相対論効果を考慮
しないと実験の傾
向を再現できない
M+
R
非相対論
P
2.60
相対論
2.50
実験
2.40
2.30
基底状態
励起状態
2.20
Cu2
Ag2
Au2
P
33.2
20.6
Pt
0.0
R
O
C
R
M+
P
20.9
-1.5
CO
P
+M
P
P
動力学的シミュレーション手法の開発
量子化学計算に基づく新しいダイナミクス手法
化学結合の解離・生成 / 電子励起状態
高速な動的手法の開発
→ ab initio Dynamics、QM/MM法、Quantal cumulant dynamics
ナノ秒オーダーの高精度シミュレーションを実現し、
大規模系における動的効果の解明
時間依存DFT法に基づく励起状態シミュレーション
従来の汎関数では扱えない：
• 擬縮重した励起状態
• 準安定構造
• 長距離電荷移動
→ LC汎関数により取り扱いが可能
レチナール分子の光異性化反応
ダイナミクス－機能という新しい概念を構築したい
UTChemプロジェクト - 理論の実用化
理論化学のソフトウエアは物質科学の共通基盤
ソフトウエアの研究開発は欧米に大きく遅れをとる
世界の分子軌道プログラム・パッケージ
•
•
•
•
•
•
•
•
DALTON
Gaussian (Gaussian, Inc.)
GAMESS (Iowa State Univ.)
Molpro
Molcas (Lund Univ.)
NWChem (PNNL)
Q-Chem (Q-Chem, Inc.)
DMol3(Accelrys)
“すべてをわれわれ自身の手で”作り上げた
日本が知的財産として世界に誇れる
物質科学のソフトウェア
なぜ新しいソフトウェアが必要か
(1) 既存のプログラムを使うことで多くの恩恵を受ける。しかし、プログラムの制
限で真のブレークスルーは成し難い。
(2) 既存のプログラムでは取り扱うことができない重要で強力な方法論が多く存
在する。
Ab initio MO, DFT, and Dynamics
Non-relativistic and Relativistic
(Two- and Four-component relativistic)
CISD, CISDT, CISDTQ,
CCSD, CCSDT, CCSDTQ,
MP2, MP3, and MP4,
are implemented into UTChem.
http://utchem.qcl.t.u-tokyo.ac.jp/
理論計算の信頼性
水素化物の分光学定数の計算
Main-group elements across the second- through fifth-period of the
periodic table J.Chem.Phys. 120, 3297 (2004)
BH
AlH
GaH
InH
CH
SiH
GeH
SnH
NH
PH
AsH
SbH
OH
SH
SeH
TeH
FH
ClH
BrH
IH
DK3-CCSD, DK3-CCSDT, DK3-CCSDTQ
Re-contracted relativistic cc-pVnZ (n=2-5)
Extrapolation
分光学定数 CH, InH

12CH (2)
115InH(1+)
Theory
Exp.
Theory
Exp.

re / Å
r0 / Å
Be / cm-1
B0 / cm-1
e / cm-1
De / cm-1
e / cm-1
xe / cm-1
0 / cm-1
1 / cm-1
2 / cm-1
D00 / eV
1.120
1.131
14.449
14.182
0.534
0.00148
2860
66
1413
4141
6743
3.47
1.120
1.130
14.457
14.190
0.534
0.00145
2859
63
1413
4146
6752
3.47
1.841
1.855
4.977
4.905
0.144
0.000229
1467
26
727
2142
3510
2.46
1.838
1.851
4.995
4.923
0.143
0.000223
1476
26
732
2157
3535
2.48

理論は実験値を次の精度で予測できる
結合距離
回転定数
振動－回転定数
調和振動数
非調和振動数
結合エネルギー
0.002 Å 以内
0.02 cm-1 以内
0.01 cm-1 以内
9 cm-1
以内
2 cm-1
以内
0.02 eV (0.4 kcal/mol) 以内
しかし、これは小さな分子に限られる
理論計算の精度
SCF: 99% of the total energy, Error 1eV/electron pair (N2-3)
MP2: 80% of the last 1%, Pair correlations (Doubles) (N4-5)
CISD: 90% of the last 1%, Singles and Doubles (N6)
CCSD: 95% of the last 1%, Coupling of pair correlations (N6)
CCSD(T): >99% of the last 1% (N7)
MRMP: 80% of the last 0.5%, Dynamical & Nondynamical (N6)
Nano-bio simulation
With the emergence of peta-scale computing platforms
we are entering a new period of modeling. The computer
simulations can be carried out for larger, more complex,
and more realistic systems than ever before.
Linear-scaling methods make this a possibility!
理論計算の問題点
系の大きさ(N)依存性が極めて大きい
200
O(N3)
計算時間はHF, DFT法で
N2.5-3.0で増加
Ｔｉｍｅ
高精度計算ではN4-7で増加
150
100
50
O(N)
0
0
1
2
3
4
5
6
ＳｙｓｔｅｍＳｉｚｅ (N )
単純にいえば計算時間がかかりすぎる
7
8
9
10
理論計算の問題点
系の大きさ(N)依存性が極めて大きい
200
O(N3)
計算時間はHF, DFT法で
N2.5-3.0で増加
Ｔｉｍｅ
高精度計算ではN4-7で増加
150
100
50
O(N)
N依存性を下げる理論
0
0
1
2
3
4
5
6
ＳｙｓｔｅｍＳｉｚｅ (N )
(Linear Scaling)が求められている！
7
8
9
10
ネックはCoulomb積分の計算
Coulomb積分とは(1/r12)を基底関数ではさんだ積分

1
 p q r  s     p (r1 )q (r1 ) r (r2 )s (r2 )dr1dr2
r12
基底関数の数をNとすると2電子積分の数はN4のオーダーで
ある。Nが大きくなると、あっという間にCoulomb積分の数
は莫大になってしまう。分子計算の悩みである
例）原子数1,000とすると基底関数の数はその約10倍
基底関数の数=10,000, 積分の数=1016 (10Peta)
Linear-Scaling Gaussian and Finite-Element
Coulomb (GFC) Method
Y.Kurashige, T.Nakajima, and K.Hirao,
J.Chem.Phys., 126, 144106 (2007)
M. Watson, Y.Kurashige, T.Nakajima, and K.Hirao,
J.Chem.Phys., 128, 054105 (2008)
DFT may be the only tool that enables
us to carry out accurate simulations for
larger systems with reasonable
computational cost. Coulomb integral
evaluation part is very often the most
time consuming step in particular with
GGA functionals.
Coulombポテンシャル
分子系の静電ポテンシャルは
• 原子核上では鋭くスパイク状に
立っている
• Coulomb力は長距離力であり、
原子核から遠く離れたところでも減
衰せずに残っている
性格の異なる2つの領域を
どうやってうまく表現するか？
ClF3
The Gaussian and Finite Element
Coulomb (GFC) method
GFC法ではCoulombポテンシャルはGauss型と有限要素の混
合基底を用いて展開され
v(r)  c  (r)  c
FE FE
i
i
i

Gauss Gauss
i
i
(r)
i
Coulombポテンシャルはポアソン方程式を解いて得られる
GFC法はCoulomb積分をLinear Scalingで計算できる
Solving the Poisson equation for V(r)
Solve the discretized Poisson algebraically
Galerkin Method
A c  b
where
Aij   i (r ) j (r )dr
bi  4  i (r )  (r )dr
Construction of the Coulomb Matrix
The Coulomb matrix
J pq    p (r )  q (r )v(r )dr
  ciGauss   p (r )  q (r ) fi Gauss(r )dr 
i
FE
FE
c

(
r
)

(
r
)
f
 i  p q i (r)dr
i
• Only overlap and kinetic-like integrals are required so far
• But we must also consider the boundary condition…
GFC - 1D Alanine -Helix Chain
CPU time
hour
O(N2.9)
75
50
25
O(N2.4)
O(N1.1)
# of basis functions/SVP
O(N0.9)
GFC – 3D Diamond
CPU time
hour
O(N2.9)
125
O(N2.5)
100
75
962 atoms (C650H312)
10,660 basis functions
Total cpu 198 hours
on IBM Power 5
50
O(N2.6)
25
O(N1.2)
# of basis functions/SVP
Human Insulin Hexamer
C257H383O77N65S6
Numbers of residues and atoms are 306 and 4,728, respectively.
ca.27,000 basis functions (SV)
1 day / 32 CPU (24 iterations)
GFC Energy Gradients
 The Coulomb force can be evaluated efficiently by using the GFC
method. The first derivative of the Coulomb energy with respect
to a nuclear coordinate X can be written as

 (r1 ) 
EJ
 
1



v(r1 )
  dr1  dr2 (r1 ) (r2 )   2 dr1 

X X 
r
X



12


 In the GFC, the Coulomb potential v(r) is expanded in the
auxiliary functions. Thus, the GFC energy gradient may also be
represented by auxiliary functions and can be evaluated
approximately by
EJ
  (r )  G
  (r )  FE
FE
 2 ciG  dr 
f
(
r
)

2
c
dr
f i (r )

i
i 



X
 X 
 X 
i
i
  
 2 Dpqci  dr  p q fi (r )
pqi
 X 
 Computational cost for gradient evaluation is equivalent to that of
one SCF iteration and scales with almost linear scaling.
GFC Energy Gradients
Efficiency of GFC Energy Gradient
48
CPU time/ hour
åvéZéûä‘ÅihourÅj
Analytical
GFC
O(N3.2)
36
24
 3D diamonds/ SVP
12
O(N1.3)
0
4000 # of basis
3000
2000
1000
0
 The GFC is much faster than analytical
gradient evaluation.
äÓíÍä÷êîÇÃêî
 The GFC scales with O(N1.3), while the analytical method scales as O(N3.2).
The GFC achieves linear scaling for gradient calculation
分子シミュレーションの2つのアプローチ
 波動関数法（ab initio MO法）
体系的理論、近似を上げれば厳密解に収束
小さな系では極めて精度の高い結果を与える
計算時間の系の大きさ(N)依存性が大きい
（大きな系では計算時間がかかりすぎる）
課題－効率的なアルゴリズムの開発
 電子密度法（密度汎関数法）
簡潔で概念的、大きな系に適用可能
半定量的理論、精度は汎関数に依存
2次の分子物性、van der Waals力や励起状態の記述が悪い
課題－汎関数の開発
Hybrid GGA functional based on the
long-range correction (LC)
J.Chem.Phys., 128, in press.
J.Chem.Phys., 127, 154109 (2007)
J.Chem.Phys., 126, 154105 (2007)
J.Chem.Phys., 120, 8425 (2004)
J.Chem.Phys., 115, 3540 (2001)
Conventional GGA has problems
1. Barrier heights in chemical reactions underestimated.
2. Van der Waals interactions repulsive
3. Excitations using time-dependent DFT for Rydberg
and CT states underestimated
4. Band gaps of insulators too small
5. Optical response function too large
6. And,…
The failure arises from the wrong long range behavior due
to the local character of the approximate xc functionals.
DFTの高精度化
LC-DFT: Long-range corrected DFT
長距離補正を施した新しい汎関数（LC-DFT) を開発し、DFTの課題を解決
電子間反発 1/r12 をEwald分割法を利用して短距離、長距離部分に分け、短
距離部はDFT汎関数、長距離部はHartree-Fock交換で表現
1

1
 L(r12)    L(r12) 
r12
 r12


erf r12  erf r12  2 (1/ a)  2r122 
 L(r12) 
,

e
,
r12
r12



6
1 1  erf r12  erf r12 


r12
r12
r12
5
4
1  erf (r12 )
r12
1
r12
3
erf (r12 )
r12
2
GGA 交換
HF 交換
 は分割を制御するパラメータ
1
0.5
1
1.5
LC-DFTで2次の分子物性の高精度な記述に成功
2
LCgau交換汎関数
近距離での柔軟性を確保するために長距離項を次のように
表現する
22
e
r
f(

r
) 2


(
1
/)
a


k e r
r

1
r12 は次のように表わされる
1 erf(r12)
2 (1/ a) 2r2 1  erf(r12)
2 (1/ a) 2r2

k
e

k
e
r12
r
r


もっとも優れた交換汎関数の１つ
2次の分子物性の改善
励起エネルギー計算の平均誤差(eV)
Hyperpolarizabilityの改善
27 valence excitations (N2, CO, H2CO, C2H4, C6H6)
B3LYP
BLYP
B3LYP
LC-BLYP
SAC-CI
41 Rydberg excitations (N2, CO, H2CO, C2H4, C6H6)
MP2
LC-BOP
HF
BLYP
B3LYP
LC-BLYP
SAC-CI
1.54
0.89
0.41
0.19
われわれの理論によってDFTのナノ・バイオ系
への応用が可能になった
2.0
CT (R ) - CT (5.0Å) (eV)
0.36
0.37
0.32
0.37
LC-BOP
BOP
AC-BOP
LC-PBEOP
PBEOP
LC-BLYP
BLYP
B3LYP
SVWN
SAC-CI
1.5
1.0
CT励起状態の記述
C2F4
R
C2H4
ΔECT(eV)
0.5
0.0
5
6
7
8
Intramolecule distance R (Å)
9
10
BLYP
B3LYP
LC-BLYP
Exptl.
5.40
7.49
12.49
12.5
LC-GGAによる化学反応の障壁 (kcal/mol)
BOP B3LYP LC-BOP Exp.
C2H6+H→C2H5+H2
NH3+H→NH2+H2
H2O+H→OH+H2
H2O+OH→OH+H2O
CH3F+H→CH2F+H2
CH3F+H→CH3+HF
CH3Br+H→CH3+HBr
CH3Cl+Cl-→Cl-+CH3Cl
CH3Br+Br-→Br-+CH3Br
1,2,3,4-C2N4H2→N2+2HCN
CH3Cl+H→CH3+HCl
….
誤差 (100以上の反応)
2.8
4.7
10.9
-3.8
3.4
16.8
0.5
-3.8
-5.5
27.3
3.4
4.9
7.5
12.8
3.8
5.9
21.0
1.9
-2.1
-3.9
39.7
5.2
6.8
12.3
18.0
5.9
9.5
24.9
3.8
3.8
2.6
50.5
8.0
8.7
5.6
2.6
7.3
11.4
18.5
8.6
9.0
30.1
5.8
2.9
1.7
51.8
10.4
希ガス2量体のポテンシャル・エネルギー
Basis set: aug-cc-pVTZ, BSSE-corrected
0 .3
0 .3
H e N e (L C - B O P + A L L )
H e - H e (L C - B O P + A L L )
B o n d en erg y (kca l/m o l)
N e - N e (L C - B O P + A L L )
N e - N e (E x p t.)
0 .1
A r- A r (L C - B O P + A L L )
A r- A r (E x p t.)
0
-0 .1
-0 .2
B o n d en erg y (kca l/m o l)
H e - H e (E x p t.)
0 .2
H e N e (E x p t.)
0 .2
H e A r (L C - B O P + A L L )
H e A r (E x p t.)
N e A r (L C - B O P + A L L )
0 .1
N e A r (E x p t.)
0
-0 .1
-0 .2
-0 .3
2 .4
2 .8
3 .2
3 .6
4 .0
4 .4
4 .8
B o n d d ista n ce (Å)
Homo-nuclear
5 .2
-0 .3
2 .4
2 .8
3 .2
3 .6
4 .0
4 .4
B o n d d ista n ce (Å)
Hetero-nuclear
4 .8
Van der Waals 相互作用
ClF…He
Naphthalene Dimer
Our DFT with aug-cc-pVQZ
Collinear minimum : -76.43
(Re = 3.41Å）
Our DFT with 6-31+G**
R
8
P
C
T
PD
P CCSD(T)
C CCSD(T)
T CCSD(T)
PD CCSD(T)
6
4
E (kcal/mol)
He
cm-1
2
Ө
T-shaped minimum:
-46.53 cm-1
(Re=3.19Å,θe=70.0°)
Antilinear minimum:
-30.04 cm-1
(Re=3.91Å）
Cl F
CCSD(T) with aug-cc-pVQZ+BF(3s3p2d2f1g)
0
Collinear minimum : -63.53 cm-1
(Re = 3.54Å）
-2
T-shaped minimum:
-41.09 cm-1
(Re=3.23Å,θe=70.1°)
-4
Antilinear minimum:
-33.80 cm-1
(Re=3.95Å)
-6
-8
2.8
3.4
4
4.6
5.2
R (Angstrom)
5.8
6.4
Our DFT
6.16 kcal/mol
CBS CCSD(T) 5.73 kcal/mol
7
in cm-1
50.0
Exfoliation energy (meV/atom)
45.0
Exfoliation energies per C-atom in polycyclic aromatic
hydrocarbon trimers (LC-BOP/6-31++G**)
Experimental data 35, 43, 52 meV
Previous calculation 8 ~170 meV
40.0
34.0
35.0
36.1
35.1
34.1
29.5
30.0
37.3
32.3 (69.0)
25.0
28.5
(49.1)
17.9
20.0
Present (MP2)
dimer
trimer
17.4 (34.4)
15.0
10.0
E(n)  E(n 1)  6n  6(n 1) 
2
 ( ) 
2
( )
120 cartesiancoordinates
Atom
X
,
72cartesian coordinates
Atom
E(n)  E(n 1)
(n ).
6n2  6(n 1)2
(C6H6)3 (C24H12)3
(C54H18)3
C(1)
C(2)
C(3)
C(4)
C(5)
C(6)
C(7)
C(8)
C(9)
C(10)
C(11)
C(12)
C(13)
C(14)
C(15)
C(16)
C(17)
C(18)
C(19)
C(20)
C(21)
C(22)
C(23)
C(24)
C(25)
C(26)
C(27)
C(28)
C(29)
C(30)
C(31)
C(32)
C(33)
C(34)
C(35)
X
Y
Z
1.230
1.229
-0.001
-1.230
-1.229
0.001
2.458
2.457
1.227
-0.002
-1.233
-2.462
-2.461
-3.690
-3.689
-2.458
-2.457
-1.227
0.002
1.233
2.462
2.461
3.690
3.689
3.686
3.685
2.454
1.225
-0.005
-1.234
-2.465
-3.694
-3.692
-4.921
-4.920
-0.708
0.710
1.418
0.708
-0.710
-1.418
1.420
2.838
3.546
2.836
3.544
2.834
1.416
0.706
-0.712
-1.420
-2.838
-3.546
-2.836
-3.544
-2.834
-1.416
-0.706
0.712
3.548
4.966
5.674
4.964
5.672
4.962
5.670
4.960
3.542
2.832
1.414
-0.034
0.041
0.076
0.034
-0.041
-0.076
0.083
0.159
0.193
0.152
0.186
0.145
0.069
0.027
-0.048
-0.083
-0.159
-0.193
-0.152
-0.186
-0.145
-0.069
-0.027
0.048
0.200
0.276
0.311
0.269
0.303
0.262
0.296
0.255
0.179
0.138
0.062
(C96H24)3
C(1)
C(2)
C(3)
C(4)
C(5)
C(6)
C(7)
C(8)
C(9)
C(10)
C(11)
C(12)
C(13)
C(14)
C(15)
C(16)
C(17)
C(18)
C(19)
C(20)
C(21)
C(22)
C(23)
C(24)
C(25)
C(26)
C(27)
C(28)
C(29)
C(30)
C(31)
C(32)
C(33)
C(34)
C(35)
-0.001
-1.229
-1.228
0.001
1.229
1.228
-2.458
-3.686
-3.686
-2.457
-2.456
-1.227
0.001
1.230
2.458
2.458
3.686
3.686
2.457
2.456
1.227
-0.001
-1.230
-2.458
-4.915
-6.144
-6.143
-4.914
-4.914
-3.685
-3.684
-2.455
-1.227
0.002
1.231
(C150H30)3
Y
Z
1.418
0.706
-0.712
-1.418
-0.706
0.712
1.413
0.701
-0.718
-1.424
-2.843
-3.549
-2.837
-3.543
-2.831
-1.413
-0.701
0.718
1.424
2.843
3.549
2.837
3.543
2.831
1.407
0.695
-0.723
-1.430
-2.848
-3.555
-4.973
-5.680
-4.967
-5.674
-4.962
-0.066
-0.084
-0.018
0.066
0.084
0.018
-0.167
-0.185
-0.119
-0.035
0.031
0.114
0.132
0.216
0.233
0.167
0.185
0.119
0.035
-0.031
-0.114
-0.132
-0.216
-0.233
-0.268
-0.286
-0.220
-0.137
-0.071
0.013
0.079
0.163
0.180
0.264
0.281
Kohn-Sham Method without SCF Procedure
Nakajima and Hirao, J.Chem.Phys., 124, 184108 (2006)
The approach is based on the low sensitivity of the
density to the choice of the functional and the basis set.
The total electron density in the ground state can be
well represented in terms of the density evaluated using
the low-quality basis set and the low-cost xc functional.
The error is remedied by the second-order perturbation
theory.
Dual-Level DFT
Solve KS equation
with low-quality basis set & low-level functional
and obtain a total density
Use a frozen density approximation
and evaluate the total energy
with high-quality basis set & high-level functional
Total Electronic Energy
The reference energy of the KS total electronic energy is given by
( 0)
L
 pq PQ tEX  DpqDPQL  pQ Pq t XCEXC
EKS
 2 Dpqhpq  2 DpqDPQ
pq
pq PQ
pq PQ
Since no rotations between occupied and virtual orbitals are allowed,
Brillouin theorem is not satisfied. The correction to the KS energy
is evaluated perturbatively by
2
Fia
(1)
EKS  2
i
a i   a
occ vir
The KS total electronic energy is given by
(0)
(1)
EKS  EKS
 EKS
DL-DFT Energy Gradients
 Like the conventional DFT gradients, the DL-DFT energy
gradients can be written as
x
x
EDL
x
  pq hpq
  JpqPQ  pq PQ    KpqPQ  pQ Pq 
x
pq
pq PQ
pq PQ
x
x
 Wpq S pq
 tXC EXC
 X pq f pqx
pq
 The CP-KS (Z-vector) equation has to be solved in the evaluation
of the DL-DFT gradients, since the DL-DFT is not variational.
 However, the procedure is efficient because the CP equation is
solved only in the low-level set.
Low
x
x
A
U

B
 iajb jb ia
jb
Accuracy
Malachite green(C23H25N2)
Low-level：LDA/6-31G
High-level：B3LYP/6-31G**

High-level
Low-level
Dual-level
B3LYP/6-31G** LDA/6-31G

Total energies (au)
-1000.819997 -991.491120 -1000.815467
Errors in the total energy (au)
9.328877
0.004530
MAE
Bond length (A)
0.0095
0.0002
Bond angle
0.38
0.01
Dihedral angle
1.10
0.05

The large molecular calculation can be performed with the DL-DFT without a
loss of accuracy.
Diels-Alder Reaction
→
＋
40
30
Low-level
n-cyclohexylmaleimide
34.55
Dual-level
-1
エネルギー差/kcal･mol
triphenylene
20
11.88
10
4.82
0
High-level
B3LYP/6311++G**
0
-10
Low-level
-20.87
LDA/6-311G
-20
-30
反応物
Reactants
遷移状態
Transition
state
生成物
Product
Valinomycin (C54H90N6O18)
Valinomycin is highly selective for potassium ions over sodium
ions within the cell membrane.
Low-level:LDA/6-31G
High-level :B3LYP/6-31G**
with PCM
K+： -24.7 kcal・mol-1
Na+： -13.3 kcal・mol-1
理論化学の現状と将来
宇宙の誕生から生命の神秘まで、化学がカバーする
領域は広い。また自然の懐は深い。多様な自然に比
べるとわれわれのシミュレーションはまだまだかも
知れない。
しかし現状でも結構面白いことができる。はっきり
定義された問題を解く腕っぷしにかけては実験に決
して劣るものではない。シミュレーションは従来の
実験、理論とは異なる新しい研究手法を実現し、科
学にブレークスルーをもたらすものとして期待され
ている。謙虚さを失わずに自然に向かえば、理論化
学はますます魅力ある存在になるであろう。
理論の果たすべき役割
シミュレーションは単に経験に追随するのではなく、実験
に先行しなければならない。
今、理論化学、シミュレーションの力量が試されており、
物質科学への影響力、貢献を増大させねばならない。
理論的に難しい問題は残されているものの、シミュレー
ションを面白い現実の問題に応用した方がいいのでは
ないか。新しい方法論、コンピュータの能力がそれを可
能にしつつある。
アボガドロ定数
質量数12の炭素12Cの12gの中に存在する
原子の数と等しい数の粒子を含む系の物
質量 (単位: モル、mol)
1
N A  6.022136710 mol
23
アボガドロ定数はミクロとマクロの架け橋
陽子の質量×NA～
 水素原子の質量
mp  N A  1.672621024 g  6.02213671023 mol1
 1.0073g / mol
素電荷 (e)×NA=ファラデー定数 (F)
e  N A  1.602181019C  6.02213671023 mol1
 96485C / mol  F
ボルツマン定数 (k)×NA=気体定数 (R)
k  N A  1.38065031023 JK 1  6.02213671023 mol1
 8.31447JK 1 / mol  R
Development of Computers
アボガドロ数
FLOPS
Floating point number
operations per second
Cray-1
0.16 TFLOPS
Earth Simulator 35.61 TFLOPS
Blue Gene/L 280.6 TFLOPS
TFLOPS=1012 FLOPS
Moore's Law
“The transistor density of
integrated circuits, with
respect to minimum
component cost, doubles
every 24 months”
コンピュータの発展
Petascale computing hardware is
around the corner.
Petascale computing chemistry
should be science driven.
Petascale resources will enable
chemistry researchers to enter a new
modeling.
Dynamics in chemical
transformation
Very Accurate results
Chemistry at periodic interfaces
Coupling multiple scales
Benchmarking abilities
Heavy element chemistry
コンピュータは数学の脆弱性を補うことができるか？
(GFlops)
1.0E+08
Todai
1st in TOP500
10th in TOP500
1.0E+07
1.0E+06
超並列型スーパコン
Next Generation
Super Computer
Peta
Earth Simulator
1.0E+05
1.0E+04
1.0E+03
Tera
1.0E+02
1.0E+01
10Peta (1016) FLOPS Massively Parallel Machine
1.0E+00
Total budget 1 billion US$ (2006～2012)
Giga
20
13
20
11
20
09
20
07
20
05
20
03
20
01
19
99
19
97
19
95
19
93
19
91
19
89
19
87
19
85
19
83
1.0E-01
(Year)
Core excitation energies in eV (singlet)
singlet
Molecule
Transition
C2H2
C1s→π*
C2H4
CVCVR-B3LYP
B3LYP
LC-BOP
CAM-BLYP
LCgau-BOP
LCgau-BOP
μ:0.47
α:0.2; μ:0.4
270.4
283.1
μ:0.42;a:0.01 μ:0.42;a:0.033
1;k:-18.0
5;k:-5.9
286.1
286.2
BLYP
B3LYP
269.7
275.2
286.0
286.1
( -16.1)
268.9
( -10.6)
274.3
(+0.2)
285.1
(+0.3)
285.1
( -15.4)
269.4
( -2.7)
282.0
(+0.4)
C1s→π*
C1s→π*
( -15.8)
270.2
( -10.4)
275.2
(+0.4)
286.0
(+0.4)
286.0
( -15.3)
270.1
( -2.7)
282.6
(+0.5)
CH2O
C1s→π*
( -15.8)
271.3
( -10.8)
276.2
( 0.0)
286.9
( 0.0)
286.9
( -15.9)
271.1
( -3.4)
283.5
( -0.1)
CO
N1s→π*
( -16.1)
382.5
( -11.2)
388.5
( -0.5)
401.3
( -0.5)
401.3
( -16.3)
382.5
( -3.9)
397.3
( -0.4)
N2
O1s→π*
( -18.5)
509.4
( -12.5)
516.7
(+0.3)
531.4
(+0.3)
531.4
( -18.5)
509.6
( -3.7)
526.7
(+5.6)
CH2O
O1s→π*
( -21.4)
512.3
( -14.1)
519.8
(+0.6)
534.5
(+0.6)
534.5
( -21.2)
512.7
( -4.1)
530.0
(+11.7)
CO
MAE
Core
Excitation
Atomiztion
Energy
( -21.9)
17.9
( -14.4)
12.0
(+0.3)
0.3
(+0.3)
0.3
( -21.5)
17.7
( -4.2)
3.5
(+11.6)
4.3
(+1.0)
0.6
9.8
7.5
10.1
7.8
5.3
7.8
4.74
2.53
2.49
2.43
2.50
rms
rms
Barrier Height
285.2
285.9
287.0
406.6
542.5
545.8
Expt.
285.8
(+0.3)
285.0
284.7
(+0.3)
285.6
286.0
( -0.4)
286.5
287.4
( -0.9)
401.6
401.0
(+0.6)
531.8
530.8
(+1.0)
535.2
534.2
Core excitation energies in eV (triplet)
LC-BOP
CAM-BLYP
μ:0.47
α:0.2; μ:0.4
269.9
282.3
C1s→π*
( -15.9) ( -10.6)
268.5 273.6
( -15.4)
268.8
( -3.0)
281.2
(+0.3)
284.5
( -0.0)
284.2
284.3
C1s→π*
( -15.8) ( -10.7)
270.5 275.1
( -15.5)
270.2
( -3.1)
282.3
(+0.2)
285.9
( -0.1)
285.1
285.95
( -15.5) ( -10.9)
15.7
10.7
( -15.8)
15.6
( -3.7)
3.3
( -0.1)
0.2
( -0.8)
0.3
Molecule
Transition
C2H2
C1s→π*
C2H4
CO
MAE
Core
Excitation
BLYP B3LYP
269.4
274.7
LCgau-BOP
LCgau-BOP
μ:0.42;a:0.011;k:- μ:0.42;a:0.033 Expt.
18.0
5;k:-5.9
285.6
285.3
285.3
Radial Distribution of Au atom
4f
DK3
2.0
NR
6s
5s
1.5
Relativistic
5p
Nonrelativistic
1.0
5d
0.5
6s
0.0
0
1
2
3
4
5
6
7
8
The relativistic kinematics of the electrons leads to markedly contracted s
orbitals and to a certain extent contracted p orbitals, while orbitals with higher
angular momentum are expanded compared with their non-relativistic
counterparts.
2.80
2.70
Bond distances of Cu2, Ag2 and Au2 (Å)
Exptl.
Nonrelativistic
Relativistic
2.60
2.50
2.40
2.30
2.20
Cu2
Ag2
Au2
Disparities in properties between Au and its similar group
members, Ag and Cu arise mainly from the relativistic effects.
Changes of spectrum due to relativity
Photoionization Spectrum of OsO4
B
C
D
A
E
19
18
17
16
15
14
Ionization energy (eV)
Nonrelativistic Treatment
Relativistic Treatment
13
12
-303.4 curves due to relativity
Changes of potential energy
-19303.5
Spin-free CASSCF
-19303.6
Spin-orbit CASSCF
Spin
Au 2D
+ Si 1-Sorbit
-303.5
Au 2D + Si 1D
Energy (a.u.)
-19000 a.u.
Doublet
Doublet
Quartet
Quartet
-19303.7
Au 2D + Si 3P
-19303.8
Ω states
-303.6
Au 2S + Si 1S
3P
Au 2DAu
+ Si2D
+ Si 3P
-303.7
Au 2S + Si 1D
Au 2S + Si 3P
-303.8
-19303.9
-303.9
Au 2S + Si 3P
-19304.0
-304.0
-304.1
-19304.1
3.0
4.0
5.0
6.0
7.0
R (a.u.)
8.0
9.0
10.0
3.0
AuSi Molecule
4.0
5.0
6.0
7.0
R (a.u.)
8.0
9.0
10.0
119Sn
NMR Chemical Shifts of SnXn
Spin-orbit effect
Changes of the chemical reaction due to relativity
O
C
Pt+
R
P
Reactant
O
C
P
R
O
C
Pt+
P
P
Transition state
R
Pt+
P
P
Product
CO insertion to cationic Pt complex
Pt
Difference density map
at the transition state
calculated by relativistic
and nonrelativistic
treatments.