4P081
アセチルコリンエステラーゼと阻害剤間の特異的相互作用：
フラグメント分子軌道計算と生化学実験○村川孟 1、岡田翼 1、M. T. H. Khan2、栗田典之1(1 豊橋技術科学大学大学院、2 University of Tromsø)Specific interactions between acetylcholineesterase and its potent inhibitors:
ab initio fragment molecular orbital calculations and biochemical experiment○T. Murakawa1, T. Okada1, M. T. H. Khan2, N. Kurita1(1 Toyohashi University of Technology、2 University of Tromsø )【はじめに】	 	
  神経細胞は、神経信号の受信とともに対応する動作を決定し、隣接細胞へ新たな信号を発信する[1]。その情報処理において、アセチルコリン (ACh) は、主に運動神経細胞から筋肉細胞へ神経信号を伝える働きをしている[2]。神経細胞は、特有の信号を神経系から受け取り、シナプスから筋肉細胞へ向けて ACh を放出し、それによって筋肉細胞の受容体が開き、筋肉は収縮する。
信号伝達に使用された ACh は、筋肉細胞と神経細胞の間にあるシナプスに存在するアセチルコ
リンエステラーゼ(AChE)によって、コリンと酢酸に分解される[3]。この分解反応により、神経信号は効率的に遮断され、ACh の断片は次の神経信号に用いる新たな神経伝達物質の再構築のために利用される。AChE による ACh の分解反応を部分的に阻害する薬は、神経伝達物質の濃度上
昇、及び筋肉細胞への神経信号の強化を促進する可能性がある。
	
  本研究では、フラグメント分子軌道 (FMO)[4,5] 計算を用い、AChE と阻害剤間の特異的相互作
用を電子レベルで解析し、既存の阻害剤の中で AChE に最も強く結合する阻害剤を明らかにし、
生化学実験と比較できる結果を得た。その結果を基に、コンピュータ内で複数の新規リガンドを作成し、AChE との結合特性を解析し、ACh の分解反応を阻害する新規リガンドを提案した。
	
 
【計算手順】	 	  タンパク質立体構造シミュレーションプログラム ICM (Internal Coordinate Mechanics) を用い、AChE と阻害剤の複合体構造を作成した。Figure 1 に示す 7 種類の阻害剤 (Bx) に着目し、置換基の位置や種類が AChE と Bx 間の特異的相互作用に与える影響を解明した。まず、AChE+Bx 複合体構造を古典分子力学（MM）計算プログラム AMBER12 を用い、水中で最適化した。その際、複合体表面から 8 Å 以内の領域に水和水を付加し、水和した複合体の構造を最適化した。さらに、
AChE と Bx 間の結合特性を明らかにするため、最適化構造の電子状態を FMO 計算により解析し、
AChE と Bx 間の結合エネルギー、及び特異的相互作用を明らかにした。その結果を基に、AChE
により強く結合可能な新規リガンドを提案し、AChE との結合特性を解析した。
	
 
【計算結果と考察】	
 Figure 2 に示すように、FMO 計算により求めた AChE と Bx 間の結合エネルギー (B.E.) と生化学実験で求めた Bx の解離定数 Kd の間には、R2 = 0.85 の強い相関がある。従って、B.E.の大きさ
により、Bx が AChE の機能を阻害する効果の強さを定性的に評価できる。また、AChE のどのア
ミノ酸残基が Bx の結合に重要であるかを明らかにするため、AChE の各アミノ酸残基と Bx 間の相互作用エネルギー (I.E.) を解析した。Bx との引力相互作用が最大の 3 個のアミノ酸を Table 1に示す。多くの Bx が Asp74、Tyr124、Tyr341 と強く引力相互作用をしている。従って、荷電ア
ミノ酸である Asp や芳香族アミノ酸である Tyr が、AChE と Bx 間の結合に重要であると言える。
	
  7 種類の阻害剤 Bx の中で、結合エネルギーが最大である B7 は、Asp74、Arg296、Tyr341 と強く引力相互作用するが、Tyr124 との引力相互作用は弱い。このような特異的相互作用の原因を明らかにするため、B7 周辺の結合構造と電荷分布を調べた。Figure 3 に示すように、B7 の末端部の NO2 基及びベンゼン環が、AChE のアミノ酸残基と弱い水素結合、あるいは Stacking 相互作用

をしている。これは B7 だけの特徴ではなく、全ての Bx に対して共通の特徴である。B7 の H 原子の電荷は 0.3、Asp74 の O 原子の電荷は−0.8 であり、これらの原子間の距離は 2.4 Å であるため、B7 と Asp74 間の相互作用エネルギーは、−28.1 kcal/mol と大きくなる。また、B7 と Tyr341のベンゼン環同士が 3.6 Å の距離でスタッキング構造を形成しており、B7 と Tyr341 の間には、−11.9 kcal/mol の強い引力相互作用が生じる。さらに、B7 のみが有する NO2 基が Arg296 と弱い水素結合を形成し、その結果、B7 と Arg296 間に−11.8 kcal/mol の引力相互作用が生じる。この相互作用が、AChE と B7 間の結合エネルギーが既存リガンドの中で最大になった主な原因であると考えられる。そこで、B7 をベースにして、様々なリガンドを作成し、AChE のリガンド結合部位に結合させ、AChE と新規のリガンド間の結合特性を FMO 計算で解析し、AChE の阻害剤として機能する可能性のある新規リガンドを提案した。詳細は、当日のポスターにて発表する。

B1

B2

B3

B4

	
 

	
 

B5

B6

B7

Figure 1 Structures of ligands Bx for inhibiting AChE function

Figure 2 Correlation between binding energies
for AChE+Bx complexes and the distribution
coefficient Kd obtained by our experiment

Figure 3 Interacting structure and charge
density between B7 and AChE residues

Table 1 Interaction energies (I.E.) (kcal/mol) between AChE residues and Bx
B1
B2
B3
B4
B5
B6
B7
Tyr124 −11.0 Tyr124 −10.6 Asp74 −15.1 Asp74 −31.6 Asp74 −23.8 Tyr124 −10.8 Asp74 −28.1
Tyr341 −7.6 Tyr341−9.8 Tyr341 −10.9 Tyr341−6.2 Tyr341 −11.1 Glu294−8.6 Tyr341 −11.9
Phe295 −7.9 Phe295−9.4 Tyr124−9.9 Tyr124−6.0 Phe295−9.0 Tyr341−8.0 Arg296 −11.8[1] K. A. Phelan, et al., J. Neurosci. 1990, 10, 2699. [2] S. Camp, et al., J. Neurosci. 2008, 28, 2459.
[3] Y. Dudai, et al., Proc. Natl. Acad. Sci. USA, 1973, 70, 2476.
[4] K. Kitaura, et al., Chem. Phys. Lett.,1999, 312, 319.
[5] Y. Mochizuki, et al., Chem. Phys. Lett., 2005, 410, 247.

4P082
がん細胞受容体に特異的に結合する新規ペプチドの提案とその結合特性：
リガンドドッキング及びフラグメント分子軌道計算による解析○水島達朗 1、川野裕章 1、霞知世 1、小林浩 2、栗田典之 1（1 豊橋技術科学大学、2 奈良県立医科大学）
A proposal of new peptide inhibitor binding to the cancer cell receptor and the analysis of its binding
properties to the receptor by ligand-docking and fragment molecular orbital calculations○T. Mizushima1, H. Kawano1, T. Kasumi1, H. Kobayashi2, N. Kurita1(1 Toyohashi University of Technology, 2 Nara Medical University)【はじめに】
がん細胞は様々なタンパク質分解酵素を産生し、細胞周辺組織を分解・浸潤し、血管に進入して血管内を移動し、標的臓器の血管内皮細胞に接着・増殖し、転移巣を形成する。従って、がん転移を抑制するには、がん浸潤を抑制することが効果的である。近年の生理医学実験[1]により、がん浸潤には urokinase-type plasminogen activator (uPA)という酵素の作用が重要であることが明らかになった。がん細胞の表面には、がん細胞受容体である uPAR (uPA receptor)が存在し、この uPARに uPA が結合することにより、がん浸潤が引き起こされる。従って、uPAR と uPA の結合を阻害できれば、がん浸潤は抑制できると考えられる。実験[2]により、uPA の ATF (amino-terminal
fragment)部位が uPAR と結合し、特に、ATF 中の 19 から 31 番目のアミノ酸が、uPAR と uPA の結合に重要であることが明らかになっている。
我々は、分子シミュレーションを用いて、uPAR と uPA 間の特異的相互作用を原子・電子レベルで解析し、上記のアミノ酸以外で uPA の正荷電のアミノ酸 Lys46 と Lys61 が結合に重要である
ことを明らかにした[3]。それらは uPAR のリガンド結合ポケットの上部に存在する負荷電の Glu33、Glu34、Glu36 と強く相互作用する。Figure 1 に示すように、uPAR のポケットの周囲には、負電荷のアミノ酸である Glu が多く存在し、Glu33、Glu34、Glu36 の反対側に、Glu132、Glu134、Glu135
が存在している。これらのポケットの両側に存在する荷電アミノ酸の Glu に強く結合するペプチドを提案できれば、uPAR のポケットを塞ぎ、uPAR と uPA の結合を阻害できると考えられる。
我々はこれまでに、分子シミュレーションを用いて、uPAR のポケットを塞ぐことが可能なペプチドを提案した[4]。本研究では、これらのペプチドを基に、アミノ酸置換を行ない、錠剤として生成可能な低分子ペプチドを新たに提案し、それらと uPAR 間の結合特性を、高精度に解析した。
【計算手順】
本研究では、uPAR の初期構造として欠損部位がGlu132Gly133
Glu134
Glu135

Glu33
Glu34
Gly35
Glu36

Ligand binding pocket

最小である X 線構造(PDB ID: 3BT1)を採用し、ホモ
ロジーモデリングを用いて uPAR の欠損部位を補い、uPAR の初期構造を作成し、古典分子力学(MM)計算

Fig.1 The ligand binding pocket of uPAR

により、uPAR の水中での安定構造を決定した。次に、Lys を 2 個含む[Gly-Lys-(Gly)n-Lys-Gly](n=1～6)のアミノ酸配列を持つペプチド(KGn)を提案

し、分子軌道計算プログラム Gaussian09 の MP2/6-31G(d,p)計算により安定構造を求め、電荷パラメータを決定した。そして、リガンドドッキングプログラム AutoDock4.2 を用いて、uPAR にペ
プチドをドッキングし、uPAR とペプチドの複合体の初期構造を作成し、MM 計算により水中での安定構造を決定した。さらに、第一原理フラグメント分子軌道(FMO)計算により、uPAR とペプチド間の結合特性を電子レベルで解析し、がん転移抑制剤としての可能性を検討した。
【計算結果と考察】
ペプチドが uPAR のリガンドポケットを塞ぐためには、25.0～37.0 Å の長さが必要である。ペ
プチド KGn の KG1～KG6 の長さを調べ、n が 3 以上であれば 25.0 Å 以上の長さが確保できることが明らかになった。そこで、KG3～KG6 のペプチ
ドを uPAR にドッキングし、FMO 法の MP2/6-31G
計算を用いて、uPAR とペプチド間の結合エネルギ

Table 1 Binding energies (B.E.)(kcal/mol) between uPAR and peptide

ーを求めた。その結果、Table 1 に示すように、KG5

Structure

B.E.

が uPAR と最も強く結合することが明らかになった。

uPAR + KG3

257.7

さらに、uPAR と KG5 の各アミノ酸間の相互作用エuPAR + KG4

294.2

ネルギーを解析した結果、KG5 は uPAR のリガンド

uPAR + KG5

301.9

ポケットを挟んで対称の位置に存在する Glu36、uPAR + KG6

291.2

Glu134、Glu135 と強く引力相互作用することがわか

uPAR + KG5-Glu

300.1

った。従って、KG5 は uPAR のこれら 3 つのアミノ

uPAR + KG5-Asp242.7酸と結合することで、uPAR のリガンドポケットを塞ぎ、uPAR と uPA の結合を阻害する可能性があると考えられる。
しかし、KG5 は両端の Lys しか uPAR と結合しておらず、中央部の 5 つの Glu が uPAR のどのアミノ酸とも結合していない。また、KG5 は全体で+2 の電荷を有しており、KG5 が人体に入った際にカウンターイオンが付加し、uPAR に特異的に結合できない可能性がある。また、人体に電荷を有するものを与えることは好ましくないと考えられる。そこで、ペプチドの中央部の結合の問題を改善するため、KG5 の中央の Gly を負電荷アミノ酸 Glu、及び Asp に置換したペプチド(KG5-Glu、KG5-Asp)を作成し、KG5 と同様の計算を行い、uPAR との結合エネルギーを求めた。
Table 1 に示すように、uPAR と KG5-Glu 間の結合エネルギーは、uPAR と KG5 間の結合エネルギ
ーとほぼ同等になった。また、置換した Glu 及び Asp は、uPAR の Lys139 と強く引力相互作用した。従って、KG5 の中央部の結合の問題点は、中央の Gly を負電荷アミノ酸に置換することで改善できると考えられる。さらに、KG5 の+2 電荷の問題を改善するため、KG5 の中央部の 5 つのGly 中の 2 個を負電荷アミノ酸に置換し、ペプチド全体の電荷を 0 にしたペプチドを新たに作成
し、uPAR とペプチド間の結合エネルギーを求めた。こちらの計算結果、及び uPAR とペプチド間の結合特性に関しては、当日のポスターで発表する。
【参考文献】[1] A.E. Baker, et al., J. Cancer, 39, 981(2003).
[2] R. D. GmbH, et al., WO98/46731 (1998).
[3] S. Tsuji, et al., J. Mol. Graph. Model., 29, 975(2011); T. Kasumi, et al., Mol. Simu. (2013) in press.
[4] T. Kasumi, et al., J. Phys., 433, 012034 (2013).

4P083

ラクトースリプレッサー二量体と DNA 間の特異的相互作用：
分子動力学計算及び第一原理フラグメント分子軌道計算を用いた解析（豊橋技術科学大学）○松下祐貴、大石叡人、大山達也、栗田典之

Specific interactions between the dimer of lactose repressor protein and DNA :molecular dynamics simulations and ab initio FMO calculations(Toyohashi University of Technology)○Yuki Matsushita, Masato Oishi, Tatsuya Ohyama, Noriyuki Kurita【はじめに】
生体内では、常に DNA 情報が mRNA に転写され、それを基にタンパク質が生成されている。転写制御タンパ
ク質ラクトースリプレッサー (LacR) は、通常は DNA に結合し、mRNA への転写を抑制している。しかし、様々なリガンドが LacR に結合することにより、LacR と DNA間の結合特性が変化し、転写機構が調節されることが実験で明らかになっている[1]。インデューサが LacR に結合
すると、LacR と DNA 間の結合を弱め、LacR が DNA から分離し、転写が活性化する。一方、アンチインデュー
サが LacR に結合すると、LacR と DNA 間の結合が強まり、転写を抑制する。我々[2]は、これまでに、第一原理フラグメント分子軌道 (FMO) 計算により、LacR 単量体とリガンド間、及び LacR 単量体と DNA 間の特異的相互作用を解析し、これらの相互作用に重要な LacR のアミノ酸残基を明らかにした。しかし、結合するリガンドの種類にFig. 1 Structure of LacR-dimer

よる LacR と DNA 間の結合特性の変化は解明できなかっ+ DNA + IPTG complex

