3P-101
CoH の ab initio 計算における DK 相対論補正の奇妙な振舞（産総研 1，NEC 2，お茶大理 3） 長嶋 雲兵 1，友成六美 2, 平野 恒夫 1, 3
Strange Behavior of DK Relativistic Correction on the CoH ab initio Calculation(AIST1, NEC2, Ochanomizu U.3) Umpei Nagashima1, Mutsumi Tomonari2, Tsuneo Hirano1,3
CoH の電子構造に関しては、Slater 型基底関数を用いて非相対論の範囲で MR-SDCI+Q および
MRCPA レベルの ab initio 計算を行い、実験値と良く対応する分子構造や分子定数を得て、その詳細を報告した 1)。既に、我々の研究を含めて、一般に、周期律表で Mn より右の第 1 列遷移金属
原子を含む分子について「分光学精度」の計算を行う場合には相対論補正をする必要があることが分かっている。 CoH も例外ではないと考え 3 次の Douglas-Kroll (DK3) 相対論補正を行う abinitio 計算を行ったところ、極めて奇妙な結果を得たので、その理由を考察すると共に、DK3 の下での最も信頼出来る計算方法を模索したので報告する。計算は Gaussian 型基底関数（GTF）を用い、C2v 対称の分子として、MOLPRO 2010.1 で行った。
我々の経験では、第 1 列遷移金属を含む分子の計算では、金属の 3d, 4s, 4p を active space とする no-core,full-valence の CASSCF を 行 っ て 、 そ の 軌 道 を 使 っ た
MR-SDCI＋Q の計算を行い相対論補正を加えれば、ほぼ構造、振動数、双極子能率の実験値を再現出来ることが分かっている 2)。 CoH に関して、3 種類の DK3 で最適化された GTF を使って DK3 ハミルトニアンのもとで計算した結果が表１の(1)、(2)、(3)である。実験値と較べて明らかに結合距離が短かすぎるし、従って伸縮振動数も高すぎる値になっている。つまり、DK3 法のovershooting なのである。
一般に、相対論補正を施すと核上に分布を持つ遷移金属原子の s 軌道は収縮し、他の軌道は広がる。
3

Φ CoH の電子配置は図１に示すように、C2v 対称のもとでは、B1 と B2 の 2 重縮退した状態として求Table 1.

Equilibrium structure and vibration frequency of 3Φ CoH at various levels of calculation.
Active space

Basis setsa)

CASSCF

MR-SDCI

DK or
NR

0) Exp. (cf., Tomonari et al (2007)1))r0 /Åωe /cm-1

ν / cm-1

1.516 1.526

1925

1856

re /Å

1) 5ZP ANO-RCC

3d,4s,4p3d,4s,4p

DK3

1.507 1.523

2002

1965

2) aV5Z-DK

3d,4s,4p3d,4s,4p

DK3

1.507 1.523

2004

1973

3) Sap-aQZP-DK

3d,4s,4p3d,4s,4p

DK31.506 1.533

2010

1987

3) Sap-aQZP-DK

3d,4s,4p3d,4s,4p

DK31.506 1.533

2010

1987

4) Sap-aQZP-NR

3d,4s,4p3d,4s,4p

NR

1.517 1.531

1923

1879

5) Sap-aQZP-DK

3s,3p,3d,4s,4p

3d,4s,4p

DK3

1.516 1.522

1919

1845

6) Sap-aQZP-NR

3s,3p,3d,4s,4p

3d,4s,4p

NR

1.509 1.519

1970

1843

a)

Sap: Sapporo GTF

まり、それぞれが b1 型の縮退した 2 個の配置関数の和と差として表される。Co と H を結びつけるのは、
図１に示すように、Co の 4s-3dσの混成軌道と H 1s との間の 2 個の結合性軌道であるため、相対論効果を取り入れた計算では、4s 軌道の収縮が顕著になって結局 Co-H 結合が短くなる。その上、3d 軌道が膨らんで 4s 電子を充分に遮蔽しないことも考えられる。 Hydride 以外の第 1 列
遷移金属を含む分子でも相対論効果を取り入れた計算では金属−配位子間の結合距離は短くなる
が、hydride の場合のように金属の 4s 軌道が配位子の 1s 軌道と直接相互作用して結合が出来ていない場合は、相対論効果による結合距離の短小化はせいぜい 0.001 Å 程度で、余り顕著ではない。
基底関数にも注意が必要である。一般に、相対論補正を行うためには、相対論のもとで最適化された基底関数を使用することが必要であり、非相対論計算を行うときは非相対論のもとで最適化された基底関数を使用することが必要である 3)。Sapporo-GTF 基底関数ライブラリ 4) には、両者が用意されているので、(3)と同じやりかたで、非相対論レベルの計算を行ったのが、表１の(4)である。奇妙なことに、相対論効果を無視した(4)の結果は結合距離、振動数の両者とも、実験値を極めて良く再現している。しかし、より低い近似の計算結果の方が、より正しく現実を再現しているということは、ab initio の原理に反している。つまり、偶発的な一致なのであろう。 CoCNの基底状態 3Φの計算 5) には、相対論効果の補正が必須だったので、(4)の結果を鵜呑みにしてはならない。
表１の(5)は、DK3 の近似のもとで、Co の内核の 3s、3p を CASSCF の active space に入れて、
3s と 4s 軌道を意図的に拡げた結果である。今度は ab initio の原理通り、結合距離も振動数も実験値を良く再現している。この方法を非相対論的に行ったのが(6)であるが、結合距離、振動数は実
験値とは乖離した結果になった。上記のように、非相対論計算では、full-valence の計算で充分なのである。
(5)の方法で、結合距離と振動数に関して実験値を再現出来たが、もう一つの重要な要素である
電荷分布を反映している双極子能率の値はどうなっているのであろうか。よい ab initio 計算とは
この 3 者がともに正しい実験値を再現していなければならない。MR-SDCI 計算の波動関数の期待値としても求まる双極子能率は MR-SDCI 計算での CI 展開における truncation のため、必ずしも正しい双極子能率の値を示さない。エネルギー的に効かなくても、電荷分布に寄与する配置関数が切り捨てられている可能性があるからである。実験は Stark 効果で求めるので、それに対応したエネルギーの電場微分としての双極子能率を求めて見た。
方法（3)、(4)、(5)、(6)の電場微分の双極子能率は、それぞれ、2.60 (2.68)、2.82 (3.00)、1.17 (2.74)、
13.8 (3.40) D となり（括弧内には期待値を示す）、いずれも実験値 6)の 1.88(8) D とは異なった値になっている（特に方法(6)はひどい結果になっている）。電荷分布には、荷電子の動的電子相関が効く筈なので、現時点で最も良い方法 である (5)で、MR-SDCI の active space に Co の内核である 3s、3p を加えた計算を試みている。これで、CASSCF と MR-SDCI の両者でバランスの取れた計算になる。膨大な CPU 負荷なので本要旨には間に合わないが、その結果に期待したい。
M. Tomonari, R. Okuda, U. Nagashima, K. Tanaka, T. Hirano, J. Chem. Phys., 126, 144307 (2007).
T. Hirano, R. Okuda, U. Nagashima, K. Tanaka, P. Jensen, J. Molecular Spectrosc.,250, 33 (2008).
3
平野 恒夫, 長嶋 雲兵，Per Jensen, 第 5 回分子科学討論会, 1E20 (2011).
4
http://setani.sci.hokudai.ac.jp/sapporo/Welcome.do5
T. Hirano, R. Okuda, U. Nagashima, P. Jensen, Mol. Phys., 105, 599 (2007).
6
H. Wang, X. Zhuang, T.C. Steimle, J. Chem. Phys., 131, 114315 (2009).
１

2

3P102

焼き戻しアンブレラサンプリング法による
水溶媒中のマロンアルデヒドに関する理論的研究（1 分子研，2 名大院理・物理，3 名大院理・構造生物研，4 名大院工・計算科学研）
森義治 1，岡本祐幸 2,3,4

Theoretical study on malonaldehyde in aqueous solution by simulated
tempering umbrella sampling(1Institute for Molecular Science, 2Dept. of Phys., Nagoya Univ., 3Structural Biology
Research Center, Nagoya Univ., 4Center for Computational Science, Nagoya Univ.)
Yoshiharu MORI1 and Yuko OKAMOTO2,3,4【序】
分子の性質を理解すること，特に化学結合の解離や生成の反応機構を理解することは重要である。このような反応機構を分子シミュレーションにより理解したい。しかしながら量子力学を用いた第一原理計算を含む分子シミュレーションは非常に高いコストを要する。つまり長時間のシミュレーションを実行することが難しく，反応を理解するために十分な統計量をとることが困難である。このようなことから，短いシミュレーション時間で効率よくサンプルをとる方法を用いる必要がある。
	  そのような方法として，われわれは拡張アンサンブル法とよばれる一連のシミュレ
ーション手法を開発してきた。これらの方法においては通常のシミュレーションよりも短い時間で多くの構造空間を探索することを可能とする。
	
  本研究では特定の化学反応を起こしやすくする拡張アンサンブル法，焼き戻しアン
ブレラサンプリング法，を提案し，その方法を実際の分子系のシミュレーションに適用した。この拡張アンサンブル法においては，設定した反応座標に関して複数のアン
ブレラポテンシャル間のランダムウォークを実現する。このことから，反応経路の安定状態だけではなく遷移状態も十分にサンプリングできることが期待される。
	 【方法】
反応座標に関する拡張アンサンブルシミュレーションをテストするために，水中のマ
ロンアルデヒドの分子内プロトン移動反応に関する第一原理シミュレーションを実行した。以上の計算を実行するために，計算プログラムとして CP2K (ver. 2.1) を用いた。計算の系としてマロンアルデヒドと水分子 71 個を用いた（総原子数：231）。

第一原理分子動力学シミュレーションは KS-DFT における Born-Oppenheimer ダ
イナミクスにより行った。交換相関汎関数として BLYP を用いた。系の温度は 300
K とし，温度制御法として Nosé-Hoover 法を用いた。初期条件が異なる独立な 4 つのシミュレーションを実行し，それぞれのシミュレーション時間として 50 ps を設定した（総シミュレーション時間：200 ps）。また拡張アンサンブル法で用いる反応
座標として，プロトン移動反応に関与する水素原子と酸素原子間の距離の差を用いた。
その反応座標に関するアンブレラポテンシャルを 11 個用意し，それらのポテンシャ
ル間のランダムウォークを実現させるようにした。
【結果】
本計算はアンブレラポテンシャルに関する焼き戻し法を用いた分子動力学シミュレーションであるので、シミュレーションの間はそのポテンシャル間をランダムウォー
クする。その結果として，水素原子は一方の酸素原子から別の酸素原子に移動する反応を繰り返すことになり，反応座標も対応した時間発展をするはずである。実際のシミュレーションにおいて，反応座標はランダムウォークをしており、化学反応経路に沿った構造を十分にサンプルできた。
	  以上のように反応経路に沿ったサンプルを十分に得ることができたので，その結果から反応座標に関する自由エネルギーを計算した。つまり反応座標に関して平均力ポ
テンシャル（PMF: Potential of mean force）を計算した。そして本計算手法により反応経路間における自由エネルギー地形を正確に計算できることを示した。
	
 以上のように，化学反応において重要な反応座標における拡張アンサンブルシミュレーションを実行し，エネルギー障壁が高い状態に関しても十分なサンプリングができることがわかった。さらにそれらの結果から自由エネルギー地形を計算し，化学反応における自由エネルギーを正確に得ることが可能となったこともわかった。発表においては，化学反応中における水溶媒の構造も議論する予定である。

3P-103
フェニルボロン酸アゾプローブグルコースセンサの超分子機能の理論的研究(上智大学 1 , 京都大学福井謙一記念研究センター2 )島田嶺太 1 , 石田俊正 2 , 早下隆士 1 , 南部伸孝 1

Theoretical study of Suparamolecular function of phenyl boronic acids azo probesensors for glucose(Sophia University1 , Fukui Institute for Fundamental Chemistry, Kyoto University 2)
R. Shimada1 , K.Ishida2 , T.Hayasita1 , S.Nanbu1【序】糖は代謝機能や, 栄養摂取, 細胞構造などの生物の基本的な機能と関わる分子であり,個人の表現型や分化, 免疫能の制御において主要な役割を果たす. このような重要性から単純な糖認識法が考案されてきた. 酵素法は高い選択性を有するが, 酵
素の安定性が十分でなく, 耐久性や連続的観測, 生体内分析などの面で問題がある.
そこで現在, 安定性の高い化学センサ分子の合成に注目が集まっている. 今までに,アゾ色素を骨格とするアゾプローブが多く合成されてきており, これらはアルカリ金属やアルカリ土類金属の認識を行う化学センサとして優れた能力を有している. 早下らは糖認識においてアゾプローブを利用し, これにフェニルボロン酸基を導入することにより, 糖選択性を持つプローブの合成に成功した.[1]これらの分子の一つである
p-BA-Azo [2]を図 1 に示す. 一方, 図 2 にこのプローブのボロン酸基の水酸基部分と糖がと結合したプローブ/糖複合体を示す. この複合体は強い発光スペクトルを示すことが確認されている. しかしながら現段階において,詳細な発光過程が明らかになっていない．

Phenyl boronic acids
図 1. p -BA-Azo

glucose
図 2. p -BA-Azo/glucose 複合体【 理 論 計 算 】 本 研 究 で は Molpro2010.1 を 用 い た 量 子 化 学 計 算 に よ り , p-BA-Azo/glucose 複合体の励起状態の解明を行った．計算方法は CAS(完全活性空間多配置)SCF 法を用い，基底関数には藤永らの midi を使用した．そして，電子基底状態(S0)及び励起状態(S1 及び S2)の最安定構造の決定を行った．さらに遷移双極子モーメントを求め，電子遷移確率を評価した．

【結果と考察】電子状態間の遷移確率を表す遷移双極子モーメントの絶対値を比較したところ，基底状態の最安定状態(S0 min) からの垂直遷移において S2-S0 間の遷移双極子モーメントが 3.08 .a.u と大きな値を示す(表 1)．このことから，この分子は光吸収後に S0 min から S2 状態に遷移すると考えられる．さらに，励起された分子は S2状態上で失活し，恐らく S2 min に到達する．そして図 3 に示す通り S2min のエネル
ギーと S1 のエネルギーがほぼ縮重していることから，非断熱遷移を経て，S1 状態へ移ると思われる．その後 S1 min まで失活し，そこから非断熱遷移により S0 に遷移することが予想される．従って，電子状態のみの計算からは実験事実に反したものとなった．つまり，実際の実験結果における S1からの強い発光を説明することができなか
遷移双極子モーメント / a.u.
表1
S 0 min
S 1 minった．そこで我々はさらに詳細な解析を行うために非断熱 ab initio MD [3]計算を実行し，その結果をポスターにて発表する.

S 1-S 00.01380.0604S 2-S 03.08000.04085Potential Energy / eVS2 4.32
S1

3

S2 min
3.68
3.63

3.48

3.36

S1 min
2.00
1.99

0.98

1
S0 min

S0
-1

0

0
2

4

6

8

10

12

図 3．S0, S1, S2 の最安定エネルギーと各状態における垂直遷移エネルギー．

S0 の最安定構造(S0 min)を 0 とした．

[参考文献][1] Kumai Mio; Kozuka Satoko; Samizo Mariko; et al. , ANALYTICAL SCIENCES,
2012, 28, 121-126[2] 4-{4-[4-(dihydroxy-Boranyl)-phenylcarbamoyl]-phenylazo}-phenol anion
を p -BA-Azo とした.
[3] T. Murakami, M. Nakazono, A. Kondorskiy, T. Ishida, S. Nanbu, Phys. Chem,
2012, 14, 11546 – 11555.

3P-104

環状テトラシロキサン化合物の水溶媒中での構造と励起エネルギーの密度汎関数法計算による研究(群馬大学) 篠原大明、工藤貴子

A DFT Study for the Structure and Excitation Energyof Cyclotetrasiloxane Compounds in the Water Solvent(Gunma University) Hiroaki Shinohara, Takako Kudo
【序】
環状・かご状シロキサン化合物は一般には水に不溶なため、水溶媒中での挙動は殆ど知られていない。そこで本研究では、計算化学の手法により環状シロキサンの一種であるシクロ
テトラシロキサン、(X2SiO)4 (D4)、 の水溶媒中での構造・光化学的挙動を明らかにする。まず X が H のシクロテトラシロキサン(図 1)に注目して、水分子との水和構造及びそれらの構造における最低励起一重項エネルギーについて DFT 計算を行った。さらに X を炭素やケイ素の、飽和および不飽和化合物などに置換した D4 (図 2)についても同様の計算を行い、置換基の変化が D4 の水溶媒中での構造・光化学的性質にどのような影響を及ぼすか調査した。

図 1 H 置換シクロテトラシロキサン

図 2 不飽和シリル基で置換した D4【計算方法】
D4 および溶媒中の D4 の構造最適化とそれらの基準振動解析は基本的に B3LYP/DZP+diffuse 関数レベルで行った。水溶液中での挙動を記述するため、対象分子を囲む広範囲の溶媒効果を PCM 法で考慮した。一方、D4 と相互作用する極近傍の水分子は、実際に水分子を作用させた。尚、比較のため PCM 法では水以外の幾つかの種類の溶媒でも計算を行った。
電子励起エネルギーは TD-DFT 法により求めた。この TD-DFT 計算では、CAM-B3LYP を基本として、長距離補正に優れたものや分散力の補正を考慮した汎関数を用いた。プログラ

ムは GAMESS および Gaussian09 を使用した。
【結果と考察】まず、B3LYP 汎関数を用いて H 置換 D4 の周囲に少数の水分子(1~20 個程度)を作用させ、その最適化構造と最低励起一重項エネルギーを求めた。H 置換 D4 には水分子との水素結合が可能な場所として、酸素原子とケイ素上の水素原子があるが、後者には水分子は配位しなかった。このため、最近傍の水和構造として４箇所の酸素原子に水分子が配位している構造が基本となる。以後、水分子同士が水素
結合ネットワークを形成するよう水分子の数を増加させていった。(図３)

それらの平衡構

造では、水分子同士の水素結合長はほぼ通常

図 3 水配位シクロテトラシロキサン

通りだが、D4 と水分子間の水素結合長は通常のものよりやや長めだった。また水分子と作用することで D4 の主にシロキサン結合部分に変化が見られた。最低励起は HOMO-LUMO 励
起が中心である。表 1 に示す様に、最低励起エネルギーは、水分子を配位させた場合真空中の D4 よりもやや短波長シフトし、その後水分子を増加させていってもわずかな増減があるのみで大きな変化はなかった。
H 体の D4 の H を種々の置換基に変えるとその種類に応じて最低励起エネルギーに変化が見られた。炭化水素やシリル基で置換した場合は励起の主体が置換基に移る傾向があるが、置換箇所によっては D4 のシロキサン骨格の軌道が置換基同士の軌道と相互作用し最低励起に拘わる。また、H 置換 D4 の 8 つの H のうち、例えば 2 つを置換するコンフォメーション
は幾つか考えられ、一つのケイ素から上下に置換基が出ているものは励起エネルギーに大き
図 3 水配位シクロテトラシロキサン
な変化をもたらしたが、他のコンフォメーションではあまり変化が無かった。また置換基を不飽和シリル基にしたものは溶媒中で H 体や炭化水素体とは異なる性質を示した。
表 1 水分子増加に伴う D4 の最低励起一重項エネルギーの推移(nm)
0個

1個

2個

4個

8個

12 個

19 個

27 個

B3LYP

176.6

174.4

174.2

169.9

171.8

177.6

175.4

179.5

CAMB3LYP

165.3

163.7

163.5

159.6

160.1

158.4

158.8

162.1

M06-2X

166.0

164.5

164.0

160.7

161.6

160.1

160.8

163.4

LC-ωPBE

154.2

153.9

153.8

151.0

151.1

149.7

149.4

152.4

3P105GPGPU によるエネルギー分布関数計算の高速化(金沢大院・自然)塩釜 拓仁, 齋藤 大明, 川口 一朋, 長尾 秀実
GPGPU acceleration for calculation of energy distribution function(Kanazawa Univ.)Takuhito Shiogama, Hiroaki Saito,
Kazutomo Kawaguchi, Hidemi Nagao

