PROTACの対となるか？〜脱ユビキチン化を促進するキメラ化合物〜

AmgenのKRAS-G12C阻害剤 Sotorasib(AMG510)が承認されましたね！*1 ケミカルバイオロジーの発展に端を発する化合物で、難攻不落と言われていたKRASを標的とする医薬品が誕生したのはとても素晴らしいです。*2

低分子医薬品の世界も、標的となるタンパクが広がったり、PROTACにMolecular glueといった新しいアプローチが実用化に向けてどんどん開発されていて眺めていて楽しいです。

で、コンセプトが明確で、なるほど面白い！そうなるかー！と思ったアプローチを知りました。「Targeted protein stabilization (TPS)」というもので、雑にまとめると「PROTACの逆をいくキメラ化合物」です。

bioRxivに論文が公開されていたのでご紹介。*3

www.biorxiv.org

Targeted protein stabilization

さて、こちらの論文で提案されているアプローチはとても分かりやすいです。

「ユビキチン-プロテアソーム系で分解されるタンパク質に対して、脱ユビキチン化酵素を働かせることで、 分解を抑制して安定化し、機能できる状態を維持しよう」ってな感じです。

変異などにより上手くフォールディングできないタンパク質はユビキチン化を受けプロテアソームで分解されます。 PROTACやMolecular glueはこのシステムを利用して、壊したいタンパク質のユビキチン化を誘導、分解促進するメカニズムでした。

今更感あふれる図

で、上の図には出てきませんが、 ユビキチン-プロテアソーム系を調節する酵素として脱ユビキチン化酵素(deubiquitinating enzyme、DUB)があります。 その名の通り「タンパク質からユビキチンを切断除去するプロテアーゼ」です。*4

Wikipediaの模式図を貼っておきます。上図の逆向きの反応です。

さて、著者らのアプローチはこうです。 「PROTACがユビキチンリガーゼを標的タンパクに近接させることで分解を促進するのなら、 脱ユビキチン化酵素を近接させれば分解を抑制・安定化できるだろう。そのような近接化に働くキメラ化合物を Deubiquitinase-Targeting Chimera(DUBTAC)と呼ぼう。」

bioRxiv 2021.04.30.441959 Fig 1aより

比較するとこんな感じの話？

やっつけの図

なるほど！！でもどうやってそんな化合物作れば良いんですかね？？？

DUBTACのデザイン戦略

DUBTACによるTPSを達成するためには大きく4つの課題があります。

1. DUBファミリーのうち、どの脱ユビキチン化酵素を利用するか？
2. また、その酵素に結合するユニットをどう見つけるか？
3. 安定化する標的タンパク質として何を選ぶか？
4. また、その標的タンパク質と結合するユニットをどう見つけるか？

著者らは、前２つの課題に対してはケモプロテオミクスを利用することで、 後２つ課題については嚢胞性繊維症(Csytic Fibrosis, CF)の原因CFTRを標的とすることでコンセプト検証（Proof Of Concept）を行っています。

Nomura Research Groupについて

さて、ここで著者情報です。今回のbioRxiv論文は University of California, BerkeleyのDaniel K. Nomura教授らのグループによる報告です(Nomura Research group HP)。 所属としてNovartis-Berkeley Center for Proteomics and Chemistry Technologiesとも記載されています。

こちらはNovartisとUC Berkeleyとが協力し、 「"undruggable"な標的に対する低分子化合物を開発するための次世代の技術」を開発するための組織とのことです*5。新しいモダリティとして深く製薬企業が関わった産学連携の成果なのかもしれません。*6

Nomura研究室の得意とするところは「ケモプロテオミクスを利用した共有結合性リガンドの発見」のようです。

より具体的には「活性ベースタンパク質プロファイリング（Activity-Based Protein Profiling, ABPP）」を利用することで、 「"undruggable"な標的タンパク質に、リガンドが結合可能なユニークなホットスポットがあるか？」、 「ポケットの探索およびそのリガンド（functional covalent ligand）の取得」ができるそうです。*7

今年３月、Acc.Chem.Resに以下のようなレビューも出版されています。

Reimagining Druggability Using Chemoproteomic Platforms Jessica N. Spradlin, Erika Zhang, and Daniel K. Nomura Acc. Chem. Res. 2021, 54, 7, 1801–1813 https://pubs.acs.org/doi/10.1021/acs.accounts.1c00065

では、具体的にどのようにDUBファミリーの探索を行ったのか？、見てみましょう。

DUBファミリーからのOTUB1の選別

さて、DUBファミリーのうち今回の目的に適した脱ユビキチン化酵素の候補をどのように選別すれば良いでしょうか？

著者らは66種の酵素について、ABPPによるケモプロテオミクスのデータを集めて解析を行いました。 このデータは「酵素のシステイン残基(Cys)について、プローブと反応する残基の位置をラベリングする」というプロファイリングを行ったものです。

目的に合う酵素は以下のような条件を満たすCys残基をもつものです。

1. アロステリックな位置 … 触媒ドメインのCysでは、TPSで利用したい酵素本来の活性が失われてしまう
2. 同一酵素内での選択性 … 同じ酵素の複数のCysのうち、特定のCysがよりラベル化されていればホットスポットの可能性が高い
3. 十分な反応性を示すCys残基 … 複数のケモプロテオミクスデータセットで繰り返しプローブと反応し、ラベリングされているもの

要するに、「触媒活性を邪魔せず、かつ、リガンドがアクセスしやすく反応性が高いCys残基」があればよい、という感じです。

以上を目標に解析を行った結果、OUTB1という脱ユビキチン化酵素が、 触媒C91とは別の箇所C23に適したCys残基を持つことがわかりました。

bioRxiv 2021.04.30.441959 Fig 1b,c,dより

ところで、シグナル分子としてのユビキチン修飾は多様性があり、ポリユビキチン鎖の連結の仕方によってシグナルの意味合いが異なることが知られています。 プロテアソーム分解のシグナルとなるのは、「ユビキチンの48番目のリジン」を介して連なったポリユビキチン鎖(K48鎖)です。

ユビキチン修飾の多様性とシグナル

幸いなことに、OUTB1は発現量が高く、かつK48鎖を標的として切断する脱ユビキチン化酵素でした。 まさしく今回の目的にピッタリな酵素だったわけです。

共有結合性リガンドの探索

標的のDUB、OUTB1が定まったので、この酵素に対して適した共有結合性のリガンドを探索しています。ここでもABPPが活躍しています。

探索結果、共有結合部位としてアクリルアミドユニットをもつ化合物EN523を見出しています。 EN523はC23選択的に反応する化合物で、OUTB1の触媒ドメインのCys（C91)との反応はみられませんでした。

EN523の結合によりOUTB1の脱ユビキチン化活性が失われないことは別途in vitro再構成系で確認しています。

bioRxiv 2021.04.30.441959 Fig. 2a,b,cより

これで、DUBTACデザイン戦略のうち、脱ユビキチン化酵素側の課題2つがクリアできました。 つづいてPOC取得に向けた標的タンパク質とリガンドの選択です。

変異CFTRと低分子医薬品 Lumacaftor

DUBTACを開発するモチベーションとして、疾患の中にはミスフォールドする変異タンパクしか作れなくなり、 ユビキチン-プロテアソーム系で分解されて、機能を発揮できなくなってしまうことでおこるものがある、ということがあげられます。

そのような疾患の例として、CFTRという遺伝子の変異でおこる嚢胞性繊維症（Cystic Fibrosis、CF)があります。 CFでは複数の変異が知られていますが、中でも最も頻度の高いPhe508の欠失したΔF508-CFTRでは、 タンパク質が不安定化し、K48ポリユビキチン化修飾-分解されてしまいます。

ΔF508-CFTRに対してはすでに低分子医薬品Lumacaftorが開発されています。 LumacaftorはCFTRに結合して安定化し、フォールディングを助けるケミカルシャペロンと呼ばれる種類の化合物です。

しかしながらLumacaftorを用いても完全にCFTRタンパクをレスキューできるわけではなく、多くが分解されてしまうそうです。 つまり、「タンパク質本来の機能を妨げないが結合するリガンド」があり、尚且つまだ効果が不十分というわけです。

このリガンドをDUBTACの標的タンパク側の結合ユニットとして使ってキメラ分子を作成すれば、 さらなるタンパク安定化効果が期待でき、コンセプト検証できるのでは？？、ってなわけです。

DUBTACの合成とTPS効果の検証

さて、DUBTACのコンセプト検証に向けた材料が出揃いました。

脱ユビキチン化酵素 OUTB1と結合するリガンド EN523と、標的タンパクCFTRに結合するリガンド Lumacaftorを長さの異なるアルキル鎖（C3、C5）でつなぎ、 2つの化合物NJH-2-056、NJH-2-057が合成されました。

bioRxiv 2021.04.30.441959Fig. 3a,bより

これらの化合物により変異CFTRの安定化は達成されたでしょうか？ ΔF508-CFTRを発現するヒト気管支上皮細胞株（CFBE41o-4.7）をもちいてその効果を確かめています。

bioRxiv 2021.04.30.441959 Fig. 3e,fより

上図の通り、リンカーのアルキル鎖C5としたNJH-2-057において、Lumacaftor単独では観測できなかったCFTRの安定化、増加が確認されました。

興味深いことに、アルキル鎖を短くしたNJH-2-056では効果が確認できず、リンカーの長さの影響を受けることがわかりました。

論文ではこの他に「NJH-2-057が狙い通りのDUBTACとして機能しているか?」メカニズムに関する実験や、「安定化効果が標的タンパク（ΔF508-CFTR)選択的なものか」 といったことに関する実験も行われています。

以上、bioRxivの論文の内容でした。

感想

さて、DUBTACによるTargeted protein stabilization (TPS)、コンセプトが明快で効果検証されており、非常に面白いアプローチだと思いました。

共有結合性のリガンドをうまく活用して、新しい機能を持つハイブリッドな化合物をつくってしまうあたり、 ケミカルバイオロジーの面目躍如という感じでとてもワクワクしますよね！

ただし今回の論文は細胞レベルでの効果検証までなので、実用化にはまだまだハードルはたくさんありそうです。

• 分子量800近くあるけど、膜透過や薬物動態はどうなの？
• DUB側が共有結合性だけど、効果の持続とか、DUBの半減期の影響は？
• ユビキチンを付け外しが繰り返されると小胞体ストレスかかりそうだけど大丈夫？　　
• あれ、そもそもLumacaftorのシャペロン効果が今回の系で確認できてないのは良いんですか？

などなど、私には一読では分からないことがいっぱいです。詳しい方教えてください。

期待

で、いちゃもんつけつつワクワクしているのは、アカデミックな検証だけでおわらないことに期待しているからです。

まず、Novartis社が関わっていそうという点。嚢胞性繊維症の治療薬というと私はVertex社のイメージが強いです。*8 ところが不勉強なので、今更ながらNovartis社がCFTR標的薬 icenticaftor(QWD251)を開発されていることを知りました。

Discovery of Icenticaftor (QBW251), a Cystic Fibrosis Transmembrane Conductance Regulator Potentiator with Clinical Efficacy in Cystic Fibrosis and Chronic Obstructive Pulmonary Disease Journal of Medicinal Chemistry Article ASAP DOI: 10.1021/acs.jmedchem.1c00343 https://pubs.acs.org/doi/abs/10.1021/acs.jmedchem.1c00343

多種多様な研究開発・新規技術を実用化する力のあるメガファーマが治療薬開発に取り組んでいると聞くと、新たなアプローチが可能になるのではないか、と期待してしまいます。

また、逆にベンチャー企業の例として、同様のアプローチを掲げているStablix Therapeuticsという会社があるそうです。 最近できたばかりのようです。*9（シリーズAってなんですか？？）

今回の記事で取り上げた論文の著者らとは別のグループの研究者らのようですが、HPで掲げられているアプローチや標的を見る限り類似の戦略をとっているようです。

興味深いのはCo-Founderの一人、Scott A. Kanner博士が以下の論文のファーストオーサーという点です。（課金してないので読んでません。すみません。）

Kanner, S.A., Shuja, Z., Choudhury, P. et al. Targeted deubiquitination rescues distinct trafficking-deficient ion channelopathies. Nat Methods 17, 1245–1253 (2020). https://doi.org/10.1038/s41592-020-00992-6 Targeted deubiquitination rescues distinct trafficking-deficient ion channelopathies | Nature Methods

こちらの論文はAbstractを見る限り人工脱ユビキチン化酵素の開発のようで、キメラ分子のアプローチとは異なりますが、 メカニズムに知見のある研究者が深く関わっているベンチャー企業がどんな化合物を作り出すのか？、楽しみです。*10

まとめ

以上、PROTACの逆をいくかのようなアプローチ、脱ユビキチン化酵素を利用するキメラ分子によるタンパク安定化の文献紹介でした。 私は専門家ではないので、このアプローチがどれくらい実用可能なものか？全く想像がつきません。*11

ですが、非専門家の私でも「なるほど！面白い！」と思える分かりやすい戦略を打ち出していて、それを実際の化合物構造として実現化している、という点で非常に興味深かったです。 また、異なる複数のグループが同じアプローチを目指しているようなのもびっくりしました。私が不勉強なだけでよく知られた手法なのでしょうか？どんどん競争が激しくなりそうです。

いずれにせよ、低分子（中分子？）化合物でできるアプローチもどんどんと進化していて楽しいですね！

毎度毎度、浅い理解で書いているので間違いが多そうです。ご指摘いただけると幸いです。ではでは！

*1:FDA grants accelerated approval to sotorasib for KRAS G12C mutated NSCLC

*2:Fragments in the clinic: AMG 510

*3:CC-BY-NC-ND 4.0 International license

*4:Wikipedia 脱ユビキチン化酵素

*5:Novartis プレスリリースよりNovartis and UC Berkeley collaborate to tackle 'undruggable' disease targets

*6:この辺りの米国の事情は私は詳しくないので間違っていたらすみません。

*7:ABPPについてはChem-Stationさんの記事「活性ベースタンパク質プロファイリング Activity-Based Protein Profiling」の説明と図が分かりやすいです。

*8:本当に格好いいですよね。尊敬しています。

*9: Stablix Therapeutics Launches with $63 Million Series A Financing

*10:なお、上記論文は今回のbioRxivの文献にも引用されています。むしろ私はこちらの文献からbioRxivに辿り着きました。

*11:POCもまだ初期のin vitro細胞系で、そもそもプレプリントの文献です。

RDKit にコントリビュートした

RDKitの2021.03 versionがリリースされましたね！

お気づきの方はいらっしゃるだろうか・・・

ついに私もRDKit コントリビューターになったよ！

といってもめちゃくちゃ些細な修正をマージして頂いただけです。せっかくなのでご紹介。

最新版のインストール

まずは最新版をインストール。以前と少しコマンドが変わっています。（RDKit documentation installation）

# チャンネルがconda-forgeになってる
conda create -c conda-forge -n my-rdkit-env rdkit

作成した環境「my-rdkit-env」にjupyter notebookもインストール。カーネルを選べるように追加しました。*1

# 仮想環境をactivate
conda activate my-rdkit-env

# jupyterをインストール
conda install jupyter

# カーネルを追加
ipython kernel install --user --name=my-rdkit-env --display-name=my-rdkit-env

Jupyter notebookを起動して確認してみます。

from rdkit import rdBase
print('rdkit version: ', rdBase.rdkitVersion)
#  rdkit version:  2021.03.1

インストールできました！*2

修正内容

マージして頂いた修正はMolFromSequenceに関するものです。MolFromSequenceを使うと、アミノ酸一文字表記の配列（シークエンス）からMolオブジェクトを構築することができます。

以前の記事で取り上げていますのでよろしければご参照ください。

magattaca.hatenablog.com

こちらの配列、オプションの引数「flavor=1」とすることでD-アミノ酸も扱えるのですが、なぜかD-セリンとD-チロシンだけが含まれておらず利用できませんでした。

無いなら付け加えてもらえばいいじゃない！　プルリク！プルリク！

github.com

まずは比較のためL-セリンとL-チロシンを確認しておきます。大文字のSとYでOKです。

from rdkit import Chem
from rdkit.Chem import Draw

L_Ser_L_Tyr = Chem.MolFromSequence('SY')

Draw.MolToImage(L_Ser_L_Tyr)

さて、問題のD-セリンとD-チロシンはどうでしょうか??

flavor=1として小文字のsとyを使えば指定できます。

# D-セリンとD-チロシン
D_Ser_D_Tyr = Chem.MolFromSequence('sy', flavor=1)
Draw.MolToImage(D_Ser_D_Tyr)

先ほどのL体と立体が逆になった構造が生成されました！

以前のRDKitのバージョンでは、これらの配列は空のオブジェクトとなってしまっていました。きちんと動いてよかった！

おわり

以上、重箱の隅をつつくような些細なお話でした。

しょうもなさすぎて蹴られるかと思いましたが、マージしていただけて嬉しかったので記事にしてしまいました。

RDKitコミュニティー優しくて好き！！

*1:Jupyterで複数カーネルを簡単に選択するための設定

