原子ダイナミクス

振動解析

機能の概要

PHASEには格子振動の基準モードを計算する振動解析機能があります。まず、原子を平衡位置からわずかに変位させて力計算を行います。その力から力定数行列を計算し、それから動力学行列を計算します。動力学行列の固有値問題を解くことにより,基準振動の振動数と固有ベクトルを計算します。

\(i\)番目の原子の安定位置から変位を\(\mathbf{u}_{i}\)とします。 変位が微小で二次以上の項が無視できるとき,格子系の運動方程式は

()\[m_{i}{\ddot{u}}_{\text{iα}} = - \sum_{\text{jβ}}^{}\Phi_{i\alpha,j\beta}u_{\text{jβ}}\]

と書けます。\(\Phi_{i\alpha,j\beta}\)は力の定数で,原子変位に関する系のエネルギー\(E(\mathbf{u}_{1},\mathbf{u}_{2},\ldots)\)の二階微分として定義されています。

()\[\Phi_{i\alpha,j\beta} = \frac{\partial^{2}E}{\partial u_{\text{iα}}\partial u_{\text{jβ}}}\]

力の定数はヘルマン-ファインマン力を原子変位で微分することにより, 求めることができます。

()\[\Phi_{i\alpha,j\beta} = - \frac{\partial F_{\text{iα}}}{\partial u_{\text{jβ}}}\]

本プログラムでは,この微分は中央差分近似で行われます。変位パラメターを\(a\)とすると

()\[\frac{\partial F_{\text{iα}}}{\partial u_{\text{jβ}}} = \frac{F_{\text{iα}}|_{u_{\text{jβ}} = a} - F_{\text{iα}}|_{u_{\text{jβ}} = - a}}{2a}\]

と書けます。 力の定数には結晶の対称性による制約があり,これを満たすように力の定数を補正する必要があります。第一に,\(i,j\)原子が空間群の対称操作\(\{ R|\mathbf{T}\}\)\(i^{'},j^{'}\)原子に移るとき,力定数テンソル\(\Phi_{i,j}\)は力定数テンソル\(\Phi_{i^{'},j^{'}}\)を回転行列Rで回転させたものに等しいです。つまり,

