4. 入力パラメータファイル：nfinp.data（詳細版）

4.1. 入力パラメータファイルの形式(F_INP ファイル)

入力パラメータファイルnfinp.dataは、どのようなモデル（原子配置など）に対し、どのような条件で計算するかという情報を記述するファイルです。

デフォルトのファイル名はnfinp.data ですが、file_names.dataにおいて、F_INPキーワードを使って自由に名前を指定できます。例えば、計算する系に関連した名前をつけることも可能です。

4.1.1. パラメータ設定形式

このファイルは、タグ名と中括弧{}で囲まれたブロックの階層構造で記述します。計算パラメータの設定は、タグ形式になっており、各タグに、結晶構造、計算精度、計算の制御などの情報を記入します。以下に, 入力パラメータファイルの記述方法を簡単に説明します。

関連のある入力データはまとめて一つの「ブロック」内に記述します。ブロックは, ブロック名{ ... } という形式で階層的に記述します。たとえば, 以下のようになります。基本的に計算パラメータは、タグ=値の形で指定します。ただし、パラメータごとに指定形式が異なります。各計算パラメータの指定方法の説明を参照ください。

Upper_block{
  Lower_block{
    ...
    tag_keyword = value
  }
}

入力パラメータファイルの作成・編集において、以下の点に注意ください。

  • 同じ階層に同じ名前のブロックを記述することはできません。

  • ブロック名の大文字・小文字を区別することはありません。

  • ブロック名が間違っている場合には、そのブロック全体が記述されていないものとみなされます（無視されます）。その場合は、デフォルト値が採用されます。エラーメッセージは表示されません。

  • 変数は、改行で区切るほかカンマ(,) で区切ることもできます。

  • 文字列型の変数に空白文字を含めたい場合、半角の2 重引用符(”)で囲みます。

  • 全角文字は使用しないでください。

最上位のブロックは、以下のものがあります。

control ブロック

全体的な計算条件の設定

accracy ブロック

計算精度の設定

structure ブロック

原子構造の設定

wavefunction_solver ブロック

波動関数ソルバーの設定

charge_mixing ブロック

電荷密度混合法の設定

structure_evolution ブロック

構造最適化、分子動力学法計算の設定

postproccesingブロック

後処理の設定

printlevelブロック

ログ出力の設定

4.1.2. 単位の指定

PHASE の入力ファイルのデフォルトの単位は原子単位ですが, 単位を明示的に指定することも可能です。表 4.1 の単位を利用することができます(デフォルトの単位は太字で表示されています)。

表 4.1 PHASEで利用可能な単位

長さ

bohr, angstrom, nm

エネルギー

hartree, eV, rydberg

時間

au_time, fs, ps, ns, s, sec, min, hour, day

速度

bohr/au_time, bohr/fs, angstrom/fs, angstrom/au_time, nm/fs, nm/au_time

力

hartree/bohr, hartree/angstrom, hartree/nm, eV/angstrom, eV/bohr, ev/nm, rydberg/bohr, rydberg/angstrom, rydberg/nm

圧力

hartree/bohr3, hartree/angstrom3, hartree/nm3, eV/angstrom3, eV/bohr3, eV/nm3, rydberg/angstrom3, rydberg/bohr3, rydberg/nm3

質量

au_mass, atomic_mass

単位は, 実数型のデータに直接指定する方法だけでなく, ブロック単位のデフォルト値を指定することもできます。ブロック単位でデフォルトの単位を指定するには, 次のように記述します。

block{
   #units angstrom
   ...
   ...
}

この例では, block の長さの単位のデフォルトがÅ単位になります。 複数指定する場合（長さ、エネルギーなど）, 空白で区切って指定してください。

4.1.3. 表形式データを含むブロック

ブロックは「表形式データ」を持つ場合があります。このようなブロックは対応する表形式データのみ記述することができます。 表形式データは2次元の行列形式でデータを指定することのできる形式であり、たとえば各原子の座標値のように複数の属性値をもつデータを任意の数定義する必要があるようなデータに対して採用されています。 その内容は、たとえば以下のようになります。

tabular_data{
     #units  angstrom
     #defaullt column1=1, column2=off
     #tag column1 column2 column3
          data11  data12  data13
          data21  data22  data23
          ...
}

ブロック定義直下に units によってブロック共通の単位を設定することができるのは通常のブロックと同じです。通常のブロックと違い単位を指定する方法はこれしかありません。#default によって表形式データのデフォルト値を設定することができます。この例では column1 というカラムには 1 というデフォルト値が、 column2 というカラムには off というデフォルト値が設定されます。値として * を指定するとデフォルト値が採用されます。最後のカラムの場合、何も指定がないとやはりデフォルト値が採用されます。#tag によってカラムを定義します。スペース区切りで利用したいカラム名を指定します。#tag の次の行からが表形式データ本体となります。#tag 行で定義したカラム順に対応するデータを指定します。

4.1.4. コメント

!または//ではじまる行は, コメント扱いとなります。

block{
 !  comment
 !  tag_keyword = value1 コメント
 // tag_keyword = value2 コメント
 tag_keyword = value3
}

ただし、!#と、!のあとに#が続く場合はコメントとはみなされないので注意してください。

4.1.5. 二項演算子を用いた実数指定（バージョン2020.01以降）

バージョン2020.01以降、二項演算子を用いて実数を指定することができるようになりました。 この機能を活用することによって、たとえば0.33333... を1/3と記述することなどができます。 +, -, /, * (もしくはx)のいずれかの文字列を含むような数値を検出した場合にそれぞれ足し算、引き算、除算、掛け算を行います。 ただし、+, -には下記の制約があります。

  • 先頭が+もしくは-の場合は使えない

  • e, dを含む数値指定(たとえば1e-3, 1d+4)には利用できない

また、一つの実数指定に複数の二項演算子を適用することはできません。

4.1.6. 入力パラメータファイル例

Si ダイヤモンド結晶（2原子）の電子状態計算を行う場合の基本的な計算条件を記述した入力ファイル例です。

control{
  condition = initial
  cpumax = 86400 sec
  max_iteration = 10000
}
accuracy{
  cutoff_wf = 25.0 rydberg
  cutoff_cd = 100.0 rydberg
  num_bands = 8
  ksampling{
    method = monk
    mesh{
      nx = 10
      ny = 10
      nz = 10
    }
  }
  initial_wavefunctions = atomic_orbitals
  initial_charge_density = atomic_charge_density
  scf_convergence{
    delta_total_energy = 1e-10
    succession = 3
  }
  force_convergence{
    max_force = 0.001 hartree/bohr
  }
}
structure{
  element_list{
  #tag element atomicnumber
        Si 14
  }
  unit_cell{
    #units angstrom
    a_vector = 0 2.732299538 2.732299538
    b_vector = 2.732299538 0 2.732299538
    c_vector = 2.732299538 2.732299538 0
  }
  unit_cell_type = bravais
  atom_list{
    atoms{
    #tag element rx ry rz mobile
          Si 0.125 0.125 0.125 0
          Si -0.125 -0.125 -0.125 0
    }
    coordinate_system = internal
  }
}
wavefunction_solver{
  solvers{
    #tag sol till_n prec cmix submat
        davidson 1 on 1 on
        rmm3 -1 on 1 on
  }
  rmm{
    edelta_change_to_rmm=5e-5
  }
}
charge_mixing{
  mixing_methods{
    #tag no method rmxs rmxe istr prec nbmix
      1 pulay 0.40 0.40 3 on 15
  }
}
Postprocessing{
  dos{
    sw_dos = ON
    deltaE = 1.e-4 hartree
  }
  charge{
    sw_charge_rspace = ON
    filetype = cube !{cube|density_only}
    title = "This is a title line for the bulk Si"
  }
}

4.2. 入力パラメータファイル nfinp.dataのタグ（キーワード）の一覧

入力パラメータファイル nfinp.dataのタグ（キーワード）の一覧を、 表 4.2 に示します。

表 4.2 入力パラメータファイル nfinp.dataのタグ（キーワード）の一覧

最上位ブロック

第2, 3ブロック

タグ（キーワード）

説明

control

全体的な計算条件設定ブロック

condition

preparation, -2：入力座標の表示,対称操作の生成,k点の生成までで終了します。

automatic, -1：継続可能であれば計算継続になります。そうでなければ、計算開始になります。

initial, 0 ：計算開始

continuation, 1：計算継続

(以下の2つはekcalによる計算で使用)

fixed_charge, 2：電荷を固定して計算

fixed_charge_continuation, 3：固定電荷+計算継続

デフォルト値はautomatic です。

cpumax

CPU時間の上限（デフォルト：86400 sec）

単位：{sec, min, hour, day}

max_iteration

max_total_scf_iteration

総SCF回数（イタレーション）の制限値

（デフォルト：10000）

max_mdstep

総MD計算数の制限値（デフォルトは無制限）

max_scf_iteration

1MDステップ内のSCF回数の制限値

（デフォルトは無制限）

nfstopcheck

ファイルnfstop.data に書かれた数値で、処理を停止するべき更新回数を決定（デフォルト：1）

sw_ekzaj

phaseで、ekcalの入力となる波動関数ファイルF_ZAJへの出力を行うときONにする。EKCALでそのファイルを読み込むときもONにする。ただし、 \(\Gamma\) 点の計算でしか使用できません。

デフォルト値はOFF です。

accuracy

計算精度の制御用ブロック識別子

cutoff_wf

波動関数のカットオフエネルギー

cutoff_cd

電荷密度のカットオフエネルギー

num_bands

バンド数

ksampling

method

k点のサンプリング法。

monk： Monkhorst-Pack法。

mesh：メッシュを生成します。

file：ファイルから入力

direct_in：直接記述

gamma ：\(\ \Gamma\)点のみ

デフォルトはMonkhorst-Pack法

mesh

メッシュの分割数

nx, ny, nz

x,y,z 方向への分割数

デフォルト値=(4,4,4)、上限値=(20,20,20)

kshift

Monkhorst-Pack法でのみ有効なタグ

k1, k2, k3

メッシュのずれの指定（入力値は[0.0,0.5] の範囲）

デフォルト値：

hexgonal の場合： k1 = k2 = 0, k3 = 0.5

そ れ以外の場合： k1 = k2 = k3 = 0.5

ただし0.5 はメッシュの刻み幅の半分の値を指す

kpoints

k 点の重みづけ

kx ky kz denom weight

k= (kx/denom, ky/denom, kz/denom)

k点の座標値と、その重みづけ

smearing

k点サンプリングのsmearing

method

parabolic ：Parabolic法（デフォルト）

cold ：Cold smearing 法（金属系で有効）

tetrahedral：Tetrahedral法

improved_tetrahedral : 改良tetrahedral法

tetrahedral またはimproved_tetrahedral としたときにはk点のサンプリングをメッシュ法にしなければなりません。

width

smearing幅（デフォルト値：0.001hartree）

method = parabolicとcold の時に使用

(タグなし)

xctype

交換相関エネルギー(LDA,GGA)

LDA ： LDAPW91, PZ

GGA ： GGAPBE, REVPBE

scf_convergence

自己無憧着場の収束判定

delta_total_energy

原子あたりの全エネルギーの計算誤差の上限\(\Delta E\)

（デフォルト値： \(10^{- 9}\) hartree）

succession

ここで指定するsuccession回数連続して誤差が \(\Delta E\) 以内に収まったときにSCF収束したと判断する。構造緩和している、あるいは有限温度計算している場合は原子位置の更新手続きに移る。 （デフォルト値：2）

force_convergence

力の収束判定

max_force

全原子に働く力の最大値がこの値より小さくなれば計算を停止

（デフォルト値：0.001 hartree/bohr）

ek_convergence

固有値の収束判定。ekcalによる計算専用の識別子

num_extra_bands

収束を促進するために追加する追加バンド（これらのバンドの固有値は収束が保証されない）の数。基底関数が増え、収束を促進する効果がある。 （デフォルト値：2）

num_max_iteration

k点一個当たりの最大の更新回数（デフォルト値：300）

sw_eval_eig_diff

固有値評価用スイッチ

{ 1, on, yes }：評価あり（デフォルト）

{ 0, off, no }：評価なし

delta_eigenvalue

固有値の許容誤差

（デフ ォルト値：\(10^{- 5}\) hartree）

succession

計算の繰り返し回数（デフォルト値：2）

(タグなし)

initial_wavefunctions

波動関数の初期値