*2:Homebrewでもできるそうです。( https://formulae.brew.sh/formula-linux/rdkit )

「ベイズ反応最適化」ツールで遊んだ話 ~EVML, auto-QChem, EDBO~

前回、化学合成のための反応条件最適化ツールの文献をご紹介しました。

Shields, B.J., Stevens, J., Li, J. et al. Bayesian reaction optimization as a tool for chemical synthesis. Nature 590, 89–96 (2021).

doi.org

こちらの文献、データ・コードともに公開してくださっています。

• auto-QChem : Quantum Chemistry
• EDBO : Experimental Design via Basian Optimization
• EVML : Expert versus Machine Leaning

せっかくなので今回はこちらで遊んでみたいと思います。

EVMLのゲームを自分で遊んだ後、EDBOに従って遊んだらどうなるか検証してみます。

反応最適化ゲーム（EVML）

概要

まずは、EVML(Expert versus Machine Leaning)（GitHub）です。

これは、開発された反応条件最適化ツール（EDBO）の性能を比較するためのWebゲームアプリです。
化学者に以下のような条件最適化ゲームをしてもらい、その結果をソフトウェアと比較し、ベンチマークとしています。

Reaction optimization game
ボスから司令のあった化合物がパラジウム触媒による直接アリール化で合成できることがわかった。
パラメータ（base, Ligand, solvent, ｃoncentration, temperature）を変えて、反応収率を最適化したい。
ボスから与えられた期限は1ヶ月。１日にこなせる実験のキャパシティーは5実験（１ batch）。
つまり最大100実験のうちに最適化しなければならない（全組み合わせの約6%）。
使えるのは、自分の知識と、各実験から得た情報のみ（ネットや本をみたらダメ）。
さあ、やってみよう！*1

EVMLの反応探索空間

このゲーム、実際の実験データをつかっているのがすごいですよね。

正直、配位子(Ligand)は色々開発されすぎててよくわかりません。・・・トリフェニルホスフィンの安心感。

塩基(Base)も溶媒(Solvent)も研究室ではみなかったものが多いです。製薬企業さんではよく使われるものなんでしょうか？？？

やってみる

早速あそんでみましょう！

GitHubのファイル（aws_setup.TXT）には、アマゾンウェブサービス（AWS）でのインストール方法が書いてありますが、 よくわからないのでMacで試します。

RのShinyパッケージでつくられたWebアプリのようなので、Rで実行すれば良さそうです。

以下のRのパッケージが必要そうなのでインストールしました。*2

1. shiny
2. dplyr
3. ggplot2
4. reshape2
5. rmarkdown
6. DT

EVMLのGitHubからフォルダをダウンロードしてきたら、「arylation_appフォルダ」の「app.R」ファイルを実行すればOKです。*3

実行すると、アプリが起動し以下のようなページが開きました！

Rで実行するとアプリが開く

左側のパネルに、ゲーム利用者の情報をいれたり、反応条件をいれたりします。
プルダウンから条件を選んで「add」し、５つ条件を選んだら「Run Experiments」をクリック！

選んだ条件の収率(yield)は「Experimet Results」タブに表示されます。 収率を見ながら条件を変えて最適化を目指しましょう！

GitHubの「results」フォルダ「analysis.ipynb」をみると、50人の化学者による結果がまとめられています。

こんな感じです。下図、上が全体で、下３つは属性で分けたサブセットです。

Expertの結果（EvML GitHubの図表より）

全体では、５バッチ前後で最高収率(100%近く)に到達している方が多い印象でしょうか？

所属の違いをみるとfacultyやEngineerは比較的直線的な改善を見せているのに対し、 Med ChemistやProcess Chemistは横ばいとなっているところがあるのが面白いですね。 条件検討の進め方にも考え方の違いがありそうです。

ベイズ最適化ソフト(EBDO)で50回行った結果も紹介されています。以下の通りです。

EDBOの結果（EvML GitHubの図表より）

こちらも５バッチ前後で最高収率に到達しているものが多そうです。ソフトは毎回Maxに到達するまで検討しつづけるのが、ヒトとの違いが見えて面白いですね。

詳細な分析と結果の考察は文献本文をご参照ください。

上から目線でコメントしましたが、私は100%近く(>99%)になるまで7バッチでした。完全にソフトに負けています。
文献で一度結果を読んでいたはずなのに。。。

合成反応の記述子DB（auto-QChem）とベイズ最適化ツール（EDBO）

続いて、ベイズ最適化を実施するためのパッケージauto-QChemとEDBOです。 前者は、関連する化学物質を記述子に変換・保管するデータベースで、後者はベイズ最適化パッケージです。

EDBOのGitHub、examplesやexperimentsディレクトリのipynbファイルには、 論文中のデータや解析がコードとともに掲載されています。ありがたいですね！

以下の記事はだいたいそこからのコピペ切り貼りです。

この記事の目標は、上記の反応最適化ゲームのデータを題材にして各ツールを使ってみることです。 論文の記述子は量子化学計算(Gaussian)によるものですが、Gaussianは使えないしよくわからないので、 より一般的なケモインフォの記述子（Mordred）を用いたいと思います。

auto-QChem

まずは、合成反応を入力に変換し保管するauto-QChemです。MongoDBタイプのデータベースで、ユーザーガイドやインターフェースについては以下で確認できます。

1. auto-QChem GitHub
2. auto-QChem GitHub documentation
3. auto-QChem Database User Guide
4. auto-QChem landing page *4

auto-QChemは名前の通り量子化学（Quantum Chemistry）記述子向けで、計算の実行を簡単にする（入力ファイルの作成と出力の解析）というものです。 Mordredを使うのであれば不要ですが、せっかくなのでさわってみます。

インストール

インストールの仕方はGitHubのinstall.mdに書いてあります。

condaで仮想環境を作って、関連するパッケージをインストールした後、GitHubから落としてきたauto-QChemのsetup.pyを使えばOKです。

# 新しい仮想環境を作ってactivate
conda create --name autoqchem python=3.7
conda activate autoqchem

# パッケージのインストール
conda install jupyter pandas scipy matplotlib pymongo pyyaml fabric xlrd appdirs openpyxl
conda install -c conda-forge openbabel=2.4.1 # v.3.0.0はよくないらしい
conda install -c rdkit rdkit
python -m pip install imolecule==0.1.13   #3dで描画するのに使う

ついで、git cloneなどでauto-QChemをローカルにもってきたあと、保存したディレクトリに移動し、setup.pyを実行します。

python setup.py install

おしまい！*5

分子の取り扱い（moleculeクラス）

auto-QChemでは分子はmoleculeクラスで扱います。OpenBabelのOBMolクラスを継承したクラスになっています。

デフォルトではmoleculeクラスの入力は分子のSMILESで、３次元構造(conformer)を構築します。 コンフォメーションはOpenBabelのOBConformerSearch()で、遺伝的アルゴリズムを使って生成されているようです。

試しにリガンドを一つ読み込んでみましょう！

論文でフィーチャーされてたCgMe-PPhというリガンドです。PubChemではこちらの名前では引っかからずmeCgPPhがシノニムにありました（PubChem CID: 53427790）

from autoqchem.molecule import molecule

# Cg-MePPhのSMILESから分子を構築
smiles_str = "CC12CC3(OC(O1)(CC(O2)(P3C4=CC=CC=C4)C)C)C"
mol = molecule(smiles_str)
print(type(mol))
# <class 'autoqchem.molecule.molecule'>

draw()という関数が用意されており、jupyter notebook上でぐりぐり動かせる3D描画が見られます。

# mol.draw()

かわいい。PubChemと描画スタイル同じですね。PubChemもimolecule使っているんでしょうか？

コンフォマーを複数発生させることもできます。また、loggingでログを表示させられるので、動作確認しておきたいときは表示させても良いかもしれません。

# ログの表示
import logging
logging.basicConfig(level=logging.INFO)

# 最大30のコンフォマー作成
mol_conf = molecule(smiles_str, max_num_conformers=30)

# INFO:autoqchem.molecule:Initializing molecule with canonical smiles: CC12CC3(C)OC(P2c2ccccc2)(CC(O1)(O3)C)C
# INFO:autoqchem.molecule:Creating initial geometry with option 'best'.
# INFO:autoqchem.molecule:Initial geometry created successfully.
# INFO:autoqchem.molecule:Conformer Search generated 6 conformations of CC12CC3(C)OC(P2c2ccccc2)(CC(O1)(O3)C)C molecule

６つのコンフォマーができたみたいです。comformer_numで各コンフォマーにアクセスできます*6。

違いがわからない。。。

その他の使い方や、Gaussianのインプットファイルの作成方法等はFunctional documentationをご参照ください。 また、データベースのインターフェースはFlaskを使っているようですがよくわからないので飛ばします。

Mordredで記述子計算

ベイズ最適化パッケージに進む前に、入力の記述子を用意しておきたいと思います。必要な分子（Ligand、Base、Solvent）のSMILESは以下です。

Ligand_smiles = ["CC(C)C1=CC(C(C)C)=C(C(C(C)C)=C1)C2=C(P(C3CCCCC3)C4CCCCC4)C(OC)=CC=C2OC", "CC(C)(C)P(C1=CC=CC=C1)C(C)(C)C", "CN(C)C1=CC=CC(N(C)C)=C1C2=CC=CC=C2P(C(C)(C)C)C3=CC=CC=C3", "P(C1CCCCC1)(C2CCCCC2)C3CCCCC3", "P(C1=CC=CC=C1)(C2=CC=CC=C2)C3=CC=CC=C3", "CC(C1=C(C2=CC=CC=C2P(C3CCCCC3)C4CCCCC4)C(C(C)C)=CC(C(C)C)=C1)C", "P(C1=CC=CO1)(C2=CC=CO2)C3=CC=CO3", "CP(C1=CC=CC=C1)C2=CC=CC=C2", "CC(OC1=C(P(C2CCCCC2)C3CCCCC3)C(OC(C)C)=CC=C1)C", "FC(F)(F)C1=CC(P(C2=C(C3=C(C(C)C)C=C(C(C)C)C=C3C(C)C)C(OC)=CC=C2OC)C4=CC(C(F)(F)F)=CC(C(F)(F)F)=C4)=CC(C(F)(F)F)=C1", "C[C@]1(O2)O[C@](C[C@]2(C)P3C4=CC=CC=C4)(C)O[C@]3(C)C1", "CP(C)C1=CC=CC=C1"]
Base_smiles = ["O=C([O-])C.[K+]", "O=C([O-])C(C)(C)C.[K+]", "O=C([O-])C.[Cs+]", "O=C([O-])C(C)(C)C.[Cs+]"]
Solvent_smiles = ["CCCCOC(C)=O,CC1=CC=C(C)C=C1", "CCCC#N", "CC(N(C)C)=O"]

それぞれMordredで記述子を計算した後、エラーを取り除きcsvで書き出します。auto-QChemのnotebook Fast featurize with Mordred.ipynbを参考にして進めていきます。

まず、必要なライブラリをインポート。

import pandas as pd
import numpy as np
import mordred
from rdkit import Chem
from mordred import Calculator, descriptors

数が少ないSolventを例に進めます。PandasのDataFrameに記述子を格納したいので、SMILESのlistをDataFrameにしておきます。 ついでにrdkitのrdmolオブジェクトに変換して列を加えておきます。

solvent_df = pd.DataFrame(Solvent_smiles, columns={'smiles'})
solvent_df['rdmol'] = solvent_df['smiles'].map(lambda x: Chem.MolFromSmiles(x))

print(solvent_df.shape)
# (4, 2)

ModredのCalculatorで記述子を計算します。

calc=Calculator(descriptors, ignore_3D=True)
md=calc.pandas(solvent_df['rdmol'])
print(md.shape)
# (4, 1613)

1613個の記述子が計算されました。エラーが出ている部分をNaNに変換して、dropnaで列を削除します。

# エラーをNaNで置き換え
md=md.applymap(lambda x: np.nan if type(x) in [mordred.error.Missing, mordred.error.Error] else x)

# NaNを含む列の削除
md=md.dropna(axis=1)

print(md.shape)
# (4, 1275)

1275個残りました。あとはcsvとして保存するだけです。

SMILESのDataFrameと記述子のDataFrameをconcatenateで結合します。rdmolオブジェクトの列は不要なので削除しておきます。 また、他の要素と区別するため、列名にprefix 「solvent_」をつけておきます。

# 保存用のDataFrameを作成。rdmolは削除しておく。
df_to_save = pd.concat([solvent_df.drop('rdmol', axis=1), md], axis=1)

# prefixをつける
df_to_save = df_to_save.add_prefix('solvent_')

# csvに保存
df_to_save.to_csv("solvent_desc.csv", index=False)

無事「solvent_desc.csv」ファイルが作成できました！

同様にして「ligand_desc.csv」「base_desc.csv」を作成しました。これで記述子の準備はおしまい。

ベイズ最適化パッケージ（EDBO）

ではいよいよ本命、EDBO(Experimental Design via Basian Optimization)です。

これはベイズ最適化を化学合成に適応するための実用的なパッケージです。

• GitHub : EDBO GitHub
• ドキュメント : edbo documentation

インストール

まずはGitHubのREADMEに従ってインストールします。

EDBO自体はpipでインストールできます。必要なライブラリはrdkit、Mordred、PyTorchのようです。

指示に従ってcondaで新しい仮想環境（edbo）を作成し、諸々インストールしました。

# https://github.com/b-shields/edbo
# (0) Create anaconda environment
conda create --name edbo python=3.7.5

# (1) Install rdkit, Mordred, and PyTorch
conda activate edbo
conda install -c rdkit rdkit
conda install -c rdkit -c mordred-descriptor mordred
conda install -c pytorch pytorch=1.3.1

# (2) Install EDBO
pip install edbo

# Running Notebooks
conda install jupyterlab

おしまい。

計算の流れとモデルの中身

EDBOで反応を最適化する流れは以下の通りです。GitHubの例（Bayesian Optimization of a Mitsunobu Reaction）を参考にしています。

1. 反応をエンコード、入力の記述子を取り込み：edbo.utils
2. 反応空間（reaction space）のパラメーター設定: components``encodings``descriptor_matrices
3. ベイズ最適化の設定（ebdoインスタンス化）: edbo.bro.BO_express
4. 計算ループの実施　　
   1. 初期化（最初の実験セットを提案）: init_sample()
   2. 実験結果の追加 : add_results()
   3. モデルの適合、獲得関数の最適化、次の実験の提案 : run()
   4. 提案実験の出力 : export_proposed()
   5. 提案実験の結果を追加し、収束するまで 2 -> 4を繰り返す

上記手順「3. ベイズ最適化の設定」でedbo.bro.BO_expressというオブジェクトを利用しています。 これは、パラメーターの自動的な読み込みや、全実験空間の組み合わせの計算、必要な前処理の実施といった作業を行ってれるもので、 edboをよりユーザーフレンドリーなプログラムにしてくれています。

BO_expressはベイズ最適化のモデルもデフォルトで設定してくれています。デフォルトの設計は以下となっています。(class edbo.bro.BO_express)

• サロゲートモデル : ガウス過程(Gaussian Process)
  • カーネル : Matérn5 (カーネルパラメータ ν = 5/2 (2.5))
  • ハイパーパラメータの決め方 : 最尤推定(MLE)
• 獲得関数 : EI (Expected Improvement)

ガウス過程はGPyTorchを使って実装されています。*7関連するやつを貼っておきます。*8*9

カッコイイやつ

では順番に計算していきましょう！

EDBOでEVMLをプレイ

今回は「EDBOの指示に従ってEVMLをプレイする」という体で、計算を実行したいと思います。つまり、EDBOの計算で提案された実験をEVMLで実施し、結果をEDBOにフィードバックして次の実験を再度提案させ、収率の向上を目指します。

ではまず反応のエンコードです。 入力として使う分子の記述子を読み込みます。edbo.utilsを使って、Mordredの計算結果を含むcsvファイルを読み込みます。

計算実行

まず反応のエンコードです。 入力として使う分子の記述子を読み込みます。edbo.utilsを使って、Mordredの計算結果を含むcsvファイルを読み込みます。

# インポート
import pandas as pd
from edbo.utils import Data

# 記述子計算のファイルを読み込み
ligands_desc = Data(pd.read_csv('data/ligand_desc.csv'))
solvents_desc = Data(pd.read_csv('data/solvent_desc.csv'))
bases_desc = Data(pd.read_csv('data/base_desc.csv'))

読み込んだデータはdataで確認でき、PandasのDataFrameとなっています。

print(type(ligands_desc))
print(type(ligands_desc.data))
# <class 'edbo.utils.Data'>
# <class 'pandas.core.frame.DataFrame'>

次に、反応空間(reaction space)のパラメーター設定を行います。今回の合成反応条件で最適化するパラメータ（concentration、temperature）を含めて要素を指定する辞書を作成します。

# 反応空間のパラメーター
components = {'Ligand':'mordred', 'Solvent':'mordred', 'Base':'mordred',
              'Concentration' : [0.057, 0.100, 0.153],
              'Temperature' : [90, 105, 120]}

# Encodingの指定。指定しない場合、自動的にOne-Hot-Encoding (OHE)になる。
encoding = {'Concentration' : 'numeric', 
            'Temperature' : 'numeric'}

# 記述子の設定 (descriptor_matrices) 
mordred = {'Ligand': ligands_desc.data,
           'Solvent' : solvents_desc.data,
           'Base' : bases_desc.data}

descriptor_matricexの設定はEDBO外部で計算した記述子を使う場合に設定する項目で、Gaussianで計算したDFT記述子の場合などに使います。Mordredの計算はEDBOでも行えるので、実際にはここで指定しなくてもよかったのですが、今回は前もって計算したデータを使うことにしました。

次に、EDBOの BO_expressオブジェクトをインスタンス化します。上で指定した要素やモデルの設定など、もろもろを書き込みます。

from edbo.bro import BO_express

# BO object

bo = BO_express(reaction_components = components,
                encoding=encoding,
                descriptor_matrices=mordred,
                acquisition_function='EI',    # 獲得関数はEI
                init_method='rand',            # 初期化はランダム
                batch_size=5,                    # １ラウンドあたりの実験数は５
                target='yield' )                   #  yieldを最適化

# 事前分布の指定。文献の設定に合わせる

from gpytorch.priors import GammaPrior

bo.lengthscale_prior = [GammaPrior(2.0, 0.2), 5.0]
bo.outputscale_prior = [GammaPrior(5.0, 0.5), 8.0]
bo.noise_prior = [GammaPrior(1.5, 0.5), 1.0]

下準備が完了したので計算ループに進みます。

まずは初期化！最初にする実験をランダムに選択します。

選択された実験はbo.export_proposed()でcsvファイルとして書き出すことができ、またnotebook上でもbo.get_experiments()で確認できます。

bo.init_sample(seed=0)             # 初期化
bo.export_proposed('data/init.csv')     # CSVファイルに書き出し
bo.get_experiments()               # ノートブック上で確認

ランダムに選ばれた初期反応条件

５つ実験が選ばれました！この条件で実験を行い、結果を追加して次のループに進むことになります。

実際の環境では、実験に時間がかかるのでPCを付けっ放しにしておくのはもったいないです。そのような場合に設定を保存しておけるよう、EDBOではpickel化するメソッドedbo.bro.BO.save()も用意されています。

bo.save()  # optimizerの保存

作業ディレクトリに「BO.pkl」が作成されました。再度読み込みたいときはedbo.bro.BO.load()を使えばOKです。

from edbo.bro import BO_express

bo = BO_express()
bo.load()

さて、話戻して、出力されたcsvファイルを見ると、以下のようにyield列があり「＜Enter Response＞」となっています。ここに収率を書き込めば、次のループで使用する入力ファイルとできます。

出力csvにyieldを書き込む

今回取り上げている反応の結果は、EBDOのGitHub 「experiments/data/direct_arylation/experiment_index.csv」に記載されています。もしくはEVMLで確認することもできます。

今回は「EDBOの指示に従ってEVMLをプレイする」ということにしているので、EVMLに条件を入力し、得られた結果をcsvに追記し、次のラウンドに回しました。

結果を記入したファイルはedbo.bro.BO_express.add_results()メソッドで読み込みます。ついでrun()を実行すれば、ガウス過程モデルのfitと獲得関数(EI)の最適化を行い、次にすべき実験を提案してくれます。

# 結果をboオブジェクトに追加
bo.add_results('data/init_results.csv')

# 計算の実施
bo.run()

# 次にする実験の取り出し
bo.export_proposed('data/round0.csv')  
bo.get_experiments()  

新しい提案の実験条件が手に入りました（round0.csv）。これをEVMLに入れて結果を得て、再度EDBOで計算、次の条件を提案を繰り返します。

・・・結果・・・・

round6の提案実験で、収率100%の反応条件が２つ見つかりました！

最初のinitバッチを含めて計7バッチ（35実験）です。だめ押しでround7も実験しましたが、収率は下がったのでここで止めました。

EVML上で確認できる「条件探索の過程」は以下の図のようになりました。

右から左に向かって新しい実験です。赤色がEBDOの提案した実験の収率で、緑色はランダムに選んだ場合の実験収率です。

最初は10%にも満たなかった収率が徐々に上がってゆき、最終的に収率100%を達成している様が見て取れます。

結果の確認

EDBOは結果解析用のツールも提供しています。

例えば、最適化ループを回す過程でoptimizerのconvergenceがどのように変化したかプロットすることができます。

bo.plot_convergence()

順調に収束している様子がわかりますね！

GitHubのnotebookで使われている解析も適用してみましょう！

実施した全40実験（8バッチ）の中で、最も収率の良かった上位5実験について反応条件を提示してみます。

全ての結果を読み込んで統合したDataFrameを作成します。カラムは反応条件（Ligand, Solvent, Base, Concentration, Temperature）と収率（yield）で、yieldで降順に並べ替えます。

results = pd.DataFrame(columns=bo.reaction.index_headers + ['yield'])
for path in ['init', 'round0', 'round1', 'round2', 'round3', 'round4', 'round5', 'round6', 'round7']:
    results = pd.concat([results, pd.read_csv('data/' + path + '_results.csv', index_col=0)], sort=False)

results = results.sort_values('yield', ascending=False)
results.head()

できました！

上位５つのリガンド、塩基はどれも同じですね。温度はより高温、濃度は濃い方が良さそうですが、溶媒の影響もありそうです。

DMAc溶媒では薄め(0.057)でも収率96.6%なのが興味深いですね。

棒グラフも書いて見ます。

import matplotlib.pyplot as plt

fig, ax = plt.subplots(1, len(results.columns.values[:-1]), figsize=(30,5))

for i, feature  in enumerate(results.columns.values[:-1]):
    results[feature].iloc[:5].value_counts().plot(kind="bar", ax=ax[i]).set_title(feature)
plt.show()

また、EDBOは構造式描画のためのChemDrawというモジュールも用意されています。

上位５反応に出てきた分子の構造式を描画しましょう。

from edbo.chem_utils import ChemDraw

for col in results.iloc[:,:3].columns.values:
    print('\nComponent:', col, '\n')
    cdx = ChemDraw(results[col].iloc[:5].drop_duplicates())
    cdx.show()

Component: ligand_smiles_index

Component: solvent_smiles_index

Component: base_smiles_index

なるほど便利！！

やっぱり構造式でみるのがが一番しっくりくる！！

感想

以上、今回は前回の記事で取り上げたNatureの文献のツールで遊んで見ました。データとツールが公開されているというのはとてもありがたいですね！

特にベイズ最適化ツールEDBOは変更が必要な設定項目も少なく、コピペするだけで簡単に使えてとても便利な印象でした。

EDBOの実践として、記述子にMordredを使い、EVMLをプレイして見ました。最初は低い収率で類似の実験を検討しており、「大丈夫かなー、本当にこれでいいのかなー」と半信半疑でしたが、収率上昇の兆候が見え始めると一気に条件検討が進んでいき最終的に収率100%の条件を見つけることができました！

私は自分でEVMLをプレイした際、99%までしか行かなかったので、100%の結果を見てびっくりしました。・・・完敗です。

ところで、文献中では最終的に記述子として量子化学計算によるものが用いられていました。今回はMordredでも良い結果に到達できることが確認できたので良かったです。Gaussianなんて一生触れそうそうにないですからね！

溶媒の沸点を考慮して反応温度に拘束をかけるとか、追加の設定を加えたらより探索効率上がるかもなー、とか思いましたがその辺りどうなんでしょう？？？

いずれにせよ、面白いツールがどんどん出てきて公開されるのは見ていて楽しいですね！

ほとんどGitHubのコードそのままコピペした記事ですみませんでした。。。ではでは！

*1:READMEをそれっぽっくかえてみました

*2:参考「aws_setup.TXT」の「5. Install required packages」

*3:私はRStudioで実行しました。Rはほとんど触ったことがなかったのでRStudioの使い方に手間取りました・・・

*4:ランディングページ（LP)というのは、ユーザーがWebで最初にアクセスするページのことらしいです。知らなかった。