1 はじめに
近年の計算機性能の向上によりタンパク質、その複合体、脂質膜などの大規模な生体系に対してのコン
ピューターシミュレーションが広く行われるようになってきた。しかし大規模系の計算では粒子数が大きいため計算コストも高い。したがってシミュレーションを高速化する手法やハード ウェアの開発が続けられている。現在開発されている次世代スパコンとして、GPU をアクセラレータとして使ったスーパーコン
ピューターが挙げられる。したがって CUDA や OpenCL 、OpenACC などの GPU を利用するための APIの開発が進められている。分子シミュレーションの中でも自由エネルギーの計算は特に計算コストが高い。
自由エネルギーを求める手法に自由エネルギー摂動法や熱力学積分法がある。これらの方法では多数の中間状態を作成しそれぞれに自由エネルギーを求めるための計算を行わなければならない。そこで松林らは溶液の分布関数理論に着目し自由エネルギーを低コスト、高精度で計算するエネルギー表示法を開発した[1-3]。この方法では溶液系のエネルギー分布関数を求める。このエネルギー分布関数を求める際に粒子間相互作用を計算する必要がある。そのため粒子数の二乗に比例して計算コストが増大する。そこで本研究では相互作用計算を CUDA Fortran を用いて高速化し 、性能の評価を行う。

2 エネルギー表示法
ここではある分子を溶質として溶媒に溶かす場合を考える。その時の溶質-溶媒間相互作用の 2 体ポテン
シャル v を予め設定しておく必要がある。v は溶質配置 ψ と溶媒配置 x の関数である。xi を i 番目の溶媒分子の座標とすると、溶質の周りの溶媒の瞬間的なエネルギー分布関数 ρ̂() はρ̂() =
で定義される。(1) 式中の∑N∑

δ(v(ψ, xi ) − )(1)iは N 個の全溶媒分子について和を取る。ここで溶質-溶媒間相互作用がある系の ρ̂() の統計平均を h・・・i 、無い場合の統計平均を h・・・i0 と表す。するとそれぞれの系の平均分布 ρ() 、ρ0 () はρ() = hρ̂()i, ρ0 () = hρ̂()i0(2)と書かれる。そしてエネルギー表示の密度汎関数理論を用いて定式化を行えば溶媒和自由エネルギー ∆µ は∫∫
ρ()∆µ = dρ()−kB T d[(ρ()−ρ0 ())−ρ()log()−α()F () + (1 − α())F0 ()(ρ()−ρ0 ())] (3)ρ0 ()
ここで α() は重み関数である。F () 、F0 () はシミュレーション上 ρ,ρ0 を厳密に求めることができないために導入された関数である。本研究ではエネルギー表示法の最初の段階であるエネルギー分布関数 ρ() 、ρ0 () を求めるプログラムについて高速化を施す。対象とする系は水分子 647 個を溶媒として水分子 1 個をその中に溶かすものを考える。水分子のモデルは TIP3P である。

3 高速化手法
ρ 、ρ0 を求める計算ではエワルド の方法によるクーロン相互作用の計算を含む。その部分が計算時間の90%を占めている。そのため、この部分に対して CUDA Fortran を使って高速化を施す。GPU はグリッド ,ブロック, スレッド といった階層構造を持っている。最小の単位はスレッドでありスレッドごとに同じプログラムが実行され 、データを GPU 側のメモリとホストの CPU 側のメモリと通信することによって計算結果を得る。従ってこの場合 do ループを GPU 上で総スレッド 数に分割し 、求めたデータをホスト側へコピーするといった手法をとる。 本研究では NVIDIA 社の Tesla C1060 を搭載した計算機、Tesla C2075 を搭載した計算機で性能評価を行った。Tesla C1060 のアーキテクチャは GT200 、Tesla C2075 のアーキテク
チャは Fermi である。なお MPI,cuda 並列いずれの場合でも通信に掛かる時間がボトルネックになりやすい。 MPI の場合ノード 間の通信はイーサネットを介して行われるため低速である。GPU の場合、ホスト

GPU 間のメモリ上のデータのコピーは、PCIExpress を介して行われるため低速である。従ってホストデバイス間のコピーは極力避けるべきである。なるべく通信の回数を減らすべきであり、通信されるデータの数も減らすべきである。したがって do ループの構造を変更しカーネル関数の呼び出しを最小限にするよう最適化を施した。並列化する手法としてはど ちらも粒子分割であり、単純に do ループの回数を全 cpu 数ま
たはスレッド 数で分割した。
100

ewald reciprocal
other

90

70

4000

80

3500

70

3000

ewaldreciprocal
other

60
60

40

2500

50

Time

Time

Speed up factor

50

40
30
20

1500

20

1000

10

10

2000

30

500

0
MPI8

0
0

10

20

30
40
Number of cpu

50

60

MPI16

MPI32

C1060

C2075

70

0
single

C1060

C2075

図 2: ρ を計算するプ ログ ラムの高速化。
図 3: ρ0 を 計算す るプ ログ ラ ムの 高速
図 1: 相互作用計算の高速化率。図中の直線 MPI8,16,32 は MPI 並列を 8, 16, 32 並列化。single は 1core での計算時間の内訳、は高速化率の理論値であり実際の値は▲で示 施した結果、C1060,C2075 は各々のアーキ
C1060,C2075 はそれぞれのアーキテクチャした。
テクチャを持った GPU を使用した結果をを持った GPU での計算時間の内訳である。
示す。

4 結果
最初に MPI による並列化に対する結果を示す。高速化率を 高速化率 = T1 /Tn で評価する。T1 は 1 プロ
セッサでの計算時間 Tn は n プロセッサでの計算時間である。図 1 は ρ() を計算するプログラムにおけるプロセッサ数に対する計算時間と高速化率である。次に GPU を用いた結果について説明する。図 2 は ρ() 、
図 3 は ρ0 () を計算するプログラムについての計算時間の内訳について MPI 並列と CUDA による並列を比較した結果である。図中の赤で示した部分はエワルド の方法によるクーロン相互作用の逆格子空間での計算時間である。緑で示した部分はその他の部分である。図 2 の結果から、MPI の並列数が増えるとその他の部分の計算時間が増大している事がわかる。これは並列数が増えることにより通信にかかる時間が増加するためである。図 2 、図 3 の結果から GPU のアーキテクチャの違いが大幅に計算性能を向上させる事が
わかる。

5 まとめ
MPI のみを使って高速化した結果から並列度を上げてもかなり早い段階で高速化率は飽和してしまうことがわかった。これは通信のオーバーヘッドが大きく影響しているからだと考えられる。GPU を使って高速化した場合について、アーキテクチャの違いが速度に大きく影響した結果となったが 、ρ() の計算は約60 倍、ρ0 () の計算については約 42 倍の高速化に成功した。

参考文献[1] N. Matsubayashi and M. Nakahara J. Chem. Phys. 113, 6070 (2000)[2] N. Matsubayashi and M. Nakahara J. Chem. Phys. 117, 3605 (2002)[3] N. Matsubayashi and M. Nakahara J. Chem. Phys. 119, 9689 (2003)

3P106多重極子とダンピング関数を用いた QM/MM 自由エネルギー計算（京大院理 1 、ワシントン大 2 ）○井上 雄介 1 、小杉貴洋 2 、林重彦 1

生体内の化学反応では、タンパク質の熱揺らぎが重要な役割を果たすと考えられており、その過程は自由エネルギー曲面上で議論されるべきものである。分子動力学 (MD) 法を用いてタンパク質を含む系の自由エネルギーを計算する際、数十ナノ秒以上のシミュレーションを実行して構造をサンプリングする必要があり、計算コストは非常に大きい。また、化学反応
は量子力学的な方法によって計算されることが必要であるが、量子計算によってタンパク質を含む系のサンプリングを取ろうとすれば、その計算コストはさらに膨大なものとなる。このような問題を解決する方法の一つとして、重要な部分を量子的に扱う QM 領域とし、それ以外の部分を分子力場を用いて古典的に取り扱う MM 領域とする QM/MM 法がよく知られている。QM/MM 法を用いれば、MM 領域内部の相互作用は低い計算コストで求めることが可能である。しかし、QM/MM 法を用いて QM-MM 間の静電相互作用を計算するとき、QM
領域の波動関数を用いたクーロン積分を求める必要があり、数十ナノ秒以上のシミュレーショ
ンを行う場合には、依然として計算コストは非常に大きなものになってしまう。この問題は、
QM 領域の波動関数のつくる電場を、相互作用点に置いた点電荷のつくる電場によって近似

する Restrained ElectroStatic Potential(RESP) 法を用いることによって解決される。RESP法を用いれば、QM-MM 間の静電相互作用も MM-MM 間静電相互作用と同様の低い計算コストで求めることができるため、十分なサンプリングをとることが容易になる。ところが、点電荷のつくる電場は等方的なものであるため、酸素原子や窒素原子等が持つ孤立電子対のつくる非等方的な電場を、点電荷近似によって表現することはできない。したがって、点電荷近似を用いた QM/MM 計算では、水素結合の方向を正しく予測することは困難である。また、電磁気学の基本法則によって予測できるように、非局在化した電子雲を点電荷で近似す
れば、実際よりも静電ポテンシャルを大きく見積もってしまう。電場の異方性の問題は、点電荷の代わりに点多重極子を相互作用点に置くことによって解決する。点多重極子を用いれば、異方的な電場を記述することができ、孤立電子対の向きを考慮した静電相互作用を低いコストで計算することが可能となる。電子雲の広がりの問題は、静電ポテンシャルの式にダンピング関数を導入することによって解決する。ダンピング関数は電子雲の広がりを表現す
るものであり、点電荷や点多重極子の場合よりも、QM 原子核近傍の静電ポテンシャルを弱める働きをする。以上により、QM-MM 間の静電相互作用は次の式によって計算される。
!
T
T
X
X
r
Q
r
q
p
r
i
ij
i
ij
ij
QM−MM
EES
=
qj fi (rij )+ i3 +
5
rij
rij
2rij
i∈QM j∈MM

ここでダンピング関数 fi (rij ) は、0 ≤ fi (rij ) ≤ 1 かつ fi (∞) = 1 を満たすものであり、本研
究では fi (rij ) = 1 − e−κi rij の形を採用した。また、多重極子は四重極子までの項を取り入れた。電荷 qi 、双極子 pi 、四重極子 Qi およびダンピング関数に含まれるパラメータ κi は、QM
領域の波動関数によって決まるものである。上の式を用いれば、点電荷近似と同程度の低い
計算コストで、電子雲の広がりと異方性を考慮した QM-MM 間静電相互作用の計算を実行することが可能である。
上述の方法によって QM-MM 間相互作用を効率良く計算することができるが、QM/MM法で MM 構造のサンプリングをとるとき、MM 構造の分布は QM 領域の状態に依存するた
め、QM 領域の波動関数や構造が変わる度に MD 計算を実行しなければならず、タンパク質等を含む巨大な系では、膨大な計算時間が必要になる。そこで、我々は Reweighting による
平均場 QM/MM 自由エネルギー計算法を開発した。この方法では、QM 領域の状態が変化したことによって変わった MM 構造の分布を、次の Reweighting の式によって補正することによって得る。

h noi
exp −β E QM−MM (d, r; R) − E QM−MM (d0 , r0 ; R)oiE ρMM (d0 , r0 ; R)
ρMM (d, r; R) = D h n
exp −β E QM−MM (d, r; R) − E QM−MM (d0 , r0 ; R)0

ここで ρMM (d, r; R) は、QM 領域の密度行列 d、QM 構造 r の下での MM 構造 R の分布を表す。したがって、理論的には最初の一回の MD 計算のみで、任意の QM 状態に対する MM 状態の分布を得ることが可能となり、非常に効率良く自由エネルギーを計算することができる。
多重極子とダンピング関数を用いた QM/MM 自由エネルギー計算法を用いて、水溶媒中のプロピオン酸などの計算を行った。プロピオン酸を QM 領域、周りの水を MM 領域とし、プロピオン酸の周りの水の分布を、点電荷で計算した場合と比較した。点電荷で計算した場
合には等方的な分布となったが、多重極子とダンピング関数を用いて計算すれば、酸素原子の孤立電子対の向きに高い分布を持つ結果となり、静電ポテンシャルの精度も大きく改善された。
発表では、計算方法の詳細と、他の計算結果についても述べる予定である。

3P-107

量子化学計算によるスクアレン分子の蛍光スペクトルの研究（慶大院・理工）小澤仁嗣，八代和徳，山本拓磨，藪下聡

Quantum chemical study on the fluorescence spectrum of squaraine molecule(Keio Univ.・Science and technology) Hitoshi Ozawa, Kazunori Yashiro, Takuma Yamamoto, Satoshi Yabushita・序
スクアレン分子は、長波長領域（630nm～750nm）の蛍光スペクトルを持つ色素分子で、生体分析、
色素増感太陽電池などに応用されている。Law は、この蛍光スペクトルが３つのピーク（短波長側か
らα、β、γ蛍光）をもつことを観測し、α蛍光は溶媒和されていない平面構造からの発光、β蛍光
は溶媒和されている平面構造からの発光、γ蛍光は TICT(twisted intramolecular charge transfer)励起状態からの発光と帰属した。我々は、LC-(TD)BOP/6-31G(D)法を用い、この分子のねじれ構造の励起状態の理論計算から、実際に TICT 励起状態は存在するが、S1 状態ではなく、その励起エネルギーは実験値に比べて高すぎること、さらに TICT の原因となる内部回転に必要な反応熱の計算から、TICT 励起状態の関与を否定した。他方 Rettig らは、この蛍光スペクトルがシングレットであると主張し、長波長側の蛍光スペクトルの肩は、溶媒和されている平面構造の 0-n 遷移によるものと帰属した。Rettig らのモデルにしたがうと、平面構造の吸収、蛍光スペクトルのストークスシフト（SS）の大きさが重要な意
味を持つが、(PCM)-(TD)BOP/6-31G(D)法でこの大きさを評価したところ、実験値をかなり過大評価した。このため、まず MRMP2/6-31G(D)法、CASSCF(11,10)/6-31G(D)法を用いた一連の計算により、スクアレン分子の平面構造の基底電子状態はジラジカル構造（図１右）よりも双極イオン対構造（図１左）が支配的であることを確認した。次に、SS の過大評価は、密度汎関数法が、電子基底状態と励起状態における平衡構造の違いを過大評価していると考え、MRMP2/6-31G(D)法を用いて中央四員環部分の構
造最適化計算を行ったところ、基底状態と励起状態の平衡構造の違いは、密度汎関数法で評価した場
合よりも大幅に減少し、実験結果と対応することが分かった。今回は、この一連の研究結果を発表し、更にこの構造を用いた吸収・蛍光スペクトル、無輻射遷移速度の温度依存性の検討を行う。

図１．スクアレン分子の構造式[1],[2]・構造最適化計算
電子基底状態(Ag)、最低励起状態(B1u)の平衡構造を(TD)BOP/6-31G(D)法、LC-(TD)BOP/6-31G(D)法、MRMP2/6-31G(D)法で最適化した場合の、中央四員環部分の結合パラメータ及び励起エネルギーを表 1に示す。ここで R1, R2 は図２のように定義した。MRMP2 法では、四員環以外の部分は、
R2

R1
図２． 中央 4 員環部分の結合パラメータ（置換基は省略）(TD)BOP/6-31G(D)法の平衡構造に固定し、R1 R2 をともに 0.01Åずつずらしたポテンシャルエネルギー面のスキャンにより極小点を決定した。さらに基底状態の平衡構造における垂直励起エネルギーΔEabs、励起状態の平衡構造における垂直励起エネルギーΔEemi、およびそれらの差から評価した SS を示す。
表 1．

基底状態(Ag)、励起状態(B1u)の構造最適化計算（最初の R1、R2 が基底状態の、次の R1、R2 が励起状態の構造パラメータ。SS: ストークスシフト）
計算方法

R1[Å]

R2[Å]

R1[Å]

R2[Å]

ΔEabs[eV]ΔEemi[eV]SS [eV](TD)BOP/6-31G(D)

1.047

1.065

1.073

1.034

2.14

1.98

0.16

LC-(TD)BOP/6-31G(D)

1.023

1.050

1.055

1.032

2.42

2.27

0.15

MRMP2/6-31G(D)

1.029

1.054

1.042

1.031

1.97

1.89

0.08

実験値

－

－

－

－

1.97

1.91

0.06

表 1 を比較すると、基底状態ではいずれも R1< R2 である。また、励起状態ではいずれも R1 > R2 とな
るが、MRMP2 法では、R1 の値が小さくなり、R1 と R2 の値が近くなることが分かる。さらに密度汎関数法のものに比べ、MRMP2 法の垂直励起エネルギーから求めた吸収・蛍光スペクトルのピーク位置お
よび SS の値は実験値とよく対応した。以上のことから、スクアレン分子の基底状態と励起状態の構造のずれ方が、相関エネルギーの評価方法の違いに大きく依存していることを示すことができた。今後、吸収・蛍光スペクトルの計算と無輻射遷移確率の温度依存性の計算などでは、密度汎関数法よりも相関エネルギーを実質的に含む計算方法で検討を行う予定である。
量子収率及び蛍光寿命の測定から、温度が上がるにつれて無輻射遷移確率は増加することが観測されている。我々はこれまで温度依存型 Engleman-Jortner の式を用いて、無輻射遷移確率の大よその値を評価してきたが、温度依存性など定量的な振る舞いは再現できていない。今後は無輻射遷移確率を、
CASSCF(11,10)/6-31G(D)法で非断熱相互作用項の計算を行い、フェルミの黄金律を直接用いて、無輻射遷移確率を評価する予定である。

・参考文献[1] K. Y. Law, J. Phys. Chem., 91, 5184, (1987).
[2] C. Comelissen-Gude, W. Rettig and R. Lapouyade, J. Phys. Chem. A., 101, 9673, (1997).

3P­108
量子古典分子動力学シミュレーションによるデコヒーレンスの評価(高等師範学校 1, カンザス大 2) 城塚

達也 1, Damien Laage1, Ward H. Thompson2

デコヒーレンスは遷移速度やスペクトルに大きく影響するが、半古典近似では考慮されていない.また、デコヒーレンスの振動遷移への影響も明らかというわけではない.そこで、半古典的な結果への補正は摂動論（フェルミの黄金則）に基づき導入する.ここで、１つの自由度に対してシュレーディンガーを解き、他の自由度（核）の波動関数に対してガウス近似を用いることにより表式を導いている.[1, 2]この手法は電子スペクトルの解析などでしばしば用いられたが、本研究では水（重水も含めて）の振動遷移（スペクトル）に応用し、[3]非断熱遷移へのコヒーレンスの影響を量子分子動力学（MD）シミュレーションを用いて研究する.この近似を用いると遷移速度やスペクトルが半古典の結果とデコヒーレンス関数の積の形になる.よって、まず MD によりデコヒーレンス関数を正確に評価する.また、デコヒーレンス関数への様々な近似[4]を考え、各成分の物理的意味を探る.更に、短時間近似を用いてデコヒーレンス関数の簡単な近似を調べる.最後に、本手法を IR スペクトルに適用しデコヒーレンスの影響を検証する.
分子間に SPC/E ポテンシャルを用いて 256 個の水分子のシミュレーションを行った.MD
シミュレーションの毎時間ステップに 1 個の OH に対して時間に依存しないシュレーディンガ
ー方程式を解き、Hellman-Feynman 力により時間発展させた.このとき、PO-DVR を用いて計算で用いるグリッドの数を減らした.[5]また、OH 伸縮にはモースポテンシャルを用いた.[6]デ
コヒーレンス関数の計算で用いる非平衡 MD の初期構造は振動基底状態の平衡化シミュ
レーションから 480 個をサンプルし振動基底・励起状態で 400 fs 時間発展させた.デコヒーレ
ンス関数を並進と回転成分に分解する際、Cartesian 座標から分子内座標に変換した.この時、回転座標は水分子の慣性モーメントテンソルが対角的になるよう選んだ.
得られたデコヒーレンス関数を図１に示す.デコヒーレンス関数の実部をガウス関数でフィッ
ティングすると約 33 fs の緩和時間が得られた.座標変換によりデコヒーレンス関数を成分分解するとこの緩和に並進と回転の両方が関与していることがわかった.また、デコヒーレンス
関数への様々な近似はいずれも定量的に十分とは言えないが、短時間近似はデコヒーレン
ス時間を粗いが簡単に見積もれることがわかった.最後に、本手法を IR スペクトルに応用す
ると半古典の結果と比べて線幅が大きくなった.

Figure 1: デコヒーレンス関数の実部と虚部.実部への並進と回転の寄与も示す.凡例中の ET と HT はガウス幅の正確な表式と高温近似を用いた結果を表す.
[1] Eric J. Heller, J. Chem. Phys. 75, 2923 (1981).
[2] Oleg V. Prezhdo and Peter J. Rossky, J. Chem. Phys. 104, 5863 (1997).
[3] D. Kraemer et al, PNAS 105, 437 (2008).
[4] Laszlo Turi and Peter J. Rossky, J. Chem. Phys. 120, 3688 (2004).
[5] Ward H. Thompson, J. Chem. Phys. 118, 1059 (2003).
[6] Kahled Toukan and Aneesur Rahman, Phys. Rev. B. 31, 2643 (1985).