選択肢：{ random_numbers, matrix_diagon, atomic_orbitals, file}

random_numbers

:乱数で初期化

matrix_diagon:小行列 対角化で初期化

atomic_orbitals:原子軌道で初期化

file:F_ZAJで指定されるファイルにより初期波動関数を与える

matrix_diagon

波動関数の初期値を小行列対角化法で与える

cutoff_wf

波動関数のカットオフエネルギー

(タグなし)

initial_charge_density

電荷密度の初期分布

選択肢：{Gauss, atomic_charge_density, file}

Gauss: ガウス分布関数の重ね合わせで初期化

atomic_charge_density: 原子の電子密度の重ね合わせで初期化

file: ファイルF_CHGTから入力

precalculation

nel_Ylm

予め計算してメモリー上に保持しておく球面調和関数の最高次数（デフォルト値は9）

structure

構造設定用ブロック識別子

unit_cell_type

単位胞の型。選択肢： {primitive, Bravais }

unit_cell

a_vector

b_vector

c_vector

a, b, c

alpha,

beta,

gamma

unit cell 単位胞の指定。以下のいずれかの方法で与える

各格子ベクトルの

\((x,y,z)\) 成分

デフォルトの単位はBohr

格子定数 a, b, c

b–c軸、c–a軸、a–b軸のなす角

(角度のデフォルト単位はdegree)

symmetry

method

選択肢:{manual, automatic}

automaticを選択すると自動的に対称性を決定します。

crystal_structure

選択肢：

{diamond,hexagonal, fcc, bcc, simple cubic}

tspace

柳瀬章「空間群のプログラムTSPACE」(裳華房)および、ABCAPのマニュアルを参照

lattice_system

{rhombohedral, trigonal,r,t,-1}

{hexagonal,h, 0}、{primitive ,simple,p,s,1}

{facecentered ,f,2}、{ bodycentered,b,3}

{bottomcentered ,basecentered ,onefacecenter ed,bot,ba,o,4}

num_generators

生成元の数 (\(1 \sim 3\) の整数値)

generators

生成元

af_generator

磁性空間群の生成元

(第3 タグなし)

sw_inversion

反転対称性の有無

(第2 タグなし)

nspin

(2023.01以降)

スピン自由度

入力値：{1,2,-2}から選択。 2D版では4も選択可能。 4はnoncollinear計算に対応。

magnetic_state

上のnspinが設定されている場合はこの設定は読み飛ばされる。

スピン自由度

入力値：{nonmgnetic,antiferro,magnetic}から選択。 nonmagneticはnonmagあるいはnoneと省略も可。antiferroはafと省略可。 magneticはmagと省略可。また、プログラム開発の歴史的経緯から、para、ferroの設定も可能。paraはnonmagneticとferroはmagneticと同じである。nonmagnetic、nonmag、noneあるいはparaはnspin=1に、mgnetic、magあるいはferroはnspin=2に対応する。また、antiferroあるいはafはnspin=-2に対応する。 2D版ではnoncollinearの設定も可能。 noncollinearはnspin=4に対応する。

noncollinearについては7.13 章 を参照。

atom_list

原子構成

coordinate_system

選択肢 ：{cartesian, internal}

atoms

rx, ry, rz

座標

element

元素名

mobile

可動性

入力値 は{1,0}、{on,off}、{yes,no} のどれでも可

weight

重みづけ

weight = 2 は sw_inversion = onの時のみ有効

このとき、反転対称の位置にも原子を生成

element_list

element

元素名(atomsのelementの入力値と一致させる)

atomic_number

原子番号

mass

質量

zeta

スピン分極 s=(nup-ndown)/(nup+ndown)

deviation

初期電荷をガウス関数の和で与えるときの各ガウス関数の偏差。

タグ名には devやstandard_deviationも使用可

wavefunction_solvers

solver

波動関数ソルバー（詳しくは4.6.2 章 と4.8 章 を参照）

sol

ソルバーの種類

MatrixDiagon ：行列対角化法

lm+MSD： lm(一次元探索)+ MSD(改良 型最急降下法)

RMM2P, RMM3：RMM 法

MSD： 修正最速降下法

pdavidson: 分割Davidson法

pkosugi: 分割Kosugi法

till_n

何回の更新まで、sol で指定された波動関数の

更新方法を適用するかを指定

dts

計算開始時の時間刻み幅

dte

itr で指定された更新の回数における時間刻み

幅。dtsの値のみが入力された場合にはdteにも同じ値を適用

itr

時間刻み幅を変化させる回数の指定

var

補間の形式。選択肢：{linear, tanh}。既定値はlinear

prec

前処理の有無。選択肢：{on,off}

cmix

電荷混合法の指定用変数。 charge_mixing タグのmixing_methods で指定されている、各方法に割り振られた番号を使って指定する。

submat

onのときsubspace_rotationの指定に従ってsubspace rotationを行う。選択肢：{on,off}

line_minimization

一次元探索に関係した制御

dt_lower_critical

dt_upper_critical

一次元探索をおこなう時の時間刻みの下限と上限

(デフォルト値はそれぞれ、0.005と2.0)

delta_lmdenom

rmm

残差最小化法

imGSrmm

RMM法で更新した波動関数に対して、Gram–Schmidtの直交化法を適用する頻度(デフォルト値は、毎回実行のimGSrmm = 1)

rr_Critical_Value

バンド毎の収束判定条件。波動関数の残差のノルムがここで指定された値以下になれば、そのバンドはそれ以降更新されない

edelta_change_to_rmm

波動関数のソルバーをRMM法に変えるときの、全エネルギー収束判定条件。ここで指定する値より全エネルギーの収束が悪いときは、その前のソルバーを続けて使う。デフォルト値は1e-3/natmhartree；ここでnatmは原子数

subspace_rotation

subspace対角化に関する制御

subspace_matrix_size

デフォルトの入力値はバンドの数 (num_bands)

num_bandsよりも大きな値が入力された場合には、強制的にnum_bandsの値を入力値に設定

damping_factor

非対角要素のダンピング係数。[0.0,1.0] の範囲外の値が入力された場合には、入力値を強制的に1.0に設定

period

solverタグのsubmatがONになっている場合、periodに1回subspace_rotationを行います。

例えばperiod=3のときiteration(i)のうち、i=1,4,7,10 ,...がsubspace rotation を行う対象になります。デフォルト値は1。

critical_ratio

非対角項の要素の値(1要素あたり)と対角項の要素の値(1要素あたり)の比がいったんcritical_ratioより小さくなった点に対しては、それ以後subspace rotationを行いません。

デフォルト値は\(10^{-15}\)

charge_mixing

電荷混合法。（詳しくは4.7 章 と 4.8 章 を参照）

mixing_methods

電荷密度の混合法。

method

選択肢：{ simple, broyden2, pulay }

デフォルトはpulay

rmxs

計算開始時の電荷密度を混ぜる割合

デフォルト値は0.4

rmxe

itr 回の更新の後に電荷密度 を混ぜる割合。

デフォルト値は 0.4。rmxs の値のみが入力された場合には、rmxe にも同じ値を適用。

itr

電荷密度の混合比(rmx)を変化させる回数

var

rmxを変化させる方法(itr回のSCF回数の間にrmxsからrmxeまで）。選択肢：｛linear, tanh}

prec

前処理の有無。選択肢：{on, off}

istr

method がsimple 以外の場合に、istr回の更新後に、指定した方法で電荷を混ぜる

nbmix

蓄えておくべき電荷密度データの回数を指定

update

nbmix回分用意されている電荷密度の配列を使い切った時の処理の選択法。

選択肢：{anew,renew}

anewはそれまでのデータを全て棄却して新規に開始。

renew は最も古いデータを最新のデータと入れ換える。

charge_ preconditioning

amix

前処理変数a

bmix

前処理変数b

structure_evolution

構造緩和計算用ブロック識別子

method

選択肢：{sd, quench, gdiis, bfgs, cg, cg2, fire, velocity_verlet, temper ature_control}

dt

時間刻み幅

stress

ストレス計算

sw_stress

ストレス計算の有 無。選択肢：{ on,off }

gdiis

(GDIISおよびBFGS を選択す る場合のタグ)

initial_method

GDIIS (BFGS)へ移行する前に利用する最適化アルゴリズム。選択肢：{ quench, cg, sd }デフォルト値はcg

gdiis_box_size

ここで指定するイオン座標更新回数分のデータをgdiis(bfgs)用配列に蓄える

gdiis_hownew

gdiis_box_sizeで指定した回数分のイオン座標のデータ配列を使い切った時の処理法の選択

選択肢：{anew, renew}

c_forc2gdiis

GDIIS (BFGS)への切替えの判定条件

デフォルト値は0.05 (hartree/bohr)

postprocessing

dos

状態密度の出力

sw_dos

状態密度出力の有無。選択肢：{ on,off }

method

選択肢：{ tetrahedral, Gaussian }

deltaE_dos

状態密度出力のエネルギー精度

variance

mehtodがGaussianの場合のガウス関数の分散

nwd_dos_window_width

出力時のエネルギー幅 \(\Delta\) Eを次式で指定： \(\Delta\) E=nwd_window_width x deldos

charge

電荷の出力

sw_ charge_rspace

電荷出力の有無。選択肢：{ on,off }

filetype

電荷出力ファイルの形式

選択肢：{ cube, density_only }

title

電荷の出力ファイルの見出し

filetype = cube の時のみ有効

printoutlevel

標準出力への出力レベルの制御

0：出力なし

1：情報を出力

2：デバッグ 用の情報を出力

base

他の変数に入力値が指定されていない時は、この値がデフォルト

pulay

Pulay電荷混合法

timing

時間指定情報

solver

電子状態解法

evdff

固有エネルギー差

rmm

残差最小化法

snl

非局所ポテンシャル

gdiis

GDIIS 法

eigenvalue

固有値

spg

空間群

kp

k 点

matdiagon

行列対角法

vlhxcq

ローカルポテンシャル

totalcharge

電子密度

submat

部分空間回転法

strcfctr

構造因子

parallel

並列化のための前処理の 結果の出力制御

input_file

入力ファイルF_INPの解析結果の出力

parallel_debug

1に設定するとゼロ番ノード以外のプロセスからもoutput00x_xxxといったファイルに出力を行う。

jobstatus

計算の進行状況をjobstatus00xに出力

jobstatus_option

状況ファイルの出力制御

jobstatus_format

tag, tag_line, tableが選択可能。既定値はtag。

jobstatus_series

ON またはOFF

4.3. 全体的な計算条件設定（Control）

計算をはじめから実行するのか継続計算を実施するのか、最大どれくらいの時間計算を継続するのか、など、計算全体に関わる条件の設定をcontrol ブロックで行います。たとえば、以下のように記述します。

control{
  condition = initial
  cpumax = 1 day
  max_iteration = 1000000
}

control ブロックにおいては、次の変数を指定することができます。

condition

初期計算か継続計算かなどの指定です。"initial" とすると計算は始めから行われ、"continuation" とすると、波動関数、電荷密度分布などの計算結果を引き継い だ継続計算が行われます。また、"automatic" とすると継続計算に必要な複数のファイル（これらは、前のジョブが正常終了した場合には自動的に生成される）が存在する場合は"continuation"、存在しない場合は"initial"と設定したのと同じ動作をします。"preparation" を指定すると、前処理（使用配列の大きさの評価、ｋ点生成など ）のみ行います。デフォルト値は"automatic" です。 ほかに、（収束した） 電荷密度分布を読み込んで、それを固定したまま波動関数のみを収束させる（バンド分散を計算する場合など）には、"fixed_charge", "fixed_charge_continuation", "fixed_charge_automatic"のいずれかを設定します。それぞれ、最初から計算するか、継続計算するか、自動判定するかの指定です。 これらの"initial"、"continuation"、 "automatic"、"preparation"、fixed_charge" 、"fixed_charge_continuation"、"fixed_c harge_automatic"は、それぞれ、整数0、1、- 1、-2、2、3、-3で代用することができます。

cpumax

PHASE計算を実行する時間を、実数＋単位、の組合せで指定します。ここで指定した時間を超えると、収束に達していなくても、継続計算用のファイルなどが出力され、計算が停止します。デフォルト値は86400 s (1 日) です。単位は必須です。使える単位は、"sec"、"s"（これは"sec"と同じ）、"min"、"hour"、および"day"です。継続計算になる可能性がある場合には、ジョブの指定時間よりも小さな値に設定しておくのがよいでしょう（例えば、ジョブの制限値が６時間で、収束後の状態密度計算などの後処理がなければ、5.8hour程度に設定する）。

