開発
このパッケージ nap はオープンソースプロジェクトなので，自由に修正および利用できる．自己責任での使用のこと．
通常，MD プログラムを特定の研究対象に適用する際にはソースコードを修正する必要がある．このパッケージのソースコードには多くの行が含まれているが，使用されているアルゴリズムは基本的なものが多く，分子動力学（MD）に精通していれば容易に理解できる．
それではコーディングと MD シミュレーションを楽しんで 😉
数式
ソースコードを読み，修正しようとする方のために，pmd がプログラム内でどのような処理を行っているかを詳しく示す． 基本的に pmd は速度ベルレ（velocity Verlet）アルゴリズムによる標準的な MD を実行する． したがって，まずは標準的な MD の数式を学ぶことを推奨する【Frenkel 参照】．
しかしながら，標準的な MD を実装するためには，多くのコード固有の定義が必要． このセクションでは，pmd で使用されるコード固有の定義について説明する．
シミュレーションセル
シミュレーションセルは 1 つのスカラー値 と 3 つのベクトル で構成され，原子構造ファイル内では以下のように記述される．
5.472 <--- l
0.500 0.500 0.000 <--- ax, ay, az
0.500 0.000 0.500 <--- bx, by, bz
0.000 0.500 0.500 <--- cx, cy, czコード内で実際に使用されるセル情報は行列 であり，次のように定義される．
\begin{equation} \mathbf{h} = l \times (\mathbf{a},\mathbf{b},\mathbf{c})= l \times \begin{pmatrix} a_x & b_x & c_x \\ a_y & b_y & c_y \\ a_z & b_z & c_z \end{pmatrix}. \end{equation} \end{aligned}$$ --- ### 位置座標 コード内では原子位置は `ra(1:3,1:natm)` として定義されており，MD シミュレーション中は 0.0 から 1.0 の範囲で正規化される． 絶対位置を求めるためには，以下のようにセル行列 $\mathbf{h}$ を乗じる．\mathbf{r}’ = \mathbf{h} \cdot \mathbf{r} = r_a\mathbf{a} +r_b\mathbf{b} +r_c\mathbf{c}.
この処理は Fortran では次のように記述される． ```fortran xi(1:3) = h(1:3,1)*ra(1,i) + h(1:3,2)*ra(2,i) + h(1:3,3)*ra(3,i) ``` --- ### 速度 原子の速度 `va(1:3,1:natm)` もセル行列を用いてスケーリングされる． さらに，時間によってもスケーリングされるため，コード内の速度には距離のスケールは含まれず，速度スケールのみが適用される． そのため，速度ベルレアルゴリズムにおいては，時間間隔 $\Delta t$ を掛けることなく速度を位置に直接加算できる． --- ### 加速度 原子の加速度 `aa(1:3,1:natm)` もセル行列によってスケーリングされ，さらに $\Delta t^2$ を乗じて，力計算サブルーチンの最終的な出力が位置単位となるようにする． これにより，速度ベルレループ内で $\Delta t^2$ を掛けることなく直接加速度を位置に加算できる． --- ### Velocity Verlet Velocity Verlet法によるMDプログラムはとてもシンプル： 1. 初期構造から原子にかかる力を計算する 2. Velocity Verlet ループの開始 1. 力を使って半分の時間刻み（$\Delta t/2$）だけ速度を更新． 2. 速度を使って時間刻み $\Delta t$ で座標を更新． 3. 新しい座標を使って原子にかかる力を計算． 4. 新しい力を使って半時間（$\Delta t/2$）だけ速度を更新． 数式で書くと次のようになる． $$\begin{aligned} \begin{eqnarray*} \mathbf{v}_i^{(n+1)} &=& \mathbf{v}_i^* +\frac{\mathbf{f}_i^{(n)}}{m_i} \frac{\Delta t}{2}, \\ \mathbf{r}_i^{(n+1)} &=& \mathbf{r}_i^{(n)} +\mathbf{v}_i^{(n+1)}\Delta t, \\ \text{Compute}\ &\ & \mathbf{f}_i^{(n+1)}\left(\left\{\mathbf{r}^{(n+1)}\right\}\right), \\ \mathbf{v}_i^* &=& \mathbf{v}_i^{(n+1)} +\frac{\mathbf{f}_i^{(n+1)}}{m_i} \frac{\Delta t}{2}, \end{eqnarray*} \end{aligned}$$ ただし下付き *i* は原子インデックスで，上付き *n* は MD-step． Fortranコードでは上記アルゴリズムは次のように書かれている． ```Fortran call get_force(namax,natm,tag,ra,nnmax,aa,strs,h,hi & ,tcom,nb,nbmax,lsb,lsrc,myparity,nn,sv,rc,lspr & ,mpi_md_world,myid_md,epi,epot0,nismax,acon,avol & ,cforce) ... do istp=1,nstp ... va(1:3,1:natm)=va(1:3,1:natm) +aa(1:3,1:natm) ... ra(1:3,1:natm)=ra(1:3,1:natm) +va(1:3,1:natm) ... call get_force(namax,natm,tag,ra,nnmax,aa,strs,h,hi & ,tcom,nb,nbmax,lsb,lsrc,myparity,nn,sv,rc,lspr & ,mpi_md_world,myid_md,epi,epot,nismax,acon,avol & ,cforce) ... va(1:3,1:natm)=va(1:3,1:natm) +aa(1:3,1:natm) ... enddo ``` ------------------------------------------------------------------------ ## Force implementation 原子間力およびポテンシャルの実装は，MD プログラミングの核となる部分． もし **pmd** に実装されていない元素の組み合わせを含むシミュレーションを実行したい場合は，自分で力計算ルーチンを実装する必要がある． 各力計算ルーチンは個別のモジュールとして定義されており，例えば `force_SW_Si.F90` は Si 系の Stillinger-Weber 型の力ルーチンを示する． これらのルーチンは `force_common.F` ファイル内の `get_force` サブルーチンを通じて呼び出される． そのため，`get_force` および `force_SW_Si` のサブルーチンを参考にすれば，独自の力ルーチンを作成することができる． シミュレーションで使用される **force-type**（力の種類）は，変数 `cforce` によって決定される． この変数は _in.pmd_ ファイル内で文字列として指定され，対応する力ルーチンが **force-type** の名前に基づいて呼び出される． 各力計算サブルーチンは，同じ引数と順序を持たなければならず，`get_force` ルーチンの冒頭で `use` を使用して提供される必要がある． --- ## ASE インターフェース Python スクリプトを使用して **pmd** を **ASE（Atomistic Simulation Environment）** と接続できる． これにより，*pmd* を用いた小規模な計算を簡単に実行できる． 以下のコードは `nappy/interface/ase/pmdrun.py` を ASE とともに使用する方法を示している． ```python import os, sys from ase.io import read from nappy.interface.ase.pmdrun import PMD atoms = read('POSCAR', index=0, format='vasp') os.system('cp /path/to/in.*.NN ./') calc = PMD(label='pmd', command='/path/to/nap/pmd/pmd > out.pmd', force_type='NN') atoms.set_calculator(calc) print(atoms.get_potential_energy()) print(atoms.get_forces()) ``` このスクリプトを実行すると，`atoms.get_potential_energy()` を呼び出した際に *pmd* がバックグラウンドで実行され，`out.pmd` ファイルから計算結果を取得する． ------------------------------------------------------------------------ ## Versioning and tagging **nap** は Git を使用してバージョン管理を行っている．開発者はコードを修正する際に，新しいブランチを作成し，変更をマスターブランチにマージすることを推奨する．タグの命名は `rev170605` のような形式になっており，これは _2017年6月5日_ のリビジョンを意味する．同じ日に新しいタグを作成する必要がある場合（通常は発生しませんが），タグ名に `rev170605a` のようにアルファベットを追加することで識別できる． **nap** のバージョン管理は（厳密ではありませんが）[セマンティックバージョニング](https://semver.org) に従っています．ただし，API が定義されていないため，メジャーバージョンが `1` になることはなく，バージョンは `v0.x.x` の形式になる． [^Frenkel]: Dann Frenkel and Berend Smit, *Understanding Molecular Simulation*, Academic Press