*5: 仮想環境でjupyter notebookの立ち上げにエラー（Bad config ~~）が出るときは、「pip install environment-kernels」をすれば動くようになるかもしれません。参考

*6:コンフォマーの数以上の数字を与えると、最後のコンフォマーが選ばれるそうです

*7:GPyTorchについてはHELLO CYBERNETICSさんの記事「GPyTorchでガウス過程を見てみよう」がとても参考になります。

*8:参考 講談社「ガウス過程と機械学習」

*9:参考ガウス過程回帰の基礎 赤穂昭太郎 システム/制御/情報 2018(62)390

「ベイズ反応最適化」(Nature 2021)がわからないのでベイズ最適化をお勉強したメモ

先日（2021/2/3）公開のNatureにとても気になる文献がでました。

Shields, B.J., Stevens, J., Li, J. et al. Bayesian reaction optimization as a tool for chemical synthesis.
Nature 590, 89–96 (2021).

doi.org

News&Viewsは無料で読むことができ、研究の概要をつかむことができます。

Machine learning made easy for optimizing chemical reactions Hein, J.E. Nature 590, 40-41 (2021)

doi.org

要するに、「有機合成の反応条件の最適化に、機械学習のツール（ベイズ最適化）を適用することで、 実験化学者を上回る結果がでた」という内容のようです。

反応条件の最適化は、実験条件（試薬・温度・濃度・・・）の組み合わせが無数にありえ、化学者の経験と感に寄るところも大きいです。 この問題を改善しうる、より効率的なソフトウェアが構築され（オープンソース！）、しかも得られた結果の中には、 経験のある化学者がだれも思いつかなかった条件でより高い収率を達成したものもあった、とのこと。

「こりゃ、おまんまの食い上げだよ！」ってことで、思わず課金してしまったのですが、読んでもさっぱりわからなかったのです。

というわけで、必要そうな周辺知識のお勉強メモ。。。

論文の簡単な紹介

まずはざっと論文の概要をご紹介します。

著者グループの過去文献(Science 2018(360)186)

さて、今回の文献はPrinceton大学のAbigail G. Doyle教授のグループとBristol-Myers Squib社によるものです。

DoyleグループのHPはこちら。 おもに有機合成反応開発を研究テーマとされているようですが、 ハイ・スループット実験（HTE）や機械学習(ML)といった手法を積極的に取り入れて成果を挙げられている、格好良いグループです。

関連する報告として、2018年Science誌にMerck社との共著を出されています。（本文公開されています。）

Predicting Reaction Performance in C-N Cross-Coupling Using Machine Learning.
Ahneman, D. T.; Estrada, J. G.; Lin, S.; Dreher, S. D.; Doyle, A. G. Science 2018, 360, 186-190.

science.sciencemag.org

Science文献では、有機合成反応の中でも必要な試薬が多く、パラメーターの多いクロスカップリング反応（Buchwald-Hartwig cross-coupling) を題材として、反応阻害物質の収率への影響を予測するモデルを構築しています。

モデルのアルゴリズムはランダムフォレスト(RF)ですが、入力となる化学物質の特徴量として、DFT計算で得た値（HOMO、LUMO etc.）や振動モード (vibrational mode)を記述子に取り込んでいる点が興味深いです。

この量子化学計算による記述子というアプローチはNature文献でも引き継がれています。

Nature文献の概要

今回の論文のポイントはだいたいこんな感じ。。。

1. ベイズ最適化の有機合成実験への実践的適用
2. 化学者が利用可能なオープンソースフレームワークの作成
   • コンピューター上でのハイスループット実験と反応の特徴量化(auto-QChem)
   • ベイズ最適化パッケージ(EDBO)
3. BMS社のハイスループット実験(HTE)データを利用した性能比較（ベンチマーク）
   • Webゲームアプリを作成し、化学者 vs. ソフトウェアの比較実験(EVML)

上記、データ・コードともに公開されており、文献ページでもリンクを確認できます*1。

それぞれ名称が何の略か見ると機能がわかりますね。

• auto-QChem : Quantum Chemistry
• EBDO : Experimental Design via Basian Optimization
• EVML : Expert versus Machine Leaning

反応の特徴量化(featurization)としては、ケモインフォマティクスの記述子(Mordred)やOne-hot-encoding(OHE)の結果と比較した上で 量子化学(DFT)を選択しています。*2

また、性能評価のために、合成反応最適化ゲームを作って化学者とソフトを競わせてしまう（Expert versus Machine Leaning）というのは面白いアプローチですね。

具体的検証として、複数の合成反応における比較結果が記載されています。 一例として、クロスカップリング反応において、機械学習は化学者よりも高い収率に到達し、なんとその条件で使われている配位子は同種の反応で報告例のないもの（!?)、だったそうです。

感想

ベイズ最適化自体は、実験計画法（Design of experiments, DOE）との比較文脈で、 化学以外の他の分野でも多く使われているそうなので、数理分野に親しみのある方からみれば、あまり複雑なことはしていないのかもしれません。 ですが、企業と組んで実際の実験データを利用し、実践的なシステムの構築、ゲーミフィケーションのとりいれといったアプローチを実現しつつ、 実際にヒト以上の性能のソフトウェアを作ってしまっているのは、本当にすごいなーと思います。

感想文はここまでにして、、、ベイズ最適化のお勉強に進みたいと思います。

目標は、文献でわからなかった以下の用語たちのお気持ちを知ることです。

1. gaussian process surrogate model
2. Matérn52
3. explore / exploit
4. aquistion funtion
5. Thompson Sampling
6. Kriging believer algorithm

ベイズ最適化について

ベイズ最適化の雰囲気

ベイズ最適化についてはSlideShareの「機械学習のためのベイズ最適化入門」が雰囲気をつかみやすかったです。*3

まとめスライドから引用します。

• ベイズ最適化は高コストなブラックボックス関数の最適点を効率よく求める手法
• ガウス過程を用いて目的関数を表現
• 獲得関数が最大となる点を次の観測点に選ぶ

だいたいこんな感じの話みたいです。*4

まず、ベイズ最適化は「形状のわからない関数（ブラックボックス関数）」について、「効率的(少ない評価回数)」に「最大値または最小値」を求めるための手法です。

ここで解きたい問題は「大域的最適化問題」。つまり「局所解」(local)に陥らずに「大域解」(global)を求めるのが目標です。

「効率的」とありますが、逆に効率が悪い例はグリッドサーチのように探索範囲のすべての組み合わせを検討するアプローチや、ランダムサーチのように完全にランダムに探索するものです。

これに対して、ベイズ最適化では「評価値の良そうなエリアを優先的に観測」しつつ、 「局所解に陥らないように（観測点が）スカスカのエリアもたまに観測」することで、より少ない回数で最適解に到達することができます。

この「次に観測する観測点」は、「これまでの結果をもとに」決められます(逐次最適化法)。人間が経験に基づいて直感的に判断するのを再現するようなものとなっています。

その判断の戦略には複数あり、例として① PI戦略 (Probability of Improvement) 、 ②EI戦略 (Expected Improvement) 、③UCB戦略 (Upper Confidence Bound) といったものがあります。

いずれの戦略においても、「何らかの関数を最大化する点を次に観測する」という選択をします。 このような関数を「獲得関数(acquisition function)」とよびます。*5

より具体的には、獲得関数は期待値(μ)や偏差(σ)の組み合わせで定義されます。 これらの値は、「最適化したい関数がガウス過程に従うと仮定する」ことで求めることができます。 つまりガウス過程の導入により獲得関数の最大化が計算可能となります。

ガウス過程では、形状のわからない関数（ブラックボックス関数）のだいたいの形を観測点から予想します。すると「未観測点の期待値(μ)と分散(σ²)」が算出可能となります。

• 期待値(μ)が大きい → 周囲の観測点が大きい → 評価値の良さそうなエリア
• 偏差(σ) が大きい　→　周囲が観測されていない → スカスカのエリア

となるので、先に見たベイズ最適化の探索の仕方に対応していることがわかります。

こうして決めた次の観測点のデータが得られたら、そのデータを取り込んでモデルを更新し、どんどんと精度を高めていきます。*6

以上がベイズ最適化の雰囲気です。

Nature文献の例では？

今回のNature文献の例では、以下のようなベイズ最適化の問題だったといえるでしょうか？

• 反応条件(xi)の組み合わせ(x1, x2,… xn)と収率(y)の関係は具体的な形がわからない関数
• 求めたいのは最大の収率(y)を与える最適な反応条件(xi)の組み合わせ
• 膨大な反応条件の組み合わせから、効率的に大域解をあたえる組み合わせを探索したい。

また、「実験化学者が見出せなかった条件にソフトウェアが到達した」という結果は、

• 化学者は経験や類似の反応の報告例をもとにあたりをつけた、良さそうな条件（期待値大）の周りを探索するため、 少ない試行回数である程度良い結果を得ることができる。
• 一方で、ソフトウェアは期待値が低く偏差が大きい組み合わせも探索に含めるため、報告例のない大域解にも到達することができた。 また、そこに至る試行回数も現実的に許容可能な範囲であった。

と、いった内容になりそうです。

ベイズ最適化とガウス過程(Gaussian Process)

ベイズ最適化の雰囲気がわかりました。気になるのは、未知の関数の計算のためにとりいれられた「ガウス過程の仮定」という、天丼みたいな用語です。もう少しつっこみましょう。

サロゲートモデル？

さて、先にベイズ最適化は逐次最適化を行うとかきました。 このように「反復的に関数評価と目的関数のモデルの更新を繰り返す手法」をSequential Model-based Optimization(SMBO)と呼ぶそうです。