max_iteration

max_total_scf_iteration

SCF計算の総イタレーション数の最大値を指定します。SCF計算の総イタレーション数がここで指定した数に達すると継続計算用のファイルなどが出力され、計算が停止します。デフォルト値は10000です。

max_scf_iteration

構造緩和計算や分子動力学計算における各MDステップ内での電子状態の更新回数（SCFイタレーション数）の最大値を指定します。例えば、構造緩和計算の最初期に構造が不安定で、電子状態に対する収束判定条件を満たすまで数百回に及ぶような大きな回数のSCFイタレーションが必要になることがあります。その場合には、SCF計算を途中で打ち切って、力を計算してより安定な原子構造に更新してから、次のSCF計算をすすめた方が有利になります。しかし、あまり小さな値（１０程度）に設定すると、計算される力の誤差が大きくなり、逆に収束を難しくすることがあるので注意が必要です。正確な力の計算が重要な場合には使用しないでください。

4.3.1. ブロックサイズ

control ブロックにおいて演算に用いるブロックサイズを指定することができます。特に大規模な系の場合デフォルトのブロックサイズは適切でない可能性があるため、適切なパラメーターを探ることによって高速な計算が実現できます。 設定可能な主なブロックサイズは下記の通り。

表 4.3 主なブロックサイズとデフォルト値

nblocksize_mgs

修正グラムシュミット

8

nblocksize_betar_dot_wfs

波動関数と射影演算子の積

32

nblocksize_vnonlocal_w

波動関数と非局所ポテンシャルの積

1000

nblocksize_submat

部分空間対角化

0 (ブロッキングをしない)

以下の記述を行うとブロックサイズの探索を行ってくれます。（3次元並列版のみ）

control{
  sw_optimize_blocking_parameters = on
}

この設定を施すと、前処理の際各ルーチンをブロックサイズを変えながら小数回実施し、最もよかったブロックサイズを報告します。Printlevelブロックにおいてipriblsizeを2以上に設定している場合、ログファイルには以下のような情報が記録されます。

!** MGS blocksize and elapsed time   525     10.0516
!** MGS blocksize and elapsed time   262      5.3681
!** MGS blocksize and elapsed time   131      3.2394
!** MGS blocksize and elapsed time    65      2.7701
!** MGS blocksize and elapsed time    32      3.5466
!** MGS blocksize and elapsed time    16      3.7317
!** MGS blocksize and elapsed time     8      5.3605
!** MGS blocksize and elapsed time     4     10.3687

最終的に得られた最適なブロッキングパラメーターは以下のように出力されます。

nblocksize_mgs           =    65
nblocksize_betar_dot_wfs =  1024
nblocksize_vnonlocal_w   =  4096
nblocksize_submat        =   700

決まったパラメーターは計算でそのまま活用されます。この文字列は入力パラメーターファイルのcontrolブロックの下にそのまま貼り付けて利用することができます。

なお、この機能は三次元版並列版においてのみ有効です。また、とk点並列とScaLAPACKをすべて同時に利用することはできないのでご注意ください。

4.3.2. 高速計算のオプション (2023.01以降)

ControlブロックにおいてSCF計算を高速にするための以下のようなオプションを利用することができます。

表 4.4 SCF計算高速化オプション

オプション

説明

sw_keep_hloc_phi

onにすると「ローカルポテンシャルに波動関数を作用した配列」をメモリーに保存し再利用することによって高速化を実現します。デフォルト値はonです。offにすると使用メモリーを削減することができます。

sw_precalculate_phase_betar_dot_wfs

波動関数と射影演算子の積の演算において位相をあらかじめ計算しておくかどうかを指定するスイッチ。有効にすると計算が高速化されるがメモリー消費が多くなります。デフォルト値はoff.

sw_precalculate_phase_vnonlocal

波動関数と非局所ポテンシャルの積の演算において位相をあらかじめ計算しておくかどうかを指定するスイッチ。有効にすると計算が高速化されるがメモリー消費が多くなります。デフォルト値はoff.

sw_reduce_fft_for_charge

RMM3ソルバー利用の際に高速フーリエ変換の回数を減らすオプション。有効にするとFFTの回数が減るが、収束性に若干の影響を及ぼす可能性がある。デフォルト値はoff.

sw_keep_hloc_phi は特に効果が高いため、デフォルト値がonになっています。メモリーが枯渇してしまうような問題に対してはoffに設定するようにしてください。

4.4. 計算精度の指定（Accuracy）

4.4.1. カットオフエネルギー

カットオフエネルギーは平面波基底を利用した計算においては計算の信頼性を決める重要なパラメーターです。

カットオフエネルギーは以下のように指定します。

accuracy{
   cutoff_wf = 25 Rydberg
   cutoff_cd = 225 Rydberg
}

cutoff_wf

波動関数のカットオフエネルギーをエネルギーの単位で指定します。

cutoff_cd

電荷密度のカットオフエネルギーをエネルギーの単位で指定します。

カットオフエネルギーは充分な精度が得られる値を事前に勘案することが理想的ですが、以下のような指針も有用です。

  • cutoff_wfはおおよそ25 rydberg

  • cutoff_cdは、ノルム保存型の擬ポテンシャルのみを利用している場合はcutoff_wfの4倍、ウルトラソフト型擬ポテンシャルを利用している場合は9倍。バージョン2020.01以降はこのような値がデフォルト値になりました。

4.4.2. バンド数

バンド数は、以下のようにaccuracyブロックの下のnum_bands変数によって指定します。

accuracy{
   num_bands = 12
}

num_bands

バンド数

バンド数は、最低限価電子数の半分+1は必要です。通常最低必要な数の2 割程度多めの数を採用します。設定値が価電子数の半分以下の場合には、自動的に設定が増やされます。またこの値を設定していない場合には自動的にバンド数が設定されます。

4.4.3. k点サンプリングとスメアリング

4.4.3.1. 基本の設定

カットオフエネルギーと同様に、k点サンプリングも計算の信頼性を決める重要なパラメーターです。 k点サンプリングは, accuracy ブロックの下にksampling ブロックを作成し、ksampling ブロックの下で設定を行います。たとえば下記のようになります。

accuracy{
   ksampling{
     method = monk
     mesh{
       nx=4
       ny=4
       nz=4
     }
   }
}

ksampling ブロックでは、下記の変数/ブロックを定義することができます。

method

k点サンプリングの方法を選びます。monk, mesh, file, gamma、directinのいずれかです。monkはMonkhorst-Pack法によるサンプリングで、通常推奨される方法であり、デフォルト値です。meshは単純なメッシュで逆空間を分割します。 四面体法により電荷密度を構成する場合や状態密度の計算を行う場合にはこれを指定します。fileはファイルから読み込みます。 バンド分散をみるために対称線に沿って多くのｋ点を入力する必要がある場合やフェルミ面の計算において大量のk点を考慮する必要がある場合などに利用します。gammaを指定すると\(\Gamma\)点のみをサンプリングします。充分大きな単位胞を使っていて、\(\Gamma\)点のみでも充分な精度が得られる場合には、これを指定します。directinは直接ｋ点の組（個数と座標）を指定します。いずれの方法でも、サンプリングｋ点に\(\Gamma\)点が含まれていて、系に反転対称中心がなければ（設定されていなければ）、\(\Gamma\)点の波動関数に関する計算は、この点の対称性を利用して他のｋ点のものに比べて３倍程度高速に実行されます（後述のとおり、これを抑制する、つまり他のｋ点と同じ演算法を適用する手段もあります）。

mesh

逆空間の分割数を指定します。以下の変数が利用できます。

nx 1番目の逆格子ベクトルの分割数を指定します。

ny 2番目の逆格子ベクトルの分割数を指定します。

nz 3番目の逆格子ベクトルの分割数を指定します。

スメアリングは、フェルミ準位付近の状態を“ぼやかす” 操作です。 これによって、フェルミ準位付近で状態を持つ金属系においても少ないk 点数で高い精度で計算ができるようになる場合があります。 スメアリングは、以下のようにaccuracy ブロックの下のsmearing ブロックにおいて指定します。

accuracy{
   smearing{
     method = parabolic
     width = 0.001 hartree
   }
}

smearing ブロックでは以下の変数を利用することができます。

method

スメアリングの方法を指定します。parabolic、tetrahedral、 cold、improved_tetrahedral、methfessel_paxtonのいずれか を指定します。通常利用するのはparabolicで、２次関数の組み合わせによって フェルミ準位付近をぼやかします（後述）。tetrahedralとimproved_tetrahedralは四面体法で、主に四面体法による状態密度計算を行う場合に利用します。coldはColdスメアリングで、金属において有効とされている方法です。

width

スメアリングの幅をエネルギーの単位で指定します。デフォルト値は0.001 hartreeです(method=parabolicの場合)もしくは0.01 hartree (method=methfessel-paxtonの場合)

この設定はmethod=parabolicもしくはmethfessel_paxtonの場合に有効です。 気を付けなければならないのは、method=tetrahedralの場合、このwidthは別の意味を持つことです。この場合、widthは"縮退しているとていると見なされる準位間のエネルギー差の閾値" になります（指定がない場合の既定値は1.e-5 hartree）。 method=parabolicのときに設定したwidthをmethod=tetrahedralの場合にもそのまま残しておかないように注意しましょう。

Widthで指定された値をwとします。method=parabolicの場合、各バンドの固有エネルギーεを下の図で示す状態密度分布を持つものとして扱います。ε-wからε+wまでの間は上に凸な二次関数、ε-2wからε-wまでの間とε+wからε+2wまでの間はそれぞれ、下に凸な二次関数で表され、それらが連続に接続しています。この状態はε-2wからε+2wまで積分して１になるように規格化されています。これにより各固有状態の占有割合が0から1の間の値を取るようになります。フェルミエネルギー±2wの間にある状態がこのスメアリングの影響を受けます。フェルミエネルギー近傍に状態が複数あり、SCF計算（および構造緩和計算）中に値の順序が入れ替わるような場合には、これらの固有エネルギーがフェルミエネルギー±2wの範囲に入るように設定すると（基底状態の電荷密度分布の変化が抑制され）収束性が改善されることがあります。しかし、wが大きい程、DFTの基底電子状態からのずれも大きくなるので注意が必要です。

../_images/parabolic.svg

図 4.1 Smearing widthと状態密度の関係。横軸はエネルギー、縦軸は状態密度。

4.4.3.2. k点サンプリングを“密度”で指定する方法（バージョン2020.01以降）

バージョン2020.01より、k点サンプリングをメッシュ数ではなく“密度”で指定することができるようになりました。以下のように設定します。

accuracy{
  ksampling{
    density = 4 bohr
  }
}

densityは長さの単位で指定します。meshブロックにおけるメッシュの指定があればそちらが優先されます。デフォルト値は4 bohrであり、これはSi結晶の場合に4×4×4のメッシュを指定することに相当します。このデフォルト値で問題ないのであれば、ksamplingブロックを丸ごと省くことも可能です。

4.4.3.3. k点サンプリング指定のデフォルトの振る舞いの変更（バージョン2020.01以降）

これまではksamplingのmethodのデフォルト値は常にmonkでしたが、バージョン2020.01以降smearingブロックのmethodがtetrahedralの場合に限りデフォルト値がmeshとなるように変更されました。これは、スメアリング手法としてtetrahedral法を利用する場合k点サンプリング手法がmeshであることが必要なためです。この仕様変更のため、「smearingのmethodがtetrahedralでksamplingのmethodとしてmeshを明示的に指定していない」場合以前のバージョンと等価な計算にならない点には留意が必要です。

4.4.3.4. Methfessel-Paxton法によるスメアリング (バージョン2021.01以降)

概要

PHASE/0にはparabolic法、四面体法などのスメアリング手法が備わっています。バージョン2021以降、これにMethfessel-Paxton法 [Methfessel89] が加わりました。Methfessel-Paxton法は、負の占有数を許容するかわりに大きなスメアリング幅を利用しても全エネルギーなどの重要な計算結果がほとんど変化しない、という特徴を持ったスメアリング手法です。

Methfessel-Paxton法について

基本の理論

スメアリング手法とは、 \(\delta\) 関数および階段関数を数値的に安定でなめらかな関数で近似する手法ということができます。Methfessel-Paxton法においては、 \(\delta\) 関数および階段関数を以下のように近似します。

