Microsoft PowerPoint pptx

Size: px

Start display at page:

Download "Microsoft PowerPoint pptx"

きみえうとだ
5 years ago
Views:

1 平成 23 年 5 月 25 日分子モデリングと分子シミュレーションポテンシャルエネルギー東京大学大学院農学生命科学研究科アグリバイオインフォマティクス教育研究プログラム寺田透

2 講義予定 1. 5 月 25 日 ( 水 ) ポテンシャルエネルギー 2. 6 月 01 日 ( 水 ) 分子動力学法とモンテカルロ法 3. 6 月 08 日 ( 水 ) 分子動力学法の応用 4. 6 月 15 日 ( 水 ) 複合体構造モデリング参考図書 : 岡崎進コンピューターシミュレーションの基礎化学同人

3 本日の講義内容ポテンシャルエネルギーとは分子軌道法実習分子力学法課題 1 力場パラメータの決定エネルギー最小化分子力学法実習課題 2

4 立体構造とエネルギー物体に力をかけて変形させると物体の持つポテンシャルエネルギーが大きくなる力同様に分子も変形するとその分子が持つポテンシャルエネルギーが変化する

5 分子のポテンシャルエネルギー分子のポテンシャルエネルギーは Schrödinger 方程式を分子軌道法を用いて近似的に解くことで計算できる分子軌道 : 分子における 1 電子の波動関数分子軌道は分子を構成する各原子を中心とした原子軌道の和で表される原子軌道は基底関数を用いて表される

6 参考 :Schrödinger 方程式 Hˆ Hˆ E elec elec = = Φ = E elec N A i + i= 1 E elec M Φ M + + A= 1 B> A N M Z r i= 1 A= 1 ia Z A Z B RAB N N i= 1 j> i 1 r N: 電子数 M: 原子数 Z A : 原子 A の原子番号 ij 分子のポテンシャルエネルギー

7 分子軌道法実習 (1) 本実習では量子化学計算ソフトウェア Gaussian 09W を用いるデスクトップにあるアイコンをダブルクリックしてこのソフトウェアのグラフィックユーザーインターフェイスGaussView 5.0を起動 Current Fragment Control Panel Molecule View Window

8 分子軌道法実習 (2) 1. Control Panel の Ring Fragment をクリックし Current Fragmentが benzeneになっていることを確認して Molecule View Window の中を左クリック 2. Control Panelのメニューから Calculate Gaussian Calculation Setup を選択 3. Job typeを Energy MethodのBasis setを6-31g(d) に設定し Submit 4. インプットファイルを保存するか聞かれるので Save しデスクトップに benzene.gjf として保存する 5. Run Gaussianウィンドウが出るので OK

9 分子軌道法実習 (3) 6. 計算が終わったら Gaussian windowを閉じるか聞かれるのではい 7. Gaussian Job Completedウィンドウでは benzene.log ファイルを選択し OK 8. Control Panelのメニューから Results Summary を選択 E(RHF) の欄に分子のポテンシャルエネルギー * が表示されている * 1 a.u. = kcal/mol

10 用語の解説 Method:Schrödinger Shödi 方程式を解くのに用いる近似法 Hartree-Fock:ab initio 法における基本 Semi-empirical: ii 半経験的 ( 大規模な系を小さな計算コストで扱えるが結果の信頼性は低い ) DFT: 密度汎関数法 ( 電子相関の効果を比較的小さな計算コストで取り込むことができる ) Basis Set: 各原子におく基底関数系 STO-3G 3-21G 6-31G 6-311G の順により複雑な軌道が表現できる必要に応じて diffuse 関数 (+ ++) 分極関数 ((d) (d,p) など ) を加える

11 分子を変形してみよう Control Panel の Modify Bond をクリックした後ベクした後ンゼンの隣り合う2つの原子をクリック次のように結合長を変える同様にエネルギーを計算してみよう ( ファイルはデスクトップにbenzene2.gjfとして保存 )