3P109Theoretical investigation of the structures and dynamics of moleculargyroscopes encapsulating halogen-substituted phenylene rotators(1Department of Chemistry, Graduate School of Science, Tohoku University,
2
Faculty of Pharmaceutical Sciences, Tokushima Bunri University)
Anant Babu Marahatta1, Manabu Kanno1, Wataru Setaka2 and Hirohiko Kono1

I. Introduction: The crystalline molecular
Statorgyroscopes encapsulating phenylene rotators have
Rotatorreceived much attention in terms of their uniquestructures, dynamics and functions. Especiallyhalogen-substituted phenylene-bridged molecularF
F
gyroscopes are expected to have dipolar units thatF
can be reoriented under the influence of electricF
fields [1]. Recently, Setaka et al. have reported theSpinsynthesis and structure of the novel molecularaxis
gyroscopes having halogen-substituted phenylene
rotators encased in three long siloxaalkane spokes Fig.1. An X-ray crystallography of the(Stator) [2]. They observed two different stable siloxaalkane molecular gyroscope having a
positions for the difluorophenylene rotator at difluorophenylene rotator at 273 K. Two
273 K (Fig. 1) and only one position for the stable positions of the rotator differ by πdichlorophenylene rotator at 153 K. Similarly, a flipping are shown.
strong deformation of the siloxaalkane spokes caused by the steric interaction with the halogen
atoms on the central phenylene ring has also been reported. During the course of theoreticallyinvestigating crystal structures and rotational dynamics of experimentally synthesized crystalline
molecular gyroscopes, we are interested in a crystalline molecular gyroscope having a
difluorophenylene rotator (ROT-F). The notable results explored currently are reported here.
II. Method: Recently, we applied the Density Functional based Tight Binding (DFTB) method[3] to investigate the molecular packing structure and rotary dynamics of the crystalline
siloxaalkane molecular gyroscope having a phenylene rotator (ROT-H) [4]. We achieved good
agreements with the experimental results. Therefore, in our present study, we have also applied
the same theoretical method for characterizing crystalline ROT-F. In accordance with theexperimental condition, a periodic boundary condition (PBC) is imposed in the calculations. By
using the force constants computed by DFTB, the motions of all nuclei are evaluated based on
the classical mechanics with constraints on the lattice constants of the unit cell.
III. Results and Discussions: The X-ray geometries of theunit cells of two degenerate molecular structures of ROT-F are
semiquantitatively reproduced by the DFTB method except
the Si-O-Si angles of the highly flexible siloxaalkane spokes(see Table 1). These angles are more acute in ROT-F than in
ROT-H derivative. It indicates that the siloxaalkane-cage of
the ROT-F expands outward with respect to the oxygen atom
for accommodating more bulky rotator. In order to confirm it,
we measured the intervening space that exists between an
Fig.2. A schematic representation
oxygen atom of each siloxaalkane arm and the nearest carbonof the siloxaalkane molecularatom of the phenylene ring (O-C distances: d1, d2 and d3 in Fig. gyroscope. The encircled atoms
2) for ROT-F and ROT-H and found that ROT-F possesses
define the dihedral angle. X = H
for ROT-H and X = F for ROT-F.

Table 1. X-ray vs. DFTB optimized parameterscomparatively large amount of
Molecular
Si-O-SiO-C distance (Å)
Phenylene dihedral
clearance around the rotating gyroscopes
angle (π)
angle ( π)
X-ray DFTB
X-ray
DFTB
X-ray
DFTB
unit (see Table 1). It assuresthat the shape of the
7.7
0.96
0.78
d1
8.4
siloxaalkane-cage in ROT-F
ROT-F
0.84
0.75
d2
5.9
5.5
0.56
0.44
5.9
5.8
0.84
0.74
d
3
deforms from that of the ROT5.3
0.96
0.81
d1 5.2H derivative to avoid steric
4.1
0.55
0.35
ROT-H
0.90
0.77
d2 4.3
effects with the fluorine on
4.1
0.96
0.81
d3 3.9phenylene. To estimate arotational energy barrier (Ea) of
the ROT-F, we scanned thepotential energy surface (PES)
under PBC (Fig. 3) and found thatEa
the difluorophenylene rotator–1experiences as high as 1700 cm
α′energy barrier (Ea of theαphenylene rotator in ROT-H
250 cm-1 [4]). The PES clearly
displays two stable degenerate
molecular structures, α and α′ andFig.3. Potential energy as a function of the difluorophenylene angle
theirrelationbyπobtained by DFTB under PBC.
difluorophenylene flipping. It
confirms the two-fold rotational axis of the ROT-F molecule with respect to thedifluorophenylene rotator. These theoretical findings are in reasonable agreements with theexperimental observations. Moreover, we performed DFTB / molecular dynamics (MD)simulations of the ROT-F at 1200 K and 300 K, and calculated the concerned dihedral angle thatdesignates the difluorophenylene orientation at each MD step. At each kinetic temperature, theinitial angle of the difluorophenylene (~0.44π) corresponds to its equilibrium position α asshown in Fig. 4 and from this position; difluorophenylene ring is observed to flip by π andreaches to its degenerate position α′ (~1.44 π). The timescale of the π flipping at hightemperature is found to be shorter than at moderate temperature, as expected.
At 1200 K

At 300 Kα′α′
π flipping

π flippingααTime (ps)Fig.4. Rotational angle of difluorophenylene as a function of time at 1200 K and 300 K.

IV. Conclusion: We successfully located two degenerate molecular structures of the ROT-F
under crystal condition. They are found to be related to each other by π difluorophenylene
flipping and have to cross ~1700 cm-l energy barrier to undergo structural exchange between
them. The flipping motion of the difluorophenylene inside a cage of the siloxaalkane spokes athigh and moderate temperatures indicates a promising function as a molecular gyroscope.
References : [1] Garcia-Garibay et al. J. Am. Chem. Soc., 124 (2002) 7719. [2] W. Setaka et al. Chem.
Lett. 39 (2010) 468. [3] D. Porezag et al. Phys. Rev. B, 51 (1995) 12947. [4] A. B. Marahatta et al. (to besubmitted).

3P3P-110
透明導電酸化物としての
透明導電酸化物としての、としての、F ドープした
ドープしたアナターゼ
したアナターゼ型酸化アナターゼ型酸化チタン
型酸化チタンのチタンのGGA+U 計算による計算による検討による検討（東大院工 1，東大院理 2）水口菜々子 1，神坂英幸 2，山下晃一 1

First-principles GGA+U Calculations of F-doped anatase titanium dioxides
for transparent conductive oxides applications(Graduate School of Eng., The University of Tokyo1，Graduate School of Sci., The University of Tokyo2)
Nanako Mizuguchi1，Hideyuki Kamisaka2，Koichi Yamashita1

1.

緒言
透明導電酸化物(Transparent Conductive Oxides : TCOs)は、透明かつ電気を通すという特異な性質を持った化合物であり、太陽電池やディスプレイの部材としてその需要を拡大している。しかし、現在主流の TCOs は酸化インジウムスズ(Indium Tin Oxides : ITOs)と言う化合物であり、希少金属のインジウムを使用するため、価格の高騰及び枯渇が懸念されている。
そこで、
ITOs に変わる TCOsとして注目を集めているのが、アナターゼ型 TiO2 系の TCOs である。これは、バンドギャップが可視光の最大エネルギーである 3.2eV より大きく可視光を吸収しないアナターゼ型 TiO2 に、ドー
プをすることによって導電性を高めて実現する。現在、Ti を Nb で置換することによって ITOs に匹敵する性能を得られるという研究報告が存在する。1)一方で、TCOs への応用を目的とする O の置換についての研究は、未だ大きな成果を挙げていない。この O 置換によるアナターゼ型 TCOsの実現を目指すにあたり、本研究では、周期表で O の１つだけ右側に位置する F をドープ元素として検討した。F は O に置換し、系内に余剰電子を 1 つ放出すると予想され、この余剰電子の拡散により導電性が高まることが期待される。しかし現実には、合成された F ドープ TiO2 は期待されたほどの性能を示さないという実験結果が報告されている。
この理由には、大きく次の 3 つが考えられる。それは、(i)誤った位置へのドープ、(ii)電子を補足するような欠陥の生成、(iii)電子を補足するような状態の存在、である。これらについて、密度汎関数理論(Density Functional Theory : DFT)に基づく第一原理計算を用いて検討した。

2.

計算手法
アナターゼ型 TiO2 の 8 倍セル（Ti16O32）に、1 つの F 原子を置換または添加した構造（ドープ率 6.25%）を計算セルに用い、一般化された密度勾配近似(Generalized Gradient Approximation :GGA)の計算レベルで DFT に基づく第一原理計算を行った。汎関数は、エネルギー計算に PW91
汎関数を、その他の計算に PBE 汎関数を用いた。PAW 型偽ポテンシャルを用い、Kohn-Sham 軌
道は平面波基底によって展開した。
また、特に要素(iii)の電子を補足するような状態の存在を正確に検討するために、Ti 3d 軌道、O2p 軌道、F 2p 軌道に GGA+U 法を適用した。これは、局在的な性質を持つ原子軌道に入っている電子の反発をよりよく表現するために取り入れられるものであり、適用により自己相互作用誤差を削減することが出来る。既往の研究において、このパラメタの適用により電子を補足する準位を再現できたという報告があり 2)、それらに基づいて本研究でもこの適用を行っている。+U の値を出す際に、今回は 2 種類の方法によるものを用いた。１つは Koopmans の定理を満たすよう+Uの値を定める方法で 2)、もう 1 つは原子のイオン化エネルギーに基づいて+U の値を算出する方法
である 3)。このうち、前者の方法による値は既往の研究のものを引用し、+U(Ti 3d)=3.40eV4)、+U(O2p)=5.25eV2)とした。後者の方法について下に図示する。
図の左側の場合、電子は他の４つの電子から反発を受ける。従ってこの状態
での電子反発によるエネルギーは 4U-である。同様に、右側の場合は 3U で

3U

4U
図１ +U 算出の算出のスキームある。この差を取ることで、問題とする軌道の+U の値を算出する。その結果、(F 2p)=4.28eV@Osite と算出された。

3.

結果・
結果・考察
まず、各ドープ位に F ドープした TiO2 の、安定構造のエネルギーに関する結果を表１に示す。

ここで、欠陥生成エネルギーΔE は欠陥（ドープ）が出来るときの全エネルギー変化であり、熱力 学 的 安 定 性 の 見 積 も り に 利 用 出 来 る 。 Δ E=[E(doped TiO2)+E(replaced atom)]-[E(pureTiO2)+E(doped atom)]であり、この値が小さいほど有利である。表より、期待していた通り O 位に置換するのが最も熱力学的に有利であることが確認された。また、スペースの関係上構造は O 位置換のものしか掲載していないが、純粋な TiO2 と比べての構造変形も O 位置換が最も小さかった。
O 位置換

Ti 位置換

格子間添加

全エネルギーE-384.6-365.1-381.5

欠陥生成エネルギーΔE-0.152

+21.0

+1.53

（単位は全て eV）

表１ 各ドープ位
ドープ位のエネルギー比較
エネルギー比較

続いて、O 位 F 置換 TiO2 の電子密度(Density Of States : DOS)図を図２に示す。適正な Hubbard +Uの値を用いて spin 分極させて計算することにより、伝導帯下端付近に不純物準位が確認出来た。
また図３に、確認された不純物準位での電子密度図を示す。中央が F 原子であるが、その周りに電子がやや局在している。以上より、F ドープ以外の要素によるフェルミ準位の変動が無視できる場合、不純物準位が電子を捕捉し、F ドープ TiO2 の導電性を低下させる可能性がある。これが、実験的に合成された F ドープ TiO2 が期待通りの性能を示さないことの一因と考えられる。

図２ O 位 F 置換 TiO2 の DOS 図

図３ 不純物準位付近の不純物準位付近の電子密度図

参考文献
1) Furubayashi, Y.; Hitosugi, T.; Hasegawa, T. et al. App. Phys. Lett. 2005, 86, 252101.
2) Morgan, B. J.; Watson, G. W. Phys. Rev. B 2009, 80, 233102.
3) Janotti, A.; Segev, D.; Van de Walle, C. Phys. Rev. B 2006, 74, 045202.
4) Kamisaka, H.; Mizuguchi, N.; Yamashita, K. J. Mat. Sci., DOI 10.1007/s10853-012-6491-2 (online 27 April2012).

3P-111
ドナーアクセプター型ポリチオフェン誘導体の電荷移動と失活過程（上智大学 1、京都大学福井謙一記念研究センター2）
小林 理 1，赤間 知子 1，石田 俊正 2，南部 伸孝 1
Theoretical study of charge transfer and deactivation process ofdonor-acceptor type thiophene derivatives
1(Sophia Univ. , Fukui Institute for Fundamental Chemistry, Kyoto Univ.2)
Osamu KOBAYASHI1, Tomoko AKAMA1, Toshimasa ISHIDA2, Shinkoh NANBU1【序論】有機薄膜太陽電池は成型をポリマー溶液の塗布などにより容易に行うことが可能であり、無機材料に比べて軽量であるため新時代の電子材料および光機能材料として注目されている。一方で発電効率の点からは無機太陽電池が優れており、有機薄膜太陽電池にはさらなる発電効率の向上が期待される。一般に太陽電池は正孔が過剰な部分である p 型半導体と電子が過剰な部分である n 型半導体を接合した構造の素子を含み、有機薄膜太陽電池では p 型部位の材料としてポリチオフェンのような π-共役高分子、n 型部位の材料としてフラーレン誘導体などが用いられる。
特にポリチオフェンは化学的、熱的安定性に優れるという材料の実用化において重要な特長を有す材料として、近年その合成法の検討とデバイスへの応用が進められている。太陽電池の発電は素子が光を吸収することで生じた励起子が接合界面において電荷が分離し電極まで運ばれることによる。励起子には寿命が存在し、その間に電荷が分離されることが発電において重要である。
電荷移動は分子内で起こる方が分子間で起こるよりも有利だと考えられ、電荷移動を分子内で起こすことが出来ればさらに発電効率を向上させることが期待できる。
我々は以前の研究においてドナー部位とアクセプター部
位を併せ持つ分子を ab initio 計算により探索した。計算に使用した分子が分子内電荷移動を起こすことは計算により得られる分子軌道より判別することが可能である。詳細にはHOMO がドナー分子上に、LUMO がアクセプター分子上に局在化することにより確認できる。2-(2’-(4’-cyano)thienyl)-pyrrolate anion (CNTP )が分子内電荷移動を起こすことを理Fig. 1 Structure of CNTP-.

論的に明らかにしている。
光電素子の性能は材料の光励起後の挙動に影響され、分子内の失活過程が効率的な電荷分離を阻害する一方、分子間の電荷移動が電荷分離を促進する。原子核の運動により電子状態が変化す
る非断熱現象は二つの電子状態のポテンシャルが近接したときに起こる。従って電子状態がバンド的に存在する材料では非断熱遷移が重要な役割を果たしていると考えられ、特にアクセプター
部位上の分子軌道からドナー部位上への分子軌道への遷移は失活過程に、異なる分子上の分子軌
道間の非断熱過程は分子間電荷移動に寄与していると考えられる。従って、これらの電子状態間の非断熱過程を調べることは重要である。
【方法論】非断熱 ab initio 分子動力学シミュレーションを行うために、Zhu-Nakamura 公式[1]に基づく Tully’s Surface Hopping (TSH)法[2]を用いた。電子基底状態及び励起状態のポテンシャルエ

ネルギー及び傾き、そして非断熱結合要素を complete active space self-consistent field (CASSCF)法により求めた。また、基底関数には Dunning の cc-pVDZ を用いた。さらに、得られたデータを用い、速度 Verlet 法に基づき原子核の座標および運動量の時間発展を行った。
次に、電荷移動反応と周囲の分子との影響を考慮するために QM/MM 法を用いた。QM 層には-

2-(2’-(4’-cyano)thienyl)pyrrole (CNTP0)を 1 分子と CNTP を 1 分子割り当て、MM 層は多数の CNTP0
及び Li+により構成した。初期構造は古典的な分子動力学シミュレーションを 120.0 fs 行って作成した。QM 層における分子内及び分子間電荷移動状態を求める量子化学計算の方法は上記と同様であり、分子力場法においては AMBER 力場 ff03 を用いた。QM/MM 計算によりポテンシャルエ
ネルギー及び傾き、そして非断熱結合要素を同様に求め、QM/MM-ZN-TSH 計算を実施した。
【結果】単分子の分子動力学シミュレーションの結果からは非断熱過程は発見されなかった。これは S1 状態が極めて安定であることを意味し、従って無輻射過程による失活は起こらないと結論できる。
QM/MM モデルにおいて、分子内電荷移動状態に対応する電子状態は S3 状態であった。S0 状態と S3 状態への遷移双極子モーメント 8.280 Debye であり、S3 状態への励起は起こりやすいと考えられる。したがって光吸収により S3 状態へ励起されると考えられる。励起エネルギーは 3.414 eV-であった。S1 状態と S2 状態はそれぞれ CNTP の HOMO から CNTP0 の LUMO と LUMO+1 への励起に対応する電子状態であり、これらの状態への非断熱遷移が分子間電荷移動を表すと考えられる。S3 状態から開始した QM/MM モデルに基づくシミュレーションからは開始数十 fs の時点において分子間電荷移動に対応する非断熱現象が確認できた。

(a)

(b)

(c)

(d)

Fig.2

(a) HOMO of CNTP in QM/MM model. (b) LUMO of CNTP0 in QM/MM model--(c) LUMO+1 of CNTP0 in QM/MM model (d) LUMO of CNTP in QM/MM model.
【参考文献】[1] Hiroki Nakamura, Nonadiabatic Transition (2nd Edition) (World Scientific, 2012)[2] John C. Tully, Richard K. Preston J. Chem. Phys. 1971, 55, 562-572

3P-112
PtCO, PdCO 分子のマトリックス振動分光における特異的希ガス効果の理論的解明（北大院理 1, 埼玉大院理 2）小野ゆり子 1、中山哲 1、高柳敏幸 2、武次徹也 1

Theoretical study of noble-gas binding effecton vibrational spectroscopy of MCO (M = Pd, Pt)(Hokkaido Univ1., Saitama Univ2.) Yuriko Ono1, Akira Nakayama1,
Toshiyuki Takayanagi2, Tetsuya Taketsugu1【 序 】 我々はこれまで NiCO 分子について高精度 ab initio 計算に基づきポテンシャル曲面の非調和性を考慮した振動準位を求め、分光測定結果との比較を行ってきた。その結果、特に変角振動モードの振動数の計算値が Ar マトリックス赤外分光法による実験値とずれることを見出し、その原因を探求する中で、Ar が NiCO に約 9 kcal/mol のエネルギーで結合し、変角振動数が 10%程度青方遷移して実験値と非常に良い一致を見せることを明らかにした。この計算に基づき、マトリックス分光法で測定された実験値は NiCO 分子ではなく希ガス化合物（Ar-NiCO）に帰属されるべきことを指摘した[1]。一方 Ni と同族金属の錯体である PdCO,
PtCO についても同様の計算を行ったが、変角振動数のマトリックス分光による実験値は ab
initio 計算で得られる値(273, 419 cm-1)の約 2 倍(615.7, 916.8 cm-1)であり、希ガスと結合する可能性を考慮しても依然として実験値と計算値の大きな不一致を説明することはできなかった[2]。2011 年 3 月に、NiCO, PdCO, PtCO に対する気相での分光実験の結果が報告されたが[3]、この論文では我々の理論計算の結果が引用され、振動数について実験値と計算値がよい一致を示すことが述べられ、希ガスマトリックス分光実験では変角振動の倍音が観測されているのではないか、との議論が展開された。本研究では、希ガスマトリックス分光実験の結果を原理的に理解するため、振動 CI 法を用い理論計算により振動数と強度を求めることとした。
【 計 算 】対象分子の各座標を変化させながら電子状態計算を繰り返すことにより得られたポテンシャルエネルギー曲面をもとに振動の Schrödinger 方程式を解けば、量子化された振動準位が得られる（ここで、振動 SCF 法、振動 CI 法を適用する）。この振動準位のエネルギー差
から、基本振動数や倍音、結合音が得られる。平衡構造におけるエネルギーの 2 次微分から得られるのは調和振動数であり、実験値との直接比較には、3 次以上の効果である非調和性を考慮して振動数を算出する必要がある。振動 CI 法には対象分子の正確なポテンシャル曲面が必要となるが ArPtCO 分子に対し CCSD(T)計算または QCISD(T)計算を行うと Pt-C の距離を伸ばした位置で収束性が悪化するので、この手法による滑らかなポテンシャル曲面の算出は困難である。また MR-SDCI(+Q)計算はコストの面から現実的とは言えない。そこで本研究においてはまず PtCO 及び PdCO 分子に対し、DK3-QCISD(T)法による全電子計算を行い平衡構造及び振動数を算出して実験値との比較を行い、次に約 40 種類の汎関数と 10 種類の基底関数を組合せて行った DFT 計算結果から平衡状態近傍での DK3-QCISD(T)法による結果及び実験値をよく再現する汎関数と基底関数の組合せを抽出し、選んだ手法を用いポテンシャル