\[D \left( x \right) = \sum_{n=0}^N [ A_n H_{2n} {\rm e}^{-x^2} ]\]
(4.1)\[S \left( x \right) = \frac{1}{2} \left(1 - {\rm erf}(x)\right)\]

\(A_n = \frac{(-1)^n}{n!4^2\sqrt(\pi)}\), \(H_{n}\) は次数 \(n\) のエルミート多項式です。 実際の計算においては、準位 \(ik\) における固有値を \(\varepsilon_{ik}\) フェルミエネルギー \(\varepsilon_F\) スメアリング幅を \(\sigma\) とすると \(x\) は \(x_{ik}=\frac{\left(\varepsilon_{ik}-\varepsilon_F\right)}{\sigma}\) という形で用いられます。また、エントロピー項は以下のように評価することができます [Kresse96]

(4.2)\[S_N = \sigma \sum_{ik} \frac{1}{2} A_N H_{2N} \left(x_{ik}\right) {\rm e}^{-x_{ik}^2}\]

図 4.2 に、 N =6までの \(D(x), S(x)\) を示します。

../_images/mp1.svg

図 4.2 N = 6までの \(D(x)\) および \(S(x)\)

図 4.2 から明らかなように、Methfessel-Paxton法の場合フェルミエネルギー付近（ \(x=0\) 付近）の電子の占有数が負になる場合があります。これは明らかに物理的ではありませんが、最終的な全エネルギーの計算値などに影響を及ぼすことはありません。

フェルミエネルギーについて

通常フェルミエネルギーは電荷の数が合うように二分法によって求めますが、Methfessel-Paxton法においては負の占有状態が発生するため通常の二分法ではフェルミエネルギーが一意に求まらない可能性があります。そこで、最低エネルギーから少しずつエネルギーを増やし、電子数が超えるエネルギーを見つけてから通常の二分法に移行する、という特殊な手続きでフェルミエネルギーを求めます。具体的な手続きは下記の通り。

  1. スメアリングに利用するスメアリング幅とMethfessel-Paxton法の次数からエネルギー幅 \(\delta\) を設定します。\(\delta=a \sigma/N\). \(a\) は入力パラメーターファイルで設定可能な調整値で、PHASE/0のデフォルト値は10です。

  2. 探索を開始する最低エネルギー値を設定します。求まっている固有値の最低値でもよいですが、保守的すぎるので \(e_{\rm min} + \left(e_{\rm max}-e_{\rm min}\right)/2n\) という値を採用します。

  3. \(\delta\) を最低エネルギーに加算していきフェルミエネルギーとし、結果得られる仮の電荷が本来の電荷を上回った段階でそのエネルギーを上限値とします。

  4. 以上の処理によって下限と上限がエネルギー幅の範囲で決定するので、ここからは通常の二分法に移行し、フェルミエネルギーを求めます。

使い方

スメアリングの設定は、入力パラメーターファイルのaccuracyブロックの下のsmearingブロックにおいて行います。Methfessel-Paxton法の設定例を以下に記述します。

accuracy{
  smearing{
    method = methfessel_paxton
    width = 0.01 hartree
    methfessel_paxton{
      order = 2
    }
  }
}

methodにmethfessel_paxtonを指定すると(methでも可)スメアリング手法がMethfessel-Paxton法になります。widthにおいてメアリング幅をエネルギーの単位で指定します。Methfessel-Paxton法の詳細設定はmethfessel_paxtonブロック（もしくはmethブロック）において行います。ここでは以下のような設定を施すことができます。

表 4.5 Methfessel-Paxton法の詳細設定

タグ名

説明

order

Methessel-Paxton法の次数を指定します。デフォルト値は2です。

例題

samples/smearing 以下に例題が配置されています。 FCC-AlとBCC-Feの例題が格納されています。それぞれの結晶について、mp2 mp4 mp6 parabolic の4つのディレクトリーがあり、 それぞれ N = 2, 4および6の場合のMethfessel-Paxton法を適用した入力例と従来のparabolic法を適用した例題に対応します。 さらにその下のサブディレクトリーには、sigma0.0001(スメアリング幅を0.0001 hartreeにした場合)とsigma0.001(スメアリング幅を0.001 hartreeにした場合) の例題が配置されています。これらの例題から得られる全エネルギーの計算結果を次の表にまとめます(表中のエネルギーはすべてhartree単位)

表 4.6 Methfessel-Paxton法および従来法によって得られる全エネルギー。

\(\sigma=0.0001\) hartree

\(\sigma=0.001\) hartree

Al, mp2

-0.2552794

-0.2552805

Al, mp4

-0.2552794

-0.2552801

Al, mp6

-0.2552794

-0.2552801

Al, parabolic

-0.2552797

-0.2552832

Fe, mp2

-21.993246

-21.993248

Fe, mp4

-21.993246

-21.993248

Fe, mp6

-21.993246

-21.993247

Fe, parabolic

-21.993247

-21.993259

表の数値から、Methfessel-Paxton法はparabolic法と比較して全エネルギーがスメアリング幅に依存しにくいことが分かります。

[Methfessel89]

Methfessel and A. T. Paxton, Phys. Rev. B 40 (1989) 3616.

[Kresse96]

G. Kresse and J. Furthmüller, Computational Materials Science 6 (1996) 15.

4.4.4. 交換相関エネルギー

交換相関エネルギーは、LDAとGGAの２種類があります。LDAはLDAPW91, PZ、GGAはGGAPBE, RPBE, REVPBEが利用できます。

accuracy{
  xctype = ggapbe
}

xctype

交換相関エネルギー(LDA, GGA)

LDA ： LDAPW91, PZ

GGA ： GGAPBE, RPBE, REVPBE

4.4.5. 収束判定

収束判定は、電子状態計算の収束判定と構造最適化の際の原子に働く力の収束判定の2 種類があります。以下のように指定します。

accuracy{
   scf_convergence{
     delta_total_energy = 1.0E-8 Hartree
     succession = 3
   }
   force_convergence{
     max_force = 2.0E-4 Hartree/Bohr
   }
}

収束判定に関わるブロック/変数は以下の通りです。

scf_convergence

SCF計算の収束判定を指定するブロックです。

delta_total_energy

全エネルギーの差の閾値をエネルギーの単位で指定します。

現在の全エネルギーと1ステップ前のエネルギーの差がここで指定した値よりも小さい場合収束判定を満たしたとみなされます。デフォルト値は1e-9です。

succession

delta_total_energyを何回連続でみたせば最終的に収束したと見なすかを指定します。ここで指定した値の回数連続で収束判定を満たせば収束が得られたと判定されます。デフォルト値は2です。

force_convergence

原子に働く力の最大値に関する閾値を設定するブロックです。

max_force

原子に働く力の最大値の閾値を力の単位で指定します。デフォルト値は1e-3です。

4.4.6. 初期波動関数と初期電荷密度

初期波動関数と初期電荷密度の設定を適切に行うと、電子状態計算を少ない回数で収束させることができます。初期波動関数および初期電荷密度は、以下のように設定することができます。

accuracy{
   initial_wavefunctions = atomic_orbitals
   intial_charge_density = atomic_charge_density
   matrix_diagon{
     cutoff_wf = 5 rydberg
   }
}

初期波動関数および初期電荷密度の設定に関わるブロック/変数は下記の通りです。

initial_wavefunctions

初期波動関数の設定方法を指定します。

random_numbers, matrix_diagon, file, atomic_orbitalsを選択することができます。random_numbers は乱数による初期化です。matrix_diagon は行列対角化によってもとめます。この際、初期波動関数作成時にのみ利用するカットオフエネルギーを採用することもできます。その指定方法は後述のmatrix_diagonブロックにおいて行います。fileは、波動関数ファイルから読み込みます。すでにある程度収束した波動関数データファイルを持っている場合はこのオプションを指定し、読み込ませることができます。atomic_orbitalsは、擬ポテンシャルファイルに記録された原子軌道データをもとに初期化を行います。デフォルト値はrandom_numbersです。

initial_charge_density

初期電荷密度の設定方法を指定します。Gauss, file, atomic_charge_density のいずれかを選択することができます。Gaussは原子を中心とした単純なガウス関数による初期化です。fileはファイルから読み込みます。すでにある程度収束した電荷密度データファイルを持っている場合はこのオプションを選択し、読み込ませることができます。atomic_charge_densityは擬ポテンシャルファイルに記録された原子の電荷密度をもとに初期化を行います。デフォルト値はGaussです。

matrix_diagon

initial_wavefunctionsにmatrix_diagonを指定している場合に、その振る舞いを制御するためのブロックです。

cutoff_wf

初期波動関数作成時に利用するカットオフエネルギーの値を指定します。デフォルト値は、通常のカットオフエネルギーの半分です。

4.4.7. カットオフ/格子定数の異なる計算から出力された波動関数/電荷密度を読み込む方法（バージョン2020.01以降）

initial_wavefunctions, initial_charge_densityでfileを指定するとファイルから波動関数や電荷密度を読み込むことができますが、2020.01未満のバージョンではカットオフエネルギーおよび格子定数が等しい場合のみ読み込むことができました。 バージョン2020.01以降は、カットオフが異なる/格子定数が異なる（すなわちメッシュが異なる）場合は新しいメッシュ上に値を補間して読み込むことができるようになりました。カットオフエネルギーもしくは格子定数が異なる波動関数データ/電荷密度データを読み込む場合以下のような設定を施します。

accuracy{
  sw_read_pwbs_info = on
  initial_wavefunctions = file
  initial_charge_density = file
}

このようにすると、nfpwbs.dataファイル(file_names.dataにおいて識別子F_PWBSを利用して変更可能)から必要な情報が読み込まれ、波動関数や電荷密度が補間をしながら読み込まれるように動作します。nfpwbs.dataファイルは計算終了時に自動的に書き出されますが、この動作を抑制したい場合は次の指定を行います。

accuracy{
  sw_write_pwbs_info = off
}

sw_read_pwbs_infoのデフォルト値はoff, sw_write_pwbs_infoのデフォルト値はonです。

4.4.8. 実空間法

PHASEは、非局所ポテンシャルの演算を逆空間で実行しますが、これを実空間で行わせることも可能です。この機能を利用するためには、以下のように設定します。

accuracy{
  nonlocal_potential{
    sw_rspace = on
    r0_factor = 1.9
  }
}

実空間法は、 [King-Smith91] および [Wang01] の方法で実現されています。 逆空間法はその演算量が \(O(N^3)\) であるのに対し実空間法は \(O(N^2)\) なので、大きな系においては実空間法の方が有利となります。 ただし、逆空間法では厳密解が得られるのに対し、実空間法は近似解しか得られない点には注意が必要です。 nonlocal_potentialブロックでは以下のような設定を施すことが可能です。

sw_rspace

実空間法を利用するかどうかを指定します。

デフォルト値はoffです。

projector_optimization

実空間法を適用するためにはプロジェクターの最適化を行う必要がありますが、その方法を指定します。このパラメーターにprefittingを指定すると[1]の方法で、mask_functionを指定すると文献[2]の方法でこの最適化が行われます。デフォルト値はmask_functionです。

r0_factor

「最適化されたプロジェクター」のおよぶ範囲を、もとのプロジェクターの何倍にするかを指定する実数。デフォルト値は1.9。
[King-Smith91]

R. D. King-Smith, M. C. Payne, and J. S. Lin, “Real-space implementation of nonlocal pseudopotentials for first-principles total-energy calculations”, Physical Review B 44 13063 (1991).

[Wang01]

Lin-Wang Wang, “Mask-function real-space implementations of nonlocal pseudopotentials”, Physical Review B 64 201107 (2001).

4.5. 原子構造（Structure）

計算に利用するモデルの指定は、structure ブロックの下で行います。たとえば、以下のようになります。

structure{
    unit_cell_type = Bravais
    unit_cell{
        #units angstrom
        a_vector = 4.914100000 0.000000000 0.000000000
        b_vector = -2.457050000 4.255735437 0.000000000
        c_vector = 0.000000000 0.000000000 5.406000000
    }
    atom_list{
        coordinate_system = Internal
        atoms{
            #units angstrom
            #tag element rx ry rz
             O 0.413100000054 0.145400000108 0.118930000000
             O 0.854599999943 0.267699999886 0.452263333333
             O 0.732300000003 0.586900000006 0.785596666667
             O 0.267699999946 0.854599999892 0.547736666667
             O 0.145399999997 0.413099999994 0.881070000000
             O 0.586899999939 0.732299999879 0.214403333333
             Si 0.530000000000 0.000000000000 0.333333000000
             Si -0.000000000072 0.529999999857 0.666666333333
             Si 0.469999999954 0.469999999908 0.999999666667
        }
    }
    element_list{
        #tag element atomicnumber mass zeta deviation
         O 8 29164.9435 * *
         Si 14 51196.4212 * *
    }
    symmetry{
          method = automatic
          sw_inversion = off
    }
}

