見出し画像

（概要編）生成AIが拓く新たな可能性：自然言語でDemand-Drivenな量子化学計算

QunaSys

はじめに

QunaSysのCRSチームでエンジニアをやっている山口です。
私は学生時代、分子動力学を用いた研究に従事していました。
現在、QunaSysには量子化学計算に精通しているメンバーが多数在籍していますが、私自身は量子化学計算について多少心得てはいるものの、リサーチチームでもないため、本業としては取り組んでいません。
しかし、QunaSysの一員である以上、量子化学計算に対する理解は不可欠だと考えています。
ただ量子化学計算は非常に複雑で、理解を深めながら計算を進めることは容易ではありません。例えば、密度汎函数理論(DFT)での計算をしようとした時にどの汎関数を選択するか、計算結果をどうやって解釈するか、などは難しい判断が求められます。
そこで、私のように量子化学計算に詳しくない人でも、やりたいことベースで簡単に計算を実行し、結果を容易に解釈できるようにするため、生成AIを用いた量子化学計算の自然言語インターフェイスの作成に挑戦してみました。
生成AIが「私たちのリクエストをどのように解釈して計算を実行するのか」を確認することで、量子化学計算に対する理解度も大きく深まることが期待できます！

Demand-Drivenな量子化学計算の実行デモ

まずはデモ動画でどんなものか体感してみてください。
3本あります。楽しんでみてくださいね！

計算前提

今回のデモ動画では以下の前提で計算されています。

  • 計算手法：PySCFでDFT計算

  • 基底関数：sto-3g

    • 最適化計算などは6-31Gもしくは6-31G(d)などを使うと思いますが、今回はサッと計算できるように一番小さいsto-3gにしました。

では早速デモ動画を見ていただきましょう！

構造最適化デモ動画

  • ざっくりデモ内容：

  1. リクエスト「CC=Cの最安定構造が知りたい。」をLLMが理解して最適な計算を方法を選択 → 構造最適化計算

  2. リクエスト内容から計算する対象分子を用意 → リクエスト内容からCC=Cを抽出

  3. 分子の内容から汎関数を選択 → B3LYP

  4. 1 ~ 3までの内容を使って分子の構造最適化を実行 & 可視化

  5. 構造最適化の結果からLLMが解釈し、回答を生成してます。

  • LLMが構造最適化した結果を解釈した内容：

質問：CC=Cの最安定構造が知りたい。
回答：CC=C（単純なアルケン、エテンを指すと思われる）の分子の最も安定な構造に関する量子化学計算の結果を解釈するには、通常、量子化学計算の出力からいくつかの重要な情報を見ます。 しかし、./optimization_steps.xyz のような外部ファイルに直接アクセスしたり、解析したりすることはできないので、このようなファイルで通常何を調べ、エテンの最も安定な構造に関する結果をどのように解釈するかを説明します。 軌道ファイルで考慮すべきポイント エネルギーの最小化 軌跡ファイル(optimization_steps.xyz)には、分子の最小エネルギー配置を見つけるために行われたジオメトリ最適化ステップに関する情報が含まれている可能性があります。エネルギーが最小になるステップを探すと、最も安定な構造に対応します。 ジオメトリの収束： ジオメトリが完全に最小に収束した兆候を最適化ログで確認します。これは通常、勾配が小さく、変位が小さ いことを確認し、その後の最適化ステップで形状やエネル ギーが大きく変化していないことを示します。 結合の長さと角度： エテンの場合、最も安定なジオメトリは2つの炭素原子間の二重結合を特徴とします。C=C結合長は約1.33Å（オングストローム）です。各炭素上のH-C-H結合角は、sp^2混成の特徴である約120度であるべきである。 振動数： これは、構造が最小であり、鞍点（安定分子ではなく遷移状態を示す）ではないことを確認する。 可視化： 分子可視化ソフトウェアを使って.xyzファイルを開き、構造を視覚的に調べることができる。これにより、構造が妥当であり（例えば、正しい結合長と角度）、エテンの予想される形状に対応していることを確認できます。 結果の解釈 エネルギー： 最適化の軌跡の中で最もエネルギーが低い構造は、 モデル化された条件下（気相、 溶媒モデルが適用されている場合など）で最も安定なエテンの形を表します。 形状： エテンの予想される形状は、sp^2混成による各炭素原子の周りの三角平面配位で、結合角は約120度です。 結合： 炭素原子間の二重結合は明らかであり、結合次数は2、結合長は単結合のC-C結合より短い。 軌跡ファイルのこれらの側面を注意深く調べることで、エテンの最も安定な構造を確信を持って特定することができます。量子化学計算の真価は、単に最も安定な構造を見つけることだけでなく、構造、エネルギー、物性の関係を理解することにあることが多いことを忘れないでください。

