LAMMPSによるFEP 自由エネルギー摂動法の計算

分子動力学法で自由エネルギーを直接計算することは難しいのですが、いろいろ工夫をすることで計算できるようになります。その一つが自由エネルギー摂動法（Free Energy Perturbation, FEP）です。

溶媒和自由エネルギーを計算する例がJ-OCTAの事例として紹介されていますので、細かい理論はこちらを参照してください。

水への溶解性評価 | J-OCTA解析事例 | CAEソリューション：株式会社JSOL

マルチスケールシミュレーションソフトウェアJ-OCTA による粒子法を用いたプロセスシミュレーション「水への溶解性評価算」の事例をご紹介します。

このFEPの計算はLAMMPSでも実行することができますが、細かいことは置いておいてとりあえず何をどうすれば動くのかを解説してみます。

LAMMPSのサンプルデータ
λで変化する相互作用について
1. pair_coeff
2. 静電相互作用
fix ADAPT all adapt/fep
compute FEP
ステップ数について
結果の解析について

LAMMPSのサンプルデータ

LAMMPSのインストーラを展開し、その中の lammps-29Aug2024/examples/PACKAGES/fep というディレクトリにFEPのサンプルが格納されています。

シンプルな例として、水中のメタン分子の計算 CH4hyd を見ていきます。
（lammps-29Aug2024/examples/PACKAGES/fep/CH4hyd）
GitHubにも置いてあります。

lammps/examples/PACKAGES/fep/CH4hyd/fep01/in.fep01.lmp at develop · lammps/lammps

Public development project of the LAMMPS MD software package - lammps/lammps

CH4hydのFEPデータはfep01とfep10の2つあります。
0：ない状態（λ＝0）
1：ある状態（λ＝1）
として、01は「メタン分子がない状態からある状態に変化する」、10は「メタン分子がある状態からない状態に変化する」を意味しています。

各フォルダにはin.fep01.lmpのようなファイル名のinファイルがあります。inファイルにはrunコマンドが2か所あり、1ファイルで2段階の計算に分けています。①平衡化の計算（10万step）、②Production run（FEP、200万step）という流れです。

λで変化する相互作用について

FEPではパラメータ λ （ラムダ）によってλ=0（ない状態）からλ=1（ある状態）を連続的に変化させますが、ターゲット分子の存在が「ある」「ない」はpairやCoulombの相互作用があるかないかということに置き換えられます。ポテンシャルの強さが0％なら無いのと同じだし、100％なら普通に「ある」ということです。

pair_coeff

pair_style名にsoftと書かれているものはこのラムダの調整がやりやすいです。またサンプルデータではpair_style が hybrid によって2つ設定されていますが、hybridが必ず必要というわけではないです。全部softポテンシャルを設定していいならhybridは要らないはずです。

softが付いているpair_style一覧はこちらにあります。

pair_style lj/cut/soft command — LAMMPS documentation

CH4hydで使用されているのは lj/cut/tip4p/long/soft で、inファイル内では
lj/cut/tip4p/long/soft 3 4 2 2 0.125 1 0.5 10.0 10.0 10.0
という設定になっています。
これは args に otype htype btype atype qdist n alpha_LJ alpha_C cutoff (cutoff2) が書かれています。TIP4P水モデルに特化したものなのでargsが多いですが、上記ページに載っているpair_styleの多くはn alpha_LJ alpha_C cutoff cutoff2という4つのパラメータを指定するようにできています。

n は λⁿ でスケーリングする際のnで1か2を指定、alpha_LJやalpha_Cはそれぞれレナードジョーンズ関数や静電相互作用の関数に入るパラメータα（アルファ）で、alpha_LJは0.5、alpha_Cは10がよく使われるようです。式の解説はこちら

softがないタイプのpair_coeffと比較して、softがあるタイプのpair_coeffには λ のパラメータが追加されています。計算の初期状態で「ない扱い」であれば0.0、「ある扱い」であれば1.0です。

pair_coeff    1    3  lj/cut/tip4p/long/soft   0.1036   3.3279  0.0  # C4H Ow
pair_coeff    1    4  lj/cut/tip4p/long/soft   0.0000   1.0000  0.0  # C4H Hw
pair_coeff    2    3  lj/cut/tip4p/long/soft   0.0699   2.8126  0.0  # H Ow
pair_coeff    2    4  lj/cut/tip4p/long/soft   0.0000   1.0000  0.0  # H Hw

これはメタン-水の相互作用を定義している部分ですが、右端の 0.0 はλ=0.0という意味です。初期状態ではメタンは「ない扱い」になっています。上記はfep01の場合ですが、fep10のinファイルではこれが1.0になっています。

静電相互作用

fep01のinファイル内の1回目のrunの直前に set type 1*2 charge 0.0 というものが書かれていますが、これは原子タイプの1番と2番の部分電荷をゼロにしています。データ上は分子はあるけど、相互作用の上では「分子がない状態」を作るわけですが、そのために静電相互作用も働かない状態にしておく必要があります。書き方の1*2は「1番から2番」という意味です。