曲面を求めた。DK3-QCISD(T)法において Pd 及び Pt 原子の基底関数には、3 次の Douglas-Kroll法で最適化された土屋らの相対論用基底系を縮約し、関谷らの相関関数を加えた基底系を用い、C と O に関しては、関谷らの sapporo-qzp を用いた。本研究においてはスペクトル強度の算出が必須であるため、双極子モーメントの曲面も同時に作成する。
【 結 果 及 び 考 察 】種々の汎関数と基底関数の組合せの計算手法（約 420 種類）を用いて PtCO
分子に関して構造最適化及び調和振動数計算を行い、算出された値を実測値 r(Pt-C) = 1.760 Å,r(C-O) = 1.144 Å, ν(Pt-C-O) = 420 cm-1, ν(Pt-C) = 600 cm-1 と比較し、これをもっとも再現する汎関数と基底関数の組合せを抽出した。PBE0MOL(0.75PBEx+0.25HF+1.0PW91c)/def2-TZVPP に
よる PtCO 分子の計算値は r(Pt-C) = 1.760 Å, r(C-O) = 1.143 Å, ν(Pt-C-O) = 416 cm-1, ν(Pt-C) =599 cm-1 であり DK3-QCISD(T)法による計算値 r(Pt-C) = 1.760 Å, r(C-O) = 1.146 Å, ν(Pt-C-O) =419 cm-1, ν(Pt-C) = 605 cm-1 と比較しても遜色ない。Table1 に PBE0MOL/def2-TZVPP により算出されたポテンシャルに VCI 法を適用して算出した PtCO, ArPtCO 分子の調和振動数・基音・倍音を示す。表より、PtCO 分子のν(Pt-C-O)は 412 cm-1 であり、その倍音は 826 cm-1 であるが倍音の強度は基音の約 1/30 であり、実験において基音と見間違う程の強度はない。比較して
ArPtCO 分子の場合ν(Pt-C-O)は 453 cm-1 であり、その強度は 0.1 と非常に弱い。一方ν(Pt-C-O)の倍音 906 cm-1 の強度は基音の 10 倍と倍音においては異常なほど強い。この振動数は既にPtCO のν(Pt-C-O)として報告されている 916 cm-1 とも近い。以上の結果より、Ar が結合し新たにフェルミ共鳴が生じることにより強度の異常が起こり、ν(Pt-C-O)の振動数の同定ミスが起こったと考えられる。また同時に、PtCO 分子は Ar マトリックス中で ArPtCO 分子を生成し、赤外分光測定では ArPtCO を観測していたとの我々の推定をも裏付けている。同様に PdCO、
ArPdCO 分子についても解析を行った。詳細な結果は当日発表する。

Table 1. Calculated harmonic, fundamental and overtone frequencies for PtCO and ArPtCO.
Frequency / cm-1 (Intensity)

ν(CO)
PtCO

ν(PtC)

ν(Pt-C-O)

Harmonic

2153 (

632.2 )

599 (

7.4 ) 416 ( 3.1 )Fundamental

2124 (

617.4 )

589 (

6.9 ) 412 ( 2.9 )Overtone4225 (

3.6 )1174 (

0.0 ) 826 ( 0.1 )ArPtCO Harmonic

2154 (

816.9 )

587 ( 15.2 ) 460 ( 0.1 )Fundamental

2125 (

797.6 )

580 ( 15.2 ) 453 ( 0.1 )Overtone4228 (

0.0 )1155 (

0.0 ) 906 ( 0.9 )【参考文献】[1] Y. Taketsugu, T. Noro, and T. Taketsugu, J. Phys. Chem. A, 112, 1018 (2008).
[2] Y. Taketsugu, T. Noro, and T. Taketsugu, Chem. Phys. Lett., 484, 139 (2010).
[3] T. Okabayashi, T. Yamamoto, E. Y. Okabayashi, and M. Tanimoto, J. Phys. Chem. A, 115, 1869
(2011).

3P113分子集団の相対配置の最適化(和歌山大システム工 1、豊田理研 2) 山門英雄 1、澤田裕 1、大野公一 2

Optimization of relative arrangements of molecules(Wakayama Univ.1, Toyota Phys. Chem. Res. Inst.2） Hideo Yamakado1,
Yu Sawada1, Koichi Ohno2【序】 筆者等は、これまで超球面探索法（SHS 法）1) を結晶構造の予測に適用し 2)、その有効性を実証してきた。 本研究では、超球面探索法に要する計算時間を大幅に短縮することを目的として、原子が集まって形成する分子集団の形を固定し、その重心の並進移動と回転の自由度を与えることにより分子集団の相対配置を最適化し、また複数の局所的安定構造を自動的に探索することを試みた。このことは、超球面探索法 1) に分子骨格を固定する束縛条件を導入することに相当する。
例 え ば 、 束 縛 条 件 や 並 進 以 外 の 対 称 性 を 一 切 導 入 し な い 場 合 、 図 １ に 示 す 実 在 す る TTF(tetrathiafulvalene: C6S4H4) – TCNQ (tetracyanoquinodimethane: C12N4H4) 結晶の場合、ユニットセル内には TTF 分子と TCNQ 分子が 2 個ずつあり全原子数は 68個、そして a、b、ｃ軸の自由度も含めて考えると、重心の並進と回転を差し引いても、3×(68+3)-6= 207 変数について最適化する必要が生じる。一方、各分子の構造を固定して、その相対配置とユニットセルの最適化を行う場合であれば、全自
由度は 27 変数まで抑えられ、これは気体分子で言えば 11
原子分子の構造最適化に相当し、SHS 法を用いて現実的な時間内に探索を行うことが可能になる可能性がある。
図１.

TTF-TCNQ の結晶構造【方法】 安定構造の探索には多原子系の基準座標を使用するように開発された
超球面探索法を、多変数関数に対する Hessian の固有ベクトルを使うように一般化した方法 3)を用いた。2 個のホルムアルデヒド分子(HCHO)について、その相対配置を探索した例を以下に示す。エネルギーの計算には、Gaussian09 を用い、計算レベルは MP2/6-31G で行った。各分子の形は、単独のホルムアルデヒド分子について上記のレベルで構造最適化を行った後に固定し、1 分子は空間に固定、もう一つの分子はその重心に 3 個、回転に 3 個の自由度を与えた上で、これら 6個の独立変数を上記の一般化された超球面探索法により安定構造の自動探索を行った。なお、初期構造としては、2 つのホルムアルデヒド分子を分子面に垂直な方向に 5 Åずらして積み重ねた構造(図 2)から出発した。
図２. 初期構造の配置【結果と考察】 図 3 に、分子骨格を固定して自動探索された二つの安定構造((a),(b))を示す。(a)は、 1
個目のホルムアルデヒド分子の酸素原子は、2 個目のホルムアルデヒド分子の 2 個の水素原子と近く、2個目のホルムアルデヒド分子の酸素原子は、1 個目のホルムアルデヒド分子の 1 個の水素原子と近い形をしており、また(b)は 2 個のホルムアルデヒド分子が同一平面上に存在し、それぞれの 1 個の水素原子が、相手分子の酸素原子と接近した形をしている。この段階では、(a)の方が(b)よりも 0.532 kJ/mol だけ不安定である。そこで、これらの構造から出発し、分子骨格に与えていた束縛条件を外して同じレベルで構造最適化を行った結果を図 4 に示す。

このプロセスに伴うホルムアルデヒドの分子構造の緩和は僅かであるが(C=O 間隔で 0.001～0.004Å程度伸びている)、分子間の相対距離は大きく縮んでいる(C…H 間で 0.19～0.48 Å程度減少)ことが判る。なおこの段階では、(a)の方が(b)よりも、1.256 kJ/mol だけ不安定であるが、これらの構造を出発点として再度計算レベルを MP2/6-31G++(2d,2p)に上げて構造最適化を行うと、(a)では 2 分子間隔が接近し（C…H 間隔で、0.05～0.15 Å程度）、(b)では 2 分子間隔が広がり(C…H 間隔で 0.03 Å程度)、
結果的に(a)の方が 3.866 ｋJ/mol だけ安定となる。得られた構造を図 5 に示す。これら(a)、(b)の構造はそれぞれ、従来理論的に存在が予言されている形 4)： (a)…Cs 型、 (b)…C2h 平面型 と対応し、束縛条件をつけて超球面探索を行う本手法の有効性が確認できる。

図 3．

(a)
(b)
自動探索された 2 種類のホルムアルデヒド 2 量体
（分子内に束縛条件付き）(a)(b)
図 4．自動探索された 2 種の構造について
束縛条件を外して最適化した構造(a)(b)
図 5． 図４の構造について、更に計算レベルを
MP2/6-31G++(2d,2p)に上げて、最適化した構造【References】
1) K. Ohno and S. Maeda, Chem. Phys. Lett., 2004, 384, 277; S. Maeda and K. Ohno, J. Phys.
Chem. A, 2005, 109, 5742; K. Ohno and S. Maeda, J. Phys. Chem. A, 2006, 110, 8933
2) 山門英雄、時子山宏明、前田理、大野公一、分子科学討論会 2009、2P133; H. Tokoyama, H.
Yamakado, S. Maeda, and K. Ohno, WATOC2011 (17-22 July 2011, Sandiego de compostela,
Spain) PIII-065; Yu Sawada, Hiroaki Tokoyama, Hideo Yamakado, Satoshi Maeda, and Koichi
Ohno, 14th ICQC (25-30 June, 2012, Boulder, Colorado, USA), IV.63 他
3) 大野公一、長田有人、前田理、分子科学討論会 2010、１E15
4) Jose M. Hermida-Ramo'n and Miguel A. Rios, J. Phys Chem. A 1998, 102, 10818

3P-114鎖状構造結晶セレンの電荷移動とバンド構造: 平面波 DFT 法での孤立系 Bloch 軌道線形結合を用いた研究（京大院理） 松井 正冬
Effects of inter-chain interactions on band structure of trigonal-Se:
A DFT study using linear combination of Bloch orbitals of isolated systems(Kyoto Univ.) Masafuyu Matsui
セレン (Se) は分子性物質と金属の中間的性質を持ち、その最安定構造は鎖状構造 (trigonal
Se; t-Se) である。 t-Se では圧力に応じて金属性が増すが、それは非結合性軌道 Lone Pair (LP)から反結合性軌道 σ* への鎖間の電荷移動の増大に起因すると説明されてきた。このような機構を理論的に明らかにするには、物性科学・分子科学の両面からの理解が必要となる。過去の半経験的手法によるバンド計算の研究から、鎖状構造を金属単純立方構造からの Peierls 転移系とみなすと
1
よく理解されることが示されている 。一方で、電荷移動錯体など分子においては孤立分子軌道の線
2形結合 (LCMO) を用いて電荷移動を評価する研究があるが 、結晶系の第一原理計算では適用な
されていない。それは無限周期系である結晶でよく用いられる平面波 DFT 法では、電荷移動の記
述のための局在化軌道とその占有状態を求める方法が無かったからである。よって本研究では、結晶における電荷移動を記述するために孤立系 Bloch 軌道線型結合 (LCBO) を導入し、t-Se に適用してその電荷移動を伴う鎖間相互作用がバンド構造に与える影響を明らかにする。
結晶系が F (= A, B, ...) で構成される unit-cell からなるとき、その Bloch 軌道波動関数は、格I
子 r にある F のみからなる孤立系の Bloch 軌道波動関数の線形結合によって、
F

N I A ,B⋯ N bandII

i ,k r=∑ ∑ ∑ Fn , k r−r I ei k K ⋅r C Fn i, k / N I
I

rI

F

n

のように表される。ここで線形結合係数 C は、孤立系軌道間の重なり行列のユニタリー変換、及び孤立系軌道と結晶系軌道間の重なり行列から求まる。このようにフーリエ変換とユニタリー変換を用いることで結晶系軌道は孤立系軌道の線型結合の形式で表現できる。また、静電・分極・反発・電荷移動の各相互作用がバンド構造に与える影響を、LCBO を用いて簡便に評価する方法として相互作用
軌道エネルギー解析を導入する。これは孤立系軌道を基底とした重なり行列、ハミルトニアン行列の、各相互作用に対応する成分を考慮・無視して固有値を求めて、相互作用に対応した軌道エネルギーとして評価する方法である。
LCBO を用いて t-Se の電荷移動の挙動を電荷分解解析により評価したところ、ほぼ全ての電荷移動が図１に示した６本の最隣接鎖間で起きていることがわかった。この電荷移動による鎖間結合は配位結合的で、図１で△と▽に配列したそれぞれ逆の配向の鎖に対して、△には電子供与側、▽には受容側と、異方的に働いていた。ただし供与と逆供与はつり合うので電荷の偏りは生じない。また
電荷移動の主な成分は LP バンドから σ* バンドへのものであることが確かめられた。鎖間距離と電荷移動及びバンドギャップの関係を示したのが図２ (a) である。電荷移動によりバンドギャップが閉じ
金属性が増すことが明らかに示された。図２ (b) は逆空間内における電荷移動と直接バンドギャップ
1 T. Oka and H. Fukutome, Prog. Theor. Phys. 79 (1988) 608
2 H. Fujimoto, S. Kato, S. Yamabe, and K. Fukui, J. Chem. Phys. 60, 572 (1974)

の関係を示したものである。電荷移動が極小化する位相が結晶の直接バンドギャップを与えている。
相互作用軌道エネルギー解析により、これは図３で示したように模式的に理解できるとわかった。すな
わち t-Se では静電、分極の影響は無視できるほど小さく、鎖間の反発相互作用によりバンド幅がひろがりバンドギャップが閉じる。そして鎖間の電荷移動により、電荷移動が大きい位相はバンドギャッ
プが開き結晶の安定化に寄与し、電荷移動の小さい位相が結晶のバンドギャップを与える。このように逆空間内の電荷移動の挙動が結晶のバンドギャップを決めるが、t-Se ではこの挙動は、HOMO,
LUMO バンドの位相と結晶の三回螺旋対称性によって決まる孤立鎖軌道間の鎖間の重なりに基いてよく理解できる。
図１. 鎖軸方向から見た t-Se の構造。矢印は電荷移動の方向を示す。△と▽で示した鎖は中央の鎖に対して逆の配向を持つ。unit-cell あたりの総電荷移動量は 0.414 でうち最隣接鎖間が 0.409 (99%)。LP から σ* への成分は0.347 (84%)。１本鎖あたりの電荷移動の成分は、供与 =0.063 (93%), 逆供与 = 0.005 (7%)

(b)

(a)

図２. (a) 鎖間距離に対するバンドギャップと直接バンドギャップ（左目盛）、電荷移動（右目盛）(b) 逆空間内での直接バンドギャップ（左目盛）と電荷移動（右目盛）Isolated systems

Repulsion

Charge transfer

CT 大:● バンドギャップ開く● 結晶安定化CT 小:● 結晶のバンドギャップ

図３. 電荷移動がバンドギャップを決める模式的説明。孤立系のバンドギャップは大きい。反発相互作用によりバンド幅が拡がりバンドギャップは閉じる。電荷移動の大きい位相ではバンドギャップが開き結晶の安定化に寄与する。バンドギャップが小さい位相が結晶のバンドギャップを与える。

3P115
量子化学に基づくⅢ－Ⅴ族六方晶の安定性と反応性(横浜国大院・工) 佐々木正和，佐藤浩太
Quantum chemical study on the structure and reactivity of Ⅲ-Ⅴ hexagonalcompounds(Yokohama National University) Msakazu Sasaki, Kota Sato
【序】
Ⅲ－Ⅴ族化合物である窒化ホウ素(BN)は六方晶平面構造をとることができ、その物理的・化学的性質から様々な用途で使用される。しかし、同じⅢ－Ⅴ族化合物であるリン化ホウ素(BP)はいまだ六方晶平面構造は発見されていない。今回の研究では六方晶平面構造の BN 及びBP の電子構造について密度汎関数法(B3LYP)による構造最適化の結果を比較研究し、さらに水素ラジカル、BH2、NH2、PH2 の各反応種を求めた安定構造の各サイトに反応させ、その反応性の違い及び構造への影響を密度汎関数法により解明することを目的とした。
【方法】
プログラム：Gaussian09 近似法：B3LYP 基底関数：6－311G(d)この計算条件において B12N12H12 及び B12P12H12 をそれぞれ BN、BP モデルとして計算した。
まず、それぞれのモデルで平面性を厳密に保ったものと平面性を大きく崩したものを初期構造とし、構造最適化計算、振動解析計算を行った。次に得られた安定構造の各サイトに水素ラジカル、BH2、NH2、PH2 を条件を変えながら反応させたモデルの構造最適化を行った。
【結果及び考察】BN においては平面性を厳密に保ったものと崩したもの
ともに平面構造に収束した。振動数計算の結果からも実数の振動数のみをもつことから安定構造であると考えられ、Fig.1 六方晶窒化ホウ素(hBN)安定構造
これは実際の六方晶窒化ホウ素の構造に一致する。
BP では平面性を厳密に保ったモデルでは平面構造に収束したが、虚の振動数をもつため安定構造ではないことが分かった。平面性を大きく崩したモデルでは、平面構造に向 Fig.2 六方晶リン化ホウ素平面構造
かって計算が進行したが、収束した構造は完全な平面構造
から少し崩れたものであった。
非 平 面 構 造 の BP の エ ネ ル ギ ー は 平 面 構 造 の も の よ り
0.161eV 低くなっており、振動数も実数のもののみ

Fig.3 六方晶リン化ホウ素非平面構造となった。また、HOMO－LUMO 差を比較すると平面構造は 3.264eV で非平面構造の方は3.315eV と非平面構造の方が少し大きくなっており、最大ハードネスの原理からも非平面構造の方が安定構造であると考えられる。
得られた安定構造に水素ラジカルを反応させた。中心および外周の水素終端させた部分とそれ以外の部分それぞれのホウ素原子、
窒素原子、リン原子合計 6 ヵ所で調べた。

Fig.4 BN と水素ラジカルとの反応

BN は 6 ヵ所のうち、5 か所でエネル
ギーが増大し、熱力学的に不安定になった。これは実際の BN が化学的に安定なことに一致する。
BP は 6 ヵ所すべてでエネルギーが減
少し、熱力学的に有利に反応が進んだ。

Fig.5 BP と水素ラジカルとの反応

また、BN が反応した箇所のみ平面構造が崩れたのに対して、BP は分子全体の平面性が崩れていった。
BN は水素ラジカルとの反応性が低かったが、BP は高かったため、さらに水素ラジカルが反応した場合について計算すると Fig6.のようにさらに大きく平面性が崩れた構造となった。このとき ΔE＝－3.170eV となり、かなり熱力学的に有利な反応であると判断できる。
水素ラジカルと同様の反応位置の条件でBN に対しては BH2 と NH2 を反応させて、

Fig.6

BP と複数の水素ラジカルとの反応

BP には BH2 と PH2 をそれぞれ反応させた。
どちらも B－N、B－P 結合をつくるように反応が進みやすく、また水素ラジカルが面に対して垂直になるように結合したのに対してこれらの反応種は反応位置の隣接原子とも
population が存在し、面に対して斜めに結合した。
また、さらに複数反応した場合についても計算すると、反応した BH2 と NH2 または BH2 とPH2 同士が結合をつくる場合があった。ただし、Fig.7 BN と BH2 の反応面が平面ではない BP は、BN ほど反応種同士が結合するモデルは多くなかった。このため BNの方が六方晶の積層成長がしやすいと考えられる。
【結言】
六方晶平面構造の BP は BN ほど完全な平面構造ではないが存在することは可能であろう。
その方法は同じように六方晶構造で平面性が崩れることがヘキサシラベンゼンの計算によって報告されている 1)ケイ素分子でも格子定数をあわせた基板上に単層膜を形成できる 2)ことか
ら同じような可能性があるだろう。ただし、その反応性は高く、BN が反応位置のみでの局所的な平面性の崩れであるのに対し、反応種と結合すると分子全体で崩れていってしまう。そのため成長させるためには成長反応の際の気相中の条件を厳しく選択する必要があると考えられる。
【参考文献】

1)Shigeru Nagase, Hiroyuki Teramae, Takako kudo J.Chem.Phys.,86(1987) 4513
2)Boubekeur Lalmi, Hamid Oughaddou, et al Applied Physics Letters 97 (2010) 223109

