pfpoly：uMLIPと反応加速を組み合わせた即時利用可能な高分子反応シミュレータ

タイトル：pfpoly：uMLIPと反応加速を組み合わせた、すぐに使えるポリマー反応シミュレーター

はじめに

こんにちは、Preferred Networksの森穂高です。普段はMaterials&DrugDiscovery事業本部のマテリアルリサーチグループでエンジニアリングマネージャを務めています。

この記事では、先日出版された私たちの論文「Ready-to-Use Polymerization Simulations: Combining Universal Machine Learning Interatomic Potential with Time-Dependent Bond Boosting for Polymer and Interface Design」で提案した、重合・硬化反応を扱う分子動力学（MD）の枠組みを紹介します。

ポリマー材料は、最終的にできあがる分子構造（鎖長、分岐、架橋、界面での結合など）によって物性が大きく変わります。実験では反応の進行や平均的な構造は観測できても、反応中間体や界面近傍での結合再編成を原子レベルで追跡し切るのは難しいことが多いです。こうした状況では計算化学の活用が考えられますが、反応を含むMDを現実的なコストで実行するのは依然として困難です。

今回の論文では、反応を含むMDを「とりあえず実行して、反応機構と相対的な傾向を確認できる」状態にまで持っていく手法を提案しています。絶対的な速度定数を求めることよりも、材料設計で必要となる比較や探索に使えることを目指しています。この論文のコンセプトムービーがYouTubeで公開されていますので、こちらもご覧ください。

重合・硬化のMDが困難な理由

重合・硬化をMDで扱う際の難しさは、主に2つあります。

1つ目は、反応を扱えるポテンシャルの準備が煩雑であることです。ReaxFFのような反応力場は便利ですが、対象系や反応の種類が変わると、パラメータセットの選定や妥当性確認が必要になります。特に、接着や摩擦などの実用的な課題で頻繁に登場する有機-無機界面では、有機物の反応や分子の凝縮構造、無機物の結晶構造系といった幅広い物理化学的性質を両立できるポテンシャルが必要です。このようなポテンシャルを調査・検証し、利用可能な状態にまで整えるには時間がかかりがちです。

2つ目は、反応がレアイベントであり、通常のMDでは反応が進まないことです。一般に反応の発生頻度は活性化エネルギーに対して指数的に減少します。例えば、ラジカルカップリングは非常に速い反応として知られていますが、そのような反応であっても発生頻度はナノ秒オーダーであり、全原子シミュレーションの典型的な時間刻みであるフェムト秒と比べると、100万ステップが必要となるレアイベントとなります。重合は逐次反応の連続なので、ナノ秒程度のMDでは反応イベントが十分に起こらず、統計を取る前に計算が終わってしまうことが多いです。反応を起こすために加速法を導入すると、今度は反応ごとにパラメータを調整する手間が生じたり、正逆両方向の反応が加速されてしまい重合のような一方向性のプロセスとの相性が悪かったりします。

1つ目のポテンシャルに関しては、ここ数年でDFTデータで学習した汎用機械学習ポテンシャル（uMLIP）が現実的な選択肢になってきました。Matlantisで提供されているPFP、M3GNet、CHGNet、MACEなどが代表例です。これらは「系ごとに力場を作り込まない」という点で非常に魅力的ですが、uMLIPに置き換えただけではレアイベントの問題は解決しません。そこで今回の論文では、uMLIPの汎用性と、反応を起こすための加速を一体的に扱うことにしました。

割り切ったアプローチ：Time-Dependent Bond Boost

基本方針は、uMLIPとしてPFPを用い、反応が進まない部分はボンドブーストで促進するというものです。ここで従来のボンドブーストと異なるのは、バイアスの深さを固定しない点です。Vashisthらの先行研究では、例えば次のような形で深さ \( f_1 \) が固定されています。

$$

V(r) = f_1{1-\exp[-f_2(r-r_0)^2]}

$$

ここで、\( V(r) \) は対象原子対の距離 \( r \) に応じて付加するバイアスポテンシャル（エネルギー）で、\( r \) は対象原子対の原子間距離です。\( f_1 \) はバイアスの振幅（エネルギーの次元を持ち、バイアスの強さ＝最大の寄与スケールを規定する）、\( f_2 \) は距離依存性の幅を決める正のパラメータ、\( r_0 \) は距離の基準（閾値）です。しかし、反応の障壁は反応タイプごとに異なるため、適切な \( f_1 \) を事前に見積もって調整する流れになりやすく、反応の種類が増えるほど準備が面倒になります。そこで我々は、バイアスの振幅を時間依存にし、反応が起きない間は単調に大きくしていくことにしました。このような時間依存のバイアスポテンシャルをTime-Dependent Bond Boost (TDBB)と呼んでいます。論文では

