03_geomopt - YoshioNishimoto/sandbox GitHub Wiki

構造最適化

構造最適化（geometry optimization）をします。 構造最適化は、ポテンシャルエネルギー面（PES; potential energy surface）の極値を見つけることです。 PESは、教科書とかでは三次元で書かれますが、実際には3Natom-6次元（自由度に一致。直線分子の場合は？）の超曲面になります。 イメージしづらいですが。 「極値」だと、極小値・極大値のどちらを計算しているのか分かりませんが、とりあえずここでは極小値を計算することを考えています。 極大値の計算はまた別の機会に（できれば･･･PySCFではできない？）。 極大値は一般的に鞍点（あんてん; saddle point）と呼ばれており、3Natom-6つの基準振動の中で一つの振動モードのみが上に凸の場合に、遷移構造・遷移状態（TS; transition structure, transition state）と呼ばれています。

PESの最小値を求めるというのは、前回の演習問題で行いました。 ただ、この場合は二原子分子で自由度が1であったため、自由度に対応する原子間距離を変化させながらエネルギーを計算し、横軸に原子間距離と縦軸にエネルギーをプロットすると、極小値を同定することができました。 ですがこれ以上大きな分子、三原子分子であったとしても自由度は3つとなり、例えば自由度の軸に沿ってそれぞれ10個の一点計算をするとしても、三つの軸に沿って原子を動かしていくと、最終的に1000回の一点計算が必要になります。 しかも、この計算を行ったとしても必ず極小値が見つかるわけではないですし、しかもそれぞれの軸に沿って10個の計算を行ったとしてもかなり荒い近似的な極小値（せいぜい精度は0.1 Angstrom程度と考えられる）しか見つかりません。 そこで、構造最適化をして自動かつ正確に極小値を求めます。

構造最適化するためには、勾配（gradient）を計算します。エネルギーを原子の核座標に関して一次微分したものです。 計算した勾配は、構造最適化中に原子を動かす方向を決めるのに必要で、そして収束条件の判定（極値では勾配がゼロになっているはずなので、ある閾値よりも勾配が小さくなっていれば、構造最適化が収束したと判断する）にも用います。 計算の詳細は置いておいて･･･この勾配は多くの場合、解析的（analytic）に計算します。 これは高校で勉強したような、f(x) = xnからf'(x) = nxn-1を計算するような方法です。 高校で勉強したと思われる導関数の定義を用いて、数値的（numerical）な計算をしているわけではありません。 つまり、勾配を計算する際に、少し原子を動かしてエネルギーを計算し、そのエネルギーの差分をとる･･･ということはしていません。 原子は動かしますが、別の理由です。

量子化学計算では、多くの場合まず構造最適化計算をして、その構造でのエネルギーやプロパティを用いて議論を行います。 というのも、エネルギーが低い構造ほどその分子構造（配座？）をとる確率が高くなるためです。 一点注意すべき事は、ここで求められる極小値は「局所的な」極小値（local minima）であり、必ずしも「大域的な」極小値（global minimum）であるとは限りません。 後者は、3Natom-6次元のPESの全域を調べないと発見できないため、新しいglobal minimumを見つけたと主張する際には注意が必要です。

構造最適化・エネルギー微分の概要としては、Szabo: Appendix Cが簡単かと思います。

実際の計算

詳しいことはマニュアルをご参照ください。 とりあえず、こんな感じ。

  1 from pyscf import gto, scf
  2 from pyscf.geomopt.berny_solver import optimize
  3 mol = gto.M(atom='N 0 0 0; N 0 0 1.2', basis='ccpvdz')
  4 mf = scf.RHF(mol)
  5 mol_eq = optimize(mf)
  6 print(mol_eq.atom_coords())

2,5,6行目を新しく追加しました。 2行目は構造最適化のドライバーを利用するための行、5行目は実際に構造最適化を行う行です。 私の環境だと、PyBernyのバージョンが低いとか言われました。 その場合は指示に従ってアップデートしましょう。 geomeTRICという外部のドライバーを使っての構造最適化も可能とのことなので、興味があればやってみてください。 上の入力ファイルは、N2分子を(R)HF/cc-pVDZという計算レベルで構造最適化を行っています。

最適化が終了した構造