12 エネルギー最小化分子を変形するとエネルギーが大きくなるエネルギーが小さくなるように変形させることでもとの構造に戻してみよう Job Type を Optimization に設定し Submit ( ファイルはデスクトップにbenzene3.gjf として保存 )

13 分子力学法量子化学計算は計算コストが大きい電子数の4 乗に比例して計算コストが増大分子の変形によるエネルギーの変化を小さな計算コストで予測できるようにしたい分子の変形を 1 結合長 (r b ) 2 結合角 (θ a ) 3 二面角 (ϕ d ) 4 共有結合を介さない原子間距離 (r ij ) の変化で表し全体のエネルギーの変化をそれぞれの変化からの寄与の和で表す E total = E bond + E angle + E dihedral + E nonbond

14 ポテンシャルエネルギー関数構造変化の成分ごとにエネルギーの変化をの変化を関数を用いて再現する例えば結合長のエネルギー関数は分子の結合長を変化させた時のエネルギーの変化を計算しこれを再現できるように決めるポテンシャルエネルギー関数の例 1 結合長 2 結合角 3 二面角 E ( ) = k ( r r ) + k ( ) + k [ + ( n )] r θ θ 1 cos φ δ + b ij b b b σ ij 4ε ij r ij 12 a a σ ij r ij a 6 4 非共有結合相互作用 a + q q i r ij j d d d d d

15 結合長 (1) 1. 講義のページから H2.gjf をダウンロードしデスクトップに保存する 2. スタートメニューからすべてのプログラム Gaussian 09W Gaussian 09W を選択し起動 3. メニューの File Open で H2.gjfを開く 4. Existing Job Edit ウィンドウが現れるのでこのメニューから File Exit & Run 5. Output File 名を聞かれるのでデスクトップに H2.outとして保存する

16 結合長 (2) 6. 計算が終了したらメニューから File Exit 7. GaussView 5.0を起動し Control Panelのメニューの File Open でH2.outを開く ( ファイルの種類を Gaussian Output Files (*.out *.log) にする ) 8. Control Panelのメニューの Results Scan を開く 9. Scan plot ウィンドウの内部を右クリックしウ Save Dataを選択しデスクトップに H2_ scan.txtとして保存

