4Bp01 エネルギー表示の方法による溶媒和自由エネルギーの解析（京都大学化学研究所） ○松林伸幸・中原勝 1. はじめにじめに溶液内過程を記述する最も重要な量は、その過程に対応する自由エネルギー（変化）である。特に、ある過程の始状態と終状態にある化学種の化学ポテンシャルが分かれば、その過程の自由エネルギー変化は分かるので、溶液内にある溶質の溶媒和自由エネルギー（化学ポテンシャル）を評価する方法の確立が、溶液の統計力学の最も重要な課題である。自由エネルギーの計算は、簡単な溶質については自由エネルギー摂動法などで可能であるが、より複雑なものを目指すとき、近似手法を定式化する必要がある。さらに、溶媒和自由エネルギーは、分子レベルの溶媒和の情報を反映する量であるため、溶媒和を統計力学的に記述する分布関数の汎関数として表すことが重要である。分布関数を定義するためには、対象とする分子の座標を設定することが最初のステップである。ある剛体分子の配置を完全に指定するためには、その分子の並進及び回転の自由度を表す一般に 6次元の座標を設定すればよい。しかし、6次元の座標空間上での分子配置の表現は、概念的にも実際運用上でも困難を伴う。この困難を解消する方法としてよく用いられるのが、相互作用点表示である。分子を相互作用点の集合として表し、相互作用点のペアの相対配置はその間の距離という 1 次元の変数で記述される。RISMは相互作用点表示に基づく溶液理論であり、高密度の分子性液体について有効であるが、現実的な 2 体分布のレベルにとどまる限り、低中密度領域で良い結果を与える溶媒和自由エネルギーの汎関数の定式化は困難である。ところが、超臨界水やアルコール溶液を用いるときは、実験上の困難もあって、低中（溶媒）密度領域を取り扱うことが多いので、低中密度領域にある分子性溶液についての近似理論を構築することが必要となる。本研究の目的は、超臨界状態を含む広い温度・密度範囲で、様々な種類の溶質に対して適用可能な、溶媒和自由エネルギーの汎関数を構築することである。この目的のために、溶質―溶媒相互作用エネルギーの分布関数によって溶媒和自由エネルギーを記述する、エネルギー表示の方法を開発した。本発表では、方法の概略を述べ、常温常圧および超臨界状態の水に溶けた代表的な無極性・極性・イオン性の溶質についてのモデル計算の結果を報告する。また、柔軟な構造を持つ分子系への拡張も容易であり、その結果についても報告したい。 2. 方法論方法論無限希釈条件における溶液系を取り扱う。x を、対象とする溶質分子と溶媒分子の相互作用を決定するに足る座標の組とし、v(x) を問題とする溶質―溶媒相互作用ポテンシャルエネルギーとする。一般には、xは溶質に対する溶媒の相対配置であり、分子振動などの内部自由度があるときはそれを x に含めてもよい。エネルギー表示では、溶媒の溶質に対する相対座標 ε を単純に v(x) の値とする。そのとき、対象とする溶質分子の周りの溶媒分子の瞬間的な分布 ρˆ (ε ) を ρˆ (ε ) = å δ (v (x ) − ε ) （i 上の和は全溶媒分子についての和） i で導入し、⋅ ⋅ ⋅ および ⋅ ⋅ ⋅ 0 を、それぞれ、溶液系と純溶媒系での平均としたとき（「純溶媒系」では溶質分子はテスト粒子として仮想的に挿入されている）、 ρ (ε ) = ρˆ (ε ) ρ 0 (ε ) = ρˆ (ε ) 0 χ 0 (ε ,η ) = ρˆ (ε )ρˆ (η ) 0 − ρˆ (ε ) 0 ρˆ (η ) 0 で定義される 3 つの分布関数によって、溶媒和自由エネルギー ∆µ を表す。エネルギー表示でも、密度汎関数理論（溶質−溶媒相互作用と分布関数の一対一対応）に基づく定式化を ∆µ に対して行うことができ、ρ(ε), ρ0(ε), χ0(ε,η) による ∆µ の表式は、2 体レベルの近似的記述に相当する。計算機シミュレーションから求めた ρ(ε), ρ0(ε), χ0(ε,η) を用いることとし、PY 及び HNC 近似を組み合わせて ∆µ の汎関数を近似的に構成した。これらの近似は、それぞれ、斥力的および引力的な相互作用領域で良い結果を出すことが、単純液体の動径分布関数の理論で経験的に知られている。エネルギー表示での近似は、動径分布関数理論での近似表式を、ρˆ (ε ) から生成した分布関数に対して形式的に書き直したものである。本研究で構成した ∆µ の汎関数は、溶質−溶媒相互作用および密度の 2 次まで厳密であり、低中密度領域で良い結果を与えることが期待できる。 3. モデル計算モデル計算下の表に、典型的な無極性・極性・イオン性の溶質の ∆µ について、常温常圧および超臨界状態の SPC/E モデル水中での近似値を、厳密に自由エネルギー計算した値と比較する。用いた ∆µ の汎関数は単一であり、低中密度領域から高密度領域まで、および、全タイプの溶質で、良い結果を与えていることが分かる。ただし、近似計算は、厳密な計算の数十分の一以下の計算量である。分子構造に flexibility のある系への拡張も容易である。調和振動子系および二面角系の溶質に対して、厳密な自由エネルギー計算が現実的に可能であり、近似的方法論のテストができる。その結果、上と同じ汎関数が、常温常圧から超臨界領域の水で良い結果を与えることが分かった。今後、いくつかの分子内自由度を持つ溶質系とその会合の解析を行う。 Solvation free energy ∆µ in the unit of kcal/mol.a) thermodynamic state solute 1.0 g/cm3, 25 ℃ 1.0 g/cm3, 400 ℃ 0.6 g/cm3, 400 ℃ 0.2 g/cm3, 400 ℃ methane 3.0±0.2 / 2.8 9.4±0.1 / 10.7 3.1 / 2.8 0.6 / 0.6 ethane 2.6±0.3 / 1.9 11.3±0.1 / 13.2 3.2±0.1 / 2.8 0.6 / 0.5 water -8.2±0.7 / -6.9±0.7 -1.5±0.4 / -0.7±0.3 -3.3±0.2 / -3.1±0.1 -2.4±0.2 / -2.3±0.2 methanol -4.4±0.5 / -3.9±1.3 4.5±0.4 / 6.0±0.3 -0.6±0.3 / -0.8±0.1 -1.2±0.1 / -1.3±0.1 ethanol -4.4±0.4 / -3.0±1.2 8.0±0.3 / 10.9±0.3 0.4±0.1 / 0.0±0.2 -0.9±0.1 / -1.0±0.1 a) Na+ -96.0±1.9 / -100.5±0.7 -85.9±0.9 / -90.2±0.1 -80.2±0.7 / -83.1±0.3 -67.5±0.9 / -74.5±0.6 Cl− -69.3±1.6 / -64.2±1.4 -54.0±0.9 / -51.5±0.3 -53.9±0.6 / -54.4±0.2 -45.5±0.9 / -46.7±0.4 NaCl -82.9±1.8 / -78.6±1.3 -61.4±0.9 / -59.7±0.4 -58.7±0.9 / -59.1±0.4 -47.5±1.2 / -50.5±0.4 Each ∆µ is expressed in the form of (approximate value) / (exact value). The error is given at 95% confidence level, and is smaller than 0.1 kcal/mol when not shown. 4Bp02 量子化学計算による超臨界水中における水分子のイオン解離反応の自由エネルギー計算（阪大院基礎工）○高橋英明 *、佐藤亘、新田友茂 1. 緒言水中の水分子のイオン解離は、溶液中で起こる最も基本的な反応であり、多数の分子が静電場として反応に関与するので量子化学計算や分子シミュレーションの対象としても興味深い。溶液中の化学反応を研究する上で、反応する溶質を量子化学的に扱うことは必須である。我々は、これまでに密度汎関数法( DFT )に基づく量子化学計算により溶質分子の電子状態を第一原理的に計算し、溶媒としての水分子を古典的なモデルで記述するハイブリッド型の第一原理分子動力学法( Quantum Mechanical /Molecular Mechanical ( QM/MM )法）を開発してきた[1]。この手法によれば、遷移状態や反応中間体などの分子の電荷分極を、溶液の熱力学条件に関わらず非経験的に決定することが可能である。また、QM/MM法は通常の古典分子動力学計算では考慮できない電子密度の広がりやその揺らぎなども再現できる。溶液中の反応経路を決定する最も重要な量は、言うまでもなく反応に対する自由エネルギー変化であり、第一原理的に自由エネルギーを決定することは量子化学においても重要な課題の一つとなっている。最近、Matubayasi らは、溶質ー溶媒間相互作用エネルギーの分布関数を変数とする溶媒和自由エネルギーの表式（エネルギー表示の理論）を提案し[2]、古典分子動力学シミュレーションで得られたエネルギー分布関数から様々な分子の溶媒和自由エネルギーを高い精度で計算することに成功した。エネルギー表示の理論はシミュレーションとの融合が容易であるだけでなく、RISM( Reference Interaction Site Model )理論で現れる相互作用点という概念を必要としないので、量子化学的な方法に馴染みやすい。超臨界水中の水分子のイオン解離反応は、超臨界水の持つ大きな反応性と関係があると考えられており、解離反応の自由エネルギーの密度依存性やその物理化学的背景を第一原理的に研究することは重要である。本研究ではQM/MM計算とエネルギー表示の理論をカップルさせた方法により、様々な熱力学条件下での水のイオン解離の自由エネルギーを計算する。 2 . 理論と計算方法理論と計算方法 QM/MM 法では、全系のハミルトニアン H を、以下のように表す。 H = H QM + H MM + H QM/MM 上式右辺の第 1 項、第 2 項はそれぞれ、量子化学的(QM)に取り扱う溶質分子のハミルトニアン、古典的(MM)に取り扱う溶媒水分子間の相互作用を表す項である。また、第 3 項は QM 系と MM 系の相互作用を表すハミルトニアンであり、次式で表せる。 H QM/MM = − ∑∫ s, i ρ (r) ⋅ qi dr + r − ris Zn ⋅ qi ∑∑ R n s, i n − ris + VLJ ここで、ρ(rr )は座標 r での QM 系の電子密度であり、rr is、qi はそれぞれ s 番目の水分子の i 番目 R n は、それぞれ、n 番の点電荷の座標及び水分子の i 番目のサイトの電荷である。また、Zn、R 目の原子核の電荷と位置座標を表す。この式の右辺の第1項はQM系の電子雲と古典水分子の点電荷とのクーロン相互作用を、第 2 項は QM 系の原子核と点電荷との相互作用を表し、第 3 項は QM - MM 間の Lennard-Jones ポテンシャルである。本研究では、H2O とその解離生成物である OH-、H3O+ を QM 系の分子とし、これらの電子状態を、BLYP を交換相関汎関数とする Kohn-Sham の密度汎関数法により計算した。また、255 個の溶媒水分子を MM 系の分子とし、これらを TIP4P モデルで記述した。溶液系 100 ps、純溶媒系 200 ps の QM/MM 計算の統計平均により溶質ー溶媒間相互作用エネルギーの分布関数 ρ(ε)、ρ0(ε)を求め、これらをエネルギー表示の近似的な汎関数に代入することにより活性化自由エネルギーを求めた。密度 0.1 ~ 0.6 g/cm 3、温度 600 K( Tr = 1.07 )の超臨界水における各溶質の溶媒和自由エネルギーを計算し、イオン解離の自由エネルギーの密度依存性を計算した。 3. 結果と考察図 1 に、溶液の温度一定( Tr = 1.07 )で、密度を 0.1 から 0.6 g/cm3 まで変化させたときの反応物及び生成物の溶媒和自由エネルギーの密度依存性を示した。この図より、密度の増加とともに水分子の溶媒和自由エネルギーは緩やかに減少するのに対して、イオン種のそれは急激に減少することが分かる。結果として、図 2に示すようにイオン解離の自由エネルギーも密度とともに単調に減少する。解離の自由エネルギー変化∆Gを内部エネルギー項∆Hとエントロピー項(-T∆S)に分けたところ、∆H は密度に依存せずほぼ一定であるのに対して、エントロピー項は密度に対して顕著に減少することが分かった。これは、イオン周りの溶媒水分子の凝集に伴うエントロピーの減少が、低密度の方が高密度領域に比べて極めて大きいという事実を反映している。解離自由エネルギーの密度依存性は、傾向において実験結果を良く再現しているが、エネルギーの絶対値は合わない。今後、OH- の電子密度の揺らぎの効果や、hydronium イオンが Zundel 型( H5O2+)や Eigen 型( H9O4+)を形成する可能性について考察する必要がある。 Reference: [1] H. Takahashi, T. Hori, H. Hashimoto, T. Nitta, J. Comput. Chem., 2 22, 1252( 2001 ) [2] N. Matubayasi, M. Nakahara, J. Chem. Phys., 1 1 33, 6070( 2000 ) 0 -40 130 -5 -50 -60 H3O+ -15 -70 -20 ∆G 22.6 MPa 120 ∆G [ kcal/mol ] H2O -10 ∆Gsol [ kcal/mol ] ∆Gsol [ kcal/mol ] 140 110 30.7 100 40.9 90 56.2 80 97.9 -80 -25 - OH 70 -30 0 0.1 0.2 0.3 0.4 0.5 0.6 -90 0.7 ρ [ g/cm3 ] 図 1 Tr = 1.07 での反応物及び生成物の溶媒和自由エネルギーの密度依存性 35.0 60 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 ρ [ g/cm3 ] 図 2 Tr = 1.07 での解離の自由エネルギーの密度依存性。矢印は、対応する圧力[MPa]を表す。 4Bp03 Theoretical Investigation of Structures, Thermodynamics and Reactivity of Dipeptide Molecules Based on Multidimensional Free Energy Maps (Rutgers University)○Tateki Ishida・Ronald M. Levy The structure, thermodynamics and reactivity of conformational changes of a series of dipeptides protein is investigated based on multidimensional free energy ( potential of mean force ) profiles. The free energy surfaces are generated by using umbrella sampling and the weighted histogram analysis method. The analytical generalized Born model including nonpolar hydration is employed to take solvent effects into account. Two representative torsion angles ( for example, phi and psi for alanine dipeptide ) are selected as reaction coordinates and two-dimensional free energy surfaces are estimated. As an exmaple, we find that the free energy profile of alanine dipeptide is in very good accord with high resolution experimental data extracted from the PDB ( Protein Data Bank ) of phi and psi dihedral angle distributions observed for alanine residues in proteins. Also, we estimate the existence probability around a minimum point for each conformation, and discuss its dependence on temperature considering the possibility to observe a conformation with small lifetime experimentally in single molecule experiments. On proline dipeptide, we observe that the results of the relative free energy between minima and the free energy surface profile depend on the selection of reaction coordinates. Based on this fact, the way of choosing reaction coordinates for a variety of amino acid types will be discussed. We will also discuss the dynamics of conformational change. 4Bp04 ab Initio 分子軌道法による分光学的精度を目指した分子定数の高精度予測（産総研グリッド研 a ・お茶大理 b ・筑波大数理 c） ○平野恒夫 a,b ・糸野幸子 b・小鷹恵利香 b ・三井由香里 b ・天野倫子 b ・福井玲 c ・長嶋雲兵 a 高精度な ab initio 分子軌道法による理論計算は、未知分子の理論的予測、実験サイドへの相補的な示唆などの面で重要である。もし、分子分光学実験に匹敵する精度の ab initio 分子軌道法計算が可能となれば、実験に先駆けて、観測されるであろうスペクトルのシミュレーションが出来るし、実験では得がたい情報、例えば電子状態の帰属、双極子能率、遷移双極子能率、Renner 定数などを得る事もできる。かかる分野を、我々は Computational Molecular Spectroscopy と呼ぶ。ここでは、電子相関を考慮した高精度な分子軌道法によって、いかに精度よく分子分光学定数を予測できるかに関して、MgNC/MgCN、NCS、FeC、FeS、FeN、FeCO、FeCN/FeNC、CoCO、CoH など、我々のグループで手掛けてきた例を挙げて、議論する。これらの計算は、ほとんど全ての場合、多配置の SDCI または ACPF 法で行ったものであるが、分子分光学で要求される精度を出すために工夫した点、注意すべき点について述べるとともに、core-valence 電子相関や、相対論的補正の必要性についても議論することにする。１） MgNC/MgCN 我々の分光学精度を目指した計算は、Guélin らが IRC +10216 で観測した星間分子の未同定ライン 1 を、MgNC によるものと帰属することに成功したことに始まった。 2,3 実験 3 と同時進行で進められた計算 2 では表１に示すように、回転定数 B0 を 0.04 %の誤差で予測することができた。また、ＭｇＮＣの励起状態a 2ΠのLIFスペクトルの実測4 Table 1 MgNC (X 2Σ+) ACPF/TZ2p+f Microwave exp. (core-valence)a (Kagi, et al.b) B0 /MHz DJ /MHz ω2 /cm-1 5969.3 0.0029 90.4 α2B /MHz -78.5 5966.8969 0.0042338 86 -70.2 a) Ref. 2; b) Ref. 3 と計算（MR-SDCI+Q/[TZ3P1f (Mg), aug-cc-pVQZ (N and C)]）5 の比較から、Renner定数ε は0.259（実験では0.05） MgNCの励起状態 2ΠのLIF スペクトルであること、および図１に示すように、実験でµΠ1/2, µΠ3/２と帰属されたラインはκΠ1/2, κΠ3/２の間違いであることを明らかにした。振動数の予測精度は数 cm-1であり、スピン-軌道相互作用定数ASOも 34.9 (実測 36.926) cm-1と、精度よく求まった。 Our assignments ２） NCS Renner効果でよく知られているNCOと等電子である NCSについては、B0 はよく知られているものの、結合距離については実験値がなかった。また、従来のかなり高度な多配置CIの計算でも、B0 の誤差は0.8 %もあった。ごく最近の我々の計算6でも、core-valenceの電子相関 Figure 1 MgNC a 2Π を無視する限り、誤差をこれ以上小さく出来なかったが、core-valenceの電子相関を取り込んだ大規模な計算（MR-SDCI+Q/[aug-cc-pCV(Q+d )Z (S), aug-cc-pCVQZ (N, C)]）を行ったところ、NC32SおよびNC34SのB0を誤差0.03 %で再現し、ｒe（N-C）=1.1784 Å、ｒe（C-S）=1.6320 Å と決定することができた。6 ３） Fe を含むラジカル Fe を含むラジカルは、開殻系であり、かつ 3d 軌道の擬縮退が多かれ少なかれ残っている場合が多いので、低い位置に沢山の励起状態があって、その電子状態および分子構造を分光学精度で計算することは困難であった。我々は、未知星間分子からの興味で、これまで FeC、FeN、FeS、FeCN、FeNC について ab Initio 分子軌道計算を行い、MCSCF 軌道を State-averaged MCSCF 法で正しく作り上げることによって、 MR-SDCI+Q または MR-ACPF 法の計算で、基底状態および低い励起状態の分子定数や電子構造を分光学実験に耐える精度で予測することに成功した。例として回転定数の計算値を表 2 に示す。分光学精度を出すには、相対論的な補正が必要で、その効果としては結合距離が 0.005 Å Table 2 Rotational Constants (B0) 程度短くなる。これだけの精度で計算できると、例えば、表 2 で 3 %の誤差を Exp. Unit in cm-1 Previous Calcs. Our Calc. FeN 12∆5/2 0.602793(17) 0.60280246(25) 0.60284 (0.0%) DFT DFT FeS 12∆i 0.20368 0.20246 (0.6%) MR-ACPF 0.1911 (-6.2 %) DFT 0.2011 (-1.3 %) FeC 13∆i 0.67291212(6) 0.66966 (-0.5%) MR-SDCI 0.6754 (0.4 %) MR-SDCI 0.6623 (-1.6 %) 23∆ 3 0.55321(15) 0.5521(10) 0.54955 (-0.5%) MR-SDCI 0.5298 (4.0 %) MR-SDCI 0.5388 (-2.4 %) 33∆2 0.5954 0.59285(6) 0.57498 (-3.0%) MR-SDCI 0.668 (3.7 %) 0.56442(15) 0.56153 (-0.5%) 0.5693 (-5.6 %) 0.6099 (1.2 %) 示す FeC の 3 3∆では、この状態には強い摂動が掛かっていることを示唆している、と読み取ることが出来る。双極子能率 µ の計算であるが、これらの Post Hartree Fock 計算ではエネルギー的には最適化されるものの、波動関数自体は必ずしもそうでもなく、例えば FeC の基底状態 3∆のµ の期待値としての値は 1.30 D （MR-SDCI）である Strong perturbation from nearby 1∆ 43∆2 MR-SDCI 0.5417 (-4.0 %) → The error of predicted B0 : 0.5 – 0.6 % が、エネルギーの電場微分としてµ を求めると、MR-SDCI+Q+Erel で 2.24 D7 となり、実測値 2.36 D8 とよく一致した。 3926 3889 1∆2 3894.27 スピン-軌道相互作用の計算の例として ∆E = T0(1∆2 – 3∆2)Unperturbed 3528 FeC の場合 7 を図 2 に示す。この結果は電子状態間の摂動を考える上で実験に貢献し、計算で予測した ASO = -180.7 cm-1 は、正しい実験値 -179.89cm-1 を導く切っ掛けとなった。また、FeCO に関しては、MR-SDCI+Q＋ 3 ∆1 Erel（相対論補正）、MR-ACPF+Erel 法により、初めて正しく基底状態 3Σ-を励起状態 5Σ- より低く算出するとともに、結合距離の実測値を 3 ∆ 再現することができた。9 (ASOΛ)2/∆E = 37 361 324 0 Calc. Figure 2 719.56 359.8 329.809 ASOΛ = 361 ∆3 CoN、CoF、CoCN/CoNC について計算中で論会で別報として発表する。 ∆2 3 ４） Co を含むラジカルあり、CoH に関しては今回の分子構造総合討 3 722 0 Exp. Spin-Orbit Interaction for FeC (Ref. 7) (Unit in cm-1) 1) M. Gu´elin, et al., Astron. Astrophys, 157, L17 (1986). 2) K. Ishii, T. Hirano, U. Nagashima, B.Weis, and K. Yamashita, Astrophys. J., 410, L43 (1993). 3) K. Kawaguchi, E. Kagi, et al., Astrophys. J., 406, L39 (1993). 4) R.R. Wright and T.A. Miller, J. Mol. Spectro., 194, 219 (1999). 5) T.E. Odaka, et al, J. Chem. Phys., 115, 1349 (2001). 6) T. Hirano, et al., Ohio State Univ. Internat. Sympo. on Mol. Spec., RG02 (2003). 7) S. S. Itono, et al.,, J. Chem. Phys., 115, 11213 (2001). 8) T. C. Steimle and W. L. Virgo, J. Chem. Phys., 117, 1511 (2002). 9) 平野恒夫ら, 分子構造総合討論会, 4C01 (2002) 4Bp05 電荷平衡法で用いられるパラメターの再決定（1 豊橋技科大、２科技団、３産総研グリッド）○中山尚史 1,2、後藤仁志 1、長嶋雲兵 3 【序論】分子力学法などのシミュレーション計算において、静電的な相互作用の見積もりは極めて重要な寄与をすることが知られている。従来のシミュレーション計算では、分子内あるいは系における静電的なエネルギーを、各原子上に点電荷を置きその電荷同士のクーロン相互作用として求めるのが一般的である。しかしながら各原子上に置く部分電荷を一定とする場合、構造変化による静電場の変化を記述することは不可能である。 Rappé と Goddard は、分子内の各原子上の点電荷 Original (Gaussian98) を簡易かつ定量的に計算する手法として、電荷平 10.0 衡法（Charge Equilibration, QEq）を開発した。この initio 計算で求めたものとほぼ同等であることが報告さ 8.0 Calc. 手法によって得られた電荷および双極子能率が、ab 6.0 4.0 2.0 れている[1]。 0.0 0.0 しかしながら、従来の QEq 法で用いられている 2.0 4.0 6.0 Exptl. 8.0 10.0 パラメターでは双極子能率の実験値を定量的に再現することが不可能である（図１参照、構造は図１：Rappé らのパラメターを用いた QEq 法によって求められた、双極子能率の実験値と B3LYP/6-31G**によって最適化したものを用いての比較いる）。よって得られる部分電荷の精度を向上させるためには計算に用いられているパラメターの最適化が不可欠である。そこで本研究では、QEq 法で用いられているパラメターを、双極子能率の実験値を再現するように最適化したので報告する。【理論】 QEq 法による系の全エネルギーは式１および２で、また QEq-GB 法による系の全エネルギーは、式３および４によって表される[2]。 E vac (Q1 LQN ) = ∑ E A (QA ) + ∑ ∑ QA QB J AB A E A (Q A ) = E A0 + χ A0 Q A + (1) A B>A 1 0 2 J AA Q A 2 (2) ここで QA は部分電荷、JAB は QEq における２中心電子間反発をそれぞれ表している。最適化は、式２ 0 に含まれる各原子固有のパラメターである χ A0 と J AA に対して行った。２中心電子間反発の計算において、Rappé らは各原子上に Slater 型関数を置き、それらの積分計算によって求めている。またこの Slater 型関数は、水素原子の場合にのみ軌道指数部分が電荷に依存する形式を採用している。しかしながらこの計算方法では、取り扱う系中の水素原子が増えるにしたがい繰り返し計算の回数が増加するため、生体高分子などでは計算時間が大幅に増加することが考えられる。そこで本研究では、JAB の計算に西本—又賀式[2]および大野—Klopman 式[3]を用いることとする。この２つを用いて双極子能率を計算した結果（パラメターは、Rappé らが決定したもので図１の計算と同じもの）を図２に示す。これを見ると計算結果の全体的な傾向は、電子間反発の計算式には依存していないことがわかる。 Ohno-Klopman 10.0 8.0 6.0 4.0 2.0 0.0 Calc. Calc. Nishimoto-Mataga 0.0 2.0 4.0 6.0 8.0 10.0 8.0 6.0 4.0 2.0 0.0 10.0 0.0 2.0 4.0 Exptl. 6.0 8.0 10.0 Exptl. 図２：西本—又賀式および大野—Klopman 式により２中心電子反発を計算した場合の双極子能率の計算値と実験値との比較（構造は、図１と同様に B3LYP/6-31G**によって最適化したもの）【原子タイプおよびパラメター最適化】表１：分割した原子タイプ計算に用いるパラメターは、従来のように各元素に対して最適化するのではなく、分子力場における原子タイプを参考として、表１のように分割し最適化を行った。パラメターの最適化には、最小二乗法を用いた。 Atom Atom type H -C, -N, -O, -P, -S C Sp3, sp2, aromatic, sp, etc., N Sp3, sp2 (aromatic), sp O Sp3, sp2 【結果】 Nishimoto-Mataga 最適化したパラメターを用いて、再度双極子能み）を図３に示す。図１および２の結果と比較して、計算値は大幅に改善された。また電荷分布も、分子軌道計算によって求められたものと比べて Calc. 率を計算した結果（西本—又賀式を用いた場合の 6.0 5.0 4.0 3.0 2.0 1.0 0.0 0.0 1.0 定性的な差は見られなかった。結果の詳細は、当日報告する。 2.0 3.0 4.0 5.0 6.0 Exptl. 図３：最適化したパラメターを用いて計算した双極子能率【謝辞】本研究は、科学技術振興事業団計算科学技術活用型特定研究開発事業「GRID テクノロジーを用いた創薬プラットフォームの構築」の援助を受けて行われました。【参考文献】 [1] A. K. Rappé, W. A. Goddard III, J. Phys. Chem., 95, 3358 (1991). [2] K. Nishimoto, N. Mataga, Z. Physik. Chem. Neue Folge, 12, 335 (1957). [3] K. Ohno, Theret. Chim. Acta, 2, 219 (1964); G. Klopman, J. Am. Chem. Soc., 86, 4550 (1964). 4Bp06 Theory of Electron Transfer from the Excited State of an Adsorbed Dye to the Conduction Band of a Semiconductor K.L. Sebastian1,2 and M. Tachiya2 1 Indian Institute of Science 2 National Institute for Advanced Industrial Science and Technology We consider the process of photoinduced heterogenous electron transfer [1-5]. A molecule, attached to a semiconductor surface, is excited to an upper electronic state by a laser pulse. The electron in the upper excited orbital can tunnel into the conduction band of the semiconductor, leaving the molecule ionized and perhaps vibrationally excited. Experimentally, electron injection times are of the order of tens to hundreds of femtoseconds. One can produce the molecule in a vibrationally coherent state and study its evolution, after electron injection. Such vibrational wave packet states, were found to continue their motion for several hundred femtoseconds in the product state, implying that electron donation occurs from a non-relaxed vibrational population of the donor. We first study of the purely electronic problem, using an Anderson like Hamiltonian, which takes only the excited orbital of the adsorbate into account. The solid is represented as a single band of finite bandwidth. Using a Green’s function approach, we calculate the exact survival probability of the excited state, for arbitrary bandwidths. We show that that the existence of a split-off state can play an important role in the total injection probability. In the wide band limit, the survival probability decays exponentially, but for finite band widths it does not. We also introduce a “pole approximation”, which is easy to use, and leads to good approximate results for the adsorbate density of states and survival probability. It consists of looking for the poles of an appropriately analytically continued, adsorbate Green’s function and using them alone, to calculate the quantities of interest. In the wide band limit the method is exact, but for finite bandwidths too, it leads to good agreement with exact calculations. We also study the influence of vibrational motion, on the electron transfer. As the electron hops into the solid, the molecule is ionized, causing a change in the vibrational Hamiltonian. We use the simplest model that includes this effect too (see Ramakrishna et. al. [1]). Further the model is of great theoretical interest as it is the simplest Hamiltonian, that has a continuum of electronic states, coupled to vibrational states. We show how one can eliminate the electronic states exactly and reduce the problem to one of just vibrational motion. Any calculation that needs to be performed can, in principle proceed as if it is just for vibrational motion of the molecule. However, the vibrational motion is now governed by a non-Hermitian Hamiltonian and the calculation, naturally, is more involved. The method works for arbitrary density of states for the conduction band, into which electron injection occurs and can be used to calculate the survival probability in any given vibrational state. We show that in this case too, the quantities of interest are all contained in an “adsorbate Green’s operator”, which acts on the space of vibrational wave functions. Its analytic continuation in the complex plane has poles and one can again use the pole approximation. In the wide band limit, this is exact, and leads to the interesting result that the survival probability is independent of vibrational motion in this limit, in agreement with the conclusions of [1]. Further, we show that in this limit, the vibrational dynamics is particularly simple. All that one has to do is, to assume that each one of the vibrational state has a finite lifetime, which is just the lifetime of the electron in the adsorbate state. In the case of finite bandwidth, one has to look for the poles of the Green’s operator and use them in the calculations. We report results of calculations using this approximate method too. [1] S. Ramakrishna, F. Willig and V. May, J. Chem. Phys. 115, 2743 (2001). [2] W. Stier and O. Prezhdo, J. Phys. Chem. B 106, 8047 (2002). [3] A. Persson, M. Ratner and H. Karlsson, J. Phys. Chem. B 104, 8498 (2000). [4] S. Ramakrishna and F. Willig, J. Phys. Chem. B 104, 68 (2000). [5] S. Ramakrishna, F. Willig and V. May, Phys. Rev. B 62, R16330 (2000).