6行目は、最適化した構造の座標（coordinate）を表示しています。 計算の最初の構造（N-N間の距離が1.2 Angstrom）はmolに入っていますが、5行目で構造最適化をした際に、最適化した構造はmol_eqに入ります。 上の例だと、このような出力になると思います。

[[0.         0.         0.11593315]
 [0.         0.         2.1517382 ]]

なお、単位は原子単位（a.u. or Bohr）です。 ゼミや学会での発表では、多くの場合オングストローム単位で結果を示すため、単位を変換して結果を示すようにします。 「N2をHF/cc-pVDZで構造最適化すると、結合距離が1.077 Åになった」とでも表現できるでしょうか。 次の計算で使いたい場合には、単位が異なることに気をつけてください。 unit="Bohr"を用いると、原子単位で構造を指定することができます。

mol = gto.M(atom='N 0 0 0.11593315; N 0 0 2.1517382', basis='ccpvdz', unit="Bohr")

構造最適化の詳細

エネルギーの一次微分は割と簡単に計算できるのですが、高校数学のように、一次微分がゼロになる点を直接計算できるわけではありません。 このため構造最適化は、Hartree--Fock法（あるいはSCF法）の計算と同じように、繰り返し計算により極小値を計算します（cf. ニュートン法）。 一次微分がゼロになる方向に原子を動かして、その後でエネルギーと一次微分を計算して、さらに原子を動かして･･･ということをしています。 構造最適化のそれぞれのサイクルの中で、エネルギーを計算するときにHatree--Fock(--Roothaan)方程式を繰り返し計算で行っています。 繰り返し計算というおおざっぱな言葉が出てきていて分かりづらいですが、やっていることは割と違います。

まずは構造最適化サイクルの1サイクル目を見てみます。

  2 Geometry optimization cycle 1
  3 Cartesian coordinates (Angstrom)
  4  Atom        New coordinates             dX        dY        dZ
  5    N   0.000000   0.000000   0.000000    0.000000  0.000000  0.000000
  6    N   0.000000   0.000000   1.200000    0.000000  0.000000  0.000000
  7 converged SCF energy = -108.914051975052
  8 --------------- SCF_Scanner gradients ---------------
  9          x                y                z
 10 0 N    -0.0000000000    -0.0000000000    -0.3121321959
 11 1 N     0.0000000000     0.0000000000     0.3121321959
 12 ----------------------------------------------
 13 cycle 1: E = -108.914051975  dE = -108.914  norm(grad) = 0.441422

New coordinatesのすぐ下の行が、最初に入力した構造です。 7行目にconverged SCF energy = -108.914051975052と表示されており、 ここではverbose=0のため何も表示されていませんが、前のページでご説明したとおりのSCF計算が行われています。 そして10・11行目に、その分子構造での（解析的）勾配が表示されています。 a.u./bohrという単位を使っています。 なお、勾配はSCF計算が収束した後（エネルギーが得られた後）でないと得られないことに注意してください。 なので、SCF計算が収束しない場合は、その時点でプログラムが異常終了することが多いです。

で、上のサイクルが終わった後に、次の構造最適化サイクルに移ります。

 15 Geometry optimization cycle 2
 16 Cartesian coordinates (Angstrom)
 17  Atom        New coordinates             dX        dY        dZ
 18    N   0.000000   0.000000   0.079377    0.000000  0.000000  0.079377
 19    N   0.000000   0.000000   1.120623    0.000000  0.000000 -0.079377
 20 
 21 WARN: Large deviations found between the input molecule and the molecule from chkfile
 22 Initial guess density matrix may have large error.
 23 
 24 converged SCF energy = -108.95047828901
 25 --------------- SCF_Scanner gradients ---------------
 26          x                y                z
 27 0 N     0.0000000000    -0.0000000000     0.1553466138
 28 1 N    -0.0000000000     0.0000000000    -0.1553466138
 29 ----------------------------------------------
 30 cycle 2: E = -108.950478289  dE = -0.0364263  norm(grad) = 0.219693

