Abstract for the Annual meeting of Japan Society for Molecular Science 2012

2E01

Molecular calculations with the Slater-type orbitals(No affiliation1 and Universidad Autónoma de Madrid2)
Kazuhiro Ishida,1 Rafael López,2 and Jaime Fernández Rico2
Introduction In recent molecular physics and theoretical chemistry, most people use
basis sets consisting of Gaussian-type orbitals (GTO). However, it is well known that
2
GTO, exp(−α r ) , have, at least, two serious defects; i.e., (a) they have no cusp at
r = 0 and (b) decay too fast at r → ∞ . In contrast with GTO, the Slater-type orbitals(STO) have the right cusp and long distance behavior. It is clear that STO should be thenatural choice for molecular calculations, but the fact that they lead to difficult integrals
have seriously hindered their use in this field. Recently, two of us (R.L. and J. F. R)
showed that the difficult integrals (which are the four-center electron-repulsion integrals,
ERI) can be calculated with using the so-called STO-NG expansion with an enoughprecision [1]. They developed a computer program named as SMILES2003 for themolecular calculations with using STO as basis functions [2]. The accompanying
coordinate expansion (ACE) algorithm developed by one of us (K. I.) is capable tocompute such integrals extremely fast [3]. We revised the SMILES2003 to
SMILES2007 with using the ACE algorithm for the four-center ERI.
Table I. Speeding up ratio for the individual ERI by the ACE algorithm
ERI type

SMILES 2003

SMILES 2007

Parallel

Total effect(1s1s|1s1s)

1 .62 ms

992 μs117 μs

14 times faster(2p2p|2p2p)

20 .8 ms

4 .74 ms

212 μs98 times faster(3d3d|3d3d)

454 ms

18 .6 ms

614 μs

740 times faster(4f4f|4f4f)

18 .8 s50 ms

4 .0 ms 4700 times faster(5g5g|5g5g)

68 .7 s

84 ms

27 ms

2500 times faster(6h6h|6h6h)

260 s

236 ms

149 ms1700 times faster

Computation time was measured with a HITAC sr11000 computer.
The 64 processors were used for the parallel computation
For all four-center computation, the STO-12G expansion is used for 1s-, 2p-, and 3d-type
STO, the STO-11G is used for 4f- and 5g-type STO, and the STO-10G is for 6h-type STO.

Table I shows that the speeding-up ratio in the computation time from as in
SMILES2003 to as in SMILES2007. The parallel version of SMILES2007 is extremely
efficient, as seen in Table I.
Table II shows the total computation time for several molecules with using the parallelversion of SMILES2007. We are now tractable to the computation of middle-sizedmolecules like adenine with the STO as basis functions, as seen in Table II.
Table II. Total computation time of the four-center ERIs over STOs for several molecules
Molecule

SMILES 2007

Parallel computation

Basis (# of STOs)

single processor

with 64 processors

Methanol CH3OH (18 electrons; Cs)
VB1 (62 STOs)

79 sec.

3.6 sec. (22 times faster)VB2 (130 STOs)

990 sec.

35 sec.

(28 times faster)

VB3 (244 STOs)

7928 sec.

270 sec.

(29 times faster)Benzene C6H6 (42 electrons D2h)
VB1 (150 STOs)

789 sec.

44 sec.

(18 times faster)VB2 (300 STOs)

8699 sec.

410 sec.

(21 times faster)

Adenine C5H5N5 (70 electrons; Cs)
VB1 (220 STO2)

19490 sec.

1095 sec. (19 times faster)VB2 (425 STOs)

201128 sec.

8698 sec. (23 times faster)

The basis set of VB1 is a “triple-zeta + polarization”
The VB2 is a “quadruple zeta + double polarization”
The VB3 is a “quintuple zeta + triple polarization”
We may conclude that we have now an efficient tool to compute molecules with STO.
The next will be the geometry optimization of molecules. We need derivatives ofmolecular integrals for such purpose. We derive the derivatives of one-electron integrals
and present them at the meeting.
References[1] J. Fernández Rico, R. López, A. Aguado, I. Ema, and G. Ramírez, J. Comput. Chem.
19, 1284-1293 (1998)[2] J. Fernández Rico, R. López, I. Ema, and G. Ramírez, J. Comput. Chem. 25,
1987-1994 (2004)[3] K. Ishida, J. Chem. Phys. 109, 881-890 (1998)

2E02
非調和ポテンシャルに対する振動座標の最適化	  	
  （イリノイ大学化学科）	
  八木清，平田聡	
  
	
  
	
  

Optimal	
  vibrational	
  coordinates	
  for	
  anharmonic	
  potential	
  (University	
  of	
  Illinois	
  at	
  Urbana-­‐Champaign)	
  Kiyoshi	  Yagi	  and	  So	  Hirata	
  

	
  
[序]	
  非調和振動状態計算は，厳密解を与える Full	
   CI 計算を除き，座標系の取り方に依存する．現在，座標系には基準座標が最もよく用いられるが，基準座標が必ずしも厳密解への速やかな収束を与えるとは限らない．基準座標は 2 次のモード間カップリングを消すが，非調和カップリングを小さくする保証はないからである．例
えば，非局在化した基準座標より局在化座標の方が X-­‐H 伸縮振動モード（X	
  =	
  O,	
  N,	
  C）
の振動状態を良く記述できることが昔から知られている[1,2]．しかし，任意の系に対し最適な振動座標を求める方法はこれまでになかった．	
  	
  本研究において，我々は vibrational	
   SCF	
   (VSCF)エネルギーを座標変換に対して最小化することで振動座標を求める新しい方法を開発した．この方法では，調和ポテ
ンシャルにより相互作用する振動子に対して基準座標が得られ，一方，相互作用のない非調和振動子（2 原子分子）の集合に対して各振動子に局在化した振動座標が得られる．すなわち，本方法は 2 つの極限において厳密解が容易に得られる座標系を与える．従って，一般の分子系でも良い座標系が得られることが期待できる．得
られた座標を用いて vibrational	
   CI	
   (VCI)計算を実施し，従来の基準座標に基づく計算よりも速やかに厳密解へと収束することを示す．	
  
	
  
[方法]	
  振動回転相互作用を無視すると，振動ハミルトニアンは基準座標 Q により以下の式で与えられる，	
  
f
1
"2
(1)	
  
Ĥ = ! ! 2 +V (Q) ,	  
2 i=1 "Qi
但し f は自由度の数，V(Q)は Born-­‐Oppenheimer ポテンシャルエネルギー曲面であり，平衡点の周りでのテイラー展開で表す．基準座標を変換し，新しい座標系を定義す
る，	
  
f

Q! i = !U siQs ,	
  

(2)	
  

s=1

! に対し，全波動関数を 1 座標関数の直積の形
但し	
   U は直交行列である．得られた Q
で表し，	
  	
  
f

! = " ! (i) (Q! ) ,	  ! nVSCF (Q)
ni
i

(3)	
  

i=1

変換行列と 1 座標関数を同時に最適化する．1 座標関数と VSCF エネルギーは，	
  ff%" 1 "2( j)( j)$!
+ ! !n j V ! !n j '!n(i)i = ! n(i)i !n(i)i ,	
  
(4)	
  
2
!
'&
$# 2 "Qi
j!ij!i

EnVSCF = ! nVSCF Ĥ ! nVSCF ,	
  

(5)	
  

: Q2 : &'(&)*$+,-'$

HFHF
stretchingmodes
aredelocalized
shown
in Figure3.
norma
to eachmolecule.
λmolecule.
derivedAppendix
B,The
is calto
each HF
λin
AppendixB,coordinates
overinthe
twoHFmo
0 , are0 , derived

-'$
4
7over
4
coordinatesaretwomoleculesre: Q̃2 : &'(/)&$+,
culatedto be
9.0×10for9.0×10
Q̃Q̃2 ,Q̃theand
4.6×10forQ̃3symmculatedto7delocalized
be
for
and
Q̃HF
, system:
and4.6×10f
flecting
theC1 andsymmetry
of
the

1
2
2v
: Q2 : &'(&)*$+,-'$: Q̃2 : &'(/)&$+,-'$
flecting
thesymmetryof thesystem:andindicating
advantageof localsymmetric
HF
and
Q̃C
indicating
thecoordinates
advantage
localedge
HF
asymmetricstretchingof coorthe
4 ,2valso the: Q4 : &011)*$+,-'$ : Q4 : &011)*$+,-'$ : Q̃4 : &'.*)&$+,-'$: Q̃4 : &'.*)&$+,-'$and Q̃4 , also

asymmetric stretching coordinates of the edge molecule

