3. 入力パラメータファイル:nfinp.data(詳細版)

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

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

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

3.1.1. パラメータ設定形式

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

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

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

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

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

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

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

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

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

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

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

control ブロック

全体的な計算条件の設定

accracy ブロック

計算精度の設定

structure ブロック

原子構造の設定

wavefunction_solver ブロック

波動関数ソルバーの設定

charge_mixing ブロック

電荷密度混合法の設定

structure_evolution ブロック

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

postproccesingブロック

後処理の設定

printlevelブロック

ログ出力の設定

3.1.2. 単位の指定

PHASE の入力ファイルのデフォルトの単位は原子単位ですが, 単位を明示的に指定することも可能です。

表 3.1 の単位を利用することができます(デフォルトの単位は太字で表示されています)。

表 3.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 の長さの単位のデフォルトがÅ単位になります。 複数指定する場合(長さ、エネルギーなど), 空白で区切って指定してください。

3.1.3. コメント

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

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

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

3.1.4. 二項演算子を用いた実数指定(バージョン2020.01以降)

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

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

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

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

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

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"
  }
}

3.2. 入力パラメータファイル nfinp.dataのタグ(キーワード)の一覧

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

表 3.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ステップ内のS CF回数の制限値

(デフォ ルトは無制限)

nfstopcheck

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

sw_ekzaj

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

デフ ォルト値はOFF です。

accuracy

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

cutoff_wf

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

cutoff_cd

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

num_bands

バンド数

ksampling

method

k 点のサ ンプリング法。

monk: Monkhorst-Pack 法。

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

file:フ ァイルから入力

direc t_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 法( 金属系で有効)

tetrahedron:Tetrahedron法

improved_tetrahedron : 改 良tetrahedron 法

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

width

smearing 幅(デフォ ルト値:0.001 hartree)

method = parabolic とcold の時に使用

(タグなし)

xctype

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

LDA : LDAPW91, PZ

GGA : GGAPBE, REVPBE

scf_convergence

自己無憧着場の収束判定

delta_total_energy

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

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

succession

誤差が \(\Delta E\) 何回連続して収 まった時に計算 を停止させるか を指定する回数

(デフ ォルト値: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 タグなし)

magnetic_state

入力値: {para, antiferro, ferro} から選択

antiferro はaf と省略も可

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