fep01ではターゲット分子の電荷をゼロにしておきますが、fep10ではそのようなことはせず、そのままの電荷で計算をスタートします。

fix ADAPT all adapt/fep

FEP計算を行う部分の設定です。

variable lambda equal ramp(0.0,1.0)
variable q1 equal -0.24*v_lambda
variable q2 equal  0.06*v_lambda

fix ADAPT all adapt/fep 100000 &
  pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_lambda &
  atom charge 1 v_q1 &
  atom charge 2 v_q2 &
  after yes

変数lambdaはrampという関数を使っています。この関数は
ramp(x, y) = x + (y-x) * (timestep-startstep) / (stopstep-startstep)
というような動きをするものです。startstep=0とすると、
ramp(x, y) = x + (y-x) * (timestep) / (stopstep) で、計算のスタート時がx、終了時がyになるようにtimestepの進行に合わせて線形変化します。したがってlambdaは計算開始時に0.0、計算終了時に1.0になり、FEPにおけるターゲット分子が「ない状態」から「ある状態」に変化することを表します。

q1とq2はそれぞれ炭素の部分電荷と水素の部分電荷で、この後に続くfixコマンドで電荷を変化させる際に使用します。

fix ADAPT all adapt/fep の部分では、lambdaの値に応じて分子のpairやCoulombをコントロールすることが書かれています。
pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_lambda は、
”pairのlj/cut… の「1から2」と「3から4」のlambdaをv_lambdaに代入する”
という意味です。CH4対H2Oの相互作用のことで、先ほど 0.0 と書いてあった部分の値が変化します。（v_lambdaとはvariable lambdaのことです）

atom charge 1 v_q1 は ”原子タイプ1番の電荷をv_q1にする”という意味です。q2も同様。

compute FEP

FEP計算にはもう一つ、compute FEPというコマンドが必要です。

variable dlambda equal 0.05
variable dq1 equal -0.24*v_dlambda
variable dq2 equal  0.06*v_dlambda

compute FEP all fep ${TK} &
  pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_dlambda &
  atom charge 1 v_dq1 &
  atom charge 2 v_dq2 &
  volume yes

fix FEP all ave/time 20 4000 100000 c_FEP[*] file fep01.fep

ここでも変数を定義していますが、このdλ＝0.05の意味は次の段落で解説します。

書き方はfix ADAPT all adapt/fepに似ています。温度の変数${TK}がありますが、これはinファイルの前半のほうで300.0が指定されています。また末尾のvolume yesは、体積変化があるアンサンブルで計算しているのでyesです。

最後、fixコマンドでcompute FEPの値を fep01.fep というファイルに書き出します。

ステップ数について

FEPの制御に関連するステップ数について見ていきます。関係するのは以下の3か所。

fix ADAPT all adapt/fep 100000 

fix FEP all ave/time 20 4000 100000 c_FEP[*] file fep01.fep

run 2000000

fix ADAPT all adapt/fep 100000 では10万step毎にpairやCoulombの値を更新します。
λ＝0.0からスタートする場合、1から9万9999 stepの間はλ＝0.0です。
variable lambda equal ramp(0.0,1.0) の計算では、最終stepである200万と更新間隔10万の関係からλの値は10万step毎に0.05ずつ変化します。したがって10万から19万9999step間はλ＝0.05です。これを200万stepまで繰り返していくことでλ＝1.0まで変化します。

なお、compute FEPに使用するdλ＝0.05 はこの0.05のことです。

fix FEP all ave/time 20 4000 100000 c_FEP[*] file fep01.fep ではFEP計算の平均値を書き出します。
Nevery = 20: 20ステップごとにデータを収集。
Nrepeat = 4000: 4000回収集したデータを平均化。
Nfreq = 100000: 平均化した結果を10万ステップごとに出力。
という設定ですが、20×4000は8万で、平均値の出力間隔の10万とは一致しません。図で描くと、

|-----------------------------------------------> 時系列
|*-----------*------------*------------*-------- ラムダが変化するタイミング
|*--*--*-----*--*--*------*--*--*------*--*--*-- サンプリングするタイミング
|-----------*------------*------------*--------- 平均値を記録するタイミング

という感じです。1から8万stepのデータから平均値を記録しますが、そこから10万stepまでのデータは捨てることになります。この理由は、FEPがラムダ変化直後の応答を見ることを目的としているからです。応答、ラムダの変化直後のデータが重要なのです。ラムダが変化した直後のデータを記録し、変化からしばらく経ったデータは捨てる、というのがFEPの場合は適切なデータ取得方法というわけです。

結果の解析について

READMEに書いてあるPythonの実行コマンド例が間違っています。

READMEには
fep.py 300 < fep01.fep
と書かれていますが、
fep.py real 300 fep01.fep
だと思います。units real, 300K, ファイル名の順番です。