構造最適化の閾値は、何を意味しているのか?

Gaussian での構造最適化時にlogファイルに表示される 4 つの値。これらが全て閾値 (threshold) 以下になると構造最適化が完了します。

         Item               Value     Threshold  Converged?
 Maximum Force            0.000000     0.000015     YES
 RMS     Force            0.000000     0.000010     YES
 Maximum Displacement     0.000030     0.000060     YES
 RMS     Displacement     0.000006     0.000040     YES
 Predicted change in Energy=-1.233934D-12
 Optimization completed.
    -- Stationary point found.

では 4 つの値はそれぞれ何を意味しているのでしょうか?

また、構造最適化を効率的に行うにはどうしたらよいのでしょうか?

英語版What Does the threshold in geometry optimization mean?

以前の記事で説明したように、gaussian では Berny アルゴリズムで構造最適化を行います。すなわち、ポテンシャルエネルギー局面 (PES) の一次微分により傾きを求め、限りなく 0、すなわち平らになるところを探索していきます。一方、Berny アルゴリズムの欠点としては収束点に近づくにつれ、一回ごとの step 幅が小さくなっていってしまい遅くなるということが挙げられます。簡単に言い換えれば、谷底に近づくほどスピードが落ちるということになります。

以上のことを踏まえて、構造最適化を数学的に考えると、以下のような概念図に行き着きます。
optimizetion_threshold2

PES 上の一次微分値は Force となります。この一次微分値は、PES の傾きを表します。
Force: F = –∇E
PES は、3N-6 次元なので、F は 3N-6 個あります。上図の Maximum force を表す式 max(F) < threshold は、全ての振動方向でPES の傾きが 0 に近いということを表しています。つまり、local minimum または local maximum にいるということを表します。

RMS というのは root mean square の略語で、二乗平均平方根を意味しています。RMS の大きさが閾値以下というのは、3N-6 個の F の平均値を出しています。構造最適化を行う際 Maximum Forceと RMS Forceでは、RMS の方が先に閾値以下になり YES と表示されるはずです。

Displacement は PES 上での探索においてのどれくらい座標が移動するかの指標となります。構造最適化計算時に PES 上のある点 r から次の点 r+1 に移動する時、その座標は eigenvalues と eigenvectors と F の関数で表されます。この計算により得られる r+1 と r の差が閾値以下になると PES で local minimum または local maximum にいると言うことができます。これは、谷底に近づくとステップ幅が小さくなっていくという Berny 法の特徴とも一致します。
(r+1 の計算方法などについては、後日また記事にします。)

この Displacement が収束しにくい分子というのは、PES 上の山または谷が急であるケースが考えられます。つまり、ある点 r と移動後の点 r+1 の間に収束点がきてしまうために、収束しづらいということになります。このような場合には、構造最適化のステップ幅を小さくすると改善します。
また、異様に構造がフレキシブルな場合なども収束が難しいです。そのような場合は、基底関数のレベルを下げると PES の傾きがなだらかになり問題が解決する場合があります。

構造最適化のコツ

以上述べてきたように、構造最適化では PES 上の傾きをもとに判断して構造を少しずつ変化させていきます。そのため、計算の初期座標から最も近い安定な構造(local minimum)を目指して動いていきます。以下に図を示します。

上図では、最も安定な点(global minimum)は、F ですが、初期構造として A や C を指定してしまうと B が計算結果として得られてしまいます。計算機は、どこが global minimum かは判断せず、傾きのみで判断してエネルギーが小さくなる方向に動いていくからです。

F にたどり着きたい場合は、D や E を初期構造として指定する必要があります。この場合、E を初期構造とした方が D を初期構造とする場合よりも計算が早く終わります。
量子化学計算は、計算コストが高いため、なるべく結果に近い構造を初期構造とした方が良いです。

では、どうすれば良い初期構造を得られるのでしょうか?

一つの方法としては、計算コストの低い計算手法を使うことが挙げられます。

手書きで E や F のような構造を書くのは非常に時間がかかるため、まずは D の構造を書き、半経験的手法(PM6など)で E に近づけます。さらに、B3LYP/6-31G(d) などで最適化し、その結果をもとに B3LYP/6-31+G(d,p) などで構造最適化すると、効率的に計算を進めることができます。

量子化学計算では、初期構造が非常に重要であり、計算時間・結果の良し悪しを左右します!

どのようにして、結果を判断するか?

得られた構造が本当に global minimum かどうかを判断するには、有機化学の知識が必要です。こう書いてしまうと、計算化学って思ったより主観が入ると思われるかもしれませんが、その通りです!

もちろん教科書からすると例外的な構造が最安定構造になる場合もありますが、そのような例はごくわずかです。計算結果を目視で確認し、置換基同士の相互作用や環のコンフォメーションなどおかしいと思ったら、構造を少しずつ変化させて global minimum を探すのが良いです。

もっと客観的な方法はないのか?

主観を入れずに客観的に判断するとしたら、全自動探索(GRRM など)、または MM 計算と組み合わせるのが良いかと思います。

しかし、QM での全自動探索は計算コストが非常に高くて現実的ではない、MM 計算と組み合わせると数百から数千通りのコンフォマーが出てくるため、それらを再度 QM 計算して絞り込んでいくなど、非常に面倒です。

参考QM でのコンフォメーション探索

参考MD 計算を利用したコンフォメーション探索【gromacs】

関連する記事



コメントを残す(投稿者名のみ必須)