dinates. Therefore,
theoptimizationof themolecules
VSCF
energydinates.
Therefore,
theoptimizationof theVSCF
e(Q1 and
Q
the
inner
(Q
2 ) and
4 and
-'$
-'$
(Q1 andQ2 )otherand theinnermoleculesQ). areOn
4 andfor3weakly&'.*)&$+,: Q3Normal: -'$&'./)&$+,
automaticallygenerates
local
coordinatesfor(Qweaklyand: Q̃:-'$optimized
automaticallygenerates
local
coordinates
the
hand,
the
optimized
coordinates
3
FIG. 3. (ColorFIG.
online)(a)and(b)
optimized
coor3.
(Coloronline)(a)Normaland(b)coor: Q3 : &'./)&$+,: Q̃3 : &'.*)&$+,
the other
the molecule.
optimized
coordinates
areAppendix
localized
non-interacting
anharmonic
oscillators.
non-interacting
anharmonic
tohand,each HF
λderived in
dinates of weaklyinteracting
HFinteractingaligned inHF
onealigned
dimension.
dinatesof weakly
in one dimension.
0 ,oscillators.
to
each
HF
molecule.
λ,derivedinAppendixB,
is cal70
The zero-point
energy(ZPE)
and(ZPE)thefundamental
freThezero-point
energyEach HF is separated
byis3separatedÅ.
culatedto be
for Q̃
andtheQ̃2fundament, andEach HF
by 3 Å.
1and
7 9.0×10
4 4.6×culatedto be by9.0×10forQ̃compared4.6×10for thoQ̃
1 and Q̃
2 , and-'$-'$quencies
obtained
VSCF
are
with
those
by
quencies
obtained
by
VSCF
are
compared
with
and
Q̃
,
also
indicatingtheadvantageof
local3&011)*$+,&'.*)&$+,
4: Q̃4 :: Q4 :and Q̃4 , also indicating
the advantage of local HF coor
: Q̃4 : &'.*)&$+,-'$: Q4 : &011)*$+,-'$

により与えられる．VSCF エネルギーは座標に依存するため，これを座標変換に対しVCI[4]-(8) in
Table
I. in
VCI[4]-(8)basedon thetwoon
coorVCI[4]-(8)Table I. VCI[4]-(8)
based
the
dinates.
Therefore,
the optimizationofthe twoVS
dinates.
Therefore,
optimizationof theVSCF energyて最小化することで座標を最適化できる．変換行列に対するエネルギーの最小化はdinates
is foundto isbefoundinthe
agreement
withinaccuracyofdinatestobe in agreement
withinaccuraTABLE I. The TABLE
zero-point
energy(ZPE)(a)andNormalthe(ZPE)
fundamental
I.(ColorThezero-point
energyand(b)
theoptimized
fundamental
automaticallygenerates
local coordinates
for w
FIG.
3.
online)andcoorautomaticallygenerates
local
coordinatesforweakly−1−1 Normal and (b) optimized coor−1 selection of coordinates is not an is- andFIG.(in3.frequencies(Coloronline)frequenciescm) ofofweaklyinteracting
HF
calculated
bycalculated
(in
cm(a)
) of weaklyinteracting
by−1 . Thus,cm
0.1non-interacting
cmthe. Thus, theanharmonic
selection of oscillators.
coordinates is notdinates
weaklyinteracting
HFalignedHF
in
one 0.1
dimension.
non-interacting
anharmonic oscillators.
Jacobi	   sweep
アルゴリズムにより数値的に実行できる．アルゴリズムの詳細は当日
dinatesofweaklyinteracting
HFalignedinone
dimension.
VSCF based on
normalandcoordinates,
the
VSCF
based
onoptimized
normalbyand
optimized andcoordinates,
and forthethe high
sue
level
calculation.
However,thelowthelevel
sue forthezero-point
high level energy
calculation.
However,
the
low
The
(ZPE)
and
fundamEach
HF
is separated
3 Å.
The
zero-point
energy
(ZPE)
and
the
fundamental fre
Each
HF
is separated
3 Å.
deviations
from
VCI[4]-(8)in byparentheses.
deviations
from
VCI[4]-(8)in parentheses.
method (=VSCF)
does
dependonthecoordinates.
The witmethod(=VSCF)
doesoncoordinates
quencies
obtained
bydepend
VSCF
arethecompared
発表する．	
  quencies obtained by VSCF are compared with those by
VSCFb energies
based
on normal
coordinates
are found
VSCF
energies
based
normal coordinates
are
VCI[4]-(8)Table onI.
VCI[4]-(8)oncoorthe
VSCFa
VCI[4]-(8)b VCI[4]-(8)VSCFa
VCI[4]-(8)in Table in
I. VCI[4]-(8)based on based
the two
得られた座標を用いて，振動状態計算を実行する．VSCF
はあくまで平均場レベル
with
large
deviations
witha meanabsolute
withlarge isdeviations
witha agreement
mean deviation
absolute
dev
dinates
found
to
be
in
withinTABLE
I.NormalThezero-point
energy(ZPE)
and the fundamental
Normalcoord.
Optimized
coord.
coord.
Optimized
coord.
dinates is foundagreement within accuracy ao
TABLE I. The zero-point energy(ZPE) and the fundamental−1to be in −1−1(MAD) byof 117.4cm.. Thus,
In strike
those
based
(MAD)
cm
In strike ofcontrast,
thoseis(in) of(0.4)
weaklyHF calculated−1
cmof−1117.4
the. contrast,
selection
coordinates
ZPE frequencies
8419.5frequencies(80.9)
8338.98338.5
ZPE
8419.5(80.9)
8338.9interacting(0.4)
8338.5であり，分光学的精度を得るには相関理論へ進む必要がある．本研究では
VCI計算(in cm) ofcmweaklyinteracting
HF
calculatedby
0.1 cm−10.1. Thus,theselectionofcoordinates
is not an is
VSCF
on andnormalandoptimized
coordinates,
the oncoordinates
on andoptimized
are inlevel
excellentwith thoptimized
coordinates
are inagreement
excellentagreemen
4125.74075.6(1.2)
4074.4
ν1 VSCF4125.7(51.3)
4075.6(1.2)
4074.4
ν1(51.3)
sue
for
the
high
calculation.
However,based
onbasednormaloptimized
coordinates,
and thesue for the−1high level calculation.
However, the low leve
−1
from
VCI[4]-(8)in parentheses.
を実行し，基準座標と最適化座標の比較を行った．以下では最適化座標に基づく
VCI
MAD of 1.5MADcm
This
clearlydemonstratetheν2 deviations
4136.2deviations
4075.6(1.0)
4074.6of. 1.5
cmresult. Thisclearlydemonstraνfrom4136.2(61.6)
4075.6(1.0)
4074.62(61.6)
method(=VSCF)
doesresult
dependon thecoordin
VCI[4]-(8)in parentheses.
method (=VSCF)
does dependonthe coordinates.
Theν34278.4 ν(193.5)
4087.04084.94278.4(193.5) (2.1)
4087.0(2.1)
4084.9impactofcoordinate
optimization
which
dramaticallyenimpactofcoordinate
optimization
which
dramatica
3ab
VSCF energies
based
oncoordinatesnormalcoordinates
a VSCFb
VSCF
energies
basedonnormalare
found
VCI[4]-(8)を oc-­‐VCI と表記する．	
  VSCF
VCI[4]-(8)ν44289.3 ν(199.8)
4087.0(2.5)
4089.54289.3(199.8)4087.0(2.5)
4089.54
hances the accuracy
ofaccuracy
VSCF.
hances
the
with
large
deviations
with absolutea mean absolute
Normalcoord.
coord.
with large
deviations
withof aVSCF.
mean
deviationcNormalcoord. 117.4
Optimized
coord. 1.5
MADc
117.4
1.5 Optimized
MAD
−1
	
  
−1
(MAD)
of
117.4
cmstrike. In
strike contrast,
th
(a)
ZPE
8419.5(80.9)
8338.9(0.4)
8338.5
(MAD)
of
117.4
cm
.
In
contrast,
those based
ZPE
8419.5 (80.9) 8338.9(0.4)
8338.5on optimized
coordinates
are in excellent
a Optimized for atheν1ground(51.3)
4075.6(1.2)4074.4
4074.4on optimized
coordinates
are in excellent
agreementagreewithstate.
fundamentals
theOptimized
for量体の分子内振動モー
theThegroundstate.
The
fundamentals
are the[結果]	
  開発した方法を水2(51.3)4125.7 4125.7
