LoginSignup
6
5

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 1

Last updated at Posted at 2019-07-12

タイトルのとおりです。ソースコード（md_water_ewald.cpp）は
https://github.com/poinco-gogo/md_water
で配布しています。学生時代に作ったものなので色々とツッコミどころ満載ですが、動くことは動きます（重要）。
力場は TIP3P で、 NVE/NVT のみの対応となります。
位置 Verlet 法を採用1し、強制的に SHAKE します。2
温度制御は Langevin 方程式に従っています（BBK integrator）。3
式を explicit に書くなら

r(t+\Delta{t})=\frac{1}{B}\left\{2r(t)-Ar(t-\Delta{t})+\frac{\Delta t^2}{m}\left[F(t)+R(t)\right]\right\}\\
A=1-\gamma\frac{\Delta t}{2}\\
B=1+\gamma\frac{\Delta t}{2}\\
<R(0)R(t)>=2k_{B}T\gamma m \delta(t)

$R(t)$ はガウスランダム過程で、上式の最終行を満たします。 $T$ は熱浴の温度、 $\delta(t)$ はディラックのデルタ関数です。摩擦係数 $\gamma=0$ とすると通常の位置 Verlet の式に戻ることに注意してください。インプットの config ファイルで gamma 0とすることで NVE のシミュレーションをすることができます。

以下、コードの解説をしていきたいと思います。
（commit: 8b01725 に基づく）

冒頭

#include "vector3.h"

このプログラムでは 3 次元ベクトルの計算に
http://cms.phys.s.u-tokyo.ac.jp/~naoki/CIPINTRO/CIP/SRC/vector3.h
を使用します（ただし、 abs を vabs に変更して下さい）。今なら Eigen を使う方があとで行列演算をしたくなった時に便利でしょう。

クラス定義（Atom）

class Atom
{       
        public: 
        Vector3 position;
        Vector3 velocity;
        Vector3 fold, fnew, vnew, rold, rnew;
        string  atomname;
        double  charge, mass, R; 
        Atom(double x, double y, double z)
        {       
                position.x = x;
                position.y = y;
                position.z = z;
        }
};

原子を記述するクラスです。R は温度制御の時に使用する定数です（下記参照）。

クラス定義（System）

class System
{
        public:
        int nstep, natom, kmax;
        double dt, T, gamma, boxsize, cutoff, volume, ewcoeff, lj, es;
        Vector3 origin, a1, a2, a3, g1, g2, g3, b1, b2, b3;

        System()
        {
                dt = 9999;
                T = 9999;
                gamma = 9999;
                nstep = 0;
                natom = 0;
                boxsize = 1;
                origin.x = 9999; origin.y = 9999; origin.z = 9999;
                cutoff = 1.;
                lj = 0;
                es = 0;
        }

        void setup_box()
        {
                volume = boxsize * boxsize * boxsize;
                a1.x = boxsize; a1.y = 0; a1.z = 0;
                a2.x = 0; a2.y = boxsize; a2.z = 0;
                a3.x = 0; a3.y = 0; a3.z = boxsize;
                double rvol2pi = 2. * M_PI / volume;
                g1 = (a2 % a3) * rvol2pi;
                g2 = (a3 % a1) * rvol2pi;
                g3 = (a1 % a2) * rvol2pi;
        }
};

系の情報を扱うクラスです。gamma は Langevin 方程式の friction coefficient というやつです。

main 関数

if (argc < 4)
{
        cout << "usage: ./a.out trjout psf pdb config\n";
        return 1;
}

引数が足りないときに使い方を表示します（return 1 でいいのか？）。

引数 意味
trjout トラジェクトリを保存するアウトプットファイル名。トラジェクトリは XYZ 形式で保存します。
psf 系に存在するすべての分子の結合（bond, angle,...）の情報を記述するファイル（psf 形式）
pdb 初期座標（PDB 形式）
config MD のパラメータを記述するファイル
Github に、 psf, pdb, config ファイルのサンプルがあります。 
srand( (unsigned int)time(NULL) );

現在時刻をもとに乱数の発生系列を定義します。あとで rand() によって乱数を発生させていますが、 C++11 ではメルセンヌ・ツイスターが標準ライブラリに実装されているので、そちらを使う方がよいでしょう（参考）。

const int nfree = natom * 3 - nwat * 3;

系の自由度を定義しています。強制的に SHAKE しているので、水の結合の自由度を全体から引いておきます。

const double ps2asu = 20.455;

本プログラムでは原子の質量を原子質量単位（atomic mass unit）、エネルギーの単位を kcal/mol、長さの単位をオングストロームで表現しています。この場合、時間の単位はピコ秒に ps2asu をかけたものになります。

const double dt_ps = system.dt * 0.001; // in ps

インプットの config ファイルに dt 2 と書いた場合、それはタイムステップが 2 フェムト秒であることを意味します。なのでいったんピコ秒に変換します。

const double A = 1. - gamma * dt * 0.5;
const double B = 1. + gamma * dt * 0.5;
const double inB = 1. / B;

冒頭の式の $A$ と $B$ を計算しています。

for (int i = 0; i < natom; i++)
{      
        Atom& at = atomVector[i];
        at.R = sqrt(2. * kb * T * gamma * at.mass / dt);
}

各原子ごとに上式の $R$ の平均的な大きさを計算します。これに正規分布 $N(0, 1)$ に従う乱数を乗じたものが、ランダム力になります。

const double rOH  = 0.9572;
const double aHOH = 104.52 / 180. * acos(-1.0);

TIP3P のジオメトリを定義しています。