波動関数ソルバ ー(詳しくは 3.6.2 章3.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のときsub space_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/natm hartree;ここ で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_rotati onを行います。

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

critical_ratio

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

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

charge_mixing

電荷混合法 。(詳しくは 3.7 章3.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

出力時のエネルギー幅ΔEを次式で指定: ΔE=nwd_window_width × 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 に設 定するとゼロ番 ノード以外のプ ロセスからもo utput00x_xxxと いったファイル に出力を行う。

jobstatus

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

jobstatus_option

状況ファ イルの出力制御

jobstatus_format

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

jobstatus_series

ON またはOFF

3.3. 全体的な計算条件設定(Control)

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

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

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

condition

初期計算か継続計算かなどの指定です。"initial" とすると計算は始めから行われ,"continuation" とすると,波動関数,電荷密度分布などの計算結果を引き継い だ継続計算が行われます。また,"automatic" とすると継続計算に必要な複数のファイル(これらは、前のジョブが正常終了した場合には自動的に生成される)が存在する場合は"continuation",存在しない場合は"initial"と設定したのと同じ動作をします。"preparation" を指定すると,前処理(使用配列の大きさの評価,k点生成など )のみ行います。デフォルト値は"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"です。継続計算になる可能性がある場合には、ジョブの指定時間よりも小さな値に設定しておくのがよいでしょう(例えば,ジョブの制限値が6時間で,収束後の状態密度計算などの後処理がなければ,5.8hour程度に設定する)。

max_iteration

max_total_scf_iteration

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

max_scf_iteration

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

3.4. 計算精度の指定(Accuracy)

3.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以降はこのような値がデフォルト値になりました。

3.4.2. バンド数

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

accuracy{
   num_bands = 12
}

num_bands

バンド数

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

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

3.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点を入力する必要がある場合やフェルミ面の計算において大量のk点を考慮する必要がある場合などに利用します。gammaを指定すると\(\Gamma\)点のみをサンプリングします。充分大きな単位胞を使っていて,\(\Gamma\)点のみでも充分な精度が得られる場合には,これを指定します。directinは直接k点の組(個数と座標)を指定します。いずれの方法でも、サンプリングk点に\(\Gamma\)点が含まれていて,系に反転対称中心がなければ(設定されていなければ),\(\Gamma\)点の波動関数に関する計算は,この点の対称性を利用して他のk点のものに比べて3倍程度高速に実行されます(後述のとおり,これを抑制する,つまり他のk点と同じ演算法を適用する手段もあります)。

mesh

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

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

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

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

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

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

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

method

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

width

スメアリングの幅をエネルギーの単位で指定します。デフォルト値は0.001 hartreeです。

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

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

../_images/parabolic.svg

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

3.4.3.2. バージョン2020.01以降の拡張

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

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

accuracy{
  ksampling{
    density = 4 bohr
  }
}

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

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

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

3.4.4. 交換相関エネルギー

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

accuracy{
  xctype = ggapbe
}

xctype

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

LDA : LDAPW91, PZ

GGA : GGAPBE, RPBE, REVPBE

3.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です。

3.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

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

3.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です。

3.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
    1. 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).

3.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
    }
}

3.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

3.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となっている必要があります。

3.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エネルギーはこの余剰電荷を打ち消すように背景に一様電荷分布をおいた状態にして計算します。

3.5.4. 原子種の指定(バージョン2020.01以降)

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

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

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

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

例:additional_charge の使用例

structure{
  charged_state{
    additional_charge = -1.0d0
  }
}

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

原子種の初期磁気モーメント (分極) については、 3.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
  }
}

3.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を参照してください。この変数を利用することによって,入力ファイルでは指定のしやすいブラベー格子で単位胞を指定しつつ,実際の計算は負荷の少ない基本格子で実行することが可能となります。lattice_system変数を利用して格子を変換させる場合,変換されるのは単位胞のみです。したがって,原子配置の定義などは基本格子の場合の定義の仕方を採用してください。たとえば,面心立方格子の計算を行う場合面心位置の原子は指定せず,原点の原子のみ指定するようにしてください。また,\(k\)点サンプリングは変換後の格子に合わせて設定してください。

generators

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

3.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とは別ファイル名を)与える必要があります。

3.6. 波動関数ソルバー(Wavefunction_Solver)

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

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

../_images/input_param_fig1.svg

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

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

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

を解いています。

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

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

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

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

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

3.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です。

3.7. 電荷密度混合法(Charge_Mixing)

3.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 は原子座標や格子が変化したあとに1から数え直す数)における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

図 3.3 電荷密度混合比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
    }
}

3.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です。

3.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 (スピンを考慮している場合).

3.9. 構造最適化 (Structure_evolution)

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

3.9.1. 構造最適化

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

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

method

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

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

dt

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

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

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

ブロック名は,GDIIS, 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)

3.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の指定が優先される。

3.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).

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

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

structure_evolution{
  stress{
    sw_stress=1
  }
}

stress

ストレス計算

sw_stress

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

3.10. 後処理(Postproccesing)

3.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 法による状態密度計算が行われてしまうのでご注意ください。

3.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つのファイルにそれぞれスピンアップ・ダウン に対応する電荷密度データが出力されます。

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

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

postprocessing{
  ...
  frequency = 5
}

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

3.11. ログレベル(PrintLevel)

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

printoutlevel{
  base = 1
}

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

base

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

timing

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

input

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

solver

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

spg

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

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