従来のSTATEではGGA等で収束させた電荷密度と原子配置を用いてvdW-DFの非局所相関エネルギーを計算していましたが、STATE 5.5.4では直接vdW-DFを用いて電荷密度と原子配置の収束計算が行えます。オリジナルのvdW-DFでは二重積分の計算にO(N^2)の計算コストがかかり、特に大きな系への適用が困難でしたが、Román-PérezとSoler*1のアルゴリズムを用いると計算コストをO(NlogN)に抑えることが可能です。ここではRomán-PérezとSolerの方法を簡略化したWuとGygi*2の方法を用いて実装を行いました。
プログラムからRomán-Pérez-Solerアルゴリズム用のvdW kernelのテーブルを読み込むため
~/STATE/src/
等にvdwdphi.datをおいて、計算を行うディレクトリにvdwdphi.datのリンクを作って下さい。
GGAで電子密度と原子構造を求めたあとicond=26に変更して計算を行うと、 従来通りnon-self-consistentなvdW-DFの計算が行えます。 この場合、実空間で非局所相関エネルギーを計算するので計算コストがかかりますが、 入力ファイルの下に
&VDW-DF FFTVDW &END
を追加することで、Román-Pérez-Solerアルゴリズムを用いることで計算が高速化されます。
vdW-DF1以外のvdW密度汎関数を用いる場合には交換汎関数と非局所相関汎関数の種類を別々に指定します。 例えばvdW-DF1に対して交換汎関数を改良したoptB86b-vdWを用いる場合には、&VDW-DFのブロックを
&VDW-DF FFTVDW x optb86b &END
のように変更します。xの次の文字列は交換汎関数を表します。
さらに交換汎関数に加え非局所相関汎関数も変更したrev-vdW-DF2を用いる場合には
&VDW-DF FFTVDW VDW-DF2 x b86r &END
とします。 VDW-DF2を省略するとVDW-DF1の非局所相関汎関数が指定されます。 交換汎関数と非局所相関汎関数の組み合わせは下表を参考にして下さい。
STATE 5.5.4以降では交換汎関数と非局所相関汎関数の組み合わせを&VDW-DFのブロックで指定しなくても、 下表のvdW密度汎関数名をxctypeに指定すればPost-processing vdW-DFの計算が出来るように実装されていますが、 この方法では非局所相関エネルギーが常にvdW-DF1で計算されるバグが判明していますので、 vdW-DF2およびrev-vdW-DF2を用いる際には注意して下さい。 最新版のSTATE 5.6.0でもこのバグは修正されていないので、とりあえず上の指定方法を用いるのが無難です。
一方、icond=0(初めから) or 1(続きから)のとき、xctypeとして以下のいずれかを選択するとself-consisitentなvdW-DFの計算が行えます。 Self-consistentの場合は&VDW-DFのブロックでFFTVDWを指定する必要はありません。 表中で非局所相関汎関数のvdW-DF1とvdW-DF2の違いは、パラメータZabが前者では-0.8491で後者では-0.8491x2.222=-1.887となる点です。
xctype | 非局所相関 | 交換 | 文献 |
vdw-df または drsll | vdW-DF1 | revPBE | M. Dion et al., Phys. Rev. Lett. 92, 246401 (2004) |
vdw-df2 または lmkll | vdW-DF2 | PW86R | K. Lee et al., Phys. Rev. B 82, 081101(R) (2010) |
c09 または c09-vdw または drsllc | vdW-DF1 | C09x | V. R. Cooper, Phys. Rev. B 81, 161104(R) (2010) |
c09-vdw2 または lmkllc | vdW-DF2 | C09x | |
optpbe または optpbe-vdw | vdW-DF1 | optPBE | J. Klimeš et al., J. Phys.: Cond. Matt. 22, 022201 (2010) |
optb88 または optb88-vdw または kbm | vdW-DF1 | optB88 | J. Klimeš et al., J. Phys.: Cond. Matt. 22, 022201 (2010) |
optb86b または optb86b-vdw | vdW-DF1 | optB86b | J. Klimeš et al., Phys. Rev. B 83, 195131 (2011) |
rev-vdw-df2 または lmkllh | vdW-DF2 | B86R | I. Hamada, Phys. Rev. B 89, 121103(R) (2014) |
STATE 5.5.4および5.6.0ではB86交換ポテンシャルの計算のバグ修正が反映されておりませんでしたので、 optB86b-vdWおよびrev-vdW-DF2を用いる際には修正版を入手して下さい。
Román-Pérez-SolerアルゴリズムではvdW kernelをq0に関して分割するため、想定するq0の最大値(qcut)と分割数(nqs)を指定する必要があります。 例えばbulkの軽元素の場合、q0の最大値は3程度なので
&VDW-DF qcut 5.d0 nqs 10 &END
とします。一方、bulkの遷移金属元素ではq0の最大値が7程度になるので
&VDW-DF qcut 10.d0 nqs 20 &END
とします。軽元素であっても孤立系や表面系ではq0の分布が長いtailを持つようになるので、後者の設定が無難のようです。
分割数nqsは大きければ大きいほど良いというわけでもなく、kernelのFourier成分を格納する配列はnqsの自乗に比例して増大するため、 特に大きな系では計算機のメモリに応じてnqsの大きさを制限する必要があります。 経験的には原子数50個程度の系でnqsを20より大きくしても 全エネルギーの変化は数meV程度に収まるようです。
何も指定しない場合はqcut=3.d0、nqs=10が設定されますが、次期バージョンではqcut=10.d0、nqs=20をデフォルトにする予定です。