LoginSignup
2
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@poinco_gogo

水（TIP3P）の古典 MD プログラム part 3

Last updated at Posted at 2019-07-24

これは https://github.com/poinco-gogo/md_water にアップロードした、水（TIP3P）の古典 MD プログラム（md_water_ewald.cpp）の解説です（commit: 8b01725 に基づく）。
水（TIP3P）の古典 MD プログラム part 1
では integrator を含む main 関数を、
水（TIP3P）の古典 MD プログラム part 2
では主に SHAKE （拘束）関数を紹介しました。

本稿ではいよいよエネルギーと力の計算について解説します。

なお、本プログラムでは、なぜかエネルギーと力の計算の関数を分離しています（そういう仕様にしたかった動機は忘れてしまいました）。よくないことに、エネルギーと力で似たような計算をしているのに、共通の関数を作らないで済ましてしまっています。「動けばいいや」の精神ですね。メンテナンスが大変ですが。

TIP3P 力場

TIP3P 力場1は以下の関数形、

E_{mn}=\frac{e^2}{4\pi\epsilon_0}\sum_{i\in m}\sum_{j\in n}\frac{q_i q_j}{r_{ij}}+\frac{A}{r_{OO}^{12}}-\frac{B}{r_{OO}^{6}}

とパラメータ（$q_i$, $A$, $B$）で定義されます。ただし $E_{mn}$ はモノマー m と n の間の相互作用エネルギーです（本プログラムでは強制的に SHAKE しているので、モノマーの内部自由度に由来するエネルギーは発生しないことに注意してください）。第 1 項が点電荷間の静電相互作用を、第 2 項が酸素原子間の LJ 相互作用を表現しています。本来の TIP3P モデルでは、水素原子には LJ パラメータが存在しないことに注意してください。これに対して、 CHARMM 力場で使用される TIP3P モデルでは水素原子に LJ パラメータを割り当てており、混乱の元になっている感があります。

論文によるとパラメータの値は、

パラメータ 値
$A$ (kcal Å^12/mol) 582.0 x 1e-3
$B$ (kcal Å^6 /mol) 595.0
$q(O)$ -0.834
$q(H)$ 0.417

$A$ と $B$ の値はエネルギーと力の計算を行う関数（ calc_pot, calc_frc ）の冒頭にハードコーディングしてあります（）。電荷の値は psf ファイルから読んでいます（ load_psf ）。

本プログラムでは、 config ファイルで指定したカットオフ距離（ cutoff ）を、近距離相互作用と長距離相互作用の境界として採用します。デフォルトではこれは 8 Å となっています。