4.5.1. ユニットセル

unit_cell_type

単位胞の指定方法を設定しています。

prmitiveかbravais を指定することができます。デフォルト値はbravais です。後述するように、単位胞を格子定数で指定する場合はこの変数をbravaisとする必要があります。また、bravaisを指定している場合、symmetryブロックの下のtspaceブロックにおいて定義できるlattice_system変数によって格子を変換させることが可能です。lattice_system変数についてはあとの説明も参照してください。

unit_cell

単位胞を指定するブロックです。セルベクトルを指定する方法と格子定数を指定する方法があります。格子定数によって指定する方法は、unit_cell_typeがbravaisの場合のみ有効です。

  • セルベクトルを指定する方法

この方法を利用する場合、ベクトル型データを利用して以下のように記述します

unit_cell{
  #units angstrom
  a_vector = a1 a2 a3
  b_vector = b1 b2 b3
  c_vector = c1 c2 c3
}

a_vector, b_vector, c_vectorによってそれぞれ\(a\)軸, \(b\)軸, \(c\)軸をベクトルで指定します。この指定方法の場合、長さの単位はブロック単位で指定する方法のみ利用できる点に 注意してください。この例では、unit_cellブロックの先頭に#units angstromとすることによって長さの単位をÅ単位に変更しています。

  • 格子定数によって指定する方法

この方法を利用する場合、以下のように記述します。

unit_cell{
  a = a0
  b = b0
  c = c0
  alpha = alpha0
  beta  = beta0
  gamma = gamma0
}

a, b, c, alpha, beta, gammaという変数を利用することによってそれぞれ格子定数\(a,b,c,\alpha,\beta,\gamma\)を指定します。 この方法で指定すると、セルベクトルは計算開始時に以下のような“下三角”形式で定義されるようになります。

a_vector = a1 0.0 0.0
b_vector = b1 b2  0.0
c_vector = c1 c2  c3

  • 簡易入力 (バージョン2022.01以降)

対称性によっては単位胞の記述を簡略化することができます。基本的なルールは

  • 格子定数b,cのデフォルト値はaに指定した値

  • 格子定数 \(\alpha \beta \gamma\) のデフォルト値は90度。ただし lattice_system の値が hexagonal の場合は \(\gamma\) のデフォルト値は120度

というものです。いくつか具体例を挙げます。

例 1.

structure{
  unit_cell{
    #units angstrom
    a = 3.5
  }
}

立方晶のため格子定数aの値のみ指定しています。

例 2.

structure{
  symmetry{
    tspace{
      lattice_system = hexagonal
    }
  }
  unit_cell_scale = 3.5 angstrom
  unit_cell{
    a = 1
    c = 1.61
  }
}

lattice_system = hexagonal によって六方晶であることを指定しています。さらに unit_cell_scale = 3.5 angstrom によって単位胞のすべての要素を3.5倍することを指定しています。3.5Åが格子定数aと等しいとすると、a = 1となり、またcの値としてはc/a比を指定すればよいことになります。

例 3.

structure{
  unit_cell_scale = 3.5 angstrom
  unit_cell{
    a_vector = 1 0 0
    b_vector = 0 1 0
    c_vector = 0 0 1
 }
}

unit_cell_scale= 3.5 angstrom として単位胞全体のスケールを指定しています。立方晶で格子定数と unit_cell_scale が等しいのであればセルベクトルを単位行列とすれば正しい指定となります。

4.5.2. 原子座標

atom_list

各原子の座標値の指定などを行うブロックです。

以下の変数/ブロックを定義することができます。

coordinate_system

原子座標を、カルテシアン座標によって定義するかフラクショナル座標によって定義するかを指定します。internalとするとフラクショナル座標によって、cartesian とするとカルテシアン座標によって指定します。デフォルト値はinternalです。

atoms

実際に各原子の座標値を指定する表形式データを記述します。代表的な属性値は以下の通りです。

element

元素名を指定します。元素名は、後述の'element_list'に おいて定義されている必要があります。

rx

\(x\)座標を指定します。

ry

\(y\)座標を指定します。

rz

\(z\)座標を指定します。

mobile

構造最適化や分子動力学シミュレーションにおいてこの原子が"可動か否"かを指定する真偽値です。可動にしたい場合onとします。デフォルト値はoffです。

mobilex

構造最適化や分子動力学シミュレーションにおいてこの原子のx, y, z座標を個別に"可動か否か"を指定する真偽値です。可動にしたい場合onとします。デフォルト値は、mobileに指定された値です。

mobiley

mobilez

weight

"重み"を設定します。この属性値に2という値を与えた場合、原点を中心とした反転対称の位置にコピー原子を配置します。デフォルト値は1です。2を与える場合、系に反転対称性があり、後述のsw_inversionパラメーターがonとなっている必要があります。

4.5.3. 原子種の指定

element_list

元素情報を指定するための表形式データを記述するブロックです。

代表的な属性値は以下の通り。

element

元素名を指定します。指定は必須です。

atomicnumber

原子番号を指定します。指定は必須です。

mass

質量を指定します。

zeta

スピンを考慮している場合の、初期スピン分極の値を設定します。

qex

電子を追加/削減する場合に指定します。

擬ポテンシャルファイルは、file_names.data ファイルにおいて、ファイルポインターF_POT(n) によって指定します。ここでn は入力における元素指定の順序に対応する整数です。たとえば、以下の要領でO とSi の元素指定が入力ファイルにおいて成されていて、

structure{
   ...
   ...
   element_list{
     #tag element atomicnumber mass zeta deviation
         O 8 29164.9435 * *
         Si 14 51196.4212 * *
   }
}

対応する擬ポテンシャルファイルがOがO_ggapbe_us_01.pp, Si がSi_ggapbe_nc_01.pp,だった場合、file_names.dataを以下のように記述します。

&fnames
F_INP = './nfinp.data'
F_POT(1)=’./O_ggapbe_us_01.pp'
F_POT(2)=’./Si_ggapbe_nc_01.pp
/

擬ポテンシャルファイルの指定は、交換相関ポテンシャルの計算方法の指定にも対応しています。公開している擬ポテンシャルファイルの交換相関ポテンシャルの計算方法はggapbeかldapw91 のいずれかですが、どちらなのかはファイル名から判定することができます。すなわち、ggapbe あるいはldapw91 という文字列が擬ポテンシャルファイルのファイル名に含まれています。なお、ggapbe とldapw91 を混在させた計算を行うことはできませんのでご注意ください。また、利用できる原子種は16 種類までです。

電子を追加/削除したい場合の設定もここで行います。qexパラメータに追加したい/削除したい電子数を指定します。たとえば、以下のように指定します。

structure{
  ...
  ...
  element_list{
    #tag element atomicnumber mass zeta deviation qex
          O 8 29164.9435  * * +1
         Si 14 51196.4212 * * 0
  }
}

追加しているのは電子なので、1つ追加すると電荷としては1減る点に注意が必要です。指定は原子種に対して行いますが、実際は系全体におよぶ設定です。SCF計算がすすむに従って各原子に分配される余剰電荷は変化します。具体的には波動関数から電荷密度をつくるときに余剰電子1個分を含むようにフェルミレベルを調整します。また、Ewaldエネルギーはこの余剰電荷を打ち消すように背景に一様電荷分布をおいた状態にして計算します。

4.5.4. 原子種の指定（バージョン2020.01以降）

バージョン2020.01以降、structureブロックにおけるelement_listの指定は非必須となりました。原子配置指定における元素名が通常の元素名もしくは元素名+数値という文字列の場合(C, C1, O2など)、element_listが存在しなくともデフォルトの設定が採用されます。qexやzetaなどの指定が必要な場合はelement_listを作成してください。

4.5.5. 系全体への電荷の付加

系に付加する電荷については、qex 以外にも、additional_charge による指定も可能です。 additional_charge で指定する数字が、正の場合は正電荷、負の場合は負電荷 (電子) を付け加えます。 なお、Ewaldエネルギーはこの余剰電荷を打ち消すように背景に一様電荷分布をおいた状態にして計算します。

下記の例では、系に -1.0 の電荷、すなわち1 個の電子を付け加えます。

例：additional_charge の使用例

structure{
  charged_state{
    additional_charge = -1.0d0
  }
}

4.5.6. 磁気モーメントの初期値の指定

原子種の初期磁気モーメント (分極) については、 4.5.3 章 のzeta 以外にも、moment による指定も可能です。 単位はボーア磁子 [ \(\mu_B\) ] です。 なお、momentによる指定を行うと、プログラム内部で自動的に zeta 値に変換されます。

下記の例では、Cr1 及び Cr2 の初期磁気モーメントに、それぞれ \(3\mu_B\) および \(-3\mu_B\) が設定されます。

例：moment による初期磁気モーメントの指定

structure{
  element_list{
    #tag element atomicnumber moment
        Cr1 24 3.0
        Cr2 24 -3.0
        O 8
  }
}

4.5.7. 対称性の指定

symmetry

系の対称性を指定するブロックです。

対称性を利用することによって、計算量を大幅に減らすことができる場合があります。以下のブロック/変数を利用することができます。

method

対称性指定の方法を指定します。manualとautomaticを選ぶことができます。manualを選択すると、生成元を直接入力することによって対称性を指定することができます。automaticを選択すると、PHASEが指定のモデルから自動的に対称性を検知し、計算に反映させます。デフォルト値はmanualです。

sw_inversion

系に反転対称性が存在する場合に、それを活用して計算量を減らすかどうかを指定する真偽値です。onの場合反転対称性を利用します。反転対称性の中心は原点(0,0,0)です。このオプションは反転対称性のある系の計算を行う場合は有効にすることが推奨されますが、反転対称性のない系で有効にすると前処理で検知し、終了するのでご注意ください。

tspace

TSPACEを利用して生成元を直接指定するためのブロックです。以下の設定を行います。

lattice_system

unit_cell_typeがbravaisの場合に、"格子の型"を指定します。 選択肢はfacecentered, bodycentered, basecentered, rhombohedralです。この変数を指定すると、指定に応じて格子が変換されます。このそれぞれがどのように格子を変換するかについては 表 5.2 を参照してください。この変数を利用することによって、入力ファイルでは指定のしやすいブラベー格子で単位胞を指定しつつ、実際の計算は負荷の少ない基本格子で実行することが可能となります。lattice_system変数を利用して格子を変換させる場合、変換されるのは単位胞のみです。したがって、原子配置の定義などは基本格子の場合の定義の仕方を採用してください。たとえば、面心立方格子の計算を行う場合面心位置の原子は指定せず、原点の原子のみ指定するようにしてください。また、\(k\)点サンプリングは変換後の格子に合わせて設定してください。

generators

生成元を指定するテーブルです。生成元は、3つまでしか定義できないという制約があります。このテーブルでどのように生成元を指定するかについては 5.2.1.2 章 を参照してください。

4.5.8. 原子配置を別ファイルで指定する方法（バージョン2019.02以降）

原子配置を、F_INPファイルではなく別のファイルによって指定することもできるようになっています。

座標データファイルを別ファイルにするためには、structureブロックおよびstructureブロックの下に定義するfileブロックにおいて関連する設定を施す必要があります。関連する変数/ブロックは下記の通りです。

変数名/ブロック名

変数名

説明

method

座標データファイルの指定の仕方を指定する文字列。directinとするとF_INPファイルから、fileとすると別ファイルから読み込む。デフォルト値はdirectin。

file

外部ファイル読み込みの設定を行うブロック。

filetype

methodの値がfileだった場合に座標データを読み込むファイルの形式を選択する。以下が利用できる。

phase0_input : F_INP形式

phase0_output : F_DYNM形式

デフォルト値はphase0_input。

frame

F_DYNM形式の場合のフレーム番号を選ぶ。0以下の値の場合最後のフレームが採用される。デフォルト値は-1。