た。
本研究では、Fig. 1 に示すように、LacR 二量体と DNA から構成される複合体を考え、そこにリガンドが結合した際の構造変化を、古典分子動力学 (MD) 計算により解析した。更に、MD 計算で求めた幾つかの構造に対し FMO 計算を実行し、結合するリガンドの種類により、LacR と DNA間の特異的結合状態がどのように変化するかを電子レベルで解析した。
【計算手順】
1.

MD 計算による LacR 二量体+DNA+リガンド複合体の水和構造の変化の解析
結晶水を含む LacR 二量体、DNA 及びアンチインデューサ ONPF の複合体の結晶構造 (PDB ID:

1EFA) を、LacR+DNA+ONPF の初期構造として採用した。この構造の DNA 骨格のリン酸基にカウンターイオンとして Na+を付加した。また、LacR にインデューサ IPTG が結合し構造 (PDB ID:
2P9H) を参考に、LacR と IPTG 間の相対位置を求め、その情報を基に、LacR+DNA+ONPF 中の
ONPF を IPTG に置換し、LacR+DNA+IPTG の初期構造を作成した。さらに、LacR+DNA+ONPFの ONPF を削除した構造を、LacR+DNA の初期構造とした。これら 3 つの構造の周囲に水分子を置き、古典 MD 計算プログラム GROMACS により、300 K で 100 ns の MD 計算を実行し、複合体の水中での構造変化を解析した。

2.

第一原理 FMO 計算による LacR と DNA 間の特異的相互作用の解析
リガンド結合による LacR 二量体と DNA 間の特異的相互作用の変化を明らかにするため、MD計算で得た複数の構造の電子状態を FMO 計算プログラム ABINIT-MP ver. 6.0 の MP2/6-31G 法を用いて解析した。その際、LacR と DNA 間、及びカウンターイオン周囲の水分子のみを考慮した。
【結果と考察】リガンド結合が LacR に与える影響を解析するため、MD 計算開始時の LacR 構造と各時間における LacR 構造のアミノ酸の Cα 原子の RMSD を解析した。Fig. 2 に示すように、リガンドが結合した LacR+DNA 複合体の構造に対する RMSD が大きくなる。これは、リガンド結合の影響がリ
ガンド周囲の LacR アミノ酸残基を中心に LacR 全体に及ぶことを示している。
IPTG が結合した LacR 二量体は、Fig. 3 に示すように、0 ~ 30 ns 時に LacR が DNA に沿った方
向に大きく揺らいでおり、その後も同様の揺らぎが観測された。LacR 二量体は、Fig. 4 に示すよ
うに、単量体それぞれの DNA binding domain (DBD) で DNA を挟む形で、DNA との複合体構造を形成している。
IPTG は、疎水基のイソプロピル基を持ち、これが LacR に結合すると LacR の Asn125
や Asp149 と反発する。この反発により、α-helix が DBD の方へと押し出され、それぞれが DBDと水素結合を形成する。これにより LacR 全体が傾いたと考えられる。
一方、ONPF は親水性のフコースを持つため、周囲のアミノ酸と反発せず、上記のような傾きは見られない。現在、MD 計算によって得られた複数の構造に対し、FMO 計算により、LacR とDNA 間の特異的相互作用を解析中であり、その結果は当日のポスターで発表する。

Arg197

RMSD (Å)

5
4
3

Asp1492

Asn1251a-helix

0
0

20

40

60

80

100

MD simulation time (ns)DBD

Fig. 2 Change in RMSD of amino acid residues
of LacR during MD simulations (black:no-ligand, red: IPTG, blue: ONPF)

Fig. 4 Schematics view of the complex with
LacR + IPTG dimer and DNA

12.8 deg.

1.5 ns

5.0 ns

7.7 ns15.0 ns30.0 nsFig. 3 Change in relative conformation between LacR + IPTG dimer and DNA
in the MD simulation: at 1.5, 5.0, 7.7,15.0 and 30.0 ns【参考文献】[1] R. Daber, et al., J. Mol. Biol. 2007, 340, 609. [2] T. Ohyama, et al., J. Comput. Chem. 2011, 32, 1661.

4P084人工塩基対を含んだ DNA 系に対するオーダーN 法 DFT 計算
1（理研 QBiC , 物材機構 2, ロンドン大学 3）○大塚教雄 1, 有田通朗 2, David Bowler3, 宮崎剛 2

Theoretical study on hydrated DNA systems with unnatural base pairby linear-scaling DFT calculations（RIKEN QBiC1, NIMS2, UCL3）Takao Otsuka1, Michiaki Arita2, David Bowler3, Tsuyoshi Miyazaki2【序】 近年の有機合成化学の技術は目覚ましく、生命システムを構成する要素を人工的に創出し、生命現象の理解とともにその機能を向上させる事が可能となってきている。核酸分野では、人工的に設計した非天然型塩基を DNA に導入し、この非天然型塩基を遺伝子情報の拡張コードとした複製・転写・翻訳シ
ステムの理解や再構成、更に人工アミノ酸を含んだ機能性タンパク質の創製を目指した研究が行われている。理研 CLST の平尾等は、独自の設計指針から人工的に作成した非天然型塩基とその塩基対（Ds-Pa, 7-(2-thienyl)-imidazo[4,5-b]pyridine (Ds)と pyrrole-2-carbaldehyde (Pa)）を開発し、この非天然型塩基対を骨格内に封入した人工 DNA は、複製と転写の機能を有する事に成功している[1]。最近では、
この非天然型塩基を天然の DNA に組み込んだ DNA アプタマーを作成し、従来の DNA アプタマーと比べ、標的タンパク質との結合能力機能を飛躍的に向上させる事を可能とした[2]。今後は、効率よく、より洗練された人工塩基対を創製する分子設計指針が必要とされている。
本研究では、我々が開発してきたオーダーN（order-N または linear scaling）法第一原理 DFT 計算プロ
グラム CONQUEST[3]を用いて、超大規模生体分子系に対する計算シミュレーションから、生体分子の機能理解や関連する分子設計技術の開発を行っている。これまでに我々は、天然型の DNA 系に対するオーダーN 法第一原理 DFT 計算の数値検証から、我々のオーダーN 計算手法が極めて高精度・安定である事を示してきた[4]。また、スーパーコンピュータ「京」上での開発と数値検証から、超高並列計算機下でもオーダーN 計算手法の堅牢さを実証している。
これまで、非天然型塩基対を含んだ DNA 系に対しては、まず非天然型塩基対のみに着目した量子化学計算から塩基対間の相互作用評価を行った[5]。次いで超高並列計算機下での超大規模計算に向けたオーダーN 法第一原理 DFT 計算の数値検証を行った。今回は、経験的な van der Waals 相互作用補正を含めた、更なる詳細な解析を行った結果を報告する。
【理論的背景】 CONQUEST で用いられているオーダーN 計算手法は、密度汎関数法における一体の密度行列を最適化する手法である。密度行列の非対角項が局所的であることから、非対角項に対する
cutoff 半径 RL を導入する事でオーダーN 法を実現している：

ρ (r, r ′) = ∑ φiα (r ) K iα , jβ φ jβ (r ′),iα , jβK = 3LSL − 2 LSLSL,

Liα , jβ = 0 for | R i − R j | > RL

ここで、行列 L は補助密度行列であり、電子数一定の条件下で全エネルギーを最小にする密度行列 Lが求められる。また φiα (r ) は、support function と呼ばれる各原子に局在した関数である。今回も擬原子
軌道（PAO）を使った計算を行っている。系の全エネルギーは、cutoff 半径に対して変分的であるという利点を持つことから、オーダーN 法を導入することによって生じる誤差を評価することが可能である。また、原子座標の自由度に対する最適化（構造最適化計算）も安定に行えるという特徴を持つ。

【計算と結果】 図 1 に示すように計算に用いた DNA モデル系は、5'-GGTAAC-Ds-ATGCG-3' と 5'-CGCAT-Pa-GTTACC-3' と した Ds-Pa 塩基対の 1 対を含んでいる。AMBER によって水分子を加える事によって、全原子数 11,912 原子（DNA: 763 原子、Na: 22
原子、H2O: 3709 分子 = 11,127 原子）を作成した。古典分子動力学計算による平衡状態計算後のスナップショットの構造を用いて、オーダーN 法 DFT 計算を行った。これは以前報告した構造と同じものである。
ここでは、人工塩基対 Ds-Pa の 1 対を含んだ DNA12 塩基対モデル系の Force に関する結果を示す。計算条件として、PBE 汎関
数と DZP 精度を持つ縮約 PAO を用いた。また、DFT 計算の vander Waals 補正として DFT-D2 法を用いた。
図 2 は、CONQUEST の SCF 計算による非天然型塩基対 DNA
モデル系の DNA 部分を構成する各原子に働く Force のノルムを示している。茶色は PBE 汎関数による Force、青色は DFT-D による補正を示している。DFT-D による各原子の Force のノルムの補Figure1. DNA system including
one unnatural base pair, Ds-Pa.

正効果は、最大でも 10-3 Ha/a0 のオーダーであり、かなり小さい事が分かる。図 3 は、DNA を構成する塩基単位で算出した Force の各成分を示している。白抜きの記号は、DFT-D2 による補正結果である。このスナップショット構造からは、塩基対単位での Force の大きな差は見られない。更に他のスナップショット
構造の比較を踏まえ、議論を進めて行く予定である。

x-comp.

y-comp.

z-comp.

Figure2. Norm of atomic force of DNA partPa

Ds

Figure3. Atomic force components of
DNA bases【参考文献】 [1] I. Hirao et al., Nature Methods., 3, 729-735 (2006). [2] M. Kimoto et al., Nature
Biotechnology, 2013 doi:10.1038/nbt.2556 [3] D. R. Bowler, T. Miyazaki, Rep. Prog. Phys. 75, 036503(2012). [4] T. Otsuka, et al., J. Phys. Condens. Matter, 20, 294201 (2008). [5] T. Otsuka, T. Miyazaki,
Int. J. Quantum Chem. 113, 504 (2013).

4P085

DNA 塩基対と塩基対とラジカルとのラジカルとの反応機構との反応機構に反応機構に関する密度汎関数計算する密度汎関数計算(豊橋技術科学大学)○清水栄佑、奥津尚子、徳山雄生、栗田典之

DFT study on the reaction mechanism between DNA base-pair and radicals(Toyohashi University of Technology)◯Eisuke Shimizu, Naoko Okutsu, Yuki Tokuyama, and Noriyuki Kurita【はじめに】
DNAは、生物の遺伝情報を記憶する物質であり、細胞内の環境や紫外線などの様々な要因によって損傷する。DNA損傷の要因の一つに放射線がある。放射線はDNAに直接作用し、DNAを構成する分子を電離し、損傷する。さらに、DNAの周囲に存在する溶媒水に作用し、OHラジカルなどのラジカルを生成する。OHラジカルは反応性が強く、糖質や脂質、タンパク質など細胞中の多くの物質と反応する。このOHラジカルがDNAと反応すると、DNAから水素原子が引き抜かれ、DNAに損傷が生じ、細胞の老化や死、がん細胞の発生、遺伝性疾患などを引き起こす原因となることが、これまでの実験により明らかになっている[1,2]。一方、OHラジカルとDNAの反応機構に関する理論的研究は広く行われている[3]が、その多くが古典分子力学法を用いた計算であり、この機
構を電子レベルで明らかにした研究例は少ない。この反応が起る機構を明らかにできれば、OHラ
ジカルによるDNA損傷を抑制する手段を提案できる可能性がある。
我々のこれまでの研究[4]では、密度汎関数法(DFT)に基づく高精度分子シミュレーション手法
を用い、DNA 塩基対の様々な位置に OH ラジカルを配置した構造を最適化し、G-C 塩基対の G のNH2 基、及び A-T 塩基対の A の NH2 基の水素原子が引き抜かれ、OH ラジカルが H2O になる構造が生成されることを明らかにした。そこで、本研究では、同様の分子シミュレーション手法を用い、DNA 塩基対と OH ラジカルの反応機構の遷移構造を、真空中、及び連続溶媒近似の水中で探索する。その結果を基に、この反応の反応障壁エネルギーを求め、DNA 塩基対に OH ラジカルが反応した時に、どのような反応が起り易いかを、電子レベルで解明する。この結果は、放射線の間接的効果により、DNA 塩基に損傷が起る現象の初期段階の解明に役立つと考える。
【計算手法】
まず、核酸構造作成プログラム NUCGEN を用い、G-C、あるいは A-T の 1 塩基対を含む二重鎖を作成した。その構造から、バックボーンを削除し、塩基対のみを切り出した。塩基対のバックボーンとの結合部位は、CH3 基でキャップした。この塩基対の構造を、第一原理分子軌道計算プ
ログラム Gaussian09 の M06-2x/6-31G(d,p)法を用い、真空中で最適化した。次に、これらの構造のNH2 基の周辺に OH ラジカルを配置し、反応前、反応後、及び反応の中間状態に対する初期構造を作成した。その際、OH ラジカルの酸素原子と DNA 塩基の水素原子間、及び OH ラジカルの水素原子と DNA 塩基の窒素原子間の距離は、アデニン塩基と OH ラジカルの反応機構を計算した
Cheng らの論文[5]を参考に設定した。作成した反応前及び反応後の構造を、真空中、及び連続溶
媒近似の水中で最適化した。その際、Gaussian09 の UM06-2x/6-31++G(d,p)法を用い、全電荷を 0、スピンを Doublet に指定した。また、連続溶媒近似には、Gaussian09 の IEFPCM 法を用いた。
最後に、これらの構造を初期構造として、Gaussian09 の STQN 法を用い、反応の遷移状態を探索した。その際、DFT 計算の条件は最適化の場合と同じ条件を用いた。また、得られた構造が遷移構造であることを確かめるため、振動解析を実行し、虚数振動数が１本のみ存在することを確認した。さらに、得られた遷移構造を初期構造にして、IRC 計算を実行し、その構造が反応前と反応後の構造を結ぶ構造であることを確認した。
【結果と考察】Gaussian09 を用い、真空中、及び連続溶媒近似の水中で求めた G-C、A-T 塩基対と OH ラジカルの反応パスと自由エネルギーの変化を Figure 1 に示す。真空中と水中での結果を比較すると、G-C、
A-T 塩基対の両者に対し、水中の方が反応後の安定化が大きくなり、反応後の構造の存在確率が水和の影響により大幅に上昇する。また、G-C と A-T 塩基対の反応障壁自由エネルギーを比較す
ると、A-T 塩基対の方が真空中では 0.3 kcal/mol、水中では 1.6 kcal/mol 小さくなり、OH ラジカル
は A-T 塩基対とより反応し易いと予測できる。現在、他のラジカルと塩基対との反応機構を解析中であり、それらの結果は当日のポスターにて発表する。

(a) G-C+•OH(b)A-T +•OH
Figure 1 Schematic free energy diagram in vacuum for the reactions between (a) G-C and •OH and (b) A-T
and •OH. The results in water approximated by continuum solvation model are shown in parentheses. Free
energies (kcal/mol) are evaluated by UM06-2x/6-31++G(d,p) at 298.15 K and 1 atmosphere of pressure.
[1] J. Cadet, et al., Free Radical Bio- logical Medicine, 2010, 49, 9.
[2] J. R. Wagner, et al., Accounts of Chemical Research, 2010, 43, 564.
[3] H. Tomita, et al., Radiation and Environmental Biophysics, 1998, 36, 235.
[4] E. Shimizu, et al., Journal of Modern Physics, 2013, 4, 442.
[5] Q. Cheng, et al., Chemistry - A European Journal, 2010, 39, 11848.