LJ 相互作用は主に近距離で作用するため、本プログラムではカットオフ距離以遠の相互作用を無視します。本プログラムは強制的に周期境界条件を課しているので、 [nearest image convention] (https://web.northeastern.edu/afeiguin/p4840/p131spring04/node41.html) で相互作用を計算します。この条件の下では、シミュレーションボックスの一辺のサイズを $L$ とするとき、粒子間距離の最大値が $L/2$ になることに注意してください。従って、 $ L $ とカットオフ $r_{cut}$ の関係、

\frac{L}{2}>r_{cut}

が必要となります。デフォルトの config ファイルでは、ボックスサイズを 16.9461011959 Å、カットオフを 8 Å にしています。

これに対して静電相互作用は長距離成分が無視できないため、カットオフ以遠の相互作用も含める必要があります。周期境界条件下では、 $N$ 個の点電荷 {$q_i$} の静電相互作用の和は

E^{elect}=\frac{e^2}{4\pi\epsilon_0}\frac{1}{2}\sum_{\vec{n}=0}\sum_{i=1}^N{\sum_{j=1}^N}^{\prime}\frac{q_i q_j}{|\vec{r}_i-\vec{r}_j+\vec{n}L|}

と書けます。ただし $\vec{n}=(n_x,n_y,n_z)$ は整数成分のベクトル、 $L$ はシミュレーションボックス（orthorhombic で辺の長さは全て等しいとする）の一辺の長さ、和の $^\prime$ は $\vec{n}=0$ の時に $i=j$ を除くことを意味します。

この無限和を計算するために、通常 Ewald の方法が採用されます。

Ewald 和

以下に Ewald 和のエッセンスを記述します。詳細な導出は他書に譲ります。2$^,$3

$\vec{r}_i$ に存在する点電荷 $q_i$ が任意の場所 $\vec{r}$ に作る静電ポテンシャルは

\phi_i(\vec{r})=\frac{e}{4\pi\epsilon_0}\frac{q_i}{|\vec{r}-\vec{r}_i|}

周期境界条件下で全ての電荷を考慮すると、

\phi(\vec{r})=\frac{e}{4\pi\epsilon_0}\sum_{\vec{n}=0}\sum_{j=1}^N\frac{q_j}{|\vec{r}-\vec{r}_j+\vec{n}L|}

ここで $i$ 番目の粒子の寄与を除いたポテンシャルを考えると2

\phi_{[i]}(\vec{r})\equiv\phi(\vec{r})-\phi_i(\vec{r})=\frac{e}{4\pi\epsilon_0}\sum_{\vec{n}=0}{\sum_{j=1}^N}^{\prime}\frac{q_j}{|\vec{r}-\vec{r}_j+\vec{n}L|}

すると上の $E^{elect}$ は

E^{elect} = \frac{1}{2}\sum_{i=1}^Neq_i\phi_{[i]}(\vec{r}_i)

と書けます。

さて、場所 r_i に存在する点電荷 $q_i$ の電荷密度は

\rho_i(\vec{r})=q_i \delta(\vec{r}-\vec{r}_i)

と書けます。ここで $\delta(\cdot)$ はディラックのデルタ関数です。実際、一般の電荷密度分布に対する $\phi_{[i]}(\vec{r})$ の表式、

\phi_{[i]}(\vec{r})=\frac{e}{4\pi\epsilon_0}\sum_{\vec{n}=0}{\sum_{j=1}^N}^{\prime}\int\frac{\rho_j(\vec{r}^{\prime})}{|\vec{r}-\vec{r}^{\prime}+\vec{n}L|}\text{d}\vec{r}^{\prime}

にこのデルタ関数を代入すると、点電荷密度に対する $\phi_{[i]}(\vec{r})$ の表式が復元されます。

さてここで以下のようなトリックを導入します。
まず、粒子 $i$ の元の電荷密度分布に、反対電荷のガウス電荷密度分布を足します。すると合成後の電荷分布は、

\rho_i^S(\vec{r})=q_i\delta(\vec{r}-\vec{r}_i)-q_iG_\alpha(\vec{r}-\vec{r}_i)\\
G_\alpha(\vec{r})=\left( \frac{\pi}{\alpha^2} \right)^{-3/2}\text{exp}\left(-\alpha^2|\vec{r}|^2\right)

となります。反対電荷の分布を足すことで、元の電荷の影響を近距離に抑えることが期待できます。
電荷分布と静電ポテンシャルの関係は、ポアソン方程式

\nabla^2\phi_i^S(\vec{r})=-\frac{\rho_i^S(\vec{r})}{\epsilon_0}

で記述されます。合成後の電荷分布に対するポアソン方程式を解析的に解くことができ、結果は

\phi_i^S(\vec{r})=\frac{e}{4\pi\epsilon_0}\frac{q_i}{|\vec{r}-\vec{r}_i|}\left[ 1-\text{erf}\left( \alpha|\vec{r}-\vec{r}_i| \right)\right]\\
=\frac{e}{4\pi\epsilon_0}\frac{q_i}{|\vec{r}-\vec{r}_i|}\text{erfc}\left( \alpha|\vec{r}-\vec{r}_i| \right)

となります。第 1 行は、ガウス電荷密度分布に由来する静電ポテンシャルが誤差関数 $\text{erf}(\cdot)$ になることを示しています。第 2 行の $\text{erfc}(\cdot)$ は相補誤差関数で、パラメータ $\alpha$ （Ewald 係数）をうまく調整することで $\phi_i^S(\vec{r})$ の影響をカットオフの範囲内に抑えることができるというわけです。本プログラムのデフォルトのカットオフは 8 Å で、 Ewald 係数の値は 0.39467 1/Å です（どちらも config ファイルで指定できます）。
合成された電荷密度分布が作るポテンシャルは periodic image も含めると、

\phi_{[i]}^S(\vec{r})=\frac{e}{4\pi\epsilon_0}\sum_{\vec{n}=0}{\sum_{j=1}^{N}}^{\prime}\frac{q_j}{|\vec{r}-\vec{r}_j+\vec{n}L|}\text{erfc}\left( \alpha|\vec{r}-\vec{r}_j+\vec{n}L| \right)

ここまでで、近距離に抑えられたポテンシャルを得ることに成功しました。現実のポテンシャルは、 $\phi^S(\vec{r})$ からガウス電荷密度分布の影響を取り除いたものになります。
上で書いたように、ガウス電荷密度分布に由来するポテンシャルは誤差関数を用いて

\phi_i^L(\vec{r})=\frac{e}{4\pi\epsilon_0}\frac{q_i}{|\vec{r}-\vec{r}_i|}\text{erf}\left( \alpha|\vec{r}-\vec{r}_i| \right)

と書けます。周期境界条件下の全てのイオンを含めると、

\phi_{[i]}^L(\vec{r})=\frac{e}{4\pi\epsilon_0}\sum_{\vec{n}=0}{\sum_{j=1}^{N}}^{\prime}\frac{q_j}{|\vec{r}-\vec{r}_j+\vec{n}L|}\text{erf}\left( \alpha|\vec{r}-\vec{r}_j+\vec{n}L| \right)

従ってエネルギーは

E^{elect}=\frac{1}{2}\sum_{i=1}^Neq_i\phi_{[i]}^S(\vec{r}_i)+\frac{1}{2}\sum_{i=1}^Neq_i\phi_{[i]}^L(\vec{r}_i)\\
=\frac{1}{2}\sum_{i=1}^Neq_i\phi_{[i]}^S(\vec{r}_i)+\frac{1}{2}\sum_{i=1}^Neq_i\phi^L(\vec{r}_i)-\frac{1}{2}\sum_{i=1}^Neq_i\phi_i^L(\vec{r}_i)\\
=E^S+E^L-E^{self}

となります。ただし

E^S\equiv\frac{1}{2}\sum_{i=1}^Neq_i\phi_{[i]}^S(\vec{r}_i)

は実空間で、

E^L\equiv\frac{1}{2}\sum_{i=1}^Neq_i\phi^L(\vec{r}_i)

はあとで見るように逆空間で和を計算します。

E^{self}\equiv\frac{1}{2}\sum_{i=1}^Neq_i\phi_i^L(\vec{r}_i)

は r_i を中心とするガウス電荷密度分布と r_i に存在する点電荷 $q_i$ の相互作用を表現しており、その値は極限値として

E^{self}=\frac{e^2}{4\pi\epsilon_0}\frac{\alpha}{\sqrt{\pi}}\sum_{i=1}^N q_i^2

と計算できます（自己相互作用項と呼ばれる）。2$^,$3これは粒子の位置に依存しないので、動力学計算の前に 1 度だけ計算しておけば良いですね。

上で述べたように、 $\phi^L$ に含まれる $\text{erf}(\cdot)$ は近距離的ではないので実空間では計算できませんが、逆空間（reciprocal space）における和に変換すれば収束が期待できます。結果を書くと

\phi^L(\vec{r})=\frac{e}{V\epsilon_0}\sum_{\vec{g}\neq 0}\frac{\text{exp}(-g^2/4\alpha^2)}{g^2}\sum_{j=1}^N q_j \text{exp}\left[-i\vec{g}\cdot(\vec{r}-\vec{r}_j)\right]

となります。ただし $V$ はシミュレーションボックスの体積で、 $k_x, k_y, k_z$ を整数として $\vec{g}=2\pi\left(k_x/L,k_y/L,k_z/L\right)$ は逆格子ベクトルです。$|\vec{g}|$ を十分大きくとればこの和は収束します。$E^L$ は、

E^L=\frac{e^2}{4\pi\epsilon_0}\frac{2\pi}{V}\sum_{\vec{g}\neq 0}\frac{\text{exp}(-g^2/4\alpha^2)}{g^2}\sum_{i=1}^N\sum_{j=1}^N q_iq_j \text{exp}\left[-i\vec{g}\cdot(\vec{r}_i-\vec{r}_j)\right]

となります。ここで $i$ と $j$ の和にプライム（$\prime$）がついていないことに注意してください。実空間では、 $\vec{L}=0$ のときに $i=j$ とすることで自己相互作用を除くことができました。逆空間では、一度全部足して閉まってから、あとで補正する（$E^{self}$）のが簡単なようです。3

ポテンシャル・エネルギーの計算

calc_pot 関数の解説です。

for (int i = 0; i < lj_pair_list.size(); i++)
{
        Atom& iat = atomVector[lj_pair_list[i]];
        for (int j = i + 1; j < lj_pair_list.size(); j++)
        {
                Atom& jat = atomVector[lj_pair_list[j]];
                Vector3 del = iat.position - jat.position;
                del.x -= boxsize * floor(del.x / boxsize + 0.5);
                del.y -= boxsize * floor(del.y / boxsize + 0.5);
                del.z -= boxsize * floor(del.z / boxsize + 0.5);
                double roo2 = norm(del);
                if (roo2 > cutoff2) continue;
                double roo6 = roo2 * roo2 * roo2;
                system.lj += A / (roo6 * roo6) - B / roo6;
        }
}

LJ 相互作用を、 nearest image convention で計算しています。

// ewald intra energy
for (int i = 0; i < atomVector.size() / 3; i++)
{
        int j = 3 * i;
        Atom& at1 = atomVector[j];
        Atom& at2 = atomVector[j + 1];
        Atom& at3 = atomVector[j + 2];
        Vector3 del1 = at1.position - at2.position;
        del1.x -= boxsize * floor(del1.x / boxsize + 0.5);
        del1.y -= boxsize * floor(del1.y / boxsize + 0.5);
        del1.z -= boxsize * floor(del1.z / boxsize + 0.5);
        double adel1 = vabs(del1);
        Vector3 del2 = at1.position - at3.position;
        del2.x -= boxsize * floor(del2.x / boxsize + 0.5);
        del2.y -= boxsize * floor(del2.y / boxsize + 0.5);
        del2.z -= boxsize * floor(del2.z / boxsize + 0.5);
        double adel2 = vabs(del2);
        Vector3 del3 = at3.position - at2.position;
        del3.x -= boxsize * floor(del3.x / boxsize + 0.5);
        del3.y -= boxsize * floor(del3.y / boxsize + 0.5);
        del3.z -= boxsize * floor(del3.z / boxsize + 0.5);
        double adel3 = vabs(del3);
        ew_intra += at1.charge*at2.charge*erfl(ewcoeff*adel1)/adel1;
        ew_intra += at1.charge*at3.charge*erfl(ewcoeff*adel2)/adel2;
        ew_intra += at3.charge*at2.charge*erfl(ewcoeff*adel3)/adel3;
}

冒頭で述べたとおり、TIP3P 力場では分子内の静電相互作用は計算してはいけません。この下のコードで、実空間の和ではこのルールが守られているのですが、逆空間の和は分子内静電相互作用を含んでしまっています。そこで、逆空間の和の元の表式、

E^{intra}=\frac{e^2}{4\pi\epsilon_0}\sum_{(i,j)\in I}\frac{q_jq_i}{|\vec{r}_j-\vec{r}_i|}\text{erf}\left( \alpha|\vec{r}_j-\vec{r}_i| \right)

に従って分子内静電相互作用を計算し、エネルギーの合計値から差し引くことにします。ただし和は分子内の粒子の全てのペアについてとります。

// ewald direct space summation
for (int i = 0; i < el_pair_list.size() / 2; i++)
{
        Atom& at1 = atomVector[el_pair_list[2 * i]];
        Atom& at2 = atomVector[el_pair_list[2 * i + 1]];

        Vector3 del = at1.position - at2.position;
        del.x -= boxsize * floor(del.x / boxsize + 0.5);
        del.y -= boxsize * floor(del.y / boxsize + 0.5);
        del.z -= boxsize * floor(del.z / boxsize + 0.5);
        double r = vabs(del);
        if (r > cutoff) continue;
        ew_direct += at1.charge * at2.charge * erfc(ewcoeff*r)/r;
}

実空間の和です。 el_pair_list は分子間のペアを保存した vector です。実空間も nearest image convention で計算していることに注意してください。

// ewald reciprocal space summation
for (int i = 0; i < g.size(); i++)
{
        double ag = vabs(g[i]);
        double agag = ag * ag;
        double re, im, dtmp;
        re = 0;
        im = 0;
        for (int j = 0; j < atomVector.size(); j++)
        {
                Vector3 del = atomVector[j].position;
                del.x -= boxsize * floor(del.x / boxsize + 0.5);
                del.y -= boxsize * floor(del.y / boxsize + 0.5);
                del.z -= boxsize * floor(del.z / boxsize + 0.5);
                double dot = g[i] * del;

                re += atomVector[j].charge * cos(dot);
                im += atomVector[j].charge * sin(dot);
        }
        g[i].z ? dtmp = 1. : dtmp = 0.5;
        ew_recipro += dtmp * exp(-agag * factor) / agag * (re * re + im * im);
}

ew_recipro *= const_recipro;

逆空間の和です。g[i].z ? dtmp = 1. : dtmp = 0.5; によって、 z 成分が 0 である逆格子ベクトルの計算結果を半減します。 part 1 における逆格子ベクトルの取り方も参照してください。

system.es += ew_direct + ew_recipro - ew_intra;
system.es *= C;

静電相互作用エネルギーの合計値を保存します。 C は $\frac{1}{4\pi\epsilon_0}$ だったような気がします。

力の計算

頑張って微分の計算を行います。
まず

\frac{\partial}{\partial \vec{r}_i}\text{erfc}\left( \alpha r_{ij} \right)
=-\frac{2\alpha}{\sqrt{\pi}}\frac{\text{exp}\left( -\alpha^2 r_{ij}^2 \right)}{r_{ij}}\vec{r}_{ij}

ただし $r_{ij} = |\vec{r}_{ij}|$ で

\vec{r}_{ij}=\vec{r}_i - \vec{r}_j

従ってある電荷ペア $q_i$, $q_j$ について,

U_{ij}^{direct}=q_iq_j\frac{\text{erfc}\left( \alpha r_{ij} \right)}{r_{ij}}

ならば

\frac{\partial U_{ij}^{direct}}{\partial \vec{r}_i}
=-q_iq_j\left[
\frac{\text{erfc}\left( \alpha r_{ij} \right)}{r_{ij}}+\frac{2\alpha}{\sqrt{\pi}}\text{exp}\left( -\alpha^2 r_{ij}^2 \right)
\right]
\frac{\vec{r}_{ij}}{r_{ij}^2}

かつ

\frac{\partial U_{ij}^{direct}}{\partial \vec{r}_j}=-\frac{\partial U_{ij}^{direct}}{\partial \vec{r}_i}

さらに、ある電荷ペアについて、

U_{ij}^{intra}=q_iq_j\frac{\text{erf}\left( \alpha r_{ij} \right)}{r_{ij}}

ならば

\frac{\partial U_{ij}^{intra}}{\partial \vec{r}_i}
=-q_iq_j\left[
\frac{\text{erf}\left( \alpha r_{ij} \right)}{r_{ij}}-\frac{2\alpha}{\sqrt{\pi}}\text{exp}\left( -\alpha^2 r_{ij}^2 \right)
\right]
\frac{\vec{r}_{ij}}{r_{ij}^2}

かつ

\frac{\partial U_{ij}^{intra}}{\partial \vec{r}_j}=-\frac{\partial U_{ij}^{intra}}{\partial \vec{r}_i}

逆空間の和について、構造因子

\sum_{i=1}^N\sum_{j=1}^N q_iq_j \text{exp}\left[-i\vec{g}\cdot(\vec{r}_i-\vec{r}_j)\right]=|S(\vec{g})|^2\\
S(\vec{g})\equiv\sum_{j=1}^N q_j\text{exp}\left(i\vec{g}\cdot\vec{r}_j\right)

の微分が必要になります。
まず

\frac{\partial S(\vec{g})}{\partial \vec{r}_j}=iq_j\text{exp}\left(i\vec{g}\cdot\vec{r}_j\right)

従って、

\frac{\partial |S(\vec{g})|^2}{\partial \vec{r}_j}=S^*(\vec{g})\frac{\partial S(\vec{g})}{\partial \vec{r}_j}+S(\vec{g})\frac{\partial S^*(\vec{g})}{\partial \vec{r}_j}\\
=S^*(\vec{g})\frac{\partial S(\vec{g})}{\partial \vec{r}_j}+c.c.\\
=2\text{Re}\left[ S^*(\vec{g})\frac{\partial S(\vec{g})}{\partial \vec{r}_j}\right]\\
=2\sum_{i=1}^Nq_iq_j\vec{g}\text{sin}\left( \vec{g}\cdot\vec{r}_{ij} \right)

後の計算は、基本的にポテンシャル・エネルギーのときと同様です。

まとめ

TIP3P （剛体）モデルの分子動力学プログラムを作成しました。
改善点は、

  • 内部自由度（結合・結合角）の運動も考慮する
  • 速度 Verlet 法にする
  • 圧力制御を可能にする
  • 一般のタンパク質も計算可能にする
  • PME を導入する
  • 並列化する

などが挙げられます。このプログラムは Github で配布しているので自由に改造して遊んでみてください。

続き→水（TIP3P）の古典 MD プログラム 番外編

  1. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. (1983) "Comparison of simple potential functions for simulating liquid water", J. Chem. Phys., 79, 926-935. ↩

  2. Lee, H.; Cai, W. (2009) "Ewald summation for Coulomb interactions in a periodic supercell" ↩ ↩2 ↩3

  3. Smith, W. (2014) "Elements of molecular dynamics" (PDF file) ↩ ↩2 ↩3

2
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
2
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?