4075.6(1.2)are
ν1
B.
Water
dimer
Water
1.5
cm−1
. Thisresultdemonstrate
clearly demon
ν2the modals.
4136.2(61.6)(1.0)4074.6 B.4074.6
virtual energies
ofvirtualenergies
modals.4075.6
MADdimerof MAD
1.5
cmof−1. Thisresultclearlytheν24136.2(61.6)of the4075.6
(1.0)
b The
b The
ドへ適用し
たresults
．νポ
テ
ン
シ
ャ
ル
エ
ネ
ル
ギ
ー
曲
面
は
ν
4278.4(193.5)4087.0(2.1)
4084.9basedonnormalandoptimizedareinagreement
results
basedonnormalandoptimizedareinagreement
impactofcoordinate
optimization
which drama
34278.4(193.5)4087.0(2.1)
4084.9impact of coordinate optimization
which dramatically
en
3
−1
−1
within 0.1. ν4
0.1
cm . (199.8)
4289.3(2.5)4089.5 4089.5
hancesaccuracyof
VSCF.
-1 theν4cm c within
4289.3(199.8)4087.0 4087.0(2.5)
hancesaccuracy
of
VSCF.
QQ
Q-112
Q1 Q
Q
Q48Q
Q
: 3783cm:Q3939cm
c Mean
1Q
2
59Q
610
7Q
8Q
9 Q10 Q11 Q12
2 Q3 Q
4the
5
6Q37Q
11
MP2/cc-­‐pVDZ
レベルで
次テイラー展開ポテンシ
absolute
deviation.
Meanc4absolutedeviation.
c MAD
117.4
1.5
MAD
117.4
1.5
Q̃3 Q̃4 Q̃5 Q̃6 Q̃7 Q̃8 Q̃9 Q̃10 Q̃11 Q̃12
ャル関数を生成し，振動状態計算に用いた．図1に
a Optimized for the ground state. The fundamentals
Q̃3 Q̃4normal
Q̃5 Q̃6 Q̃7andQ̃8 Q̃optimized
4 shows
coordinates ofa Optimized
9 Q̃10 Q̃11 Q̃12
thefor the ground state. The fundamentals are the are Figure
B. dimer
Water
dimernormal andWatervirtualenergiesof isthemodals.
water
dimer
for
the
modes.
The coordina
norFigure4 intramolecularshows
optimized
virtual energiesof themodals.
(b) B.
The HamiltonianwithQFF
invariantwithrespect
The
Hamiltonian
with
QFF
is
invariantwithrespect
水素結合ドナー水分子のOH
伸縮振動モードを示す．
bb The results
results
onand
normaloptimized
are in agreement
based
onbasednormalareinthatagreement
mal coordinates
be forcharacterized
as the bending,water candimer
the intramolecularmodes. The
to the coordinatetransformation,
so optimized
thatand
VCI
based
on
to Thethecoordinatetransformation,
so
VCI
based
on
−1
−1
withinwithin 0.1
cm 0.1. cm .
基準座標は水分子の対称・逆対称伸縮振動となる．
andsymmetric
stretchingof bencoordinatescanbecharacterized
as12 the
Q74Q
Q
Q96OH
Q10
QQ
c Mean andQ1 Qmal
Q41Qasymmetric
Q
Q
7Q
811
9Q
11 Q
normal candoptimized
coordinatesareequivalent
thenormaloptimized
coordinates
areinequivalent
in the2 Q3and
52Q63Q
85Q
1210 Qmodes
absolute
deviation.
Mean
absolute
deviation.
each watermolecule,
thoughsymmetricOH
OHstretching
stretch- moand symmetric
andtheasymmetricfull VCI limit.full
However,truncated
VCI
does dependon depend
VCI limit.
However,truncated
VCI
does
on
また，対称伸縮振動はアクセプター水分子が
5%ほ
-1 Q̃
-1 Q̃ symmetricing modes
foundbedelocalized
water
OH st
Q̃slightlyQ̃:Q̃each3807
cm:Q̃3918cm
Q̃thoughQ̃Q̃3arethe coordinates,
thus theand
behaviorof behavior
convergence
the andcoordinates,
thus theof convergence
7 Q̃
811
9 Q̃
11 Q̃to12the other3Q̃to6
74 molecule,
85Q̃
96Q̃
10
1210the4 Q̃5 Q̃water
moleculeby∼6%
inweight.
The
optimization
ing
modes
are
found
to
be
slightly
tocoorthe
ど混ざり，非局在化している．一方，最適化座標はshowsandnormalanddelocalized
optimized
may be different.
Infollowing,
wefollowing,onthe
con-respect
shows 4normal
optimized
coordinatesobe thedifferent.
the
focus
on
the respect
con- Figure 4Figure
The
Hamiltonian
with
QFF
is we
invariantwith
The may
Hamiltonian
withInQFF
isfocusinvariantwithhastransformedtheOH
stretchingmodes
to localOH
watermoleculeby∼the
6%intramolecularin weight.
Themodes.
optimiz
water
dimer
for
vergenceof
VCI[m]-(k)withrespect
to
m.
water
dimer
for
the
intramolecularmodes.
The
norvergenceof
VCI[m]-(k)withrespect
to
m.
to
the
coordinatetransformation,
so
that
VCI
based
on
to
the
coordinatetransformation,
so
that
VCI
based
on
図 1.	
  水 2 量体の水素結合ドナー水分子
ドナー水分子の OH 結合へ局在化する．これは OH stretching
coordinates,
as shown
in Fig.
4 (c),aswhile
hasmaltransformed
thecharacterizedOH
stretchingmodes
tocoordinates
canbethe leavbending
coordinatescanbecharacterized
as locathe
All theelectronic
structure
calculations
werecarriedAll optimized
theelectronic
structure
were
normalandcoordinates
arecalculations
equivalent
in thecarriednormalandoptimized
coordinates
are equivalent
inmal
the
の OH
伸縮振動モードと調和振動数．
35,36
35,36
ing
the
bending
coordinates
unchanged.
present
re-while
stretching
coordinates,
asshown
in Fig.
4stretching(c),andsymmetric
andasymmetric
OHThestretchingmodeso
伸縮振動の強い非調和性を反映している．	
  andsymmetric
andasymmetricOH
out at fullthe VCI
MP2/cc-pVDZ
levelof theory
using
the
out
at
the
MP2/cc-pVDZ
levelof doestheorythelimit.
However,truncated
VCI
dependon
full
VCI
limit.
However,truncated
VCI
doesusingdependon(a)	
  基
準座標と(b)最適化座標．	
  
37
37
sult
indicates
strongof OH
ingthethebending
coordinates
unchanged.
The
prese
each watermolecule,
thoughthesymmetricOH
stretcheachwateranharmonicity
molecule,
thoughthestretching
symmetricO
Gaussian09
program,andthevibrational
structurecal-structure
Gaussian09
program,andthevibrational
calthe coordinates,andthus
thethusbehaviorof convergence
thecoordinates,
andthebehaviorofconvergence	  表 1 に VCI計算の結果を示す．VCI-­‐(n)は量子数を3838drivesthearecoordinates
be localized
toothesulting
indicates
thestrongof theOH
stretingwhich
modes
are
found
to
be
slightly
delocalized
to
modes
found
toanharmonicity
betoslightly
delocalizedtoculations
bybe
MaVi.
culations
byInMaVi.
maydifferent.
the following,
we focuswe
onfocusthe conmay
be different.
In the following,
onmodes,the contheOHbonds.
modes,
which
drivesthecoordinatestobelocalizwatermoleculeby∼6%
inweight.
The
optimizationwatermoleculeby∼6%
inweight.
The
opvergenceofVCI
VCI[m]-(k)with respect
to m. to m.
n 個まで励起を許した計算である．基準座標と最適化座標による
VCI-­‐(8)は 0.1	  vergenceof
VCI[m]-(k)with respect
theOHcalculations
bonds.
VCI
six-dimension
have
been
hasandtransformed
the OHinstretchingmodes
tomodes
local OH
has
transformedtheOH stretchingto
All theAllelectronic
structure
calculations
were carried
the electronic
structure
calculations
were VSCF
carried計算となる．表1から明らかなように，
cm-­‐1 の精度で一致し，従って，benchmark
35,36
carried
out
with
the
above
coordinates
forshown
VSCF
and VCI
calculations
inthe
six-dimension
havew
stretching
coordinates,asshown asinFig.
4intramolec(c),
while
leav
35,36
stretching
coordinates,
in Fig.
4 (c),
out at IV.
the MP2/cc-pVDZ
level of level
theoryusing theout
atthe MP2/cc-pVDZ
of theoryusing theIV. APPLICATIONSAPPLICATIONS
37
ular modes.
The
ZPE
andthefundamentalareingbending
coordinates
unchanged.
Theoutbending
with
above
coordinates
forpresent
theintram
37 vibrational structure caling thecoordinates
unchanged.
TherepGaussian09
program,and theGaussian09
program,and the vibrational structure
cal-thecarried
基準座標を用いた
VCI-­‐(1)はVCI-­‐(8)と比較し大きな誤差を生じている．特に，OH
逆frequencies
38
sult
indicates
thestronganharmonicityofOH
stretchingpresented
in
Table
II.
The
VCI-(8)results
based
on
the
ular
modes.
The
ZPE
andfundamental
frequenci
38
sult
indicates
thestronganharmonicityofOH
culationsby MaVi.
culationsby MaVi.fluoride
A. Weakly-interactinghydrogenA.
Weakly-interactinghydrogenmodes,presented
which
drivesthecoordinates
toa belocalizedtotwo coordinates
are ingood
agreement
withmaximum
in
Table
II.
The
resultso
対称伸縮振動は誤差が大きく，
ν5 と
ν6 fluoride
modes,
which
drivestheVCI-(8)coordinates
tobased
be lo
−1
the
OH
bonds.
deviation
of
0.1
cm
.
Thus,
VCI-(8)servesasa
reftwocoordinates
are
in
good
agreement
with
a
max
表
1.	
  