また、この「目的関数のモデル」のことをサロゲートと呼びます。これらの用語を使うと、「ベイズ最適化は、サロゲートをベイズ的に構築するSBMO」ともいえます。*7

サロゲートは先の「形状のわからない関数」のモデルにあたりますが、その種類はいくつかあり、ガウス過程、ランダムフォレスト、Tree Parzan Estimator(TPE)などがあります。

Nature文献中のGaussian process surrogate modelという記載はこれだったんですね！

ガウス過程を５分で理解！？

では、ガウス過程(Gaussian Process, GP)はどのようなものなのでしょうか？

本棚見たらこんな書籍をもってました！なんと冒頭、第0章のタイトルは「たった５分でガウス過程法がわかってしまう」！！！ まじかーー！！

ガウス過程と機械学習 (機械学習プロフェッショナルシリーズ)

ガウス過程は「観測データをもとに、元の関数（f(x)）の形を推定」する方法の一つです。 「回帰問題（データから誤差を織り込んだ入出力(x-y)関係を推定する問題）」の解法の一つ、というわけですね！

教科書でよく見る回帰問題の標準的な解法は最小二乗法(OLS)です。最小二乗法では誤差を最小とするように、最適なパラメーターをビシッ!と決めてやります 対してガウス過程回帰法では、ビシッと一つの形として推定するのではなく、あやふやさを織り込んだ広がりを持ったもの（「確率分布の雲」）と推定するのが特徴です。

   

また、この関数の推定（モデル化）においてガウス分布を道具として使います。複数のガウス分布を組み合わせることで、多様な関数の形を表現できるので「ある程度の滑らかさを持つ関数だったら何でもアリ」という広い条件のもとで探索することができます。  

念のため、「ビシッ = 直線」「あやふや = 曲線」という意味ではないです（上図、OLSが１次関数なのは手頃な図がなかったからです）。ここで言いたいのは、「ベイズ主義でいくぜ！」ってことです。

頻度主義では「知りたいパラメータは1つの真の値」（ビシッ）でした。一方でベイズ主義では「知りたいパラメータは確率変数である」（雲）と考えます*8。ベイズ推定によりパラメータを推定し得た関数の形は確率的なものなので、このモデルによる予測も広がりのあるものとなります。

観測データはおおよそ予測値（期待値）のあたりだけどズレるのは織り込み済み。予測の確度が高い部分は広がり（分散）が小さく、自信がなければ大きくなる。と、いうのが「あやふやさを織り込んだ広がり」とした推定のようです。

ベイズに話が逸れましたが、何はともあれガウス過程のガウスがガウス分布に由来しそうなことがわかりました。数学っぽくするためにWikipediaを引用しておきます。*9*10

では、ガウス過程の残りの「過程」とは何でしょうか？

これはガウス過程が「「確率過程」(stochastic process)の一種である」ことに由来するそうです。 確率過程とは、「入力の集合(x_1, x_2, ..., x_N)に対応する確率変数の集合(y_1, y_2, ..., y_N)に同時分布p(y_1, y_2, ..., y_N)を与える確率モデル」です。*11

また、観測Yの確率モデル（あるいは確率的生成モデル (probablistic generative model)）とは、「現実世界における観測Yは、何らかの確率分布p(Y)からのサンプリング Y ~ p(Y)によって得られたものだ」とする仮説のことです*12。サンプリングと考えることで、「観測される値の生成過程を確率分布で表現」することができます*13。

ちなみに、もともと「時系列に対する理論」として生まれたので「過程」という名前がついていますが、時系列以外の一般に成り立つ理論だそうです。 *14また、確率過程が仮説であることから、「ガウス過程を仮定」という表現の「仮定」の意味もわかりますね。

格好いい表現を書籍からお借りします。

さて、だいたい必要そうな用語をひろいました。これらを踏まえて「ガウス過程と機械学習 第０章」の表現をお借りしましょう。

簡潔にいうと、ガウス過程とは「関数f(x)を確率変数と見立てた確率分布」であり、ガウス過程回帰とは「データから関数f(x)の確率分布をガウス過程の形で求める方法」です。 また、求めたいものは「関数f(x)の「事後分布の雲」」であり、ベイズ推定で求めます。*15

ガウス過程の定義を以下に引用しておきます。

上記でN次元としていますが、Nは無限に大きくても良いのでガウス過程は「無限次元のガウス分布」です。

定義の用語を使えば、ガウス過程回帰は、平均 μ と共分散行列 K からfの事前分布であるガウス過程を定め、観測データによってfの事後分布を求める、ということになります。

以上、ガウス過程の用語でした。「ガウス過程と機械学習」の色々な章から継ぎ接ぎしたので、正確なお話は書籍を買いましょう！図表が多く語り口も優しいので数学嫌いでも取りつきやすいですよ！*16

ガウス過程とカーネル法

さて、ベイズ最適化の関数のモデルにおいて使われるガウス過程の雰囲気がわかって来ました。ガウス過程の定義を参照すると、関数の形を求めていく作業は、平均 μ と共分散行列 K からfの事前分布であるガウス過程を定め、観測データによってfの事後分布を求める、作業になります。

この問題を解くために、より具体的な形に入っていくにはカーネルという言葉を押さえておく必要があります。サポートベクターマシーンでよく見る格好いいやつですね！

具体的には共分散行列がカーネル関数から計算されるカーネル行列となっています。カーネル関数をつかうことで、計算量をずっと減らすことができます(カーネルトリック)。*17

カーネル関数の定義には様々あり、以下の書籍には複数のアプローチからの導入やアプローチ間の関係性が議論されています。*18

カーネル多変量解析―非線形データ解析の新しい展開 (シリーズ確率と情報の科学)

このうち、確率モデルからの導入において、線形モデルのヘイズ推論から出発して、ガウス過程とカーネル関数の関係性が説明されています。また、「ガウス過程と機械学習 Chaptor 3」において、線形モデルからパラメータ(重み w)の消去によりガウス過程を導く流れからも説明されています。*19

だいたいこういう雰囲気の話です。

上の流れで見たことから、カーネル関数のグラム行列を分散共分散行列とみなすと、カーネル関数とガウス過程は等価となります*20。カーネル関数をつかうと、特徴空間における内積をデータから直接計算することができるので、特徴ベクトルを明示的に使った計算が不要となり、計算量を減らすことができます。*21

さて、カーネル関数の概略がだいたいわかりました。では、具体的にはどのような形の関数なのでしょうか？

カーネル関数の形

カーネル関数の具体例に入る前に、その意味合いについて触れておきます。

カーネル関数は「サンプル間の類似度を評価するもの」ということができます。カーネル関数は共分散行列の要素を与えており、多変量ガウス分布において2つの変数の間の共分散となります。「共分散が大きい = 似た値をとりやすい」ということになるので、「入力 xとx'の近さ」（類似度）をあらわす意味合いがある、ということになります。*22

カーネル関数の満たすべき条件は書籍を参照していただくとして、具体的な例を以下に挙げます。*23

一番下、Matérnカーネルという複雑な形を引用しました。Matérnカーネルのパラメータ ν は関数の滑らかさを表すもので、カーネルの性質を決めます。 「ν = ∞」のとき、RBFカーネルに一致し、「ν = 3/2」はMatérn3、「ν = 5/2」はMatérn5と呼ばれます。

今回題材にしているNatureの文献ではMatérn52という記載があり、カーネル関数としてMatérn5（ν = 5/2）を使用していると言う意味だと思われます。

以上、カーネル関数の形についてでした。

これで、ガウス過程の概要と、具体的な式表現をざっと拾い終わりました。 先に見た通り、ガウス過程はベイズ最適化のサロゲートモデルの一つで、目的関数の表現に関わるものでした。 では続いて、ベイズ最適化の残り、効率的な探索部分に関して見て行きましょう。

探索と活用と獲得関数

さて、ベイズ最適化の概要でみたように、べイズ最適化では効率的な探索のため、獲得関数を最大化するように次の観測点を選んでいくのでした。

探索と活用のトレードオフ

このような、観測を繰り返して最適解を探す問題において効率化のために考えるべき点として「探索と活用のトレードオフ」(Exploration-exploitation trade-off)があります。

（しらみ潰しやランダムでない）効率的な探索には、すでにある情報（過去の経験）の活用が鍵となります。 やってみて良い結果になった行動を続けていけば、同様に良い結果が得られる可能性が高くなります（経験の活用(Exploitation)）。

一方で、活用ばかりでは、離れた場所にある、もっと良い結果につながる行動を見つけることができません。

そこで、不確実な未経験の行動についても試してみることが大事になります（探索(Exploration)）。ですが、探索ばかりでも、今度は過去の経験を生かせない、と言う問題があります。*24

このように、探索と活用にはトレードオフがあるので、両者のちょうど良いバランスを取ることが目指されます。このバランスをどのように取るか？バランスの取り方（戦略）を反映するのが獲得関数(acquitson function)となっています。

獲得関数

では、獲得関数には具体的にどのようなものがあるのでしょうか？ 以下の３つが有名みたいです。*25 *26 *27

上記、①PI戦略の獲得関数は改善確率に基づくもので、改善量を考慮していません。 従って、「改善量は大きいが不確実性が高い」選択肢があったとしても、「改善量は小さいが確度は高い」選択肢が選ばれる可能性があります。つまり活用(exploitation)重視なので、局所解に陥る可能性があります。

対して、②EI戦略では改善量も考慮に入れられており、imporovement functionの期待値が獲得関数として使われます。ベイズ最適化ではEIが一般的によく使われるそうです。

③UCB戦略は多腕バンティット問題でよく出てくるバンディットアルゴリズムです。 ベイズ最適化を多腕バンディット問題と捉えると、あらかじめ決めた試行（観測-モデル更新）回数で最適解を見つける問題と考えられます。 最適化の試行回数をT、現在の試行回数をtとすると、ハイパーパラメーターν(> 0)を使ってκを置き換えて、上図のGP-UCBのように表すこともできます。

以上見てきた３つの獲得関数（戦略）は、解析的な形をとっており、決定的(deterministic)なアプローチとなっています。これに対してランダムなアプローチもあります。

Thompson Sampling

ランダムなアプローチの例としてThompson Samplingというアルゴリズムがあり、Nature文献中でも他の獲得関数と比較検証されています。

Thomson Samplingについてはバンディットアルゴリズムとしての解説が多いので、多腕バンディット問題に従って内容をおっていきます。

Thompson Samplingはベイズ戦略の一種で、アーム(腕)を「そのアームが最適（期待値が最大）である事後確率」に従ってランダムに選択します。 この事後確率は以下の図中のように表せますが、実際には求める必要はありません。同等となる動作が知られています。*28*29*30

上図中のような手続きを行うことで、「期待値最大の事後確率に基づくアームの選択」が可能になります。 *31 Thompson SamplingはUCB戦略と比べて余分な探索が少なくなることが知られており、有限の試行回数でより良い性能を達成できるそうです。

ベイズ最適化（ガウス過程）にThompson Samplingを利用する際のアルゴリズムを引用しておきます。*32

右側の日本語は本文を適当に訳したものです。

以上、ベイズ最適化の獲得関数を見てきました。最後のThompson Samplingは、他の３つと異なり決定的ではなく、ランダムな手法という違いがありました。 加えて、Thompson Samplingには「次の観測点の候補をサンプリングにもとづいて決める」という手順から、ベイズ最適化の並列化（バッチ処理）と相性が良いと予想されます。

この記事の最後の話題として並列化に触れておきましょう。

ベイズ最適化の並列化（バッチ処理）

今回のNatureの論文の目的は、有機合成反応条件の探索へのベイズ最適化の応用でした。

ベイズ最適化は観測をもとに意思決定を行うサイクルを回す逐次的最適化です。実際の問題への適応を考えると、このサイクルで最も時間がかかるのは実際に実験を行う部分です。

実験作業の効率化を考えると、「1つ実験して、結果を踏まえて条件を変更して次の実験をおこなう」というサイクルよりも、 「同時に並べて複数の実験を行って、それらの結果を踏まえて、また新しく複数の実験を行う」といった、並列化した作業（バッチ処理）の方が好ましくなります。

先に、Thompson Samplingは並列化と相性がいいと書きました。 １サイクルあたりN個の実験を行うバッチ処理に適応したいのであれば、事後確率分布からN個サンプリングを行えば良いと考えれば、なるほどそうかも、と言う感じですね。

では、「決定的なアプローチ（PI、EI、UCB）の場合、どうやってバッチ処理にあわせるか？」というと、Nature文献中ではKriging beliver algorithmと言うのを用いています。

クリギング

クリギング(Kriging)は鉱山技術者ダニー・クリーグにより開発された、採石場全体に含まれる鉱物の総量を推定するための方法です。

サンプリング（ボーリング検査）のコストを減らすため、すでにある座標 x1,..., xN での実測値 fx1,... fxNの重み付き線形和を使って 任意の空間座標 (x)における計測値を推定すると言うものです。現象の生じる位置関係をモデルに含めて説明能力を増すために空間統計学などで使われています。

クリギングとガウス過程回帰法は原理的な方法論としての内容は同じです。カーネル関数でみたMatérnカーネルは空間統計学で広く使われているそうです。*33

クリギング信者アルゴリズム

クリギング信者アルゴリズム(Kriging believer algorithm) は、クリギングに基づく逐次的最適化の問題において、 特に複数の点を同時に最適化(multi-points optimization)する際の戦略です。用語は違いますが、ガウス過程によるベイズ最適化で並列最適化を考えているのと同じですね。

複数の点の同時最適化において、決定的な獲得関数例えばEI）を使っている時、難しくなるポイントは多変量ガウス分布の積分の計算が出てくることです。 この複雑な計算を回避するヒューリスティックな戦略の一つがKriging Believerです。*34

ここでは最適化の各イテレーションで考慮する条件的知識(conditional knowledge)を弱めることを考えます。 具体的には、「一つ前のイテレーションにおける、選んだ箇所での応答(response)に関する条件的知識」を、 「クリギング予測器(krigind predictor)の期待値」で置き換えます。 決定的な値(deterministic value)に置き換えることで、計算を楽にしてやる、という感じですね。*35

Natureの文献中では、「チェスプレーヤーが、複数の手を先読みして、各ステップにおける駒の操作で起こりうる結果を推測する」ようなものだと例えられていました。 Supplementary Information Figure S29. S30.が、Thompson Sampling、Kriging believer algorithmそれぞれの場合のアルゴリズムを記載しているので引用しておきます。*36

以上、ベイズ最適化の並列化(Parallel Optimization)に関する用語でした。

ベイズ最適化まとめ

ベイズ最適化について関連する用語を拾ってきましたが、話題が散らかって細かくなりすぎた感じもあるので、最後にもう一度大枠を見直しておきます。

結局、ベイズ最適化は「高コストなブラックボックス関数の最適点を効率よく求める手法」であり、「関数のサロゲートモデルとしてガウス過程」が出てきました。 ガウス過程の式においては、「カーネル関数（共分散行列）として何を選ぶか？」という選択がありました。 また、最適化のための探索の過程では、効率的な探索のため複数の戦略があり、それぞれ獲得関数を最大化するように次の観測点を選ぶ、という手順が取られました。

最適化のイメージとして、arXiv:1012.2599 (2010). Figure 1がわかりやすいので引用しておきます。

また、解くと言う観点から、ベイズ最適化に関与するハイパーパラメータとしては、ガウス過程の平均関数とカーネル関数のパラメータがあります。これらは最尤推定(MLE)や事後分布最大化(MAP)推定によって求められるそうです。 *37

最後に発展的な話題として、最適化の各イテレーションで観測点を複数とする並列最適化、という点についてもふれました。ベイズ最適化の並列化は難しいらしく、色々と研究がされているみたいです。*38

　まとめ

以上、Nature文献の感想と、読んでもわからなかったベイズ最適化の雰囲気を知るためのお勉強記録でした。

機械学習で新しい反応条件が見つかるというのは個人的にはかなり衝撃でした。学生時代、後追いロースループット実験しかしていなかったショボい身としては、 「そこそこ良い収率出せる配位子あるのに、わざわざ報告例のない配位子にチャレンジするの謎。そこから最適条件までもっていけるとかどう考えても予測不能でしょ。・・・」って感じです。

お仕事なくなっちゃうね！

正直な話、Gaussain（高価そう）でゴリゴリ計算して、ハイスループット実験マシーンをバシバシ動かせる先端機関、数える程しかないのでは？？と思ってしまいました。 ソフトがオープンになっても、入力データを準備できる環境を手に入れるのが大変そうです。

また、本記事は用語まとめといいつつ、書籍やインターネット上の資料の引用記事のようになってしまいました。ごめんなさい。 正しい内容が知りたい方は書籍「ガウス過程と機械学習」を購入しましょう！色々なところでもオススメされているので買って損はないはず！

後半、Thompson SamlingやKriging believer algorithmはどこを参照にすれば良いかわからず、検索して出てきた英語の資料をDeepLに投げ込んだりしているので、かなり内容が怪しくなっています。 オススメの日本語資料、書籍があればご教示いただけると幸いです。

ではでは。

やったね！

*1:購読不要

*2:Supplementary Informationでも比較を確認可能です。