振動解析デモ動画

  • ざっくりデモ内容：

  1. リクエスト「c1ccccc1の振動モードを知りたい。」をLLMが理解して最適な計算を方法を選択 → 振動解析

  2. リクエスト内容から計算する対象分子を用意 → リクエスト内容からc1ccccc1を抽出

  3. 分子の内容から汎関数を選択 → B3LYP

  4. 1 ~ 3までの内容を使って分子の振動解析を実行 & 可視化

  5. 振動解析の結果からLLMが解釈し、回答を生成してます。

  • LLMが振動解析した結果を解釈した内容：

質問：「c1ccccc1の振動モードを知りたい。」
回答：化学式からベンゼンを表すc1ccccc1という分子の振動数は、典型的な振動スペクトルの値から解釈すると、非常に珍しいものです。一般的に分子の振動周波数はcm^-1（波数）で示され、その範囲は数百から数千cm^-1です。あなたが挙げた値は、予想される値よりも数桁小さく、データの誤解や誤った表示の可能性を示唆しています。

量子化学計算の文脈を考えると、計算された振動数の単位と大きさが適切であることを確認することは極めて重要です。ベンゼンの場合、予想される振動周波数は、特徴的な環の伸縮、屈曲、およびCH伸縮モードを強調するもので、通常、以下の範囲で見られる：

CH伸縮モード：約3000～3100cm^-1
CC伸縮および環変形モード：約1000～1600 cm^-1
面外屈曲モード：約600～900cm^-1
この食い違いを考えると、提供された振動数のスケールが間違っているか、物理的な振動数に直接変換されない計算過程の予備段階である可能性がある。また、これらの値は周波数に正しく変換される前のヘシアン（力定数）マトリックス計算の固有値を表している可能性もある。

振動モードを有意義に解釈するためには、上記の典型的な範囲内の周波数が表示されることが期待される。振動モードファイル（vibrational_mode_1.xyz, vibrational_mode_2.xyz, vibrational_mode_3.xyz）は、各振動モードに対応する分子内の原子の変位ベクトルを表していると思われる。これらのファイルは、各振動中に分子がどのように変形するかを視覚化し、各振動モードの性質（例えば、結合の伸張、角度の屈曲、またはより複雑な運動を伴うかどうか）を理解するのに役立ちます。

要約すると、計算で得られたベンゼンの振動モードを正しく解釈するためには、周波数値の見直しと潜在的な補正が必要です。モードファイルは分子の物理的な振動に対応している必要があり、これを解析することで分子動力学を理解し、ベンゼンの予想される周波数範囲内で振動モードの同一性を確認することができます。

一点計算デモ動画

  • ざっくりデモ内容：

  1. リクエスト「CC=CC=CCの最安定構造が知りたい。」をLLMが理解して最適な計算を方法を選択 → 一点計算

  2. リクエスト内容から計算する対象分子を用意 → リクエスト内容からCC=CC=CCを抽出

  3. 分子の内容から汎関数を選択 → B3LYP

  4. 1 ~ 3までの内容を使って分子の一点計算を実行

  5. 一点計算の結果からLLMが解釈し、回答を生成してます。

  • LLMが一点計算した結果を解釈した内容：