17 結合長 (3) [a.u.] Energy [ H 2 の共有結合長の変化に伴うエネルギーの変化 y = x x x x Bond length [Å] H2_scan.txt をExcelで開きグラフを書くと 4 次関数でよく近似できることがわかる構造変化に伴うエネルギー変化をあらかじめモデル化しておくことで低い計算コストでエネルギーを求めることができる

18 エネルギー関数の近似 [a.u.] Energy [ H 2 の共有結合長の変化に伴うエネルギーの変化 y = x x x x r E E Bond length [Å] ( r + Δr) = E( r ) 0 0 E + r ( r + Δr) E( r ) k( r r ) 実際にはエネルギー最小状態のまわりを熱ゆらぎしている生体分子のシミュレーションで扱う温度は300 K 程度 kt = 0.6 kcal mol 1 = 10 3 a.u. 実際の結合長の変化は 0.04 Å 程度 2 ( r) 1 E( r) 2 ( 3 Δr + Δr + O Δr ) r r 2 r 0 Δr が小さければ 2 次式で近似できる

19 参考 : テイラー展開関数 f(x) の x = p のまわりでの展開 x = p からの微小な変位を Δx とおく f ( p + Δx) = f ( p) 1 + L+ k! = f ( p) df + dx k d f dx df + dx 2 ( x) 1 d f ( x) 2 x= p ( x) k Δx + x = p Δx k 2! dx + L 2 x= p Δx 2 ( x) 1 d f ( x) 2 ( 3 Δx + Δx + O Δx ) 2! dx 2 x = p x= p 2 次までの展開 + 誤差項 (Δx 3 のオーダー )

20 運動の比較 E( ( ) ( ) 2 E r = k r r0 E () r F( r) = = 2k ( r r ) r μ & r = 2k r ( r ) 2k ω =, μ = μ 0 m1m2 m + m () t = Acos( ω t + δ ) r0 r 黒実線 :abb initoit MD 赤破線 : 古典近似解析解ポテンシャルエネルギーを量子化学計算を用いて求める代わりにエネルギー関数を用いて近似してもほぼ同じ軌道が得られる

結合角のエネルギー関数 -75.998 H 2 O の結合角の変化に伴うエネルギーの変化モデル系 : 水分子 -75.999-76 [a.u.

006 90 95 100 105 110 115 120 E E Angle [degree] ( θ Δθ ) = E( θ ) 0 Gaussian

21 結合角のエネルギー関数 H 2 O の結合角の変化に伴うエネルギーの変化モデル系 : 水分子 [a.u.] Energy [ θ 0 θ E E Angle [degree] ( θ Δθ ) = E( θ ) 0 Gaussian job file: H2O.gjf 2 ( θ ) 1 E( θ ) 2 ( 3 Δθ + Δθ + Δθ ) E O 2 θ 2 θ ( θ + Δθ ) E( θ ) k( θ θ ) θ 0 0 θ 0

分子間相互作用の計算 (1) 1. GaussView 50 5.0 を起動し Control Panel の Element Fragmentをクリックして O を選択 Current Fragment がOxygen Tetravalentになっていることを確認して Molecule View Window の中を左クリックク 2.

22 分子間相互作用の計算 (1) 1. GaussView を起動し Control Panel の Element Fragmentをクリックして O を選択 Current Fragment がOxygen Tetravalentになっていることを確認して Molecule View Window の中を左クリックク 2. ややはなれた別の位置を左クリック 3. Control Panel の Element Fragment をクリックして C を選択 Current FragmentがCarbon Tetrahedralになっていることを確認して Molecule l View Window の H 2 Oの水素原子から1つ選びクリック CH 3基に置換 4. 同様にもう1つのH 2 O 分子の水素原子をCH 3 基に置換

23 分子間相互作用の計算 (2) 5. Control Panel メニューの Calculate l Gaussian Calculation Setup を開き Job typeを Energy MethodのBasis set を 6-31G(d) に設定し Submit ( ファイルはデスクトップにmethane2.gjfとして保存 ) 6. 同様に CH 3 OH 1 分子についてもエネルギーを計算 ( ファイルはデスクトップにmethane1.gjfとして保存 ) 7. 以下の式を用いて相互作用エネルギーを求めるを求める AB ( E E ) ΔE = E + A B ΔE = [2 ( )] = au = 1 15kcalmol a.u kcal mol

24 分子間相互作用の成分 van der Waals 相互作用斥力 +δq δq +δq 引力 δq 静電相互作用

25 静電相互作用エネルギー関数量子化学計算では電子密度 ρ(r) が計算される電子密度 ρ(r) から分子のまわりの静電ポテンシャルϕ(r) が計算できる静電ポテンシャルを再現するように原子に部分電荷をおく静電相互作用を部分電荷どうしの相互作用として計算する E ij q q i j δq = +δq R R i j

26 van der Waals 相互作用 MP2/cc-pVTZ HF/cc-pVTZ メタンメタン分子間に引力が働いている r 無極性分子間に働く引力 ( 分散力 ) は電子相関を考慮した高精度な量子化学計算によって初めて現れる

27 van der Waals 引力の起源無極性分子であっても電子雲のゆらぎによって瞬間的に双極子が生じ近くの分子に双極子を誘起する双極子誘起双極子相互作用は常に引力引力項の起源 <μ>=0 <E(r)>=0

28 van der Waals エネルギー関数 r r=σ で E =0 r=σ で E vdw =0 r=2 1/6 σ で E vdw = ε ( 最小値 ) () σ EvdW r = 4ε r 12 σ r 6 斥力項引力項

29 二面角のエネルギー関数 (1) モデル系 : ブタンブタンの二面角の変化に伴うエネルギーの変化 φ Energy [a.u.] Gaussian job file: butane.gjf Dihedral angle [degree]

] -157.275-157.28-157.285-157.29 van der Waals -157.295-157.

30 二面角のエネルギー関数 (2) 周期関数の重ね合わせと 1 番目と4 番目の原子間の van der Waals 相互作用の和で表される E ( ) [ + ( )] 0 φ k 1 cos n φ δ = i i i i σ + 4ε r σ r Energy [a.u.] van der Waals 周期関数 Dihedral angle [degree] 1 番目と 4 番目の原子間のみ van der Waals 相互作用を考慮

31 課題 1 ブタンの二面角の変化に伴うエネルギーの変化を表すポテンシャルエネルギー関数のパラメータk n δ ε σを決定せよ講義のページにある butane_scan.xlsx をダウンロードして利用すること 1 番目と 4 番目の原子間距離は Gauss View 5.0 のScan Plotウィンドウのメニューの Plot Plot Molecular l Property を開き Bond 1 4 を指定して求めること

32 力場パラメータの決定力場パラメータとは? ポテンシャル関数で用いられるパラメータ ( 平衡結合長ばね定数部分電荷など ) 非経験的パラメータ量子化学計算の結果からパラメータを求める経験的パラメータ構造や熱力学量などの実験値を再現するようにパラメータを決める

33 問題点と解決法 (1) タンパク質や核酸など生体高分子は多数の原子からなる生体高分子全体について量子化学計算を行うのは困難同じアミノ酸同じヌクレオチドなど同じ構成単位には同じ力場パラメータを使うさらに異なる構成単位の間でも化学的に類似した環境にある原子は同じ原子種とみなし ( 点電荷を除いて ) 同じ力場パラメータを割り当てる同じ原子種を含む小さなモデル化合物についてパラメータを決定する

34 原子種とモデル化合物の例 H H1 O O HC HA HC N CT C HC C CT HA CA HA CT N HC CA CA HC CT HC HC HC H CA CA HA CA HA HA CA HA CA CA H1 HA CA CA H1 CT OH HA C HA H1 HO OH HO モデル化合物 Jorgensen & Tirado-Rives, J. Am. Chem. Soc. 110, 1657 (1988)

35 問題点と解決法 (2) 凝縮相 ( 液相など ) では分子が接近しているため第 3の分子の位置が 2つの分子の相互作用に影響を与える気相で決めたポテンシャルエネルギー関数をそのまま適用できない気相凝縮相

E 有効ポテンシャルエネルギー ( r, r, r ) = E ( r, r ) +

2 3 = E 1 eff 2 1 3 ( eff eff r, r ) + E (

36 E 有効ポテンシャルエネルギー ( r, r, r ) = E ( r, r ) + E ( r, r ) + E ( r, r ) + Δ E ( r, r, r ) = E 1 eff ( eff eff r, r ) + E ( r, r ) + E ( r, r ) , E(r 1,r 2,r 3 ):3 分子系の相互作用エネルギー E(r 1,rr 2 ):22 分子系の相互作用エネルギー E eff (r 1,r 2 ): 有効 2 体間相互作用エネルギー厳密な多体間相互作用の効果を 2 体間相互作用のエネルギー関数に取り込む凝縮相このエネルギー関数のパラメータを実験値を再現するように決める 1 2 3

水分子のモデル (1) r(oh) SPC 10 1.0 TIP3P 0.9572 HOH 109.47 104.52 r* 1.7766 1.

41 0.417 q H r(oh) ( ) [Å], HOH [degree] r* [Å 12 kcal mol -1 ], 2 r* =

37 水分子のモデル (1) r(oh) SPC TIP3P HOH r* ε q q H r(oh) ( ) [Å], HOH [degree] r* [Å 12 kcal mol -1 ], 2 r* = 6 2σ ε [Å 6 kcal mol -1 ] r* van der Waals 相互作用は酸素原子間のみ計算する q O = 2q H Jorgensen et al. J. Chem. Phys. 79, 926 (1983)

38 水分子のモデル (2) 密度蒸発熱定圧比熱膨張率圧縮率 SPC TIP3P 実験値密度 [g cm 3 ] 蒸発熱[kcal mol 1 ] 定圧比熱 [cal mol 1 K] 膨張率 [10 5 K 1 ] 圧縮率[10 6 atm 1 ] いずれも 25 1 atmにおける値 Jorgensen et al. J. Chem. Phys. 79, 926 (1983)

39 生体高分子の力場パラメータポテンシャルエネルギー関数のパラメータ ( 力場パラメータ ) は分子シミュレーションのソフトウェアと共に配布されている AMBER / CHARMM h / GROMOS, GROMACS th h/ /

40 エネルギー最小化 (1) 立体構造 ( 座標 ) を変化させてエネルギー関数の値が最小になるようにすること立体構造最適化とも呼ばれる分子動力学シシを行う際には原分子動力学シミュレーションを行う際には原子同士のぶつかりを排除したり構造のゆがみを正すために最初に必ず行う

41 エネルギー最小化 (2) 1 次のアルゴリズム Steepest descent( 最急降下 ) 法最も単純収束までに多段階を要することがある Conjugate gradient( 共役勾配 ) 法エネルギー関数が N 次元の2 次形式で近似できる場合 N 回の操作で極小に到達する 2 次のアルゴリズム Newton-Raphson 法収束は早いが Hessian( 2 E) の計算に膨大な時間がかかる

42 Steepest descent 法 1. 初期構造をr 0 とする 2. 現在の構造 r i における勾配 g i = E(r i ) を計算する 3. 勾配方向にエネルギーが最小になる構造を求め r i+1 とする 4. 上記 2 3をエネルギーがこれ以上小さくならなくなるまで繰り返す f ( x, y) x ( x, y ) = ( = x 10 ) 2y2 y 2 の最小化 f +

43 Newton-Raphson 法以下に従って E(r+Δr) ) が 0 となる ( すなわち E が極小値をとる ) 変位 Δrを求める 1 ( ) ( ) ( ) 2 E r + Δr = E r + E r Δr + Δr E( r) Δr + L 2 2 E r + Δr = E r + E r Δr + L = 0 Δr Δr ( ) ( ) ( ) [ ( )] 2 1 E r E ( r ) エネルギーがこれ以上小さくならなくなるまで繰返す収束が非常に早いので次元が低い場合に有用

44 エネルギー最小化法の問題点いずれもエネルギー関数が初期構造極小値をとる立体構造を探索する手法生体高分子には多数のエネルギー極小構造が存在するためエネルギー最小構造の探索は困難分子動力学法モンテカルロ法を利用したsimulated annealing 法が使われることが多いエネルギー最小化アルゴリズムを用いて到達可能な構造エネルギー最小構造

45 Discovery Studio 3.0 Client の起動 Discovery Studio 3.0 Clientのアイコンをダブルクリックして起動メニューの View Explores Tools を選択し Tools タブを表示させておくメニューの View Toolbars Sketching を選択し Sketching tool barを表示させておく

46 低分子化合物の生成 (1) 1. メニューの File New Molecule l Window を選択し新しいウィンドウを開く 2. Sketchingツールバーから Ring を左クリック Molecule Windowの真ん中付近を左クリック 3. ViewツールバーのRotate Sketching t モードを解除を左クリック 4. メニューの View Hierarchy を選択し Hierarchy Window を表示 5. 炭素原子を1つ左クリックし選択黄色でマークされる 6. Hierarchy Windowのツリーを展開し対応する原子がマークされていることを確認

47 低分子化合物の生成 (2) 7. Molecule l Window の中で右クリックしてメニューを出し Attributes of C5 を選択( 原子名は選択した原子によって異なる ) HybridizationがSp3になっていることを確認 8. Molecule Windowの中で何もないところを左クリックして選択を解除 9. Ctrl キーと A を同時に押して全原子を選択 10. メニューの Chemistry Bond Aromatic を選択 Attributes の Hybridization が Sp2 になっていることを確認せよ

48 低分子化合物の生成 (3) 11. Simulation ボタンを左クリック 12. Tools タブの Change Forcefield を左クリックして展開する 13. Forcefieldを CHARMm Partial Chargesを Momany-Rone とし Apply Forcefield を左クリック 14. 水素原子が付加され Forcefield Statusが Molecule 1 typed with CHARMm となっていることを確認 15. Attributesで部分電荷や原子種がアサインされていることを確認せよ

49 参考 :Molecule Window の操作 (1) Molecules Windowの中を左クリックしてアクティブにしてから以下の操作を行う回転をクリックしてからクしてから 3D Window の中で左ドラッグ並進をクリックしてから 3D Window の中で左ドラッグズームをクリックしてから 3D Window の中で左ドラッグ

50 参考 :Molecule Window の操作 (2) 選択 Molecule Window 上で原子をクリックその原子が選択され黄色い四角でマークされる Molecule Window 上で原子をダブルクリックその原子を含む残基が選択されマークされる Hierarchy Window でもチェイン残基原子グループ (backboneなど) 単位で選択できる何もないところをクリックすると選択を解除できる Hierarchy WindowではCtrlキーを Molecule WindowではShiftキーを押しながらクリックすると複数選択ができる属性 (attribute) t 選択した後右クリックで表示できる

51 参考 :Molecule Window の操作 (3) Home 最初の向き位置に戻す Fit to Screen ( 選択した ) 構造を Window にフィットするように並進拡大縮小 Center Structure ( 選択した ) 構造の中心が Window の中心に来る ( 選択した ) 構造の中心が Window の中心に来るように並進

52 参考 : 力場パラメータ Discovery Studio には CHARMM が統合されているシミュレーションの対象に応じて力場パラメータを使い分けるのが良い CHARMm:General-purpose Momany and Rone all-atom atom forcefield that also provides automatic parameter estimation charmm22:academic all-atom forcefield used for simulating protein systems charmm27:academic all-atom forcefield used for simulating DNA and protein systems

53 エネルギー最小化 (1) 1. メニューの Window Close All で現在出ているWindowを閉じる (Saveするか聞かれるが No で良い ) 2. メニューの File New Molecule Window を選択し新しいウィンドウを開く 3. Sketchingツールを使って dibenzo-p-dioxin O を作れ ( 原子の変更はメニューのの Chemistry Element ) O

54 エネルギー最小化 (2) 4. Forcefield には CHARMm を用いる 5. Toolsタブの Run Simulation のツリーを展開し Minimization をクリック設定ウィンドウが現れるウ 6. 右下のように設定して Run 7. メニューの View Data Table でエネルギー値を確認せよ

55 SMILES による化合物の表現低分子化合物の構造の表現方法の 1つの SMILESがある新しいMolecule Windowを開きメニューの File Insert From SMILES を選択 Smiles stringに入力すると対応する分子が生成される (aromatic は小文字 ) benzene: c1ccccc1 dibenzo-p-dioxin: i c13ccccc1oc2ccccc2o3 1O 2O3 alanine: [N+][C@@H](C)C(=O)[O-]

56 課題 2 以下の分子を作成しエネルギー最小化計算をせよ Oseltamivir: インフルエンザ治療薬 ( 商品名 :Tamiflu) Force fieldはcharmmとすることエネルギー最小化後の構造の図を PNG 形式で保存すること

57 課題の提出課題 1については作成したExcelファイル課題 2についてはエネルギー最小化後の構造の図のファイルを添付してメールで寺田宛 u-tokyo.ac.jp) に送ること課題 2についてはエネルギー最小化前と後のエネルギー値を本文に記載すること件名は分子モデリング課題とし本文に氏件名は分子モデリング課題とし本文に氏名と学生証番号を必ず明記すること

Microsoft PowerPoint pptx

Microsoft PowerPoint pptx 平成 24 年 5 月 21 日分子モデリングと分子シミュレーションポテンシャルエネルギー東京大学大学院農学生命科学研究科アグリバイオインフォマティクス教育研究プログラム寺田透講義予定 1. 5 月 21 日 ( 月 ) ポテンシャルエネルギー 2. 5 月 28 日 ( 月 ) 分子動力学法とモンテカルロ法 3. 6 月 04 日 ( 月 ) 分子動力学法の応用 4. 6 月 11 日 (