filetypeがphase0_outputだった場合、外部ファイルから取得できる情報は原子座標（と場合によっては原子速度）のみです。その他の原子配置に与えられる属性値(mobileなど)は別途F_INP形式のファイルから読み込む必要があります。これは、デフォルトの振る舞いではF_INPそのものになります。後述のように、ファイルポインターF_COORD_ATTRによって別のF_INP形式のファイルを指定することも可能です。

分子動力学シミュレーションのF_DYNMファイルには各原子の速度が記録されており、これを初期速度として採用することも可能です。この場合、以下のような設定を施すことによって初期速度をプログラムが割り当てることを抑制する必要があります。

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

file_names.dataファイルにおいて座標を読み込むファイルの指定を行います。関連するファイルポインターは下記の通りです。

ファイルポインター

説明

F_POS

原子座標データファイルを読み込むファイル。

デフォルト値は、fietypeがphase0_inputである場合F_INPが指す値、phase0_outputである場合F_DYNMが指す値。

F_COORD_ATTR

filetypeがphase0_outputである場合に原子の属性値を読み込むファイルの指定。 デフォルト値はF_INPが指す値。

file{filetype=phase0_output}としたとき、file_names.dataにF_POSの設定がないとF_DYNMが指すファイルから 原子座標を読み込みますが、有限温度計算や構造緩和を行った場合、F_DYNMデータを上書きしてしまうので、注意が必要です。 上書きを避けるためには、原子座標を読み込むファイルをF_POSで（F_DYNMとは別ファイル名を）与える必要があります。

4.6. 波動関数ソルバー（Wavefunction_Solver）

4.6.1. PHASE における計算フロー

PHASE における計算フローを 図 4.3 に示します。

../_images/input_param_fig1.svg

図 4.3 PHASE における計算フロー

波動関数の更新の過程で、Kohn-Sham方程式

(4.3)\[(H_{\text{KS}}(\rho_{\text{inp}}) - \epsilon_{i})\psi_{i} = 0,\]

を解いています。

ある試行の波動関数が与えられ、

(4.4)\[\Delta_{i} = (H_{\text{KS}}(\rho_{\text{inp}}) - \epsilon_{i})\psi_{i}\]

の演算を繰り返し行うことに より、 (4.3) の解が得られることになります。 その際、\(\Delta_{i}\) はエネルギー\(\epsilon\) の波動関数 \(\psi\) に対する勾配と解釈する事 ができるので、この勾配が0に近づくように波動関数が更新されます。 フローチャート 図 4.3 中の電荷密度の作成の過程では、 更新した波動関数から、新しい電荷密度 \(\rho\) が、以下の処方で与えられます。

(4.5)\[\rho_{\text{out}} = 2\sum_{\text{occ.}}^{}|\psi_{i}|^{2},\]

図中の内側のループでは、入力の\(\rho_{\text{inp}}\)と新しい\(\rho_{\text{out}}\)が 一致するまで計算が行われます。 この作業は SCF(自己無動着場)計算と呼ばれています。 外側のループでは、与えられた原子配置に対して力の計算が行われ、この力が0(閾値以下)になるような 原子配置に到達するまで、計算が続行されます。

4.6.2. 波動関数ソルバー

SCF計算において、「波動関数ソルバー」によって波動関数の更新を繰り返し行います。

波動関数ソルバーは、wavefunction_solverブロックで指定します。

wavefunction_solver{
  solvers{
    #tag sol      till_n prec cmix submat
      pdavidson 2     on 1    on
      rmm3     -1     on 1    on
  }
  rmm{
    edelta_change_to_rmm = 1e-3
  }
}

wavefunction_solverブロックで利用できるブロック/変数は以下の通りです。

solvers

どの波動関数ソルバーを、

どのタイミングで適用するかを設定する重要なテーブルです。以下の属性値を利用することができます。

sol

利用するソルバーのアルゴリズムを指定します。msd, lm+msd, cg, pkosugi, pdavidson, rmm3 から選択します。 msdは、修正再急降下法です。 1SCFあたりの計算負荷は最も軽い手法ですが、この手法のみで収束を得るのは困難です。主に初期波動関数に対して利用します。lm+msdはmsd法に一次元探索を備えた手法です。1回あたりの計算負荷は比較的軽い手法で、msd法と比較して収束の速い手法です。cg法は共役勾配法です。 計算負荷はlm+msdよりも多いですが、収束性は良い場合が多いです。davidson法は、1回あたりの計算負荷は高いですが信頼性の高い手法です。pkosugi, pdavidsonはdavidson法を並列化した手法であり、通常はdavidson法よりも推奨されます。rmm3法は1回あたりの計算負荷はDavidson 法よりも軽く、信頼性もDavidson法に劣らない場合の多い手法です。ただし、ランダムな波動関数に適用すると正しい解へ収束しない場合があるので、rmm3法を利用する場合は他のソルバーである程度波動関数を収束させてから移行する設定にする必要があります。

till_n

SCF計算の何ステップ目までsolを適用するかを指定します。上の例では、till_n値は、pdavidsonは2なので2回目まで使用する設定となっています（ただし後述のedelta_change_to_rmmで指定する条件を満たさない限り2回目以降もpdavidsonが利用され続けます）。負の数を指定すると収束するまでそのソルバーを使い続けます。したがって、rmm3は最後まで利用される設定となっています。

prec

前処理の有無を真偽値で指定します。デフォルト値はonで、通常変更する必要はありません。

cmix

採用する電荷密度混合法を指定します。電荷密度混合法は電荷密度混合法 において解説します。

submat

部分空間対角化を行うかどうかを真偽値で指定します。デフォルト値はonで、通常変更する必要はありません。

davidson

davidson法の詳細な振る舞いを設定したい場合に利用するブロックです。以下の変数を利用することができます。

max_subspace_size

davidson法で利用する部分空間の最大サイズを指定します。デフォルト値はバンド数の4倍です。

ndavid

davidson法はすこしずつ部分空間を拡大しながら波動関数を更新しますが、その回数を指定する変数です。デフォルト値は5です。

rmm

rmm法の詳細な振る舞いを設定したい場合に利用するブロックです。

edelta_change_to_rmm

rmm法は、ある程度収束した波動関数に適用しないと正しく動作しない場合があります。そこで、ここで指定した値よりも全エネルギーがよくなった時点でrmm法へ移行します。デフォルト値は、1e-3/natm hartree；ここでnatmは原子数。

line_minimization

lm+msd法やcg法は1次元探索を行い、最適なきざみ幅をもとめます。その1次元探索の詳細設定を行うブロックです。

dt_lower_critical

1次元探索の下限の刻み幅を指定しま す。デフォルト値は0.1です。

dt_upper_critical

1次元探索の上限の刻み幅を指定しま す。デフォルト値は2.0です。

4.6.2.1. pkosugi, pdavidsonソルバーの高速化オプション (バージョン2023.01以降)

Pkosugiおよびpdavidsonソルバー利用の際、収束しているバンドは処理をスキップすることによって高速化を実現することができます。以下の要領でバンドが収束しているかどうかの判定の閾値を設定します。

wavefunction_solver{
  davidson{
    delta_eig_occup = 1e-8
    delta_eig_empty = 1e-4
  }
}

davidson ブロックの下の delta_eig_occup で占有準位に対する収束判定条件を、 delta_eig_empty で非占有準位に対する収束判定条件を設定することができます。デフォルト値は両方とも1e-8 hartreeですが、非占有準位は全エネルギーに直接貢献しないため、バンドギャップがある場合はその収束判定は甘くても問題ない場合が多いです。非占有準位の閾値を甘く設定することによって全体の収束性を高めることができるかもしれません。

4.7. 電荷密度混合法（Charge_Mixing）

4.7.1. 電荷密度混合法

SCF 計算において、前回のSCF ステップで得られた電荷密度を一定程度混合することによって計算を進行させます。ここでは、この“電荷密度の混合方法” について説明します。電荷密度混合法の設定は、例えば、下記のようにcharge mixing ブロックで指定します。

charge_mixing{
    mixing_methods{
        #tag method rmxs rmxe  prec istr nbmix
              pulay  0.4 0.4   on   3  15
    }
    charge_preconditioning{
        amix = 0.9
        bmix = -1
    }
}

charge mixing ブロックにおいては、以下のブロック/変数が利用できます。

mixing_methods

電荷密度混合法を指定するためのテーブルです。

電荷密度混合法はいくつでも定義することができます。実際に利用する混合法は、前節のsolversテーブルのcmix 属性値によって指定します。cmix属性値では、利用したい電荷密度混合法を1始まりの定義順で指定します。このテーブルは、以下の属性値をもちます。

no

正の整数。cmixと対応付けられる番号。 この指定がない場合、methodなどを指定する行に1から順にnoが付加されます。

method

電荷密度混合のアルゴリズムを選びます。simple, broyden2, pulay, DFP のいずれかを選択することができます。simpleは単純混合法です。broyden2はBroydenの2番目の方法、pulayはPulayによるRMM-DIIS法です。 DFPはDavidson-Fletcher- Powell法です。broyden2法、pulay法、DFP法は、いずれも準ニュートン法の一種です。mixing_methodsブロックを定義していない場合にはpulayが既定値になります。

rmxs

混合比パラメータrmxの初期値を0から1の間で指定します。既定値は0.5（mixing_methodsを定義してい ない場合の既定値は0.4）。

rmxe

混合比パラメータrmxの最終値を0から1の間で指定します。指定がない場合、rmxsが使われます。

itr

時間刻み幅を変化させるSCFの回数。既定値は100。

var

補間方法。選択肢:{linear, tanh}。既定値はlinear。

prec

前処理の有無を真偽値で指定します。通常onにします。また既定値もonです。

istr

method=simple以外のときに有効なパラメータです。0以上の整数を指定します（既定値は3）。 broyden2、pulay、あるいはDFPが採用されている場合でも、SCF回数がistr+1に達するまではsimple法を使って電荷密度を混合します（少なとも一回はsimple法の適用が必要）。

nbmix

broyden2、pulay、あるいはDFPを選択している場合、過去の電荷密度の履歴を利用します。その履歴SCF回数を指定します（既定値は15）。

charge_preconditioning

前処理の係数を設定します。 amix, bmixというパラメータを同名の変数によってこのブロックの下で設定することができますが、通常設定しなくても問題ありません。 複数の電荷密度混合法を使う場合でも共通したamix、bmixが使われます。

amix

0と1の間の実数を指定。既定値は1.0。

bmix

正の実数を指定。 負の値を指定した場合、あるいは未設定の場合（-1が既定値として採用される）、後述の式の中の \(G_0^2\) を0.63として処理。

prec=offとした場合（charge_preconditioningが無効な場合）、単純混合法が働いているSCF回（method=simple以外の場合も、SCF回数がistr+1に達するまでは単純混合法が働きます）には次の式に従って電荷密度を混合します。

\[\rho_{\rm new} \leftarrow \left(1-{\rm rmx} \right) \rho_{\rm old} + {\rm rmx} \times \rho_{\rm new}\]

補間方法の指定varがない場合（あるいはlinearが指定されている場合）は、 rmxsとrmxeの間を線形補間した値を混合比パラメータrmxとしますが、varにtanhが指定されている場合、SCFの回数 i （ここで i は原子座標や格子が変化したあとに１から数え直す数）におけるrmxは次の式に従って設定されます。

\[{\rm rmx} = {\rm rmxs} + \frac{1}{2} \left( {\rm rmxe} - {\rm rmxs} \right) \left( tanh( 10 \times \frac{i-1}{\rm itr} - 5 ) + 1 \right)\]
../_images/extract_rmxt.png

図 4.4 電荷密度混合比rmxの推移。設定はitr=30、rmxs=0.1、rmxe=0.8。

charge_preconditioningが有効の場合、以下の要領で \(G\) の成分ごとの混合比を変えます。

\[\begin{split}\rho_{\rm new}\left(G\right) \leftarrow \left(1-f\left(G\right)\right) \rho_{\rm old} \left(G\right) + f\left(G\right) \rho_{\rm new} \left( G \right), \\ f\left( G \right) = \frac{ {\rm rmx} \times {\rm amix} } {1+ \left( \frac{G_0^2}{G^2} \right)}, \\ G_0 = {\rm bmix} \times G_{\rm min}\end{split}\]

ここで、 \(G_{\rm min}\) は原点以外の最小の逆格子ベクトルの長さです。

charge mixingグロックは次のように表記することも可能です（複数の収束法を利用する場合）。 "*"は既定値を使うことを意味します。