V
CI-­‐(1),	
  oc-­‐VCI-­‐(1),	
  および
VCI-­‐(8)による水2
theOHbonds.
の誤差はそれぞれ
54.0,	
  115.0	
  cm-­‐1
であandVCI
calculations
insix-dimension
have been
erence toVSCF
evaluate
theof. VSCF
In order to
test
the developed
and program,deviation
ofaccuracy
0.1 VCI
cm−1
Thus,
VCI-(8)servesas
In orderto testthetheorydeveloped
theory and program,
VSCF
andcalculations
inlower-rank
six-dimension-­‐ and
IV. APPLICATIONS
量体の分子内振動モードに対する基本振動数(cmIV.
APPLICATIONS
carriederenceout
withthe above
coordinates
forVSCFthe intramolec
VCI.HFVSCF
andVCI[1]-(8)withnormalcoordinates
yield
we have first studied
modelsystem,in whichfourin HF
to evaluate
theaccuracy
of
and
lower
we have afirststudieda modelsystem,which four
る．一方，最適化座標を用いた
oc-­‐VCI-­‐
carried
out
with
the
above
coordinates
for
the
in
1)．括弧内に VCI-­‐(8)からの誤差，MAD に絶対平ular modes.
The
ZPETheandZPE
fundamental
frequencies
areular
modes.
andfundamental
frequ
A. Weakly-interactinghydrogenhydrogen
fluoride fluoride(1)は精度が格段に改善し，絶対平均誤A. Weakly-interactingin Table in
II. Table
The VCI-(8)
based
on ba
the
均誤差を示す．	
   presentedpresented
II. The results
VCI-(8)results
two coordinates
are in good
agreementwith a maximum-­‐1
two coordinates
are
in
good
agreement
withaまで減少した．
差は 37.5	
   から 4.6	  cm
deviationdeviationofoc-­‐VCI-­‐(1)	
  0.1 cmof−10.1. Thus,In ordertest tothetestdeveloped
theory andprogram,−1 VCI-(8) serves as a ref
VCI-­‐(8)	
  cm. Thus, VCI-(8) serves
In toorderthe developed
theoryandVCI-­‐(1)	
  program,	  	
  evaluate
the accuracyof VSCF ofandlower-rank
havewefirststudieda modela system,
in whichinfour
HF fourerenceoc-­‐VCI-­‐(1)の CI we
次元は系の大きさに対し
to evaluate
the accuracy
VSCF
and l
have
first studied
model system,
which
HF toerenceν 	
  1603.3	  (6.6)	
  
1602.3	  (5.6)	  
1596.7	
  て線形であるため，大規模な分子系に対ν 	  1645.5	  (5.7)	
  
1644.6	  (4.8)	  
1639.8	  
して非常に有効な方法である．	
  
ν 	
  
3640.4	  (21.1)	
  3621.7	  (2.4)	
  3619.3	
  
	
  
ν 	
  
3715.6	
  (22.7)	
  3695.1	  (2.1)	
  3692.9	  	
  [結論]	
  振動状態計算において，座標系
3836.8	  (54.0)	
  3788.4	  (5.6)	
  3782.8	  ν 	  をどのように選べば良いか，という古くν 	
  3909.6	
  (115.0)	
  3801.5	  (6.9)	
  3794.6	
  からの大きな問題を解決した．VSCF エ
MAD	
  
	
  	
  
37.5	
  
	
  	
  
4.6	
  
	
  	
  
ネルギーを最小化することで，座標は適切に局所化する．詳細は当日発表するが，座標が局所化することにより，ポテンシ
ャルエネルギー曲面はコンパクトな形へ変換され，VSCF 波動関数が VCI 計算の良い出発点となる．そのため，座標を最適化することで，相関レベルの振動状態計算は厳密解へより速い収束を示す．座標を最適化することで，振動状態計算の適用範囲が大きく広がる．	
  	
  [参考文献]	
  	
  
1.	
  	
  M.	
  S.	
  Child	
  and	  L.	
  Halonen,	
  Adv.	
  Chem.	
  Phys.	
  57,	  1	  (1984).	
  
2.	
  	
  B.	
  R.	
  Henry,	
  Acc.	
  Chem.	
  Res.	
  10,	
  207	  (1977).	
  
1
2
3
4
5
6

2E03
Na クラスターのプラズモン励起における集団性とラマン増強機構への寄与（分子研 1 ，総研大 2 ，国立台湾大 3 ）安池 智一1,2 ，野田 真史 1 ，信定 克幸 1,2 ，林 倫年 3Collectivity of plasmonic excitations in sodium clusters and its relation to Raman enhancement(IMS1 , GUAS2 , NTU3 ) T. Yasuike1,2 , M. Noda1 , K. Nobusada1,2 , M. Hayashi3【序】金属ナノ粒子は一般に可視領域に著しく強い吸収を持ち，それはバルク金属におけるプラズ
マ集団振動からの類推によって粒子内の局在プラズモン励起に起因すると考えられている．局在
プラズモンの強い光応答特性は，表面増強ラマン分光へ応用され単一分子分光を可能にするなど
その有用性に注目が集まっている [1]．一方で，構成原子数が数個から数十個であるような金属ク
ラスターにおいても，プラズモンに対応する励起状態の存在が知られてきた [2]．この状態への遷
移双極子モーメントは著しく大きく，和則の観点から粒子内のほとんどの電子がその励起に関与するとされている．本研究では，小さな Na クラスターにおけるプラズモンに対応する励起状態について，励起の集団性に注目した解析を行い集団性の発現メカニズムを明らかにする．さらに，
プラズモン励起波長近傍での近共鳴ラマンとしてのクラスター表面増強ラマン効果について，プ
ラズモン励起による電子運動と分子振動の相関の観点から議論を行う．
【Na クラスターにおけるプラズモン励起】密度汎関数理論 (DFT) に基づく実空間グリッド法計算プ
ログラムに Kohn-Sham 基底状態に対する周波数領域の線形応答理論を実装し，励起状態および光吸収スペクトルの計算を行った．Na 原子はノルム保存
250

擬ポテンシャルによって 3s 電子のみを露に取り扱い，

200

交換相関汎関数には Perdew-Zunger の局所密度近似

150を用いた．直線状の Nan (n = 6, 10, 14) に対して得る状態はごく少数に留まっていることが分かる．また，クラスターの構成原子数を増やすと，灰色で示したクラスター全体に非局在化した双極子的な遷移密度を持つ Σ, Π 励起の遷移モーメントが原子数に比
例して大きくなることが見て取れる．さらにクラス
ターサイズを大きくするとこの傾向はより顕著にな

5.18

Π

0
250

1

2

3

4

Na10

200

1.01

150
100

4.25

Σ

Π

50
0
250

1

2

3

4

Na6

200

ることが知られ [3]，Σ, Π 励起は連続的にバルクプラ

150

ズモンへと繋がると考えられる．遷移モーメントの100

集中現象は，対応する励起における電子運動の集団50性を示唆する．

1.02

50

Intensity ( Arb. Units )84 の双極子許容状態が存在するが，光吸収に関与すΣ

100られた吸収スペクトルを図 1 に示す．図示したエネル
ギー領域には n = 6, 10, 14 に対してそれぞれ 40, 63,

Na14

0

【プラズモン励起状態における集団性】各励起状態に

1.01

2.55

Σ

Π
1

2

3

4

Excitation Energy ( eV )おける particle-hole および hole-particle 対の励起振 図 1: 吸収スペクトルと遷移密度（図中数
幅 Xia ，Yia から，各励起空間における遷移密度行列を 字は集団性指標）

XiaYy∑ ia 2
Txi,a = ∑2 , Ti,a =X
ia ia
ia Yiaで定義する．T 行列の非ゼロ特異値を並べたベクトル λ は自然軌道表示の励起演算子ベクトルと∑見なすことができる [4]．λ に於ける有効励起対の数 n∗ = 1/ i |λi |4 は励起の集団性の指標として有用である．Na14 の Σ, Π 励起に対する n∗ は（2 つの励起空間での重み平均として）それぞれ

1.02, 5.18 となり，これらの双極型運動は集団性の観点において大きく異なることが分かった．Σ
励起の大きな遷移モーメントは集団運動の結果ではなく，一つの励起電子対が分子軸に沿った長距離に渡る双極型運動をすることにその原因がある．一方で Π 励起においては，励起電子対の大部分が関与した集団的な運動によって遷移モーメントの集中が起こっている．Π 励起の集団性は，縮重した複数の素励起が同程度の強さで結合することによって生じる [5]．