3P-116
アセチルアセトン法によるルチジン誘導体の生成反応に関する理論的研究(1 城西大理、2NTT 環境エネルギー研) 寺前裕之 1*、丸尾容子 2【はじめに】近年シックハウス症候群と呼ばれる、住宅内装材などの化学物質が原因で起こると考えられている目や喉の痛み、あるいはアトピー性皮膚炎などが社会的な問題になってきている。ホルムアルデヒドはこのシックハウス症候群の主な原因物質と考えられている。またホルムアルデヒドは発がん性や変異原性を持ち、人間の健康に多大な影響を及ぼす。
ホルムアルデヒドは建築材料、壁紙、塗装材料、家庭用品などに広く使用されており、WHO
では 30 分での被爆量の安全値として 0.08ppm を設定するなど、低減のための対策が進められつつある。問題はどのように室内でのホルムアルデヒドの濃度測定をおこなうかである。
ホルムアルデヒドの測定には通常は溶液中ではアセチルアセトン法が用いられる。これはアセチルアセトン（β ジケトン）２分子とアンモニアイオンがホルムアルデヒドと反応してルチジン誘導体が生成する反応を利用する。
このルチジン誘導体は黄色に呈色し 410nm 付近に吸収極大を持つ。この波長の吸収強度を測定することでホルムアルデヒドの濃度を決定する。しかし、アセチルアセトン法は加熱を必要とし、溶液中では有用であるが気相での測定には不向きである。
近年、丸尾らは pentane-2,4-dione (R1,R2 = CH3)およびその置換体、1-phenyl-1,3- butanedione(R1=CH3,R2 =Ph)と 1,3-diphenyl-1-3-propanedione (R1=Ph,R2 =Ph)の３種類の β ジケトン類とアンモニウム塩を多孔質ガラス中に存在させることにより、気相での測定に使用できることを示した。
水溶液中ではルチジン誘導体の光吸収強度が時間経過で減衰してくるが、多孔質ガラス中では減衰しない(1-phenyl-1,3-butanedione ;フェニル体)。1,3-diphenyl-1,3-propanedione（ジフェニル体）は水溶液中では反応しないが、多孔質ガラス中では反応し、また HCHO の量を増やすと、強度が減衰するといった興味深い性質がいくつか明らかになってきている。
我々は以前の研究で、これらの性質を明らかにする第一歩として、ルチジン誘導体の生成反応について ab initio 分子軌道法を用いて反応機構を検討したが、多孔質ガラス中の反応については検討していなかった[1]。
本研究では ab initio 分子軌道法を用いて多孔質ガラス中での特異的な反応性を解明することを目指し、ルチジン誘導体生成スキーム中で重要と考えられる、FLUORAL-P 生成反応についてガラスのモデルとしての H2Si=O や(OH)2Si=O が存在した時の反応がどのように変化す
るかについて調べたので報告する。

【計算方法】分子軌道計算には Gaussian09 プログラムを使用した。基底関数には 6-31G**基底を使用し HF および MP2 法によりエネルギー勾配法を用いて、構造最適化を行った。求まった最適化構造は振動数計算により安定点および遷移状態であることを確かめた。
【結果と考察】FLUORAL-P 生成の素反応について図 1 に示す。

図 1 FLUORAL-P(5)生成の素反応

この反応の活性化障壁は全体の反応中で 2 番目に高く、また TS は反応課程中で最もエネルギーが高くなる。真空中での活性化障壁は MP2/6-31G(d,p)+ZPC レベルで 47.7 kcal/mol であり、水脱離反
応としては値が大きい。ただし、溶液中でのみ反応が進行するという事実と矛盾しない。
ここで、多孔質ガラスのモデルとして、H2Si=O と(OH)2Si=O を、水溶媒のモデルとして水分子を考えて、FLUORAL-P との複合体を計算した。非常に興味深いことに、図 2 に示したよ
うに、これらのモデルの Si 原子は FLUORAL-P の酸素原子と結合を作って 5 配位となってしまう。一方、炭素原子類似の H2C=O では結合は作らない。これは一見奇妙ではあるが、CO2は結晶を作らず、SiO2 は結晶となることと関連していると考えられる。

図 2 FLUORAL-P と(a)H2Si=O、(b)(OH)2Si=O、(c)水、(d) H2C=O との複合体の構造
またさらに興味深いことに、この Si 原子との間に結合を作ることにより、図 1 に示した反応の活性化エネルギーが表 1 に示したように、低くなることがわかった。
Table 1. Summary of Environmetal Effect
HF/3-21G

HF/6-31G(d,p)

MP2/6-31G(d,p)

Gas phase

50.5

58.4

47.5

H2O

51.6

58.0

48.9

H2O catalyzed

38.6

60.9

42.3

H2SiO

42.9

54.8

40.0

(OH)2SiO

39.2

52.5

36.9

参考文献
[1] Teramae, H.; Maruo, Y, Y., Intern. J. Quantum Chem., doi: 10.1002/qua.24105

3P117
ベンゾポルフィリン/C60 誘導体界面における電荷分離過程の配向依存性（東大院工）

永野智也，藤井幹也，城野亮太，山下晃一

Orientation dependence of charge-separation process at the interface
tetrabenzoporphyrin/C60 derivatives(The University of Tokyo)
Tomoya Nagano, Mikiya Fujii, Ryota Jono, Koichi Yamashita
１．序論
近年、エネルギー問題が重視される中で半永久的に利用できる太陽エネルギーを利用した太陽電池が盛んに研究開発されている。特に原料が安価で製造が容易とされている有機薄膜太陽電池が２１世紀に入ってから盛んに研究開発されるようになり、今後も新規物質開発の可能性が注目されている。
有機薄膜太陽電池において高い変換効率を達成する大きな要因の一つとして、Donor と Acceptor の界面において生成した励起子を再結合させずに電荷分離させることがあげられ、この電荷分離速度を向上させるために Donor と Acceptor の界面を制御する研究が様々になされている。そして新しい分子設計による自己組織化制御や界面での分子配向制御により、今後更に電荷分離速度が向上すると考えられる。しかしこの電荷分離機構の特性はまだ十分解明されていないため、電荷分離に適した界面の様
子はまだ分かっておらず、理論的なアプローチが期待されている。そこで本研究では、有機薄膜太陽電池における電荷分離機構の特性を解明し、高い変換効率達成のための要因を探ることを目的とする。
具体的に Donor として tetrabenzoporphyrin（BP）、Acceptor として最近の研究で新しいフラーレン
誘導体として開発された bis(dimetylphenylsilylmethyl) [60] fullerenes（SIMEF）を用いた系について、電荷分離に適した界面配向の示唆を量子化学計算を用いて行う[1]。界面配向のモデルとして、界面における二分子系に注目し、この二分子の重心間距離と相対角度という分子配向の変化に伴う電荷分離速度の変化を評価する。

２．計算手法量子化学計算のパッケージには Gaussian09 を用いた。分子のエネルギー計算、構造最適化には密度汎関数法（DFT）を用
い、汎関数 B3LYP、基底関数 6-31G(d)、6-31+G(d)を用いた。
励起状態計算には TD-DFT を用い、汎関数 B3LYP、基底関数
6-31+G(d)を用いた。遷移双極子モーメント並びに双極子モーメントの計算には CIS を用い、基底関数 6-31G(d)を用いた。
電荷分離を評価する二分子の配向には、二分子系を様々な初期構造から構造最適化したものの中で最安定だった構造及びその重心間距離、相対角度を変化させたものを用いた（図１）。
電荷分離速度の計算方法としては Marcus 理論の中間温度域

図１：分子配向のモデル図（左：重心間距離、右：相対角度）（分子内振動モードのフォノンエネルギーが熱揺らぎより十分小さく、溶媒振動モードのフォノンエネルギーが熱揺らぎより十分大きい温度域）の式（１）を用いた。ここで Vif は始状態（励起状態）と終状態（電荷分離状態）の電子カップリング、λ s は媒体の再配置エネルギー、T は温度（300K を用いた）、Sq は核のスペクトル関数、ΔG0 は始状態と終状態のエネルギーギャップ、ω q は核の平均角周波数である。物理量の算出法は以下の論文[2][3]を参考にし、電子カップリングの計算には Generalized
Mulliken Hush 法を用いて書き下した式（２）を用いた。ここでμ if は始状態から終状態への遷移双極子モーメント、Δμ if は始状態と終状態の双極子モーメントの差、ΔEif は始状態と終状態のエネルギー差である。
∞
 (∆G 0 + λ s + nℏω q )2 
S qn
2π
1
2
k=
Vif

∑ exp(−Sq ) n! exp−ℏ
4π λ s Τ n =04λ s k B T



µ if

Vif =

(∆µ if )2 + 4(µ if )2

（１）

（２）

∆E if

３．結果と考察
計算によって得られた BP の UV-VIS スペクトルの 590nm、580nm の光吸収ピークについて、電荷分離速度の重心間距離、相対角度依存性を算出した結果を示す（図２）。電荷分離速度は各励起状態によって配向依存性が異なり、各励起状態毎に電荷分離に適した構造が存在することが示唆された。またこのような配向依存性を示す原因として、本系では電荷分離速度の項の中で電子カップリングの項の遷移双極子モーメントが最も配向依存性が強く支配的であることが分かった。つまり二分子系の最安定構造が電荷分離に最も適した配向ではなく、分子配向制御することによって更なる電荷分離速度の向上が期待できると考えられる。特に、今回の結果では相対角度については一概に言えないが、重心間距離は安定構造の重心間距離 7.0Åより短くすることで電荷分離速度を数倍向上させられると考

14

14

13

13

590[nm]580[nm]12
11

log(k) [1/s]log(k) [1/s]えられる。計算手法、結果の詳細は当日発表する。

12

590[nm]580[nm]

11
10

10
6

6.5

7

7.5

重心間距離[Å]

9
5

10

15
20
相対角度[°]

25

図２：各励起状態からの電荷分離速度の配向依存性

参考文献[1] Matsuo, Y., Sato, Y., Niinomi, T., Soga, I., Tanaka, H., and Nakamura, E., J. Am. Chem. Soc.
131, 16048 (2009)[2] Bredas, J., Beljonne, D., Coropcoanu, V., and Cornil, J., Chem. Rev. 104, 4971 (2004)[3] Yuanping, Y., Coropceanu, V., Bredas, J., J. Am. Chem. Soc. 131, 15777 (2009)

3P118カーボンナノチューブ/ポリチオフェン
界面構造の理論的研究(東大院・工) 西村 亮彦、藤井 幹也、 城野 亮太、 山下 晃一

Theoretical analysis on carbon nanotube - polythiopheneinterfaces(Univ. of Tokyo, Eng.) Nishimra Katsuhiko, Fujii Mikiya,
Jono Ryota, Yamashita Koichi

序論
近年、有機薄膜太陽電池はその大量生産可能性や製造時のエネルギー消費が比較的少な
いこと、また素材の柔軟性などから、次世代の太陽電池として開発が進められている。その中でも、P3HT/PCBM バルクヘテロ接合太陽電池は比較的確立した成果が得られており、様々な研究がこれまでになされている 1 。PCBM と同じナノ炭素材料であるカーボン
ナノチューブは、PCBM と比べてキャリア移動度が高い 2,3 ことや、その高い比表面積な
どから、バルクヘテロ接合太陽電池材料として PCBM よりも適していると考えられる。
そのため、P3HT/SWNT バルクヘテロ接合太陽電池の実現可能性について多くの研究が成されてきた。しかし、その多くは P3HT/SWNT 接合太陽電池は P3HT/ PCBM 接合太陽電池と比べて著しくキャリア寿命が短く、太陽電池として実用できるレベルではないという結果に終っている。
このような状況の中で近年、これまでの結果に比して非常に長寿命で、P3HT/ PCBM
接合と比較できる程度のキャリア寿命が P3HT/SWNT 接合において得られたとする結果が報告された 4 。このことは P3HT/SWNT 接合面における P3HT がバルクの P3HT とは異なる分子構造・電子状態にあることが原因であると考えられている。
本研究の目的は、P3HT/SWNT 接合面における P3HT の分子構造・電子状態を調査し、より高効率の有機薄膜太陽電池の開発へとつなげることである。

方法
本研究ではモデル系として数百原子を含むポリチオフェンと SWNT を用意し、半経験的分子軌道法を用いて解析を行なった。

P3HT/SWNT 接合面においては P3HT は π − π スタックによって SWNT と分子間結合
し、SWNT に螺旋状に巻きついた構造をとることが知られている 5 。π − π スタックには

P3HT の側鎖は寄与しないと考え、P3HT の側鎖を除いたポリチオフェンと SWNT について、様々な初期構造から構造最適化を行なった。また、π − π スタックの引力は分散力に相互作用に起因することが知られているため、PM3 に分散力補正を加えた PM3-D 6,7を用いて構造最適化を行なった。
また、ZINDO/S による励起状態計算により、分子構造に依存して光学的性質がどのように変化するかについても当日報告する。

結果

図 1: ポリチオフェン/SWNT の安定構
造 (赤:ポリチオフェン骨格炭素)

図 2: P3HT/SWNT の螺旋構造 (赤:ポリチオフェン骨格炭素 緑:側鎖炭素)

ポリチオフェン/SWNT 系について様々な螺旋状初期構造から PM3-D により構造最適化を行なった結果、全ての初期構造について図 1 のようにポリチオフェンが直線状構造に収束した。また、エネルギー面はポリチオフェンが SWNT 表面上を滑る方向には非常に平坦であることがわかった。このことから側鎖の存在が接合面における分子構造の決定に大きく寄与しているものと思われた。そこで側鎖を含む P3HT/SWNT 系の構造最適化を行なったところ、図 2 のような結果が得られ、構造の決定において側鎖の重要性が示唆された。
また、図 1 の安定構造ではポリチオフェンの共平面性が崩れており、共平面性の崩れた構造は π 共役系の切断を引きおこし、キャリア移動度を低下させ再結合を起こりやすくすると考えられるので、側鎖の存在がポリチオフェン誘導体/SWNT 接合太陽電池を作成する上で重要であると思われる。

参考文献[4] Stranks, S. D. et al. Nano Lett. 2011,
11, 66–72.

[1] Hauch, J. A. et al. Sol. Energ. Mat.
[5] Goh, R. G. S. et al. Appl. Phys. Lett.
Sol. C. 2008, 92, 727 – 731.
2006, 88, 053101.
[2] Fuhrer, M. S. et al. Nano Lett. 2002, [6] McNamara, J. P. et al. Phys. Chem.
2, 755–759.
Chem. Phys. 2007, 9, 2362–2370.
[3] Mihailetchi, V. et al. Adv. Funct. [7] Morgado, C. A. et al. J. Chem. Theory
Mater. 2003, 13, 43–46.
Comput. 2007, 3, 1656–1664.

3P-119

異分子間電子移動の解析
筑波大・化

鞆津

典夫、相川

小春、守橋

健二

Theoretical investigation on hole transfer between different molecules(Univ. of Tsukuba)

Norio Tomotsu, Koharu Aikawa, Kenji Morihashi[序]有機 EL は数種の分子層を積層して陰極および陽極に電圧をかけることにより電子と正孔がそれぞれの電子輸送層・正孔輸送層を通過し、発光層で結合する。層内の同一分子間の電子移動および層間での異分子間の電子移動が有機 EL の高効率化には重要な因子となる。二分子間の電子移動においては結晶構造が明確な分子において多くの研究例があるが、結晶をとらない分子や、異分子間の電子移動についてはいまだ研究が行われていない。
本研究では、Wu ら[1]が開発した制約密度汎関数理論(Constrained

DFT（以下 C-DFT)）を用いて Alq3 または正孔輸送材（以下 HT、具体的には BCP,CBP,NPD）の一方に正孔を置いた 2 分子の計算を行い、Marcus の理論に基づいた正孔移動速度を算出し、考察した。
[計算]
GAUSSIAN09 を用い、B3LYP/6-31G*で fac 体 tris(8-hydroxyquinolinato)aluminium (以下fAlq3), mer 体 tris(8-hydroxyquinolinato)aluminium (以下 mAlq3), BCP, CBP,α-NPD,の各分子の中性分子とカチオン分子の安定構造を求めた。Fig-１に BCP, CBP,α-NPD,の中性分子(黒)とカチオン分子(青)の構造を示す。

BCP

αNPD

CBP

Fig-1 正孔輸送材の分子構造、青色はカチオン構造
2 分子の配座は中性分子またはカチオン＋

BCP

分子の構造を固定し、いくつかの配座を初期構造とし、最安定構造を求めた。Fig-2に mAlq3 と BCP の構造を示す。黒の分子構造は mAlq3+と CBP（始状態）、青の分子
構造は mAlq3 と CBP+、(終状態)の安定構

Alq

3

BCP

造である。CBP, α-NPD と mAlq3 の構造も同様にして求めた。

Fig-１