*3:この資料を後悔してくださっている牧山幸史氏はHOXO-M Inc.という企業の方らしいです。O'REILLYの「機械学習のための特徴量エンジニアリング――その原理とPythonによる実践の翻訳をされた企業さんらしい。不勉強で知らなかった。

*4:以下の図の表記は講談社MLP 「ガウス過程と機械学習」の記載も参考にしています

*5:UCB戦略はオープンソースの逆合成解析ソフト AiZynthFinderでも使われた気がします。

*6:こちらの記事（「ベイズ最適化って何がベイズ？」）では、「ベイズ的に獲得関数を更新していく」のでベイズ最適化というのでは、と書かれていました。

*7:グリー株式会社　尾崎嘉彦氏の資料「機械学習モデルのハイパパラメータ最適化」の記載を元にしています。SildeShare

*8:Wikipedia ベイズ確率

*9:Wikipedia 正規分布

*10:Wikipedia 多変量正規分布

*11:「ガウス過程と機械学習 Chapter 3 ガウス過程」 記載より

*12:「~」は確率的生成モデルを表現するのに用いられる記号です。

*13:同上「Chapter 4 確率的生成モデルとガウス過程」 記載より

*14:同上 「Chapter 3 ガウス過程」 記載より

*15:Nature文献中に"Gaussian process represents a distribution over functions"という記載がでてきましたが、この辺りの事を言っているのだと思います。

*16:まあ、私は0章よむのに１週間かかってますけどね。読み終えられる自信がない。。。

*17:「ガウス過程と機械学習 Chapter 3 ガウス過程」記載より

*18:岩波書店「カーネル多変量解析」1,2章ではカーネル多変量解析の概要を難しい数式なしに把握できるよう説明してくださってます。リプレゼンター定理とか出てくる用語が格好いいです。数式の間を言葉でつなぐように書かれているので、式変形を追ううちに目標を見失う私にとってはありがたかったです。同書４章をよんでサポートベクトルマシンの意味するところがやっとわかった気がします。6章,7章は再性核やらヒルベルト空間やら私には全く理解できない話でした。

*19:この流れについてはHELLO CYBERNETICSさんの記事「ガウシアンプロセスを使ってみたい人への最低限の知識」でもわかりやすく説明してくださってます。

*20:また、詳細は省きますが、カーネル関数と等価なガウス過程を事前分布として用いた場合のMAP推定の結果と、カーネル関数の重み付き和のモデルを正則化付きで求めた関数近似の結果は一致します。「カーネル多変量解析 第２章」をご参照ください

*21:Wikipedia カーネル法

*22:明治大学　データ化学工学研究室（金子研究室）HP 「カーネル関数って結局なんなの？→サンプル間の類似度と理解するのがよいと思います！」の記載が参考になります。

*23:「ガウス過程と機械学習 Chapter 3」の式より

*24:参考　Hiromichi Okazaki氏の資料 SlideShare「探索と活用の戦略　ベイズ最適化と多腕バンディット」

*25:参考 牧山幸史氏の資料 SlideShare「機械学習のためのベイズ最適化入門」

*26:明治大学　データ化学工学研究室（金子研究室）HP 「ベイズ最適化(Bayesian Optimization, BO)」

*27:E.Brochu, V.M.Cora, N.de Freitas : A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning, arXiv:1012.2599 (2010).

*28:ALBERT Official Blog バンディットアルゴリズム 基本編

*29:IBIS2014 　本多淳也 氏　資料「多腕バンディット問題の理論とアルゴリズム」

*30:参考 Qiita 「Thompson Samplingで広告配信を最適化してみた」

*31:「アームが最適である確率」に着目した方法を確率一致法といい、Thompson Samplingはベイズ統計の枠組みを確率一致法に適用した方策と言えるそうです。

*32:Kandasamy, K., Krishnamurthy, A., Schneider, J. & Poczos, B. Parallelised Bayesian Optimisation via Thompson Sampling. PMLR 84:133-142, 2018.

*33:「ガウス過程と機械学習 Chapter 6 ガウス過程の適用」

*34:Constant Liar と言う戦略もあるそうです。

*35:参考 Ginsbourger D., Le Riche R., Carraro L. (2010) Kriging Is Well-Suited to Parallelize Optimization. In: Tenne Y., Goh CK. (eds) Computational Intelligence in Expensive Optimization Problems. Adaptation Learning and Optimization, vol 2. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-642-10701-6_6

*36:「Kriging believer algorithm」という用語は出てきませんが、明治大学　金子研究室HPの記事（「ベイズ最適化において一度に複数の実験をするときに候補を選択するシンプルな方法」）で解説してくださっている方法は、ここで引用したアルゴリズムと同じことだと思います。間違っていたらすみません。

*37:創薬・材料探索のための機械学習　「ベイズ最適化(その２)：獲得関数」

*38:Slideshare 「機械学習モデルのハイパパラメータ最適化」

Macで独習！量子化学計算入門！！ 〜インストールでつまずいた話〜

電子書籍「独習　量子化学計算」を購入し、一通りポチポチ遊んでみました。

ゼロからわかる!! 独習 量子化学計算: 理論からはじめない新しい量子化学計算の本

「ゼロからわかる！！」の謳い文句の通り、ど素人でもソフトウェアを動かすことができ、とてもオススメです！

WindowsでもMacでも実践できるように解説してくださっているのですが、 私のMacではソフトウェアのインストールなど、環境構築でいくつかつまずきました。

備忘録としてメモを残します。

尚、本記事は以下の環境での事例です。

• MacBook Pro 2017
• macOS High Sierra version 10.13.6
• プロセッサ名: Intel Core i5
• コア数　２

書籍の紹介とオススメポイント

まず書籍について、初学者からみた感想とオススメポイントをご紹介します。

「ゼロからわかる！！独習　量子化学計算」は、めちゃくちゃわかりやすくて素晴らしいサイト「PC CHEM BSICS.COM」さんが出されている電子書籍です。

こちらのサイト、量子化学の用語を検索すると何度もトップにヒットしました。いずれの記事も簡潔でわかりやすく、初学者の私にとって非常に参考になりました。

また、理論・用語の解説だけでなく、「無料で楽しむ量子化学計算サイト」とのことで、 自宅のPCで試すことができるソフトウェアとその利用方法についてもたくさん紹介してくださっています。

気軽に試せるのは初心者にとってとてもありがたいですよね！

で、この素晴らしいサイトの電子書籍「独習　量子化学計算」ですが、 サブタイトルは「理論からはじめない新しい量子化学計算の本」！！！

「難しいことよくわかんないけど、とりあえずやってみたい！」という、私にとってピッタリすぎる謳い文句です！

早速購入して一通りやってみましたが、「初学者が独学でつまずくポイントがよくわかっていらしゃる・・・」という感じでとても分かりやすかったです。

具体的には、以下のような構成となっているのが良かったです。

1. 実際に計算を経験することを通して学習できるように構成
2. ソフトウェアのインストールから使い方、結果の解析まで解説
3. 一つ一つの章が簡潔で、初学者でも力尽きずに取り掛かかれる
4. ステップバイステップでより複雑な計算の設定がわかるようになる

OSはWindows、Mac両方に対応しており、一方でしか使えないソフトについては、その旨と別の方法の解説がのっています。（使われているソフトはちょっとWindowsより？）

また、操作方法は実際の画面の図で解説されているのでわかりやすくなっています。

個人的にとても助かったのは、操作方法を繰り返し説明してくださっている点です。

私は記憶力が貧弱なので、書籍をよんでいるとよく「基本忘れた！また、前の章に戻って確認しなきゃ。面倒くさいからもういいや。」と投げ出してしまうのですが、今回はそうならずにすみました！*1

ちょっと気になった点は、操作画面の写真の解像度が少し低く、 拡大してもボタンや設定の文字がよくみえないということがありました。*2

あと、「文字検索できたら便利なのになー」と思いました。*3

また、ソフトウェアのインストールや様々なサイトの紹介に際してリンクが記載されていますが、 Kindleから直接Webページに飛んだり、文字列のコピーをしたり、ということができませんでした。ただし、この点は「書籍のサポートページ」にリンク集を用意してくださっているので、そちらで解決しました。

さらにリンクだけでなく、サポートページには「サンプルファイル」も用意してくださっています。「自分で作った初期構造だと計算が収束しない！振動数めっちゃ負！もう嫌だ！！」となった時には、サンプルファイルを利用して計算を進めていくことができます。

「できの悪い人間の諦めポイント」がよくわかっていらっしゃる！ありがとうございます！！

というわけで、「量子化学計算興味はあるけど、どこからスタートすれば良いかわからない」という方は この書籍から始めるのがとてもオススメだと思います！

Macで環境構築する際につまずいた点

では本題、Macで環境構築する際につまずいた点とその回避策の備忘録です。

Avogadroの英語化

まず、Avogadroについてです。 Avogadroはさまざまな量子化学計算ソフトウェアの入力ファイルを作成することができ、 本書でも繰り返し利用方法が解説されています。

MacでAvogadroをダウンロードすると、メニューが日本語化されているのですが、 ちょっと微妙で、英語・日本語の重複があったりして混乱します。

こんな感じ・・・

Avogadroの日本語がちょっと微妙

書籍の表記と合わせるためにも英語化した状態で利用するのがオススメです。

Avogadroのメニューを英語にする方法は、PC CHEM BASICS.COMさんの記事 「Avogadroを使ってみよう」 でも解説されていますが、私の環境ではうまくいきませんでした。

サイト記載の方法は以下の通り。

1. Avogadroのアプリケーションをcontrolボタンを押した状態でクリック
2. 「パッケージ内容を表示」を選択
3. Contents⁩>Resources⁩内の「ja.lproj」をフォルダごと削除

この通りにしても日本語のままでした。

そこで、Avogadroのメールフォーラム(?)（こちら）の記述を参考にして、さらに「xxxxx.qm ファイル」を削除しました。

こちらのファイルもAvogadroを右クリックして「パッケージの中身を表示」すれば確認できます。

このqmファイルが翻訳に関与しているそうなので、削除すれば自動的にデフォルトの英語となります。

こんな感じ・・・

英語にしてしまった方がわかりやすい

メニューの重複も無くなってわかりやすくなりました。

Fireflyのテスト失敗（Wineの設定、prefix）

次にFirefly(PC GAMESS)(version 8.2.0)の計算を実行しようとしてつまずいた点です。*4

Fireflyインストール後、以下のコマンドでテストを実行しましたが失敗しました。

/Applications/Firefly/BENCHMARKS-THREADED-TESTJOBS-RUNSCRIPT-DualCore-Only.sh

ターミナルで実行すると、以下のようなメッセージがでて計算が走りませんでした。

wine: 'ホームディレクトリ/.wine' is a 64-bit prefix, it cannot be used with 32-bit Wine.

「テストファイルが問題？」と、別に用意した「水(H2O)の構造最適化計算」を試しましたが、 こちらもアウトプットがファイルが作成されず、Fireflyがすぐに閉じてしまいました。

エラーメッセージのWINEですが、「Windows用のソフトウェアをMacで利用するための環境をつくる」ためのものだそうです。

Mac用のFireflyをインストールすると、内部に「WINE」というディレクトリもあり、Fireflyパッケージの中に含まれているようです。

で、ソフトウェア起動時に、WINEの設定がホームディレクトリ下の「.wine」ディレクトリにつくられて読み込まれるようなのですが、 どうもこのprefix(?)のバージョンの違いが問題となったようです。

そこでこの「.wine ディレクトリを一度全て削除」してしまったあとに、再度Fireflyのテストを試したところ無事実行することができました！

Firefly実行後に再度確認すると、同じ場所に新しい「.wine」ディレクトリができていました。

注) この方法だとWINEに依存する既存の他のパッケージに影響があるかもしれません。もっと良い方法をご存知の方はご教示いただければ幸いです。

私の場合、以前に別のWindows用ソフト（WinMOPAC）を利用しようとして、Wineを個別にインストールしたことがあり、 その時に作られた設定が今回のFireflyと合わなかったのが原因だったようです。FireflyのWINEはwine 1.1.33 ソースに基づいているようですが、10年くらい前のもののようなので、以前のWINEが64-bitで、Fireflyが32-bitだったのかなぁと推測しています。

尚、Firefly付随のインストールマニュアルに書かれているテスト用のコマンドは、現在のバージョンとずれているので修正が必要です。

マニュアルでは、Dual coreの場合ターミナルで以下を実行するように記載されています。

/Applications/Firefly/BENCHMARKS-THREADED-TESTJOBS-RUNSCRIPT.sh

実際には、配布されているファイルは以下でした。

/Applications/Firefly/BENCHMARKS-THREADED-TESTJOBS-RUNSCRIPT-DualCore-Only.sh

テストされる際にはご注意を！

Fireflyでインストールが推奨されているアプリ

Firefly 8.0.0(for Mac)で推奨されているソフトウェアとして以下のものがあります。

1. TextWrangler
2. wxMacMolPlt

Firefly 7.1.G update 以降、これらのソフトウェアへの依存性(dependency)はなくなっているようなので、 別に入れなくても問題はないようです。

ですがFireflyのインストールマニュアルに「便利だし無料だからマジでおすすめ！（“strongly recommended”）」と書かれていたので導入しました。 デフォルトの「/Applications」ディレクトリに入れておくことがオススメされています。

TextWranglerは無料のテキストエディタです。 Firefly計算実行した際に出力の「.outファイル」がTextWranglerで自動的に開かれ、 逐次内容が更新されるので計算が走っているかどうかの確認に便利です。・・・が、邪魔といえば邪魔です。。。

無料のテキストエディタは色々あるので、別にTextWranglerに拘らなくても良いのかなーという感じです。 オススメな理由をご存知な方がいらっしゃったら教えていただければ幸いです。

もう一つのwxMacMolPltは、計算結果の解析に使うことができ、「独習　量子化学計算」でも何度もでてきます。こっちはインストールしておきましょう！

GamessQ立ち上げ時のエラー

2021年現在、「GAMESS for MacOS X」にはGamessQが含まれていません。別途インストールする必要があります。

私のMacOS環境では「GamessQ (v1.2.2)」を使おうとすると以下のエラーが出ました。バックエンドの立ち上げに失敗しているみたいです。

GamessQのエラー

「Add」などのボタンがグレーになっていて使うことができません。・・・困った。

回避策（？）はこちらのGitHubページに書かれていました。（→Error launching backend on fresh install under OS X #19）

1. GamessQのディレクトリにあるコマンドラインバイナリ「gamessqd」をダブルクリックして実行
2. アプリ「GamessQ.app」を立ち上げ

黒い四角のアイコン gamessqdをクリックするとターミナルが立ち上がります。 この状態でGamessQのアプリを立ち上げれば先のエラーが出ず、正常(?)に起動できました！

「GamessQ.app」の前に「gamessqdd」を立ち上げる

「Add」や「Clean up」といったボタンが使えるようになっています（アクティブ？）。

アプリを閉じるとターミナルにもlogoutと表示されてプロセスが終了します。

Gamess計算実行時のエラー（HOST_NOT_FOUND）

GamessQが使えるようにできたので、計算を実行したところ結果が得られずログに以下のエラーが出ました。

 Error: Gethostbyname(XYZXYZ) returned HOST_NOT_FOUND.
 ddikick.x: Fatal error detected.

Gethostbyname(XYZXYZ)のXYZXYZの箇所には、実際には私のMacBookに設定した名前が入っています。*5

上記のエラーはGamessQを使わずに、Gamess単独で以下のテストを実行しても再現しました。（GAMESS(US)は少し古い2018年のビルドです）*6

1. ターミナルでGamess実行ファイルを格納したディレクトリに移動
2. ./gms tests/standard/exam01
3. output file name? [tests/standard/exam01.log]と聞かれるのでReturnキーを押す
4. 出力されたexam01.logを確認

Gethostbyname()は指定したホスト名から、ホストに関する情報を得るための命令のようです。
つまり「PC(MacBook）に設定した名前のホストが見つからないよ!!（HOST_NOT_FOUND）」みたいなエラー？？？

色々調べた結果、macOSのhostsファイルの設定を追加してやればよいのではないか、となりました。

hostsは「./etc」（あるいは「./private/etc」）にあります。テキストエディタで開くとこんな感じでした。

##
# Host Database
#
# localhost is used to configure the loopback interface
# when the system is booting.  Do not change this entry.
##
127.0.0.1 localhost
255.255.255.255 broadcasthost
::1             localhost

「よくわからないけどローカルだからローカルホストでしょ！」ということで、 7行目の後ろにPCの名前（XYZXYZ）を追加して以下のように書き換えました。（書き換えには権限が必要です）

127.0.0.1 localhost XYZXYZ(←MacBookにつけた名前)

hostsファイル変更後に、Gamessのテストを再度実行すると以下のようなポップアップが出てきました。

ddikick? OK!OK!・・・知らんけど。

「ネットワーク受信接続」を許可すると「exam01.dat」ファイルが作成され、無事計算実行できました！

「exam01.log」の中身も 「EXECUTION OF GAMESS TERMINATED NORMALLY」となっていました。

こちらのページによると、 「127.0.0.1はローカル・ループバック・アドレスと呼ばれ、自分自身を指す特別なIPアドレスである。「localhost」という名前でも参照できる。自分自身の上で動作しているサービスへ接続する場合は、このIPアドレスを利用できる。」だそうです。

GamessはローカルPCで計算を実行する場合でも、内部でネットワーク接続処理（TCP/IP?)を行う構成になっているのでしょうか？？？