質問：CC=CC=CCが安定かどうかを知りたい
回答：化学式CC=CC=CCで表され、原子位置でさらに詳細が示される分子を、B3LYP汎関数とsto-3g基底セットを用いた密度汎関数理論（DFT）を用いて解析した。得られた結果を解釈してみよう： エネルギー 分子の全エネルギーは-231.78152928244455原子単位（au）である。この値は量子化学計算では絶対的な値ですが、主に分子構造の安定性や反応のエネルギーを評価する際の比較指標となります。 通常、エネルギー状態が低いほど、より安定な配座を示す。参考文献や比較対象がなければ、このエネルギーが高いレベルの安定性を示すかどうかを明確に言うことは難しいが、類似の分子やコンフォーメーションのコンテキスト内で、相対的な安定性を判断することができる。 HOMO-LUMOギャップ HOMO（最高被占分子軌道）エネルギー：-0.15028384231415 au LUMO（最低非占有分子軌道）エネルギー：0.07488075286715457 au HOMOエネルギーとLUMOエネルギーのギャップは、分子の化学反応性と安定性を知る上で重要である。一般にギャップが大きいと、HOMOからLUMOへ電子を励起するのに必要なエネルギーが大きくなるため、分子がより安定で反応性が低いことを意味する。 この場合、ギャップは約0.225165auで、比較的中程度である。このことは、この分子がそこそこの安定性を持ちながらも、特に適切な条件下では反応に参加する可能性があることを示唆している。 解釈 CC=CC=CC分子は、エネルギーとHOMO-LUMOギャップから、この計算の条件下では安定な配置を持つように見える。しかし、適度なHOMO-LUMOギャップのため、その反応性は無視できない。 その安定性をより包括的に理解するには、この分子のエネルギー論を類似構造のものと比較したり、潜在的な反応経路における挙動を評価したりすることが有益であろう。 さらに、ラジカルの形成の可能性や外部条件（触媒の存在、温度など）の影響など、他の要因も考慮することで、この分子の安定性と反応性についてさらなる洞察が得られるかもしれない。 結論 DFT計算から、この分子はB3LYP関数とsto-3g基底セットの制限内で安定な電子構造を持つことが示された。 適度なHOMO-LUMOギャップは、ある程度の反応性を示唆しており、特に有機合成や材料科学における潜在的な用途という観点から、さらなる化学的研究の興味深い候補となる。

QunaSysの社員からのコメント

QunaSysの社員にもコメントをもらいました！

  • 社員1
    　振動解析では、”この食い違いを考えると”というような出力がありますが、これは実験事実と比較した上での計算結果の食い違いを指しているのでしょうか？DFTの精度上、実験値と計算値の振動数は合わないことのほうが多いので、食い違いが生じるのは仕方ないような...という気持ちの上で確認したいと思っています
    　構造最適化では”ジオメトリの収束： ジオメトリが完全に最小に収束した兆候を最適化ログで確認します。これは通常、勾配が小さく、変位が小さ いことを確認し”と書いているけど何に対して小さいかまで書いて欲しい。「指定した収束しきい値以下まで各値が小さいことを確認する」とかだといい感じ。優・良・可・不可の4段階評価を付けるなら良くらいの感覚です。
    　一点計算に関しては違和感は無いです！
    総じてこれが学生の頃にあったらなと思います。

大体はいい感じみたいですね。
これのMD版があると私も嬉しかったかもしれません。

「Demand-Drivenな量子化学計算」でできることの全体像を見る

今回だと「Demand-Drivenな量子化学計算」では以下のような構成になっています。

使ってる技術スタック

  • 言語：Python

  • フレームワーク

    • クライアントサイド：Streamlit

    • サーバーサイド：FastAPI + LangServe

  • 量子化学計算用ライブラリ：PySCF

  • LLMフレームワーク：LangChain

  • LLM：OpenAI

できる計算の種類

計算としてできるのは以下です。基本的にPySCFで計算を実行しています。

  • 一点計算

  • 構造最適化

  • 振動解析

振動モードや最適化の可視化はPy3Dmolを使用しています。