4P086
Elongation 法による天然・人工核酸の超効率的電子状態計算とその機能設計（九大院・総理工 1，JST-CREST2） ○折本 裕一 1，Liu Kai1，青木 百合子 1, 2

Highly efficient O(N) calculations of natural and artificial DNAs by elongation method（Kyushu Univ.1，JST-CREST2） ○Yuuichi Orimoto1，Liu Kai1，Yuriko Aoki1, 2【序】 DNA をはじめとする核酸は、バイオ分野に限らず電子材料等のナノテクノロジー領域における応用展開でも期待されており、特に近年、図１a に示す金属錯体型人工 DNA1)のような人為的に化学修飾を施した人工核酸が盛んに合成され、その多様さから新規デバイス応用の観点から注目されている。天然核酸や人工核酸の物性・機能を理解し、さらにはそれらを自在に制御して新規機能を創出するためには電子論的立場から分子構造と電子状態の関係を正確に知ることが本質的な第一歩になると考えられる。
非経験的な量子化学計算の適用が本目的に有効と考えられるが、核酸のような巨大系の電子状態を得ることはその莫大な計算コストにより、計算機が発展した現在でも容易ではない。系のサイズ N の３～４乗に比例した計算コストがかかるため、扱う系が大きくなるとともに電子状態を得るための時間（及びメモリ等の計算資源）が爆発的に増大してしまい、いずれ計算不可能となる。
当研究グループが開発してきた高分子の理論的重合法（３次元 Elongation 法）2)は、巨大系に対してもその電子状態を化学的精度を保持する高精度で、なおかつ N の１乗に比例した計算コスト（オーダーN）
で超効率的に得ることができる。本研究では、Elongation 法による天然核酸、及び人工核酸の電子状態の効率的獲得とその物性抽出を行い、生体材料の機能設計手法としての有効性を検証する。
【方法】 Elongation 法 2)は、計算機上での高分子の重合反応と考えることができ、反応末端に任意のモノ
マーを付加させつつ反応末端近辺のみ解くことで、逐次的に系の電子状態を伸長させていく方法である。
その具体的な手順を、金属錯体型
人工 DNA を例として図 1b に示す。
(a)
最初に、適当なサイズの出発クラス
Ag
タに対する電子状態計算を行い、
2Cu
1
系全体に広がった正準分子軌道~
Frozen Active（CMO）を獲得する。その後、CMO
CMO
RLMO RLMOのユニタリ変換によって、反応末端
(b)
伸長
伸長
側（攻撃ユニットが付加する側）のCu
Zn
Active 領域局在化軌道(RLMO)と、
Cu
Cu
攻撃
攻撃ユニットから遠い側の Frozen
Zn
ユニット
相互作用領域
RLMO を得る。後者は攻撃ユニット
相互作用領域
Frozenと相互作用がないほど離れている
RLMO
ため、計算から外すことができる。
伸長…
Active RLMO と攻撃ユニットからな…
る相互作用領域に対する部分空間
金属錯体型人工DNA の電子状態の Fock 行列を対角化することで、
攻撃ユニット分が伸長した全系の
図１ （a）金属錯体型人工 DNA の例 11a), 21b)、（b）Elongation
電子状態を精度を落とさず得ること
法による電子状態計算の手順ができる。ユニタリ変換による軌道局在化と、反応末端に相当する相互作用領域だけの固有値問題を繰り返すことで、Frozen RLMO が次々と蓄積され、電子状態が伸長される。各ステップで相互作用領域はほぼ一定であり、さらに軌道のtailing が及ばなくなった基底関数を随時計算から外す AO Cutoff 法によって、Fock 行列の作成と変換（AO→RLMO）にかかる時間も一定となる。これにより、巨大系に対しても高精度を保ったまま超高効率な+

2+

オーダーN 計算が達成される。
B-poly(dA)・poly(dT) DNA
伸長方向【 結 果 と 考 察】 核 酸 の電 子 状態 計 算へ のElongation 法の有効性を検証するため、天然
核酸 B-poly(dA)・poly(dT) DNA に対するベン
チマークテストを行った（図２）。グラフ中、従来
攻撃ユニット600000法では、計算時間が系のサイズとともに急激に増大しているのに対し、Elongation 法では AO
500000 各サイズに対する
従来法
CPU時間 (sec)Cutoff 開始以後、各ステップにかかる計算時間
400000がほぼ一定になっていることが分かる。また、計算精度は、全系をダイレクトに解く従来法と比AO cutoff
300000
べて、全エネルギーの誤差（全 E（Elongation
開始
20ユニット目における法）－全 E（従来法））が塩基対 20 ユニット目で従来法との全エネルギー誤差
200000-9
2.59×10-9 Hartree/atom
わずかに 2.59×10 Hartree/atom であった。全
100000系に換算しても 0.002 kcal/mol の誤差であり、Elongation法
高精度に電子状態が得られている。
0
次に、様々な核酸に対して Elongation 法が
4 6 8 10 12 14 16 18 20 22 24 26 28
安定的に電子状態を算出できるか、非周期性
塩基対ユニット数の塩基配列、二重らせんのピッチの違い、カウンターイオンの有無などの様々なモデルを用
図２ Elongation 法の計算効率（天然核酸によるベいて検証を行った（図３）。図の中で、3～6 はランチマークテスト(HF/STO-3G)2b)）
ンダムな配列（T-C-G-C-T-G-T-A-T-G-T-T-G-A-A-C-C-T-C-T）を持っており、ピッチは 6 だけが A 型（他は B 型）を持つ。また、4 と 5 ではカウンターイオ
ン（Na+と K+）を全てのリン酸部に付けた。その計算結果から、全てのケースで塩基対 20 ユニット目における全エネルギー誤差のオーダーが 10-8 Hartree/atom を保持できていることが分かった。これは、１万原子の系を考えても 10-1 ～10-2
kcal/mol 程度の誤差であり、非周期系の巨大な核酸に対しても十分に化学的精度
6
3 B-DNA (random)を満たした解析が可能であ
A-DNA
Na+
(random)
ることが分かった。
このように天然核酸に対して高効率、高精度、加え20ユニット目における従来法との全エネルギー
4 B-DNA (random) + Na++て様々な系の安定計算をK
誤差（Hartree / atom）
確認できた Elongation 法を、
3
2.29×10-8
今後、人工核酸の電子状
44.47×10-8
態計算にも応用し、さらに
5
8.46×10-8は分子構造と物性（導電性、
5 B-DNA (random) + K+
6
3.10×10-8強磁性、非線形光学特性
図３ 様々な核酸（塩基配列：T-C-G-C-T-G-T-A-T-G-T-T-G-A-A-C-C-Tなど）の関係解析を通して-C-T）に対する Elongation 法の計算精度（HF/STO-3G）
天然核酸及び人工核酸の機能設計を展開していく。
1)2)(a) K. Tanaka, A. Tengeiji, T. Kato, N. Toyama and M. Shionoya, Science, 299, 1212 (2003); (b) K.
Tanaka, Y. Yamada and M. Shionoya, J. Am. Chem. Soc., 124, 8802 (2002); (c) G. H. Clever, C. Kaul
and T. Carell, Angew. Chem. Int. Ed., 46, 6226 (2007).
(a) A. Imamura, Y. Aoki and K. Maekawa, J. Chem. Phys., 95, 5419 (1991); (b) Y. Aoki and F. L. Gu,
Phys. Chem. Chem. Phys., 14, 7640 (2012).

4P-087

FATM-iExg 法によるシュレーディンガー解の計算:
含窒素有機化合物への応用( 量 子 化 学 研 究 協 会 , J S T - C R E ST )○ 黒川 悠索，中嶋 浩之，中辻 博

Solving the Schrödinger equation with the FATM-iExg method: Application tonitrogen-containing organic compounds(QCRI, JST-CREST) ○Yusaku I. Kurokawa, Hiroyuki Nakashima, Hiroshi Nakatsuji【緒言】
原子・分子のシュレーディンガー方程式の厳密解に限りなく近い波動関数を求めることが可能な理論として FC 法（Free Complement, 自由完員関数法）が提案されている [1]。order n の FC 波動関数   は order n までの完員関数{φi}の線形結合   n    iMi Cii n  と定義される。ここで完
 n

n

  



員関数系{φi}は、任意の初期関数φ(0) を用いて、 i   1  g  H  E     として得られ、この右辺を展開した時に現れる線形独立な関数からなる。係数{Ci}の決定方法としては積分法（変分nn0法）とサンプリング法の二つが提案されているが、いずれの方法でも Order n が十分に大きいと FC 波
動関数は正確な波動関数に収束する。この方法をそのまま分子系に適用(direct 法)することもできる
が、精度を落とすことなく効率よく分子系の波動関数を表現できることが望ましい。我々はまず有機分子を対象に考え、これに適した波動関数の構築を目指した。詳しい理論については 4E06 及び4E07 を参照されたい。
【方法と結果】
有機分子の構成要素は主に水素、炭素、窒素、酸素である。FATM 法は、個々の構成元素の正確な波動関数から分子全体の波動関数を求めるために開発された。本研究では、まず、窒素と水素の波動関数を LSE サンプリング法によって用意し、その関数を用いて FATM- iExg 法
によって窒素化合物の波動関数を求めた。
まずは窒素原子の波動関数を求めた。窒素原子の初期関数を N0  1s 2 2s 2 2 px 2 py 2 pz とすると、Order が２(n=2)の時、674 個の完員関数が生成された。一般に 1 電子軌道の動径部分は

r n exp   n r  の形をとるが n が大きくなると関数が外側に広がる。今電子密度が大きい領域を記述したいので高い r の指数を持つ 1 電子関数に対しては  n   c n として電子の広がりを抑えた(η=1.04)。ここで、軌道指数αc は Clemenｔi を用いた。さらに、2s 軌道を 1s 軌道と直行させるために、
Table 1. FATM-iExg 法による NH 分子のエネルギー
エネルギー (au)

Exact との差(au)H square error

FATM iExg 法-55.146 055

7.57×10-2

0.835

FATM iExg 法 (+混成配置)-55.162 201

5.95×10-2

0.677

FATM iExg 法 (R=10au)-55.092 357

3.15×10-3

0.097

正確なエネルギーの参照値-55.089 2

Table 2. FATM-iExg 法による N2 分子のエネルギー
エネルギー (au)

Exact との差(au)H square error

FATM iExg 法 (Exchange=1)-109.231 118

3.11×10-1

3.37

FATM iExg 法 (Exchange=2)-109.547 167-4.87×10-3

2.01

正確なエネルギーの参照値-109.542 7

2s : r exp   2 s r   2s ' : 1  cr  exp   2 s r  の置き換えを行った。サンプリング点は原点周りに球対称に Ns=3x106 個とった。この時各サンプリング点における H square error の値が 2  r    H  r    r   E    r    を満たすサンプリング点 rμのみ選択し評価に用いた(ε22=10-8)。新たなサンプリング点のセットで波動関数を求め、再び閾値より小さな H square error を与える点を選択する、という操作を収束するまで行った。すると E=-54.592 319 au (ΔE= -0.00311 au),
σ2= 0.09308 が得られ、波動関数は、

  0.2470 1s 2 2s '2 2 px 2 p y 2 pz  0.2270  r71s 2  2s / r2  2 px 2 p y 2 pz  0.1958r221s 2 2s 2 2 px 2 p y 2 pz 2が得られた(各 cf は規格化・反対称化されている)。サンプリング点でシュレーディンガー方程式が成立しているという条件を満たすように係数を決定しただけにも関わらず、   がメイン0配置となり、化学的直観と合致している。 NH 分子の計算にはこの波動関数を用いた。
次に NH 分子の波動関数を FATM 法により求めた。原子核近傍の記述には上で得た窒素原子と水素原子の波動関数を用い、結合領域の関数は order=1 とした。すると７４７個の完員関数が生成された。Table 1 に示したように、平衡核間距離においてエネルギーは E=-55.146 055au (ΔE= 7.57×
10-2)、R=10au においては E=-55.092 357au (ΔE= 3.15× 10-3)と求まった。次に 2s 電子と 2p 電子の混成の効果を露わに含めるため、初期関数に 2s->2p の励起配置も含めた。するとエネルギーは改善され、E=-55.162 201au (ΔE=5.95× 10-2)となったが、まだ不十分である。これは分極関数やπ結合を表現する関数などが足りないためである。これらも含めた結果など、詳細は当日に発表する。
一般に波動関数は、パウリの原理の要請により全電子を反対称化しなければならない [2]。しかし遠く離れた原子間の複数の電子の交換の寄与は指数関数的に減少する。この原理に基づいた高速
反対称化法を iExg 法と呼ぶ。Table 2 に示したのは、iExg 法を N2 分子に適用した結果である。二つの窒素原子間の電子の交換を 1 つまでにした場合(Exchange=1)、エネルギーは E=-109.231au(ΔE=0.311au)であったが、電子の交換を 2 つまで行うと(Exchange=2) E=-109.5472au (ΔE=-4.87×
10-3)となり、結合エネルギーをより正しく表現できた(窒素原子内の電子の反対称化は全て行っている)。ハイトラー・ロンドンの研究にあるように結合エネルギーのほとんどは交換エネルギーからくるものであり、本結果と合致している。発表当日は計算法の詳細も併せ発表する。
以上のように、本研究では窒素原子の精密な波動関数を求め、FATM-iExg 法を用いることで原子間の結合も正しく描写できることが示された。この方法を一般の有機分子のシュレーディンガー解の計算に応用する。
Reference:[1] H. Nakatsuji, J. Chem. Phys., 113, 2949 (2000), H. Nakatsuji, H. Nakashima, Y. Kurokawa, A.
Ishikawa, Phys. Rev. Lett. 99, 240402 (2007), H. Nakatsuji, Acc. Chem. Res., 45, 1480 (2012).
[2] H. Nakashima and H. Nakatsuji, J. Chem. Phys. 139, 044112-1 (2013).

4P088キラル理論分子技術、キラサクの生体分子への応用（量子化学研究協会）宮原友夫, 中辻博

Theoretical chiral molecular technology, ChiraSac applied to biological molecules（QCRI）Tomoo Miyahara, Hiroshi Nakatsuji【序】円二色性(CD)スペクトルは、一重結合の回転のような低エネルギ
ー自由度や、水素結合・スタッキング
等の弱い相互作用に鋭敏に反映する。SAC-CI 法[1-3]を用いれば、CD
スペクトルの持つこれらの分子情報を、全て克明に解析することができ、予言することもできる[4]。そこで我々は、キラル分子の CD スペクトルの持つ分子情報の理論解析・予測を目的として、SAC-CI 法と Gaussian 中の有用な方法論を統合したキラサクを構築している（図１）[5]。本発表では、キ

図１. キラル理論分子技術「キラサク」

ラサクを生体分子へと応用し、CD ス
ペクトルとその周辺アミノ酸の効果について考察する。
【方法】
本研究では、
Gaussian[6]に搭載されている DFT 法、QM/MM
法、SAC-CI 法を統合して用いた。エストラジオールを DFT 法により構造最適化し、その構造を用いて SAC-CI 法により CD スペクトルを計算し、実験スペクトルと比較検討することにより、キラル分子の溶液中の安定構
造、CD スペクトルの立体配座依存性、相互作用情報などについて考察