に比べて非弾性ラマン散乱の強度が著しく増大する．直線(a)Na4 クラスターの末端にピラジンの窒素原子が吸着した仮想的な構造を考え，Σ 励起に近共鳴な励起光を照射した際の電子の時間発展を，実空間 DFT 法によって追跡した．Na4の長軸方向を x とし，入射光は x 方向の直線偏光とした．
図 2(a) はこのときの誘起双極子の x 成分の時間発展 dx (t)dx(t) (eÅ)

小さな Na クラスターに吸着した場合にも，分子は孤立分子

Induced dipole moment (eÅ)

【Pyrazine-Na クラスターにおけるプラズモン励起と分子振動の相関】金属ナノ粒子と同様に
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
0

10

20

30

40

50

Time (fs)
Time
(fs)

(b)

（青は全体，赤はピラジン部分の寄与）で，Na4 のプラズモ
ン励起によってピラジン部にも同位相の誘起双極子モーメ
ントが生じていることが分かる．対応する遷移密度の ω 成 図 2: (a) 誘起双極子の時間発展 (b)対応する遷移密度分 ∆ρω (r) [図 2(b)，ω は入射光の角振動数] を見ると，これは Σ モードのプラズモン運動の空間分布と重なる位置に(a) Enhanced mode

ピラジンが吸着したためであることが理解できる．つまり，ピラジンの分子振動は Na クラスターのプラズモンの電子運動と空間的に相関し得る．非弾性ラマン散乱の強度は分極率の基準振動微分に比例し，これは例えば xx 成分について以下のように書き直すことができる．
∫ ()∂αxx (ω)∼∆ρωQ(0) +ϵQi (r) − ∆ρωQ(0) −ϵQi (r) xdr∂Qi(b) Non-Enhanced mode

右辺被積分関数の括弧内は平衡点から基準振動 Qi の方向に ±ϵQi ずらした構造で求めた遷移密度の ω 成分の差であ

図 3: 基準振動方向の変化に伴う遷移密度の変化 (a) 増強モード：1557
る．図 3(a),(b) はそれぞれ増強モード，非増強モードに対 cm−1 (b) 非増強モード：1260 cm−1する差遷移密度分布であり，増強モードではクラスター全り，この差遷移密度分布によってラマン散乱の強度が決ま

体に非局在化していることが分かる．つまり，ある振動モード方向の構造変化がプラズモンの電子運動に大きな影響を与える場合にそのモードは著しい増強を受けると理解することができる．吸着分子の振動モードに対するプラズモン励起の敏感な応答性は，プラズモン励起がクラスター全体に非局在化していることに起因している．
【参考文献】[1] Y. Fang, N.-H. Seong, D. D. Dlott, Science 321 (2008) 388. [2] C. Yannouleas,

R. A. Broglia, PRA 44 (1991) 5793. [3] J. Yan, Z. Yuan, S. Gao, PRL 98 (2007) 216602. [4] I.
Mayer, CPL 437 (2007) 284. [5] T. Yasuike, K. Nobusada, M. Hayashi, PRA 83 (2011) 013201.

2E04
周期的分子性結晶に対する QM/MM 法の開発：
SCF 点電荷による結晶効果のモデル化(京大福井謙一研究セ 1) 青野 信治 1、榊 茂好 1*

Proposal of QM/MM method for Periodic Molecular Crystal：
Crystal：Consideration of Crystal Effects by SCF point charge(Fukui Institute, Kyoto Univ.1) Shinji Aono1, Shigeyoshi Sakaki1*【序論】

理論化学における重要な課題の１つに、溶媒や結晶といった周囲の環境からの影響を考慮して注目する分子の構造や電子状態を検討する事が挙げられる。特に近年、分子デバイスの開発に対する要請を背景に、結晶およびアモルファスなど固体中での電子励起や発光、スピン状態などの分子的性質が実験、理論双方で研究対象として多くの興味を持たれている。
周期的分子性結晶を扱う際に有用な手法として、主に次の２つの理論的手法が挙げられる。１つは周期的密度汎関数理論であり、もう１つは QM/MM 法である。前者においては、注目する分子に周期的並進操作を施す事で無限サイズの結晶を直接計算する事ができ、特に注目する分子とその周囲からの集団的な影響を等価な計算レベルで記述できる点が長所に挙げられる。しかしながら密度汎関数理論に基づくため、多参照性が求められる電子励起状態などの記述において原理的に限界がある。一方、後者は周囲の分子集団を古典的モデルによって記述する点において、注目する分子に掛る結晶効果の記述精度に限界があるが、注目する分子の励起状態に対して高精度な電子状態計算を適用する事が原理的に可能である点が長所として挙げられる。
本研究では周期的分子性結晶中での励起状態の理論化学研究を目標として、
今回まず QM/MM 法の枠組みにおいて、基底状態下での結晶分子の構造および格子ベクトルの最適化法を開発した 1。
この周期的分子性結晶の最適化を、水分子が O サイト向きで Pt 金属に配位した
Pt(CN)4・2H2O 2 および水分子が H サイ
ト向きで Pt 金属に水素結合した transPtCl2(NH3)(N-glycine)・H2O 3 という２

図１基本単位格子内の分子構造 (a) Pt(CN)4・2H2O (C2 対称性)

つの系へ適用し（図１）、また孤立した(b) trans-PtCl2(NH3)(N-glycine)・H2O (Ci 対称性)

クラスター系での最適化構造と比較して、得られた構造と電子状態の違いについて検討した。
【計算方法】 本研究における QM/MM 法は周期的分子性結晶の基本単位格子を QM 領域、あるいは
基本単位格子内に対称性がある場合は対称操作を行う事で基本単位格子を再現できるような QM 領域を選ぶ事とした。これは注目する QM 領域に対して並進および対称操作を行う事で周期的分子性結晶全体を再構成できるからである。また QM 領域以外の結晶部分を MM 領域と見做し、QM エネルギー
EQM と QM/MM 相互作用エネルギーEQM/MM の 1/2 の和（式(1)）を変分原理によって最小化するように、注目する QM 領域の電子波動関数ΦQM を決定する。

E crystal ∝ E QM +1
1
E QM/MM = Φ QM Hˆ gas Φ QM + (E ES + E LJ )
2
2

式(1)

本研究は式(1)において QM/MM 相互作用を静電項と Lennard-Jones 項の和で表現し、特に静電項に含まれる MM 分子電荷に関しては QM 分子電荷（本研究では ESP 電荷）と並進および対称操作を介して等価になるように Fock 行列を定めた。QM および分子 MM 電荷を SCF 法で求め、更に解析的
gradient を用いて分子座標と格子ベクトルおよび対称操作に関する最適化を行った 1。
（Scheme）

Scheme : 周期的分子性結晶における分子電荷 Q、分子座標 r、格子および対称操作 t の決定

この基底状態下の結晶分子および格子の構造最適化法を B3PW91 汎関数による DFT レベルで上記の Pt 錯体へ適用した。Pt は 60 個の内殻電子を Hay-Wadt ECP に置き換えて（11111/11111/111/1）
基底を、その他の原子には 6-31++G**基底を用いた。Lennard-Jones 変数には AMBER および SPC
モデルの値を用いた。
【結果】 テトラシアノ錯体、グリシン錯体での基本単位格子は４分子により構成されているものの、
2 の２分子である。またテトラシアノ錯体の 1、22対称性を持つため独立な分子は、図１に示した 1、2間には水配位子の H サイトと CN 配位子の N サイトを介した水素結合が２つ、11、22’ 間には１つ存在、グリシン錯体の分子間では水分子の H サイトとグリシンの COOH 基の O サイトおよび水し（図２）分子の O サイトとグリシンの COOH 基の H サイトを介した水素結合が存在する（図３）。孤立したテトラシアノ錯体の単量体では最適化構
造は結晶構造中と一致するが、二量体、四量体では Ow(1A)-Hw 距離（図２）が伸び、プロトン移動性が現れる。一方、結晶を考慮した QM/MM 計算ではプロトン移動性は現れない。また、孤立したグリシン錯体の単量体では水分子が Cl 配位子と NH3配位子の間で水素結合を作り、結晶構造中と大きく異なる。二量体、四量体となるに

図２：テトラシアノ錯体 (a) 基本セル内 (b) 基本セル間水素結合

従って結晶構造に近付くが、水分子の相対位置や Pt(11)-Pt(11’) 距離は依然として結晶 構 造 中と は 異な る 。結晶 を 考 慮し た
QM/MM 計算によって得られた最適化構
造はこれら２つの錯体共に実験結果を再現し、本 QM/MM 法は周期的分子結晶中の基底状態をとる分子のモデル化に有用と考えられる。

