振動解析で何がわかるの?
振動数計算を行うことで、主に以下の3つの情報が得られます。
1) 計算で得られた構造が確かな安定構造か、又は遷移構造か
構造最適化計算では計算が終了したからと言って必ず自分の求めた安定構造であるとは限りません。例えば、ソフトが計算終了の閾値に達したために計算を止めてしまうこともあれば、求めようと思った安定構造とは違う局所的な安定構造が出力されることもあります。振動解析をおこなうことで、振動数がすべて実数であれば安定構造、虚数振動が一つあるときは遷移状態の構造であることが判断できます。
2) 熱力学的諸量 (エンタルピー、エントロピー、Gibbs自由エネルギーetc.)
振動数計算をおこなうことでエンタルピーやエントロピー、ギブス自由エネルギーなど熱力学的諸量が得られます。これらを利用して、さらに反応速度や反応熱、平衡定数などの計算ができます。これで、反応が進行しやすい・しにくいなどを定量的に議論できるようになります。Outputファイルに出力される熱力学的諸量の見方はこちらを参照してください。
3) IR(赤外吸収分光法)やラマンスペクトル
HF法(ハートリック・フォック法)やDFT法(密度汎関数理論)により得られる基準振動とその振動数は、赤外線やラマン分光装置で計測される値と直接比較することができます。量子化学計算で得られるIRスペクトルは、試料測定時のノイズを除いた振動に由来する純粋な吸収ピークを示し、化合物の同定や各物質間で強度の相対比較を行うことができます。
計算の準備をしよう
振動解析を行うことで、上記の3つの情報が一度に得られます。今回は、IRスペクトルをGUIソフトで表示させ、実測値との比較を行ってみましょう。計算にはFirefly(PC GAMESS)を用います。振動数計算は構造最適化計算と異なり計算負荷が大きいので気長に待つ必要があります。それでは早速、計算に必要な入力ファイルを作成しましょう。
1) 今回は、アセトフェノンを計算してみましょう。MolView又は他の分子モデリングソフトでアセトフェノンのMOLファイル(拡張子.mol)を作成してください。
2) 保存したMOLファイルに分かりやすいようにAcetophenonと名前を付けて、MoCalc2012に読み込みます。計算条件はHF/6-31G(d)(理論/基底)に指定してください。赤枠①『Method』の部分を「6-31G*」に変えるだけでOKです。
(一応「Job Type」が「Geometry Optimaization(構造最適化)」になっていることを確認してください。振動数を計算するには、まず対象とする化合物の安定構造を得るために構造最適化計算が必要です。)
3) 次に「More Keyworda」の所に「HSSEND=.TRUE.」と打ち込んでください(赤枠②)。このキーワードを付け加えることで、構造最適化計算後に振動計算を連続実行することができます。
※別の方法として、「HSSEND=.TRUE.」を付けずに構造最適化計算を行った後、「Job Type」の「Vibrational Frequencies」や「IR Spectrume」で計算する方法もあります(計算内容は同じです)。構造最適化の計算を確認してから、振動数計算を実行したいときなどに使います。構造最適化がきちんと行われず連続で振動計算を実行すると、時間のムダになってしまうことがあるので、この操作は覚えておくとよいでしょう。
計算を実行しよう
「Run Firefly!」ボタンを押して、計算を実行します。最近のパソコンであれば、30分程度で計算は終了すると思います。
※ ステータス表示に「EXECUTION OF FIREFLY TERMINATED ABNORMALLY~」と表示され、アウトプットファイルに「ERROR: MEMORY REQUEST EXCEEDS AVAILABLE MEMORY」と出力されている場合は、MWORDSを増やしてみてください。キーワードのMWORDSは計算に使用するメモリ量を指定します。先の画面では、「MWORDS=250」に設定してます(Keywords欄)。単位はMW(メガワード)で、1MW=8MBです。
指定メモリ数を大きくし過ぎると、ステータス表示に「EXECUTION OF FIREFLY TERMINATED ABNORMALLY~」と表示され、アウトプットファイルには「◯◯◯◯ WORDS OF MEMORY UNAVAILABLE」とメモリ不足のエラーが出力されます。その際は、MWORDSを下げてみてください。
計算結果を確認しよう
1) ステータス表示が「GAMESS/Firefly-calculation successfully terminated. Have fun! 」になったら計算終了です。Resultsボタンの赤枠①「IR Spectrume」を押すと計算で得られたIRスペクトルが表示されます。
2) 赤枠②のFileから「Import Reference Spectrume」を選ぶと実測値との比較が行えます。Chemspider, NIST, SDBSにデータが存在する化合物であれば比較ができます。また、「Add Reference Spectrume」から自分で測定したスペクトルデータも読み込むことが可能です。その際、拡張子は.jdx、.dx、jcm(JCAMP-DXファイル形式)である必要があります(大半のIR測定装置に対応してると思います)。また、一応画像でも追加できるようになってます。今回はNISTを選んで、アセトフェノンのCAS番号「98-86-2」を入力して「Start」を押してください。
3) 赤色でNISTのアセトフェノンのIRスペクトルデータが表示されます。赤枠の「Scale Freqs」にチェックマークを入れると計算で得られたスペクトルのスケールファクター※が適応されます。
4) スペクトルを比較すると2000cm-1付近のC=O伸縮のスペクトルが高波数側に過大評価されているのがわかります。完全一致とはいきませんが、実測のスペクトルの帰属には使える結果が得られたと思います。
時間がある方はB3LYP/6-31G*で計算して結果を比較してみても面白いでしょう。
※ちなみにどのスペクトルのピークが、どの官能基の振動モードに対応しているかを確認するには、IR SpectrumウインドウのツールバーからVibrations >Interactive correlation of animated vibrations bandsを選択します。Webブラウザが起動し以下のような画面が表示されるので、確認したいピークを選択すれば、対応する振動モードのアニメーションが表示されます。
スケールファクターについて知ろう
先程、スケールファクターという言葉が出てきました。ちなみに量子化学計算では、調和振動数を計算しています。実験のスペクトルと直接比較できる基本振動数を計算で得るには非調和性というものを考慮する必要があるんですが、計算負荷がかなり大きくなります。そこで近似(調和近似)を用います。一般的にこの近似による誤差が3~5%程度と言われているため、これに対する補正としてスケールファクターを用います。
スケールファクターは、計算に使用した理論(計算手法)/基底関数で値が異なります。例えば、今回使用したHF/6-31G*のスケールファクターは「0.8985」が使用されてます。MoCalc2012では以下の文献のスケール因子を採用しています。
A.P. Scott, L.Radom, J. Phys. Chem., 100, 16502 (1996). DOI: 10.1021/jp960976r
参考にいくつかの理論/基底関数に対するスケールファクターを示します。
ちなみに、MoCalc2012フォルダ内のParamsフォルダ中にある、”FreqScalingParams.txt" に値を追加することで設定にないスケールファクターを追加することも可能です。
HF/3-21G:0.9056, HF/6-311G**:0.9085, B3LYP/6-31G*:0.9603
CCCBDBでもスケールファクターを確認できます。以下のURLにアクセスしてください。
http://cccbdb.nist.gov/vibscalejust.asp
MacでIR計算を行いたい時の対処法
MoCalc2012はWindows環境でしか実行できません。しかし、他のソフトを利用すればMacでも同様にIRスペクトルを表示させることが可能です。例えばMacの場合、AvogadroやMacMolPltなどでInputファイルを作成し、Mac版のFireflyで計算します。結果(Outputファイル)の可視化は、GaussSumを使うとこんな感じで表示できます。ただし、どの振動モードがどのスペクトルに対応しているかはアニメーションで表示することはできません。また、GaussSumをMacで使う場合は、別途Wineなどの導入が必要になります。
どうしても、振動モードのアニメーションを表示させたい場合は、WebMOのデモ版を使用する手があります。これを利用すると、以下のようなスペクト表示と振動モードのアニメーションが表示できます。
WebMOはWebブラウザがあれば良いので、新たにソフトを導入する手間はありません。詳しい使い方は各ソフトの関連記事を参照してみてください。また、この方法はWindowsでも利用できるので是非活用してみてください。
おわりに
今回、IRスペクトルの計算を例に振動数計算の方法を説明しました。実験のIRスペクトルと比較するためにはスケールファクターが必要なこと、そのためには理論/基底関数に適したスケールファクターの選択が必要であることを述べました。
先に述べたように振動数計算からは、他にも非常に有用な情報が得られます。重要な計算ですので、今後詳しく説明していきます。
この記事で使用したOutputファイルを以下に置いておきます。
Output: http://pc-chem-basics.blog.jp/AcetophenonIR(out).out
※) 当ブログの連載が書籍になりました!!興味がある方は、お手に取っていただければ幸いです!!
(Amazon購入ページ)(書籍の紹介記事)
【関連記事】
NMRスペクトルを計算してみよう
紫外可視吸収スペクトルを計算してみよう
熱力学的諸量をみてみよう
構造最適化について知ろう
コメント
コメント一覧 (23)
構造計算からスペクトルの帰属をシミュレーションしてみたくてたどり着きました。
構造計算からIRスペクトルまではできたのですが、
特定の官能基とピークの帰属を
デモンストレーションするようなことは可能でしょうか。
可能でしたら、その手法をご教示いただければと思います。
不躾な質問で恐縮ですが、よろしくお願いいたします。
どうもこのブログの管理人です。
MoCalc2012の場合は、IR Spectrumウインドウのツールバーから操作できます。
Vibrations >Interactive correlation of animated vibrations bandsを選択すればOK
です。
Webブラウザが開き、表示されるスペクトルのトップピークを選択すると対応する振動
モードがアニメーションで表示されます。
表示方法は、後日記事に追加いたします。
IRスペクトルを見たくて、指定のやり方に則って計算を行おうとしているのですがRun Firefly!を押してから一瞬でステータスにEXECUTION OF FIRELY TERMINATED ABNORMALLY AT (現在時刻)と出て動作が終了してしまいます。Outputを見るとなにやら計算はしているようですがそれ以外のメニューが灰色になっており使えません。ご指導のほどよろしくお願い致します。
どうもこのブログの管理人です。
「EXECUTION OF FIRELY TERMINATED ABNORMALLY」とはプログラムが異常終了したメッセージです。
Outputファイルの最後の方にエラーの内容が記載されているはずなので、それを確認してみてください。
エラーの内容に因って対処方法が異なります。
対処方法についてはこのブログの以下のページを参考にしてみてください。
■GAMESSの基本的なエラーの対処法
http://pc-chem-basics.blog.jp/archives/3317285.html
エラーメッセージ?を見ると
ADDRESS 0x0059D7B5 HAS INITIATED PROGRAM ABOUT BECAUSE OF FATAL ERROR(S)
と出ています。初心者なので右も左もわからず状態なのですが先程の返信にあったURLから対策を探して見ます。
構造最適化を飛ばして初期配置からIRスペクトルを予測するにはどうしたら良いでしょうか?
測定したIRの簡易的な帰属にfireflyを利用したいと考えています。
よろしくご教示ください。
どうもこのブログの管理人です。
初期配置からIRスペクトルを計算するにはMoCalc2012で、Job TypeをMolecular Energy
に変更し一点計算を行います。一点計算とは、座標を動かさないで行う計算のことです。
また直接Inputファイルから設定する場合は、RUNTYP=ENERGYにすればOKです。
ご教示いただいた方法で一度トライしてみます。
たびたびのご相談で申し訳ありません.
ご回答いただいたようにMoCalc2012でMolecular Energyで試してみましたが
計算終了後にIR Spectrumeがでてきませんでした.(当該部にはSCFが表示)
他条件は本ブログにあるように6-31G*に変更し,HSSEND=.TRUE.と
打ち込みました.まだ何か条件が足らないのでしょうか.
たびたび恐縮ですが再度ご指導いただければ幸いです.
どうもこのブログの管理人です。
申し訳ございません、失念しておりましたm(_ _;)m。
RUNTYP=ENERGYとHSSEND=.TRUE.のキーワードだけでは計算できませんね。
改めて、
MoCalc2012で、Job TypeをVibrational Frequencies
Inputファイルから設定する場合はRUNTYP=HESSIAN
(※この場合、$STATPT~$ENDにHSSEND=.TRUE.は必要ありません。)
で計算してみてください。
Inputファイルはこんな感じです。
$CONTRL COORD=UNIQUE RUNTYP=HESSIAN SCFTYP=RHF $END
$CONTRL UNITS=ANGS MAXIT=100 MULT=1 $END
$SCF DAMP=.TRUE. $END
$BASIS GBASIS=N31 NGAUSS=6 NDFUNC=1 $END
$SYSTEM TIMLIM=3000 MWORDS=250 $END
$STATPT NSTEP=200 $END
$DATA
Acetophenon
C1
O 8.0 0.0003790000 -1.1699990000 2.2991100000
C 6.0 -0.0003200000 -0.0451800000 0.2052120000
C 6.0 -0.0001220000 1.1912940000 -0.4404930000
C 6.0 -0.0002170000 -1.2225500000 -0.5426350000
C 6.0 -0.0003220000 -0.1065370000 1.6579580000
C 6.0 0.0001800000 1.2502980000 -1.8341450000
C 6.0 -0.0000150000 -1.1635460000 -1.9362880000
C 6.0 0.0001830000 0.0728280000 -2.5820930000
C 6.0 0.0000740000 1.1893160000 2.4510110000
H 1.0 -0.0001240000 2.1335040000 0.0982450000
H 1.0 -0.0003160000 -2.1966000000 -0.0606530000
H 1.0 0.0002790000 2.2127030000 -2.3376580000
H 1.0 0.0000870000 -2.0801430000 -2.5188600000
H 1.0 0.0003840000 0.1188340000 -3.6671770000
H 1.0 -0.9025260000 1.7679690000 2.2373060000
H 1.0 -0.0003260000 0.9552490000 3.5207100000
H 1.0 0.9031740000 1.7670720000 2.2376060000
$END
無事に計算できました.
ご丁寧に解説いただき大変助かりました.
ありがとうございました.いろいろと遊んでみます.
質問がありまして書き込みさせていただきます。
>別の方法として、「HSSEND=.TRUE.」を付けずに構造最適化計算を行った後、「Job Type」の「Vibrational Frequencies」や「IR Spectrume」で計算する方法もあります
とありますが、具体的にどのようにするのでしょうか。
一度構造最適化したアウトプットファイルを使って、IR計算を行うのでしょうか?
お時間ありましたら、ご回答いただけると幸いです。
よろしくお願いいたします。
どうもこのブログの管理人です
ご質問の通り一度構造最適化したOutputファイルを使って、振動数の計算を行います。
具体的には、 MoCalc2012を使う場合
1)構造最適化後にxyzファイルが出力されるので、それをMoCalc2012に読み込みます。
2) そうすると何番目の構造を使用するか?と尋ねるサブウインドが開くので一番最後の
構造最適化されたModel番号を指定します。
3) あとは、「Vibrational Frequencies」や「IR Spectrume」に設定して、計算を実行
すればその構造に関する振動数の計算が行なえます。
他のソフトを使う場合、
Outputファイルから、「EQUILIBRIUM GEOMETRY LOCATED」で検索すればその文字
列の下にxyz形式で構造最適化後の構造が出力されています。
それを使用して、新たに振動数計算用のInputファイルを作成すればよいかと考えます。
記事を参考に他の分子でも構造最適化→IR計算を行っているのですが、虚数の振動数が出てきてしまいます。
そこで解決策となる糸口を考えているのですが、ご教授いただければ幸いです。
1 アウトプットファイルから、具体的にどの原子をいじれば良さそうだ、わかりますか。
2 SCFとエネルギー閾値を、厳しく設定する(それぞれ10-7乗、10-6乗)ことで、解決しますか。
3 そのほかにアドバイスがありました、ぜひよろしくおねがいします。
どうもこのブログの管理人です。
お役に立てていれば幸甚です。
1. 振動のアニメーションで虚数振動がどの振動モードに対応しているかで
確認できると思います。
2. 収束判定を厳しくして時間をかければ求まることもありますが、一概に
そうとは言えません。
まずは、出力ファイルから収束の状況を確認して、繰り返し計算や判定値
を厳しくするだけでよいのか見極める必要があります。既におかしな構造
に収束しつつある計算を続けても意味がありません。
初期構造で類似のものがデータベースにあるのならそれを用いるか、一度
低レベルの計算で構造の見通しを立てるのも良いかと思います。
3. なにより初期構造が重要ですので、既知や類似の座標を利用するのが
早道です。それでも収束しない場合は、方法論を見直すためにキーワード
を追加・変更することを検討すると良いと思います。
以下の記事などを参考にしてください。
http://pc-chem-basics.blog.jp/archives/1736191.html
管理人さま、お返事ありがとうございます。
MacMolPitで虚数振動(値は-23.03です)を、view⇒show nomarl modeで見た際に、全74原子のうち、7つの原子から赤矢印が出ているのが確認できました。またアニメーションモードで分子全体が上下に振動しているように見えるのですが、どの振動モードかわかりません。
そこで、追加質問なのですが、
1. 振動モードはどのように判別すればよろしいでしょうか。
2. また、いじるのはこれら7つの原子ということになるでしょうか。その場合、どのようにいじればよいか、もしアドバイスがありましたら、ぜひご教授お願いします。
具体的な化合物がわからないのでこれ以上コメントは難しいのですが、
非直線分子の場合3N-6個の振動モードがあるのことをご存知ですか?
まず基準振動について、少し勉強してみるとよいかと思慮いたします。
虚数振動の変位はエネルギーの低下をもたらすので、より安定な構造へ
の変化を意味します。虚数振動モードの方向に変位の大きさにしたがっ
て原子座標をずらした構造から,最適化を再度行う方法があります。
その場合、振動解析の力の定数を読み込ませてリスタートするとよいで
しょう。
リスタートの方法は以下を参照ください。
http://pc-chem-basics.blog.jp/archives/3807643.html
ただし、対象系によっては初期構造を見直して再計算した方が効率的な場合
のことが多くありケースバイケースです。
難しい事がわからない場合は、先の回答の通りもう一度初期構造を見直して
既知の類似構造を使うか、計算レベルを落とした安定構造で再度計算する方
が手取り早いと考えます。
管理人様
いつもお世話になっております。返信が遅くなりましてすいません。
お忙しいところご丁寧に回答頂きまして、ありがとうございます。
無事、構造最適化後、振動数計算ができました。
ありがとうございました。
お世話になります。
ブログを参考にさせていただきまして、なんとかIRスペクトルを出力することまではできるようになりました。
リファレンスのIRスペクトルをImportする際に、"リモートサーバーがエラーを返しました:(407)プロキシ認証が必要です"という警告がでて、上手くいきません。プロキシの突破法をご存じでしたら、ご教示いただけますと幸いです。
よろしくお願い申し上げます。
A
会社員Aさんこんにちは!!
例えば、社内ネットワークなどを使っておられて、プロキシを使用していない自宅など他のネットワークでは407エラーが出ないにであれば、管理者に問い合わせてみるのがよいかもしれません。
407エラーは私も遭遇したことがないので、お力になれず申し訳ありません。
別のソフトウェアでも同様なことが頻繁にあるので、社内ネットワークが原因かなあと考えています。プロキシの件は社内のネットワーク担当に聞いてみることとします。
また、別件で下記の分子の最適化計算を行おうとしているのですが、スピン多重度を1重項に指定して計算を行うと
エラーが出て計算が止まってしまいます。 .gamに記載のWarningにはMultiplicity was changed 1 to 2 for job consistency.とあるのですが、、
解決策がございましたら、ご教授いただければ幸いです。
SMILES:FC(F)(C(F)(F)C(F)(F)F)C(F)(F)C(F)(F)C(F)(F)CC
会社員A
おそらく、初期構造、又は電荷、スピン多重度、閉殻系などの入力ファイルで指定する電子状態に矛盾があるためかもしれません。入力ファイルを一度見直して初期構造や正しい電子状態を指定すると解消するかもしれません。
ご参考まで。
Inputファイルに矛盾があったことが原因でした。ChemSketchから.molを出力した際に水素の情報が抜けていたためでした。