図 2. エストラジオール

した。また、QM/MM 法で求めたエストラジオール―エストロゲン複合体の構造をもとに、UV・CD スペクトルを計算し、エストラジオールと周辺アミノ酸の相互作用について考察した。
【エストラジオール】エストラジオール(図 2)の二面角を回転させたときの CD スペクトルの変化を図 3 に示す。計算による最安定構造の SAC-CI CD スペ
クトルは、実験スペクトルと一致しなかった(図 3(a))。一

Table 1. Model of estradiol-estrogen
receptor complex

方、Δ=330º のときの SAC-CI CD スペクトルの結果が、

Model

QM

MM

実験スペクトルと良く一致した(図 3(l))。回転障壁は 2.5

A

EstradiolNonekcal/mol と小さく容易に回転できる。
また Δ=0º と Δ=330º

B

EstradiolAll amino
acids (239)C

Estradiol +18 amino acids

None

D

Estradiol+
18 amino acids

Other aminoacids (221)のエネルギー差は 0.64 kcal/mol とさらに小さい。従って、PCM ではなく、溶媒を顕わに加えると安定構造が変わる
可能性がある。

図 3. エストラジオール(E2)の SAC-CI CD スペクトル(赤)と実験スペクトル(黒:[6])【エストラジオール―エストロゲン複合体】エストラジオール―エストロ
ゲン複合体の UV・CD スペクトルについて 4 種類のモデルを用いて計算した(Table 1)。エストラジオールのみのモデル(A,B)では、アミノ酸の電荷の有無によって CD の符合が変わるなど、電荷が結果に大きな影響を与えている。UV もスペクトルの形状に変化が現れる。一方、エストラジオールを取り囲むアミノ酸 18 個を加えたモデル(C,D)では、アミノ酸の電荷が結果に与える影響は小さい。取り囲む 18 個のアミノ酸がエストラジオールと相互作用し、その性質を決めていて、それ以外のアミノ酸は、エストラジオー
ルの性質に直接影響を与えないことが分かる。
当日は、ロドプシンの cis-trans 異性化と CD スペクトルの関係についても発表する予定である。
【参考文献】

図 4. 蛋白質中の E2 の SAC-CI スペクトル[1] H. Nakatsuji, K. Hirao, J. Chem. Phys. 1978, 68, 2053, H. Nakatsuji, Chem. Phys. Lett.1978, 59, 362.;
1979, 67, 329, 334; Bull. Chem. Soc. Jap. 2005, 78, 1705. [2] M. Ehara, J. Hasegawa, H. Nakatsuji, Theory
and applications of Computational Chemistry, The First 40 Years, Elsevier Oxford, 2005; p1099. [3] SAC-CI
homepage. http://www.qcri.or.jp/sacci/ (16/12/2012). [4] T. Miyahara, H. Nakatsuji, H. Sugiyama, J. Phys.
Chem. A, 117, 42 (2013). [5] T. Miyahara, H. Nakatsuji, J. Phys. Chem. A, submitted. [6] M. J. Frisch, et at.
GAUSSIAN 03, Gaussian, Inc. Pittsburgh PA, 2003. [6] T. Takakuwa, Jasco Report, 1995, 37, 38.

4P089巨大活性空間へ適用可能な多参照配置間相互作用法の開発：DMRG-MRCI(総研大 1 ，分子研 2 ) ○齋藤雅明 1 ，倉重 佑輝 1,2 ，柳井 毅 1,2
Highly scalable multireference configuration interaction theory with internal contractionof density matrix renormalization group wave function(The Graduate Univ. for Advanced Studies1 , IMS2 )Masaaki Saitow1 , Yuki Kurashige1,2 and Takeshi Yanai1,2[序]多参照電子状態理論は, Hartree-Fock 法或は密度汎関数（DFT）法をベースとする単参照理論が苦手とする, 静的電子相関の正確な記述を可能とする. 中でも多参照配置間相互作用（MRCI）法は静的相関のみ
ならず, 動的電子相関の記述にも優れたハイブリッドとして知られている. MRCI は CASPT2 などの多参照摂動論の弱点である, intruder state の存在に対しても数値的に安定であり, 且つ多参照結合クラスター（MRCC）法よりも定式化がシンプルである為に, MRCC に見られる特異性 1 を持たない. MRCI 法は, 基底・励起状態を問わず, 高精度に波動関数を構築する事が可能であり, 結合の解離・生成を伴う反応過程の記述にも適している. 通常, MRCI 法では参照関数からの一電子及び二電子励起のみを考慮する為に, 打ち切られた CI 法としての性質から, 大きさに関する無矛盾性を欠く. ただ, 四電子励起の効果を a priori あ
るいは a posteriori に取り込む事で, この欠点は容易に補正可能である事が知られる.2

MRCI 法は 1970 年代に Bünker らにより, 参照関数からの励起行列式を変分基底として定式化された（非縮約表現）. 3 ただこのアプローチでは, 活性軌道数に対して指数関数的にハミルトニアンの次元が増大する為に, 非常に小規模な系に対してのみ適用可能であった. その後, 1980 年代に Werner らによって,よりコンパクトな変分基底である内部縮約基底を用いた内部縮約型（IC-）MRCI が導出されたものの, セミインターナル励起に関するハミルトニアン行列要素の計算に 5 体の縮約密度行列（RDM）の評価が不可避となる.4 故に, Werner らによる初期 IC-MRCI の計算コストは系のサイズに対して 11 乗のオーダーでスケールされ, 実用的な方法論と呼ばれるに至らなかった. 今現在, 最もスケーラブルであると言われている部分内部縮約型（PIC-）MRCI では, 5-RDM の評価を回避する為に, セミインターナル励起の記述に非縮約表現を用いる.5,6 PIC-MRCI は高効率な計算プログラムに実装され, 小活性軌道系に対しては非常に高速である事が知られているが, 非縮約表現を用いている為に, 活性空間の大きさに対して計算コストは指数関数的に増大する. 従って, 適用可能な系のサイズは限られる. そこで我々は, 非縮約表現を用いない完
全内部縮約型（FIC-）MRCI の導出及び, 計算コードへの実装を行った.

[理論]
我々のアプローチは, セミインターナル励起に関するハミルトニアン行列要素をダイレクトに評価するのでは無く, 多重交換子を経由して導出する事に依って, 5-RDM の出現を数学的に厳密に回避可能である事に立脚する. セミインターナル行列要素に関して次の恒等式が成立する：
000hΨ0 |E SI HESI |Ψ0 i = hΨ0 |[E SI , [H, ESI ]]|Ψ0 i + E0 hΨ0 |E SI ESI |Ψ0 i
ここで, ESI はセミインターナル励起演算子, Ψ0 は参照関数を表す. 5-RDM は右辺第一項の交換子の中で打ち消されて寄与しない. 全体として, ハミルトニアンは 4 体までの RDM を用いて構築可能となるが,我々は更に, 4-RDM を低次の RDM を用いて近似する（キュムラント近似）. RDM の生成には, 従来法で成されているように Full-CI を行うのでは無く, ab initio 密度行列繰り込み群（DMRG）法を用いる.7,8 また, プログラマブルな FIC-MRCI 方程式の導出及び, 計算コードへの実装は, 我々が開発したテンソル生成
ライブラリ FEMTO を用いて行った.9,10 結果として得られる計算コードは, 並列化 API である Message

Passing Interface 及び Global Arrays Toolkit に基づいており, 大規模に並列実行可能である.

[結果]図に我々の FIC-MRCI と, Werner らによる部分内部縮約型（WK-, CW-MRCI）とのポリエン分子（C6 H8 ∼ C24 H26 ）に対する計算時間の比較を示す. 全ての C 2pz 軌道を活性軌道に設定し, H 及び C 1s
軌道は frozen とした. 部分内部縮約型 MRCI 計算には Molpro2012 を用いた. また, FIC-MRCI ハミルト
ニアンの構築には, キュムラント近似を用いた（cu(4)）. WK-及び CW-MRCI は, 活性軌道数が少ない場
合には非常に高速であるが, 活性軌道数の増大に従って, 指数関数的に計算時間が増大する. その一方で,我々の FIC-MRCI は非縮約基底を用いていない為に, 活性空間の増大に対して多項式的にスケールされる.
その結果として, 既存のアプローチが適用不可能なサイズの活性軌道系に対しても, 現実的な計算時間で展開可能である.
精度に関するベンチマークとして, 表にポルフィリン分子の S0 – T0 状態間のエネルギー差を示す. CW-

MRCI 計算では, ポルフィリン分子が持つ 24 の π 軌道の内, 8 つを活性軌道に設定した. 一方で, DMRG をベースとした FIC-MRCI（DMRG-MRCI）計算では, 全ての π 電子及び π 軌道を活性空間に設定した. H,

C, N 1s 軌道は frozen とした. 8 軌道モデルでは, 実測値を 0.2 eV 程度過大評価する事が判る. 一方で, 24
活性軌道の場合では, 大きさに関する無矛盾性を補正する（MRCI+Q 及び ACPF）事で, 実測値との非常に良い一致が見られる事が判る. ここでは示していないが, 我々の DMRG-MRCI の傾向は, 高精度理論である拡散モンテカルロの傾向とも良く一致する事を確認した.10 なお, 今現在, 鉄ポルフィリン（CompoundI）に関する計算を実行中である.
MethodBasis set

S0 / Eh

T0 / Eh∆E(S0 -T0 ) / eV

CW-MRCI

6-31G*-985.264 078-985.199 718

1.75

CW-MRCI+Q

6-31G*-986.023 128-985.958 101

1.77

FIC-cu(4)-MRCI

6-31G*-985.410 545-985.358 663

1.41

FIC-cu(4)-MRCI+Q

6-31G*-986.504 297-986.447 988

1.53

FIC-cu(4)-ACPF

6-31G*-986.581 052-986.524 167

1.55

–

–

–

1.58

CAS(8e, 8o)CAS(26e, 24o) with the DMRG

Experiment[文献]
F. A. Evangelista, and J. Gauss, J.

Chem. Phys. 134, 114102 (2011).
2 P. G. Szalay, T. Müller, G. Gidofalvi,

H. Lischka, and R. Shepherd, Chem.
Rev. 112, 108 (2012).
3 K. J. Bünker, and S. D. Peyerimhoﬀ,

Theor. Chim. Acta. 35, 33 (1974).
4

H.-J. Werner, and E.-A. Reinsch, J.

Chem. Phys., 76 3144 (1982).
5

H.-J. Werner, and P. J. Knowles, J.

Chem. Phys., 89 5803 (1988)
6 K. R. Shamasundar, G. Knizia, and H.-

J. Werner J. Chem. Phys 135, 054101(2011).
7 G.K.-L. Chan, and M. Head-Gordon,

J. Chem. Phys., 116 4462 (2002).

Computational time per iteration / sec.

1

1.0×106
Polyene;CnHn+2 / CAS(ne,no)1.0×105

1.0×104

FIC-cu(4)-MRCI
㻌

1.0×103

CW-MRCI㻌
WK-MRCI㻌

1.0×102

10.0

1.0
6

8

10

8 Y. Kurashige, and T. Yanai, J. Chem.

12
14
16
18
20
22
Number of the active orbitals (n)Phys., 130 234114 (2009).
9 FEMTO ::An Integrated Toolset

For The Automated Tensor Generation, Version 0.1.0, M. Saitow
10 M. Saitow, Y. Kurashige, and T. Yanai, J. Chem. Phys., 139 044118 (2013).

24

4P090Density matrix renormalization group method for prediction of
hyperfine structure
Tran Nguyen Lan,1 Yuki Kurashige,1,2 and Takeshi Yanai1,2
1

The Graduate University for Advanced Studies, Okazaki, Aichi 444-8585
2

Institute for Molecular Science, Okazaki, Aichi 444-8585
E-mail: lantran@ims.ac.jp

Despite all recent progress, the accuracy of hyperfine coupling constants (HFCCs) is stillchallenging for quantum chemistry. In this work, we have performed the density matrix
renormalization group (DMRG) calculations to predict the HFCCs of several π−radicals: BO,
CO+, CN, and AlO. The DMRG algorithm used herein was considered as a near fullconfiguration interaction (FCI) method and the electron correlation effects weresystematically investigated using the complete active space (CAS) procedures, i.e CASSCF
and CASCI. The geometry of radicals has been adapted from experiments. The EPR-III basis
set was used for C, O, and N elements, while the IGLO-III basis set was used for Al element.
Although the HFCCs of BO and CO+ calculated by conventional methods, i.e. density
functional theory (DFT) and coupled cluster (CC), were already in the acceptable agreementswith experimental values, our DMRG calculations have pushed the accuracy of HFCCs
further. For CN radical, our DMRG calculation has captured the high-order effect ofcorrelation treatment. Regarding AlO radical, our DMRG calculation properly described themulti-reference problem resulting from the balance between ionic states. Therefore, the
HFCCs of C center in CN and Al center in AlO are in the good agreements with experimental
values.
Concerning the CAS-type procedure, most CASCI calculations provided the poor agreementswith experimental values, while the quite accurate HFCCs of less electronegative centers canbe found when the orbital optimization procedure was employed with the same CAS. Thesesituations can be attributed to the nature of one-particle bases in which the DMRG
calculations were performed. The canonical Hatree-Fock (HF) orbitals were used as theorbital bases for the DMRG-CASCI calculation. Whereas, the DMRG calculation in
conjunction with CASSCF procedure was performed on much more compact orbitals whichwere obtained through orbital optimization procedure.
In summary, our results have confirmed that it is critical to correlate the core electrons tocorrectly obtain the spin density at the nucleus. Additionally, the inclusion of polarization
shells is necessary for describing the dynamical correlation, which makes the proper spinpolarized effects. Finally, the assessment for the convergence of HFCCs with respect to thenumber of renormalized states M indicates that the acceptable result can be obtained even if
total energy is not converged with respect to M.

Table 9: Hyperfine coupling constants (in MHz) for 2 AlO molecule. The IGLO-III and EPR-III
basis sets were used for Al and O, respectively. The total number of AOs is of 84. For comparison,
we also adapted the CCSD and CCSD(T) results calculated by Kossmann and Neese.
27 Al

Methods

17 O

A(K;c)

A(K;d)
11

A(K;d)
22

A(K;d)
33

A(K;c)

A(K;d)
11

A(K;d)
22

A(K;d)
33

DMRG-CASCI(21e,36o)

710.91

46.45

46.62

93.06

1.82

44.98

44.89

89.87

DMRG-CASSCF(21e,36o)

722.73

54.19

54.09

109.28

15.06

50.73

50.43

101.16

B3LYP

512.21

59.97

59.97

119.93

8.17

66.22

66.22

132.43TPSS

656.79

56.10

56.10

112.21

9.52

59.91

59.91

119.83

BP

653.71

56.86

56.86

113.72

14.21

59.60

59.60

119.20

CCSD

482.40

57.20

57.20

114.30

18.10

63.80

63.80

127.70

CCSD(T)

565.30

56.20

56.20

112.40

19.30

58.90

58.90

117.80Exp

gas-phase

738

56

56

112

Exp

Ne-matrix

766

52

52

104

50

100

n/a
2

50