図３：グリシン錯体 (a) 水分子の相対位置 (b) 分子間水素結合[1] : S. Aono, S. Sakaki, Chem. Phys. Lett. 2012,2012 544, 77[2] : C. Mühle, J. Nuss, M. Jansen, Z. Kristallogr. New Cryst, Struct. 2009,2009 224, 9.
[3] : S. Rizzato, J. Berges, S. A. Mason, A. Albinati, J Kozelka, Angew. Chem. 20102010, 49, 7443.

2E05
分子錯体形成における電荷移動項と分散項の競合・共同効果(慶大･理工、豊田理研)

岩田末廣

Cooperative and Competitive Roles of Charge-Transfer Terms in Molecular Cluster
Formation(Keio Univ. , Toyota Phys. Chem. Res. Inst.) Suehiro Iwata
序

この十年ほどの間、多くの閉殻分子(原子)から構成されているクラスターの結合エネルギーを｢基底関数の不均衡」に伴う誤差を避けながら効率よく精度高く計算する理論の構築とプログラム開発を進めてきた。不均衡に伴う誤差は、一般的には BSSE(basis set superposition error)と一まとめて呼ばれている。しかしながら、一電子基底関数の不均衡(Orbital basis inconsistency, OBI)と多電子関数の不均衡(Configuration basis inconsistency, CBI)という二つの性質の異なる不均衡からなっており、当然のことながら、適切に克服する理論手法には違いがある。これまでの研究で[1]、局所射影
分子軌道法(LP MO)による 1 電子励起 3 次摂動論(3SPT)によって、非常に時間のかかる Counterpoise 手順を使わずに、BSSE を含まない結合エネルギーを計算すること出来るまでに至った。分子間相互作用で重要な分散項(dispersion)の計算には多電子関数を多くのスレーター関数の重ね合わせで表さなければならず、CBI を考慮しなければならない。最近利用されている近似法は、構成分子(原子)の周波数依存分極率から見積もる方法であるが、CBI がないけれども、分子軌道に重なりがある領域との接続法に論理的な問題があり、経験的なパラメータの導入が避けられない。LP MOに基づく摂動法(PT)では、被占軌道も励起軌道もクラスターを構成する個々の分子に局在している。
その結果、励起電子配置を選別することによって、CBI を避けることが出来る。

方法 LP MO PT では、分子クラスターの波動関数を次のように展開する。

摂動法によって係数 A を求めることによって、各項からの寄与を計算することが出来る。 さらに、CT(電荷移動)と分散項は分子対 X-Y に、LE(局在励起)と分子内相関項は各分子 X に分割して計算することが出来る。相互作用項の各項は｢観測量｣ではなく、したがって成分分割は一義的ではなく、その方法、用いる多電子波動関数、一電子基底関数に強く依存する。したがってその意味づけには常に注意を払い、また、論文間の比較にも制限がある。

図 1 電荷移動項と分散項。 H 授与と H 供与型分子の型による分類