localhostに追加するのが正しいのか全くわかりませんがとりあえず動いたのでヨシ！！

IPアドレスをもう少し検証

せっかくなのでもう少し実験してみました。

先ほどのページには以下の記載がありました。*7

一般的には「127.0.0.1」というIPアドレス（IPv4の場合）が利用されるが、実際にはIPアドレスの最上位のバイト（最上位の8bit）の内容が「127」でありさえすればよいので、「127.0.0.1～127.255.255.254」の範囲内ならばどのIPアドレスでも利用できる（127.0.0.0と127.255.255.255の2つはブロードキャスト・アドレスのため除外される）。例えば127.0.0.2でもよいし、127.1.2.3でもよいが、一般的な用途では127.0.0.1だけが利用される。

別のアドレスでもOKなら、localhostの後ろに加えずに「127.0.0.2 XYZXYZ」としてもよいのでは？？？

ということでhostsファイルを以下のようにして再テストしてみました。

127.0.0.1 localhostt
127.0.0.2 XYZXYZ(←MacBookにつけた名前)

127.0.0.2を使っても計算は終了し、「exam01.log」は EXECUTION OF GAMESS TERMINATED NORMALLYとなりました。

ですが、データファイル「exam01.dat」（PUNCH出力）が作成されず、代わりに「exam01.F05」というファイルができました。 また、ターミナルもコマンド実行中のままで次に移行しませんでした。

logには計算結果が書かれていたものの、「.datファイル」（PUNCHファイル）への出力、計算プロセスの終了がうまくいっていない感じでしょうか？

正常な動作のためにはlocalhostと同じIPアドレスにしておいた方が良さそうです。

ソフトウェアの通信やコマンドラインの処理、何が行われているのかさっぱり意味がわからないです。。。。

wxMacMolPltによる結果の可視化の違い（FireflyとGamessQ)

「独習　量子化学計算　Chapter 8 フロンティア軌道で反応を予測しよう」では、MacMolPltにより分子軌道を可視化する方法が解説されています。

この操作で、書籍では「Fireflyで得た計算結果（.outファイル）」を解析していましたが、私の環境では再現できませんでした。

一方で、「GamessQで同じinputファイルを使って得た計算結果（.logファイル）」では、可視化操作がうまくいきました。

行いたかった作業は以下のように「Select Orbital Set : Molecular Eigenvectors」で、 「Select Orb : Energy」から目的の分子軌道を選択し描画することです。

やりたかった作業

上図は、GamessQで得たlogファイルを用いたものです。対して、Fireflyのoutファイルでは以下となりました。

Fireflyの出力ファイルだとできない

上図のように、そもそも選択できる描画「Surfce Type」が少なく、「Atomic Orbitals」しか選択できませんでした。

Fireflyの計算結果（.outファイル）の中身を確認すると、以下のように「MOLECULAR ORBOTALS」の情報が記載されていたので、 おそらく描画ソフト側（wxMacMolPlt）の問題ではないかと推測しています。

Fireflyでも分子軌道は出力されている

私のPC環境（or MacOS）特有の問題かもしれませんが、同じ問題に出くわした方のためにご参考までに紹介しておきます。

おまけ Gatekeeperへのアプリ許可の追加

GamessQのbackend エラー回避で参照したGitHubにGateKeeperのパーミッションについても言及されていました。 「administratorの権限がないと、先の回避策が使えないのでは？？？」って感じの話みたいです。

こちらのページによると、「Gatekeeperとは、Mac App Storeからインストールするかユーザーの右クリック操作により許可されたアプリだけを起動できるようにする、セキュリティ強化の仕組み」だそうです。

「開発元本当に信頼できる？？？」って聞いてくるあれですね！

この許可をターミナルのコマンド「spctl -add　パス」を使って処理することもできるとのことです。

例えば、GamessQの場合なら以下のようなコマンドとなります。（パスは適宜修正してください）

spctl --add "/Applications/GamessQ/GamessQ.app"
spctl --add "/Applications/GamessQ/gamessqd"

許可リストから外す場合は「spctl -remove　パス」とすれば良いそうです。

spctl --remove "/Applications/GamessQ/GamessQ.app"
spctl --remove "/Applications/GamessQ/gamessqd"

インストールの際にsudoと組み合わせておけば、administrator権限がない場合でも問題を回避できるようになるみたいです。（GitHubでは実習で生徒が使うPCの設定が議論されてました）

まとめ

以上、「独習　量子化学計算」の感想と、Macで量子化学計算の環境構築につまずいた点のご紹介でした！

繰り返しになりますが、とりあえず計算を始めてみたい初学者にとてもオススメな書籍という感想です。*8

ハートリー・フォックやDFTという用語をなんとなく知ったものの、実際にはどんな感じか全くわからなかったので非常に勉強になりました。

動かしてみると、「SCF計算ってわかったつもりになってたけど、具体的にどの部分の何を計算してるんだっけ？？？」とか、 「DFTで計算したら、HFでできたアセトンの構造最適化失敗するやん！負の振動とは一体？？？」などなど、新たな疑問がたくさん出てきました。

すこし初期構造を変えるだけで計算がうまくいかなくなったりして、計算化学の実験科学的な側面をちょっと体感できたのでよかったです。たぶん論文に載っているような結果は、膨大な試行錯誤の結果なのでしょうね。格好いい図くらいに思っててすみませんでした。

あと、途中で北浦-諸熊解析という用語が出てきて、おおっよくわからないけど格好いい！！！これが噂のONIOM？？FMO？？とテンション上がりましたが全然関係なかったです。

水分子１つの計算でつまずいているのに、タンパク質の計算を理解できる日はくるのでしょうか。。。。

ゆれるアセトン

*1:レベル低い感想ですみません

*2:外付けモニタと私の視力のせいかも・・・

*3:私のKindleアプリが古いせいかも

*4:Fireflyのライセンス認証のメールが中々来ないと思っていたら、Gmailの迷惑フォルダに入っていました。受信日、翌日だったのでレンスポンスめちゃくちゃ早くてびっくりです。

*5:恥ずかしいので隠した

*6:テスト方法参考　量子化学計算ソフト　GAMESS のエラー

*7:ローカル・ループバック・アドレス（127.0.0.1）とは？

*8:回し者ではないですよ！ノーアフィリエイトですよ！

雰囲気量子化学入門（後編） ~ハートリー・フォック法の課題と密度汎関数理論〜

雰囲気量子化学入門・後編！前回に引き続き、明日知ったかぶれる量子化学を目指して用語をひろっていきます！

前編ではシュレーディンガー方程式からハートリー・フォック方程式までを辿りました。 後編ではハートリー・フォック法の問題点とその解決策、そしてDFTについてお勉強！

全体の流れ

まず、全体の流れを再掲します。だいたいこんな感じの話です。

今回は下半分、ハートリー・フォック方程式の誤差から出発したいと思います。

ハートリー・フォック法の誤差（電子相関）

さて、シュレーディンガー方程式から近似で導いたハートリー・フォック方程式で得られるエネルギーは、全体のエネルギーの99%以上を見積もることができます *1。残る誤差は1%程度ですが、この大きさは化学に適応する点で非常に問題となります。

なぜかというと、化学で考慮したい課題は、化学反応の活性化エネルギーやイオン化エネルギーといった問題ですが、これらは全エネルギーの絶対値ではなく、 そのエネルギー差（相対値）によって結果が左右されます。このエネルギー差の大きさが、ハートリー・フォック法の誤差と同程度となるため、 化学的精度にとって無視できない誤差となってしまいます。*2

この誤差は電子相関と呼ばれ、レフディンにより電子相関エネルギーは以下のように定義されています。

大まかにいうと、電子相関は「多電子系において1個の電子の位置を定めると、他の1個の電子の存在確率が第１の電子の位置の影響を受ける」 というようなことを意味するようです。*3

ではなぜハートリー・フォック法ではこの相関を一部しか取り込めず、誤差となってしまうのでしょうか？要因と解決策を辿ってみましょう。

電子相関の起源とポスト・ハートリー・フォック法

ここではより具体的な電子相関の起源と、その解決策として考案されたハートリー・フォック法の延長線上にある手法を辿ります。

ざっとこんな感じです。

動的電子相関と静的電子相関

まず、電子相関は大きく動的電子相関と静的電子相関の２種類に分類されます。*4

動的電子相関はハートリー・フォック法が独立粒子モデルに基づいていることに由来します。 このモデルでは「1つの電子に着目した平均場近似」となっているため、電子間の位置の相互作用を平均的なものに置き換えています。 つまり、2つ以上の電子の衝突や散乱といった効果が含まれません。この、「平均化により表現できなくなった部分」が動的電子相関の起源となっています。

格好いい表現でいうと、「ハートリー・フォック波動関数は相関カスプ条件を満足しないため、動的電子相関が欠如しており、 短距離電子反発を余計に取り込んでしまう」（？？？）ということになるようです。

要するに、実際には電子同士は互いに反発するのでゼロ距離まで近づく確率はとても低いけど、ハートリー・フォックはそんな状況が存在する可能性も含めてしまう*5。余分な反発まで考えてしまうので、結合エネルギーの過小評価や、結合長の過大評価につながったりしてしまう、ってことでしょうか？？

一方、静的電子相関は１つのスレーター行列式（電子配置）を使って波動関数を近似していることに由来します。

格好いい表現では「いくつかの電子状態（電子配置）がエネルギー的に近く、擬縮退している場合にみられる 配置間の長距離相互作用（擬縮退効果）」（？？？）となります。

良くわからないのでもう少し具体的に、例として水素分子を水素原子に解離することを考えます。

この時、閉殻ハートリー・フォック（RHF）法では「１つの空間軌道に2つの電子をペアで詰める」ので、2つの水素原子の分子軌道のうち、一方に電子2個、もう一方に電子が0個のイオン化した状態しか得られません。ですが、実際にはそれぞれ1個ずつ電子を持つ水素原子として解離した状態の方がエネルギー的に安定です。従って"不自然"な状態を過大評価してしまうことになっています。

つまり、「１つよりも同時に２つを考えた方がいい結果が得られる時があるのに、それが上手くできないので問題だ！」って感じでしょうか？

・・・良くわかりませんがとにかく誤差があるんだ！！！どげんかせんといかん！

ポスト・ハートリー・フォック法

ハートリー・フォック法で取りこぼしてしまう電子相関を考慮するために、様々な手法が考案されています。

上図に示したように大きく次の3つの分類が可能です。

• ① 動的電子相関を考慮するためのポスト・ハートリー・フォック法
• ② 静的電子相関を考慮するための多配置 (multi-configuration) SCF法
• ③ どちらも考慮する多参照 (multi-reference)電子相関法

このうち、ポスト・ハートリー・フォック法では、仮想軌道を用いることで励起状態の電子配置(励起配置)を考慮します。 原理的には「全ての励起配置を考慮すれば、用いている基底関数の空間内で、電子相関を厳密に取り扱える」ことになります。 一方で、励起配置の数が増えるほど計算コストが増大するという課題があります。

ポスト・ハートリー・フォック法には大きく① 配置間相互作用(CI)法、② 摂動法、③ クラスター展開(CC)法があります。 それぞれこんな感じです。

CI法は線型結合で表現されており変分法による近似でCI係数を決定します。

一方、摂動法はその名の通り近似に摂動法をもちいています。 ハミルトニアンを無摂動ハミルトニアンと摂動項に分割し、できるだけ少ない逐次近似で真の解に近づくような分割の仕方が目指されます。

これらに対して、CC法では演算子が指数関数となっている点に特徴があります。CI法と同じ数の変数でより多くの励起配置を考慮でき、 より高い精度の結果が得られると言う利点があります。

いずれも「電子の励起配置をどこまで考慮するか？」によって呼び方が変わります。

MP2とかCCSD(T)とかなんとなく聞いたことありましたけどポスト・ハートリー・フォック法だったんですねー。

多配置SCF法、多参照電子相関法

ポスト・ハートリー・フォック法以外も少し触れておきます。

まず、静的電子相関を考慮する方法として、多配置SCF法があります。 基本的考えは「複数のスレーター行列を使って波動関数を表現する」というもので、配置関数によるCI展開を用います。 配置関数の選び方に複数あり、代表的なものはCASSCF法です。*6

また、両方（動的＋静的）の電子相関を考慮する方法として多参照電子相関法があります。 多配置SCF法で静的な部分を考慮しつつ、ポスト・ハートリー・フォック法のいずれかを用いて動的な部分を取り込む、組合せのアプローチです。*7

以上、ハートリー・フォック法の拡張を見てきました。これらはいずれも「波動関数に基づいて電子相関を考慮するアプローチ」ですが、計算コストが大きくなるという問題点があります。

これに対して、より低い計算コストで電子相関を取り込むことができる手法があります。それが密度汎関数理論(DEnsity Functional Theory, DFT)です。

ようやく本題！

密度汎関数理論

DFTでは、波動関数を求めるのではなく、代わりに「電子密度から出発して電子相関を考慮」します。 この「波動関数 → 電子密度」というのがDFTの最大の特徴です。

これまで、解けない「波動関数に関するシュレーディンガー方程式」を解くため、近似を頑張ってハートリー・フォック方程式へと導きました。 対して、DFTではシュレーディンガー方程式と等価な「電子密度に関するコーン・シャム方程式」を用います。

どうしてこれがOKになるのか、順番に見ていきましょう！

お話の概要を再掲しときます。

ホーヘンベルグ・コーンの定理とエネルギー汎関数

まず、ホーヘンベルグ・コーン（Hohenberg-Kohn）の定理と呼ばれる基本定理から出発します。

以下の2つの定理からなります。

① 存在定理は要するに、「系のすべての物理量は基底状態の電子密度を与えれば一意的に決まる、すなわち電子密度の汎関数である」ということを意味するそうです。

上図に記載の通り、エネルギーも電子密度の汎関数（エネルギー汎関数）となっています。右辺はそれぞれ、T[ρ]：運動エネルギー、V_ext[ρ]：外部ポテンシャルエネルギー、V_ee[ρ]：電子間相互作用エネルギーの項を表します。

また、②　変分原理は、「種々の電子密度 (ρ)のうち、エネルギー汎関数(E[ρ])に代入して得られるエネルギーが最小となるような電子密度（ρ_0）が、 系の基底状態の電子密度であり、「正解」である」といった意味になるそうです。つまり、「電子密度(ρ)を色々と変化させて、とにかく条件に合うものを見つけてしまえば解けたことになる」、と！ *8

たぶん凄いことなんですが、私の能力ではイマイチ凄さがつかみきれない！悔しい！

要するに、「等価だから、シュレ-ディンガー方程式を解く代わりにDFTのアプローチを採用した方が計算楽かもよ！」ってことですね。知らんけど

エネルギー汎関数の形

さて、等価な問題と書きましたが、正確な電子状態が求まるのは、「エネルギーと電子密度の間の対応関係を正確に表す汎関数がわかっている場合」に限られます。

しかし、残念ながらそのような汎関数の形はわかっておらず、結局DFTにおいても近似に基づいた汎関数を使うことになります。

近似として以下のようなものが提案されましたが、精度が低く*9、実用性に問題がありました。

これらに対し、コーンとシャムが導入した近似により、DFTは実用上有用な方法となりました。

コーン・シャム近似

コーンとシャムが導入した近似はだいたい以下のようなものです。

厳密な表式がわからず近似の精度の悪さに関係していた運動エネルギーに着目し、正確に表せる部分とそれ以外の残りの部分に分けました。 前者はコーン・シャム軌道と呼ばれる軌道を導入して表し、残りの表しきれない部分を交換相関汎関数に押し込めてしまいました。

また、加えて、古典的な電子間クーロンエネルギーで近似できない量子力学的な電子間相互作用についても、この交換相関汎関数に押し込められました。

「運動エネルギーの”正確に表せる部分”」にもう少し踏み込みます。ここでは、現実の複雑な電子同士の相互作用のある系の代わりに、「相互作用のない仮想的な系（Kohn-Shamの補助系）」を導入します。この補助系は「現実の系と同じ電子密度（+同じ全エネルギー）」をもつように設定されており、 補助系の運動エネルギー汎関数は厳密に表現することができます。 *10

こうして得られたコーン・シャム法のエネルギー汎関数では、独立粒子運動エネルギー項(T_KS)と 原子核-電子間ポテンシャルエネルギー汎関数(V_ne)、古典的な電子-電子間のクーロン相互作用によるエネルギー項(J)、 そして交換相関汎関数(E_XC)とに分割されて表されるようになりました。

交換相関汎関数の形

さて、コーン・シャム近似でよくわからない部分を押し込められた交換相関汎関数ですが、やはり厳密な形がわからず、 様々な近似で表すことがアプローチされています。

こんな感じ。

いよいよB3LYP出てきました！！

さて、交換相関汎関数の最も単純な近似は一様電子ガスモデルから導出された「局所密度近似(local density approximation、LDA)」でした。 この精度を高めるために、電子密度の勾配(∇ρ)の効果（電子密度の変化に対するエネルギーの依存性）を考慮してLDAを補正した、「一般化された密度勾配近似(generalized gradient approximation, GGA)」が考案されました。