HFCCs (in MHz) for AlO radical. The IGLO-III and EPR-III basis sets were used for Al
and O, respectively. The total number of AOs is of 84.
Table 11: The convergence of total energy, Fermi contact term, and spin-dipole term at C centerwith respect to number of renormalized states M. For the spin-dipole term we only present the firstelement. The active space is CAS(13e,30o). The convergent MOs of DMRG-CASSCF(13e,30o)
calculation with M = 512 was used as the reference for the other DMRG calculations 0.7
and onlyAl
DMRG-CASSCF(15e,33o)
CASCI procedure was performed for these calculations.
0.6
Energy (Eh )A(K;c) (MHz)A(K;d)
11 (MHz)128

92.570835

576.8053.18

25692.575539

563.7952.97

51292.577510

561.95

52.92

1024

92.578114

559.4452.86

0.5

ρα -β (a.u.)M

DMRG-CASSCF(21e,31o)
DMRG-CASSCF(21e,36o)

0.4
0.3
0.2

O

0.1

The convergence of total energy, Fermi
0contact term, and spin-dipole term at C
-1
0
1
2
3
4
center (CN) with respect to number ofz-axis (a.u.)
renormalized states M. The active space is↵Figure1: Spin density distribution of SNO with largest ni of AlO radical. The results were calcuCAS(13e,30o). The EPR-III basis set
wasSpin density distribution of spin natural
lated by DMRG-CASSCF(15e,33o), DMRG-CASSCF(21e,31o), and DMRG-CASSCF(21e,36o).
used.	
  orbital (SNO) with largest spin
occupation number of AlO radical. The
geometry of AlO radical (in Bohr) are:O(0.000, 0.000, 0.000) and Al(0.000,
0.000, 3.057)

34

37

4P091
相対論効果と電子相関効果を同時に取り込んだ大規模計算法の開発(早大先進理工 1，早大理工研 2，JST-CREST3，京大 ESICB4)○中野匡彦 1，清野淳司 1，中井浩巳 1-4

Large-scale Computational Method with Simultaneous Inclusion ofRelativistic and Electron Correlation Effects(Advanced Science and Engineering, Waseda Univ.1, Rise, Waseda Univ.2, JST-CREST3, ESICB, Kyoto Univ. 4)○Masahiko Nakano1, Junji Seino1, Hiromi Nakai1-4【序】周期表の全元素を扱える高精度な量子化学計算を実現するためには，相対論効果と電子相関効
果を同時に考慮した方法論の開発が必須となる．そこで我々は，高精度な 2 成分法の一つである無限次 Douglas-Kroll-Hess (IODKH) 法に基づく多電子ハミルトニアン (IODKH/IODKH ハミルトニアン)[1]に注目して理論開発を行ってきた．その結果，IODKH/IODKH ハミルトニアンを用いた Hartree-Fock(HF) および電子相関計算において，原子番号 100 を超える超重元素に対しても，4 成分法の結果を高精度に再現できることを確認してきた[1,2]．さらに，相対論効果の局所性を利用した局所ユニタリー変
換 (LUT) 法を開発し，相対論的ハミルトニアン生成に要する計算時間を大幅に短縮するとともに，ス
ピンフリーレベルでは非相対論と同等のコストで計算できることを示してきた[3-5]．
本研究では，スピン－軌道相互作用項などを含んだスピン依存ハミルトニアンに基づき，主要な相対論効果を全て考慮した大規模計算法を提案する．そのためには，①スピン依存ハミルトニアン生成の高速化，②あらゆるスピン依存演算子を表現できる HF および電子相関理論への拡張，③大規模計算のための分割統治 (DC) 法を統合した理論の開発が必要となる．①と②の開発はこれまで独立的に報告しているため[2-5]，本研究では，これらの理論を③と統合可能な形式へ拡張した．
【理論】本研究では，IODKH/IODKH ハミルトニアンに含まれるスピン依存項による𝛼𝛼/𝛽𝛽スピンの混成の効果を適切に表現するために，一般化スピン軌道 (GSO) を用いた一般化非制限 HF (GUHF) 法およ
び GUHF 波動関数を参照とした電子相関法を基盤とする．本稿では，電子相関法についてのみ述べる．
GSO を用いた DC 電子相関法では，通常法と同様に，相関エネルギー𝐸𝐸corr を持つ分子全体を複数の

𝑠𝑠
小さな部分系𝑠𝑠に分割して，相関エネルギー𝐸𝐸corrの計算を行う[6-8]．

𝐸𝐸corr ≈{subsystem}

�
𝑠𝑠

𝑠𝑠
𝐸𝐸corr(1)

𝑠𝑠
ここで，GUHF 波動関数を参照とした一般化非制限 MP2 (GUMP2) 相関エネルギー𝐸𝐸corr
は，
𝑠𝑠
𝐸𝐸corr={occ(𝑠𝑠)} {vir(𝑠𝑠)} {𝛼𝛼,𝛽𝛽} 𝒞𝒞(𝑠𝑠)

1
�
2
𝑖𝑖𝑖𝑖

∗

𝜔𝜔
� � ��𝐶𝐶𝜇𝜇𝜇𝜇
� �𝜇𝜇𝜇𝜇�𝑔𝑔�𝑖𝑖𝑖𝑖 �𝑗𝑗𝑗𝑗��𝑡𝑡̃𝑖𝑖𝑖𝑖,𝑎𝑎𝑎𝑎 − 𝑡𝑡̃𝑖𝑖𝑖𝑖,𝑏𝑏𝑏𝑏 �
𝑎𝑎𝑎𝑎

𝜔𝜔

𝜇𝜇

(2)

と表される．ただし，𝑔𝑔�𝑖𝑖𝑖𝑖 は IODKH/IODKH ハミルトニアンの二電子項， 𝒞𝒞(𝑠𝑠)は部分系𝑠𝑠に属する基底
̃
関数の集合，𝑡𝑡𝑖𝑖𝑖𝑖,𝑎𝑎𝑎𝑎= �𝑎𝑎𝑎𝑎�𝑔𝑔�𝑖𝑖𝑖𝑖 �𝑏𝑏𝑏𝑏���𝜀𝜀𝑖𝑖 + 𝜀𝜀𝑗𝑗 − 𝜀𝜀𝑎𝑎 − 𝜀𝜀𝑏𝑏 �は二電子励起振幅，𝐶𝐶および𝜀𝜀はそれぞれ GUHF 法から得られる MO 係数および軌道エネルギーである．

式(2)に含まれる MO 積分�𝑖𝑖𝑖𝑖�𝑔𝑔�𝑖𝑖𝑖𝑖 �𝑗𝑗𝑗𝑗�は，非相対論と同様に AO 積分�𝜇𝜇𝜇𝜇�𝑔𝑔�𝑖𝑖𝑖𝑖 �𝜌𝜌𝜌𝜌�の積分変換によって与えられる．ただし，IODKH/IODKH ハミルトニアンには Pauli 行列𝝈𝝈を含む様々な形式のスピン依存演

算子が含まれており，それぞれの形式に応じた異なる積分変換の様式を要する．例えば，スピン－軌道相互作用に対応する演算子𝝈𝝈𝑖𝑖 ⋅ 𝜴𝜴𝑖𝑖 の変換は，{𝛼𝛼,𝛽𝛽} {AO}1𝜔𝜔 ∗ 𝜔𝜔′𝜏𝜏 ∗ 𝜏𝜏′(𝑖𝑖𝑖𝑖|𝝈𝝈𝑖𝑖 ⋅ 𝜴𝜴𝑖𝑖 |𝑗𝑗𝑗𝑗) =� � �𝐶𝐶𝜇𝜇𝜇𝜇� 𝐶𝐶𝜈𝜈𝜈𝜈 �𝐶𝐶𝜌𝜌𝜌𝜌� 𝐶𝐶𝜆𝜆𝜆𝜆 (𝜇𝜇𝜇𝜇|[𝜔𝜔|𝝈𝝈𝑖𝑖 ⋅ 𝜴𝜴𝑖𝑖 |𝜔𝜔′ ]|𝜌𝜌𝜌𝜌)𝛿𝛿𝜏𝜏𝜏𝜏′
2 ′ ′

(3)

𝜔𝜔𝜔𝜔 𝜏𝜏𝜏𝜏 𝜇𝜇𝜇𝜇𝜇𝜇𝜇𝜇と表される．ここで， 𝜴𝜴は IODK 法に固有の軌道角運動量に相当する演算子である．

さらに，式(3)に現れる IODKH/IODKH に特有の�𝜇𝜇𝜇𝜇�𝑔𝑔�𝑖𝑖𝑖𝑖 �𝜌𝜌𝜌𝜌�については，AO 積分�𝑝𝑝𝑝𝑝�𝑔𝑔�𝑖𝑖𝑖𝑖 �𝑟𝑟𝑟𝑟�のユニタ

リー変換 (相対論的変換) によって生成される．相対論効果の局所性に着目した LUT 法の下では，
�𝜇𝜇𝜇𝜇�𝑔𝑔�𝑖𝑖𝑖𝑖 �𝜌𝜌𝜌𝜌� = �

𝑈𝑈𝜇𝜇∗𝐴𝐴𝑝𝑝𝐴𝐴 𝑈𝑈𝜈𝜈𝐵𝐵𝑞𝑞𝐵𝐵 �𝑝𝑝 𝐴𝐴 𝑞𝑞 𝐵𝐵 �𝑔𝑔�𝑖𝑖𝑖𝑖 �𝑟𝑟 𝐶𝐶 𝑠𝑠 𝐷𝐷 �𝑈𝑈𝜌𝜌∗𝐶𝐶 𝑟𝑟 𝐶𝐶 𝑈𝑈𝜆𝜆𝐷𝐷𝑠𝑠𝐷𝐷
�𝜇𝜇 𝐴𝐴 𝜈𝜈 𝐵𝐵 �

1 𝐶𝐶 𝐷𝐷�𝜌𝜌 𝜆𝜆 �
𝑟𝑟𝑖𝑖𝑖𝑖(𝐴𝐴 = 𝐵𝐵 = 𝐶𝐶 = 𝐷𝐷)(the others)(4)と定義される．ここで， 𝑈𝑈はユニタリー変換行列要素， 𝐴𝐴, 𝐵𝐵, 𝐶𝐶, 𝐷𝐷は基底関数の属する原子を表す．すなわち，相対論的変換は一中心積分にのみ適用し，複数原子にまたがる二電子積分は非相対論的な
Coulomb 項で近似する．式(4)のように，ユニタリー変換は系内の各原子に対して独立に定義されるため，相対論的変換の計算コストは，系内に含まれる原子数に対して線形となっている．さらに式(1)のように部分系ごとの電子相関エネルギーを定義することで，電子相関計算全体に要する計算コストの線形スケーリングを達成できる．
【結果】Fig. 1 に， (HF)n (n = 10, 20, …, 100) の(DC-)GUHF 計算および (DC-)GUMP2 計算に要した (n 4.74 ) (n1.16 )CPU 時間を示す．この際，LUT-IODKH/IODKH ハミル
トニアンを用い，DC-GUHF 法のバッファ領域は左右 5分子，DC-GUMP2 法のバッファ領域は左右 2 分子とした．また，基底関数には，DK3 basis set [9]を採用した．
GUHF および GUMP2 計算のコストはそれぞれ𝒪𝒪(𝑛𝑛3.47 ) (n3.47 ) (n1.76 )お よ び 𝒪𝒪(𝑛𝑛4.74 ) で あ る の に 対 し ， DC-GUHF お よ び

DC-GUMP2 計算のコストはそれぞれ𝒪𝒪(𝑛𝑛1.76 )および𝒪𝒪(𝑛𝑛1.16 )となった．さらに，DC-GUHF 法の計算コスト

は，高速多重極法 (FMM) を採用することで𝒪𝒪(𝑛𝑛)まで削減することができる．すなわち，LUT-IODKH/IODKH-

DC 法を用いることで，post-HF 法を含めた全ての計算過程において，高次相対論効果を考慮した線形スケー

Fig. 1. CPU time for calculating GUHF andGUMP2 energies of zig-zag (HF)n chains with andwithout DC method (in second).

リングの達成が可能となる．
【参考文献】[1] J. Seino, M. Hada, Chem. Phys. Lett. 461, 327 (2008). [2] 中野匡彦, 清野淳司, 中井浩巳,
日本化学会第 92 春季年会, 慶應義塾大学 (2012). [3] J. Seino, H. Nakai, J. Chem. Phys. 136, 244102 (2012).
[4] J. Seino, H. Nakai, J. Chem. Phys. 137, 144101 (2012). [5] J. Seino, H. Nakai, J. Chem. Phys. 139, 034109(2013). [6] M. Kobayashi, Y. Imamura, H. Nakai, J. Chem. Phys. 127, 074103 (2007). [7] R. Zaleśny, M. G.
Papadopoulos, P. G. Mezey, J. Leszczynski, in Linear-Scaling Techniques in Computational Chemistry andPhysics: Methods and Applications, (Springer, Dordrecht, 2011), pp.97-127. [8] M. Kobayashi, H. Nakai, Phys.
Chem. Chem. Phys. 14, 7629 (2012). [9] T. Tsuchiya, M. Abe, T. Nakajima, and K. Hirao, J. Chem. Phys. 115,
4463 (2001).

4P092重元素化合物に対する相対論的構造最適化計算*(早大先進理工 , 早大理工研**, JST-CREST***, 京大 ESICB****) ○中嶋

裕也*, 清野

淳司*, 中井浩巳*,**,***,****

Relativistic geometry optimization for heavy-element compounds(*Waseda Univ., **JST-CREST, ***ESICB, Kyoto Univ.) ○Yuya Nakajima*, Junji Seino*, Hiromi Nakai*,**,***【緒言】相対論効果は、量子化学計算を用いた重元素化合物の構造決定において、重要な要素である。この効果を簡便に考慮できる手法として、相対論的なポテンシャルを用いた有効内殻ポテ
ンシャル（RECP）法やモデルコアポテンシャル（MCP）法がある。これまでにこれらの手法に基づいた重元素化合物の構造最適化が広く行われてきた。しかし、ポテンシャル法による相対論効
果の取り扱いの適用限界は明らかでない。これを明らかにするには、高信頼性を有する相対論的計算手法を用いた重元素化合物に対する系統的な数値検証が不可欠である。我々はこれまで無限
次 Douglas-Kroll-Hess（IODKH）法[1]に局所ユニタリー変換（LUT）を組み合わせることで、高精度かつ高速なエネルギー計算手法を開発してきた[2, 3]。さらにスピン非依存（SF）ハミルトニア
ンを用いた解析的微分法の開発も行い、数値検証の結果、4 成分相対論法と同程度の精度かつ非相対論と同等の計算時間で計算可能な手法であることを示してきた[4]。本研究では様々な軽・重元素を含み、種々の結合様式を有する化合物に対して、LUT-IODKH 法に基づく構造最適化計算を行う。これにより構造における相対論効果に関して詳細な知見を得ることを目標とする。さらに重元素ではスピン依存（SD）効果も重要であると予想されるため、SD ハミルトニアンに対する解析的微分計算手法の開発も行う。この手法を用いた系統的な数値検証により、SD 効果の構造に対する寄与を明らかにする。
【理論】SD-IODKH ハミルトニアンに対する Hartree-Fock エネルギーの解析的微分表現は以下のような形になる。
 ,   ||  '  '  V  '  1  ,  '  (1)E  ,
 ' h '  ' ' S   nuc    i D  D
D
D 

 





 RA  , '  , iRA  , '  ,RA RA  2  , ', , '  , , ,  RA 