$$

f_1(t)=\min(\gamma t, f_1^{\max})

$$

という形を用い、初期は弱く、反応が起きなければ徐々に強くして「いずれは起こす」設計にしました。ここで \( t \) はシミュレーション時間、\( f_1(t) \) は時刻 \( t \) におけるバイアス振幅、\( \gamma \) は振幅の増加率、\( f_1^{\max} \) は振幅の上限値です。\( \min(\cdot,\cdot) \) は2つの引数のうち小さい方を返す関数です。固定された深さに頼るのではなく、時間方向に自動でスキャンしていくという感覚です。もう一つ重要なのは、結合形成と結合解離で異なる関数形を使い、生成物側（あるいは解離側）に到達したらバイアスが消えるようにしたことです。結合形成側には

$$

V_f(r,t) = f_1(t){1-\exp[-f_2(r-r_0)^2]}

$$

を、結合解離側には

$$

V_d(r,t) = f_1(t)\exp[-f_2 r^2]

$$

を用いています。ここで \( V_f(r,t) \) は結合形成を促進するための時間依存バイアスポテンシャル、\( V_d(r,t) \) は結合解離を促進するための時間依存バイアスポテンシャルです。いずれも距離 \( r \) と時間 \( t \) に依存し、振幅は \( f_1(t) \) で与えられます。こうすることで、反応が起きた後にもバイアスが残り、結合を解離し続ける状況を避けています。距離の閾値 \( r_0 \) は、反応ごとに設計し直すのではなく、対象原子のファンデルワールス（vdW）半径の和にスケール因子 \( \lambda \) を掛けた

$$

r_0 = \lambda \sum_{a\in{\text{pair}}} r_a^{\mathrm{vdw}}

$$

で決めています。ここで \( r_a^{\mathrm{vdw}} \) は原子 (a) のvdW半径（ファンデルワールス半径）で、上付きの \( \mathrm{vdw} \) はvdW半径であることを示すラベルです。\( \lambda \) は無次元のスケール因子、\( \sum_{a\in{\text{pair}}} \) は対象原子対 \( \text{pair} \) を構成する2つの原子 (a) についての和を表します。ここも「反応ごとの手作業を減らす」ための工夫です。この方法の立ち位置は割り切っています。バイアスは時間とともに大きくなり、最終的に障壁を超えることもあり得るので、ハイパーダイナミクスのような物理時間の再構成（タイムリスケーリング）は前提にしません。したがって、この枠組みは絶対的な速度定数を求める用途ではなく、反応機構の理解と相対比較を主目的としています。論文中でも、その理由（タイムリスケーリングの成立条件が一般に満たされないこと）を明示しています。

実装とシミュレーションフロー

uMLIPにはPFP v6.0.0（CRYSTAL PLUS D3 mode）を用い、MDエンジンにはOpenMM 8.1.1を使用しました。バイアスポテンシャルはOpenMM-Torchを通じてTorchScriptとして組み込み、力はPyTorchの自動微分で評価しています。このように自動微分可能な形でエネルギーを記述することで、独自のバイアスポテンシャルを実装しやすくしています。自動微分で力を求める実装は、一般に解析微分による実装よりも低速ですが、ポテンシャルとしてuMLIPを用いる場合はほとんどのケースでそちらがボトルネックになるため、今回の用途ではバイアスポテンシャルがパフォーマンスに与える影響は軽微です。実際のシミュレーションフローを図1にまとめています。最初に官能基に基づいて反応性原子のグループを定義し、距離条件で反応候補の組を作り、原子が重ならない候補だけを選んでバイアスポテンシャルを付与します。バイアス付きMDで反応を誘起した後はバイアスなしで緩和し、結合状態の変化を判定して反応性グループを更新し、次の反応へ進みます。結合判定はvdW半径に基づくルールを採用し、距離がvdW半径和の60%未満になったら結合とみなす、という形に統一しています。

image図1 シミュレーションフロー

ケーススタディ

この枠組みは、特定の反応一つを精密に再現することよりも、化学的に異なる系を同じ手順で実行し、傾向と反応機構を把握できることを重視しています。そこで論文では、ラジカル重合、縮合重合、界面でのエポキシ硬化の3つの系で検証しました。ここではその中でも代表的な例として2つを紹介します。

連鎖重合：ラジカル重合における濃度依存性と分子量成長の傾向