LLMがサポートしてくれている部分

  • リクエストの内容を解釈し計算方法の選択 :

    • ユーザーの自然言語による指示を理解し、適切な計算方法を自動的に選択します。

  • 分子の計算可能な形式への変換 :

    • ユーザーが入力からSMILESや化合物の名前を抽出し、量子化学計算に適した形式に自動的に変換します。

  • 構造最適化の収束判定 :

    • 構造最適化計算の過程を監視し、収束条件を満たしているかどうかを判断します。

  • 汎関数の自動選択 :

    • 計算対象の分子に適した汎関数を自動的に選択します。

Agentという技術を使ってます

「Demand-Drivenな量子化学計算」ではLLMをつかったAgentという仕組みを使っています。
今回は少しだけAgentについて解説させていただきます。

Agentは自然言語処理と推論を組み合わせて、ユーザーからの指示を理解し、適切なアクションを選択・実行することで、タスクを自動的に完了させる仕組みです。
ざっくりいうと、Agentとはモデルに意思決定の機能を与えて機能しているシステムなんですね。

Agentを使っている様子

このAgentを使えば、これまでは人間が手作業で行っていたドメインに特化した専門的なタスクを、AIに自動化させることができるんです。
例えば、化学の分野で言えば、特定の化合物を計算するのに適した密度汎関数を探したいとしましょう。Agentを使えば上記の図ような感じで、目的の論文や汎関数を自動で見つけ出してくれます。

解説すると、以下のようになります。
ユーザーリクエスト：化合物Aを計算するのに密度汎関数理論を使いたい。もっとも適した汎関数のchemRxiv論文を探して、教えてください。

Agentの行動 :
1. 化合物A+DFT+functionalでchemRxiv検索
2. 情報取得してユーザーリクエストでベクトル検索
3. ユーザーリクエストに合わせて、Agentが最適な汎関数の情報を回答として生成

精度に関しては別途考える必要がありますが、Agentではこのようなことができるようになります。
実際に今回作った「Demand-Drivenな量子化学計算」のシステムでもこのようなAgentの仕組みを使用しています。
ドメインに特化したことをLLMにやらせたいならAgent機能は必須になってきますね。

現状の課題と将来の展望

デモ動画の内容からはわかりづらいですが、「やりたいことの解釈と計算方法の選択」や「汎関数の自動選択」などの機能において、まだ精度の面で課題があります。
しかし、これらの課題は、RAG（Retrieval-Augmented Generation）などの技術を導入することで、大幅な改善が見込まれます。

例えば、「汎関数の自動選択」については、計算対象の化合物に関する論文を検索し、そこで使用されている汎関数を選択するという方法が考えられます。
また、「やりたいことの解釈と計算方法の選択」については、ユースケースと計算方法の対応表を作成し、RAGを用いてユースケースを検索しLLMに適切な計算方法を判断させるなどの方法が有効だと考えられます。

さらに、最近提唱されたLLaMP（Large Language Model Made Powerful for High-fidelity Materials Knowledge Retrieval and Distillation）の仕組みを取り入れることで、結晶構造の生成などにも対応でき、より簡単かつ精度よく計算できる範囲が広がると期待されます。

今後は、これらの新しい技術を積極的に取り入れて私がより便利に量子化学計算をできるように「Demand-Drivenな量子化学計算」プロジェクトの精度と利便性をさらに高めていきたいですね！

最後に

もちろんまだまだ拡張の余地がありますが、今回はデモ用に作ったものなのでこの辺でとどめています。
例えば、計算結果の可視化とか、軌道の可視化とかいろいろいれたいなとか思ったり思わなかったり…
次回は「Demand-Drivenな量子化学計算」で使っている自然言語のインターフェイスをどう作るかに関して解説していきたいと思います。よかったらフォローして待ちください！

CRSチームでは今後も化学の研究活動のサポートに向けて調査・開発を進めて行きます。
今回のような「Demand-Drivenな量子化学計算」の仕組みをより知りたい方はお気軽に下記までお問い合わせください！