ORCAのIGLO法でNMR化学シフトの高速計算
GAMESSを使ってNMRの化学シフトを計算することができますが(記事)、この計算にはかなりの時間とメモリ量を要し、量子化学計算で算出する利点の数々もこのリソース消費との比較で霞んでしまうことがあります。
商用プログラムでは化学シフトの高速計算を実現しているものもありますが、無償で利用可能なプログラムの中にも高速で計算できるものがあります。
ORCA(http://www.thch.uni-bonn.de/tc/orca/)は、ドイツ・ボン大学のFrank Neeseらによって開発・メンテナンスが行われているプログラムですが、このプログラムにはIGLO法によるNMR遮蔽定数計算が実装されています(しかもLinux版では並列化されています)。ORCAのIGLO法によるNMR計算では、GAMESSのGIAO法による計算に比べて圧倒的に計算時間が短く、しかもメモリ消費量も格段に少なくなっています。ppm換算で小数点以下1桁までしか出力されないため、1H-NMRの計算には不向きですが、13C-NMRの計算には十分活用可能です。
果たしてその精度はどれくらいなのか、実際にいくつかの有機分子で計算を行ってみました。
代表的な有機低分子での計算
一般的な13C-化学シフトの範囲は0～220 ppmほどですので、その範囲をだいたい網羅するように9種類の有機分子を選び、CDCl3中での実測13C-化学シフト値※1をどれくらい再現できるかを調べました。9種の有機分子と、基準となるTMS(Tetramethylsilane)の構造を以下に示します。
構造最適化と振動解析はB3LYP/Def2-SV(P) (溶媒和はCOSMO法でCHCl3を考慮)で行い、得られた構造を使い、HF/IGLO-II, BLYP/IGLO-II, B3LYP/IGLO-IIの3種類のモデル化学でIGLO法による遮蔽定数の計算を行いました。代表的な入力ファイルの一例を以下に示します。
(構造最適化と振動解析)
遮蔽定数計算の例を見るとわかるように、Gaussianのように複数のジョブをつなげて一つの入力にすることができます。構造最適化→振動解析→遮蔽定数計算を一つの入力ファイルで実施することも可能です(上記では2つに分けていますが)。また、構造最適化計算で出力されるxyzファイルを参照することで、座標の入力を省略することもできます。
このような連続ジョブの利点は他にもあって、最初のジョブで得られたSCF軌道を次のジョブのinitial guessとして使うことが出来、これによってSCF計算のステップ数を削減できる可能性もあります。
さて、実測値と計算値の比較が以下のグラフになります。
このグラフは、生の計算値と実測値の相関をExcelで計算し(切片0の線形近似)、その傾きの逆数でスケーリングしたものを計算値(scaled)として使っています。遮蔽定数の計算結果は、化学シフト値の増大に伴って実測とのずれが広がる傾向にあり、振動計算と同様にスケーリングすることで実測値との一致を良くすることができます(scaling factorの研究例もあるようです)。
こうして眺めてみると、HF(WFT)とBLYP(DFT)でその傾向が真逆であることがわかります。HFが30～100 ppmで過小評価、100～170 ppmで過大評価するのに対し、BLYPは40～100 ppmで過大評価、140以上で過小評価する傾向にあります。両者をハイブリッドしたB3LYPではその中間的な値となり、全体的によい一致を示しています。また、DFTでは塩素が結合した炭素の化学シフトを過大評価し、特に2,2-DichloropropaneのC2は大きく逸脱しています。実測値との誤差という意味ではDFT(特にhybrid DFT)が優れているようですが、大小関係の正確さという意味ではHFも優れています。
GIAO/GAMESS vs IGLO/ORCA
さて、計算時間がGIAO/GAMESSとIGLO/ORCAでどれくらい異なるかを比べたのがTable 1です。TMSの遮蔽定数計算における計算時間をGIAO/6-31G(d)とIGLO/6-31G(d)で比較しています(SOSCF, SCF convergence threshold=1.0d-7に共通設定)。
言うまでも無く、圧倒的にORCAのIGLO法が速いです(100倍以上)。使用メモリ量についてはORCAの出力ファイルに記載されないため正確な数値はわかりませんが、明らかにORCAの方が少なくて済みます。
| Wall Clock Time (sec) | ||
| 1 proc. | 4 proc. | |
| GIAO (GAMESS(US) 2009R1) | 4383.9 | --- |
| IGLO (ORCA 2.6.35) | 36.7 | 13.9 |
大きめな分子での計算
今までに計算してきた分子はいずれも小さなものでしたが、実際に計算したい分子というのはもっと大きいものです。 IGLO/ORCAが非常に効率よく13C-化学シフトの計算を行えることを先ほど示しましたので、より「現実的」なサイズの分子で計算を行ってみました。
Triamcinolone Acetonide(以下TAと略)は、分子量434.5の合成コルチコステロイドで、中時間作用型ステロイド系抗炎症薬として利用されています(参考：英語版wikipedia)。この分子の13C-NMRスペクトルは16.26～209.75 ppmと広範囲にピークが散在し、また配座が固定されているため、ベンチマークとして非常に都合の良い分子です。TAの13C-化学シフトの計算を、小さな分子の時と同じモデル化学で実施しました。(※ただし構造最適化はRI-mPWLYP/SV(P)で実施し、振動解析は時間がかかるため省略)
傾向は先ほどの例とあまり変わりません。HFとBLYPが逆の傾向を示し、その中間を行くB3LYPがもっとも精度が高い結果を与えています。また、先ほど塩素が結合した炭素でDFTの誤差が大きかったわけですが、こちらでフッ素が結合した炭素(C8)の精度は悪くありません。
計算時間は構造最適化(RI-mPWLYP/SV(P))が63.6 min, 化学シフト計算がHF/IGLO-IIで113.5 min, BLYP/IGLO-IIで11.6 min(!!!), B3LYPで80.8 minでした(Linux版で4プロセス並列計算の場合)。ORCAはDFT計算が非常に高速で、GGA functionalを利用した場合、計算内容によってWFTよりも4～40倍高速とのことです。TAの最適化をRI-mPWLYPで行ったのは、そこに理由があります(B3LYPでは時間がかかる)。もちろんpure DFTの問題点はいろいろあるでしょうが、それらを理解した上で利用すれば、かなり計算時間を短縮することができます。
実測データ出典
SDBSWeb : http://riodb01.ibase.aist.go.jp/sdbs/
(National Institute of Advanced Industrial Science and Technology, 2009/06/07)
[本稿で行った計算の入出力ファイル(orca_iglo.zip), 1.35 MB]
[本稿で行った計算の結果をまとめたエクセルシート(orca_iglo.xls), 76.0 KB]