計算結果 水クラスター内の水素結合系に適用する。本報告で計算する系の構造は、MP2/aug-ccpVDZ によって最適化されている。クラスター内分子の構造も最適化してあるが、counterpoise 補正は適用されていない。ほとんどは Xantheas とその共同研究者によって報告されているクラスターである。
図１は、水クラスター(H2O)11 の 11 個の異性体(10kJ/mol 以内）の中の水分子対の電荷移動項 ECTと分散項 Edisp の関係を示している。顕著なことは、強い結合(両項が-10kJ/mol より強い）とき、両項はほとんど直線的な関係を持っている。絶対値の値も-10kJ/mol 程度で逆転している(縦軸と横軸の尺度が異なることに注意)。クラスター内の水分子は、水素供与結合している水素の数 m と水素授与している水素の数 n で分類することが出来る(DmDn)。左図では各水素結合対を水素授与分子のタイプで、右図では水素供与分子のタイプで区別して色づけしてる。D2A1 水分子が水素授与し、D1A2
水分子が水素供与している時、強い水素結合が形成されていることが明らかになる。このことは、OH
伸縮振動数の低周波数移動や OH 結合距離の計算値で知られていたことではあるが、LP MO PT
計算によって明瞭に示すことが出来た。ECT の値は、水素結合の強さの尺度になり、水素結合の有無の判断にも利用できる。クラスター全体の結合エネルギーへの寄与では、遠距離相互作用である
Dispersion 項が最大で、ついで本報告の定義での電荷移動項、ついで静電的な(分極を含む）相互作用となっている。図 1 から明らかなように、全体としては電荷移動項と Dispersion 項は共同してクラスター形成に寄与しているが、異性体毎に詳しくみると、大小に差があり、競合して異性体間エネル
ギーの差をつくっている。図 2 に異性体 X434 の主要水素結合に対する両項の値を示す。D1A3 タイプの水分子をもつ特異な構造となっている。
水クラスターの異性体の数は、水分子の数ととも急速に大きくなる。例えば、n=11 では、区別可能な digraph(水素結合を方向性もつ矢印で表す)は 460,066,726 と数え上げられている。もちろんポテンシャルエネルギー曲面上の極小値は少なくはなる。極小構造の探索にはいろいろな手法が用いられている。安定な異性体の構造と、その中に見いだされる水分子のタイプ、
水素結合対のタイプとの相関を調べ、表にした。相対エネル
ギーの順は MP2 と LP MO 3SPT+Disp では一部異なるが、その差は一つを除いては 1kcal/mol 以内である。

図２ 最安定異性体。
W8 は D1A3[1]Iwata, J.Phys.Chem.A 114 (2010) 8697; [2] Iwata, PCCP 14(2012) 7787
ClusterThenumberof HBsRelative/kJ mol-1
aug-cc-pVDZ

H donor H acceptor
D1A2-D2A1

D1A1

D2A1

D1A2

D2A2

5
4
3
4

0
1
2
1

6
5
4
5

b

5
5
4
5

0
0
1
0

17
16
16
16

0.00
1.45
1.70
2.69

3SPT+
Disp
0.00
-1.38
-2.07
-2.85

16

3.04

5.91

3

2

4

4

1

16
16
16

3.94
6.26
6.32

6.68
-0.67
10.60

4
2
3

1
3
3

5
3
3

5
3
3

0
2
2

16

7.69

10.74

3

2

4

4

1

17

7.81

9.15

2a

2

3

3

3

MP2

X434
X515
X443
X551
PX318
9
PI1226
X4412
PI2476PX420
5PI2877the number of waters

2E06
分極モデルを用いる分子動力学法と溶液理論の結合による溶媒和自由エネルギーの計算（東北大院・理）鈴岡大樹、高橋英明*、石山達也、森田明弘

Computation of the solvation free energies by the molecular simulations
with a polarizable force field combined with a theory of solutions(Graduate School of Science, Tohoku University)
Daiki Suzuoka, Hideaki Takahashi*, Tatsuya Ishiyama, Akihiro Morita【はじめに】 凝縮系の化学反応の経路を理論的に予測するには、その過程に伴う自由エネルギ
ー変化を計算することが本質的に重要である。反応の自由エネルギーや活性化自由エネルギーは
熱力学サイクルによって、反応物や生成物、遷移状態の溶媒和自由エネルギーから計算可能である。最近、溶質－溶媒間相互作用ポテンシャルの分布関数を基本変数とする分布関数理論（エネルギー表示の理論（ER 法））が開発された[1]。ER 法によれば、溶液系と純溶媒系におけるエネ
ルギー分布関数を用いて、効率良くかつ正確に溶媒和自由エネルギーが計算できる。
溶液において溶質が大きく電荷分極する系では溶媒の分極が強く誘起されるので、溶質の溶媒
和自由エネルギーが著しく変化する。したがって、溶質だけでなく溶媒の電荷分極を自由エネル
ギー計算に取り入れる必要がある。本研究では、古典の分極モデルに基づく分子シミュレーショ
ンと ER 法を結合した新規な方法（MMpol－ER 法）を開発し、その有効性を検証した[2]。
【方法論と計算の詳細】 オリジナルの ER 法は溶質－溶媒間相互作用が二体的であることを前提として構築された理論である。従って、分極モデルに基づく分子シミュレーションに直接適用することはできない。MMpol－ER 法においては、ある中間状態を導入して溶媒和自由エネルギー
を二体的な相互作用による寄与  と残りの多体的な寄与に分割する。
    (1)

ここで導入する中間状態とは、溶質－溶媒分子のみが二体的に相互作用するような系である。純溶媒系からこのような中間状態に変化する際の自由エネルギー変化が  である。  の計算には通常の ER 法を直接適用することが可能である。一方、この中間状態から、すべての分子が分極
力場によって相互作用する溶液系に変化するときの自由エネルギー変化が  である。  は
QM/MM-ER 法[3,4]での多体の寄与の取り扱いを、分極する溶媒の系に拡張することで簡潔かつ厳密に計算される。具体的には、溶液系と中間状態の全相互作用エネルギーの差によって定義される分極エネルギーをエネルギー座標として導入する。すると、溶液系と中間状態のの分布関数
P(), P0()によって、溶媒の分極電荷に応答して溶質の電荷が揺らぐことによる自由エネルギー変化が次のように厳密に定式化される。
 P      W  
 P   

 0


   d  kB T ln (2)

ここで、kB はボルツマン定数、T は絶対温度であり、W()は規格化された任意の重み関数である。また、は原理的にはに依らず一定となるので、(2)式の右辺の角括弧内を R()と定義し、

その一定性を確認する。
本研究では、異なる中間状態を経由する二つの経路 A, B を導入し、それらから計算される溶媒
和自由エネルギーを比較することにより本方法論の有効性を検証する。経路 A では二つの中間状態を導入し、溶媒和自由エネルギーを三つの項に分割する。一方、経路 B では中間状態を一つだけ導入することで溶媒和自由エネルギーを二つの項に分割する。いずれも二体的な寄与  は一項のみで、残りは多体的な寄与として計算される。
本計算ではベンチマークテストとして、水およびメタノール分子の水への溶媒和自由エネルギ
ーを計算する。また、溶質および溶媒分子の分極は CRK モデル[5]で記述した。溶質の内部座標は固定し、溶媒は内部自由度を持つ 500 個の水分子で構成した。比較のため、数値的に厳密な熱力学的積分法（TI 法）を用いて溶媒和自由エネルギーを計算した。分子シミュレーションの統計
量は TI 法と経路 A が 3 ns、経路 B は 2 ns とした。
【結果】 図１に経路 B の多体の寄与の計算に用いた種々の関数を示す。R()は P()と P0()の重なりが大きい領域でほぼ一定となっている
ことが分かる。また、経路 A, B に伴う自由エ
ネルギー変化を表１に示す。比較のために、TI 法から得られた溶媒和自由エネルギーTIの値をあわせて示した。水およびメタノールの溶媒和自由エネルギーについて、本方法で求められる自由エネルギーが経路にほとんど依らないことが確かめられた。また、これらの値は TI 法から得られる値ともよく一致していることが分かり、本方法論の正確さが確認できた。さらに、誤差論から MMpol-ER 法が
TI 法の 1/10 以下の統計量で同程度の統計誤差を達成する事が予測され、これについても実際に確かめられている。今後、本方法を拡張することによって、溶質を量子化学的に記述する

図１．経路 B における、溶質と溶媒分子の電荷
揺 ら ぎ に 関 わ る エ ネ ル ギ ー 分 布 関 数 P(),P0()と重み関数 W().

新規な方法(QM/MMpol－ER 法)を開発する計画である。
表１．水およびメタノールの水に対する溶媒和自由エネルギー。括弧内の数値は誤差。
（単位：kcal/mol）

Species
H2O

CH3OH

参考文献

Path



1

2



A

2.6 (0.0)-3.9 (0.2)-4.0 (0.2)-5.3 (0.4)

B

0.7 (0.1)-6.5 (0.3)―-5.8 (0.3)

A

2.3 (0.0)-3.5 (0.2)-2.2 (0.1)-3.5 (0.2)

B

0.7 (0.1)-4.3 (0.2)―-3.7 (0.2)

TI-5.7 (1.3)-2.4 (1.1)[1] N. Matubayasi et al., J. Chem. Phys. 113, 6070 (2000). [2] D. Suzuoka and H. Takahashi et al.,J. Chem. Phys. (2012) submitted. [3] H. Takahashi and N. Matubayasi et al., J. Chem. Phys. 121, 3989 (2004). [4]H. Takahashi et al., J. Chem. Phys. 136, 214503 (2012). [5] A. Morita et al., J. Am. Chem. Soc. 119, 4021(1997).

2E07

励起状態と溶媒効果:界面錯体を用いた理論化学的研究
1.東大院・工，2.CMSI 拠点研究員，3.東大・先端研
城野 亮太 1, 2，藤沢 潤一 3，瀬川 浩司 3，山下 晃一 1

Solvent Effect on Excited States of Surface Complex
1. Department of Chemical System Engineering, The University of Tokyo
2. Computational Materials Science Initiative (CMSI)
3. Research Center for Advanced Science and Technology, The University of Tokyo
JONO, Ryota1,2; Fujisawa, Jun-ichi3; SEGAWA, Hiroshi3; YAMASHITA, Koichi1【序】
光誘起電子移動反応は，光合成・光触媒・太陽電池などで主要な反応であり，電荷分離状態の生成が効率に強く影響する．この反応を計算科学的に明らかにすることでより高効率な材料の設計が可能になるが，電荷分離状態について実験を再現し予測する標準的な方法は未知である．
我々は酸化チタンと TCNQ からなる界面錯体が，電荷移動遷移によって可視光から近赤外領域の光を幅広く吸収する新物質であることを実験的に報告し，hybrid 汎関数を用いた密度汎関数計
算によって電荷移動遷移を帰属した．[1] その一方で我々は，人工光合成分子を用いた研究で，電荷移動遷移を高精度 ab initio 計算と同等の精度で求めるためには長距離補正を施した密度汎関数法(LC-DFT)を用いる必要があることを確認している．[2] 本講演では，界面錯体の可視光領域の吸収を LC-DFT を用いて再現し，溶媒が電荷分離状態に及ぼす効果について報告する．
【方法】
酸化チタンナノクラスターはアナターゼ型の特性を有するように設計した Ti9O18[1]およびTi14O28 を用いた．界面錯体には TiO2-TCNQ を用いた．LC-PBE は人工光合成分子での計算で最も SAC-CI と励起エネルギーの一致が良かった．そこで，界面錯体について，溶媒を連続誘電体モデル(IEFPCM)で考慮し LC-PBE で構造を最適化し，LC-PBE を用いて時間依存密度汎関数法(TD-DFT)による励起状態計算を IEFPCM 環境下で計算した(Fig1a)．また，LC-PBE を用いて，真空中で構造最適化および励起状態計算を行った(Fig1b)．さらにアセトニトリル分子を露に配置した系の励起状態計算を行った(Fig1c, Fig1d)．
a

b

c

d

Figure 1 HOMO of TiO2-TCNQ at (a) LC-PBE+PCM//B3LYP+PCM, (b)LC-PBE//LC-PBE,
and (c)(d) HOMO of TiO2-TCNQ + AN at LC-PBE//LC-PBE level of theory.

【結果】Fig1 に示したように，励起状態計算時のLC-PBE による最高占有軌道(HOMO)はい
ずれの構造でも TCNQ に分布していた．一方でその励起スペクトルは Fig2 のようにそれぞれの計算条件で大幅に異なる結果を示した．特に IEFPCM で溶媒を近似した場合，
電荷移動遷移は酸化チタン内励起よりも長波長側には観測できなかった．真空中で最適化した構造における吸収は可視光領域で電荷移動遷移を示し，定性的には実験事実と一致した．ここで Fig1c のようにアセトニトリ
ル分子を配置すると，その吸収は短波長側にシフトした．これは電荷移動錯体の持つ双極子とアセトニトリル分子の双極子が，基底状態で安定化・電荷分離状態で不安定化する向きに配置された結果である．一方 Fig1d のように配置した場合はそれぞれの双極子が基底状態では不安定化・電荷分離状態では安定化するような相互作用をするため，長波長側にシフトした．実際の実験環境では溶媒分子
は多数存在しており，その構造を一意に特定することはできない．したがって実用的な電荷分離状態の計算を行うためには PCM のよ
うな連続誘電体モデルの改良が必要である．
当日は Ti14O28 を用いたアナターゼ(101)面の結果，詳細な吸収の帰属，及び他の溶媒についての結果を追加して報告する．

Figure 2 Calculated absorption spectra[1] R Jono, J Fujisawa, H Segawa, K Yamashita J.Phys.Chem.Lett 2, 1167-1170 (2011)[2] R Jono, K Yamashita J.Phys.Chem.C 116, 1445-1449 (2012)

2E08
電荷移動状態を経由した励起エネルギー移動に対する計算手法の開発(神戸大院システム情報 1) 藤本和宏 1
Transition-Density-Fragment Interaction Combined with Transfer Integral Approach for
Excitation-Energy Transfer via Charge-Transfer States(Kobe Univ.1)

Kazuhiro J. Fujimoto1【序】 励起エネルギー移動(EET)は分子ペアや分子集合体において観測される現象であり、その特性は植物の光合成に代表されるように生体内での反応において広く用いられている。生体内での EET の本質的理解は、高効率の太陽電池の開発に対しても役立つと考えられている。
EET 機構を解明するために、これまでに数多くの理論的研究が行われてきた。代表的なものとしては、Förster 機構と Dexter 機構が挙げられる。
Förster 機構とは分子間の擬クーロン相互作用が起源となる EET であり、Dexter 機構とは分子間における電子スピンの交換によって生じる EET である。
これら二つの機構は EET 研究の中核を担ってきたが、1994 年に Harcourt らは電荷移動状態を経由した EET という新たな機構を提案した[1]。Fig. 1 に示すように、Harcourt らの提案した EET は、ドナー分子が励起状態となる始状態とアクセプター分子が励起状態となる終状態以外に、各分子が電荷Fig. 1. Schematic illustration of the EET via CT states.
Blue and red lines show the HOMO and LUMO of themonomer (I or J), respectively.

分離した中間状態を考慮したものである。Förster 機構を取り扱うための高精度な計算手法が近年開発されるようになってきたのに対し[2, 3]、取扱いが複雑な Harcourt らの機構に関する研究は十分に行われていないのが現状である。
【結果】 本研究では、Transition-Density-Fragment Interaction (TDFI)法と Transfer Integral (TI)法を組み合わせた方法を提案し、この方法が電荷移動状態を経由した EET に対して有効であることを示す[4]。TDFI 法は擬クーロン相互作用を記述するために開発された手法であり、これまでの研究においてキサントロドプシンにおける EET[2]やエキシトン結合型の CD スペクトル[3]に対して適用されてきた。そして、従来の研究で多用されてきた双極子・双極子近似や遷移電荷を使用した方法(TrESP 法)と比較して、TDFI 法は劇的に精度を向上させることを示してきた。また、TDFI
法では分子間の self consistency を満たすために DFI 法[5]を採用しているが、この結果、電子的カップリングエネルギーの計算精度が高くなることも示してきた[3]。今回の研究では TDFI 法に二つの改良を加えた。まず、電子のスピン交換相互作用を含めた形に拡張を行った。これにより、

Förster 型の EET 以外に Dexter 型の EET への適用を可能とした。二番目の拡張は、フラグメント分割した遷移密度行列間の重なりに対する補正を加えた。さらに、これら二つの拡張を加えた
TDFI 法に TI 法を組み合わせることを試みた(TDFI-TI 法)。本研究では相互作用に対する 3 次の摂動までを考慮した以下の形で電子的カップリングエネルギーの評価を行った。
TDFI
TDFI
TDFI
TIF  TDirect  TIndirect  VCoul
 VExch VOvlp  V TI(2)  V TI(3) 

TDFI
TDFI
TDFI
、 VExch、 VOvlpは擬クーロン相互作用、スピン交換相互作用、重なり補正項をそれ
ここで、 VCoul

ぞれ示し、各項は TDFI 法によって求めた。また、 V TI(2) と V TI(3) は中間状態を経由した相互作用(2次と 3 次の摂動項)を示し、これらは TI 法によって求めた。
TDFI-TI 法の精度を確認するために、
エチレン２量体を用いたテスト計算を行った。その結果、TDFI-TI 法は以前の
擬クーロン相互作用のみしか考慮していなかった TDFI 法と比べて格段に精度が向上することが確認できた(Table I)。
この結果に基づいて、電子的カップリン
グエネルギーの分割解析を行った(Fig.
2)。その結果、エチレン２量体の距離が3.5 Å 以上の場合では擬クーロン相互作用の寄与が大きく、また距離が狭い場

Fig. 2. Electronic coupling energies between two ethylene molecules
as a function of intermolecular distance.

合には電荷移動状態の関与する間接的な相互作用の寄与が大きいことが分かった。さらに、交換相互作用による寄与は非常に小さいことが分かった。
本発表では、その他のテスト計算の結果、電荷移動状態と分子軌道の関係、さらに Dexter 機
構についても述べる予定である。
Electronic coupling energies of ethylene dimer (cm-1).

Table I.

CIS

TD-rCAM-B3LYP

Full QM

1976

1454

TDFI-TI

1997

1417

1649

1097

TDFI
a

a

Previous TDFI method [3].

References:[1] R. D. Harcourt, G. D. Scholes, and K. P. Ghiggino, J. Chem. Phys. 101, 10521 (1994).
[2] K. J. Fujimoto and S. Hayashi, J. Am. Chem. Soc. 131, 14152 (2009).
[3] K. J. Fujimoto, J. Chem. Phys. 133, 124101 (2010).
[4] K. J. Fujimoto, J. Chem. Phys. 137, 034101 (2012).
[5] K. Fujimoto and W.-T. Yang, J. Chem. Phys. 129, 054102 (2008).

2E09
Theoretical Investigation of Proton Conductivity
in Proton Exchange Membranes for Fuel Cells(Fuel Cell Nanomaterials Centre, University of Yamanashi) Raman K. Singh, Takao Tsuneda,
Kenji Miyatake, and Masahiro Watanabe
Introduction: Nafion (DuPont) is the most used proton exchange membrane in polymericelectrolyte membrane fuel cells (PEMFCs) because of its highly proton conductive properties
and excellent chemical and physical stabilities. However, some potential drawbacks such as
high production cost, high methanol permeability and poor mechanical stability at hightemperatures of these polymers have prompted for development of novel membranes.1 Inquest of alternatives to perfluorinated polymers membranes for automotive applications,multiblock sulfonated poly (arylene ether sulfone)s (SPE) (Fig. 1) containing highly
sulfonated hydrophilic component seem to be very
promising.1 The multiblock SPE membranes have shown
more developed phase separation between hydrophobic
and hydrophilic blocks than that of the random
copolymer equivalents. It has been observed that the
proton conductivity of SPE-bl-2 membrane is higher thanNafion (NRE212) over a wide range of humidity and athigh temperature (>100 oC) but fails at low humidity. The
reasons behind the poor proton conductivities of these
SPE membranes at low humidity are still unknown andunexplored either experimentally or theoretically, and Figure 1: The optimized structure
deserve utmost attention for the future membrane of SPE membrane, in whichdevelopment. This study intends to figure out why Nafion green and blue dotted circles
indicate the separated and paired
shows higher proton conductivity than those of othersulfonic acid groups.
hydrocarbon membranes at low humidity.
Computational Details: Geometry optimizations and frequency calculations have been
carried out by the generalized Kohn-Sham method using long-range corrected Becke 1988
exchange+ Lee-Yang-Parr correlation (LC-BLYP) functional in conjunction with cc-pVDZ
basis set without imposing any symmetry restriction and all structures have positive real
frequencies. The computations have been performed using GAUSSIAN 09 program.
Results and Discussion: To discuss the state of protons in the hydrated sulfonic acid groups(–SO3H), we firstly optimized the geometries of Nafion and SPE models with various numberof water molecules, i.e. –SO3H-(H2O)n for n = 0, 1, 2, …, 10. It was found that the optimized
distance between the sulfonic acid group and proton increases sharply from 1.032 to 3.084 Å
and 1.037 to 2.084 Å for one to three hydration water molecules in Nafion and SPE,
respectively, while the distance between the sulfonic acid group and the closest watermolecule are nearly constant up to
eight hydration water molecules forboth the membranes. This clearly
indicates that minimum three watermolecules are needed to detachproton from the sulfonic acid of
both the membranes. For SPE,
however, we found that the sulfonicacid groups hydrated in pairs(“paired” in Fig. 1) detach no proton
Figure 2: The optimized hydrated structures of Nafion
even with 12 hydration water (left) and SPE (right) membrane models, for which themolecules (i.e. 6 molecules per detachment of proton (red circle) requires threegroup). The presence of such hydration water molecules in vacuo.
inactive sulfonic acid groups is
supported by an IR spectrum study. In addition, we also found that the proton is always
located at the edge of the hydration water cluster after the dissociation as supported by a study
on ab initio molecular dynamics simulations and isotopic-exchange spectroscopy.2
Furthermore, we analysed the molecular orbitals of hydrated Nafion and SPE membranes todiscuss the proton detachments from the viewpoint of electronic structure. As results, it wasnoticed that the electron distribution of HOMO drastically changes from the backbone of thechain to the sulfonic acid groups in Nafion as the hydration level increases. This seems to
make it easy to dissociate protons and protonated water clusters from the sulfonic acid groups.
On the other hand, the HOMO and LUMO of SPE are hardly affected by the hydration.
To investigate the proton conductivityof Nafion and SPE at low humidity infurther detail, we calculated thedissociation energies of protonatedwater clusters from the sulfonic acidgroups. As a result, we found that theserequire more than 100 kcal mol-1 invacua. This indicates that protonatedwater clusters hardly detach from the Figure 3: The optimized structure of double Nafion
sulfonic acid groups at low humidity. models hydrated in pairs with 5 hydration waterWe, therefore, considered that the molecules. The distance of the sulphur atoms is
proton conductance proceeds through fixed at the averaged distance given in a previous
the relay. It was found that the MD calculation.
dissociation energies drastically decrease to around 40 kcal mol-1 through the relay. We arenow exploring the feasibility of this relay.
References:
1. Miyatake, K.; Bae, B.; Watanabe, M. Polym. Chem. 2011, 2, 1919-1929.
2. Buch, V. et al., Proc.. Nat. Acad. Sci. 2007, 104, 7342-7347.