しかし、依然として「交換エネルギーをポテンシャル関数で表す近似」では化学的精度を得るには粗い場合があり、密度汎関数法とハートリーフォック法とのハイブリッドと呼ばれる方法が提案されました。この「混成GGA汎関数(hybrid-GGA)」では、ハートリー・フォック交換項を使ってGGAに補正が加えられています。*11

交換相関汎関数の具体的な形として多くが提案されていますが、通例「交換汎関数+相関汎関数」と両者を組み合わせた名称で呼ばれます。 これは、「電子のスピンに由来する効果を記述する交換汎関数」と「電子相関に由来する効果を記述する相関汎関数」とにわけて それぞれ開発されてきた経緯によります。

で、Beckeによる交換汎関数(B)とLee-Yang-Parrによる相関汎関数(LYP)の組み合わせに、ハートリー・フォック交換項を加えた混成GG汎関数が B3LYPであったと！なるほど！！

コーン・シャム方程式とその解き方

さて、交換相関汎関数についてもわかりコーン・シャム法のあらましがわかってきました。あとはこの問題をどう解けば良いか？、です。

が、ちょっと用語が多すぎて目標を見失なってしまったので、もう一度振り返ります。

とにかく我々の目標は「系（分子）の状態（物理量）」を知ることです。シュレーディンガー方程式では波動関数を求めればOKでした。 対して、DFTでは基底状態の電子密度を求めます。で、この電子状態を求めるのに必要な（だけど正確な形がわからない）汎関数の近似のために、 コーン・シャム軌道や交換相関汎関数が導入されました。ここで、コーン・シャム軌道は相互作用のない独立電子系（コーン・シャムの補助系）を考えているので。 １つの電子に関する軌道です。

１電子軌道とエネルギーの問題・・・前回の記事でみたハートリー・フォック法と同じ形式の問題ですね！

ハートリー・フォック法における解法のアプローチはこんな感じでした。

コーン・シャム法においても同じアプローチを用いることができます。

変分原理よりエネルギーを最小化するものが最も良い近似なので、コーン・シャム軌道の規格直交条件を付加条件としてラグランジュの未定乗数法を使うことで、コーン・シャム方程式と呼ばれる方程式が導出されます。

こんな感じ。

導いたコーン・シャム方程式は、ハートリー・フォック方程式と同じ形の1電子方程式となっているので、同様に解くことができます。

ハートリー・フォック・ローターン方程式を導いた時のように、基底関数系を用いた展開によって行列問題に変換（コーン・シャム行列）すれば、 SCFの手続きによる繰り返し計算によって解くことができます。

前回の記事で「行列問題への置き換えはDFT計算においても主流の解き方」と書きましたが、確かにコーン・シャム法も同じアプローチで解けましたね！

また、この置き換えにおいて、DFTでも基底関数系が用いられました。前回見たように基底関数系の例として有名なものにPople型基底関数(ex. 6-31G)があるのでした。

これで論文で良く見る表現「B3LYP/6-31G」が「交換相関汎関数/基底関数系」の組み合わせを表すことがわかりましたね！！これで今日から知ったかぶれる！！

コーン・シャム法の振り返りと解の意味するもの

無事、密度汎関数理論において解くべき方程式（コーン・シャム方程式）とその解法まで辿ることができました。最後に、以上の手続きで得られたものについて振り返りたいと思います。

まず、密度汎関数理論の導入としてハートリー・フォック法におけるエネルギーの誤差の問題（電子相関)から出発しました。 ハートリー・フォック法の延長にあるポスト・ハートリー・フォック法では、電子相関を「基底状態と励起状態との電子配置間の相互作用」として取り込みました。

対して上図で示したDFTのコーン・シャム方程式において、交換相関ポテンシャルは「交換相関エネルギー(E_XC)を電子密度（ρ）で微分したもの」として定義されています。つまり、コーン・シャム法においては、電子相関はポテンシャル汎関数として取り込まれていることになります。

ハートリー・フォック法とコーン・シャムDFTの式を見比べてみます。

ほとんど同じ形ですが、コーン・シャムDFTではハートリー・フォック法のフォック演算子の交換演算子(K)が、交換相関ポテンシャル(ν_XC)に変わっています。交換相関ポテンシャルに、ハートリー・フォック交換項に加えて電子相関の効果も含まれているため、DFTではより精度良く系の記述ができるようになっているそうです。

さらに、DFTはポスト・ハートリー・フォック法と比べて計算コストが低く抑えられています。ポストHF法では複数の電子配置を考えるので、精度を高めようと考慮する励起配置の数を増やしていくと、計算コストがどんどんと高くなっていきます。対してコーン・シャム法では「基底状態の電子配置のみを考え、かつ１電子方程式を解くのみで電子相関を取り込める」ので、計算が楽！ということのようです。

また、コーン・シャム近似では運動エネルギー汎関数の表現のためにコーン・シャム軌道を導入しました。 注意すべきなのは、この軌道は問題をとくために仮想的に導入された人工の系（コーン・シャムの補助系）における電子軌道であり、 「直接には物理的な意味をもっていない」という点です。従ってコーン・シャム方程式のエネルギー固有値も実際の系のものと見做すことはできないそうです。

ただし実際には、コーン・シャム方程式のエネルギー固有値は実験値と良く比較され、その根拠はJanakの定理と呼ばれるものにあるようです。より詳しくはこちらの記事をご参照ください。

「ハートリー・フォック方程式とコーン・シャム方程式似てるし、解き方も一緒だから、結果の意味も同じでしょ！」と思ってしまいましたが、そんな簡単な話ではなかったんですね！詳しくは理解できませんでしたがとりあえずコーン・シャムに基づいてエネルギーの議論に持ち込んでもOKそうなので安心（？？）。

まとめ

以上、雰囲気量子化学入門（後編）でした。後編では、ハートリー・フォック行列の誤差(電子相関)から出発し、 その解決策として、拡張であるポスト・ハートリー・フォック法をみたあと、本命密度汎関数理論(DFT)の概要を辿りました。

イマイチ、出発点の動的電子相関と静的電子相関の具体的な意味を掴みきれていない気もしますが。。。

ぼんやりと聞いていたDFTが、波動関数ではなく電子密度に基づいたものであり、シュレーディンガー方程式ではなく、 等価な別の問題を解くものであったというのは全く理解していなかったので、びっくりしました。

別問題とはいいつつも、式の形や解くための変換の考え方はハートリー・フォックと共通する部分が多かったので、 前回の記事と合わせて全体の流れを辿ること何となく話を追うことができました。

中島隆人先生の「量子化学　ー分子軌道法の理解のためにー」、初学者にとって分かりやすい構成になっていてオススメですよー！

相変わらず、具体的な数式部分はさっぱり分かりませんでしたが、とりあえずB3LYP/6-31Gの意味がなんとなくわかったのでヨシ！！

今回もひどい誤解がたくさんあると思うのでご指摘いただければ幸いです。

参考資料

この記事は前編同様主に「量子化学　ー分子軌道法の理解のためにー 中島隆人　著」（裳華房）の内容を参考にしています。

上記書籍はハートリー・フォック法の問題点の指摘とその解決策の流れで密度汎関数理論が導入されており、流れを追いやすい一方で、 密度汎関数理論そのものについての記載量はあまり多くありません。 そこで、この記事の密度汎関数理論における記載はインターネット上に公開してくださっている以下の資料を参考にさせていただいています。 ただし、非専門家の私の独学での把握に基づいて記載しているので、誤りの責は全て私にあります。

密度汎関数理論の大まかな概要と用語についてはChem-Stationのゼロさんの一連の記事がコンパクトで非常にわかりやすかったです。

• Chem-Station
  • 計算化学：汎関数って何？
  • 計算化学：基底関数って何？
  • 計算化学：DFTって何？ part II
  • 計算化学：DFTって何？ part III

より踏み込んだ内容、数式の展開、コーン・シャム法の意味するところなどについては@dc1394さんの一連の解説記事が非常に参考になりました。

• @dc1394氏による解説
  • できるだけ簡単にHatree-Fock法を解説してみる（C++17のコード付き）
  • できるだけ簡単にHatree-Fock法を解説してみる part2
  • できるだけ簡単に密度汎関数理論（Density Functional Theory, DFT）を説明してみる（前編）
  • できるだけ簡単に密度汎関数理論（Density Functional Theory, DFT）を説明してみる（後編）

また、汎関数についての具体的な例と説明、B3LYP/6-31Gについては以下の記事が参考になりました。

• 計算科学.com
  • 汎関数の選び方「DFT計算」
  • 汎関数データベース
• PC CHEM BASICS.COM
  • B3LYP/6-31G(d)とは何なのか？

ハートリー・フォック法の誤差（電子相関）に関する記述は以下を参考にしています。動的電子相関・静的電子相関の具体的な意味合いについてはYouTubeの湯どうふさんの動画が大変参考になりました。

• 理化学研究所　常田グループ（理論化学）　講義資料
• YouTubeチャンネル 湯どうふ（化学の勉強）
  • 【量子化学講義 中級 #6】電子相関、独立電子近似を超えて【最終回】
  • 波動関数の満たすべき条件：カスプ条件【量子化学、理論化学
• wikipedia 電子相関

常田貴夫先生のこちらの本、Chem-Stationでも計算科学.comでもオススメされているのでとても気になるのですが、 懐が寂しいので買うべきか迷っています。。。

密度汎関数法の基礎 (ＫＳ物理専門書)

※ 書き終わってから気づきましたが、大阪大学　白井光雲先生の以下の資料は「原理の主張と意味」に焦点をあてて解説してくださっていて、とてもわかりやすいかったです。もっと早く知りたかった。

*「第一原理と密度汎関数理論」 11 March 2005

おまけ（メモ）

• 半経験的分子軌道法

前回、ハートーリー・フォック法では経験的なパラメーターを用いないので、HF法に基づく分子軌道法をab initio分子軌道法や非経験的分子軌道法と呼ぶとかきました。 これに対して、煩雑な計算をすべて（１から）解くのではなく、一部、実験で得られたパラメーターをつかって計算する方法が半経験的分子軌道法です。 計算コストを下げる以外にも、電子相関の効果をパラメーターを通して取り込むといった点が期待できるそうです。 手法の例としてはAM1(Austin Model 1)法や、PM3(Parametic Methot 3)法といったものがあります。

• v-表示可能性問題

DFTの基本定理「ホーヘンベルク・コーンの定理」の第１定理（存在定理）は「v-表示可能性の仮定」を用いています。これは「任意の電子密度 ρを与えると、必ず、それを基底状態の電子密度とするような外部ポテンシャル v が１つは存在する」というものです。 *12 この表現を使うと、存在定理は「基底状態の波動関数と、v-表示可能である電子密度とは１対１の対応がある」と言い換えることもできるそうです。*13

問題は、「v-表示可能性であることの必要十分条件がわかっていない」ことで、このv-表示可能性問題の解決として「レヴィ(Levy)の制限付き探索」があります。

• レヴィの制限付き探索

レヴィらは「v-表示可能な範囲」よりも広い「N-表示可能」な範囲に拡張して以下のように書き直しました。

この密度汎関数Q[ρ] （Levy-Liebの密度汎関数）を使うとv-表示可能な領域の外でもホーヘンベルク・コーンの定理の第２定理（密度変分原理）が成り立ちます。これに従って変分探索を行い、エネルギーを最小化するような電子密度（ρ）を探せばv-表示可能性問題を回避できることになります。

もう少し具体的には、２段階の最小化工程となっています。

1. まずρを固定して、そのρを与える波動関数Ψ_ρの中で「運動エネルギー（^T）と電子間相互作用（^V_ee)の期待値」を評価し、期待値を最小化するようなΨ_ρを探します。このΨ_ρをQ[ρ]と定義します。(１つ目の式)
2. 次に、今度はρを固定せずに、2つ目の式のエネルギーを最小化するようなρを探します。

このプロセスが「レヴィの制限付き探索」と呼ばれるものです。

尚、N-表示可能は「N電子系において基底状態に対応する反対称化した波動関数(Ψ)を与える電子密度(ρ)」といったような意味で、以下のような簡単な条件をみたすどのようなρに対してもそのようなΨが存在することがわかっているそうです（電子密度のN表現可能性）。

以上、よくわからなかったので本文に含めなかった事項の追加メモでした。

あと、いらすとやさんなんでもある。

*1: ハートリー・フォック極限のエネルギーで99.5%程度見積れるそうです。HF極限は基底関数を無限に増やしたHFR 計算によるエネルギーに対応するもののようです。

*2:以下の図で矢印の向きがごちゃごちゃしてしまってすみません。左半分について、「正確なエネルギーよりもハートリー・フォックエネルギーの方がエネルギーの値が高く、精度が上がるに連れてエネルギー値が低くなっていく」ように並べています。HF方程式を導くときに変分法を使ったので、「変分原理 ：近似のエネルギー ε ≧ 真のエネルギー E 」を意識すれば、このような並び順となるはず。たぶん・・・（参考 wikipedia - Electronic correlation ）

*3:参考 Wikipedia - 電子相関

*4:電子相関の説明と図の作成に当たってはWebに公開されている理化学研究所 常田貴夫先生の講義資料とYouTubeの湯どうふ（化学の勉強）さんの説明動画を参考にさせていただいています

*5:波動関数は電子の存在確率に関係する（コペンハーゲン解釈）ということを踏まえている？？

*6:価電子軌道（HOMO,LUMO）付近の分子軌道を用いて、その組の中で電子の励起で生じる電子配置を考慮（active space）します。CASSCFではこのactive spaceの全ての電子配置を考慮し、配置関数とする方法です。

*7:尚、多参照に対して単参照電子相関法はHF波動関数から出発して動的電子相関を考慮する相関法のことです。

*8:参考　できるだけ簡単に密度汎関数理論（Density Functional Theory, DFT）を説明してみる（前編）

*9:原子の全エネルギーを15%~50%程度低くみつもってしまうそう

*10: コーン・シャム軌道の項の表式（T_KS) 参照

*11:参考　「大学院講義　有機化学 Ⅰ.分子構造と反応・有機金属化学」　東京化学同人　第１版

*12:外部ポテンシャルは、物質中の電子に対しては原子核によるCoulombポテンシャルを意味するそうです。

*13:前記事のコメント欄でご紹介いただいた資料をご参照ください

雰囲気量子化学入門（前編） ~シュレーディンガー方程式からハートリー・フォック法まで〜

量子化学ってなんだか格好良くて憧れてしまいますよね！で、学生の頃疑問だったのが講義と実践の圧倒的解離。。。

講義ではいつも「シュレーディンガー方程式入門！」「水素原子解いちゃうよ！」で終わってしまうのに、学会や論文では、「ここはDFTでー、B3LYPでー」みたいな謎用語が繰り出される。。。、 「え！何それ？？何この飛躍？？？」となっていました。

で、数式わからないけど知ったかぶりたい！格好つけたい！というわけでそれっぽい用語（？）をひろってみました。

参考文献はこちら！本棚の奥から出てきた本です。

量子化学 -分子軌道法の理解のために (化学の指針シリーズ)

では早速、雰囲気量子化学入門！まずは前編！ハートリー・フォック法についてお勉強！

基本の復習

まず、基本の復習です。とりあえずシュレーディンガー方程式が解ければ、その分子がどんな感じのやつかわかるんだ、と！

で、「ハミルトニアンが決まるのが大事」ということですが、 どうも「ハミルトニアンはエルミート演算子」ということに関連しているらしい。

「固有値が実数だから観測量として意味をもつ」、ということでしょうか？

これを踏まえてもう一度定常状態のシュレーディンガー方程式を見返します。こんな感じ？

・・・エルミートってそんな物理化学的な意味合いにつながってたんですね。 線形代数の格好いい名前だけど、なんだかよくわからないやつくらいにしか思ってませんでした。。。

では、この大事なハミルトニアンをどう導くか？

「古典的なハミルトン関数をつくっておいて演算子を使って書き直す」ことで導出できるそうです。 以下のような「量子化の手続き」と呼ばれる対応規則を用いればOK！！簡単！！

分子のハミルトニアンの式は長いので省略します。（・・・LaTexにもう飽きた）

雰囲気量子化学の流れ

さて、本題。水素原子からDFTへの穴埋めです。

あやふやな雰囲気ですが、キーワードを拾っていくとこんな感じみたいです。

多粒子問題のシュレーディンガー方程式を解けないので、近似を頑張って1粒子問題のハートリーフォック方程式までもっていった。 でも、どうしても誤差（電子相関）の問題が残った。解決のためにポスト・ハートリーフォック法が考えられたが、計算コストがとても大きくなった。

で、より計算コストの低い解決策が密度汎関数法(DFT)で、「波動関数ではなく電子密度から出発する」という根本的な違いがある。 DFTが解くのはシュレーディンガー方程式そのものではなく、等価な別のもの。原理的には厳密に電子相関を見積もることができるらしい。

ただDFTにも「汎関数の正確な形がわからない」という問題があり、近似が導入される。現在のDFT計算の多くはコーン・シャム近似に基づいており、 コーン・シャム法では汎関数の運動エネルギー項のためにコーン・シャム軌道を、また交換相関汎関数と呼ばれる項を導入した。 *1