ここで E は全エネルギー、RA は原子 A の核座標、ω, ω’, τ, τ’はスピン関数、D は密度行列、h は 1
電子ハミルトニアン、 は空間軌道、Vnuc は核間反発エネルギー、εi は軌道エネルギー、S は重なり行列を表す。式(1)は SD ハミルトニアンに基づいているため、スピンに関する制限を完全に取り除き、 および  スピン関数の総和と空間軌道の総和で表される。一方、SF ハミルトニアンの場合は、スピンを変化させる演算子が含まれないためスピンの方向を制限することが可能であり、空間軌道の総和のみで表される。また、LUT-IODKH 法において、式(1)の 1 電子ハミルトニアン
 

行列の微分 h / RA は次のようになる。

 A h LUT
B
2
RA

 
 A  VNNR B


R
N A
 A
 

 A T NR   VNNR   VNSF   VNSD BRN A, BN A , B
N A , B A  A T NR   VNNR B


R
N
 A

 A  B
 A  B, RAB   
 A  B, RAB   (2)ここで T NR は非相対論（NR）の運動エネルギー、 VNNR は NR の核－電子ポテンシャル、 VNSF および
VNSD は相対論的に変換された核－電子ポテンシャル（SF：SF 部分、SD：SD 部分）、A, B は部分系

（ここでは特に原子とする）、χ は基底関数、 RAB は原子 AB 間の距離を表す。また、τ は距離の閾値で相対論的に変換する範囲を決定し、3.5 Å とした。ここで、式(2)の 2 行目の VNSD が 0 の場合にSF 計算となる。
【結果と考察】Table 1 に第 4 周期から第 6 周期までの遷移金属を含む金属錯体の構造パラメータを周期ごとに示す。X-Y は原子 X と原子 Y の結合距離を表す。本研究ではハミルトニアンとして
非相対論（NR）、RECP（SDD）、MCP、SF-LUT-IODKH を用い、計算レベルは B3LYP とした。
NR および LUT-IODKH では基底関数として非縮約型 DKH3-Gen-TK/NOSeC-V-TZP を用いた。また、RECP では第 3 周期までの元素の基底関数として cc-pVTZ を用いた。Δr は NR からの差分、
つまり相対論効果の結合距離への影響を表す。この結果から、MCP 以外、相対論効果を考慮した手法ではすべて結合距離が短くなり、族番号の増加とともに相対論効果の影響が大きくなる。特に 12 族（Zn, Cd, Hg）では、同周期の他の族と比べて顕著に結合距離が短くなる。また RECP、
MCP、LUT-IODKH では、各錯体に対して相対論効果の寄与の程度が異なる。すなわち第 4 周期元素では、結合距離の差分が 0.005 Å 程度とほぼ一致していることから、相対論効果の寄与はい
ずれの手法でも同程度である。第 5 周期元素では、族番号が大きくなるにつれて MCP と
LUT-IODKH の差分が大きくなる一方、RECP と LUT-IODKH の差分がほぼ一定である。また第 6周期元素では、MCP と LUT-IODKH の結合距離の差分が 0.01 Å 程度であり、RECP と LUT-IODKHの結合距離の差分が最大で 0.005 Å である。
当日はこれらの結果を結合様式の観点からも考察し、どのように相対論効果が影響しているのかを議論する。さらに系統的に変化させた分子種を用いて、SD 効果の構造に対する寄与についても議論する予定である。
Table 1. Geometry parameters for transition metal complexes calculated by NR, RECP, MCP and LUT-IODKH at theB3LYP level.
RECP
MCP
LUT-IODKH
Period Group Complex Parameter
NR
r
Δr1)
r
Δr1)
r
Δr1)
4
4
TiOCl2
Ti-O
1.598
1.593
(-0.005)1.601 ( 0.003)
1.595(-0.003)
Ti-Cl
2.239
2.232
(-0.007)2.237 (-0.002)
2.236(-0.003)8
Cr(CO)6
Cr-C
1.926
1.920
(-0.006)
1.931 ( 0.005)
1.920(-0.006)
12
ZnMe22)
Zn-C
1.958
1.945
(-0.013)1.933 (-0.025)
1.942(-0.016)
5
4
ZrBr4
Zr-Br
2.495
2.495
( 0.000)
2.494 (-0.001)
2.488(-0.006)
8
RuO4
Ru-O
1.693
1.682
(-0.011)1.687 (-0.006)
1.681(-0.012)12
CdEt23)
Cd-C
2.201
2.153
(-0.048)
2.171 (-0.030)
2.158(-0.043)67
Re(CO)5Br Re-C (Eq.)4) 2.077
2.023(-0.054)2.004 (-0.073)
2.022(-0.055)
Re-C (Ax.)5) 1.995
1.944(-0.051)
1.929 (-0.066)
1.944(-0.051)
Re-Br
2.703
2.656
(-0.047)
2.633 (-0.070)
2.652(-0.051)
8
OsO4
Os-O
1.737
1.697
(-0.040)1.684 (-0.053)
1.700(-0.037)
12
HgMe22)
Hg-C
2.247
2.123
(-0.124)
2.111 (-0.136)
2.118(-0.129)1)

Difference in bond length from NR. 2) Me = CH3. 3) Et = CH2CH3. 4) Equatorial. 5) Axial.

【参考文献】[1] M. Barysz and A. J. Sadlej, J. Chem. Phys. 116, 2696 (2002). [2] J. Seino and H. Nakai, J. Chem. Phys. 136, 244102.
(2012) [3] J. Seino and H. Nakai, J. Chem. Phys. 137, 144101 (2012). [4] Y. Nakajima, J. Seino, and H. Nakai, APCTCC6,
PS5(S-B), July, 2013.

4P093
局所応答分散力法の改良と汎関数依存性の検証
1(早大先進理工 ，早大理工研 2，JST-CREST3，京大 ESICB4) ○五十幡 康弘 1，中井 浩巳 1-4
Refinement of local response dispersion method and its functional dependence(Advanced Science and Engineering, Waseda Univ.1, RISE, Waseda Univ.2, JST-CREST3, ESICB4)
Yasuhiro Ikabata1, Hiromi Nakai1-4【緒言】
近年，密度汎関数理論による弱い相互作用の記述を目的として様々な分散力補正法が提案されている。当研究室で開発してきた局所応答分散力(LRD)法[1-3]は，系の電子密度から分散力係数を
非経験的に見積る分散力補正法である。LRD 法を LC-BOP 汎関数と組み合わせると，水素結合，π/π スタッキング，ラジカルダイマー[2]，励起状態[3]などにおける弱い相互作用を高精度に再現できる。しかし，LRD 法による高次の分散力係数は過大評価の傾向が確認されている[1]。さらに，
LC-BOP 汎関数との組み合わせでは，分子内の弱い相互作用や電荷を有する系に対する精度が不十分であると報告されている[4,5]。そこで本発表では，分散力係数をより正確に見積るための改良を示し，それに基づいて汎関数依存性を検証した結果を発表する。
【理論】
LRD 法による分散力補正エネルギーは，二原子間の寄与の和として次式のように与えられる。
(n)
Edisp = − ∑∑ Rab− n Cnab [α ] f damp
( Rab , [α ])(1)n a <b

ここで，n = 6, 8, 10, …であり，Rab は原子間距離， Cnab は分散力係数，fdamp はダンピング関数である。分散力係数は分子中の有効原子分極率 α に依存する。

α lma (iω ) = ∫ drwa (r ) wa' (r )

ρ (r )
ω (r ) + ω 2
20