charge_mixing{
    mixing_methods{
        #tag no method rmxs rmxe  var  itr istr prec nbmix
             1  simple 0.1  0.4  tanh  5   *    on   *
             2  pulay  0.4   *    *    *   0    on  10
    }
}

4.7.2. 収束を加速させるテクニック

ここでは、SCF計算がなかなか収束しない場合について試すことのできるテクニックを紹介します。

部分空間対角化

部分空間対角化の有効/無効の設定は、変数submat によって制御することができます。

wavefunction_solver{
  solvers{
    #tag sol till_n  dts  dte  itr  var     prec  cmix submat
        lmMSD -1 0.2  1.0  40 linear  on 1 on
  }
}

部分空間対角化の適用を、波動関数を更新する前に適用するか後に適用するかによって収束の振る舞いが変化します。 これは、特にRMM法を利用している場合に大きな影響を与えます。 デフォルトの振る舞いでは波動関数更新前に部分空間対角化が適用されますが、 波動関数更新後にする場合には以下のように変数before_renewalをoffとします（経験的には、before_renewal=onの方が収束性がよい場合が多いです）。

wavefunction_solver{
  solvers{
    #tag sol till_n  dts  dte  itr  var     prec  cmix submat
        lmMSD -1 0.2 1.0 40 linear on 1 on
    }
    submat{
      before_renewal=off
    }
}

また、部分空間対角化はバンド数が多い方がより有効に作用します。バンド数を増やせばそれだけ計算量も増えますが、この効果によって全体の計算時間は短くなる場合もあります。

SCF 計算をある回数で打ち切る方法

初期の原子配置が安定な原子配置から遠い場合、SCF計算を収束させるのに多くの繰り返し計算が必要となる場合があります。このような場合は、たとえ電子状態が充分に収束していなくとも構造最適化をすすめることによって結果的に正しい解へより少ない計算時間で到達することができる場合があります。そこで、入力の指定の収束条件を満たしていなくとも収束したとみなし、構造最適化を進める機能がPHASEには備わっています。この機能を利用するためには、controlブロックの下でmax_scf_iteration変数を設定します。

control{
   ...
   max_scf_iteration = 50
}

この例では、50 回のSCF 計算を行っても収束判定を満たせなかった場合、その時点で至っている電子状態を利用して原子間力を計算し、構造最適化を進行させます。

電荷密度の差の混合比を変更する方法

スピンを考慮している場合、電荷密度混合は全電荷とスピン電荷密度（アップスピンの電荷密度とダウンスピンの電荷密度の差）に分離して混合します。全電荷とスピン電荷の混合比をそれぞれ違う値に設定することが可能です。このような設定を行うには、下記の要領でspin_density_mixfactor変数を定義します。

charge_mixing{
  spin_density_mixfactor = 4
  mixing_methods{
    #tag   no  method rmxs rmxe prec istr nbmix update
        1   broyden2  0.1 0.1 on 3 15 renew
  }
}

この例の場合、spin density mixfactor は4 であり、電荷密度の差の混合比は0.1 × 4 = 0.4 という値が採用されます。全電荷とスピン電荷を混合するのではなくアップスピンの電荷密度とダウンスピンの電荷密度を直接混合する場合、以下の要領でsw recomposing 変数にoff を設定します。

charge_mixing{
  sw_recomposing = off
  ...
}

スピン電荷密度の混合に利用するアルゴリズムを変更する

スピン電荷密度に対して、強制的に単純混合法を採用することも可能です。このような設定は、以下のようにspin_density ブロックを作成し、sw_force_simple_mixing変数を定義しその値をon とします。

charge_mixing{
  sw_recomposing=on
  spin_density_mixfactor = 4
  mixing_methods{
    #tag no method rmxs rmxe prec istr nbmix update
          1   broyden2  0.1 0.1 on 3 15 renew
  }
  spin_density{
    sw_force_simple_mixing = on
  }
}

スピンを固定する方法

一定の間スピンを固定してSCF計算を行うと収束性が改善する場合があります。この設定は、下記の要領でstructureブロックの下にferromagnetic_stateブロックを作成し行います

structure{
  ...
  ferromagnetic_state{
    sw_fix_total_spin = on
    spin_fix_period = INITIALLY
    total_spin = 1.0
  }
  ...
}

ferromagnetic_stateブロックでは以下の変数を利用することができます。

sw_fix_total_spin

"on"とした場合、スピンを固定した計算を行います

spin_fix_period

スピン固定の方法を指定します。

"INITIALLY"を指定した場合、SCF計算の初期は固定し、すこしずつ拘束を外していきます。"WHOLE"と指定した場合、計算終了までスピンを固定します。整数を指定した場合、その回数だけ固定しあとは通常の計算を行います。

total_spin

アップスピンとダウンスピンの差を指定します。単位胞全体の値を指定してください。

欠損電荷を混合する方法

PAW法を利用している場合、欠損電荷の混合が行われます。DFT+U法を利用している場合、占有行列の混合が行われますが、これも実質上は欠損電荷の混合をおこなっていることと同等です。この混合に対して通常の電荷密度と同様のアルゴリズムで混合させるには、以下のようにcharge_mixingブロックにsw_mix_charge_hardpart変数を定義し、その値をonにします。

charge_density{
  ...
  sw_mix_charge_hardpart = on
  ...
}

このように設定することによって、PAW法およびDFT+U法利用時の収束性が向上する場合があります。PAW法とDFT+U法を利用する場合は、このパラメータのデフォルト値はonです。

4.8. 波動関数ソルバーおよび電荷密度混合法の自動設定

PHASEに搭載されている波動関数ソルバーには、MSD法、lm+MSD法、Davidson法、CG法、RMM法、直接対角化法などの基本ソルバーと補助ソルバーとしてのsubspace rotationがあります。さらに、電荷密度混合法として単純混合法、Pulay法、Broydenによる2番目の方法などを搭載しています。これらを、問題に応じて適切に組み合わせることによって高速な収束が期待できます。しかし、このように問題に応じて適切に組み合わせるのは非常に手間がかかる作業です。そこで、PHASEには、適切な波動関数ソルバーや電荷密度混合法をプログラムが自動的に選択する機能があります。この機能は、様々な系に対し収束させることができるようになっていますが、もしなかなか収束させられない場合は、手動で波動関数ソルバーや電荷密度ミキサーの設定を行ってください。

「ソルバーセット」は、利用したい計算機能や並列数、バンド数などに応じて自動的に適切なものが採用される仕組みになっているので、利用にあたって特に気にする必要な項目はありません。この自動選択機能は, 波動関数に関してはwavefunction_solverブロックの下のsolversブロックが、電荷密度に関してはcharge_mixingブロックの下のmixing_methodsブロックが存在しない場合に有効となるので、本機能を利用したい場合は上述の設定を削除するかコメントアウトしてください。wavefunction_solverブロックは存在していても構わないので、ソルバーに関する詳細設定が必要な場合は対応するサブブロックにおいて行います。たとえば、本機能を利用しつつrmmソルバーは収束が10-6 hartreeよりもよくなったタイミングで利用したい場合は、以下のような記述を行います。

wavefunction_solver{
  rmm{
    edelta_change_to_rmm = 1e-6 hartree
  }
}

また、本機能の関連機能として、利用したい電荷密度混合法が1種類のみの場合、以下の簡易表記が利用可能となっています。

charge_mixing{
  method = pulay
  rmx = 0.2
  istr = 4
  nbxmix = 10
}

各変数は、以下のような意味を持ちます。

method

電荷密度混合の手法を選択します。simple, broyden2, pulay

のいずれかが有効。simpleは単純混合法、broyden2法はBroydenによる2番目の方法、pulay法はPulayによるDIIS法です。デフォルト値はpulay.

rmx

混合比を指定します。デフォルト値は0.4 (スピンを考慮していない場合), 0.1 (スピンを考慮している場合)

istr

broyden2法ないしpulay法を採用している場合に、はじめ何回をsimple法で混合するかを指定します。デフォルト値は3 (スピンを考慮していない場合)、5 (スピンを考慮している場合)

nbxmix

broyden2法ないしpulay法を採用している場合に、電荷密度の履歴を保持しておく回数を指定します。デフォルト値は15 (スピンを考慮していない場合), 5 (スピンを考慮している場合).

4.9. 構造最適化 (Structure_evolution)

構造最適化、分子動力学法計算に関するパラメータは、structure_evolutionブロックで指定します。

4.9.1. 構造最適化

structure_evolution ブロックに、構造最適化の設定をします。

...
structure_evolution{
    method = quench
    dt = 50
    ...
}
...

method

構造緩和の方法を指定します。

構造緩和のオプションとして、quench (quenched MD法)、cg (CG法)、cg2法（改良CG法）、gdiis (GDIIS法), bfgs (BFGS法), fire (FIRE法), lbfgs (L-BFGS法) のいずれかが選べます。デフォルト値はbfgsです。

dt

構造緩和を行う際の時間刻みです。大きい方が早く収束へいたりますが、大きすぎると計算を正しく進行させることができなくなる場合があります。 quenchおよびgdiisの場合に意味を持ちます。デフォルト値は原子単位で100です。構造最適化の場合は、このパラメータは数値解法の都合で導入しているものであり、物理的な意味はありません。

GDIIS法, BFGS法, L-BFGS法は原子に働く力が大きい場合安定に計算できない場合があるので、力が大きい内はquenched MD法かCG法を利用し、 ある程度力が小さくなってからGDIIS/BFGS/L-BFGS法に切り替える、という動作をします。GDIIS/BFGS/L-BFGS法に切り替える前の最適化手法と切り替えの判定条件は、それぞれ変数initial_methodとc_forc2gdiisを利用して 次のように設定します.

structure_evolution{
    method = gdiis
    dt = 50
    gdiis{
        initial_method = cg
        c_forc2gdiis = 0.0025 hartree/bohr
    }
}
...

ブロック名は、GDIIS, BFGS, L-BFGS共通でgdiisです。デフォルト値はinitial_methodがcg, c_forc2gdiisが0.05 hartree/bohr です。

gdiis

GDIIS およびBFGS を選択する場合のタグ

initial_method

GDIIS (BFGS) へ移行する前に利用する最適化

アルゴリズム。選択肢：{ quench, cg, cg2, sd }, デフォルト値はcg2

gdiis_box_size

ここで指定するイオン座標更新回数分のデータをgdiis(bfgs)用配列に蓄える

gdiis_hownew

gdiis_box_sizeで指定した回数分のイオン座標のデータ配列を使い切った時の処理法の選択

選択肢：{anew, renew}, デフォルト値はrenew

c_forc2gdiis

GDIIS (BFGS) への切替えの判定条件

デフォルト値は0.05 (hartree/bohr)

4.9.1.1. limited-memory BFGS法による構造最適化 (バージョン2021.01以降)

概要

PHASE/0には構造最適化手法としてquench法、CG法、GDIIS法、BFGS法、FIRE法などが搭載されています。バージョン2021.01以降、文献 [Hjorth17] のlimited-memory BFGS法 (l-BFGS法)が選択肢に加わりました。 通常のl-BFGS法のほか、文献において提案されている前処理を施すこともできるようになっています。

l-BFGS法について

limited-memory BFGS法とは、BFGS法を効率化した最適化手法です。ヘッセ行列をメモリーに明示的に持つのではなく、iterationごとに構築していくことに特徴を持つ手法です。ヘッセ行列をあらわに持たないため、メモリー要求が小さいことが名称のlimited memoryの由来です。そのアルゴリズムの擬コードは、以下のように記述することができます。

(4.6)\[\begin{split}s_{k-1} = x_k - x_{k-1} \\ y_{k-1} = g_k - g_{k-1} \\ \rho_k = \frac{1}{s_k y_k} \\ q = g_k \\ {\rm do \quad }i=k-1, k-2, ..., k-m \\ \alpha_i = \rho_i s_i q \\ q = q-\alpha_i y_i \\ {\rm end do} \\ \gamma_k = \frac{s_{k-1} y_{k-1}}{y_{k-1} y_{k-1}} \\ H_k^0 = \gamma_k I \\ z = q H_k^0 \\ {\rm do \quad }i=k-m, k-m+1, ..., k-1 \\ \beta_i = \rho_i y_i z \\ z = z+s_i \left( \alpha_i - \beta_i \right) \\ {\rm end do} \\ x_{k+1} = x_k - \alpha_k z\end{split}\]