\[R\mathbf{r}_{i} + \mathbf{T} = \mathbf{r}_{i^{'}}\]
()\[R\mathbf{r}_{j} + \mathbf{T} = \mathbf{r}_{j^{'}}\]

ならば,

()\[\Phi_{i,j} = R^{T}\Phi_{i^{'},j^{'}}R\]

でなければなりません。 第二に,力定数テンソル\(\Phi_{i,j}\)の成分\(\text{αβ}\)を原子の番号\(j\)すべてにわたり足し合わせると,ゼロになります。つまり,

()\[\sum_{j}^{}\Phi_{i\alpha,j\beta} = 0\]

です。 第三に,力定数行列は対称でなければなりません。つまり,

()\[\Phi_{i\alpha,j\beta} = \Phi_{j\beta,i\alpha}\]

です。 換算変位\(w_{\text{iα}} = u_{\text{iα}}\sqrt{m_{i}}\)と 動力学行列\(D_{i\alpha,j\beta} = \Phi_{i\alpha,j\beta}/\sqrt{m_{i}m_{j}}\)を用いて, 格子系の運動方程式 (5.25)

()\[{\ddot{w}}_{\text{iα}} = - \sum_{\text{jβ}}^{}D_{i\alpha,j\beta}w_{\text{jβ}}\]

と書き換えます。この方程式を解くために, \(w_{\text{iα}} = Q\xi_{\text{iα}}e^{i\omega t + \delta}\) という解を仮定します。

()\[\omega^{2}\xi_{\text{iα}} = \sum_{\text{jβ}}^{}D_{i\alpha,j\beta}\xi_{\text{jβ}}\]

これは固有値が\(\omega^{2}\)で,固有ベクトルが\(\xi_{\text{iα}}\)となる, 行列\(D_{i\alpha,j\beta}\)の固有値問題です。振動解析機能ではこの固有値問題を解き,格子振動の基準モードを求めます。

入力パラメータ

振動解析を行うには,まず原子が平衡位置にあることが必要です。平衡状態にないと動力学行列の固有値が負になり,振動数が純虚数のソフトモードが現れます。平衡位置の原子座標は,構造最適化機能を用いて計算します。構造最適化計算が終了したら,nfdynm.dataの最後に記述されている最適構造での入力パラメータファイルを作成します。

振動解析の設定は、Phononブロックで指定します。

Phonon{
   sw_phonon = on
   sw_calc_force = on
   sw_vibrational_modes = on
   displacement = 0.05
}

振動解析の入力変数を以下に示します。

振動解析に関係する変数の説明

変数名またはタグ名

デフォルト値

説明

sw_phonon

OFF

格子振動解析 設定ブロックを有効にする かどうかのスイッチです。

sw_calc_force

OFF

振動 解析のための力計算を行う かどうかのスイッチです。

ON:格子振動 解析のための力計算を行い ます。(計算した力はforc e.dataに出力されます。)

OFF:sw_vibrational_modes =ONなら ファイル"F_FORCE"から力 のデータを読み込みます。

displacement

0.1

原子変位パラメーター。

sw_vibrational_modes

OFF

格子振動解析を行う かどうかのスイッチです。

ON:格子振動解析が 行われ,mode.dataファイ ルに結果が出力されます。

OFF:格子 振動解析は行われません。

norder

1

差分次数を 変更するパラメターです。

sw_polynomial_fit

OFF

ON:多項式フィッ トで力の微分を求めます。

OFF:差 分で力の微分を求めます。

  • 原子座標と対称性の入力

原子座標は反転対称があってもすべて入力する必要があります。すなわち,原子座標のweight属性値による省略は利用できません。また,sw_inversionはOFFとする必要があります。振動モードの分類と入力座標の対称性チェックに系の空間群を使用するので,結晶構造またはその空間群をsymmetryブロックで正しく指定します。ただし,対称性自動判定機能を利用することも可能です。

  • 元素の質量の指定

元素の質量はelement_listブロックの変数massで指定する。原子単位(a.u.)ではなく,原子質量単位(amu)で入力するには,#units atomic_massを#tag行の上に挿入する。

  • 原子変位の選択

原子変位はPhononブロックのdisplacementで設定します。通常,原子変位は0.1 a.u.以下にとると良いです。振動数の原子変位依存性を調べて,希望する振動数の収束が得られる原子変位に設定します。norderを2に設定することで,差分の次数を3から5に換えることができます。diplacementで設定した値をuとすれば,原子変位は-u,-u/2,u/2,uになります。sw_polynamial_fitをONにして多項式フィットにすれば,norderを2より大きく設定できます。そのときの原子変位は-u/norder,-u/(norder-1),...,u/(norder-1),u/norderです。norderを大きくすると微分精度はよくなりますが,力計算の回数が増えるので,計算時間はnorderが1の場合の2*norder倍になるので注意してください。

計算結果の出力

振動解析結果は、振動解析結果ファイルmode.data、力のデータforce.dataに出力されます。

mode.dataには振動解析の結果が記述されます。まず最初に基本並進ベクトル\(\mathbf{a}_{i} = (a_{\text{ix}},a_{\text{iy}},a_{\text{iz}})\)が次の形式で記述されます。

--- primitive lattice vectors ---
  a_1x a_1y a_1z
  a_2x a_2y a_2z
  a_3x a_3y a_3z

次に原子の数natmと各原子の座標\((x_{i},y_{i},z_{i})\)と質量\(m_{i}\)とラベルname(i)が次の形式記述されます。

--- Equilibrium position and mass of each atom---
 Natom = natm
 do i=1,natm
    i  x(i)  y(i) z(i) m(i) name(i)
 end do

次に振動解析の結果が次の形式で記述されます。

--- Vibrational modes ---
 Nmode= nmode, Natom= natm
 do m = 1,nmode
   n=   m  representation(m) acvtive(m)
   hbarW= omega_ha(m) Ha = omega_ev(m) eV; nu= omega_nu(m) cm^-1
   do i=1,natm
      i  vec(m,i,1) vec(m,i,2) vec(m,i,3)
   end do
 end do

representationは既約表現の配列です。active(m)はラマン活性なモードあれば Rになり,赤外活性なモードであればIRとなります。両活性であれば,IR&Rとなります。 サイレントモードの場合には何も表示されません。 vecは固有ベクトルの配列で,omega_haはHartree単位での振動数で, omega_evは電子ボルト単位での振動数で,omega_nuは波数です。

力データファイルforce.dataには力の定数を計算するための力のデータが記述されます。その力データは次の形式で出力されます。

num_force_data, norder, sw_ploynomial_fit
do i = 1, num_force_data
   displaced_atom, displacement(1:3)
   do ia = 1, natm
      i, force_data(ia,1:3,i)
   end do
end do

num_force_dataは力を計算する原子配置の数で,displaced_atomは変位した原子の番号で,配列displacementが原子の変位ベクトル\((u_{x},u_{y},u_{z})\)です。 norderは入力で指定したnorderの値が記述されます。sw_ploynomial_fitは入力のsw_ploynomial_fitがONのときに,ONを表す1が記述されます。OFFの場合には,0が記述されます。

sw_calc_forceをOFFに設定することで,出力された力を読み込み,振動解析をやり直すことができます。元素の質量を変更することは問題ありませんが,力計算に関連する変数は変更してはなりません。

一部の原子のみ振動解析の対象とする方法(バージョン2020.01以降)

バージョン2020.01以降より,一部の原子のみを振動解析の対象とすることができるようになりました。本計算機能を用いると,たとえば表面に吸着した分子の振動解析を行う際分子と表面数層のみを考慮することによって必要な原子間力計算回数を減らすことが可能になります。

基本的には通常の振動解析と同様の設定を施せばこの計算機能を用いることができます。各原子のmobile属性値を,変位の対象にしたい場合onに,したくない場合offに設定します。

structure{
  atom_list{
    atoms{
      #default mobile = on
      #tag element rx ry rz weight mobile
        Cu 0 0 0 1 off
        Cu 0 0.250000000015 0.090375000005 1 off
        Cu 0.250000000015 0 0.090375000005 1 off
        Cu 0.250000000015 0.250000000015 0 1 off
        Cu -0.00133991609598 -6.13618712873e-06 0.183172868018 on
        Cu 0.00243063003687 0.247502210432 0.274544823036 on
        Cu 0.248598303756 0.999991722423 0.274874965184 on
        Cu 0.249918873555 0.250806431959 0.18373526966 on
        Cu 0 0.500000000029 0 1 off
        ...
        ...
    }
  }
}

通常mobile属性値のデフォルト値はoffですが,振動解析の場合はonがデフォルト値となります。なお,mobile属性がoffの原子が存在する振動解析の場合は対称性を考慮しないようにしてください。

原子間力をファイルから読み込む方法(バージョン2020.01以降)

原子間力が計算済みの場合(force.dataファイルが存在し,原子間力データが記録されている場合),計算開始前にファイルから読み込み,その原子間力の計算をスキップさせることが可能となりました(2020.01以降この動作が規定の振る舞い)あえて原子間力を新たに計算し直したい場合は,以下の設定を施します。

phonon{
  sw_read_force_pre = off
}

phononブロックの下のsw_read_force_preによって可能な場合はファイルから原子間力を読み込むかどうかを設定します。このパラメーターのデフォルト値はonなので,offとすることによってこの振る舞いを抑制し,原子間力を計算し直すことができます。

計算例:水分子の振動解析

水分子の振動解析例を紹介します。入力データは samples/phonon/H2O 以下にあります。

構造最適化

振動解析を行うには原子が平衡状態になければなりませんので,振動解析を行うときと同じ条件で構造最適化を行います。平衡状態にないと動力学行列の固有値が負になり,振動数が純虚数のソフトモードが現れます。水分子の構造最適化の入力例を以下に示します。

control{
        condition = initial
        cpumax = 1 day ! maximum cpu time
        max_iteration = 6000
}
accuracy{
        cutoff_wf =   25.00  rydberg
        cutoff_cd =  225.00  rydberg
        num_bands = 8
        xctype = ggapbe
        initial_wavefunctions = matrix_diagon
        matrix_diagon {
          cutoff_wf = 5.0 rydberg
        }
        ksampling{
          method = gamma
        }
        scf_convergence{
          delta_total_energy = 1.e-10
          succession = 3
          num_max_iteration = 300
        }
        force_convergence{
          max_force = 1.e-4
        }
        initial_charge_density = Gauss
}
structure{
        unit_cell_type = primitive
        unit_cell{
              a_vector =  15.0               0.0                0.0
              b_vector =   0.0              15.0                0.0
              c_vector =   0.0               0.0               15.0
        }
        symmetry{
             tspace{
                lattice_system = primitive
                generators{
                   #tag rotation tx ty tz
                        C2z      0  0  0
                        IC2x     0  0  0
                }
             }
             sw_inversion = off
        }
        magnetic_state = para
        atom_list{
             coordinate_system = cartesian
             atoms{
                !#default mobile=on
                !#tag  rx             ry       rz          element
                      -1.45           0.000    1.123       H
                       1.45           0.000    1.123       H
                       0.0            0.0      0.0         O
             }
        }
        element_list{  #units atomic_mass
                       #tag element  atomicnumber zeta  dev
                            H             1       1.00  0.5
                            O             8       0.17  1.0    }
}
wf_solver{
        solvers {
        !#tag sol    till_n dts dte itr  var    prec cmix submat
              msd      5    0.1 0.1   1    tanh on   1    on
              lm+msd  10    0.1 0.4  50    tanh on   1    on
              rmm2p   -1    0.4 0.4   1    tanh on   2    on
        }
        rmm {
          edelta_change_to_rmm = 1.d-6
        }
        lineminimization {
          dt_lower_critical = 0.1
          dt_upper_critical = 3.0
        }
}
charge_mixing{
        mixing_methods {
        !#tag id method   rmxs rmxe itr var    prec istr nbxmix update
              1  broyden2 0.3  0.3  1   linear on   5    10     RENEW
              2  simple   0.2  0.5  100 linear on   *    *      *
        }
}
structure_evolution{
   method = cg
}

file_names.dataにはelement_listと同じ順番でポテンシャルファイルH_ggapbe_nc_01.ppとO_ggapbe_us_02.ppを指定します。この入力を使用して得た水分子の構造を 図 5.34 に示します。

../_images/image82.svg

水分子の構造

振動解析

構造最適化後に振動解析を行うには,入力の原子座標を最適化した座標に換えて,Phononブロックを加え,振動解析の設定をします。最適化原子座標は構造最適化計算の出力ファイルnfdynm.dataに記述されている最後のステップの原子座標です。

atom_list{
     coordinate_system = cartesian
     atoms{
        !#tag  rx             ry       rz          element
              -1.446816228    0.000    1.123327795 H
               1.446816228    0.000    1.123327795 H
               0.0            0.0      0.0         O
     }
}

振動解析の設定はたとえば以下のようにします。原子変位は0.05とします。

Phonon{
   sw_phonon = on
   sw_calc_force = on
   sw_vibrational_modes = on
   displacement = 0.05
}

PHASEを実行します。

% mpirun  ../../../bin/phase

PHASEを実行すると、振動解析結果のファイルmode.dataが出力されます。

振動数レベル図はツールfreq.plを使用して作成します。分子の場合には以下のように-molというオプションを付けてfreq.plを実行します。

% freq.pl -mol mode.data

この例題の水分子の基準モードの振動数を 図 5.35 に示します。

../_images/image83.svg

水分子の振動モードの振動数

基準振動の固有ベクトルの図を作成するための拡張trajectory形式のファイルは、ツールanimate.plで作成します。原点の移動を指定したファイルcontrol.inpを用意します。control.inpは以下のように記述します。

origin  7.5 7.5 7.5

ツールanimate.plを以下のように実行します。

% animate.pl mode.data control.inp

基準振動の固有ベクトルの拡張trajectory形式のファイルmode_*.tr2が生成されます。

この例題の水分子の基準モードの固有ベクトルを 図 5.36 に示します。生成された振動モードの拡張trajectory形式のファイルmode_7.tr2,mode_8.tr2,mode_9.tr2を可視化したものです。

../_images/image84.png

水分子の振動モードの固有ベクトル

計算例:シリコン結晶(Si2)

入力パラメータ

シリコン結晶の振動解析の例題です。計算例題は、 samples/phonon/Si2 です。

入力パラメータファイルnfinput.dataでは、element_listにシリコン原子の質量28.0855 amuが指定されています。質量の単位を原子質量単位とするため、#unitsの後にatomic_massを指定しています

element_list{     #units atomic_mass
         #tag element  atomicnumber mass
                             Si       14     28.0855
        }

振動解析のパラメータをPhononブロックで指定します。

Phonon{
   sw_phonon = on
   sw_calc_force = on
   displacement = 0.1
   sw_vibrational_modes = on
}

sw_calc_forceとsw_vibrational_modesがともにONなので,振動解析のための力計算を行い,振動解析が行われます。

PHASEを実行します。

% mpirun  ../../../bin/phase

計算が終了すると,出力ファイルmode.dataに振動解析の結果が出力されます。mode.dataの最初の部分は以下のようになっています。

--- primitive lattice vectors ---
  0.0000000000   5.0875600000   5.0875600000
  5.0875600000   0.0000000000   5.0875600000
  5.0875600000   5.0875600000   0.0000000000
--- Equilibrium position and mass of each atom---
Natom=    2
   1   1.2718900000   1.2718900000   1.2718900000    51196.42133 Si
   2  -1.2718900000  -1.2718900000  -1.2718900000    51196.42133 Si
--- Vibrational modes ---
Nmode=    6 Natom=    2
n=     1 T1u      hbarW =  0.00000000E+00 Ha =  0.00000000E+00 eV; nu=  0.00000000E+00 cm^-1
    1   0.0000000000  0.0000000000  0.7071067812
    2   0.0000000000  0.0000000000  0.7071067812

最初の二行目から三行目は基本並進ベクトルをあらわしています。六行目は原子数を表しています。その次の行からは,原子の番号,デカルト座標,質量,ラベルが一行にあらわされています。Vibrational modesというタイトル行の次の行にはモード数と原子数があらわされています。これ以降には各振動モードの既約表現を先頭行として,次行に振動数があらわされ,その次の行から固有ベクトルがあらわされています。固有ベクトルは原子の番号の後にその原子に帰属するベクトルの3成分があらわされています。

振動数レベル図

振動解析の出力ファイルmode.dataの振動数のデータから振動数レベル図を作成します。 以下のように、ツールfreq.plを実行すると、Postscript形式の振動数レベル図freq.epsが出力されます。

% freq.pl mode.data

シリコン結晶の振動解析の振動レベル図を 図 5.37 に示します。この図から振動数が517 cm-1であるモードがあることが分かります。 このモードの既約表現はT2gであるので,同じ振動数のモードが三重に縮重しています。T2gモードがラマン活性である場合,図中の規約表現の右側にRが表示されます。赤外活性である場合にはIRと表示されます。

../_images/image85.svg

バルクSiの領域中心フォノンモードの振動数

振動モードの可視化

振動解析の出力ファイルmode.dataから拡張Trajectory形式のファイルを作成することにより,固有ベクトルを矢印表示したり,原子が振動するアニメーションとして振動モードを可視化したりできます。ツールanimate.plを使用して、振動解析の出力ファイルmode.dataから振動数の情報を取り出し,拡張Trajectory形式のファイル(拡張子:tr2)を作成します。

原点の移動とセルベクトルの変更を指定したファイルcontrol.inpを用意します。control.inpは以下のように記述します。

origin  1.27189 1.27189 1.27189
vector1 10.17512 0 0
vector2 0 10.17512 0
vector3 0 0 10.17512

ツールanimate.plを以下のように実行すると,拡張Trajectory形式のファイルがモードの数だけ出力されます。

% animate.pl mode.data control.inp

この例題では切り出すセルをブラベー格子の単位胞にとり,セルの原点にシリコン原子がくるように設定しています。 たとえば,出力された拡張Trajectory形式のファイルmode_6.tr2を可視化すると, 図 5.38 のように固有ベクトルが矢印で示されます。 図 5.38 に示されているセルは,出力されたgrid.mol2ファイルを読み込ことで表示できます。また,出力された拡張Trajectory形式から、原子の振動を可視化することができます。

../_images/image86.svg

バルクSiの領域中心フォノンモードの固有ベクトル

計算例:銅 (100)面に吸着したエチレン分子の振動解析

Cu(100)面にエチレン分子を配置した系の振動解析を実施しました。サンプルデータは samples/phonon/Cu100_C2H4 の下にあります。まずはCu(100)面にエチレン分子を吸着させ,通常の構造最適化を実施しました。 結果得られた原子配置は 図 5.39 に示す通り。

../_images/image87.jpeg

Cu (100)面にC2H4分子を吸着させた系。

分子のみを考慮するケース(サンプルデータが置かれたディレクトリーはmolonly),分子と表面第一層まで考慮するケース(サンプルデータが置かれたディレクトリーはph1),分子と表面第二層まで考慮するケース(サンプルデータが置かれたディレクトリーはph2),そしてすべての原子を考慮する(サンプルデータが置かれたディレクトリーはphall)振動解析を実施しました。比較のため,孤立分子の振動解析も行いました。得られた振動数の上位5つを 表 5.13 に報告します。

振動数top 5(単位:cm-1

分子のみ

分子+1層

分子+2層

全原子

孤立分子

3041.4

3039.7

3039.4

3039.4

3057.8

3014.0

3012.2

3012.0

3012.0

3028.5

2946.3

2944.3

2944.1

2944.1

2973.7

2942.2

2940.1

2939.8

2939.8

2957.6

1475.9

1478.4

1478.7

1478.9

1618.1

分子のみ考慮する場合ではやくも3 cm-1以下の精度で振動数が計算できていることが分かります。1層取り込むとこれが0.5 cm-1以下となり,ほぼ収束していると考えてよいでしょう。2層取り込んだ結果は最大でも0.2 cm-1の誤差となっており,これは同じ結果が得られているといってもよい精度です。

フォノンバンド計算

機能の概要

PHASEには, \(\Gamma\) 点だけでなく一般のk点における格子振動解析を行い、フォノンの状態密度やバンド構造を計算する機能があります。 文献 [Parlinski97] のアルゴリズムに従って計算します。

利用方法

基本的な入力パラメータ

この機能を利用するためには, \(\Gamma\) 点の場合と同様phononブロックを作成し,設定を行います。典型的には,以下のようになります。

phonon{
  sw_phonon = on
  sw_vibrational_modes = on
  sw_calc_force = on
  displacement = 0.1
  method = dos
  lattice{
    l1 = 2
    l2 = 2
    l3 = 2
  }
  dos{
    mesh{
      nx = 10
      ny = 10
      nz = 10
    }
  }
}

基本的な入力パラメーターを以下に示します。

第1 ブロック識別子

第2,第3 ブロック識別子

タグ識別子

説明

phonon

フォノン計 算の設定を行う ためのブロック

sw_phonon

PHASEの

振動解析機能を 利用するかどう かを指定します 。Γ点のみの場 合と同様です。

sw_vib rational_modes

振動解 析を行うかどう かを指定します 。Γ点のみの場 合と同様です。

sw_calc_force

力定数の計 算を行うかどう かを指定します 。Γ点のみの場 合と同様です。

displacement

力定数の計 算を行う際に, 原子をどの程度 移動させるかを 指定します。Γ 点のみの場合と 同様です。デフ ォルト値は0.1 bohrです。

method

“手 法”を指定しま す。状態密度計 算の場合はdos, バンド構造計 算の場合はband を指定します。

lattice

一 般k点 の振動解析は, スーパーセルの 力定数を必要と します。そのス ーパーセルの大 きさを指定する ブロックです。

nx

a軸を何倍するか を指定します。

ny

b軸を何倍するか を指定します。

nz

c軸を何倍するか を指定します。

dos

状態密度の計算 方法を指定する ブロックです。

mesh

状態密度計算 に利用する k 点のメッ シュを指定する ブロックです。

nx

1 つめの逆格子ベ クトルの分割数 を指定します。

ny

2 つめの逆格子ベ クトルの分割数 を指定します。

nz

3 つめの逆格子ベ クトルの分割数 を指定します。

methodをbandと設定すると,フォノンバンドの計算になります。バンド構造の計算は,電子バンド構造と同様band_kpoint.plを利用して計算するk点の情報が記録されたkpoint.dataファイルを作成したあとに実行します。

バンド数およびk点分割数

フォノンバンド計算を行う場合,スーパーセルの作成が行われます。バンド数やk点分割数は,生成されたスーパーセルに合わせてPHASEが自動的に変更します。以下の注意が必要です。

  • スーパーセルは,ブラベー格子に対して作成されます。通常のPHASE による計算の場合unit cell type がBravais の場合は基本格子に対してk 点を定義しますが,フォノンバンドの場合はブラベー格子に対して行うようにしてください。

  • バンド数は,定義した原子にしたがって通常の計算と同じように決定してください。

イオン性結晶の場合の設定方法

イオン性結晶の場合,Γ点においてオプティカルモードの縦波と横波が異なった振動数を持ちます。この現象は,LO-TO 分裂と呼ばれます。この効果を取り入れる場合,入力ファイルにおいてさらに以下の指定を行う必要があります。

phonon{
  sw_lo_to_splitting = on
  electronic_dielectric_constant{
    exx = 2.6
    eyy = 2.6
    ezz = 2.6
    exy = 0.0
    exz = 0.0
    eyz = 0.0
  }
}

変数sw_lo_to_splitting をon とすることによってLO-TO 分裂を考慮した計算を行うことができます。electronic_dielectric_constantブロックには,電子系の誘電テンソルを指定します。electronic_dielectric_constantブロックの下のexx, eyy, ezz, exy, exz, eyz に対応する誘電テンソルの成分を指定します。誘電テンソルは,実測値もしくはUVSOR-Epsilon による計算値をご利用ください。原子の有効電荷も指定する必要があります。これは,作業ディレクトリーにeffchg.data ファイルを作成し,以下のように指定します。

2
1
1.12 0.0 0.0
0.0 1.12 0.0
0.0 0.0 1.12
2
-1.12 0.0 0.0
0.0 -1.12 0.0
0.0 0.0 -1.12

ファイルの1 行目に原子数を記述します。2 行目以降に有効電荷の値を指定します。まず指定対象の原子のID を指定し,さらに有効電荷テンソルを指定します。有効電荷テンソルは,形式電荷を利用することもできますが,

UVSOR-Berry によって得られたボルン有効電荷を利用することが望ましいです。

計算の実行

入力データが準備できたら,通常通りPHASE を実行します。まずは,PHASE は入力の指定にしたがってスーパーセルを作成します。ログファイルには以下のように報告されます。

natm_super,natm2_super= 64 64
ia,cps(3),pos(3),ityp
1 1.27189 1.27189 1.27189 0.06250 0.06250 0.06250 1
2 8.90323 8.90323 8.90323 0.43750 0.43750 0.43750 1
3 1.27189 6.35945 6.35945 0.06250 0.31250 0.31250 1
4 8.90323 13.99079 13.99079 0.43750 0.68750 0.68750 1
5 6.35945 1.27189 6.35945 0.31250 0.06250 0.31250 1
6 13.99079 8.90323 13.99079 0.68750 0.43750 0.68750 1
7 6.35945 6.35945 1.27189 0.31250 0.31250 0.06250 1
8 13.99079 13.99079 8.90323 0.68750 0.68750 0.43750 1
9 11.44701 1.27189 1.27189 0.56250 0.06250 0.06250 1
10 19.07835 8.90323 8.90323 0.93750 0.43750 0.43750 1
11 11.44701 6.35945 6.35945 0.56250 0.31250 0.31250 1
12 19.07835 13.99079 13.99079 0.93750 0.68750 0.68750 1
13 16.53457 1.27189 6.35945 0.81250 0.06250 0.31250 1
...
...

natm_super がスーパーセルの原子数です。cps は原子のカルテシアン座標, pos はフラクショナル座標です。itypは原子の種類を識別する番号です。また, スーパーセルに合わせて変化したバンド数とk点のメッシュが次のように報告されます。

num_bands will be changed.
neg,meg= 192 192
k-point mesh will be changed.
mesh= 1 1 1

neg が新しいバンド数, mesh が新しいk 点メッシュです。

出力ファイル

mode.dataファイル

振動解析の結果はmode.data ファイルに記録されます。フォノンバンドの場合のmode.dataファイルは,たとえば以下のようになります。Γ点の場合と比較して,振動モードの記述の仕方が異なります。

--- Vibrational modes ---
Nmode= 6 Natom= 2 Nqvec 120
iq= 1 q=( 0.00000, 0.00000, 0.00000) ( 0.00000, 0.00000, 0.00000)
n= 1 T1u IR
hbarW= 0.00000000E+00 Ha = 0.00000000E+00 eV; nu= 0.00000000E+00 cm^-1
1 0.0000000000 0.7071067812 0.0000000000
2 0.0000000000 0.7071067812 0.0000000000
1 0.0000000000 0.0000000000 0.0000000000
2 0.0000000000 0.0000000000 0.0000000000
n= 2 T1u IR
...
iq= 2 q=( 0.01875, 0.01875, 0.03750) ( 0.02316, 0.02316, 0.00000)
n= 1 B2 IR&R
hbarW= 0.63506708E-04 Ha = 0.17281054E-02 eV; nu= 0.13938112E+02 cm^-1
1 0.4999599615 -0.4999599615 0.0000000000
2 0.4999599615 -0.4999599615 0.0000000000
1 0.0063274755 -0.0063274755 0.0000000000
2 0.0063274755 -0.0063274755 0.0000000000
n= 2 B1 IR&R
...

モードの数と原子数の後に,k点の数が示されます。各k 点の振動数モードの記述の前に, k 点の座標が部分座標とカルテシアン座標で示されます。振動モードの固有ベクトルは一般には複素数となるので, 固有ベクトルの実部の後に, 虚部が記述されます。なお,Γ点の場合と同様に振動の対称性および赤外ラマンの活性/不活性の判定が出力されますが,この情報はΓ点以外では意味がない点にご注意ください。

phdos.dataファイル

フォノンの状態密度はphdos.data ファイルに出力されます。その内容は,典型的には下記のようなものです。

# Index Omega(mHa) Omega(eV) Omega(cm-1) DOS(States/Ha) DOS(States/eV) DOS(States/cm-1) IntDOS(States)
0 -0.00050000 -0.00001361 -0.10973732 0.00000000  0.00000000 0.00000000 0.00000000
1  0.00950000  0.00025851  2.08500903 0.00473815  0.17412390 0.00002159 0.00001500
2  0.01950000  0.00053062  4.27975539 0.01996324  0.73363561 0.00009096 0.00012976
3  0.02950000  0.00080274  6.47450174 0.04568839  1.67901746 0.00020817 0.00044927
4  0.03950000  0.00107485  8.66924810 0.08191360  3.01026946 0.00037323 0.00107853
5  0.04950000  0.00134696 10.86399446 0.24722290  9.08527497 0.00112643 0.00286860
6  0.05950000  0.00161908 13.05874081 0.37130693 13.64527929 0.00169180 0.00591423
7  0.06950000  0.00189119 15.25348717 0.49343689 18.13347292 0.00224826 0.01020273
8  0.07950000  0.00216331 17.44823352 0.67844022 24.93222060 0.00309120 0.01602478
...................
...................

1 列目は状態密度のインデックス,2, 3, 4 列目がそれぞれmHa, eV, cm-1 単位のエネルギー,5, 6, 7 列目がそれぞれstates/Ha, states/eV, states/cm-1 単位での状態密度,8 列目が積算状態密度です。積算状態密度は,最も高エネルギーの状態においては原子数×3 になります。

解析用Perlスクリプト

フォノンバンド計算の結果解析用のPerl スクリプトがPHASE には備わっています。以下の3 種類のPerl スクリプトを利用して結果の解析を行うことができます。

phonon_dos.pl

フォノンの状態密度データから,「フォノン状態密度図」を作成するPerlスクリプトです。以下のように使用します。

% phonon_dos.pl phdos.data OPTIONS

phdos.data が,PHASE が出力するフォノン状態密度データです.実行すると,phonon_dos.eps というEPS形式の画像ファイルが作成されます。下記のオプションを利用することができます。

--units=UNITS or -u UNITS

エ ネルギーの単位を指定します。mHa, meV, THz, cm-1 の いずれかです。デフォルト値はcm-1 です。

--width=WIDTH or -w WIDTH

作成される図の 幅を指定します。デフォルト値は1 です。

--erange=[emin:emax] or -e [emin:emax]

エネルギーの範囲を指定します。

--drange=DRANGE or -d DRANGE

状態密度の範囲を指定します。

--title=TITLE or -t TITLE

図のタイトルを指定します。

--font=FONT or -f FONT

グラフに利用するフォントサイ ズを指定します。デフォルト値は18 です。

--keep or -k

中 間のデータファイルを保持する場合 ,このオプションを有効にします。

--mono or -m

モノクログラフを描画したい場 合にこのオプションを指定します。

--dinc=DINC

状態密度の目盛を指定します。

--einc=EINC

エネルギーの目盛を指定します。

phonon_band.pl

フォノンバンドのデータから「フォノンバンド図」を作成するPerl スクリプトです。以下のように使用します。

% phonon_band.pl mode.data OPTIONS

mode.data が,振動解析の結果が記録されたファイルです。実行すると,phonon_band.eps というEPS 形式の画像ファイルが作成されます。オプションには,下記のようなものがあります。

--control=CONTROL or -c CONTROL

band_kpoint.pl ファイルの入力ファイルを指定 します。デフォルト値はbandkpt.in です。

--ptype=PTYPE or -p PTYPE

グラフ描画に利用 するプロット種を指定します。line を指定すると実線,circle を指定すると丸でフォノンバンド を描画します。デフォルト値はline です。

--units=UNITS or -u UNITS

エ ネルギーの単位を指定します。mHa, meV, THz, cm-1 の いずれかです。デフォルト値はcm-1 です。

--width=WIDTH or -w WIDTH

作成される図の 幅を指定します。デフォルト値は1 です。

--erange=[emin:emax] or -e [emin:emax]

エネルギーの範囲を指定します。

--title=TITLE or -t TITLE

図のタイトルを指定します。

--font=FONT or -f FONT

グラフに利用するフォントサイ ズを指定します。デフォルト値は18 です。

--mono or -m

モノクログラフを描画したい場 合にこのオプションを指定します。

--keep or -k

中 間のデータファイルを保持する場合 ,このオプションを有効にします。

--einc=EINC

エネルギーの目盛を指定します。

phonon_energy.pl

フォノンの状態密度から,振動に由来する内部エネルギーとヘルムホルツの自由エネルギーや比熱を計算するスクリプトです。振動に由来するヘルムホルツの自由エネルギーに通常のDFT計算で得られる全エネルギーを加えれば,有限温度における固体の自由エネルギーを計算することができ,温度誘起の相転移を解析することなども可能です。

波数kのフォノンのエネルギーは,その振動数\(\omega_{\mathbf{k}}\)を利用して\(\left( \frac{1}{2} + n \right)\hbar\omega_{k}\)と記述することができます。分配関数は\(Q_{k} = \sum_{n}^{}e^{\frac{- U_{\text{kn}}}{k_{B}T}}\)と記述されるので,フォノンのエネルギーを代入し整理すると以下の結果が得られます。

\[Q_{k} = \frac{\exp\left\lbrack - \frac{\hbar\omega_{k}}{2k_{B}T} \right\rbrack}{1 - \exp\left\lbrack - \frac{\hbar\omega_{k}}{k_{B}T} \right\rbrack}.\]

ヘルムホルツの自由エネルギーは\(F_{\text{vib}} = \sum_{k}^{}{- k_{B}T\log Q_{k}}\)と記述できるので,以下のように記述されます。

\[F_{\text{vib}} = \sum_{k}^{}\left\lbrack \frac{\hbar\omega_{k}}{2} + k_{B}T\log\left( 1 - \exp\left\lbrack - \frac{\hbar\omega_{k}}{k_{B}T} \right\rbrack \right) \right\rbrack.\]

振動による平均の内部エネルギーは\(U_{\text{vib}} = \frac{1}{Q_{k}}\sum_{n}^{}{U_{\text{kn}}\exp\left( - \frac{U_{\text{kn}}}{k_{B}T} \right)}\)と記述できるので,以下のように記述することができます。

\[U_{\text{vib}} = \sum_{k}^{}{\left\lbrack \frac{\hbar\omega_{k}}{2} + \frac{\hbar\omega_{k}}{\exp{\left\lbrack \frac{\hbar \omega}{k_{B}T} \right\rbrack - 1\ }} \right\rbrack.}\]

エントロピーは,\(F_{\text{vib}}\)および\(U_{\text{vib}}\)から\(\left( U_{\text{vib}} - F_{\text{vib}} \right)/T\) と計算することができます。定積比熱は内部エネルギーの温度に関する偏微分で与えられるので,以下のように計算することができます。

\[C_{v} = \frac{\partial U_{\text{vib}}}{\partial T} = k_{B}\left\lbrack \frac{\frac{\hbar\omega_{k}}{k_{B}T}\exp{\left( \frac{\hbar\omega_{k}}{2k_{B}T} \right)}}{\exp{\left( \frac{\hbar\omega_{k}}{k_{B}T} \right) - 1}} \right\rbrack^{2}.\]

phonon_energy.plを利用すると,以上のような計算を実行することが可能です。以下のように利用します。

% phonon_energy.pl mode.data

この操作によって,以下の3つのファイルが作成されます。

phonon_energy.dataファイル 内部エネルギー,ヘルムホルツ自由エネルギー,エントロピー,比熱の計算結果が以下の形式で記録されているファイルです。

# T (K) Internal Energy (eV) Free energy (eV) Entropy (eV/K) Cv (kB/atom)
0 0.125434126153072 0.125434126153072 0 0 30 0.12552700746085 0.125409486111375 3.91737831580881e-06
0.0820122071540538 60 0.126828216477476 0.124936822438767 3.15232339784872e-05
0.435633166874193 90 0.130001095247047 0.123379006005857 7.35787693465625e-05
0.787404251770626 120 0.134948880737123 0.120473935403623 0.000120624544445835 1.12444793146534
.......
.......
.......

1列目に温度が,2列目以降からは内部エネルギーとヘルムホルツ自由エネルギーがeV単位で,エントロピーがeV/K単位で,原子あたりの比熱がkB単位で記述されます。

phonon_energy.eps ファイル 内部エネルギー,ヘルムホルツ自由エネルギー,エントロピーを温度の関数としてプロットしたEPSファイルです。

phonon_Cv.eps ファイル 比熱と温度の関係をプロットしたEPSファイルです。

Si結晶の場合に得られるphonon_energy.epsおよびphonon_Cv.epsの例を示します。phonon_energy.plスクリプトは,フォノン状態密度計算を実行した際に得られるmode.dataファイルを利用する必要がある点にご注意ください。フォノンバンド計算を実行した際に得られるmode.dataファイルを利用すると,以下のようなエラーが発生し途中で終了してしまいます。なお,得られるエネルギーは,入力で指定した原子数分となります。

% phonon_energy.pl mode.data
weight undefined for q-point no. 1 at /home/user/phase/bin/phonon_energy.pl line 131, <MD> line 4450.
../_images/image91.png

phonon_energy.eps(左図)とphonon_Cv.eps(右図)の例

phonon_energy.plスクリプトのオプションは,下記の通りです。

--width=WIDTH or -w WIDTH

作成される図の 幅を指定します。デフォルト値は1 です。

--trange=[tmin:tmax] or -t [tmin:tmax]

温度の範 囲を指定します。デフォルト値は0 Kから3000 Kです。

--nT=NT or -n NT

温度の点の数を指定 します。デフォルト値は100です。

--font=FONT or -f FONT

グラフに利用するフォントサイ ズを指定します。デフォルト値は18 です。

--mono or -m

モノクログラフを描画したい場 合にこのオプションを指定します。

--tinc=TINC

温度の目盛を指定します。

--einc=EINC

エネルギーの目盛を指定します。

--cinc=EINC

比熱の目盛を指定します。

例題

シリコン結晶

最も簡単な例の1 つとして,シリコン結晶のフォノンバンドとフォノン状態密度の計算を実行した例を紹介します。 この例題の入力ファイルは, samples/phono_band/Si 以下にあります。

まずはバンド計算を行います。 samples/phonon_band/Si/band 以下の入力ファイルを利用します。

band_kpoint.pl 用の入力ファイル,bandkpt.in の内容は,以下のようになっています。

0.02
-0.8333333 0.8333333 0.8333333
0.8333333 -0.8333333 0.8333333
0.8333333 0.8333333 -0.8333333
0 0 0 1 # {/Symbol G}
1 1 0 2 # X
5 3 0 8 # U
0 0 0 1 # {/Symbol G}
1 0 0 2 # L

このbandkpt.in ファイルを利用して,以下のようにkpoint.data ファイルを作成します。

% band_kpoint.pl bandkpt.in

入力の,原子配置の指定は以下のようになっています。

structure{
  unit_cell_type = bravis
  unit_cell{
    a = 10.17512
    b = 10.17512
    c = 10.17512
    alpha = 90.0
    beta = 90.0
    gamma = 90.0
  }
  symmetry{
    tspace{
      lattice_system = facecentered
    }
    method = automatic
  }
  atom_list{
    coordinate_system = internal
    atoms{
      #tag element rx ry rz mobile
          Si 0.125 0.125 0.125 0
          Si 0.875 0.875 0.875 0
    }
  }
  element_list{
    #units atomic_mass
    #tag element atomicnumber mass
        Si 14 28.0855
  }
}

unit_cell_type をbravais とし,lattice_system パラメータによってこの系がfacecentered, すなわち面心であることを指定しています。上述したように,通常のPHASE の計算ではこのような指定がなされている場合単位胞を基本格子に変換しますが,フォノンバンド計算ではそのようなことは行われません。次に,phonon ブロックを次のように記述しています。

Phonon{
  sw_phonon = on
  sw_calc_force = on
  sw_vibrational_modes = on
  lattice{
    l1 = 2
    l2 = 2
    l3 = 2
  }
  method = band
}

スーパーセルは,a, b, c 軸それぞれを2 倍とする設定を採用しています。以上の設定のもとPHASE を通常通り実行します。計算が終了すると,その結果がmode.data に記録されます。mode.data ファイルからフォノンバンド図を得るためには,以下の操作を行います。

% phonon_band.pl mode.data --control=bandkpt.in

--control オプションでband kpoint.pl 用の入力ファイルを指定していますが,この指定がない場合はバンド図に特殊点を表す縦線などは描画されなくなります。結果は次に示す図のようになります。

../_images/image92.png

シリコン結晶のフォノンバンド

フォノンの状態密度の計算に必要な入力データは, samples/phonon_band/Si/dos 以下にあります(力定数は計算済みなので,bandディレクトリーの下にあるforce.dataファイルをコピーして利用すると力定数計算をスキップすることも可能ですが,この場合はsw_calc_forceパラメータをoffに設定してください)このサンプルの入力パラメータファイルには,以下のような記述がなされています。

Phonon{
  sw_phonon = on
  sw_vibrational_modes = on
  lattice{
    l1 = 2
    l2 = 2
    l3 = 2
  }
  dos{
    mesh{
      nx = 10
      ny = 10
      nz = 10
    }
  }
  method = dos
}

method = dos と指定することによって状態密度計算を行うことを指定しています。dos ブロックの下のmesh ブロックにおいて,状態密度計算で利用するメッシュを10×10×10 としています。入力をこのように編集し終えたらPHASE を実行します。フォノン状態密度の計算結果はphdos.data ファイルに記録されます。このデータをもとにphonon dos.pl スクリプトを利用してフォノン状態密度図を作成します。

% phonon_dos.pl phdos.data

この結果得られるフォノン状態密度図は次に示す通りです。

../_images/image93.png

シリコン結晶のフォノン状態密度。

ヨウ化カリウム

ヨウ化カリウムはNaCl 型の結晶構造をもつ,イオン性の結晶です。 図 5.43 にその結晶構造を示します。ここでは,この結晶を例にLO-TO 分裂を考慮したフォノンバンド計算を紹介します。この例題の入力ファイルは, samples/phonon_band/KI 以下にあります。

../_images/image94.png

ヨウ化カリウムの結晶構造

LO-TO分裂を考慮した計算を行うためには,電子系の誘電テンソルと有効電荷が必要です。これらは以下のようにして得ました。

  • 誘電テンソル:UVSOR-Epsilonを利用して計算しました。この際,2.2 eVのギャップ補正を施しました。結果は,xx, yy, zz方向がそれぞれ2.6となりました。

  • 有効電荷:UVSOR-Berryを利用して,ボルン有効電荷テンソルの計算を行いました。結果は,カリウムの有効電荷が1.1262, ヨウ素の有効電荷が-1.1262となりました。

これらを設定し,sw_lo_to_splittingをonとする以外はシリコン結晶の場合と同じです。に,得られたフォノンバンドを示します。比較のため,LO-TO分裂を考慮せずに計算した結果も合わせて表示しています。赤線がLO-TO分裂を考慮せずに計算した結果,青線が考慮して計算した結果に対応します。この図から明らかなように,Γ点付近ではLO-TO分裂によって考慮しない場合は縮退している状態が分裂しています。

../_images/image95.png

KI 結晶のフォノンバンド。赤線がLO-TO分裂を考慮せずに計算した結果,青線がLO-TO 分裂を考慮して計算した結果。

スズの温度誘起相転移

最後に,フォノン自由エネルギー解析の簡単な適用例としてスズの温度誘起相転移の例を紹介します。この例題の入力ファイルは, samples/phonon_band/Sn/a-Sn (αスズ)および samples/phband/Sn/b-Sn (βスズ)にあります。

スズには,αスズとβスズと呼ばれる同素体があります。αスズはダイヤモンド構造,βスズはその名の通りβスズ構造をとります。その結晶構造を, 図 5.45 に示します。

../_images/image96.png

αスズ(左図)とβスズ(右図)の結晶構造

βスズ構造はダイヤモンド構造をc軸方向から押しつぶしたような結晶構造であり,体心正方晶を取ります。常温ではβスズが安定ですが,低温下ではαスズが安定になります。これは,結晶そのものの全エネルギーはαスズの方が低いが,温度上昇に伴うフォノンの自由エネルギーの低下はβスズの方が大きいためある温度で自由エネルギーはβスズの方が低くなり,相転移するからであると考えられます。このような現象を,フォノンの自由エネルギー計算と結晶の全エネルギー計算を組み合わせて確認していきます。

まずは,格子定数の最適化を行いました。ただし,βスズ構造のc/a比は0.54614と固定して最適化しました。結果は,次の表に示す通りです。

格子定数a (Å)

格子定数 c (Å)

全エネルギー (ha./cell)

αスズ

6.6555

6.6555

-136.147884

βスズ

5.9184

3.2323

-136.144694

この結果から明らかなように,全エネルギーはαスズの方が低いので,絶対零度ではαスズが安定であると考えられます。

続いて,得られた安定な格子定数のもとで振動解析を行いました。自由エネルギーを評価する場合に必要な計算は,状態密度のみです。αスズ,βスズに対してシリコンの場合と同様の設定をPhononブロックで行い,振動解析を実施しました。計算終了後に得られたmode.dataファイルを,phonon_energy.plスクリプトで処理します。

% phonon_energy.pl mode.data

結果得られるphonon_energy.dataファイルの3列目にフォノンの自由エネルギーが記録されます。これは単位胞あたりの値なので,上述の全エネルギーの単位胞あたりのエネルギーを加え,温度の関数としてプロットすると 図 5.46 のような結果が得られます。

../_images/image97.png

αスズとβスズの自由エネルギーと温度の関係。赤線がαスズ,緑線がβスズに対応する。

図 5.36 においてαスズの曲線(赤線)とβスズの曲線(緑線)が交差する温度が転位温度と考えられます。この計算ではおおよそ510 Kとなりました。実際には290 Kなので相転移温度が高く評価されてしまいましたが,このような計算によって温度誘起の構造相転移を説明できることはお分かりいただけたと思います。

使用上の注意

  • 一般のk 点における振動解析を実行するためには,スーパーセルに対する振動解析を行う必要があります。したがって,Γ点のみの場合と比較すると非常に多くの計算時間が必要です。

  • フォノンバンド計算において最も計算量が多いのがスーパーセルに対する力の計算です。このデータは,1 度得られたら再利用することができます。たとえばフォノンの状態密度を計算したあとにフォノンバンドを計算する場合,または異なる対称線にそったフォノンバンドを計算する場合などは,以下のようにsw_calc_ forceパラメータをoff とすることによって力計算をやり直すことをさけることができます。なお,力計算の結果が保存されているファイルはforce.dataというファイルです。バンド計算と状態密度計算を異なるディレクトリで行う場合にsw_calc_force = offとする場合はこのファイルを当該ディレクトリにコピーして利用してください。

phonon{
  sw_phonon = on
  sw_calc_force = off
}
  • 通常の計算の場合,入力で指定されたBravais格子はPrimitive 格子に変換されて計算が行われます。ところが,フォノンバンド計算の場合はこの変換は実施されず,Bravais 格子のままスーパーセルが作成され,計算が行われます。k点サンプリングメッシュを検討する際などに注意が必要です。

  • スーパーセル構築のパラメータ,l1, l2, l3 はもとの対称性を保つような指定の仕方をしてください。異なる対称性の場合,意味のある計算は行われません。

参考文献

Parlinski97
  1. Parlinski, Z. Q. Li and Y. Kawazoe, Physical Review Letters vol. 78 pp. 4063 (1997)

分子動力学法シミュレーション

機能の概要

PHASEは、原子に働く力を利用して分子動力学法シミュレーションを行うことが可能です。 エネルギー一定,温度一定,圧力一定の分子動力学シミュレーションなどが実行できます。

入力パラメータ

分子動力学法シミュレーション機能と関連あるタグの一覧を表に示します。

分子動力学法シミュレーション機能に関連のあるタグの一覧

第1 ブロック識別子

第2,第3 ブロック識別子

タグ識別子

説明

struc ture_evolution

原子座標データ の更新方法を指 定するブロック

method

原子座標 の更新方法を指 定する。分子動 力学シミュレー ションの場合,

v elocity_verlet (エ ネルギー一定の 分子動力学シミ ュレーション)

tempe rature_control (温度一定の 分子動力学シミ ュレーション)

pr essure_control (エネルギ ー・圧力一定の 分子動力学シミ ュレーション)

temperature_pr essure_control (温 度・圧力一定の 分子動力学シミ ュレーション)

dt

時間刻 みを指定する。

デフ ォルト値は100 au (約2.4 fs)

tempe rature_control

温 度制御の設定を 行うブロック。

method

温度 制御の方法を指 定する。nose_ hooverかveloci ty_scalingのい ずれか。nose_ hooverの場合は Nose-Hoover熱 浴による温度制 御が,velocity _scalingの場合 は温度スケーリ ングによる温度 制御が行われる 。デフォルト値 はnose_hoover.

sw_r ead_velocities

原子の初期速 度を,PHASE/0 に自動生成させ るのではなく手 動で入力する場 合にこのパラメ ータをonとしま す。デフォルト 値はoffです。

set_in itial_velocity

原子 の初期速度をプ ログラムが自動 的に設定するか どうかを指定す るスイッチ。デ フォルト値はon

sw_sh ift_velocities

全 運動量がゼロに なるよう,MDス テップごとに速 度をシフトする かどうかを指定 するスイッチで す。デフォルト 値はoffです。

thermostat

熱浴の設定を 行うブロック。 表形式データ。

temp

温 度を指定する。

qmass

熱浴の質 量を指定する。

tdamp

熱浴の 質量を直接指定 するのではなく 、その周期を時 間の単位で指定 する。qmassに よる指定の方が 優先される。な お、qmassもtda mpも指定がない 場合,tdamp=5 0×dtに相当する 周期の質量がデ フォルト値とし て採用される。

pr essure_control

圧 力制御の設定を 行うブロック。

pressure

目的の圧 力を指定する。

mass_baro

バー ロスタットの質 量を指定する。

m11,m12,m13

m21,m22,m23

m31,m32,m33

格 子の制御に拘束 条件を加える。 たとえば,m11 = offとすると11 成分(a軸・ a軸)が変 化しなくなる。

structure

atom_list

atoms

原子配置を指定 するブロック。 表形式データ。

mobile

原子 が“可動”かどう かを指定する。 可動にする場合 onを指定する。

thermo group

原子に熱浴を割 り当てる。定義 した順に,整数 値で指定する。

デ フォルト値は0 (熱浴 に割り当てられ ていない状態)

vx,vy,vz

原子 の初期速度を手 動入力する場合 ( sw_read_veloci ties=onの場合) に各原子の速度 のx,y,z値を原 子単位で入力す る。入力が省略 される場合,0 とみなされる。

element_list

元素情 報を指定する。 表形式データ。

mass

対応する元素の 質量を指定する 。デフォルト値 は原子単位であ ることに注意。

計算結果の出力

座標データはfile_names.dataファイルのF_DYNMによって指定されるファイルに各ステップでの座標値が出力されます。 その形式は、構造最適化の場合と同様です。

  • 原子座標

原子座標は,構造緩和を行った場合と同様,file_names.data中のF_DYNM識別子によって指定されるファイル (既定のファイル名はnfdynm.data)に記述されます。入力においてprintoutlevelブロックの下のiprivelocity変数の値を2以上にしていた場合 (分子動力学シミュレーションの場合はデフォルト値),各原子の速度のデータも出力されます。速度のデータは,力のデータのあとに原子単位で出力されます。

  • 各ステップでのエネルギー

各ステップでのエネルギーは,file_names.data中のF_ENF識別子によって指定されるファイル (既定のファイル名はnfefn.data)に出力されます。サンプルによって得られる結果を以下に記します。

iter_ion,iter_total,etotal,ekina,econst,forcmx
   1      18     -7.8953179624     0.0000042358    -7.8953179624     0.0186964345
   2      30     -7.8953851218     0.0000665502    -7.8953185716     0.0183575424
   3      43     -7.8955768901     0.0002565396    -7.8953203505     0.0173392067
   4      56     -7.8958649874     0.0005418445    -7.8953231430     0.0156398790
   5      69     -7.8962052587     0.0008785990    -7.8953266596     0.0132645441
   6      83     -7.8965425397     0.0012120826    -7.8953304571     0.0102355854
   7      97     -7.8968179539     0.0014840140    -7.8953339398     0.0066063151
   8     111     -7.8969784478     0.0016420281    -7.8953364197     0.0024736141
   9     125     -7.8969875377     0.0016502900    -7.8953372478     0.0020111576
  10     139     -7.8968352058     0.0014992046    -7.8953360011     0.0066379641
  11     153     -7.8965440599     0.0012113794    -7.8953326806     0.0111430822
                          ...............................
                          ...............................
                          ...............................

一列目は原子座標の更新回数,二列目は電子のSCF計算の回数です。三列目は,系の内部エネルギー, 四列目は系の運動エネルギーです。五列目は系の内部エネルギーと運動エネルギーを足した値であり,エネルギー一定の分子動力学シミュレーションにおける保存量です。

使用方法:エネルギー一定の分子動力学シミュレーション

エネルギー一定の分子動力学シミュレーションの入力パラメータ例です。

計算例題は、 samples/dynamics/molecular_dynamics/NVE です。

accuracy{
    cutoff_wf = 9.00 rydberg
    cutoff_cd = 36.00 rydberg
    num_bands = 8
    xctype = ldapw91
    force_convergence{
        max_force = 1.0e-8 Hartree/Bohr
    }
    initial_wavefunctions = matrix_diagon
    ksampling{
        mesh{
            nx = 4
            ny = 4
            nz = 4
        }
    }
    scf_convergence{
        delta_total_energy = 1e-12 Hartree
        succession = 3
    }
}
...
...
structure{
    unit_cell_type = primitive
    unit_cell{
        a_vector = 0.0000000000        5.1300000000        5.1300000000
        b_vector = 5.1300000000        0.0000000000        5.1300000000
        c_vector = 5.1300000000        5.1300000000        0.0000000000
    }
    atom_list{
        atoms{
            #tag element rx ry rz mobile
             Si 0.130 0.130 0.130 yes
             Si -0.130 -0.130 -0.130 yes
        }
    }
    element_list{
        #tag element atomicnumber
         Si 14
    }
}
...
...
structure_evolution{
    method = velocity_verlet
    dt = 100
}
...
...

この入力は,シリコン結晶の入力を少し変更したものとなっています。 atomsでは,各原子の“mobile”変数を“yes”と設定しています。ここを“no”あるいは“0”と設定すると,その原子は 分子動力学シミュレーションを行っても動くことはありません。さらに座標値をあえて安定でない値にしています。 具体的には,Si結晶の二つの原子を(111)方向にお互いから離れるように少しだけずらしています。

手法の選択

structure_evolutionブロックでは,“method”変数を“velocity_verlet”としています。この選択によって小正準集合の分子動力学シミュレーションを行うことができます。また,各ステップでの更新量(変数dt)を,原子単位で“100”としてい 上で述べたようにこの値は \(2.418 \times 10^{-15}\) sに相当します。

初期速度の与え方

ここまで説明したサンプルの入力を利用すると,原子の初期速度は全て0と設定されます。 原子に初期速度を与える場合,下記のような入力を準備してください。

structure_evolution{
    method = velocity_verlet
    dt  = 100
    temperature_control{
        thermostat{
            #tag temp
                 300
        }
    }
}

ここで,“temp”変数で初期の温度をケルビン単位で設定します。原子の初期速度は,この温度になるように, かつ正規乱数に従って,全運動量が0になるように設定されます。

原子ごとに異なる初期温度を設定することも可能です。この場合,まず下記のような入力を作成します。

structure_evolution{
    method = velocity_verlet
    dt  = 100
    temperature_control{
        thermostat{!#tag temp
                         300
                         500
                         700
        }
    }
}

次にatomsの各原子に,“thermo_group”という変数を設定します。

structure{
    ...
    atom_list{
        atoms{
        !#tag rx ry rz  element mobile weight  thermo_group
   0.1159672611      0.1235205209      0.1215156388    Si   1   1  1
  -0.1329067626     -0.1264216714     -0.1225370484    Si   1   1  2
   0.1273740089      0.6305999369      0.6247606249    Si   1   1  3
            ...
            ...
        }
    }
    ...
}

この例では一番目の原子が300Kに,二番目の原子が500Kに, 三番目の原子が700Kになるよう初期速度が設定されます。

計算結果の例

この計算例の計算結果の内部エネルギー,運動エネルギー,全エネルギーを 図 5.47 に示します。

../_images/image98.svg

内部エネルギー,運動エネルギー,全エネルギーと時間の関係。

使用方法:温度一定の分子動力学シミュレーション

温度一定の分子動力学シミュレーションの入力パラメータ例です。

計算例題は、 samples/dynamics/molecular_dynamics/NVT です。

熱浴の設定

structure_evolutionブロックにtemperature_controlブロックを指定します。

structure_evolution{
    method = temperature_control
    dt  = 50.0
    temperature_control{
        thermostat{
            #tag temp  qmass tdamp
                 300   5000 10000
        }
    }
}

上記の入力例では,まず“method”変数をtemperature_controlとしています。この変数によって温度制御を行うように指定します。 ついで,“dt”変数を設定しています。これは,時間刻みの指定です。原子単位で入力します。例で示されている50.0という値は, 約1.2fsに相当します。

さらに,temperature_controlブロックで熱浴の設定を行っています。 “thermostat”ブロックで各熱浴のパラメタを設定します。“temp”パラメタによってその熱浴の目的とする温度(ケルビン単位), “qmass”パラメターによって熱浴の質量(原子単位)を設定します。“qmass”パラメータによって質量を直接指定するのではなく、“tdamp”パラメーター(時間の単位)によって熱浴の周期を指定し間接的に熱浴の質量を設定することも可能です。“qmass”と“tdamp”が両方設定されている場合,“qmass”が優先されます。また、いずれの指定もない場合50×dtの周期が実現するように熱浴の質量が設定されます。

熱浴の割り当て

structure ブロックの,atoms ブロックを設定する必要があります。設定例を以下に記します。

structure{
     ...
    atom_list{
        num_atoms = 8
        cooordinate_system = internal
        atoms{
        !#tag rx ry rz  element mobile weight  thermo_group
   0.1159672611      0.1235205209      0.1215156388    Si   1   1  1
  -0.1329067626     -0.1264216714     -0.1225370484    Si   1   1  1
   0.1273740089      0.6305999369      0.6247606249    Si   1   1  1
  -0.1152089939     -0.6164829779     -0.6221565128    Si   1   1  1
   0.6299472943      0.1341313888      0.6253193197    Si   1   1  1
  -0.6305720382     -0.1290073650     -0.6187967685    Si   1   1  1
   0.6151271805      0.6206113965      0.1333834419    Si   1   1  1
  -0.6276524003     -0.6268549639     -0.1175099372    Si   1   1  1
        }
    }
    ...
}

各原子に“thermo_group”パラメータを割り振っています。このパラメータに 熱浴の識別番号を設定します。なお,熱浴の識別番号は熱浴の定義順に割り振られます。また,他の属性値と同様,“#default” タグを利用することによってデフォルト値を設定することも可能です。ここの例では全ての原子に同じthermo_groupを設定していますが, 各原子が異なる熱浴に関連付けられていても問題ありません。

熱浴の“段数”の設定 (バージョン2019.01以上)

通常の熱浴をさらに熱浴で制御する,という計算手法がNosé-Hoover chain法です [Glenn92] 何段にもわたって再帰的に適用することも可能となっています。Nosé-Hoover chain法を用いることによって,より少ないステップ数で熱浴が平衡状態に至ることなどが期待できます。利用するためには,以下のようにtemperature_controlブロックにおいてnum_chain変数を使います。

structure_evolution{
    method = temperature_control
    dt  = 50.0
    temperature_control{
      num_chain = 5
      thermostat{
          #tag temp  qmass tdamp
               300   5000 10000
      }
    }
}

理論的には各段数における熱浴の質量を別途設定することも可能ですが,最初の熱浴の質量以外結果にほとんど影響を与えないため,のこりの熱浴の質量としてはPHASE/0が最初の熱浴の質量から自動的に算出した値が採用されるようになっています。

num_chainのデフォルト値は1, すなわち通常のNosé-Hoover法です。

“温度プロファイル”の設定 (バージョン2019.01以上)

シミュレーションの進行とともに入力パラメーターファイルの指定に応じてターゲット温度を変化させる分子動力学シミュレーションを行うことができます。変化のさせ方を“温度プロファイル”と呼びます。温度プロファイル機能を利用するためには,まずtemperature_controlブロックにおいてsw_temperature_profile = onとします。さらに,thermostatブロックにおいて温度プロファイルの設定を行います。温度プロファイル機能を利用していない場合thermostatブロックにおいて定義するテーブルの行は一つの熱浴を表しますが,利用している場合はテーブルの行は一つのプロファイルを表すことに注意が必要です。温度プロファイル機能を利用している場合にさらに複数の熱浴を定義するためには,属性値no (id, thermo_groupも可)によって切り分けます。具体的には,たとえば下記に示すようになります。

structure_evolution{
  ...
  ...
  temperature_control{
    sw_temperature_profile = on
    thermostat{
      #tag no tempi tempf till_n tdamp
          1 8000 8000 3000 5000
          1 8000 300 7000 5000
          2 400 500 1000 5000
    }
  }
}

上述のように属性値noによって熱浴を指定します。tempi, tempfはそれぞれ初期温度および終温度に対応する属性値です。till_nは,この温度プロファイルが適用されるステップ数です。属性値tdampによって熱浴の緩和時間を指定します。tdampのかわりに属性値qmassを使って熱浴の質量を直接指定することも可能です。以上より,この例ではthermostatテーブルの各行は次のように解釈されます。

  • 1行目:1つ目の熱浴の温度を8000Kで3000ステップ適用する。

  • 2行目:1つ目の熱浴の温度を8000Kから300Kまで7000ステップかけて降温する (プロファイル終了時点で総MDステップ数10000)。以降は300Kのまま進行する。

  • 3行目:2つ目の熱浴の温度を400Kから500Kまで1000ステップかけて昇温する。以降は500Kのまま進行する。

なお,最後のプロファイルのtill_n以降は,そのプロファイルのtempfの温度が使い続けられる仕様になっています。

ランジュバン熱浴の設定(バージョン2021.02以上)

ランジュバン熱浴の運動方程式は次の様に定義することができます。

()\[\begin{split}\dot{r}_i = p_i/m_i \\ \dot{p}_i = F_i - \gamma_i p_i + \sqrt{\frac{2m_i\gamma_ik_B T}{\Delta t}} R\left(t\right) \\ \left< R\left(t\right) \right> = 0, \left< R\left(t\right) R\left(t^\prime \right)\right> = \delta \left(t-t^\prime \right)\end{split}\]

\(\gamma_i\) は時間の逆数の単位を持つ量で,摩擦係数のような役割を果たします。これは,能勢フーバー熱浴の緩和時間(もしくは質量)に対応するものであり,デフォルト値が設定されていますが入力パラメーターファイルを通じて指定することもできます。\(R\left(t\right)\) はランダム力であり,原子に加えることによって最終的には所定の温度へ至るように動作します。 \(\left<R\left(r\right)\right>\) は分散1, 中心0の正規乱数によって評価します。

ランジュバン熱浴を利用するには,たとえば以下のように設定します。

structure_evolution{
  method = velocity_verlet
  temperature_control{
    method = langevin
    thermostat{
      #tag temp tdamp
               300   5000
    }
  }
}

structure_evolutionにおいてmethodをvelocity_verletとし,temperature_controlブロックのmethodをlangevinに設定するとランジュバン熱浴を利用することができます。ランジュバン熱浴の摩擦係数 \(\gamma_i\) はこれまでの熱浴の緩和時間指定と同じ方法で行います。すなわち,緩和時間を直接もしくは熱浴の質量を通じて間接的に指定します。 熱浴に属する原子の指定方法などはこれまでの方法と同様です。

ランジュバン熱浴は統計力学を直接適用するような熱浴であり,能勢フーバー法と比較して統計力学的により尤もらしい集団の計算が実現できる場合があります。 例として,ダイヤモンドの結晶を用いてランジュバン熱浴と能勢フーバー熱浴双方で計算を行い,結果を比較してみました。

用いた系はダイアモンド64原子の系です。カットオフエネルギーは25 Rydberg, \(\textbf{k}\) 点サンプリングは \(\Gamma\) 点のみとしました。時間刻みは1 fsとし,緩和時間としてはデフォルト値(時間刻みの50倍)を採用しました。合計1万ステップの分子動力学シミュレーションを実施しました。

図 5.48 にシミュレーションの結果得られた運動エネルギーの履歴を示しました。この図から明らかなように,ランジュバン熱浴によるシミュレーションではターゲット温度に相当する運動エネルギー(図中の黒い水平線)に はやい段階で達しています。また,能勢フーバー熱浴の結果と比較すると運動エネルギーの振幅も小さい傾向です。結晶の計算の場合ランダムな作用による減衰が少ないため決定論的な能勢フーバー熱浴を用いると人工的な振動が減衰しにくい傾向がありますが, ランジュバン熱浴の場合はそもそもがランダム力を付与しているため,「ランダムな作用による減衰」が発生しやすい影響が現れた結果と考えられます。

../_images/langevin_image1.svg

運動エネルギーの履歴。赤:ランジュバン熱浴,緑:能勢フーバー熱浴。

図 5.49 には運動エネルギーの分布を示します。カノニカルアンサンブルの場合,運動エネルギーはマクスウェルボルツマン分布に従うはずです。運動エネルギーのマクスウェルボルツマン分布は以下のように記述することができます。

()\[f\left(K\right) = \frac{2}{\pi} \frac{1}{\left(k_B T \right)^\frac{3}{2}} \sqrt{K} \exp{\left(-\frac{K}{k_B T}\right)}\]

比較のため, (5.36) 式の分布も 図 5.49 にプロットしました。一見して明らかなように,ランジュバン熱浴の方がマクスウェルボルツマン分布によく従う分布となっています。 これもやはりランダム力を取り入れているがゆえのモードの減衰が要因と考えられます。

../_images/langevin_image2.svg

運動エネルギー分布。赤丸:ランジュバン熱浴,緑丸:能勢フーバー熱浴,黒線:マクスウェルボルツマン分布

なお,現バージョンにおいては後述の圧力一定の分子動力学シミュレーションにランジュバン熱浴を用いることはできません。

使用方法:温度・圧力一定の分子動力学シミュレーション

PHASE/0は, [Souza97] , [Hernandez01] などの文献において解説されている手法によって圧力一定の分子動力学シミュレーションを行うことができます。

入力パラメーターファイルの書き方

まず,温度一定の分子動力学シミュレーションと同様の手続きで入力パラメーターファイルを作成します。通常の温度一定の場合との違いは,圧力も一定であること,また圧力浴の設定が必要なことの2点です。以下に例を示します。

structure_evolution{
  method = temperature_pressure_control
  dt = 50.0
  temperature_control{
    num_thermostat = 1
    set_initial_velocity = on
    !!method = velocity_scaling
    thermostat{
      #tag temp qmass
            300 4000
    }
  }
  pressure_control{
    pressure = 0.0
    mass_baro = 1
    m11 = on
    m22 = on
    m33 = on
    m12 = on
    m13 = on
    m23 = on
  }
}

設定方法の詳細は,下記の通り。

  • MD手法の選び方 | methodにtemperature_pressure_controlもしくはpressure_temperature_controlを指定すると温度-圧力一定の分子動力学シミュレーション,すなわちNPT分子動力学シミュレーションを実行することができます。pressure_controlとすると圧力一定の分子動力学シミュレーション,すなわちNPH分子動力学シミュレーションを実行することができます。
  • 温度制御の設定 | 温度制御は,NVTの場合と同様temperature_controlブロックにおいて行います。ただし,NVTシミュレーションでは複数の熱浴を定義することが可能であるのに対し,NPTシミュレーションの場合に定義できる熱浴は一つのみである点に注意してください。また,上記の例ではコメントアウトされていますが,method = velocity_scalingと指定することによって速度スケーリングによる温度制御を利用することもできます。Nosé-Poincaréの熱浴で温度制御を行うことが理想的ですが,安定に運動方程式が時間発展できない場合は速度スケーリングを利用してください。指定がなければNosé-Poincaré熱浴が採用されます。
  • 圧力制御の設定 | 圧力制御は,pressure_controlブロックを作成し,その下で設定します。pressureに目的の圧力を指定します。圧力の単位のデフォルト値は原子単位ですが,GPaなどを利用することもできます。mass_baroには熱浴の質量を指定します。デフォルト値は1であり,単位の指定はできません。目安としては,熱浴の周期(単位胞の体積変化の周期)が100 MDステップ相当程度の時間になるように決めるとよいでしょう。
  • 計量テンソルに対する拘束条件の設定 | mxyパラメーターによって,計量テンソルの変化に対して拘束条件を課すことができます。たとえばm11 = offとすると11成分(a軸・a軸)が変化しなくなり,m12 = onならば12成分(a軸・b軸)および21成分(b軸・a軸)が変化しなくなります。このパラメーターのデフォルト値はすべてonなので,特に指定がなければ計量テンソルのすべての成分が変化します。

NPTシミュレーション結果で得られる主な出力ファイルは,エネルギーなどの履歴を記録したnfefn.dataファイル,座標データの履歴を記録したnfdynm.dataファイル,計量テンソルの履歴を記録したnfmetric.dataファイル,そして格子定数の履歴を記録したnflatconst.dataファイルである。それぞれについて説明する。

nfefn.dataファイル

エネルギーなどの履歴を記録したファイルであり,構造最適化やNVTのMDシミュレーションの場合も得られます。NPTシミュレーションの場合,以下のような出力が得られます。なお,紙幅の都合で1行を改行して表示しています。

iter_ion, iter_total, etotal, ekina, econst, pressure
1 13 -31.8045273788 300.0000000000 0.0000000000 0.0209863336 -0.0000504658
2 25 -31.8045248143 301.4440100456 0.0000586853 0.0211839341 -0.0000573770
3 37 -31.8043935680 299.1271868912 0.0001085792 0.0214781523 -0.0000574428
4 49 -31.8041402859 293.2700519011 0.0001476826 0.0218951630 -0.0000570036
5 61 -31.8037768714 284.3267516440 0.0001733906 0.0224437768 -0.0000560637
6 73 -31.8033194235 272.9264725662 0.0001832795 0.0231125156 -0.0000546420
7 85 -31.8027870051 259.8057152864 0.0001752150 0.0238698884 -0.0000527691
8 97 -31.8022003504 245.7426304561 0.0001475268 0.0246685850 -0.0000504870
9 109 -31.8015806171 231.5028071591 0.0000992076 0.0270606007 -0.0000478460
10 121 -31.8008921515 217.8018629938 0.0000862100 0.0294678421 -0.0000449028
11 133 -31.8002640917 205.3149756742 0.0000004873 0.0315593847 -0.0000420615
...
...

一行が1タイムステップのデータに相当します。1列目が原子(単位胞)の更新回数,2列目がSCF計算の総更新回数,3列目が系のエネルギー(原子単位),4列目が瞬間的な温度(ケルビン単位),5列目がハミルトニアン(原子単位),6列目が原子に働く力の最大値(原子単位),7列目が瞬間的な圧力です(すべて原子単位)。ハミルトニアンは理想的には保存するはずですが,第一原理計算においてエネルギーに対する歪みテンソルの厳密な微分を求めることは困難なため,保存がよくない傾向にある点に注意してください。

nfdynm.dataファイル

基本的なファイルフォーマットは,通常のnfdynm.dataファイルと同じです。ただし,NPTシミュレーションの場合毎ステップ格子定数が変化するので,格子定数などが記述されるヘッダー部が毎ステップ書き込まれる点に違いがあります。

nfmetric.dataファイル

各ステップにおける計量テンソルが記録されるファイルです。典型的な内容は下記の通り。

2
105.8775163849 0.0064685429 -0.0017420207
0.0064685429 105.8774556260 0.0027757249
-0.0017420207 0.0027757249 105.8776712432
3
105.8648537936 0.0196786348 -0.0055233812
0.0196786348 105.8644796107 0.0091955626
-0.0055233812 0.0091955626 105.8652454083
4
105.8462013539 0.0389223098 -0.0123434716
0.0389223098 105.8452005929 0.0190717120
-0.0123434716 0.0190717120 105.8468474757

まずステップ数を表わす整数値が記録され,ついで3×3のテンソルが記録されます。なお,半ステップずれた数値解法を採用しているため,記録は2ステップ目からとなります。

nflatconst.dataファイル

各ステップにおける格子定数が記録されるファイルです。典型的な内容は下記の通り。なお,紙幅の都合で1行を改行して表示している

2 10.2896800915 10.2896771391 10.2896876164 89.9984979129 90.0009426965 89.9964995371 1089.4462540511
3 10.2890647677 10.2890465841 10.2890837983 89.9950232125 90.0029893382 89.9893495841 1089.2504025643
4 10.2881583072 10.2881096705 10.2881897084 89.9896762432 90.0066816443 89.9789308011 1088.9605527024
5 10.2869678503 10.2868710406 10.2870094050 89.9825817981 90.0125180935 89.9656805224 1088.5783846407
6 10.2855038285 10.2853392827 10.2855504459 89.9738915620 90.0209385359 89.9500901373 1088.1067192172
7 10.2837804786 10.2835276570 10.2838245727 89.9637869404 90.0323061920 89.9327006608 1087.5497224800
8 10.2818163114 10.2814547573 10.2818483397 89.9524812878 90.0468878988 89.9140944258 1086.9130916478
9 10.2796344063 10.2791449593 10.2796435437 89.9402204612 90.0648344761 89.8948814720 1086.2041781879

一行が1タイムステップのデータに相当します。1列目は単位胞の更新回数,2列目から7列目がそれぞれ格子定数, 最後の8列目が単位胞の体積です。長さの単位はBohr単位,角度の単位は度,体積の単位はBohr3です。

分子動力学シミュレーションにまつわるその他の設定

速度スケーリングによる温度制御の分子動力学シミュレーション

Nos\(\acute{e}\)-Hooverの熱浴ではなく,原子の速度を温度が合うようにスケールし直すことによって温度を制御することも可能です。このような計算は,temperature_controlブロックの下のmethod変数をvelocity_scalingとすることによって実現することができます。

structure_evolution{
  ...
  temperature_control{
    method = velocity_scaling
    ...
  }
}

速度は,目的の温度に合うよう毎MDステップスケールされます。

全運動量がゼロになるよう速度をシフトする方法

分子動力学シミュレーションにおいて理論上全運動量は保存します(ゼロになる)が,実際は数値誤差によりゼロとはならず,系全体が並進してしまう場合があります。これを防ぐには,以下のように変数sw_shift_velocitiesの値をonとし,MD計算中全運動量がゼロとなるようにします。

structure_evolution{
  ...
  temperature_control{
    ...
    sw_shift_velocities = on
  }
}

原子の初期速度を手動で指定する方法

原子の初期速度は通常ランダムな正規分布が得られ,対応する温度が設定した温度で,全運動量が0になるように自動的に決定されますが,手動で各原子に割り振ることも可能です。このような指定は,原子配置のvx, vy, vz属性値によって行います。また,設定した速度を上書きしてしまわないようtemperature_controlブロックの下においてsw_read_velocities = onを設定する必要があります。たとえば,以下のように設定します。

structure{
     ...
    atom_list{
        atoms{
        !#tag rx ry rz  element mobile thermo_group vx vy vz
            0.11    0.12    0.11    Si   1 1 0.001 0.0014 0.0008
           -0.13   -0.13   -0.14    Si   1 1 -0.001 -0.002 0.0001
            0.12    0.63    0.62    Si   1 1 0.0003 -0.0005 -0.00028
            ...
            ...
        }
    }
    ...
}
...
structure_evolution{
  ...
  temperature_control{
    sw_read_velocities = on
  }
}

なお,速度の単位のデフォルト値は原子単位です。

また,この方法で初期速度を指定した場合でも,デフォルトの設定では温度に合うよう,また全運動量が0になるよう速度はシフト・スケールされます。この動作を抑制したい場合,以下の要領で変数set_initial_velocityの値をoffとします。

structure_evolution{
  ...
  temperature_control{
    sw_read_velocities = on
    set_initial_velocity = off
  }
}

原子を領域に閉じ込める方法 (バージョン2019.01以上)

矩形領域および円筒領域に閉じ込めることができます。このような計算機能は,狭い領域に閉じ込められた系の振る舞いを調べたり,ESM法やdipole補正法などアルゴリズムの都合上真空層が必須である機能において真空層に原子が入ってこられなくする目的で利用することができます。

入力の記述

領域は,structureブロックのregionxブロックにおいて定義することができます。ここでxは領域の識別番号です。

structure
  region1{
    region_group = 1
    type = cylinder
    radius = 3.5 angstrom !半径
    cylx = 5 angstrom
    cyly = 5 angstrom ! 円筒の中心位置
    orientation = 3 ! 円筒の向き
    cylzmin = -1000 ! 円筒の下端
    cylzmax = 1000 ! 円筒の上端
    sw_tally = on ! エネルギーを加える
    eps = 0.001 !ポテンシャルの深さ
    sigma = 1.5 !距離のスケール
  }
}

xは1から始まります。xを2, 3, …とすることによって任意の数の領域を定義することができます。regionxブロックにおいて以下の変数を設定することによって領域を定義します。

変数名

説明

region_group

“領域のグループ”を整数で指定します。この値が同じ 領域はひとかたまりの“グループ”として扱われます。後述 の原子配置テーブルの属性値region_groupにこの数値を指 定することによって領域グループと原子を紐づけます。デ フォルト値はregionxブロックのxの値です。

type

領域の種類を指定します。

cylinder (円筒型)もしくはbox (直方体型)とします。デフォルト値はboxです。

radius

円筒の半径を指定します。デフォルト値は6.5 b ohrです。type=cylinderの場合のみ意味のある設定です。

cylx

円筒の中心位 置のx座標を指定します。デフォルト値はセルの境 界です。type=cylinderの場合のみ意味のある設定です。

cyly

円筒の中心位 置のy座標を指定します。デフォルト値はセルの境 界です。type=cylinderの場合のみ意味のある設定です。

orientation

円 筒の向きを指定します。1の場合x方向,2の場合 y 方向,3の場合z方向に向きます。デフォルト値 は3です。type=cylinderの場合のみ意味のある設定です。

cylzmin

円筒の長さ 方向の下限を指定します。デフォルト値は-10 10 です。type=cylinderの場合のみ意味のある設定です。

cylzmax

円筒の長 さ方向の上限を指定します。デフォルト値は10 10 です。type=cylinderの場合のみ意味のある設定です。

xmin

直方体のx 方向の下限値を指定します。 デフォルト値は-10 10 です。type=boxの場合のみ意味のある設定です。

xmax

直方体のx方向の上限値を指定します。 デフォルト値は10 10 です。type=boxの場合のみ意味のある設定です。

ymin

直方体のy 方向の下限値を指定します。 デフォルト値は-1010 です。type=boxの場合のみ意味のある設定です。

ymax

直方体のy方向の上限値を指定します。 デフォルト値は1010 です。type=boxの場合のみ意味のある設定です。

zmin

直方体のz 方向の下限値を指定します。デフォルト値は デフォルト値は-1010 です。type=boxの場合のみ意味のある設定です。

zmax

直方体のz方向の上限値を指定します。デフォルト値は デフォルト値は1010 です。type=boxの場合のみ意味のある設定です。

sw_tally

領域に起因するエネルギ ーを全エネルギーに加える場合この値をonに設定します。

eps

ポテンシャル\(\varepsilon\left( \frac{ \sigma}{r} \right)^{12}\)\(\varepsilon\)の 値をエネルギーの単位で指定します。デフォルト値は1e-3 hartreeです。

sigma

ポテンシャル\(\varepsilon \left( \frac{\sigma}{r} \right)^{12}\)\(\sigma\) の値を長さの単位で指定します。 デフォルト値は1 bohrです。

定義した領域は,そのregion_idを紐づけたい原子のregion_idに指定することによって割り当てることができます。たとえば,以下のように設定します。

structure{
  …
  atom_list{
    coordinate_system = cartesian
    atoms{
      #default mobile=on,thermo_group=1
      #units angstrom
      #tag element rx ry rz region_group
            H 4.231707 4.904619 6.374683 1
            H 5.716594 4.994127 6.011627 1
            O 5.118193 4.883964 6.766158 1
            H 4.167342 5.876768 8.210465 2
            H 5.481543 5.259672 8.697061 2
            O 4.627457 5.590603 9.014168 2
            …
            …
    }
  }
}

この例では,1番目,2番目,3番目の原子はregion_group=1のすべての領域に紐づけられ,4番目,5番目,6番目の原子はregion_group=2のすべての領域に紐づけられます。

計算の実行

計算は,通常の分子動力学シミュレーションと同じように実行することができます。入力パラメーターファイルが読み込まれると以下のように読み込んだ領域の情報が報告されます。

!** region statistics
!** num_regions = 1
!** status for region no 1
!** region type : CYLINDER
!** orientation : 3 (1->x, 2->y, 3->z)
!** radius : 6.6140409725
!** cylx,cyly : 9.4486299607 9.4486299607
!*cylzmin,cylzmax: -1.797693134862316E+308 1.797693134862316E+308
!** sigma, epsilon : 1.0000000000 0.0010000000
!** tally : F
!** n target atoms : 36
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
21 22 23 24 25 26 27 28 29 30 31 32 33 37 38 39

特に,n target atoms以降領域に割り当てられている原子のリストが出力されるので,想定通り領域が設定できているかどうかを確認することができます。

使用における注意点

  • 分子動力学シミュレーション機能に特別な制約はありません。ウルトラソフトおよびPAW擬ポテンシャル,並列計算,継続計算に対応しています。 ただし,分子動力学シミュレーションにおいては原子の位置は熱揺らぎによって激しく振動するので,系が厳密な意味での対称性を持つことはありません。したがって,対称性関連の設定は施さないでください。

  • 原子の質量は,通常の構造最適化においては直接系の性質を左右するものではありませんが,分子動力学シミュレーションの場合は意味のある量です。したがって,本機能を利用する場合は元素の正しい質量を指定する必要があります。PHASEにおける標準の質量の単位は,原子単位です。たとえば,陽子の質量は原子単位で1822.877333 です。

  • 運動エネルギーはハートリー単位で記述されますが, 運動エネルギーと温度との間には\(E_{\text{kin}} = \frac{3}{2} \times N_{\text{atom}} \times k_{B}T\) という関係があります。ここで運動エネルギーを\(E_{\text{kin}}\),原子数を\(N_{\text{atom}}\),ボルツマン定数を\(k_{B}\), 瞬間的な温度を\(T\)と記述しました。従って運動エネルギーから系の温度を知りたい場合,まず運動エネルギーを原子数で割り, \(3.1578 \times 10^{5}\)という値を掛けて(ハートリー単位から\(k_{B}T\)単位への変換),さらに\(\frac{3}{2}\)で割って下さい。

  • 原子座標の更新回数は入力で指定した時間刻みの値(structure_evolutionブロック内のdtという変数で指定)を掛けること によって実時間での経過が分かります。時間の単位は入力で指定することが可能ですが,原子単位系を利用した場合(デフォルト) \(2.418 \times 10^{- 17}\)という値を掛ければ「秒」に変換することが可能です。たとえば100 a.u.という時間は2.418 fsに相当します。

  • 温度一定の分子動力学シミュレーションにおける、熱浴の質量\(Q\)の値について注意点を挙げます。 小さすぎる\(Q\)を採用すると,系のダイナミックスに熱浴に起因する 人為的なモードが生じてしまい,また場合によっては計算が破綻してしまいます。他方大きすぎる\(Q\)を採用すると,系が熱平衡に達するのに 多くのステップ数を必要とするようになってしまいます。

    \(Q\)の値は,系の特徴的な振動の周期と熱浴の振動の周期がおおよそ同等か,熱浴の方が長くなるように選ぶとよいとされています。 熱浴の振動の周期と系の振動の周期の間の関係は,大雑把には次の式で評価できます [Nose91]

    \(\tau = \frac{2\pi}{\omega} = 2\pi\left( \frac{Q}{2gk_{B}T} \right)^{1/2}\)

    ここで\(\tau ,\omega\)はそれぞれ系の周期と周波数,\(g\)は系の自由度(3\(\times\)その熱浴に関連付けられている原子の数),\(k_{B}\)は ボルツマン定数,\(T\)がその熱浴の温度です。例えば,\(\tau\)を0.05 ps,原子の数を8,温度を300 Kとして上式で\(Q\)の値を見積もると,原子単位でおおよそ4600程度となります。PHASE/0では、\(Q\)の値を直接指定することも、周期を介して間接的に指定することも可能となっています。質量の指定も周期の指定もない場合、周期が時間きざみdtのおおよそ50倍となる\(Q\)の値が採用されます。

参考文献

Glenn92

Glenn J. Martyna and Michael L. Klein, and Mark Tuckerman, “Nosé-Hoover chains: The canonical ensemble via continuous dynamics”, Journal of Chemical Physics 97 15 (1992).

Souza97

Ivo Souza and JoséLuís Martins, Phys. Rev. B 55 (1997) pp 8733-8742.

Hernandez01
  1. Hernández, Journal of Chemical Physics, 115 (2001) pp. 10282-10290.

Nose91
  1. Nos\(\acute{e}\),Progress of Theoretical Physics Supplement No 103,1991,pp.1-46.