複数のビニルモノマー（メタクリル酸メチル、メチルメタクリレート、スチレン、酢酸ビニル、ジフェニルエチレン、イタコン酸ジメチル）について、モノマー転化率の時間変化を単一指数関数でフィットして見かけの成長反応速度定数 \(k_p^*\) を求め、実験値の順位と比較しています。その結果、Spearman ρ=0.66、Kendall τ=0.60、順位の平均絶対偏差（MAD）=1.0という相関が得られました。完全一致ではありませんが、チューニングなしで相対的な傾向がある程度得られるという点で、狙い通りの結果と考えています。一方で、PFPが推定する障壁の順位と実験傾向が一致しない例（スチレンの順位など）も見られます。これは率直に弱い点であり、凝縮相では活性化障壁だけでなく拡散や配座、エントロピー要因も効くため、障壁の誤差と \(k_p^*\) の誤差が一対一で対応しない点についても議論しています。そもそも時間依存バイアスによって時間軸が歪むため、絶対的な反応速度を求める設計ではなく、順位付けや比較に主眼を置いた立ち位置となります。

image図2 様々なビニルモノマーに対する見かけの反応速度定数(a)と実験値との順位比較(b)

固体界面での反応：CuO界面でのエポキシ硬化における界面近傍の反応抑制とCu–O–C結合

最後に、CuO（OH終端）基板上でのエポキシ-アミン硬化を扱っています（図3）。この系は、有機硬化反応と酸化物表面反応が共存するため、反応力場でバランスを取るのが難しい類の問題だと考えています。TDBB＋PFPはこうした問題に対応できます。シミュレーショントラジェクトリを解析して得られた深さ方向の反応密度（図3 (a)）を見ると、界面近傍はエポキシ・アミンの密度が高いにもかかわらず反応密度が低下しています。これは、界面近傍では運動性が低下して反応が抑制されるという理解と整合する結果です。スナップショットからも、CuO–OH表面での開環反応を介して、ポリマーネットワークと基板との共有結合（Cu–O–Cを含む結合）が形成される様子が確認できます（図3 (b)）。数は多くありませんが、界面接着に寄与する結合が生じるという点で重要です。

image図3 CuO界面でのエポキシ硬化における反応密度プロファイル(a)とスナップショット(b)

手法上の制約と今後の課題

TDBBは「とにかく実行できること」に重点を置いた設計であるため、代償もあります。まず、時間依存なバイアスポテンシャルによって物理時間との対応関係が崩れるため、絶対的な速度定数を再現する用途には向きません。次に、ラジカル重合については現状の実装では停止反応・連鎖移動反応を考慮していないため、分子量分布まで含めて再現したい場合は反応カタログの拡張が必要です。最後に、相対的反応性の細かい順位付けは、uMLIPの反応障壁再現精度に上限があり、数kcal/molレベルの差を安定して見分けるのは現状では容易ではありません。

今後の方向性としては、遷移状態や反応障壁の情報をより明示的に含むデータでuMLIP側の障壁再現精度を向上させつつ、副反応や停止反応も含む形で反応カタログを拡張し、必要に応じて加速と時間再構成（キネティック・リコンストラクション）を組み合わせて、傾向の観察から定量的な予測に近づけていくのが自然な流れだと考えています。

おわりに

今回の研究は、反応力場を系ごとに作り込む方向性というより、uMLIPを前提として「反応が進まない」というボトルネックだけを、できるだけ手間をかけずに解決する方向を目指しました。TDBBは、物理時間の正確さを捨てる代わりに、反応ごとの事前チューニングをできるだけ避け、反応機構の理解と相対比較を行えるようにします。重合・硬化を材料設計のワークフローに組み込む上では、まず「実行できて、結果が見える」ことが重要だと考えており、その最低ラインを超えるための提案がPFPとTDBBの組み合わせです。

Materials & Drug Discoveryでは、Matlantisを超えて材料開発を加速する様々な技術を研究しています。一緒に研究開発に取り組む研究者/エンジニアを募集しています。この記事を通じて興味を持たれた方は、以下のページをご覧ください。今回の記事のような原子シミュレーション手法に関する知見をお持ちの方も大歓迎です！

2027新卒採用 – マテリアルディスカバリーエンジニア / リサーチャー

マテリアルリサーチPE（Matlantis事例）

マテリアルディスカバリーリサーチャー

化学系R&D・企画出身者の新規事業開発

コンサルタント出身者の新規事業開発

The post pfpoly: A Ready-to-Use Polymer Reaction Simulator Combining uMLIP and Reaction Acceleration appeared first on Preferred Networks Tech Blog.