18・19行目で表示される構造が、1サイクル目の構造とは違います。 21・22行目は気にしないで良いです。 24行目では、再びSCF計算をしています。 構造最適化サイクルで、毎回SCF計算と解析的勾配の計算をしているわけです。 エネルギー・勾配ともに減少していることに注意してください。 極小値ではエネルギーが局所的に最小値となっているはずなので、エネルギーが減少するという挙動は正しいです。 また、極値では勾配がゼロに近くなっているはずなので、勾配が減少するという挙動も正しいです。 ただ、必ずしも単調減少していくわけではありません。 少なくとも3Natom-6次元でのパラメータの最適化なので、エネルギーが増えたりすることもよくあります。

そして最後の構造最適化サイクルでは、

 85 --------------- SCF_Scanner gradients ---------------
 86          x                y                z
 87 0 N    -0.0000000000    -0.0000000000    -0.0000031240
 88 1 N     0.0000000000     0.0000000000     0.0000031240
 89 ----------------------------------------------
 90 cycle 6: E = -108.955558787  dE = -4.74512e-08  norm(grad) = 4.41806e-06

というように、勾配の値がかなり小さくなります（3.12D-6）。 verbose=4とすると、追加でConvergence criteria（収束しきい値）に関する情報も表示され、

6 Convergence criteria:
6 * Gradient RMS: 3.12e-06 < 0.00015 => OK
6 * Gradient maximum: 3.12e-06 < 0.00045 => OK
6 * Step RMS: 1.55e-06 < 0.0012 => OK
6 * Step maximum: 1.55e-06 < 0.0018 => OK
6 * All criteria matched

というように、勾配のRMS (root mean square)、勾配の最大値、ステップ（原子の変位）のRMS・最大値を計算し、 それらの値全てがしきい値（0.00015, 0.00045, 0.0012, 0.0018）より小さくなった場合に、「構造最適化計算が収束」したと判定します。 構造最適化した構造を議論する際は、最後に出てきた構造を用いて議論すれば良いと思います。

基準振動解析

別ページに書きました。 極小値であるか極大値であるかは、二次微分を計算すると分かるのですが、 なんか直線分子の場合はバグっているような･･･。 PySCFディレクトリにあると思われるexamples/10-thermochemistry.pyをご参照のうえ、freq_infoの結果を表示してみましょう。 極小値（付近）の場合はfreq_auやfreq_wavenumberがすべて実数になります。

演習問題

  • NH3を、6-31Gと6-31G(d)基底関数で構造最適化してみましょう。 どちらの構造が正しいでしょうか。 そして、なぜ正しくない構造が得られるでしょうか。 基底関数に含まれている関数の角運動量をもとに、考えてみましょう。 保留
  • 上記の構造最適化1サイクル目に表示されている
13 cycle 1: E = -108.914051975  dE = -108.914  norm(grad) = 0.441422

のうち、norm(grad) = 0.441422はどのような計算式により計算されているでしょうか。

  • ベンゼン（でも何でも良いですが）の「MP2/cc-pVDZ//HF/6-31G(d)」計算をしてみましょう。 ここでは、HF/6-31G(d)の構造最適化を行い、最適化した構造を用いてMP2/cc-pVDZによる一点計算（エネルギー計算）をするという意味です。 ご参考のために、H2のMP2/cc-pVDZ//HF/6-31G(d)は、E(MP2) = -1.15478598152265 E_corr = -0.0262803424848388となりました。
  • cis-1,3-butadieneとtrans-1,3-butadieneを任意の計算レベル（手法や基底関数を自由に決めて良い）で構造最適化を行い、どちらが何kcal/mol安定か調べてみましょう。結果を示すときは、計算レベルも記載するようにしましょう。
  • ディールス・アルダー反応の反応エネルギーを計算してみましょう。とりあえず、反応物（別々に）と生成物の構造最適化を行い、左辺と右辺とでエネルギーを比較します。左辺と右辺、どちらのエネルギーがどれだけ低くなるでしょうか。またこの結果から、例えば常温で常識的な時間で反応が起こる・起こらないについて議論することができるでしょうか。時間があれば、基底関数を変化させたときにどの程度計算時間やエネルギー差が変化するか見てみましょう。分子の初期構造は可視化ソフトウェアで作成することをオススメします。
  • 時間があるなら、数値的に勾配を計算するプログラムを作成して、解析的な勾配と一致することを確認してください。
⚠️ **GitHub.com Fallback** ⚠️