mAlq3+と BCP の始状態(黒）と終状態(青)更に構造異性体である fAlq3 についても CBP, α-NPD, BCP との構造を最適化し、得られた構造の
C-DFT 計算を当研究室が開発した DFT 計算ソフトを用い[2]、B3LYP/6-31G で行い、電子状態から

Marcus 理論を用いて速度定数を求めた。
[結果と考察]計算で得られた Marcus パラメータと速度定数を Table-１にまとめた。Alq3+ + HT → Alq3 + HT+の反応を正方向、その反対を逆反応とした。
Table-1正方向の移動速度

エネルギー差 kJ/mol
活性化エネルギー kJ/mol
正移動の速度定数 s-1
逆移動の速度定数 s-1

ｍAlq3-BCP ｍAlq3-CBP ｍAlq3-NPD ｆAlq3-BCP ｆAlq3-CBP ｆAlq3-NPD
23
9
58
24
50
120
34
10
137
29
186
6086
1.91E+08
9.72E+12
1.73E-12
2.16E+09
1.20E-18
0.00E+00
3.32E+12
2.68E+09
4.17E-41
5.23E+13
1.90E-31
0.00E+00

BCP は正逆ともに移動速度が大きく、特に逆反応の方が速いという結果を得た。Fig-3 に fAlq3 と
BCP との電子移動の際の Corresponding Orbitals を示した。

Fig-3

fAlq3 と BCP の正孔輸送における移動前後の Corresponding Orbitals

CBP も同様な傾向を示しているが、正逆反応の速度定数の大きさは BCP の場合と反対になっている。
更に fAlq3 に対しては移動速度が低下している。Fig-4 に fAlq3 と CBP の正孔移動前後の構造とCorresponding Orbitals を示した。mAlq3 と fAlq3 では配座が異なり、Alq3 に CBP が接近して安定構造を取る際に fAlq3 は mAlq3 に比べ、接近しにくいために、軌道の相互作用が起こりにくくなっていると考えられる。

Fig-4

fAlq3 と CBP の正孔輸送における移動前後の Corresponding Orbitals

これに対し生成系のエネルギーが反応系のエネルギーよりかなり大きな NPD は正方向も逆方向も活性化エネルギーが高く、正孔の移動はほとんど起こらないと考えられる。
[参考文献][1] Q. Wu, T. Van Voorhis, Phys. Rev. A 72, 024502 (2005).
[2] T. Ogawa, M. Sumita, Y. Shimodo , K.Morihashi, Chem. Phys. Lett.,

511, 219 (2011).

3P-120

カーボンアロイ触媒の電⼦状態と触媒活性に及ぼす置換基効果(東⼤院・⼯) 畑 智⾏, ⽜⼭ 浩, ⼭下 晃⼀

CACs electronic state of carbon alloy catalyst andsubstituent effects on the catalytic activity(University of Tokyo, Eng) Hata Tomoyuki, Ushiyama Hiroshi, Yamashita Koichi
1 緒言
燃料電池とは、水素が酸素によって酸化される際に放出される自由エネルギーを、直接電気エネルギーに変換するデバイスである。近年では非白金系正極触媒の研究が行われており、正極においては、ZrTa酸化物や錯体添加フロン樹脂、炭素素材等が触媒活性を示すことが発見されている。今回は、炭素素材のなかでもグラフェンナノリボンを層状に配置した触媒構造であるカーボンアロイ
触媒(CACs) に注目した。
CACs は、酸素のエッジ上への吸着により、酸素還元触媒として働くことが実験から予想されており、代替正極材料として注目されている。正極での酸素還元反応は、CACs への吸着反応と酸素への電子移動反応の二つの素過程に分けられる。後者の電子移動反応が律速であり、よってCACsと吸着酸素のエネルギー準位を比較することで、触媒活性を定量的に評価できることが知られている。
また、第一原理分子動力学計算により、CACs の表面を窒素等の他元素で置換することで、触媒活性の上昇が期待される結果が得られている[1][2]。しかし、そもそも酸素がCACs 上に吸着するか否かに関しては、吸着機構が不明瞭なため明確な基準がなく、活性の確認はできても、設計への指針が示されていないのが現状である。以上より、CACs の触媒設計において、酸素還元反応の活性を総括的に評価し設計指針を開発するためには、酸素の吸着機構を知ることが極めて重要である。以上を背景として本研究では、多数の窒素置換CACs触媒モデルに対して第一原理電子状態計算を行い結果を比較することで、CACs の酸素吸着機構及び窒素置換の効果についてより深く理解することを目的とする。具体的には、窒素の置換箇所が吸着機構に与える影響を化学的に評価することを目指した。
2 方法
今回は、以下のようなCACs モデルを、置換数を変えて計9 種類作成した(図1)。以降はモデルをmodel[置換数-置換箇所] と呼称し、置換箇所の指定には図1のように定義したラインを用いた。
例えば図2 の場合は窒素を第2ラインに4置換しているため、model[4D-2L] と呼称する。置換していないCACs はmodel0 とした。先行研究から、第1 ライン置換構造は酸素が吸着しないことが分かっている。CACs 触媒モデルの計算には、密度汎関数法(DFT)、汎関数:にはGGA-RPBE、基底関数には
DZP 局在基底を採用し、計算パッケージ:にはsiesta-2.0.1を用いた。すべてのモデルについて吸着媒
単独での状態密度を比較した結果、先行研究における触媒活性の変化に伴うピークシフトが確認できたものの、吸着機構については多くの情報は得られなかった。そのため、CACs モデルに実際に酸素を吸着させ、その状態密度を計算し比較した(図2)。先行研究で発見された素過程の通り[1]、エッジ端

の炭素に酸素が吸着することで反応が進行すると仮定している。また状態密度のピークに対応する軌道を調べ、酸素の分子軌道と吸着媒との相互作用を調べることで、吸着機構とその変化を説明した。

(図1 : 置換モデル)(図2 : 酸素吸着モデル)

3 結果
モデルに対し酸素を吸着させたところ、先行研究と同様に第1ライン置換構造は酸素が吸着しなかったのに対して、第2ライン置換構造には酸素が吸着した(図2)。また、model0に対しても酸素は吸着した。計算結果の考察には、置換の影響が顕著だろうとの考えから、まず4置換モデルを用いた(図3)。酸素分子単体とのエネルギーを比較することにより、CACsのFermi エネルギー近傍の軌道と相互作用している軌道は、酸素の1πg*軌道であると分かった。状態密度のピークに対応する軌道を見ていくと、Fermi エネルギー近傍においてエッジ上の水素と炭素間の反結合性軌道や、吸着スポット
上に状態がない非結合性軌道をみることができた。また両モデルの状態密度について、酸素の軌道に相当するピークにおける電子の占有状態を調べると、吸着しなかったmodel[4D-1L]では水素炭素間の反結合性軌道(図4)を2電子が占有している一方で、吸着したmodel[4D-2L]では反結合性軌道と非結合性軌道とがFermiエネルギー上で縮退し、電子が分数占有することで、構造が安定化している
ことが分かった。以上のような酸素分子とCACsの軌道相互作用を、軌道相関図を用いて解析した。

(図3 : 状態密度図)(図4 : 反結合性軌道)[1] T. Ikeda, M.Boero, S.-F. Huang, et al.: J. Phys. Chem. C, 2011.
[2] S.-F. Huang, K. Terakura, T. Ozaki, et al.: Phys. Rev. B, 2009.

3P-121
DFT 計算を用いたニトリルの水和反応における CeO2触媒の基質特異性の検討(名古屋大学大学院・工) 吉川幸男, 沢邊恭一, 薩摩篤Theoretical study on substrate specificity of a CeO2 catalyst for nitrile hydration using DFT
calculation(Nagoya Univ.) Yukio Yoshikawa, Kyoichi Sawabe, Atsushi Satsuma【序】 塩基性触媒を用いたニトリルのアミド化は、OHの CN 基への求核攻撃による水和反応
で開始する。水溶液中で OHを供給できる CeO2 触媒はニトリルの水和反応に高活性であることが報告された [1]。NaOH 水溶液を用いたニトリル水和反応では、α位にヘテロ原子を持つ 2-シアノピリジン(以下 2CNPy)よりもγ位にヘテロ原子を持つ 4-シアノピリジン(同 4CNPy)のほうが高収率であり、これは各基質分子の LUMO による求核試薬反応性の予想と一致する。しかし、CeO2触媒を用いた反応では、Table 1 に示すように 2CNPy の収率が 99%となり、3CNPy と4CNpy はほとんど反応しない。さらに、他の基質でも(X:ヘテロ原子)の構(a)2CNPy(b)3CNPy(c)4CNPy＞99(30℃,12h)

5(100℃,24h)

4(100℃,24h)

Nitrile造を持つ基質に対し、高収率でアミド化が進行するという基質特異性が CeO2 触媒で発現

Yield
[%]

した[1]。この反応の律速段階は CeO2 表面で生成した OHの求核攻撃であることが示唆されている。したがって、基質の CN 基が表面 OH 基と隣接する吸着構造をとると、反応

Table 1 CeO2 触媒を用いたアミド化の基質適応性が進行しやすいと考えられる。我々は、基質特異性を発現する基質分子は、ニトリルの N とヘテ
ロ原子の両者が表面と二点結合する吸着構造をとることで、CN 基が触媒表面と接近可能になると予想した。そこで、反応活性と予測されている CeO2(110)表面における 2CNPy と 4CNPy の吸着構造を DFT 計算によって研究し、CeO2 触媒の基質特異性について検討した。
【計算方法】 DFT 計算には 30Ry のカットオフエネルギーを持つ平面波基底を用い、擬ポテン
シャルとして Vanderbilt の ultrasoft 型のものを採用した。CeO2 表面には 2×2 の大きさの (110)
面を 6 層または 3 層重ねたスラブ模型を用い、真空層の厚さは 14 Åとした。汎関数は PBE を使用し、1.5eV の Hubbard U 値で DFT+U 計算を行った。サンプル k 点は Monkhorst-Pack による 2×2×1 を使用した。吸着エネルギーは Ead = E(孤立系)－E(吸着系)と定義し、正の値が安定な吸着を示す。計算プログラムは pwscf を用いた。構造最適化計算では、表面 2 層と吸着分子の座標を自由にして最適化計算した。
【結果と考察】 CeO2 表面にはルイス酸点となる Ce 原子が存在するので、この酸点と基質の塩基点である N 原子が相互作用して吸着すると考えられる。そこで、最初に欠陥のない CeO2 表面と 2CNpy および 4CNpy の吸着構造を調べることで、2CNpy が二点のヘテロ原子を介して吸着するかどうかを検討した。6 層のスラブ模型上での吸着構造と吸着エネルギーを Fig. 1 に示す。
1

M. Tamura, H. Wakasugi, K. Shimizu, A. Satsuma, Chem. Eur. J. , 17(2011), 11428.

4CNpy ではニトリルまたはピリジン環の N のどちらかで吸着し、分子が直立した構造となる。
そ れに対し て、
2CNPy ではニトリルとピリジン環の二つの N が表面 Ce 原子と結合する。この構造では CN 基が表面に対して斜め横の配向となるが、N 原子
二点と表面 Ce 原子間Fig.1 4CNPy(a) (b)と 2CNPy(c)の吸着構造および吸着エネルギー(eV)．
(d)は 2CNPy 吸着の電子密度等高面(0.02 a.u.-3) ．

で共有結合していることは Fig.1(d)の電子密度で確認された。なお、構造(c)の吸着エネルギーは構造(a)と同じであり、二点で吸着しているにも関わらず吸着結合は強くなっていない。これは、二点結合をすると各 N 原子とルイス酸点の間の距離(Fig.1(c) 2.93Å,2.89Å)が、一点結合で吸着する構造(Fig.1(a) 2.73Å)より離れてしまうためである。
次に表面 OH 基が存在する CeO2(110)表面の酸点の強さについて考察する。なお、以下はすべて 3 層のスラブ模型を用いた計算結果である。CeO2 触媒の水和反応では、酸素欠陥サイトで水が解離して生成した表面 OH 基が水和反応に関与すると示唆されている[1]。
そこで、欠陥サイト上で水が解離した構造を最適化計算したところ、Fig.2 の結果を得た。表面 OH 基は表面に対する配向が異なる 2 種類の構造が得られた。したがって、これらの OH 基に隣接する Ce
原子のルイス酸点も二種類存在する。垂直な OH 基側の Ce 酸点を
A サイト、横向きの OH 基側の Ce 酸点を B サイトとする(Fig. 2)。

Fig.2 酸素欠陥上で解離吸着した H2O の最適化構造

これらの酸点と OH 基のない表面の酸点の強さを比較するため、塩基分子として HCN 分子を用いて吸着エネルギーを計算した(Table 2)。
その結果、OH 基が存在する表面の方が A、B サイトでそれぞれ 0.03、0.05 eV 吸着エネルギーが小さくなった。したがって、OH 基の存在下では、ルイス酸強度はやや弱くなっている。しかし、2CNpy の吸着
エネルギーは、OH 基が存在する表面の方が 0.13 eV 大きい。各酸点の強さは弱いが、CeO2 表面

吸着分子

with OH
adsorbates

HCN

without OH
adsorbates

2CNpy

吸着サイト

Ead (eV)A site0.23

B site0.21

A sites

0.75

HCN

0.26

2CNpy

0.62

Table 2 CeO2 表面に 3 層スラブ模型を用いたときの吸着エネルギー

OH 基の存在下で Ce 間距離が 3.88Åから 4.02Åに伸び、2CNpy の CN 側の吸着で Ce-N-C の間の角度が 121°から 126°に広がったため、CN 側の酸－塩基相互作用が強くなり、二点吸着が強まったと考えられる。
以上の結果より、次のことが示唆される。CeO2(110)表面の [1̅10]方向に並んだ Ce 間距離は約4Åである。基質分子において、塩基点となる N 原子とヘテロ原子 X が 4Åより若干短い距離で横に並ぶ構造があると、CeO2(110)表面はこの基質分子を二点で吸着できる。そのため、CN 基は表面に接近可能となり、水和反応において基質特異性が現れる。

3P122半導体に連結した金属錯体の電子状態に関する理論的研究(豊田中研 1, JST さきがけ 2) ○白井聡一 1, 佐藤俊介 1,2, 森川健志 1, 旭良司 1shirai[at]mosk.tytlabs.co.jp (送信時は[at]を@に変更)【緒言】 地球温暖化抑制と石油に替わる炭素源確保を両立しうる手段として、光エネルギーを利用した CO2 還元が盛んに研究されている。最近、我々は酸化物半導体の表面に金属錯体を連結した新奇な CO2 還元光触媒を見出した[1]。さらに、この触媒を水酸化触媒系と組み合わせることにより CO2 と H2O からギ酸 HCOOH を生成することに成功し、初の人工光合成系として提案した(Figure 1)[2]。一方、反応効率はまだ低く、その向上が課題となっている。還元側触媒は、半導体を光励起し、生じた励起電子を表面の金属錯体に移動させ、錯体上で CO2を還元する。ゆえに、触媒の活性は半導体表面と錯体の電子状態、およびそれらの相互作用に依存すると推測されるが、その詳細は明らかになっていない。そこで、実験的解析と量子化学計算を併用して詳細を解明し、触媒活性の向上につなげたいと考えている。本研究ではその一環として、酸化チタンクラスターに Re トリカルボニ
ル錯体Re(4,4’-X2-bpy)(CO)3Cl を連結した系を CO2 還元光触媒のモデルとして計算を行い、電子状態を解析した。特に錯体と半導体を架橋するアンカー部位X に着目し、その電子的性質が錯体の電子状態に与える影響を詳細に解析したので報告する。
HCOOH

Electron Transferhν

H2O Oxidationhν

CO2 + 2H+
-

M

e

e-

Metal Complex

h+
h+

H2O

CO2 Reduction

X=H, COOH, PO3H2, OH+

1/2O2 + 2H

Figure 2. Re(4,4’-X2-bpy)(CO)3Cl.

Figure 1. Total reaction of the Z-scheme system forCO2 conversion to formic acid.
【計算の詳細】 計算法として B3LYP を汎関数とする密度汎関数法、基底関数系として LANL2DZ を用いた。ただし、P 原子には cc-pVDZ の分極関数を追加した。全ての計算は、Gaussian03 を用いて行った。
［Re 錯体］ Re(4,4’-X2-bpy)(CO)3Clのトリカルボニル(CO)3 は赤外吸収スペクトルに特徴的なピークを与え、実験的にも容易に解析できる。そこで、X=H,
COOH, PO3H2, OH の各錯体について構造最適化
および振動計算を行い、アンカーの電子的性質と(CO)3 伸縮振動数の相関を解析した。さらに、錯体の電荷を±1 として同様の計算を行い、電荷状態と振動数の相関を解析した。
［複合体モデル］ Re錯体は酸化チタン表面に対し、２つのアンカー部位で吸着した構造をとることが示

ProtonFigure 3. Re complex-TiO2 cluster hybrid model.
Re(4,4’-COOH2-bpy)(CO)3Cl is presented as an example.

唆されている[3]。一方、(TiO2)n の n=10-16 のクラスターの構造および物性に関しては、詳細な理論的研究が行われている[4]。本研究では、その中で２座吸着構造を表現可能な(TiO2)15 クラスターに X=COOH, OH の錯体を連結し、半導体表面に連結された錯体を表現した(Figure 3)。X より解離した２つのプロトンはともにクラスター上の酸素原子に吸着させ、モデルの電荷を中性とした。構造最適化、振動計算および自然電子密度解析を行い、錯体単体の性質と比較した。
【結果および考察】 Re 錯体単体に関する計算結果

Table 1. Frequencies of (CO)3 stretching inを Table 1 に示す。いずれも(CO)3 伸縮振動に特徴

Re(4,4’-X2-bpy)(CO)3Cl complex (cm-1).
X的な３つの振動モードが見られた。X=H に対し、X
Hが電子吸引性置換基となる X=COOH, PO3H2 では振動数増、電子供与性置換基となる X=OH では振動数減となった。また、全ての場合で電荷+1 では中1886.81909.8
1984.5

1958.1
1961.32041.51850.1
1866.4
1957.6

PO3H2

1886.31908.1
1984.8

1964.3
1966.62040.01846.81867.01957.4

OH

1875.7
1898.01980.7

1955.7
1962.52036.51840.51854.6
1951.3反結合性となっている(Figure 4)。ゆえに、錯体が電子不足となり、HOMO から電子が抜けると CO 結合
は強く、電子豊富となって LUMO に電子が入ると弱－1
1839.41854.01949.9COOH

性状態に対して振動数増、−1 では振動数減となった。錯体の HOMO および LUMO では、CO 結合がmolecular charge
neutral+11879.1
1962.5
1902.1
1966.4
1982.42039.6

くなると考えられる。振動数は結合が強いほど高い
ため、振動数の計算結果は軌道から予測される挙

(a)

(b)

動と整合する。すなわち、半導体に連結した錯体の電子状態は、(CO)3 伸縮振動の振動数変化より推定できる。次に、X=COOH, OH の錯体を連結した複合体モデルの計算結果を Table 2 に示す。振動数は、

Figure 4. Molecular orbitals of Rebpy(CO)3Cl(a) HOMO, (b) LUMO

X=COOH, OH ともに錯体単体の場合に近い。また、錯体部位の電荷は中性に近く、錯体･半導体間の顕著な電子的相互作用は生じていないと考えられる。

Table 2. Frequencies of (CO)3 stretching and natural

実際の半導体表面には、本モデルでは考慮していない構造のゆらぎや欠陥が存在しうる。また、錯体
以外にも水分子や酸素分子が吸着している可能性

charge ρ in Re complex - TiO2 cluster hybrid.
-1

Frequencies (cm )がある。錯体の電子状態に影響しうるこれらの因子を近似的に考慮した計算も進めており、当日はそれらの結果も合わせて報告する。

ρ[Re-complex]
ρ[(TiO2)15]

X=COOH1885.71910.6
1984.3

X=OH
1873.91900.3
1980.1+0.247-0.247+0.300-0.300【参考文献】[1] Sato, S.; Morilawa, T.; Saeki, S.; Kajino, T.; Motohiro, T. Angew. Chem. Int. Ed. 2010, 49, 5101.
[2] Sato, S.; Arai, T.; Morikawa, T.; Uemura, K.; Suzuki, T. M.; Tanaka, H.; Kajino, T. J. Am. Chem. Soc. 2011, 133, 15240.
[3] Anfuso, C.; Snoeberger, R. C.; Ricks, A. M.; Liu W.; Xiao, D.; Batista, V. S.; Lian, T. J. Am. Chem. Soc. 2011, 133, 6922.
[4] Qu, Z-. W.; Kroes, G. -J. J. Phys. Chem. C. 2007, 111, 16808.

3P123
分子動力学法を用いた水素ハイドレートの研究（金沢大院•自然）東	  真史,	
 三浦	
  伸一	
 

A molecular dynamics study on hydrogen hydrate(Kanazawa Univ.) Masashi Higashi, Shinichi Miura【序】	  クラスレートハイドレート構造と呼ばれる、水分子が水素結合によりかご状の構造をとり、その中に通常はガスとして存在する分子（ゲスト分子）を取り込むことができる構造がある。ガス分子と水分子の間には化学結合は存在せず、ガス分子と水分子との相互作用でゲスト分子を閉じ込めている。このハイドレートに水素分子が取り込まれたものが水素ハイドレートである。また、水素が氷の様々な多形に取り込まれた構造が知られている。水素ハイドレートは水素の貯蔵物質として期待されており、環境に優しいエネルギー源として次世代エネルギーの有望な候補の一つとして、メタン
ハイドレートと共に注目されている物質である[1]。	
  	
  	
  	
  	
 
	
  一方、水素ハイドレート中の水素分子は一つの空隙中に静止しておらず、非常に速い速度で氷の空隙間を動き回っていることが実験より示唆されている[2]。水素分子の拡散係数は通常の固体中のものと比べて桁違いに大きく、むしろ液体に近い値を持つ。
このことより拡散過程の活性化エネルギーが低いことが予想される。つまり、この水
素ハイドレートが他の固体の物質と比べて何か特徴的な構造特性を持つことが予期される。本研究では、水素ハイドレートの経路積分分子動力学計算（PIMD）を実施し、その性質の解析を行った。	
 	
 【計算方法】	
 	
  本研究では、水分子のモデルにはq-TIP4P/F[3]を用いた。このq-TIP4P/F	
 は、分子間相互作用に剛体モデルのTIP4P/2005[4]のパラメーターを用いた4サイトのフレキシ
ブルモデルである。水素分子はフレキシブルモデルを用いる。水素分子の分子間相互作用のパラメーターには剛体モデルに採用されているものを用いた。水分子内の相互作用はモースポテンシャルをテーラー展開し４次までの式としたもので表す。水素分子のH-H結合の相互作用にも同様に修正したモースポテンシャルを用いる[5]。また、
水分子はH-O-H結合角に関するポテンシャルには調和型の関数を用いた。分子間相互作用は水分子の酸素-酸素間、水分子の酸素-水素分子の重心間でLJポテンシャルタイプの相互作用をする。部分電荷としては、水分子は水素原子にプラスの電荷がありM	
 サイトにマイナスの電荷がある。水素分子は水素原子にプラスの電荷があり重心の位置にマイナスの電荷がある。電荷のあるサイト同士でクーロン相互作用をし、クーロ

ン相互作用の計算にはエワルド法を用いた。計算方法としては経路積分分子動力学法を用いた。離散化した経路積分表示により、分配関数は	
 となる。またこのときのWMは、である。これは、温度Tにおける各分子M個のサイトからなるN個の環状高分子の系の相互作用と等しい。各サイトは最隣接のサイトとバネで結合されており、分子間では同じラベルのサイト間のみに元の相互作用の1/M倍の相互作用をする。ここで、高分子の各サイトに仮想的な運動量と質量を導入し、運動方程式を導き、分子動力学計
算（MD）を実行する。

氷の構造には氷Ihを用いた。全水分子数を９６個とした。そして、氷の中に水素分子を一つ入れた状態で計算を行った。温度は100Kとし、温度制御には能勢-Hoover 鎖法を用いた。左図は平衡化後のスナップショットを示す。計算結果、解析結果は当日示す。

図 : 平衡化後のスナップショット
参考文献[1] W. L. Mao, C. A. Koh, and E. D. Sloan, Phys. Today 42 October 2007 .
[2] H. L. Strauss, Z.Chen, and C.-K.Loong, J. Chem. Phys, 101, 7177 (1994) .
[3] S. Habershon, T. E. Markland, and D. E. Manolopoulos, J. Chem. Phys, 131,
024501 (2009) .
[4] J. L. F. Abascal, and C. Vega, J. Chem. Phys, 123, 234505 (2005) .
[5] J. Schaefer, W. A. Lester Jr, D. Kouri, and C. A. Wells, Chem. Phys. Lett. 24(1974) .

3P124

時間依存断熱状態法のナノカーボンへの応用：
近赤外フェムト秒パルス誘起振動ダイナミクスの動力学計算（東北大院理 1）山崎 馨 1，中村 尭祉 1，新津 直幸 1，菅野 学 1，河野 裕彦 1

Application of Time-Dependent Adiabatic State Approach to Nanocarbons:
Simulation of Near-infrared Excitation of Raman Active Modes(Tohoku Univ. 1) Kaoru Yamazaki1, Takashi Nakamura1, Naoyuki Niitsu1,
Manabu Kanno1 and Hirohiko Kono1
1．序論
近年，高い強度を持ったフェムト秒 (10 −15 秒, fs) 幅の近赤外レーザーパルス（近赤外フェムト秒パルス）は，ナノカーボン製分子回路の作成に利用できると期待されている[1]．近赤外フェムト秒パルスがナノカーボンに照射されると，そのパルス形状に依存して異なる分子振動にエネルギーが注入される[2,3]．この振動エネルギーが他の振動モードに再分配されて，ピコ秒 (10 −12 秒, ps)からマイクロ秒（10 −6 秒, μs）程度の時間スケールで転位や炭素フラグ
メントの脱離などの構造変形が起きる [2,3]．このため，ナノカーボンの近赤外フェムト秒パ
ルス誘起構造変形の動力学を解明するには，少なくともパルス照射から構造変形に至るまでの数十 ps から数 ns の分子動力学（MD）計算を効率良く実行する必要がある．
そこで我々は，近赤外フェムト秒パルスにより誘起されるカーボンナノチューブなどの大型のナノカーボンの動力学を高速かつ半定量的に追跡できる時間依存断熱状態 DFTB/MD 法を開発した．この手法では，レーザー電場との相互作用により時々刻々と変形する時間依存断熱ポテンシ
ャルエネルギー曲面(時間依存断熱 PES) En(R,t) [4]が原子核に及ぼす力を，半経験的分子軌道法の１つである Self-consistent charge density functional based tight binding (SCC-DFTB)法[5]により求め，
MD 計算を行う．本発表では時間依存断熱状態 DFTB/MD 法の概略と近赤外フェムト秒パルス（列）を用いた C60 フラーレンの非共鳴ラマン振動励起ダイナミクス[3]への適用例について議論する．
2．理論：時間依存断熱状態 DFTB/MD 法
まず，時間依存断熱状態 DFTB/MD 法の基礎となる時間依存断熱状態法[4]について概説する．時刻 t，核座標 R における時間依存断熱ハミルトニアン H (R,t) が，レーザー電場が印加されていないときのハミルトニアン H (R, ε(t) = 0) 及びレーザー電場と分子の電気双極子相互作
用 V(t) = − με(t) の和で書き表されるとする．この H (R,t) = H (R, ε(t) = 0) + V (t) を各時刻において定常状態のシュレディンガー方程式 H (R,t) Ψn (R,t) = En (R,t) Ψn (R,t) に代入し，時間依存断
熱状態 n= 0, 1, 2,…の波動関数 Ψn(R,t) とレーザー電場により変形した 時間依存断熱 PES
En (R,t) を求める[4]．着目する En (R,t) が原子核に及ぼす力を求めることでレーザー誘起化学
反応の MD 計算が可能である．今回我々が開発した時間依存断熱状態 DFTB/MD 法は，ナノカーボンの化学反応動力学を半定量的に表現できる SCC-DFTB 法[5]のハミルトニアン HDFTB (R,ε(t) = 0)を H (R, ε(t) = 0) として用いることで，カーボンナノチューブ等の大型ナノカーボンのレーザーパルス誘起ダイナミクスを効率良く計算できる．なお，今回の計算では時間依存断熱電子基底状態 (n = 0) のみを考慮した．

3．結果と考察
直線偏光 (z 軸偏光) のガウス関数型パt = 167.5 fs
t = − 70.0 fs
t = 107.0 fs
ルス (半値幅 Tp = 70 fs, ピーク強度 Ipeak =7.0 × 1014 W/cm2, 中心波長 λ = 1800 nm)
を C60 に 1 つ照射した時のトラジェクト
リのスナップショットを図 1 に示す．
zzz
この時 C60 が偏長・扁平構造を交互に
xxx
取る hg(1)モード (調和振動数 270 cm−1)
図 1：時間依存断熱状態 DFTB/MD 法を用いてが主に励起された (図 1)．トラジェクト計算された z 軸偏光近赤外フェムト秒パルス
リのパワースペクトル (図 2，実線) もこ(Tp = 70 fs, Ipeak = 7.0 × 1014 W/cm2,
λ = 1800の結果を支持する．また，Tp = 30 fs (図nm) による C60 の非共鳴振動ダイナミクスのス
2，点線) とした場合は，550 cm−1 付近のナップショット．パルス強度が極大になった時
ag(1)モード (調和振動数 566 cm−1) の相間を t = 0.0 fs と定義した．C60 が偏長・扁平構
対ピーク強度が Tp = 70 fs の時よりも高くなっている．以上の結果は我々が最近報
造を交互に取る hg(1)モード (調和振動数 270
告した DFT 法 (B3LYP/3-21G 法) に基づ
cm−1 ) が主に励起されている．
く同様な MD 計算 (DFT/MD 法) の結果
hg (1)[3]と良く一致した．また，表 1 に示すよ
うに，近赤外パルスとの相互作用によって C60 に 注 入 さ れ る 振 動 エ ネ ル ギ ー
E in DFTB は DFT/MD 法による結果 (E in DFT)
[3]と半定量的に一致した．さらに，時間依存断熱状態 DFTB/MD 法による MD 計
算の所要時間は DFT/MD 法よりも 2 桁ほag (1)
ど少ないことがわかった．
以上より，時間依存断熱状態DFTB/MD−1
Frequency /cm
法は低い計算コストでDFT/MD法によるナノカーボンの近赤外フェムト秒パルス 図 2：時間依存断熱状態 DFTB/MD 法によるトラ
誘起振動励起ダイナミクスの結果を半定 ジェクトリのパワースペクトル．Tp = 70 fs (実量的に再現できる有効な手法だといえる． 線) と Tp = 30 fs (点線) の双方で 250 cm−1 付近
当日はDFT/MD法の結果[3]と比較しなが の h (1) モードが強く励起されている．
g
ら，2つのガウス関数型パルスを時間差を変えて照射した場合の振動ダイナミクス 表 1: ガウス関数型近赤外フェムト秒パルスにも含めてより詳細に議論する予定である． よって C60 に注入された振動エネルギーE in：
DFTB/MD 法(E in DFTB) と DFT/MD 法 (E in DFT)
参考文献：[1] R. Stӧhr et al., ACS Nano 5, [3]の比較
5141–5150 (2011); [2] T. Laarmann et al.,Tp / fs
E in DFTB / eV
E inDFT/ eV [3]Phys. Rev. Lett. 98, 058302 (2007); [3] N.
30
16
19
Niitsu, K. Yamazaki et al., J. Chem. Phys.
70
20
30
136, 164304 (2012); [4] Y. Sato, H. Kono, S.
Koseki, and Y. Fujimura, J. Am. Chem. Soc.
125, 8019–8031 (2003); [5] M. Elstner, et al., Phys. Rev. B 58, 7260-7268 (1998).

3P-125
光励起による C60 フラーレンの解離ダイナミクスのシミュレーション(東北大院理 1, 東北大多元研 2) 中村尭祉 1, 新津直幸 1, 菅野学 1, 河野裕彦 1, 上田潔 2

Simulation of the fragmentation dynamics of photoexcited
C60 fullerene(Graduate School of Science, Tohoku University1, Institute of Multidisciplinary Research for
1
Advanced Materials, Tohoku University2) Takashi Nakamura ,Naoyuki Niitsu1, Manabu Kanno1, Hirohiko Kono1, Kiyoshi Ueda2【序】近年、高強度レーザーパルスを用いて分子の反応を制御する実験が行われている。フラーレンC60では光子エネルギーやパルス長を変化させることにより、イオン化と解離などの過程を制御できることが実験的に知られている。例えば、波長800 nm程度の近赤外レーザーを用いた場合、パルス長が5 psでは様々な解離生成物が得られるが、パルス長が25 fsと短い場
合は多価カチオンが生成し、解離は抑制される [1]。一方、X線自由電子レーザー(XFEL)を用いた場合は内殻軌道から電子がイオン化され、その後の電子再配置によってオージェイオン
化が起こり、多価イオン化すると考えられている。さらにイオン化後の親カチオンは余剰振動エネルギーを持つため、近赤外レーザーでは観測できなかった12価を超える超多価カチオンの解離が観測できるようになると考えられている。近赤外レーザーによる励起の場合とは異なり、XFELでは親カチオンの余剰エネルギーを見積もることができる。本研究では、XFEL
により高い余剰エネルギーを持ったC60の多価カチオンからフラグメントが解離する動力学
過程を理論計算により調べた。
【手法】C608+, C6010+, C6012+を構造最適化し、様々な余剰振動エネルギーに対して計算を行い、解離ダイナミクスの価数・エネルギー依存性を調べた。また、XFEL によるイオン化を想定
して C608+に 40, 80 eV、C6010+に 50, 100 eV、C6012+には 60, 120 eV の初期エネルギーを与え同様に計算を行った。これらのダイナミクス計算には密度汎関数緊密結合(Density-Functional
Tight-Binding; DFTB)法を用いた。
【結果・考察】いずれの価数においても初期エネルギーに依存した数種類の解離パターンが確認された。低いエネルギー(40~60 eV)を与えた場合、振動励起に伴Fig. 1 低エネルギー時(40~60 eV)に見られ

うStone-Wales(SW)転位(Fig. 1) [2]及びその後のC2の脱離 た Stone-Wales 転位. ①-②の結合が 90°回が観察された。C2の脱離には数百ps～数nsの時間を要 転し、構造が変化する.
した。C60n+(n=8, 10, 12)の安定構造とSW転位後の構造のエネルギー差は数eV程度であり、余剰エネルギーは十分に大きい。C2の脱離までにSW転位は複数回起こっており、頻度が高いほど脱離が早かった。通常、フラーレンの五員環は孤立五員環則により隣り合わないが、 Fig. 2 高エネルギー時(80~120 eV)に見らSW転位が起こることによりこの束縛が破れ構造

れた C60n+ (n = 8, 10, 12)の 2 パターンの解離過程

に歪みが生じ不安定化する。脱離時間とSW(a)

転位の頻度に相関が見られるのはこのため60と考えられる。
フラグメントの脱離に加え低エネルギー時

40

Count

高いエネルギー(80~120 eV)を与えた場合、

C608+_80 eV

20には見られなかった球形ケージの崩壊が確0認された。Fig. 2のように脱離→崩壊と、崩

0

5

10

壊→脱離の2種類の解離プロセスが観察され、ケージの崩壊は中性C60 に70 eV以上のエネ

30

伴い短くなっていったが、価数への依存性は

20
10
0

見られなかった。今後クーロン反発の効果を0

より正確に見積もれる計算方法を用いることで、崩壊の価数に対する依存性を評価する。
高エネルギーを与えた場合、いずれの価

C6010+_100 eV

Count

で起こった。崩壊時間はエネルギーの上昇に

5

15

20

(c)
15

C6012+_120 eV

10

Count

ルギーの上昇に伴い短くなる傾向にあった。

10

Fragment size

数・エネルギーにおいても主な解離生成物はC2であった。フラグメントの脱離時間はエネ

20

(b)

ルギーを与えた場合においても確認されている [3]。ケージの崩壊は数ps～数十ps程度

15Fragment size

5

また脱離時間の分散もエネルギーの上昇に0伴い小さくなっていくことが確認された。より高いエネルギーを与えた場合、解離がFig.
2のProcess2に偏るため、このような結果が得られたものと思われる。またFig. 3にC60n+(n=8,

0

5

10

15

Fragment size

20

Fig. 3 C60n+ (n = 8, 10, 12)に 80, 100, 120 eV を与えた際のフラグメント分布. (a) C608+に 80 eV, (b) C6010+に 100
eV, (c) C6012+に 120 eV を与えた.

10, 12)から最初に脱離したフラグメントの分布を示した。エネルギーと価数が大きくなるに
つれて、フラグメントサイズの分布に広がりが生じていることが分かる。これは球形ケージの崩壊後に脱離が起こった場合、大きなフラグメントが脱離しやすいことが原因であると思われる。C60、C70等のフラーレンに大きなエネルギーを与えた場合、エネルギーの上昇に伴い
C2～C5程度の小さなフラグメントの生成量が増大することが確認されている [4]。多段階の解離も確認されており、本研究においても解離後の親カチオンやフラグメントのダイナミクスを追跡する。
当日の発表ではエネルギーや価数を変えた結果についても報告し、それらへの依存性について詳細に議論する。また、フラグメントの持つエネルギー等についても報告する。

[1] I. V. Hertel et al., Adv. At.Mol. Opt. Phys. 50, 219 (2005).
[2] A. J. Stone and D. J. Wales, Chem. Phys. Lett. 128, 501 (1986).
[3] N. Niitsu et al., J. Chem. Phys. 136, 164304 (2012)[4] T. A. Beu et al., Phys. Rev. B 83, 024103 (2011).

3P-126
ケイ素・炭素混合ベンゼン原子価異性体に関する理論的研究(群馬大学) 召田温美, 工藤貴子

A Theoretical Study of the Valence Isomers of the Mixed Si/C-Benzene(Gunma Univ.) Atsumi Mesuda, Takako Kudo
【序論】
ある一つの分子骨格中に炭素とケイ素が混在する場合、炭素とケイ素の割合とその性質の推移との関係はどうなっているのだろうか？本研究ではこの素朴な疑問の答えを得るため、ベンゼンの原子価異性体を取り上げた。
C6H6 で表される化合物には、有機化学で有名なベンゼンを含め計 5 つの原子価異性体が存在する（下図）。それぞれ、芳香族性を持つと言われているベンゼン、二重結合を二つ持つデュアーベンゼンとビシクロプロペニル、二重結合を一つ持つベンズバレン、そして単結合のみから構成されるプリズマンである。この骨格の炭素を順次ケイ素で置換し、その構造や安定性の変化を ab initio 分子軌道法および密度汎関数法計算により調べた。

ベンゼン

デュアーベンゼン

ビシクロプロペニル

ベンズバレン

プリズマン【計算方法】
SinC6-nH6（n=1~6）で表される炭素・ケイ素混合ベンゼン原子価異性体の内、考えられる全ての組み合わせの分子構造を HF/6-31G(d)および B3LYP/6-311+G(d,p)レベルで構造最適化し、基準振動解析により平衡構造であることを確認した。また、各異性体構造の相対安定性を支配する因子を明らかにするため、ベンゼン型異性体の共鳴安定化エネルギーおよび各異性体のひずみエネルギーを homodemotic reaction energy 計算から見積もった。更に、ベン ゼ ン 型 異 性 体 に つ い て は 、 芳 香 族 性 の 指 標 の 一 種 で あ る NICS 値 の 計 算 を
B3LYP/6-311+G(d,p)レベルで行った。NICS 値は環の重心からの距離依存性を調べた結果、環の重心およびそこから 0.8 Å 離れた場所の二カ所での値を採用した。尚、プログラムは
Gaussian03 を使用した。
【結果と考察】ここでは、ベンゼン型およびプリズマン型異性体のみについて述べ、他の異性体について
は当日発表する。
ベンゼン型異性体

全ての組み合わせの構造が得られた。下の図１は、ベンゼンの相対エネルギーと構造である。また、環構造の変化としては、C の多い系では平面構造を取っていたが、Si が多く連なっているものは、非平面の構造であった。
ベンゼン系に関しては、共鳴安定化エネルギーの計算と NICS 値計算を行った。その結果、共鳴安定化エネルギー計算では、全体的に相対エネルギーと相関関係を得られた。一方、NICS 値は一般には負の値が大きいと芳香族安定化が大きいとされているが、NICS 値が小さいにも関わらず、エネルギーが不安定であったため、相対エネルギーとあまり良い相関関係は見られなかった。表１は、NICS 値の結果である。
プリズマン型異性体
プリズマン型異性体についても全ての組み合わせの構造が得られた。この系では、プリズ
マンは、三員環アルカン 2 つと、四員環アルカン 3 つの構成要素があるので、今回、三員環と四員環のひずみエネルギーを計算した。
その結果、プリズマンの安定要素は、炭素の多い系では三員環、ケイ素の多い系では、四員環と相関関係が見られた。また、炭素・ケイ素半分の、C3Si3H6 の系では、その中間の結果を得られた。これは、三員環アルカンのひずみエネルギーは、炭素を含む三員環アルカンの方がケイ素を含む三員環アルカンよりも安定であり、反対に四員環アルカンは、ケイ素を含む四員環アルカンのひずみエネルギーの方が炭素を含む四員環アルカンのひずみエネルギーより小さいため、より安定な方のひずみエネルギーに支配されていると考えられる。

0.0

0.0

0.0

2.4

8.3

6.6

11.8

13.7

9.7

図１ ベンゼン型異性体の相対エネルギー（kcal/mol）
計算レベル B3LYP/6-311+G(d,p)
表１ ベンゼン型異性体のNICS値 計算レベル B3LYP/6-311+G(d,p)
Si

Si

Si

Si

Si

構造
Si

Si

Si

Si
Si

Si

Si

Si

Si

Si
Si

Si

Si Si

Si Si

Si
Si

Si

Si

Si

Si Si

Si

Si

Si

Si

Si

Si

Si

Si

Si
Si
Si

相対エネルギー（kmol/mol）
NICS(0.0Å)(ppm)-7.942 -7.767 -7.120 -9.288 -7.719 -6.947 -10.003 -8.599 -7.109 -10.728 -9.976 -10.988-11.420

NICS(0.8Å) (ppm)-10.556 -9.102 -7.464 -9.737 -8.190 -6.280 -9.746 -8.253 -6.277 -10.010 -9.259 -10.674-10.581

0.0

2.4

11.8

0.0

8.3

13.7

0.0

6.6

9.7

3P-127
イオン液体中のイオン間ダイナミクスと動的不均一性に関する理論的研究(分子科学研究所) 石田

干城

Molecular Dynamics Study of Dynamical Properties in Ionic Liquids(Institute for Molecular Science) Tateki Ishida【序】陽・陰イオンの組み合わせからなるイオン液体はイオンのペアを変えることで様々な種類のものを合成することが可能である。したがって、様々な用途に応じたイオン液体を“設計”することを可能にするという目標に向けて、多種多様な研究が行われてきている。特に、イオン液体に関する理論的研究は実験研究の進展と共に近年、急速にその重要性を増してきている。
イオン液体の“設計”に向けての重要な点の一つと考えられるものとして、イオン液体の多様な特性を理解することが挙げられる。イオン液体の多様な特性には液体構造などを含む
静的な挙動に関するもの、また緩和過程や動的不均一性などの動的な挙動がある。これらの内で特にイオン間の相互作用、つまりイオン間ダイナミクスに影響をうけると考えられる動的挙動はイオン液体の性質を如実に表わすものであり、通常の液体や従来の溶融塩とは明らかに異なる性質の解明への鍵となる。動的な挙動を理解することはイオン液体の特性のコントロールという大きなステップとなることは明らかであり、その解析が望まれる。
我々はこれまでにも分子レベルでの理論的解釈を行うためにシミュレーションによる研究を進めてきた。イオン液体中における陽・陰イオン間の相互相関の解析や、イオン間での運動量移動とその際のイオン間ダイナミクスの追跡から、イオン液体中でのいわゆる「かご効果」についての解析を行い、さらに分極効果により「かご効果」が減少する傾向があることを初めて明らかにした。
これらのイオン間ダイナミクスは比較的短時間の領域（~ 1 ps）での相互相関関数を解析することで可能となるが、それに続く sub-diffusive region (regime) と呼ばれる領域から拡散領域へとわたる時間領域においてはさらに動的不均一性と呼ばれる動的挙動が出現し得ると言われ、この時間領域での緩和過程は時にはナノ秒オーダー以上にも及ぶ。多くの通常液
体では過冷却状態において出現するとされるこのような動的挙動がイオン液体では室温において現れ得ることが近年のシミュレーションによる研究からも示唆されてきており、系統的な研究と解析が必要とされてきている。以下では、これらの課題に対して本研究で行ったことに関して述べる。
【イオン間ダイナミクスと動的不均一性についての理論的研究】
分子動力学シミュレーションによる手法を用い、イオン間ダイナミックス、及び動的不均一性に関する解析を行った。具体的には対象とするイオン液体の系として、[BMIm][NTf2]を含む数種類の系を選び、分子動力学シミュレーションを実行した。シミュレーション結果
から、イオン液体中における陽イオンと陰イオンの速度などに関する相関関数を求め、また
角運動量相関関数、中間散乱関数や non-Gaussian parameter の計算を行った。これらの解析から、イオン間のダイナミクスに関してより分子描像を捕えた解析を可能とした。加えて、動的不均一性と空間的相関構造の関係についても研究を行った。
解析結果から、「かご効果」が顕著であるとされる短時間領域でのイオンの挙動はイオン間ダイナミックスを大きく反映するものであり、対称性なども含めたイオンの形状にも依存すること分かった。さらに、室温においてイオン液体は顕著な動的不均一性を示し、その緩和過程は数十ナノ秒のオーダーに及ぶことが示された。また、これらの遅い緩和過程は陽・陰イオンの拡散過程に大きく影響を及ぼすことが分かった。

【参考文献】(1) “Molecular Dynamics Study of the Dynamical Behavior in Ionic Liquids through InterionicInteractions”, T. Ishida. J. Non-Cryst. Solids, vol. 357, pp.454-462 (2011).
(2) “The Dynamical Properties on Ionic Liquids: Insights from Molecular Dynamics Study”, In Ionic
Liquids - New Aspects for the Future; Kadokawa, J., Ed.; InTech: Rijeka, Croatia, 2012; Accepted.

3P-128
テトラシアニド金属錯体とメチルビオロゲンからなる結晶の吸収スペクトル解析(立教大理 1,お茶大理 2)
小嶋亮平 1，森寛敏 2，望月祐志 1，上之原和佳 1，松下信之 1

Spectral analysis of crystals composed of [M(CN)4]2- {M=Pt,Pd,Ni} and Methylviologen(Rikkyo University1, Ochanomizu University2)
Ryohei Kojima 1, Hirotoshi Mori 2, Yuji Mochizuki 1, Waka Uenohara 1, Nobuyuki Matsushita 1【序】テトラシアニド金属錯体 ([M(CN)4]2- {M=Pt,Pd,Ni}) 及びメチルビオロゲン(MV2+)は、Ni
錯体を除き孤立系においては可視領域に吸収スペクトルが存在しない、しかし、無色である Pt・
Pd 錯体とメチルビオロゲンを組み合わせることで、有色結晶が得られることが報告されている。
また、同一構造ながら、色が異なる結晶が得られる場合があることも明らかとなっている。しかし、同一構造錯体が異なる光学的特性を示す原因は解明されておらず、我々は各種テトラシアニ
ド金属錯体とメチルビオロゲンを組み合わせた結晶について励起状態計算を行い色の発現機構の理論的解析を行った。この要旨では、Pt 錯体の場合についてまとめる。
【計算】単体（孤立系）である[M(CN)4]2-{M=Pt,Pd,Ni}, MV2+, MV+、及び一対の[M(CN)4]2- とMV2+の系 ([M(CN)4]2-＋MV2+) を系統的に計算した。構造は[M(CN)4]2-では LC-BLYP で構造最適化し、MV2+, MV+, [M(CN)4]2-＋MV2+では MV[M(CN)4]結晶の X 線結晶構造解析で得られた実験値を用いた。励起状態計算は SF/SO-CASPT2（スピン軌道相互作用の有無）で行い、基底関数
は、Pt,Pd,Ni は内殻省略近似基底の一種である MCPtzp を使用し、それ以外の原子は cc-pVDZを用いた。また、アニオンである[M(CN)4]2-では電子の広がりを考慮し、diffuse 関数を含んだ基底(MCPtzp+,aug-cc-pVDZ)を使用した。
【結果と考察】・[Pt(CN)4]2Fig.1 に[Pt(CN)4]2-の 8 電子 5 軌道の活性空間を用いた計算結果と吸収スペクトルを重ねた図を示す。
計算から得られた理論スペクトルと実験の吸収スペク
トルの波長領域が一致していることがわかる。この理論
スペクトルについて解析を行うと、スピン軌道相互作用
によって軌道が混合し各状態を構成していることがわ
かった。
・[Pt(CN)4]2-+MV2+
Fig.2 に MV [Pt(CN)4]の単結晶偏光吸収スペクトルと [Pt(CN)4]2-+MV2+の SO-CASPT2 計算結果との比較を示す。また、Table.2 に SO-CASPT2 の解析結果、

Fig.1 [Pt(CN)4]2-実験 1・理論スペクトル

Table.3 に SO-CASPT2 の状態を構成する SF-CASPT2 の波動関数キャラクターを示した。
SF-CASPT2 の各番号は SO-CASPT2 の各状態を構成する波動関数キャラクターに対応している。
Fig.2 から理論スペクトルと吸収スペクト
ル(a)の立ち上がりと一致していることがわかり、可視領域で得られている理論スペ
クトルの解析を行うと、吸収スペクトルは[Pt(CN)4]2-（供与体）から MV2+（受容体）への電荷移動励起であることがわかった。
また、同一構造ながら色が異なる結晶については MV+ で可視領域での吸収スペクト
ルを示す状態遷移で説明可能である。
Fig.2 MV[Pt(CN)4]スペクトル 2 と計算値(a)E∥DA カラム(b) E⊥DA カラム
Table.1 SO-CASPT2 計算結果
eV
cm-1
nm
f
0.00
0
1
2.74
22070
453 6.11E-04 6
2.74
22086
453 6.59E-06 6
2.74
22134
452 5.65E-03 6
2.80
22620
442 1.72E-02 2
3.00
24212
413 4.68E-03 7
3.01
24246
413 2.40E-03 8
3.01
24310
411 2.18E-02 7
3.03
24457
409 1.16E-02 8
3.42
27612
362 3.75E-03 7
3.43
27670
361 8.74E-03 9
3.43
27697
361 1.79E-03 9
3.48
28035
357 7.69E-04 8
3.53
28463
351 1.49E-03 9
3.53
28491
351 1.90E-02 9
3.54
28568
350 1.10E-02 9
3.57
28828
347 4.53E-04 5(c)水/DMF 混合溶液中
波動関数キャラクター
100%
78%
80%
81%
73%
36%
38%
45%
37%
33%
34%
36%
34%
47%
46%
41%
78%

7
7
4
7
8
7
8
7
9
7
7
3
8
4
8
8

7%
11%
7%
15%
25%
34%
20%
18%
33%
31%
29%
30%
29%
25%
25%
12%

Table.2 SF-CASPT2 波動関数キャラクター
波動関数キャラクター
12222000093%
22u2000d2%
222u00d0
2
22u200d0
96%
222u00d0
2%
22u2000d32u2200d090%u22200d0
9%4
222u00d0
94%
22u200d0
3%222200005u22200d090%2u2200d0
9%
222u00d0
6
22u200u099%
7
2u2200u078%
u22200u0
13%
222u00u0
8
222u00u090%
2u2200u09%
u22200u09
u22200u086%
2u2200u0
13%
10
22u2000u83%
22u20u00
17%
1 C. D. Cowman, H. B. Gray, Inorg Chem., 15, 11, (1976), 2823.
2 D. Shiota, N. Matsushita Chem. Lett. 37,(2008),399.

3
8
8
8
3
9
4
4
3
8
8
6
6
6
2
7

7%
8%
5%
9%
22%
14%
19%
16%
9%
17%
20%
15%
11%
11%
13%
4%

8

3%

2

3%

7
6
9
3
9
9
8
4
4
9
7
8
7
6

4%
2%
13%
11%
9%
15%
8%
9%
7%
8%
10%
9%
12%
3%

3

2%

4
5
5
3
6
6
2
5
4
7
3
3

2%
1%
5%
10%
7%
7%
5%
6%
2%
8%
4%
2%

※赤は 1 重項状態
2%
1%

黒は 3 重項状態

1%
1%

スピン,電子対

9%
1%

u,d,2 はαスピン,β3P129エチレンおよびベンゼン誘導体のH/D同位体化学シフトに関する理論的研究（金城大・医療健康 1，金城大・社会福祉 2，北陸先端大院・ナノセンター3）
杉森 公一 1，川辺 弘之 2，島原 秀登 3

Theoretical investigation of H/D isotope-induced chemical shift in
ethylene and benzene derivatives(Faculty of Health Sciences, Kinjo Univ.1, Faculty of Social Work, Kinjo Univ.2, Nano Material
Technology Center, Japan Advanced Institute of Science and Technology (JAIST) 3)
Kimikazu Sugimori1, Hiroyuki Kawabe2, Hideto Shimahara3【序】重水素置換に代表される同位体効果は，分子の平均構造や熱力学的諸量，反応速度に影響を与える。それらの理論的解析や帰属に，分子軌道（MO）法や密度汎関数（DF）法などによる量子化学計算は強力な方法の 1 つである。しかし MO 法や DF 法では，核は電子よりも重いため静止した状態（断熱状態）にあるとする Born-Oppenheimer（BO）近似に依っており原子核の量子性・質量は露に考慮されない。そのため，同位体効果による僅かな構造変化に伴う分子物性変化の予測には，分子振動の非調和性の考慮や，波束ダイナミクス・量子モンテカルロ法による取扱いが要求される。多成分量子化学理論により核の波動関数を含む量子化学計算の試みもあり，同位体効果の関与する磁気的分子物性へも適用されている[1-2]。
核磁気共鳴（NMR）から得られる化学シフトにも同位体効果が影響し、数 ppb のわずかな同位体シフト（isotope shift）が生じる。Jameson らは，ab initio 法によって化学シフトの同位体効果を予測した[3]。同位体置換による分子構造変化の影響はポテンシャルエネルギー曲面(PES)，磁気遮蔽定数曲面(MSS)上の変化として表すことができる。最近，我々は二原子分子
H2 / HD，HCl / DCl，NaH / NaD の 1 次元 PES を用い，Morse 非調和振動子の解析的振動波動関数のボルツマン分布から平均構造（平均結合長）を求め，有限温度下での同位体シフトの理論計算を報告した[4]。
本研究では，エチレンおよびベンゼン誘導体における重水素置換による同位体化学シフトについて理論的に研究する。
【計算方法】モデル分子として C2H2，C2H4，C6H6 の重水素置換体を対象とする。２置換体までを考えたエチレンの例を Figure 1 に示した。核の伸縮振動運動に由来する核間距離 R の偏位のみを考慮する（Figure 1 の例では、R2 および R4 が重水素置換による核間距離に対応するが、すべての R の偏位を考慮する）。まず，X-H (X=C, N, O) 原子間距離に対する PES V(R)とこれに随伴して変化する MSS σ(R)を求める。次に，得られた PES を Morse ポテンシャル
VMorse(R)にフィッティングし平衡核間距離 Re，解離エネルギーDe，曲率 α を求める。得られた Morse ポテンシャルから振動波動関数の解析解 Ψv(R) (v は振動準位)が得られ，核間距離 Rとの期待値<Ψv(R)|R|Ψv(R)>によって平均核間距離<R>が求められる。σ(Re)と σ(<R>)の差が核振動に由来する効果に相当し，σ(<RX-H>) − σ(<RX-D>)が 1 次の同位体シフト（primary isotope
effect）1Δ となる。なお，重水素置換の部位から化学結合を介してさらに１つ離れた原子の時期遮蔽定数 σ の変化は，2 次の同位体シフト（secondary isotope effect）2Δ と呼ぶ。Figure 1 のモデル 1 では，C-1 位の σ 変化は 1Δ であるが C-1 位に結合した H と C-2 位では 2Δ となる。
H
13

H

R1

D

C
H

H

H

R2

13

D

C
R3

H

H

D

H

D

2a

H

D

R4

13

R5

C

1

0

H

R4

13

D

C

R5

H

H

H

R4

13

C

R5

D

2b

2c

Figure 1. H/D isotopomers and isotopologues of ethylene molecule.

計算方法として DF 法(B3LYP 汎関数)，Coupled-Cluster Singles and Doubles(CCSD)法を用いた。NMR 計算には Gauge-Independent Atomic Orbital(GIAO)法[5]，基底関数系には三倍原子価
correlation-consistent 基底関数系 cc-pVTZ，電子状態計算には Gaussian 09[6]および CFOUR
ver.1.2[7]を用いた。
【結果】B3LYP，CCSD/GIAO による結果の一部を Table 1 に示す。平衡距離での 13C 遮蔽定数[ppm]，Morse 振動波動関数から得た平均距離での 13C 遮蔽定数[ppm]，同位体シフト[ppb]を示す。bond order の増加によって 1Δ は減少し，実験結果[8]の傾向とも一致する。ただし，エチレンでの過小評価は，高次の同位体シフトを求める際に顕著に影響する。構造変化と電荷分布に関連付けた詳細な議論は当日報告する。
Table 1. Calculated NMR shielding tensor σ (isotropic value) [ppm] at the equilibrium distance Re, theaverage distance and their isotope shift 1Δ [ppb] by using GIAO/CCSD/cc-pVTZ (T = 300 K).
<σC1>T

Isotopomer

σC1 ( Re )ΔC1

Exp.[8]0d ethylene

74.04574.967001d ethylene

74.04575.113-146-274

0d ethane

188.257

185.028001d ethane

188.257

185.294

-266

-284

0d benzene

67.23567.71600

1d benzene

67.23566.867-151-283【参考文献】[1] Y. Kita and M. Tachikawa, J. Mol. Struct. (THEOCHEM) .912, 2-4 (2009).
[2] Y. Kita, T. Udagawa, and M. Tachikawa, Chem. Lett. 38, 1156-1157 (2009).
[3] A. C. de Dios and C. J. Jameson. Annu. Rep. NMR Spectrosc. 29, 1–69 (1994).
[4] K. Sugimori and H. Kawabe, Int. J. Quantum Chem. 110, 2989-2995 (2010).
[5] R. Ditchfield, Mol. Phys. 27, 789-807 (1974).
[6] M. J. Frisch, G. W. Trucks et al., Gaussian 09, Rev. A. 02, Gaussian, Inc., Wallingford CT (2009).
[7] J. F. Stanton, J. Gauss et al., CFOUR Coupled-Cluster techniques for Computational Chemistry, version 1.2.
[8] S. Berger, Encyclopedia of Nuclear Magnetic Resonance, 2, 1168-1172 (1996).

3P130
異種金属原子による窒素分子を挟んだ逆サンドイッチ型錯体の電子構造に関する理論的研究(京大福井謙一研セ) 中垣 雅之，榊 茂好

Theoretical Study of Inverted Sandwich Type Complexes of Hetero
Transition Metals Bridged by N2(Kyoto univ. FIFC) Masayuki Nakagaki, Shigeyoshi Sakaki【緒言】β-ジケチミナト骨格を有する DDP (DDPH = 2-{(2,6diiso-propylphenyl)amino}pent-2-en; SCHEME 1)が配位した金属錯体は、有機分子を挟んだ逆サンドイッチ型二核錯体(ISTC)を形成する。これら ISTC は金属原子やサンドイッチされる分子の種類
によって多様なスピン状態をとり 1、その理解は理論化学のみならず材料科学の面からも興味深い。窒素分子を挟んだ錯体は特に興味深く、Cr 錯体では窒素分子は η2-side-on で配位し１重項状態が最安定であるが 1、Fe 錯体では窒素分子は η1-end-on で配位し７重項状態をとる 2。これまでの研究で、金属の種類による構造・スピン状態の違いが、電子配置の違いによる金属原子‐窒素分子間の結合の違いに起因することを明らかにして来た。Cr 原子と Fe 原子では窒素分子との結合の性質が大きく異なることから、その Cr-Fe 異核錯体の電子状態は興味深い。本研究では Cr, Mn, Fe 原子を組み合わせた異核錯体の理論計算を行い、スピン状態及び窒素分子の配位構造の違いを検討した。
【計算方法】DDP の置換基を水素原子に置き換えた AIP (AIPH = 1-amino-3-imino-prop1-ene; SCHEME 1)を配位子とした錯体、(μ-N2)[M1(AIP)][M2(AIP)]2 {(M1,M2) = (Cr, Mn),(Cr,Fe), (Mn,Fe)}の各スピン状態について、CASSCF 法を用いて構造最適化を行い、CASPT2 法を用いてスピン状態の相対安定性を評価した。CASSCF 法の active space には、金属原子の3d 軌 道 と 窒 素 分 子 の ２ つ の * 軌 道 を 用 い た 。 遷 移 金 属 の 内 殻 電 子 は
Stuttgart-Dresden-Bonn の ECP で置き換え、原子価軌道には triple zeta 基底関数を用いた。
AIP の N 原子には aug-cc-pVDZ を、その他の原子には cc-pVDZ を用いた。
【結果及び考察】(μ-N2)[Mn(AIP)][Fe(AIP)]及び(μ-N2)[Cr(AIP)][Fe(AIP)]のスピン状態の相対エネルギーを表１に示す。(μ-N2)[Mn(AIP)][Fe(AIP)]では 8 重項が安定である。Mn と Feの等核錯体はともに高スピン状態(Mn:9 重項、Fe:7 重項)が安定であり、窒素分子との結合も類似している。CASSCF 法の結果から異核錯体においても等核錯体と同様の電子状態である
ことが示された。また、η1-end-on 配位が η2-side-on 配位よりも 1.7 kcal/mol 安定であり、その値は Mn 錯体(−0.2 kcal/mol)と Fe 錯体(3.6 kcal/mol)のほぼ平均となっている。
(μ-N2)[Cr(AIP)][Fe(AIP)]ではエネルギーの近い２つの電子状態が存在する。B1 状態は１重項状態から９重項状態までが近い値をとり、η2-side-on 配位が η1-end-on 配位に比べて非常に安定である。一方 η2-配位の A2 状態、η1-配位の B2 状態では３重項から７重項までの３つのスピン状態のみが近い値をとる。これらの状態では η2-side-on 配位と η1-end-on 配位のエネルギー差が小さい。η2-配位の 9B1 状態が最も安定であるが η1-配位の中では 7B2 状態が安定であり、エネルギー差も小さい。
表 1 スピン状態の相対エネルギー(kcal/mol)(μ-N2)[Mn(AIP)][Fe(AIP)]
η1-end(μ-N2)[Cr(AIP)][Fe(AIP)]

η2-sideη1-endη2-sidemultiplicity A2 state B2 state multiplicity B2 state B1 state A2 state B1 state
10tet

26.5

12.1

9tet

19.8

12.8

15.1

0.0

8tet

0.0

1.7

7tet

2.7

13.8

2.6

1.0

6tet

5.7

6.4

5tet

4.1

14.4

2.2

1.8

4tet

9.6

9.7

3let

4.8

14.9

1.4

2.2

2let

20.2

16.8

1let

13.6

15.0

16.8

2.2

(μ-N2)[Cr(AIP)][Fe(AIP)]の CASSCF 法によって得られた自然軌道のうち、金属-窒素分子間の結合に関係する dyz 軌道、dxy 軌道、
N2 の*軌道から形成される５つの軌道を図１に示す。
B1 状態は Cr 二核錯体の電子状態に近く、金属原子の dyz 軌道と N2 の*z 軌道からなる(yz)結合のみをもつ。Fe-dyz 軌道は非結合性軌道として１電子占有されていることから、Cr-N2間の方が Fe-N2 間よりも結合が強い。一方、η2-A2 状態は Fe 二核錯体の電子状態に近く、(yz)結合に加え、*x 軌道 dxy 軌道からなる δ(xy)結合も存在する。非結合性軌道の Cr-dxy 軌道、が１電子占有され δ(xy)結合は Fe-N2 側に局在かすることから、Fe-N2 間の方が Cr-N2 間よりも結合が強い。軌道相関図や局在化軌道を用いた、より詳細な考察に関しては当日報告する。
9B1 状態(1.96)(1.00)(0.02)(1.00)(1.00)

Cr-dyz+N2-*z

Fe-dyzCr-dyz−N2-*z

Cr-dxyFe-dxy(1.55)(0.45)(1.26)(1.00)(0.74)

Fe-dyz+N2-*z

Fe-dyz−N2-*z

Fe-dxy+N2-*x

Cr-dxyFe-dxy−N2-*x

3A2 状態

図 1 (μ-N2)[Cr(AIP)][Fe(AIP)]の CASSCF 自然軌道。左が Cr 原子、右が Fe 原子、( )内は占有数を表す。
[1] Y. I. Kurokawa, Y. Nakao, S. Sakaki, J. Phys. Chem. A, 2010, 114, 1191, 2012, 116, 2292.
[2] W. H. Monillas et al., J. Am. Chem. Soc. 2007, 129, 8090.
[3] J. M. Smith et al., J. Am. Chem. Soc. 2006, 128, 756.