× ∇Rlm (r ) ⋅ ∇Rlm (r + R a − R a' )(2)

ここで l = 1, 2, …，− l < m < l であり，wa は空間分割関数， Rlm は体球調和関数である。ω0 はプラ
ズマ分散関係の長波長極限に対応する。これまでは van Voorhis らによる式[6]13ωq2=0 (r ) = kF2 (1 + λs 2 )2(3)を採用していた。ここで k F = (3π 2 ρ )1 / 3 ， s = ∇ρ (2kF ρ ) であり，λ は経験的パラメータである。また，ダンピング関数は

 n−4
(n)
(Rab / R [α ])−6 f damp( Rab , [α ]) = exp −
2



(4)

としていた。本研究では，高次の分散力係数の過大評価を防ぐために，次式のように ω0 を l に依存させる形に変更した。

ω

2
q =0

1  1+ l 2 (r ) = k F2 1 +
λs 
3 
2


2

(5)

C6 係数は l = 1 の結果から得られるので式(3)と式(5)は同一の結果となる。この変更に伴い，ダンピング関数を n に依存しない形，すなわち式(4)にて n = 1 としたものに変更した。

【数値検証１：分散力係数】
数値検証１：分散力係数】
希ガス二量体に対する LRD 法による C6 係数の絶対誤差が最小になるよう λ を決定した，高次の分散力係数を式(3)，式(5)に基づいて計算した。結果を Table 1 に示す。式(5)を採用することで，
文献値[7-9]に近い C8，C10 係数を与えることに成功した。
Table 1. Dispersion coefficients of rare-gas dimers obtained by the LRD method.
The LC-BOP functional and aug-cc-pVQZ basis set were used.

C 6 (a.u.)

He2
Ne2
Ar2
Kr2

LRD

Ref.

1.558
6.336
61.21
119.2

1.461
6.383
64.30
129.6

C 8 (a.u.)
LRDconventional present
31.82
16.77
174.93 100.85
3063
1931
7436
4771

Ref.
14.11
90.34
1623
4187

C 10 (a.u.)
LRDconventional present
594.1
155.8
4388
1355
13467649780
401372 154672

Ref.
183.6
1536
49060
155500【数値検証２：相互作用エネルギー】
弱い相互作用のベンチマークセットである S22 に対して，LRD 法をいくつかの長距離補正汎関数と組み合わせて相互作用エネルギーを計算した。ダンピング関数に含まれるパタメータは希ガ
ス二量体の相互作用エネルギーを再現するようフィッティングした。22 種類の系のそれぞれについて，CCSD(T)完全基底極限からの差分を Figure 1 に示す。LC-BOP 汎関数を改良した LC-BOP12が特に良い結果となったが，CAMB3LYP，ωB97X も MP2 より平均絶対誤差が小さい結果となった。また，分散力主体の系では組み合わせる汎関数ごとに大きな差が生じた。当日は分子内相互Difference from reference (%)__

作用や電荷を有する系の検証結果も発表する。
60
40

Hydrogen-bonded

Dispersiondominated

Mixed

20
0
-20
-40

LC-BOP12+LRD-60

CAMB3LYP+LRD-80

ωB97X+LRD-100

MP2/CBS

Figure 1. Percentage differences from the reference interaction energies of theS22 set. The 6-311++G(2d,2p) basis set was used for the DFT calculations.
[1] T. Sato and H. Nakai, J. Chem. Phys. 131, 224104 (2009).
[2] Y. Ikabata and H. Nakai, Chem. Phys. Lett. 556, 386 (2013).
[3] Y. Ikabata and H. Nakai, J. Chem. Phys. 137, 24106 (2012).
[4] S. Steinmann and C. Corminboeuf, J. Chem. Theory Comput. 7, 3567 (2011)[5] S. Steinmann and C. Corminboeuf, J. Chem. Theory Comput. 8, 4305 (2012).
[6] O. A. Vydrov and T. Van Voorhis, J. Chem. Phys. 130, 104105 (2009).
[7] Z. C. Yan, J. F. Babb, A. Dalgarno, and G. W. F. Drake, Phys. Rev. A 54, 2824 (1996).
[8] A. Kumar and W. J. Meath, Mol. Phys. 54, 823 (1985).
[9] A. J. Thakkar, H. Hettema, and P. E. S. Wormer, J. Chem. Phys. 97, 3252 (1992).

4P094ペア行列を用いた平均場電子状態計算による静的電子相関の取り扱い（早大・高等研）○小林正人

Static electron correlation treated by mean-field electronic structure
calculations with pairing matrix(WIAS, Waseda Univ.) ○Masato Kobayashi【緒言】
密度汎関数理論や Hartree-Fock (HF)波動関数を出発点とする単参照電子状態理論は、多くの化学現象や化学反応機構を説明し、成功を収めてきた。しかし、軌道の擬縮退に伴って生じる静的電子相関を取り込むことができない。このような電子相関は、遷移金属を含む分子や遷移状態・解離系などの平衡構造から離れた構造では本質的に大きな寄与をもつ。静的電子相関を取り込む手法としては、CASSCF 法や電子対波動関数（ジェミナル）を用いる方法[1,2]のような多配置 SCF 法をベースとした理論が現在の主流である。しかし最近では、
Hartree-Fock-Bogoliubov (HFB)法がジェミナル積波動関数と形式的に類似した波動関数を与えることから、HFB 法で導入されるペア行列を用いた静的電子相関計算に注目が集まっている[3,4]。本研究では、そのひとつである ζ-corrected HF (ζ-CHF)エネルギー汎関数を用いた平均場電子状態計算の性能を評価し、また解析的エネルギー勾配を導出したので報告する。
【ζ-CHF エネルギー汎関数とその性能】
超伝導準粒子の平均場理論である HFB 法や、密度行列汎関数の立場から導出された CHF エネ
ルギー汎関数の表式に基づいて、Staroverov と Scuseria は ζ-CHF エネルギー表式を提案した[3]。
閉殻系の ζ-CHF エネルギーは、空間軌道を用いて次式で与えられる。


*
(1)
ECHF ( )    2h P    2       P P     K 
K    
ここで P   a†a  は一体の密度行列であり、 K    a a  はペア行列と呼ばれる。ζ

は適当なパラメータであり、ζ = –1 のときは HFB エネルギーと、ζ = 1 のときは CHF エネル
ギーと一致する。このエネルギーを P と K に関して変分的に求めるのが ζ-CHF 法である。これらの行列要素は、HFB 法と同様、以下の-77.80方程式を解いて得られる準粒子軌道から例として、図 1 にエチレン C2H4 分子の二重結合回転に対するポテンシャルエネ-77.85-77.90

Energy (Eh)

求められる。
 F
 Δ   X   S 0  X  ε 0  (2)
 †

 

* 
  Δ F   Y   0 S  Y  0 ε -77.95-78.00ルギー曲線を示す。HF 法では二面角が 90°-78.05のところで微分不可能となることが知ら-78.10

れており、CASSCF 法などで静的電子相関

HF
CHF
CASSCF(2,2)
-CHF ( = 0.90)
-CHF ( = 0.75)

0

30

60
90
120
Dihedral angle (°)

150

180

を取り込むことによりこれが改善される。 Fig. 1. Potential energy curves for C H rotation obtained
2

4

by HF, CASSCF, and (ζ-)CHF calculations (6-311G**).

CHF 法は、平衡構造付近では HF 解を与えるが、構造が離れるにつれて HF 法よりも低いエネル
ギーを与える。90°で滑らかな曲線が得られるものの、静的電子相関が過剰に見積もられ、エネル
ギー障壁が正しく再現されない。ほどよい ζ の値を用いる事により、90°付近の記述のみを適切に改善することができることが確認された。ペア行列を用いた手法によるエネルギーの改善法としては、ほかに HFB 波動関数に対して射影演算子を作用させる方法が提案されている[5]。

【ζ-CHF 法に対するエネルギー勾配】(1)式の ζ-CHF エネルギーを原子座標で微分することにより、エネルギー勾配を以下のように表せることがわかる。
P
K  
ECHF ( )
 (H-F force)  2  F    

Q
Q
Q 
 

(3)

ここで第 1 項は Hellmann–Feynman (H-F)力であり、第 2 項が原子中心の基底関数展開により生じる Pulay 力である。(2)の方程式と準粒子軌道の直交条件を用いると、(3)式は以下のような簡単な式で表せることがわかった。
S
ECHF ( )
 (H-F force)  2W 
Q
Q


(4)

ここで W は HF 法のエネルギー勾配にも登場するエネルギー重み付き密度行列である。(4)式では H-F 力にペア行列と ζ に依存する項があるものの、Pulay 力の表式は HF 法に対するものと全く同じになる。
図 2 にオリゴエン C10H12 分子を HF 法、
155

CASSCF 法、ζ-CHF 法で構造最適化した
時の各 C–C 結合距離をプロットした。結150

CHF 法では結合距離がほぼ一定となってしまう。CHF 法では静的電子相関を過R(C-C) [pm]合交替は RHF 法が最も大きく、続いて
CASSCF 法、UHF 法の順に小さくなり、

RHF
UHF
CHF ( = 1.00)

-CHF ( = 0.90)
-CHF ( = 0.85)
CAS(10,10)

145

140

剰に取り込み、非局在化した構造を好む
ことがわかる。この傾向も、ζ を導入することによって改善し、1 よりも小さいζ を用いることで、UHF 法や CASSCF 法に近い結果が得られることが確かめら
れた。

135

130

0

2

4

6

8

10

PositionFig. 2. Optimized C–C bond-lengths of C10H12 oligoene
obtained by RHF, UHF, CASSCF, and (ζ-)CHF
calculations (cc-pVTZ).

[1] P.R. Surján, in Correlation and Localization (Springer, 1999), pp. 63–88.
[2] M. Tarumi, M. Kobayashi, and H. Nakai, Int. J. Quantum Chem. 113, 239 (2013).
[3] V.N. Staroverov and G.E. Scuseria, J. Chem. Phys. 117, 11107 (2002).
[4] T. Tsuchimochi and G.E. Scuseria, J. Chem. Phys. 131, 121102 (2009).
[5] G.E. Scuseria, C.A. Jiménez-Hoyos, T.M. Henderson, K. Samanta, and J.K. Ellis, J. Chem. Phys. 135,
124108 (2011).

4P095
DyF 分子の基底・励起状態の解析（DIRAC-KRCI による）（中京大国際教養*，名市大名誉教授**）○山本茂義*，舘脇洋**(Chukyo Univ.,* Nagoya City Univ.**) °Shigeyoshi Yamamoto,* Hiroshi Tatewaki**【 緒言 】
ジスプロシウムは磁性で重要なランタノイドであるが，4f 軌道の縮重，相対論効果，電子相関
効 果 の 考慮が 必 要 である ． 二 原子分 子 DyF の 基 底 状態は Kaledin ら の 研 究 [1,2]に よ り ，（は電子の全角運動量の分子軸射影．）Dy+(4f106s1)F–(2p6)
Dy+(4f96s2)F–(2p6)，=7.5 とされている．
配置が隣接するため，Dy+ 部分の 4f, 6s, 6p から生じる多数の電子配置を考慮することが重要である．4 成分相対論に基づく CI 法により，3.0 eV 以下の励起状態を求め，f-shell Omega decomposition
法[3,4]を適用して励起状態の帰属を行う．これは，CI 係数の重みを Ωf（f 電子だけで合成される
角運動量）ごとに分割し，波動関数の特徴付けを行う方法である．電子相関効果は，DIRAC-2012
プログラム[5]の中の KRCI モジュール[6]を使用した RASCI 計算により考慮する．
【 CASCI(11) 】
Dy の基底関数は Dyall の triple-zeta を縮約せず，d 関数（exponent 0.020134537）を 1 個増強した[30s 24p 17d 11f 2g]である．F の基底関数は Koga, Tatewaki, Matsuoka 開発の[13s 9p 2d]．分子全体で大成分が 396 GTFs である．(4f, 6s)を開殻とした average HF/DC (Hartree-Fock /Dirac-Coulomb
Hamiltonian)計算により得られるスピノルを CI 計算に用いた．
まず，(4f, 6s, 6p) の計 11 個の軌道（22 スピノル），11 電子の CASCI 計算を行った．Slater 行
列式の数は 14,094（=7.5 の場合）である．基底状態が =7.5 であることが確認できた．また，
=7.5 の最低 4 状態は，Dy+ 部分の主配置が，それぞれ (4f96s2), (4f106s1), (4f106s1), (4f106s1)となり，実験と符合した．これは電気的中性の分子の HF/DC スピノルを使った場合である．カチオンのス
ピノルを使った場合は 4 番目まで全て (4f96s2) となって実験と合わないため，本研究では中性分子のスピノルを一貫して使用する．
なお，(4f, 6s, 6p, 5d) の計 16 個の軌道を用いた CASCI 計算も行い，3.0 eV 以下では 5d の寄与は無視できることを確認した．
図 1 に示すように多数の電子状態が存在する．3.0 eV 付近に 6p が関与した状態があることが確認できた．これは実験事実と符合する．191 個の電子状態に対して f-shell Omega decomposition 法を適用し，31 個の「family」に分類した．family は最高の状態から派生する一連の電子状態のことで，図 1 においては，線で連結して示されている．これまでに発表した EuF[3]や GdF[4]との違いは，decomposition の表を 4f9 と 4f10 の 2-column 構造とする点である．
遷移モーメントが大きな（>4 Debye）遷移が多数見出された．これらは，低エネルギー (4f106s1)
状態と高エネルギー (4f106p1) 状態との間の遷移であった．実験からは，(4f96s2)の基底状態 =7.5(v=5)と高エネルギー (4f96s16p1) =8.5 (v=0) 状態間の遷移，その他が報告されており，若干の相違がある．今回，フランク・コンドン因子を考慮していないが，検討の余地がある．
【 RASCI(0|8|125) 】

(4f, 6s) 計 8 個の軌道（16 スピノル）を RAS2 とし，スピノルエネルギーが 20 au 以下の 125 個の仮想軌道を RAS3 とした RASCI(0|8|125)を実行し 4 根を解いた．CI の次元数は 8,550,226（=7.5の場合）である．RASCI によるポテンシャルエネルギー曲線（PES）から，(4f96s2) の基底状態に対し Re=1.976 Å, e=596 cm-1 を得た．これらは実験値 1.945 Å, 610 cm-1 [2]に対応している．また，(4f106s1) =8.5 (family c) 状態では，Re=2.086 Å, e=488 cm-1 を得た．Gotkis [7]の推定値 2.057 Å, 499
cm-1 に対応している．T0 の計算値は 0.52 eV であった．なお，実験値は 0.30 eV である．
さらに F-2s, F-2p, 5s, 5p の 8 占有軌道を RAS1 とした RASCI(8|8|125) を，R=1.976 Å では =7.5，
R=2.086 Å では =8.5 に対し実行した．両者のエネルギー差（T0 相当）は 0.63 eV で，RASCI(0|8|125)
より結果は悪くなった．なお，RASCI(8|8|125) の次元数は 483,922,171 である．

Figure 1.

The low-lying electronic states of DyF classified into families.

【 参考文献 】[1] M. C. McCarthy, J. C. Bloch, R. W. Field, and L. A. Kaledin, J. Mol. Spectr. 179, 253 (1996).
[2] L. A. Kaledin, R. T. Holbrook, and J. A. Kunc, J. Appl. Phys. 83, 3499 (1998).
[3] S. Yamamoto, H. Tatewaki, and H. Moriyama, Theor. Chem. Acc. 131,1 (2012).
[4] S. Yamamoto and H. Tatewaki, Comput. Theor. Chem. 980, 37 (2012).
[5] DIRAC, a relativistic ab initio electronic structure program, Release DIRAC12 (2012), written by H. J.
Aa. Jensen, R. Bast, T. Saue, and L. Visscher with other contributors.
[6] S. Knecht, H. J. Aa. Jensen, and T. Fleig, J. Chem. Phys. 128, 014108 (2008).
[7] I. Gotkis, J. Phys. Chem. 95, 6086 (1991).

4P-096

多中心数値積分の精度向上を目指した球面積分法の提案（東大生研） ○松田潤一，平野敏行，佐藤文俊

Proposal of the Spherical Integration Method
toward the Accuracy Improvementof Multi-center Numerical Integration(Institute of Industrial Science, The University of Tokyo)○Junichi Matsuda, Toshiyuki Hirano, Fumitoshi Sato【序】
密度汎関数（DF）法では、交換相関（XC）項の計算に数値積分法（グリッド法）を用いるのが一般的である。その精度はグリッド点数に依存し、振動計算など高い精度が要求される場合は、非常に大きなグリッドセット（GS）が必要となる。XC 項の計算は計算律速の一つのため、大規模 DF 計算においては、適切な GS を設定する必要がある。
我々のこれまでの研究から、分子に対する球面グリッドの配置角により精度のばらつきが発生
し、より大きな GS を用いても精度が向上しない場合があることがわかっている[1]。そこで、空間分割と球面グリッドとの関係に着目し、球面グリッドの配置角を適切に設定することで精度のばらつきを最小化する研究を行なっている[2]。本研究では考案した方法が、ガウス型基底関数を用いた計算で、XC ポテンシャル計算の精度向上において有効であることを確認した。また、XC
エネルギー計算についてもグリッド角の設定方法を新たに考案し、同様に有効であることを確認した。
【球面積分と空間分割】グリッド法において、全空間積分である多中心積分は、各原子を中心とした単中心積分への展開と、動径積分および球面積分グリッドによる単中心積分の組み合わせにより計算される。このときの単中心積分は𝑁
𝑁
𝐼𝐴 = 𝑤𝐴 𝒓 𝐹 𝒓 𝑑𝑟 𝑑𝛺 = 𝑖 𝑟 𝑤𝑖𝑟 𝑗 𝛺 𝑤𝑗Ω 𝑤𝐴 𝒓𝒊,𝒋 𝐹 𝒓𝒊,𝒋 − 𝛿𝐴
(1)であり、空間分割重み𝑤𝐴 によって被積分関数𝐹 𝒓 に窓をかけた形となる。𝛿𝐴 はグリッド法により発生する誤差である。式(1)の総和を全原子についてとれば、積分空間を重ね合わせたこととなり、全空間積分の値が求まる。
本研究では、式(1)のうち球面積分𝑁𝛺Ω(2)
𝑗 𝑤𝑗 𝑤𝐴 𝒓𝒊,𝒋 𝐹 𝒓𝒊,𝒋における誤差のばらつきの解消を狙う。式(2)は𝐹 𝒓 の形状に関わらず、球面グリッドの配置角に対する誤差の分布およびその幅が、ある程度類似していると分かった。よって球面積分グリッドと空間分割との組み合わせの段階で既に無視できない誤差が発生していると考えた。これは同時に、球面グリッド配置角の適切な設定によって、単中心積分への展開時に生じる誤差を球面積分で補正できる可能性を示唆している。
【球面グリッド配置角の決定方法】
XC 項の計算における𝐹(𝒓) は XC ポテンシャル𝜇XCである。その誤差分布とよく一致する既知の𝐹ref (𝒓) が存在するならば、適切な球面積分グリッドの配置角を事前に知ることができる。これを用いれば、球面積分における誤差の打ち消しによって計算精度が向上する。
検討の結果、𝐹ref 𝒓 = 1 が全体の積分誤差への影響が大きい大部分の球面積分計算において有効であると分かった。その一例を Fig.1 に示す。黒線で示される配置角では𝐹 𝒓 = 𝜇XC の球面積
分値が真値となる。紫線の角度では 𝐹ref 𝒓 = 1 の球面積分値が真値となる。両者がほぼ一致していることがわかる。すなわち、SCF 計算の前に𝐹ref 𝒓 = 1 についての計算から適切な角度を見積

もっておくことができ、この角度となるように回転を施したグリッドを用いれば、XC ポテンシャ
ルの計算について、従来の方法よりも 1 桁程度またはそれ以上の精度向上が期待できる。
【本方式の実用化】
分子のワールド座標が DF 計算の結果に影響を与えることが既に知られている。普及している
多くの DF 計算ソフトウェアは、分子の向きを一意に決定することで回転不変性を確保し、同一分子での計算結果が変わらないようにしている。しかし DFT-MD など、分子構造が変化する場合にはこれだけでは不十分となる。本提案方式では十分な回転不変性を実現できる可能性がある。
ただし、本方式はまだ完全ではない。半径 8Å 付近の球面積分の層で、何故か Fig.1 の相関関係が崩れることが分かっていて未解決である。遠方からの寄与であるため誤差は大きくはないが、なまじ他の層での結果が良好なため、現在の本方法の主な誤差となっている。したがって、この原因が究明されれば更に計算精度が上がる可能性がある。その詳細と解決方法については当日報告を行いたい。また、XC エネルギーの積分計算における本方式の有効性についても合わせて報告する。

Fig.1 球面積分値とグリッド配置角の例
x,y 軸: 球面グリッドの角度 φ,θ 分子,原子: Gly,Cα 被積分値:μXC
半径: 2.7Å NΩ=194 黒線: 誤差 0 紫線: 被積分値 1 で誤差 0【参考文献】
1. Junichi Matsuda, Toshiyuki Hirano, Fumitoshi Sato, ”A Study of the Numerical Integral
Scheme for Large-scale Density Functional Calculation”, 14th International Density
Functional Theory Conference (DFT11), P65, 2011.
2. 松田潤一, 平野敏行, 佐藤文俊, “球面グリッドの回転による多中心数値積分の精度向上”, 第6 回分子科学討論会, 4P-099, 2012.

4P098複素基底関数法における複素軌道指数最適化(慶大院) ○松崎 黎, 藪下 聡
Optimization of complex orbital exponents in the complex basis function method(Keio University) ○Rei Matsuzaki, Satoshi Yabushita【序】光電子分光は広範な実験で用いられているが、その理論的解析はいまだに不十分な場合が多い。我々は複素基底関数法[1]を用いて光イオン化の計算に取り組んでいる。複素基底関数法のデメリットとして、微分断面積の計算に非常にたくさんの基底関数が必要であることが指摘されてきた[2]。この欠点を克服するために、複素数
軌道指数まで変分的に最適化することで最良な基底関数セットを構築する手法を開発してきた。その過程で最適な軌道指数を複素平面状にプロットすると特徴的なパターンが現れたので本研究ではその原因を探った。
【理論】光イオン化の断面積は、双極子近似の元で、式(1)の解である一次摂動波動関数
できる。ここで、エネルギー の光を固有エネルギーを持つ束縛状態から計算することがに照射したときを考えている。式(1)を汎関数形式に変形すると、式(2)のようになり、汎関数形式を用いると変分原理による計算が可能になる。ここで、丸括弧はブラ側の複素共役をとらない積分を意味する。複素基底関数法では 複素関数の基底関数セットの線形結合を試行関数として用いる。本研究では式(3)で表されるような複素 STO を用いる。線形結合をとった時の展開係数を変分原理によって最適化すると、汎関数は式(4)のように、基底関数に含まれる軌道指数に依存する行
列形式で表現される。それぞれの行列およびベクトルは式(5)(6)(7)で表現される。
【計算】本研究では基底関数に含まれる軌道指数 を最適化することで、最適な基底関数セットを構築する。式(4)は試行関数の展開係数のみを最適化した結果なので、 を最適化するには、さらに式(4)に対してなる点を数値的に探せばよい。我々は解析的微分法を用いた Newton-Raphson 法により の最適化を行った。
【計算結果】図 1 に基底状態の水素原子の光イオンに対して 5 つの複素 2p-STO を用いて最適化した複素軌道指数の複素平面状での分布を示す。プロットはいくつかの光電子の波数 に対して行った。図を見ると、実軸付近のは始状態 1s の軌道指数 に集積し、虚軸付近では光電子の波数に集積する。また、全体としては円弧の振る舞いをすることがわかる。この振る舞いの原因として基底関数同士の重なりが等しくなるように振る舞うことが原因の一つとして考えられる。つまり STO の場合、軌道指数が等比数列になると基底関数同士の重なりが等しくなる
が、今回の計算では軌道指数は複素数なので、その値が複素等比数列に近い形に最適化され複素平面上で円弧を描いていると考えられる。さらに詳しく調べるために、
図4にの解析関数としての性質に注目した。図 2～の虚部を一つの引数だけを変数として他の引数を最適点に固定して作った図を載せる。解析関

数の性質より、実部および虚部はラプラス方程式を満たすので、一般に関数の安定点は鞍点になるが、今回の計算結果にもそのことが反映されている。また、図 2～4 の凡例をみると、変数を動かした時の関数の敏感さに大きな差があり、実軸に近い軌道指数は敏感であり、虚軸に近くなると鈍感になっていく。I の安定点付近の振る舞いを調べるために、次のような手順で Hessian の固有値解析を行った。
よび、を成分とする 10×10 の行列の固有値問題を解いたところ、その固有値が

とし、

お
の形に

得られた。以上の特徴は(4)式の汎関数が、複素軌道指数を変数とする多変数解析関数であることを反映したものである。ポスター発表においてさらに詳しい解析を説明する予定である。

(1)

図1. 5つの基底を用いた場合の最適化した

(2)

(5)

(3)

(6)

(4)

(7)

図2. を変数,他の を最適値に固定したときの

図3.

を変数としたとき(左図)とを変数とみなしたとき(右図)の。他の は最適値に固定

図 3.

を変数としたとき(左図)とを変数とみなしたとき(右図)の。他の は最適値に固定【参考文献】[1] T.N.Rescigno et al., Phys. Rev. A, 31, 624(1985)[2]C.W.McCurdy et al., ibid., 35, 657(1987)

4P099量子化学計算によるスクアレン分子の無輻射遷移過程の研究（慶大院・理工）○小澤仁嗣、八代和徳、山本琢磨、藪下聡

Theoretical study of the nonradiative transition process of the squaraine molecule(Keio Univ. Science and Technology) ○Hitoshi Ozawa, Kazunori Yashiro, Takuma Yamamoto, Satoshi Yabushita

・序Bis[4-(dimethylamino)phenyl]squaraine(以後スクアレン分子と略す)は、長波長領域に蛍光を

図１．スクアレン分子の構造式[1],[2]発する色素である。この分子については、複数の実験研究者により異なった形状の蛍光スペクト
ル、及び無輻射遷移の顕著な温度依存性などが観測されている[1],[2]。そしてその励起緩和機構は、実験研究者による解明が試みられたにも関わらず、未だに詳細はわかっていない。
我々は、量子化学計算によるこの分子の励起緩和機構を解明することを試みてきた。これまで、PCM 法を用いた第１励起状態の平面→ねじれ構造の反応熱の溶媒極性依存性の計算から、第１励
起状態ではこの分子の内部回転は起こらないことを示した。また、CASSCF 法、MRMP2 法による平衡構造、及び Duschinsky 効果を含めた Franck-Condon 因子を用いた吸収・蛍光スペクトルの理論計算から、この分子の蛍光スペクトルの起源が B3u 励起状態の平面構造からのものである
ことを示した。本研究では、量子化学計算によるこの分子の無輻射遷移過程の解明を試みる。特にこの分子の無輻射遷移過程の可能性として、B3u 励起状態とエネルギー的に近接する B3g 励起状態を経由するもの、及び B3u 励起状態からの Jortner の理論あるいは Marcus 理論的な電子移動反応の２通りを考え [3],[4] 、それぞれの検討を行う。その結果、この分子の無輻射遷移過程は、
Marcus 理論的な分子内電子移動反応である可能性が高いことを示す。
・B3g 励起状態の検討
ま ず 、 B3g 励 起 状 態 を 量 子 化 学 計 算 で 検 討 し た 。 こ の 励 起 状 態 は 、 真 空 中 の 計 算（LC-TDBOP/6-31G(D)）では第１励起状態であるが、溶液中（PCM-LC-TDBOP/6-31G(D)）では第２励起状態となり、その第１励起状態である B3u 状態とのエネルギー差は約 0.4～0.5eV である。また、B3u 励起状態の方が B3g 励起状態よりも溶媒効果を強く受けるため、溶液中では B3u 励
起状態から B3g 励起状態経由の無輻射遷移は、活性化エネルギーが高く起こらないと考えられる。
・Marcus 電子移動理論を応用した B3u 励起状態からの無輻射遷移過程
次に、B3u 励起状態からの電子移動反応が起こる可能性を考えた。B3u 励起状態では、基底状態からの電子励起の際に、置換基部分から中央４員環へ-0.08e～-0.2e の電子移動が起こり、さらに再配向エネルギー（溶質と溶媒の両方を含む）の溶媒極性依存性が大きい。
まず電子移動反応が、Jortner の理論に従って溶媒とスクアレン分子の両方の振動モードの寄与によって起こる場合を考え、数値計算による検討を行ったところ、スクアレン分子の振動モードが寄与した電子移動反応は起こらないことが分かった。この原因は、この溶質分子の基底状態と励起状態の構造のずれは非常に小さく、その Huang-Rhys 因子の総和は小さいため、相対的に溶媒の振動の寄与の方が大きいためと、低振動数の溶質の振動モードの Huang-Rhys 因子の相対値が大きいことから、中間温度域の振舞い（熱エネルギーk B 𝑇よりも溶質の振動エネルギーの方が大幅に高い場合）が成り立たず、高温極限（Marcus 理論）の振舞いをするためと考えられる。
Marcus 理論によれば、反応自由エネルギーΔG の溶媒極性依存性は、電子移動反応の活性化
エネルギーに、つまりその速度に影響を与える。そこで、この分子の無輻射遷移を、Marcus 理論における逆転領域での分子内電子移動反応とみなすことができるとした。このような解釈の妥
当性は、過去にカルバゾール誘導体で報告されている[5]。このとき、その無輻射遷移速度は次式で与えられる。
(1)𝑘nr = (2π|𝑉DA |2 ⁄ℏ)(4π𝜆s k B 𝑇)−1/2 exp(− 𝐸a ⁄k B 𝑇)

ここで、𝜆sは再配向エネルギー（溶質と溶媒）、𝐸a は活性化エネルギーで、また𝑉DA は電子トンネ
ル行列要素である。初めに、PCM-LC-TDBOP/6-31G(D)で計算された再配向エネルギー（溶質と溶媒）（ジクロロメタンとエタノールで 0.205eV、0.332eV）を用いて活性化エネルギーを Marcusの式 Ea = (G+B)2/(4B)（調和近似、G は反応自由エネルギーであり、ほぼ電子励起エネルギー）で評価したところ、実験値をかなり過大評価した。これは、この無輻射遷移過程が弱い電荷
分離状態間の電子移動反応であるため、近接溶媒分子のみが反応に大きく寄与し、その結果、誘電飽和の影響が活性化エネルギーに強く生じるためと考えられる[6],[7]。Rettig らのエタノール中
アレニウスプロットの実験データと再配向エネルギーの計算値を用いて𝑉DA の値を逆算したところ、約 1kJ/mol であった。そこで、この分子系の𝑉DA が 1kJ/mol で、エタノール中での𝐸a の実験
値が 0.165eV～0.172eV であることを踏まえ、𝐸a が 0.17eV と 0.24eV の場合の電子移動反応速度を(1)式を用いて数値計算して図２を得た。

8.00E+09

knr [s-1]

6.00E+09
4.00E+09
Ea=0.17eV

2.00E+09

EtOH
CH2Cl2Ea=0.24eV

0.00E+00
150

200

T [K]

250

300

図２． 電子移動反応速度の計算結果（破線が実験値、実線が計算値）
その結果、電子移動反応速度の計算値が、無輻射遷移速度の実験値と同程度の大きさになることがわかった。このことから、この分子の無輻射遷移過程は Marcus 理論的な電子移動反応であり、0.17eV、0.24eV の活性化エネルギーはそれぞれエタノール、及びジクロロメタン中のものに対応する可能性があることがわかった。
[1] Law, K. Y. J. Phys. Chem. 1987, 91, 5184.
[2] Cornelissen-Gude, C.; Rettig, W.; Lapouyade, R. J. Phys. Chem. A 1997, 101, 9673.
[3] Jortner, J. J. Chem. Phys. 1976, 64, 4860.
[4] Marcus, R. A. Rev. Mod. Phys. 1993, 65, 599.
[5] Kapturkiewicz, A.; Herbich, J.; Karpiuk, J.; Nowacki. J. J. Phys. Chem. A 1997, 101, 2332.
[6] Yoshimori, A. and Kakitani, T. J. Phys. Soc. Jpn 1992, 61, 2577.
[7] Yelle, R. B. and Ichiye, T. J. Phys. Chem. B 1997, 101, 4127.

4P100
局所制御理論による非極性直線分子の配列制御レーザー場設計（慶大院・理工）○菅原道彦

Theoretical designing of alignment control laser field for non-polar linear
molecules by local control theory(Keio Univ.) ○M. Sugawara【序】レーザーを利用した分子の振動・配向状態の量子制御は、理論・実験の両面から注目を集めている。レーザー・分子間の相互作用として双極子近似が有効である分子系
では、局所制御理論による制御場設計が可能である。一方で、非極性分子の配列制御等、レーザー場の２乗に比例する分極相互作用が重要である系では、動的シュタルクシフトを考慮する必要があり、通常の局所制御理論に基づいた定式化で扱うことが難しい。本
研究では、制御レーザー場を主要振動数で振動する部分とそれに比べて十分にゆっくり変化するエンベロープ関数の積として表現することにより、局所制御理論によるレーザ
ー場設計が可能であることを示す。
【理論】非極性分子のハミルトニアン Ĥ 0 、分極演算子を α̂ 、制御入力レーザー場を u (t )としたとき、全系のシュレディンガー方程式di
Ψ (t ) = − ( Hˆ 0 − αˆ u (t ) 2 ) Ψ (t ) ,(1)hdtと与えられる。但し、 Hˆ 0 = ∑ Ei i i （分子の回転固有状態： { i } 、固有エネルギー：{Ei } ）。ここで、 u (t ) を次の形に規定する。

u (t ) = ∑ 2uij( e ) (t ) cos[ωij t ] .

(2)

i≠ j

ωij は状態 i と j の間の２光子遷移に関わるレーザー場の角振動数成分である。また、
2uij( e ) (t ) は各成分のエンベロープ関数であり、その時間変化は ωij の時定数と比較して十分ゆっくりであるとする。(1)式中のハミルトニアンはHˆ 0 − αˆ u (t ) 2 = Hˆ 0 + αˆ ∑ uij( e ) (t ) + αˆ ∑ uij( e ) (t ) cos[2ωij t ]
i≠ j

(3)

i≠ j

と変形され、右辺第２項はエンベロープ関数の時間変化に依存する動的シュタルクシフ
トを生じさせる項である。ここで、エンベロープ変化について断熱的な描像を導入し、
各 時 刻 で Gˆ 0 ≡ Hˆ 0 + αˆ

∑ u (t ) を 対 角 化 す る ユ ニ タ リ 変 換 T% を 導 入 す る 。 演 算 子(e)ij% (t ) ≡ T% † Ψ (t ) （ 但 し 、
F%0 = T% †Gˆ 0T% 、 β% = T% †αˆT% 、 及 び 新 た な 系 の 状 態 Φ
F% φ% = λ% φ% (i = 0,1,L) ）を導入するとシュレディンガー方程式(1)は、
0

i

i

i

d %
i
% (t )
Φ (t ) = − ( F%0 − β% ∑ u%ij( e ) cos[2ωij t ]) Φ
(4)
dt
h
i≠ j
と変形できる。さらに、時間発展演算子 U% 0 (t ) ≡ exp[−iF%0t / h ] を用いて相互作用表示に移行する。すなわち、 F (t ) = U% † (t ) F% U% (t ) 、 β (t ) = U% † (t ) β%U% (t ) 、 Φ (t ) ≡ U% † (t ) Φ (t ) を0

導入すると、(1)式は

0

0

0

0

0

0

d
i
Φ (t ) = − β (t )∑ u%ij( e ) cos[2ωij t ] Φ(t )
dt
h
i≠ j

(5)

(

)

ここで、動的シュタルクシフトを考慮した共鳴条件である 2ωij = λ% j − λ%i / h を課すことによって激しく変化する振動数成分 cos[2ωij t ] を取り除くことができ、最終的に(5)式はdi
Φ (t ) = − ∑ βij u%ij( e ) Φ (t )dth i≠ j(6)へ変形される。ここで、 β ij は i → j 間の２光子遷移に関わる分極定数成分である。この運動方程式に局所制御理論を適用し、エンベロープ関数 u%ij( e ) 及び断熱的に変化する各振動数成分を設計する。制御の基準となる評価関数を制御量演算子 Ô の期待値として
y (t ) = Φ(t ) Oˆ Φ(t ) のように定義する。 y (t ) の時間微分を(6)式を用いて計算すると、

y& (t ) = 2∑ Im Φ (t ) Oˆ (t ) βij Φ(t ) u%ij( e )

(7)

i≠ j

となる。この時間微分を常に 0 以上に保つようにするためには、各パルスエンベロープをu% ( e ) (t ) = A Im Φ (t ) Oˆ (t ) β Φ (t ) ( A > 0 ) と設計すればよい。
ij

0

ij

0

【結果】等核２原子分子モデル系の回転準位の分布制御に上述の理論を適用した。基底として、球面調和関数の基底状態 Y0,0 = 0 と第２励起状態 Y2,0 = 1 をとり、初期状態の 0 から励起状態 1 への遷移を２光子共鳴過程を利用して実現する。Oˆ = 1 1 と定義し、y (t )として状態 1 の準位分布をとり、上述の理論により設計されたレーザー場の ω12 及び制御レーザー場 u (t ) を図１に示す。図２に示すように 0 、 1 の準位分布の変化を見ると、設計指針通り状態 1 の準位分布が短調増加している。一方、設計された制御レーザー場
照射下の、状態 0 及び 1 の準位分布の時間変化を図２に示す。最終状態として、分子がより配列している状態である目標状態 1 は得られているものの途中では 1 の準位分布の単調増加は保証されていない。これは、設計理論で保証されている条件があくまで 1の分布の単調増加であり、レーザー場によω12る動的シュタルクシフトが顕在化する時
間領域では 1 と 1 が異なった状態であu (t )

ることに起因する。しかし、エンベロープ
関数が断熱的に変化する限り、パルス照射
前後では、状態 0 と 0 、状態 1 とがそれぞれ一致するため、(6)、(7)式を基にした制御則を適用することが出来る。

0

図２

1

状態 0 と 1 の準位分布の時間変化

図１

0

図３

設計された制御レーザー場（赤・下）と、その中心振動数（青・上）の時間変化

1

状態 0 と 1 の準位分布の時間変化