で、この交換相関汎関数として最も有名なものにB3LYPがある。

やった！B3LYPでてきた！

さっぱり意味がわかりませんが、とりあえずこんな感じに追っていけば論文でよく見るアレにたどり着ける！

ハートリー・フォック法までの流れ

では、前半シュレーディンガー方程式〜ハートリー・フォック方程式までの流れをもう少し詳しく追って見ましょう。

こんな感じ。

ボルン・オッペンハイマー近似と分子軌道

多原子分子のシュレーディンガー方程式は厳密には解けないので近似が必要です。 近似法の一つとして分子軌道法があり、その基礎としてボルン・オッペンハイマー近似（≒断熱近似）があります。 これは「電子の運動に対して原子核の運動を固定させて考えよう」というもので、原子核と電子を分離することで、 「原子核と電子の多粒子問題」を「電子のみに着目した問題」へと簡略化することができます。

「原子マジで重いしもう止めて良くない？？」ってやつですね！

「電子のみ」となりましたが、依然として多電子系は３体以上の多体問題なのでさらに近似が必要です。 ここで導入されるのが分子軌道（Molecular orbital, MO）で、「一つの電子の座標だけを含む１電子軌道関数」です。 分子軌道の概念をもちいることで「１電子の問題」にまで近似することができます。

ちなみに、電子の座標には位置の座標だけでなく電子スピンの座標も含まれます。

MOが出てくると実験化学屋でも親しみを感じられますね！光れ！HOMO-LUMO！

パウリの排他原理とスレーター行列式

さて１電子の問題まで近似できましたが、ここでさらに「波動関数はパウリの排他原理を満たさなければならない」という問題がのこってます。 パウリの排他原理は「２つの電子の座標を交換すれば、波動関数の符号が変わる」というもので、「電子の交換に関して波動関数が反対称である」といいます。

分子軌道法の最も簡単な形、ヒュッケル分子軌道法では「n電子系の波動関数をn個の分子軌道の積」として近似します。 ですが、このハートリー積と呼ばれる近似では、電子の座標を交換しても波動関数の符号が変わらずパウリの原理を満たせません。

ここで、スレーターは行列式の性質「任意の2つの行(or 列)を入れ替えると符号が変わる」に着目し、波動関数を行列式で表しました（スレーター行列式）。 行列式波動関数 (Slater determinant)は自動的に反対称化されているためパウリの原理を満たすことができます。 *2

変分法とハートリー・フォック方程式

さて、近似と波動関数の表し方が決まりました。あとは変分法によりハートリー・フォック方程式を得ることができます。

変分原理によると「近似のエネルギーは真のエネルギーと等しいか、それよりも高く」なります。 従って「最良の近似波動関数は近似のエネルギーを最小にするように選べば良い」ということになります。

以上を踏まえて、スレーター行列式について、波動関数を構成する軌道関数に関してエネルギー期待値を最小化すれば最も良い近似式となります。 ラグランジュの未定乗数法を使って規格直交条件のもとで変分法を適用することで導くことができ、これがハートリー・フォック方程式です。 *3

ハートリー・フォック法関連用語

ハートリー・フォック(HF)方程式の導出までの流れを辿りました。すでに頭の中がぐちゃぐちゃなので、整理を兼ねて関連用語を拾ってみましょう。

• 独立粒子モデル

「一つの電子に着目して、その電子は原子核と他の電子の作る平均場の中を運動している」という描像で、HF法はこの考え方に従っています。 電子一つとする近似で、問題の簡略化に繋がっていますが、同時に「２つ以上の電子がお互いに衝突したり、散乱したりする効果を含まない」という点で誤差の要因となります。*4

• ab initio 分子軌道法

HF法では経験的なパラメーターを用いないので、HF法に基づく分子軌道法をab initio分子軌道法や非経験的分子軌道法と呼びます。

• 軌道エネルギーとハートリー・フォックエネルギー

HF法のエネルギー（ハートリー・フォックエネルギー）はエネルギーの期待値です。これに対してHF方程式のフォック演算子の固有値は軌道エネルギーとよばれます。 ヒュッケル法の場合と異なり、「HFエネルギー ≠ 全軌道エネルギーの和」であることに注意です。

• 空間軌道とスピン軌道

分子軌道の電子の座標には位置の座標と電子スピンの座標があります。 前者による空間部分を空間軌道とよび、空間軌道にスピン関数をかけたものをスピン軌道といいます。 スピン関数は規格直交条件をつかって積分できるので、スピン軌道によるHF方程式から、空間軌道だけのHF方程式を導くこともできます。

• RHF法とUHF法

制限付き(Restricted)HF法と非制限(Unrestricted)HF法の違いはスピンの取り扱いによります。
RHF法では「 αスピンの電子とβスピンの電子をペアにして、同じ空間軌道につめる」ことで得たスピン軌道を使います。 対して、UHF法では「ペアにせず、異なる空間軌道**につめる」ことで得られたものを使います。

ハートリー方程式を解くには？

さて、シュレーディンガー方程式の近似式であるハートリー・フォック方程式を得ることができました。これを解くにはどうすれば良いでしょうか？

ここまでで得たHF方程式は、フォック演算子の運動エネルギー項に微分を、クーロン演算子と交換演算子に積分を含む微積分方程式であり、 多原子分子のような多中心の場合、数値的に解くことは困難です。 そこで行列形式のハートリー・フォック・ローターン方程式に誘導し、self-consistent field (SCF)の手続きと呼ばれる手法により解く、ということが行われます。

この行列問題への置き換えはDFT計算においても主流の解き方となっているそうです。

ハートリー・フォック・ローターン方程式

ハートリー・フォック・ローターン方程式を導くには、まず基底関数を用いて分子軌道をLCAO (Linear Combination of Atomic Orbital)展開します。

LCAO展開は、分子軌道を「N個の原子軌道の線型結合」で表すもので、「分子の軌道は原子軌道の形を色濃く残している」といったタイトバインディングの考えに基づく近似です。 この時の原子軌道を基底関数、展開係数を分子軌道係数と呼びます。

HF方程式導出の時と同様、変分原理に従って、近似エネルギーが最小となるような近似式が最良の波動関数となります。 LCAO展開では分子軌道が「線型結合で近似した試行関数」となっているので、リッツ (Ritz)の変分法を用い、 変分パラメーターである展開係数(分子軌道係数)を求めれば良いことになります。

Ritzの変分法では永年方程式と呼ばれる方程式を解くことで展開係数を得ます。永年方程式は行列式なので、ここで行列問題への置き換えができていることになります。

ちなみにローターンのつづりはRoothaanです。ジョジョの擬音っぽい。ウリィィィィィ！！

self-consistent field (SCF)の手続き

さてハートリー・フォック・ローターン方程式に導くことで行列問題となりました。

ハートリー・フォック・ローターン方程式は「重なり行列をもった一般化固有値問題」の形を取っています。

単純な場合であれば、重なり行列(S)の単位行列への変換（直交化）*5により、「重なり行列をもたない普通の固有値問題」へと変換することで、簡単に計算機で解くことができるようになります。

ですが、ハートリー・フォック・ローターン方程式の場合、フォック行列に分子軌道係数からなる密度行列が含まれているため単純ではなく、 「繰り返しの手続き」を用いて計算することになります。この手続きをself-consistent field (SCF)の手続き（自己無撞着場の手続き、つじつまのあった場の手続き）と呼びます。

上図の手続きをn回繰り返し、n回目と(n-1)回目の係数行列 or 密度行列がほとんど同じものになれば、SCFが収束(converge)したことになり、 ハートリー・フォック・ローターン方程式が解けたことになります。

基底関数

無事ハートリー・フォック・ローターン方程式の解き方までわかりました。 話を戻して、先のLCAO展開で導入された基底関数の中身について触れておきます。これは後編のDFTにも関係してきます。

基底関数の形

まず、基底関数の形としては大きくスレーター型(Slater-type)とガウス型(Gaussian-type)があります。 スレーター型は水素原子の波動関数と同じ形で、より正確ですが計算コストが高く、現在では半経験的分子軌道法計算以外ではあまり用いられていないそうです。 対するガウス型は正確さに劣るものの、2電子積分を解析的に求めることができ、計算コストでメリットがあるため、基底として多用されているようです。

上図で式の記号はそれぞれ、主量子数(n)、方位量子数(l)、磁気量子数(m)、球面調和関数(Y)を表します。 また、R_Aは原子Aの座標ベクトルで、軌道指数(ζ)は関数の広がりを表します。

カスプ条件（Wikipedia 加藤の定理）とは、「電子密度が原子核の位置において尖点(カスプ)を有することを述べる」ものです。 スレーター型、ガウス型、両者のグラフを見比べれば、前者はx=0で尖っていてカスプ条件を満たし、後者は滑らか（微分が0）で条件を満たさないことがわかります。

短縮ガウス型基底関数

さて、多用されるガウス型軌道(GTO)ですが、より少ない個数の関数系で精度の高い計算を実現するため、短縮ガウス型基底という形で用いられます。

以下のような形を取っており、軌道指数の異なる原子ガウス型関数の線型結合で与えられます。 「GTOを複数組み合わせることで近似的により正確なSTOに近づける」ことが目指されています。

代表的な基底関数系として以下のものがあります。

いよいよ6-31Gがでてきました！よく見るけどよくわからないやつ！

格好いいのでもう少し用語を拾います。

pople型基底関数系

こういうX-YZGといった形で表されるsplit-valence基底系は、ジョン・ポープル（Pople)らのグループにより提案されたのでPople型基底関数系とも呼ばれるようです。
ハイフンの前半は内殻軌道、後半は原子価軌道(valence ornital)の表し方を示し、後半の数によって呼び方が変わります。*6
例えば、6-31GのようなX-YZGはダブルゼータ(DZ)型、6-311GのようなX-YZWGであればトリプルゼータ(TZ)型と呼ばれます。 *7

具体例として、炭素原子では6-31Gは以下に対応します。*8

• 内殻軌道 = 1s → ６個の原始(primitive)GTOからなる短縮(contracted)GTO(CGTO)で表す
• 原子価軌道 = 2s2p →　3個の原始GTOを短縮した短縮GTOと、1個の原始GTOの線型結合で表す

原子価軌道は短縮GTOと原始GTOの組み合わせとなっていますが、グラフ形状をみると後者の原始GTOのほうが動径方向に広がっており、 前者が軌道のより内側、後者がより外側も表現する役割を担っているようです。*9

Pople型以外でよく使われる基底系としては、Huzinaga-Dunning系、Roos ANO系、Dunning cc-pV nZ系といったものがあるそうです。 最後のダニング（Dunning）らによる「Correlation-Consistent（電子相関と整合な）基底系」（ex. cc-pVDZ）も論文でよく見かける気がします。

分極関数、分散関数

さて、6-31Gがわかりました。では、変化形の6-31G(d)や6-31+G(d)とは？？？

それぞれ、「(d)は分極関数(polarization function)」、「+ は分散関数(diffuse function)」を表すそうです。

分極関数は「一つ多く節をもつ補助関数（角運動量が高い）」で、これを追加することで基底系に柔軟性が加わり、 分子中の原子の電子密度の分極を表現できるようになります。

例えば水素原子であれば、球対称な1sに分極p関数を追加することで、非対称な分子軌道も表せる、といったことになります。

分極関数の表記は追加した軌道をそのまま書くか、アスタリスクをつけます。*10

一方、分散関数（広がった関数）は「非常に平たいガウス関数」で、原子軌道の原子から遠く離れた尾部を正確に表現することができるようになります。 例えば、アニオンや励起状態の電子といった、より広がった電子分布を表現する際に有効な補助関数となるそうです。*11

分散関数の表記は「Gの前にプラス(+)」をつけます。「++」だとより軽い元素（HとHe)にも分散関数が追加されることを示します。

・・・表記１文字にそんな違いがあったとは。。。知りませんでした。

以上、基底関数系の用語紹介でした。

前編まとめ

以上、雰囲気量子化学入門（前編）では、DFTの知ったかぶりに向けてシュレーディンガー方程式の復習からハートリー・フォック方程式とその解き方までを辿って見ました。

ざっくりとまとめると、多体問題解くために、近似して、行列表現にしたよ！って感じでしょうか？それにしてもどうしてこんな発想がでてくるのでしょう？本当すごいですね。

最後に、分子軌道を表す要素として基底関数系についても触れました。よくみる格好いいアレの正体がちょっとわかって楽しかったです。

っていうか脳みそシンプルなので分極関数の形を見て、「え、これ混成軌道と同じじゃないの？そんなちょいちょい足さんでも、もう混ぜてしまえばええんやないの？」って思ってしまいました。たぶん正確さ（定量性？）とか違うんですよね？？？そもそも分子軌道は１電子の近似だから、化学結合の原子価結合法とは別物なのでしょうか？さっぱりわからない。

あとPople型でゼータと呼ぶのがなぜかもわかりませんでした。唯一分かったのはエルミートには格好いいだけじゃない意味があったということ！

格好つけるために数式をLaTeXでコピペしてみましたが、意味はわからなかった！*12。

ひどい誤解がたくさんあると思うのでご指摘いただければ幸いです。

参考資料

この記事は主に「量子化学　ー分子軌道法の理解のためにー 中島隆人　著」（裳華房）の内容を参考にしています。 かつて優秀な研究室の先輩にDFTについて知りたいといったら、「とても簡単にわかりやすく紹介されているから」とオススメされました。 私の能力では字面を追うだけで精一杯でしたが、数式の変換も記載してあり、コンパクトに幅広い情報がまとまっているので、 入門書としてとてもよいのではないかな？と思います。ご興味を持たれた方はぜひご一読を！そして、私に解説してください。お願いします。

また、インターネット上に公開してくださっている以下の資料がわかりやすく、参考にさせていただきました。

• @dc1394氏による解説
  • できるだけ簡単にHatree-Fock法を解説してみる（C++17のコード付き）
  • できるだけ簡単にHatree-Fock法を解説してみる part2
  • できるだけ簡単に密度汎関数理論（Density Functional Theory, DFT）を説明してみる（前編）
  • できるだけ簡単に密度汎関数理論（Density Functional Theory, DFT）を説明してみる（後編）
• 理化学研究所　常田グループ（理論化学）　講義資料
• 計算科学.com
  • Pople型基底関数
• PC CHEM BASICS.COM
  • B3LYP/6-31G(d)とは何なのか？
  • 基底関数はどれを選べばいいの？
• Chem-Station
  • 計算化学：汎関数って何？
  • 計算化学：基底関数って何？
  • 計算化学：DFTって何？ part II
  • 計算化学：DFTって何？ part III
• wikipedia
  • 基底関数系（化学）)

おまけ

念のため、観測量に関連して「演算子Aの期待値」の定義を復習します。ついでに記号が似てるのでブラケット表現も。 だいたいこんな感じ。

*1: コメントでご指摘いただいた内容をもとに修正してみました 01/18 update

*2:スレーター行列式には「反対スピンを持つ電子が同じ軌道を占めることが可能になる」という欠点があり、１電子軌道関数に反します。 同じ向きのスピンを持つ電子については、行列式の「同じ成分を持つ行(or 列)が２つあると値が0になる」という性質から、同じ軌道を占めることができません。

*3:以下の図では正準(cononical)ハートリー・フォック方程式としています。これはHF方程式をラグランジュの未定乗数法で解く際に、未定乗数が対角行列となるようにユニタリ変換したものです。正準HF方程式はフォック演算子の固有値方程式であり、固有値が軌道エネルギー(εi)、対する固有関数がスピン軌道(φi)となります。「HF方程式の解をユニタリ変換してもHF方程式の解となる」こと、「フォック演算子がエルミート演算子であるため、対角化するようなユニタリ変換を選ぶことができる」ことから、このような選び方が可能になります。ここでもエルミートが活躍！（参考 Wikipedia-ハートリー＝フォック方程式 ）

*4:後編で出てくる電子相関の問題です

*5:重なり行列はエルミート行列なので、ユニタリー変換により単位行列に変換することができます

*6:この数は「異なる軌道指数を持つ関数のセットの数」を表すそうです。 参考 Chem-Station 計算科学：基底関数って何？

*7:Wikipediaに「ゼータ ζ はSTO基底の指数部を表わすために用いられることが多い」と書かれていますが、なぜそれでPople型の呼称に使われるのかはよくわかりません。どなたか教えてください

*8:使われている原子GTOの数をまとめると、1sに4個、2Sに4個、2pに4個なので、結局s軌道に10個、p軌道に4個使われていることになります。
また、短縮後の関数は1sに１個、2s2pに2個なので、結局s軌道に3個、p軌道に2個使われています。 これらをまとめて「(10s4p)→[3s2p]」といった表記もなされるそうです。

*9:参考 1. 計算科学.com Pople型基底関数
参考 2. Wikipedia 基底関数系（化学）)
参考 3. PC CHEM BASIS.COM B3LYP/6-31G(d)とは何なのか？

*10:軌道をそのまま書くのとアスタリスクで微妙に意味が違うようなのですがよくわかりませんでした。

*11:参考 PC CHEM BASIS.COM 基底関数はどれを選べばいいの？

*12:LaTeXiTとかいうのすごい便利ですよ！