ここで \(x_k, g_k`は :math:`k\) 回目のiterationにおける座標データとエネルギーの座標微分、 \(H_k^0\) は初期ヘッセ行列の推定値に相当する行列です。 その行列要素は (4.6) のように \(H_k^0 = \gamma_k I\) と計算することが一般的ですが、これを以下のような前処理行列で置き換えることが文献では提案されています。

(4.7)\[\begin{split}P_{ij} = -\mu \exp \left( -A \left( \frac{r_{ij}}{r_{nn}}-1 \right) \right) , r_{ij} < r_{\rm cut}, \\ 0 , r_{ij} > r_{\rm cut}\end{split}\]

ここで \(A, \mu, r_{nn}, r_{\rm cut}\) は任意に設定できるパラメーターです。 文献によると \(A\) はあまり結果を左右しないようです。文献でも採用されている3という値をデフォルト値としています。 \(r_{nn}\) は最近接間距離の最大値であり、 原子配置から自動的に決まる値です。 \(r_{\rm cut}\) はカットオフ距離に相当します。典型的には \(r_{nn}\) の2倍という値にります。 最も結果を左右するパラメーターが \(\mu\) です。このパラメーターは小さい方が速い収束が見込めますが、小さくし過ぎるとロバスト性が損なわれてしまいます。 PHASE/0のデフォルトの振る舞いとしては、前処理行列の行列要素の大きさが結果としておおむね通常のl-BFGS法の \(\gamma_k\) と同程度の値になるような値が採用されます。

使い方

構造最適化の設定は、入力パラメーターファイルのstructure_evolutionブロックにおいて行います。l-BFGS法を用いる場合、典型的には下記のような指定になります。

structure_evolution{
  method = lbfgs
  lbfgs{
    c_iteration2gdiis=3
    c_forc2gdiis = 0.05
    gdiis_box_size = 6
    initial_method = cg2
    maxstep = 0.1
  }
}

タグmethodにlbfgsを指定するとl-BFGS法を利用することができます。Limited-memory BFGS法（および通常のBFGS法、GDIIS法）の詳細設定は、structure_evolutionブロックの下のlbfgsブロック(bfgs, gdiisも可)において行うことができます。lbfgsブロック(もしくはbfgs, gdiisブロック)において設定できる主なパラメーターは下記の通り。

表 4.7 l-BFGS法の詳細設定

タグ名

説明

c_iteration2gdiis

lbfgs法、bfgs法、gdiis法は計算の最初期から適用されるわけではなく、はじめは閾値を満たすまでは別の最適化手法が用いられます。閾値を満たしたあと、このタグの指定の回数経てからlbfgs法 (もしくはbfgs法、gdiis法)にに移行します。デフォルト値は3.

c_forc2gdiis

lbfgs法、bfgs法、gdiis法へ移行する際の閾値を力の単位で指定します。デフォルト値は0.05 hartree/bohr

gdiis_box_size

lbfgs法、bfgs法、gdiis法の履歴の大きさ。デフォルト値は6.

initial_method

lbfgs法、bfgs法、gdiis法へ移行する前の構造最適化手法を指定します。デフォルト値はcg2.

maxstep

lbfgs法において、最適化1回で原子が動ける距離の最大値を長さの単位で指定します。デフォルト値は0.1 bohr. この設定値はlbfgs法の場合のみ利用できます。

多くのタグ名がGDIIS法にちなんだものとなっているのは、本機能がGDIIS法の履歴活用手法を用いているためです。

Limited-memory BFGS法は、上述の通り前処理行列を作用させることができます。その設定は、下記の要領で行います。

structure_evolution{
  method = lbfgs
  ...
  sw_prec = on
  prec{
    A = 3.0
    mu = -1
  }
}

structure_evolutionブロックにおいてsw_prec = onとすると前処理が有効になります。前処理の詳細設定はprecブロックにおいて行います。設定できるパラメーターは下記の通り。

表 4.8 l-BFGS法の詳細設定

タグ名

説明

A

(4.7) のパラメーター \(A\) . デフォルト値は3.

mu

(4.7) のパラメーター \(\mu\) . 負の値を指定すると、履歴に応じた最適値が自動的に計算されます。正の値の場合はその値がそのまま利用されます。デフォルト値は-1.

例題

PHASE/0の構造最適化の例題は、 samples/structural_evolution 以下に用意されています。 すべての例題についてquench, bfgs, cg2, gdiisの例題が用意されていますが、ここにlbfgsディレクトリーを追加しました。 PHASE/0付属の例題を用いて、CG2, BFGS, l-BFGS法が収束するまでに要したionic iterationの回数を以下に報告します。

表 4.9 l-BFGS法の詳細設定

系

CG2

BFGS

l-BFGS

Si(001)面

84

102

63

SiO 2 結晶

15

15

11

TiO 2 結晶

21

31

17

ジクロロシクロヘキサン

77

60

45

この表から、おおむねl-BFGS法が高速であることが分かります。必ずよい結果が得られるとは限りませんが、多くの問題でここで見たような傾向が見られており、 l-BFGS法は有力なオプションであると言えます。

[Hjorth17]

Ask Hjorth Larsen et al. J. Phys.: Condens. Matter 29 (2017) 273002.

4.9.2. 分子動力学法計算

分子動力学法計算に関するパラメータは、structure_evolutionブロックで指定します。

structure_evolution{
    method = velocity_verlet
    dt = 100
}

method

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

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

temperature_control (Nosé-Hoover熱浴 による温度一定の分子動力学シミュレーション)

velocity_scaling(速度スケーリングによる温度一定の分子動力学シミュレーション)のいずれか

dt

時間刻みを指定する。

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

分子動力学の場合、構造最適化の場合と異なり物理的な意味のある量である。

thermostat

熱浴を定義するブロック。

temp

温度を指定する。

qmass

熱浴の質量を指定する。

tdamp

熱浴の"周期"を指定する。ここで指定された周期から逆算して熱浴の質量がもとまる。qmassとtdampが混在する場合、qmassの指定が優先される。

4.9.3. 構造更新時の電荷密度、波動関数の予測更新（収束性の向上）

PHASEには、構造最適化や分子動力学シミュレーションを行っている際に、波動関数や電荷密度を原子配置の変化に合わせて“補外”することによって収束性を向上させる機能が備わっています。補外は [Arias94] で紹介されている方法によって行っています。

この機能を利用するには、structure_evolutionブロックにpredictorブロックを作成し、そこで本機能に関する設定を行います。

structure_evolution{
  predictor{
    sw_charge_predictor = on
    sw_extrapolate_charge = on
    sw_wf_predictor = on
  }
}

predictorブロックで定義できる変数です。

predictor

sw_charge_predictor

電荷の予測を行うかどうかを指定するスイッチ。

デフォルト値はon

sw_extrapolate_charge

予測の際に電荷密度の補外を行うかどうかを指定するスイッチ。デフォルト値はon

sw_wf_predictor

波動関数の予測を行うかどうかを指定するス イッチ。デフォルト値はoff

また、printoutlevelに変数ipripredictorを定義し、その値を2以上にすると補外をする際に原子配置の予測の精度やなどの情報がログファイルに出力されます。

[Arias94]

  1. Arias, M. C. Payne and J. D. Joannopoulos, “Ab initio molecular-dynamics techniques extended to large-length-scale systems”, Physical Review B 45, 1538 (1992).

4.9.4. ストレステンソル計算

ストレステンソル計算を行うには、structure_evolutionブロックのstressブロックで指定します。

structure_evolution{
  stress{
    sw_stress=1
  }
}

stress

ストレス計算

sw_stress

ストレス計算の有無。選択肢：{ on,off }

4.10. 後処理（Postproccesing）

4.10.1. 状態密度（DOS）

SCF計算が収束したのち、状態密度の計算を行うことができます。電荷密度の計算を行うには、postprocessingブロックの下のdosブロックで設定します。

postprocessing{
    dos{
        sw_dos = on
        method = gaussian
        deltaE_dos = 1e-4 hartree
    }
}

dos ブロックでは以下の設定を行うことができます。

sw_dos

状態密度計算を行うかどうかを指定する真偽値です。

状態密度の計算を行う場合onとします。

method

状態密度の計算方法を指定します。gaussianとtetrahedralのいずれかを選択することができます。gaussianを選択した場合、エネルギー準位をガウス関数によって幅を持たせた上で計算した状態密度が得られます。tetrahedralの場合四面体法による高精度な状態密度計算を行うことができます。ただしtetrahedralを利用する場合後述の四面体法が利用できる条件もご参照ください。

deltaE_dos

状態密度計算に利用されるエネルギーの幅をハートリー単位で指定します。デフォルト値は1e-4 hartreeです。

状態密度の計算方法としてtetrahedralを利用する場合、以下の条件が満たされている必要があります。

  • k点サンプリング手法としてmesh 法を採用している

accuracy{
  ksampling{
    method = mesh
  }
}

  • smearing の方法としてtetrahedral 法を採用している

accuracy{
  smearing{
    method = tetrahedral
  }
}

以上が満たされていないとgaussian 法による状態密度計算が行われてしまうのでご注意ください。

4.10.2. 電荷密度

SCF計算中は逆空間で電荷密度を扱いますが、収束した電荷密度を実空間に逆フーリエ変換し、出力させることも可能です。こうすることによってPHASE-Viewerなどを利用して電荷密度の可視化を行うことが可能です。電荷密度を実空間に出力させるためには、postprocessingの下のchargeブロックで設定を行います。

postprocessing{
  charge{
    sw_charge_rspace = on
    filetype = cube
  }
}

chargeブロックの下では以下の変数の設定を行います。

sw_charge_rspace

電荷密度を実空間で出力するかどうかを指定する真偽値です。

onにすると実空間の電荷密度が出力されます。

filetype

電荷密度データのデータフォーマットを指定します。density_onlyとcubeが選べます。density_onlyの場合電荷密度のみが出力されます。 デフォルト値はdensity_onlyです。cubeの場合、Gaussian Cube形式で電荷密度が出力されます。このパラメータは、cubeに設定することを推奨します。

title

Gaussian Cubeファイルの"見出し"指定します。空白文字を含める場合、全体を半角の2重引用符で囲みます。

また、filetypeとしてcubeを選択した場合、file_names.dataファイルにおいて電荷密度ファイルのファイル名を変更しておくことを推奨します。

&fnames
...
F_CHR = './nfchr.cube'
/

変更しない場合のデフォルト値はnfchr.dataです。

スピン分極を考慮している場合は、file_names.dataで指定したファイル名がnfchr.cubeであったとすると、nfchr.up.cubeとnfchr.down.cubeという2つのファイルにそれぞれスピンアップ・ダウン に対応する電荷密度データが出力されます。

4.10.3. 構造最適化/分子動力学シミュレーションの最中に後処理を行う方法

ここで説明した後処理は、特に設定が施されていない場合力が収束した後に実行されます。構造最適化や分子動力学シミュレーションの最中に後処理を行いたい場合、以下のような記述を行います。

postprocessing{
  ...
  frequency = 5
}

変数 frequency に正の値を指定した場合、指定した回数に 1 回の頻度で後処理が行われるようになります。結果は、状態密度の場合は dos_iterxx.data ファイル、電荷密度は nfchr_iterxx.data ファイルに記録されます（xxは原子配置の更新回数に相当する数値に読み替えてください）。なお、この機能によって出力できるのは、状態密度と電荷密度のみです。

4.11. ログレベル（PrintLevel）

PHASEはoutput000というファイル（000は計算を行うたびに1ずつ増えます）にログを記録します。そのログの詳細度の設定はprintoutlevelブロックで行います。

printoutlevel{
  base = 1
}

最上位にprintoutlevelブロックを作成し、その下にログレベルを制御する変数を定義します。ログレベルを制御するための変数は 0,1,2のいずれかの値をとり、数字が大きいほどより詳細な出力が得られます。デフォルト値はすべて1です。 ログレベルを制御する変数として主なものは以下の通りです。

base

計算全体のログレベルを指定します。特に指定のない項目はここでの指定に従います。

timing

計算時間に関わるログレベルを制御します。

input

入力に関わるログレベルを制御します。

solver

波動関数ソルバーに関わるログレベルを制御します。

spg

空間群に関わるログレベルを制御します。

base=2に設定すると膨大な量の出力が得られ、ログファイルが見づらくなってしまいます。得られる情報のほとんどはデバッグ情報なので、特別な事情がない限りbase=2は指定しないことを推奨します。