// make reciprocal vectors
vector<Vector3> g;
int kmax = system.kmax;
int sqkmax = kmax - 1;
sqkmax *= sqkmax;
double dum = 0;
for (int k = 0; k < kmax; k++)
for (int i = -kmax + 1; i < kmax; i++)
for (int j = -kmax + 1; j < kmax; j++)
{
        if (k * k + j * j + i * i > sqkmax or (i==0 && j==0 && k==0))
                continue;
        Vector3 vtmp(i * system.g1.x, j * system.g2.y, k * system.g3.z);
        if (vabs(vtmp) > dum)
                dum = vabs(vtmp);
        g.push_back(vtmp);
}

逆格子ベクトル $\vec{G}=2 \pi \left(\frac{k_{x}}{L_{x}},\frac{k_{y}}{L_{y}},\frac{k_{z}}{L_{z}} \right)$ を定義しています。 $k_{x}$, $k_{y}$ については、 config ファイルで指定した kmax の値をもとに -kmax から +kmax までとりますが、 $k_{z}$ については0 から +kmax までしかとりません。この条件の下で結果を 2 倍すれば求めるべき値を得ることができます。4
kmax は、静電相互作用エネルギーが収束するまで充分大きくとりましょう。
（エワルド係数 ewcoeff との兼ね合いで設定するべき kmax の値が変わることに注意してください。）

double ew_self = 0;
for (int i = 0; i < atomVector.size(); i++)
{
        ew_self += atomVector[i].charge * atomVector[i].charge;
}
ew_self *= -system.ewcoeff / sqrt(M_PI) * 332.0636;

エワルドの自己相互作用エネルギー項です。332.0636 は、$e^2/4 \pi \epsilon_{0}$ だったような気がしますが、記憶が定かではありません。

for (int i = 0; i < natom; i++)
{
        Atom& at = atomVector[i];
        at.velocity.x = gauss() * sqrt(kb * T / at.mass);
        at.velocity.y = gauss() * sqrt(kb * T / at.mass);
        at.velocity.z = gauss() * sqrt(kb * T / at.mass);
}

初速度をボルツマン分布からサンプリングしています。

for (int i = 0; i < atomVector.size(); i++)
{
        Atom& at = atomVector[i];
        at.fold = at.fnew;
        at.rold = at.position;
        at.position = at.rold + dt * at.velocity + dtdt_div2 / at.mass * at.fold;
}

冒頭の式で、 $r(t+\Delta t)$ を計算するために $r(t)$ と $r(t-\Delta t)$ が必要です。ここではインプットの座標が $r(t-\Delta t)$ だと思うことにして、 $r(t)$ をテイラー展開の 2 次までとることで求めます。

r(t) = r(0)+v(0)\Delta t + \frac{\Delta t^2}{2m}F(0)

for (int istep = 1; istep <= nstep; istep++)
{

この中が MD のメインループです。

calc_frc(atomVector, lj_pair_list, el_pair_list, system, g);

まず position に保存された座標に基づいて force を計算します。計算された force は fnew に保存されます。

for (int i = 0; i < atomVector.size(); i++)
{
        Vector3 noise( gauss(), gauss(), gauss());
        Atom& at = atomVector[i];
        at.rnew = 2. * at.position - at.rold + gamma * dt_div2 * at.rold + dt * dt / at.mass * (at.fnew + at.R * noise);
        at.rnew = at.rnew * inB;
}

冒頭の式に基づいて $r(t+\Delta t)$ を計算します（ rnew に保存される）。

for (int ishake = 0; ishake < 100; ishake++)
{
        if (shake(atomVector, shake_list))
                break;
        if (ishake == 99)
        {
                cerr << "error: shake does not converged.\n";
                return 1;
        }
}

SHAKE ループです。

for (int i = 0; i < atomVector.size(); i++)
{
        Atom& at = atomVector[i];
        at.vnew = 0.5 / dt * (at.rnew - at.rold);
}

$v(t+\Delta t)=\frac{r(t+\Delta t)-r(t-\Delta t)}{2\Delta t}$ で、更新された速度ベクトルを計算します。結果は vnew に保存されます。

double K = calc_kin(atomVector);
calc_pot(atomVector, lj_pair_list, el_pair_list, system, g);
system.es += ew_self;

運動エネルギー、ポテンシャルエネルギーを計算します。
運動エネルギーは vnew に基づいて計算しています。
ポテンシャルエネルギーは position に基づいて計算しているみたいです。これは変ですね。 rnew に基づく計算であるべきではないだろうか？

for (int i = 0; i < atomVector.size(); i++)
{
        Atom& at = atomVector[i];
        at.rold = at.position;
        at.position = at.rnew;
        at.velocity = at.vnew;
}

次のステップに備えて変数の値を入れ替えています。

以上で MD のループは終わりです。
ちょっと長くなってきたので、続きは別にします。
続き→水（TIP3P）の古典 MD プログラム part 2

  1. 拘束の導入が簡単なので学習目的にはよいでしょう。しかし速度 Verlet の方が velocity のエラーが少ないので良いかもしれません。有名な MD プログラムは大抵速度 Verlet か leap frog なのではないでしょうか。 ↩

  2. 時間の刻み幅を大きくしたいからですね。 ↩

  3. BBK integrator の速度 velocity バージョンは Bussi, G.; Parrinello, M. "Accurate sampling using Langevin dynamics", Phys. Rev. E, 2007, 75, 056707. で提案された integrator の、時間の刻み幅が小さい極限で一致する、気がする。 ↩

  4. Tuckerman, M. E. "Statistical Mechanics: Theory and Molecular Simulation", Oxford Graduate Texts, pp. 661 ↩

6
5
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
6
5

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?