PCで苔を育てる人

自作シミュレーションゲームPraparatを作っています。 人工生命をシミュレーションするゲームです。https://www.nicovideo.jp/watch/sm41192001

【夏休みの自由研究】LK-99のバンド計算をしてみよう！

科学

巷では常温常圧超伝導が話題ですね。

きっかけになったのは、7月22日にarXivに投稿された下記の論文です。

arxiv.org

常温常圧超伝導は、人類の悲願とも言える存在であったため、 これだけの注目が集まるのも無理はないかと思います。 一方で、𝕏(旧:Twitter)を始めとするSNSでは、特に専門家以外の方々の間で 熱狂的に議論されているように見えるのは興味深いです。



LK-99の理論について

このLK-99については、実験だけではなく理論の方面からの研究も行われています。 そのうちの1つが密度汎関数法(DFT)を使ったバンド計算です。 すでにいくつかのグループが報告していますが、フェルミ面付近にフラットなバンドが見えており、 これが超伝導体にも見られるものらしいのですが、私は専門外なのでよく分かっていません。

重要なことは、フェルミ面付近でフラットなバンドを持つことは、 超伝導体であるための十分条件ではないということです。

  • LK-99が超伝導体であることが理論的に示された ← 間違い
  • ❍ 超伝導体の特徴の1つであるフラットバンドをLK-99が有することが分かった



LK-99のバンドを計算してみる

さて、ここからが本題ですが、ここでは、LK-99のフラットなバンドを 手持ちのパソコンで計算してみることにします。

理論の詳細には立ち入らず、とりあえずこの記事の通りに作業すれば、 LK-99のバンドが計算できるというところを目指します。

夏休みに入ったばかりの大学生が、自由研究感覚で挑戦できると良いなと思っています。 理屈は分からずとも、今話題になっている論文と同じ結果を、家庭のパソコンでも 出力できるという経験は、印象深いものになるのではないでしょうか。

最初に断っておきますが、私は専門家でもなんでもないので、 内容に踏み込んだ議論は行いません。また、ここでの計算結果にも責任は持てません。 このブログは私の自由研究でもあるわけです。


参考にする論文

今回は、こちらの論文を参考にしたいと思います。

arxiv.org

理論の論文をいくつか見ましたが、計算条件が一番分かりやすかったのがこれでした。


DFTソフトのインストール

まずは必要なソフトのインストールです。

論文を見てみると次のように書いてあります。

The first-principles calculations were performed using the Vienna ab initio simulation package (VASP).

ということで筆者らはVASPというソフトを使って計算を行っています。

VASPのホームページにアクセスすると、どうやらライセンスの購入が必要なようです。 価格は書いておらず、まずはお問い合わせからとのこと。



値段が分からないままお問い合わせをするのは怖いので、他のサイトで確認します。

大学生のお小遣いで買える程度だとよいのですが...


VASPの価格 (上杉 著「第一原理シミュレーションによる添加元素の最適化設計」より)


はい...？

ちょっと、お小遣いでは厳しい額ですね...。 しかも教育機関向けと書いてあるので、個人では買えないのでしょう。

というわけで別のソフトを検討する必要が有ります。

幸いなことに、フリーで使えるオープンソースの計算ソフトがいくつか存在します。

今回は、その中でも有名なQuantum Espresso(QE)を使ってみたいと思います。

github.com

インストール自体は、下記のとおりに実行すれば難なく終わります。

% pwd
/home/username
% sudo apt install build-essential gfortran
% sudo apt install openmpi-bin
% git clone https://github.com/QEF/q-e.git
% cd q-e
% ./configure
% make all

ただ、私が検証に利用したパソコンは、すでに色んなライブラリを インストールしてしまっていたため、もしかすると上記だけでは 足りないかも知れません。

その際は、お手数ですがご連絡いただけますと、大変助かります。


必要なファイルの用意

次は計算に必要なファイルを用意します。

下記のようにフォルダを用意しておきましょう。

% pwd
/home/username
% mkdir lk99
% cd lk99


擬ポテンシャル(Pseudopotential)

まずは、擬ポテンシャル(Pseudopotential)を用意します。 これは、電子が原子核から受けるポテンシャルだと思ってもらえれば良いです。

下記のサイトからダウンロードできます。

pseudopotentials.quantum-espresso.org

LK-99は、 {\rm Pb_9 Cu (PO_4)_6 O }という組成をしているので、 Pb, Cu, P, Oの擬ポテンシャルを用意する必要があります。

今回は、下記のコマンドでダウンロードしてください。

% pwd
/home/username/lk99
% mkdir pp
% cd pp
% wget https://pseudopotentials.quantum-espresso.org/upf_files/Pb.pbe-dn-kjpaw_psl.1.0.0.UPF
% wget https://pseudopotentials.quantum-espresso.org/upf_files/Cu.pbe-spn-kjpaw_psl.1.0.0.UPF
% wget https://pseudopotentials.quantum-espresso.org/upf_files/P.pbe-nl-kjpaw_psl.1.0.0.UPF
% wget https://pseudopotentials.quantum-espresso.org/upf_files/O.pbe-n-kjpaw_psl.1.0.0.UPF


結晶構造と計算条件

次はLK-99の構造ファイルを入手します。 幸いなことに、これは論文のSupplementary Materialの中で与えられているので、 これをそのまま使わせていただきます。

論文の後ろの方にある"Pb9Cu(PO4)6O"という文字から始まるセクションを コピーして、POSCARというファイル名で保存してください。 置く場所はどこでも良いです。

▼POSCARの例

% cat POSCAR
Pb9Cu(PO4)6O
1.00000000000000
9.9306823821051555 -0.0000014818366361 0.0000000000000000
-4.9653399077953200 8.6002239607356703 -0.0000000000000000
0.0000000000000000 0.0000000000000000 7.4109945172796605
Pb Cu P O
9 1 6 25
Direct
0.9986807449642078 0.7699202447216552 0.2474756380877509
0.9977014382162520 0.2577003216124319 0.7545056960306785
0.2300797552783452 0.2287605002425670 0.2474756380877509
0.7422996783875685 0.7400011176038168 0.7545056960306785
0.7712395147574378 0.0013192690357936 0.2474756380877509
0.2599988973961879 0.0022985757837491 0.7545056960306785
0.3333333129999971 0.6666666269999979 0.0102252768912011
0.3333333129999971 0.6666666269999979 0.4962563424132084
0.6666666870000029 0.3333333429999996 0.5217140548665616
0.6666666870000029 0.3333333429999996 0.0634955699714910
0.6238378061607577 0.5940989414527301 0.2330614375419467
0.3715115645822017 0.3913852574380640 0.7499251076346213
0.4059010585472700 0.0297388347080183 0.2330614375419467
0.6086147425619355 0.9801263371441329 0.7499251076346213
0.9702611952919840 0.3761622238392447 0.2330614375419467
0.0198736928558696 0.6284884654178010 0.7499251076346213
0.4970341950086190 0.6419127807282180 0.2504819573050943
0.4722574174374839 0.3124116086150668 0.7489819492636393
0.3580872192717823 0.8551214142804011 0.2504819573050943
0.6875883913849331 0.1598458078224275 0.7489819492636393
0.1448785557195963 0.5029657759913753 0.2504819573050943
0.8401541621775773 0.5277425535625101 0.7489819492636393
0.7465890784344098 0.6974373665162791 0.0877326196589775
0.2518710331486960 0.3302579489920529 0.9082269759127742
0.3025626334837207 0.0491516819181211 0.0877326196589775
0.6697420510079468 0.9216131131566275 0.9082269759127742
0.9508483470818850 0.2534109515655927 0.0877326196589775
0.0783869158433784 0.7481289968513065 0.9082269759127742
0.2751067320564097 0.3631086026138700 0.5737407252403873
0.7132743878686607 0.6201447197159712 0.4140528338576862
0.6368913973861298 0.9119981584425314 0.5737407252403873
0.3798552802840290 0.0931296381526946 0.4140528338576862
0.0880018705574747 0.7248932979435926 0.5737407252403873
0.9068703908473114 0.2867256421313418 0.4140528338576862
0.5441603686812513 0.4191736441937305 0.1764358135764958
0.4780900223868556 0.5733419996779769 0.7722935441387523
0.5808263258062671 0.1249867244875282 0.1764358135764958
0.4266579703220141 0.9047479627088743 0.7722935441387523
0.8750132455124694 0.4558396313187486 0.1764358135764958
0.0952520072911234 0.5219099776131445 0.7722935441387523
-0.0000000000000000 0.0000000000000000 0.2917658511111268

正直、この構造ファイルさえ手に入れば、半分は勝ったも同然です。

ただし、このファイルはVASP用のフォーマットになっているので、 これをQE用に変換してあげる必要があります。

下記のサイトにアクセスしてください。

www.materialscloud.org

フォーマットをPOSCARにし、「ファイルを選択」で先程作成した POSCARファイルをアップロード、"Calculate my structure"をクリックしてください。

すると良い感じの画面が出てくるので、下の方にスクロールしてQEのインプットを見つけてください。

これをクリップボードにコピーして、下記のような感じで保存してください。

% pwd
/home/username/lk99
% mkdir inputs
% cd inputs
% cat > input.in

余談ですが、PythonのASEやPymatgenなどのライブラリを使えば、 こういった作業は自動でこなせるようになります。


実行

それではいよいよ計算の開始です。 DFTによるバンド計算の流れは、概ね次のようなものになります。

  1. 構造最適化
  2. SCF
  3. NSCF
  4. DOS
  5. バンド図

このうち、1の構造最適化が一番時間がかかるのですが、 今回はこれをスキップして、論文で与えられている構造を そのまま使うことにします。 また、4のDOSについても今回は行いません。

  1. 構造最適化
  2. SCF
  3. NSCF
  4. DOS
  5. バンド図

ちなみに、今回私が6コアのマシンで検証した際の計算時間は下記のとおりでした。

  1. 構造最適化
  2. SCF (1時間39分)
  3. NSCF (3時間5分)
  4. DOS
  5. バンド図 (4時間4分)


SCF計算

最初にSCF計算を行います。

まずは次のように作業してください。

% pwd
/home/username/lk99
% mkdir 01scf
% cd 01scf
% cp ../inputs/input.in .

ここで、input.inを次のように編集します。

▼input.inの例

% cat scf.in
&CONTROL
    prefix      = 'lk99'
    calculation = 'scf'
    tstress     = .true.
    tprnfor     = .true.
    outdir      = './out'
    pseudo_dir  = '../pp'
/
&SYSTEM
    ibrav       = 0
    nat         = 41
    ntyp        = 4
    ecutwfc     = 40
    ecutrho     = 160
    occupations = 'smearing'
    degauss     = 0.003675
    smearing    = 'gaussian'
/
&ELECTRONS
    electron_maxstep = 100
    conv_thr         = 1e-09
    mixing_beta      = 0.7
    diago_david_ndim = 4
/
ATOMIC_SPECIES
Pb         207.2 Pb.pbe-dn-kjpaw_psl.1.0.0.UPF
Cu        63.546 Cu.pbe-spn-kjpaw_psl.1.0.0.UPF
P   30.973761998 P.pbe-nl-kjpaw_psl.1.0.0.UPF
O         15.999 O.pbe-n-kjpaw_psl.1.0.0.UPF
ATOMIC_POSITIONS angstrom
Pb       1.1293232294     3.7547448097     1.8340405970
Pb      -1.3023964628     7.9497623531     5.5916375765
Pb       1.1489751864     7.7008734208     1.8340405970
Pb      -1.2681567294     3.4974336207     5.5916375765
Pb      -2.2782984159     5.7448282092     1.8340405970
Pb       2.5705531922     5.7532504658     5.5916375765
Pb       4.9653411911     2.8667410733     0.0757794710
Pb       4.9653411911     2.8667410733     3.6777530328
Pb       0.0000000000     0.0000000000     3.8664200002
Cu       0.0000000000     0.0000000000     0.4705653209
P        8.2105722025     2.2426423547     1.7272170358
P        6.7113430460     0.4992594225     5.5576908611
P        3.8832110300     5.9892429296     1.7272170358
P        6.1426395162     5.5625638601     5.5576908611
P        2.8022403406     0.3683379355     1.7272170358
P        2.0420410111     2.5383999373     5.5576908611
O        6.7139137896     2.6538520455     1.8563204123
O        3.1386080324     8.4202916340     5.5507011195
O        4.2754221980     4.4874938779     1.8563204123
O       -3.8961492874     7.1081916914     5.5507011195
O        3.9066875856     1.4588772965     1.8563204123
O        0.7575412550     1.6719631142     5.5507011195
O       -1.0142168212     3.1313758775     0.6501859633
O        0.8614076797     8.5737741449     6.7308651389
O        2.7605985430     6.1561977489     0.6501859633
O       -2.8904688646     5.0593370810     6.7308651389
O       -1.7463817217     7.9128728132     0.6501859633
O        2.0290611848     3.5673352137     6.7308651389
O        5.8943805130     0.2560738791     4.2519893691
O       -0.9612700690     2.4666418617     3.0685432816
O        6.7617256411     4.9766463243     4.2519893691
O        3.3098017112     6.5344179893     3.0685432816
O        2.2399174191     3.3675030165     4.2519893691
O       -2.3485316423     8.1993865886     3.0685432816
O        8.2878846617     0.7382457515     1.3075648471
O        6.8662625520     2.0641280221     5.7234632213
O        5.1474004762     6.8083957849     1.3075648471
O        4.7099638023     4.9142937880     5.7234632213
O        1.4607384353     1.0535816834     1.3075648471
O        3.3197972189     1.6218014097     5.7234632213
O        0.0000000000     5.7334821466     2.1622751229
K_POINTS {automatic}
4 4 5 0 0 0 
CELL_PARAMETERS angstrom
    9.9306823821     0.0000000000     0.0000000000
   -4.9653411911     8.6002232198     0.0000000000
    0.0000000000     0.0000000000     7.4109945173

次のコマンドで実行です。

% pwd
/home/username/lk99/01scf
% pw.x < input.in > out.o

複数コアが使えるのであれば、下記のようにすることで並列化出来ます。

% OMP_NUM_THREADS=1 mpirun -np 4 pw.x < input.in > out.o


NSCF計算

次はNSCF計算です。

次のように作業してください。

% pwd
/home/username/lk99
% mkdir 02nscf
% cd 02nscf
% cp ../01scf/input.in .
% cp -r ../01scf/out .

input.inを若干修正します。

% pwd
/home/username/lk99/02nscf
% diff input.in ../01scf/input.in 
3c3
<     calculation = 'nscf'
---
>     calculation = 'scf'
15c15,17
<     occupations = 'tetrahedra'
---
>     occupations = 'smearing'
>     degauss     = 0.003675
>     smearing    = 'gaussian'
71c73
< 6 6 7 0 0 0 
---
> 4 4 5 0 0 0

コピペできるように、全文も載せておきます。

▼input.inの例

% cat input.in 
&CONTROL
    prefix      = 'lk99'
    calculation = 'nscf'
    tstress     = .true.
    tprnfor     = .true.
    outdir      = './out'
    pseudo_dir  = '../pp'
/
&SYSTEM
    ibrav       = 0
    nat         = 41
    ntyp        = 4
    ecutwfc     = 40
    ecutrho     = 160
    occupations = 'tetrahedra'
/
&ELECTRONS
    electron_maxstep = 100
    conv_thr         = 1e-09
    mixing_beta      = 0.7
    diago_david_ndim = 4
/
ATOMIC_SPECIES
Pb         207.2 Pb.pbe-dn-kjpaw_psl.1.0.0.UPF
Cu        63.546 Cu.pbe-spn-kjpaw_psl.1.0.0.UPF
P   30.973761998 P.pbe-nl-kjpaw_psl.1.0.0.UPF
O         15.999 O.pbe-n-kjpaw_psl.1.0.0.UPF
ATOMIC_POSITIONS angstrom
Pb       1.1293232294     3.7547448097     1.8340405970
Pb      -1.3023964628     7.9497623531     5.5916375765
Pb       1.1489751864     7.7008734208     1.8340405970
Pb      -1.2681567294     3.4974336207     5.5916375765
Pb      -2.2782984159     5.7448282092     1.8340405970
Pb       2.5705531922     5.7532504658     5.5916375765
Pb       4.9653411911     2.8667410733     0.0757794710
Pb       4.9653411911     2.8667410733     3.6777530328
Pb       0.0000000000     0.0000000000     3.8664200002
Cu       0.0000000000     0.0000000000     0.4705653209
P        8.2105722025     2.2426423547     1.7272170358
P        6.7113430460     0.4992594225     5.5576908611
P        3.8832110300     5.9892429296     1.7272170358
P        6.1426395162     5.5625638601     5.5576908611
P        2.8022403406     0.3683379355     1.7272170358
P        2.0420410111     2.5383999373     5.5576908611
O        6.7139137896     2.6538520455     1.8563204123
O        3.1386080324     8.4202916340     5.5507011195
O        4.2754221980     4.4874938779     1.8563204123
O       -3.8961492874     7.1081916914     5.5507011195
O        3.9066875856     1.4588772965     1.8563204123
O        0.7575412550     1.6719631142     5.5507011195
O       -1.0142168212     3.1313758775     0.6501859633
O        0.8614076797     8.5737741449     6.7308651389
O        2.7605985430     6.1561977489     0.6501859633
O       -2.8904688646     5.0593370810     6.7308651389
O       -1.7463817217     7.9128728132     0.6501859633
O        2.0290611848     3.5673352137     6.7308651389
O        5.8943805130     0.2560738791     4.2519893691
O       -0.9612700690     2.4666418617     3.0685432816
O        6.7617256411     4.9766463243     4.2519893691
O        3.3098017112     6.5344179893     3.0685432816
O        2.2399174191     3.3675030165     4.2519893691
O       -2.3485316423     8.1993865886     3.0685432816
O        8.2878846617     0.7382457515     1.3075648471
O        6.8662625520     2.0641280221     5.7234632213
O        5.1474004762     6.8083957849     1.3075648471
O        4.7099638023     4.9142937880     5.7234632213
O        1.4607384353     1.0535816834     1.3075648471
O        3.3197972189     1.6218014097     5.7234632213
O        0.0000000000     5.7334821466     2.1622751229
K_POINTS {automatic}
6 6 7 0 0 0 
CELL_PARAMETERS angstrom
    9.9306823821     0.0000000000     0.0000000000
   -4.9653411911     8.6002232198     0.0000000000
    0.0000000000     0.0000000000     7.4109945173

次のコマンドで実行です。

% pwd
/home/username/lk99/02nscf
% pw.x < input.in > out.o

この計算が終わるとフェルミエネルギーが取得できるので、 控えておいてください。

% grep Fermi out.o 
     the Fermi energy is     6.6810 ev


バンド計算

最後に、お待ちかねのバンド計算です。

% pwd
/home/username/lk99
% mkdir 03bands
% cd 03bands
% cp ../01scf/input.in .
% cp -r ../01scf/out .

input.inは下記のように編集してください。

▼input.inの例

% cat input.in 
&CONTROL
    prefix      = 'lk99'
    calculation = 'bands'
    tstress     = .true.
    tprnfor     = .true.
    outdir      = './out'
    pseudo_dir  = '../pp'
    verbosity   = 'high'
/
&SYSTEM
    ibrav       = 0
    nat         = 41
    ntyp        = 4
    ecutwfc     = 40
    ecutrho     = 160
    occupations = 'smearing'
    degauss     = 0.003675
    smearing    = 'gaussian'
/
&ELECTRONS
    electron_maxstep = 100
    conv_thr         = 1e-09
    mixing_beta      = 0.7
    diago_david_ndim = 4
/
ATOMIC_SPECIES
Pb         207.2 Pb.pbe-dn-kjpaw_psl.1.0.0.UPF
Cu        63.546 Cu.pbe-spn-kjpaw_psl.1.0.0.UPF
P   30.973761998 P.pbe-nl-kjpaw_psl.1.0.0.UPF
O         15.999 O.pbe-n-kjpaw_psl.1.0.0.UPF
ATOMIC_POSITIONS angstrom
Pb       1.1293232294     3.7547448097     1.8340405970
Pb      -1.3023964628     7.9497623531     5.5916375765
Pb       1.1489751864     7.7008734208     1.8340405970
Pb      -1.2681567294     3.4974336207     5.5916375765
Pb      -2.2782984159     5.7448282092     1.8340405970
Pb       2.5705531922     5.7532504658     5.5916375765
Pb       4.9653411911     2.8667410733     0.0757794710
Pb       4.9653411911     2.8667410733     3.6777530328
Pb       0.0000000000     0.0000000000     3.8664200002
Cu       0.0000000000     0.0000000000     0.4705653209
P        8.2105722025     2.2426423547     1.7272170358
P        6.7113430460     0.4992594225     5.5576908611
P        3.8832110300     5.9892429296     1.7272170358
P        6.1426395162     5.5625638601     5.5576908611
P        2.8022403406     0.3683379355     1.7272170358
P        2.0420410111     2.5383999373     5.5576908611
O        6.7139137896     2.6538520455     1.8563204123
O        3.1386080324     8.4202916340     5.5507011195
O        4.2754221980     4.4874938779     1.8563204123
O       -3.8961492874     7.1081916914     5.5507011195
O        3.9066875856     1.4588772965     1.8563204123
O        0.7575412550     1.6719631142     5.5507011195
O       -1.0142168212     3.1313758775     0.6501859633
O        0.8614076797     8.5737741449     6.7308651389
O        2.7605985430     6.1561977489     0.6501859633
O       -2.8904688646     5.0593370810     6.7308651389
O       -1.7463817217     7.9128728132     0.6501859633
O        2.0290611848     3.5673352137     6.7308651389
O        5.8943805130     0.2560738791     4.2519893691
O       -0.9612700690     2.4666418617     3.0685432816
O        6.7617256411     4.9766463243     4.2519893691
O        3.3098017112     6.5344179893     3.0685432816
O        2.2399174191     3.3675030165     4.2519893691
O       -2.3485316423     8.1993865886     3.0685432816
O        8.2878846617     0.7382457515     1.3075648471
O        6.8662625520     2.0641280221     5.7234632213
O        5.1474004762     6.8083957849     1.3075648471
O        4.7099638023     4.9142937880     5.7234632213
O        1.4607384353     1.0535816834     1.3075648471
O        3.3197972189     1.6218014097     5.7234632213
O        0.0000000000     5.7334821466     2.1622751229
K_POINTS crystal
73
    0.0000000000     0.0000000000     0.0000000000 1 ! GAMMA 
    0.0769230769     0.0000000000     0.0000000000 1
    0.1538461538     0.0000000000     0.0000000000 1
    0.2307692308     0.0000000000     0.0000000000 1
    0.3076923077     0.0000000000     0.0000000000 1
    0.3846153846     0.0000000000     0.0000000000 1
    0.4615384615     0.0000000000     0.0000000000 1
    0.5000000000     0.0000000000     0.0000000000 1 ! M
    0.4523809524     0.0952380952     0.0000000000 1
    0.4047619048     0.1904761905     0.0000000000 1
    0.3571428571     0.2857142857     0.0000000000 1
    0.3333333333     0.3333333333     0.0000000000 1 ! K
    0.2888888889     0.2888888889     0.0000000000 1
    0.2444444444     0.2444444444     0.0000000000 1
    0.2000000000     0.2000000000     0.0000000000 1
    0.1555555556     0.1555555556     0.0000000000 1
    0.1111111111     0.1111111111     0.0000000000 1
    0.0666666667     0.0666666667     0.0000000000 1
    0.0222222222     0.0222222222     0.0000000000 1
    0.0000000000     0.0000000000     0.0000000000 1 ! G
    0.0000000000     0.0000000000     0.0666666667 1
    0.0000000000     0.0000000000     0.1333333333 1
    0.0000000000     0.0000000000     0.2000000000 1
    0.0000000000     0.0000000000     0.2666666667 1
    0.0000000000     0.0000000000     0.3333333333 1
    0.0000000000     0.0000000000     0.4000000000 1
    0.0000000000     0.0000000000     0.4666666667 1
    0.0000000000     0.0000000000     0.5000000000 1 ! A
    0.0769230769     0.0000000000     0.5000000000 1
    0.1538461538     0.0000000000     0.5000000000 1
    0.2307692308     0.0000000000     0.5000000000 1
    0.3076923077     0.0000000000     0.5000000000 1
    0.3846153846     0.0000000000     0.5000000000 1
    0.4615384615     0.0000000000     0.5000000000 1
    0.5000000000     0.0000000000     0.5000000000 1 ! L
    0.4523809524     0.0952380952     0.5000000000 1
    0.4047619048     0.1904761905     0.5000000000 1
    0.3571428571     0.2857142857     0.5000000000 1
    0.3333333333     0.3333333333     0.5000000000 1 ! H
    0.2888888889     0.2888888889     0.5000000000 1
    0.2444444444     0.2444444444     0.5000000000 1
    0.2000000000     0.2000000000     0.5000000000 1
    0.1555555556     0.1555555556     0.5000000000 1
    0.1111111111     0.1111111111     0.5000000000 1
    0.0666666667     0.0666666667     0.5000000000 1
    0.0222222222     0.0222222222     0.5000000000 1
    0.0000000000     0.0000000000     0.5000000000 1 ! A
    0.5000000000     0.0000000000     0.5000000000 1 ! L
    0.5000000000     0.0000000000     0.4333333333 1
    0.5000000000     0.0000000000     0.3666666667 1
    0.5000000000     0.0000000000     0.3000000000 1
    0.5000000000     0.0000000000     0.2333333333 1
    0.5000000000     0.0000000000     0.1666666667 1
    0.5000000000     0.0000000000     0.1000000000 1
    0.5000000000     0.0000000000     0.0333333333 1
    0.5000000000     0.0000000000     0.0000000000 1 ! M
    0.3333333333     0.3333333333     0.5000000000 1 ! H
    0.3333333333     0.3333333333     0.4333333333 1
    0.3333333333     0.3333333333     0.3666666667 1
    0.3333333333     0.3333333333     0.3000000000 1
    0.3333333333     0.3333333333     0.2333333333 1
    0.3333333333     0.3333333333     0.1666666667 1
    0.3333333333     0.3333333333     0.1000000000 1
    0.3333333333     0.3333333333     0.0333333333 1
    0.3333333333     0.3333333333     0.0000000000 1 ! K
    0.3333333333     0.3333333333    -0.0666666667 1
    0.3333333333     0.3333333333    -0.1333333333 1
    0.3333333333     0.3333333333    -0.2000000000 1
    0.3333333333     0.3333333333    -0.2666666667 1
    0.3333333333     0.3333333333    -0.3333333333 1
    0.3333333333     0.3333333333    -0.4000000000 1
    0.3333333333     0.3333333333    -0.4666666667 1
    0.3333333333     0.3333333333    -0.5000000000 1 ! H2 
CELL_PARAMETERS angstrom
    9.9306823821     0.0000000000     0.0000000000
   -4.9653411911     8.6002232198     0.0000000000
    0.0000000000     0.0000000000     7.4109945173

この"K_POINTS crystal"の下に書いてあるのが、 まさにバンド図の描き方を指定しているところです。

先程アクセスしたSeeK-pathの出力をそのまま使うと、 このバンド図の指定が細かすぎて、計算に時間が かかってしまうため、私の方で不要な部分はカットし、 また、てきとうに間引いてあります。

次のコマンドで実行です。

% pwd
/home/username/lk99/03bands
% pw.x < input.in > out.o


バンド図のプロット

計算したバンドをプロットします。 下記の作業を行ってください。

% pwd
/home/username/lk99
% mkdir 04plotbands
% cd 04plotbands
% cp -r ../03bands/out .

次のようなinput.inを作成してください。

% cat input.in 
&bands
  outdir  = './out'
  prefix  = 'lk99'
  filband = 'lk99.band'
  lsym    = .true.
/

次のコマンドで実行です。

% pwd
/home/username/lk99/04plotbands
% bands.x < input.in > out.o

計算が終了したら、結果をプロットします。

gnuplotなどを使えば、このままでもプロットできるのですが、 あまり綺麗な図にならないので、今回はpythonでプロットしたいと思います。

下記のようなpythonスクリプトを用意してください。

▼pythonコードの例

import numpy as np
import matplotlib.pyplot as plt

def read_bands(filename):
    bands = []
    band  = []
    for line in open(filename):
        line = line.split()
        if len(line)==0:
            bands.append(np.array(band))
            band = []
        else:
            band.append([ float(v) for v in line ])
    return(bands)
            
def get_p_online(b1, b2, b3, path):
    dk = 0.0
    positions = [dk]
    for i in range(1, len(path)):
        ks  = path[i-1][0]*b1 + path[i-1][1]*b2 + path[i-1][2]*b3
        kd  = path[i][0]*b1 + path[i][1]*b2 + path[i][2]*b3
        v   = (kd - ks)
        dk += np.sqrt(sum(v**2))
        positions.append(dk)
    return(positions)

def main():
    plt.rcParams["font.size"] = 18
    plt.figure(figsize=(16, 9))
    
    elec_num = (14*9 + 19*1 + 5*6 + 6*25 - 1)//2 + 1
    fermi    = 6.6810
    lat_a    = 9.9306823821
    bandfile = "./lk99.band.gnu"
    
    b1 = np.array([0.6327042861, 0.3652919899, 0.0000000000])
    b2 = np.array([0.0000000000, 0.7305839798, 0.0000000000])
    b3 = np.array([0.0000000000, 0.0000000000, 0.8478194516])
    
    # (G, M, K, G, A, L, H, A), (L, M), (H, K, H2)
    vG  =  np.array([0.0000000000, 0.0000000000,  0.0000000000]) 
    vM  =  np.array([0.5000000000, 0.0000000000,  0.0000000000]) 
    vK  =  np.array([0.3333333333, 0.3333333333,  0.0000000000]) 
    vA  =  np.array([0.0000000000, 0.0000000000,  0.5000000000]) 
    vL  =  np.array([0.5000000000, 0.0000000000,  0.5000000000]) 
    vH  =  np.array([0.3333333333, 0.3333333333,  0.5000000000]) 
    vH2 =  np.array([0.3333333333, 0.3333333333, -0.5000000000]) 

    s1  = get_p_online(b1, b2, b3, (vG, vM, vK, vG, vA, vL, vH, vA))
    s2  = get_p_online(b1, b2, b3, (vL, vM))
    s3  = get_p_online(b1, b2, b3, (vH, vK, vH2))

    s1  = [ s/((2.0*np.pi)/lat_a) for s in s1 ]
    s2  = [ s1[-1] + s/((2.0*np.pi)/lat_a) for s in s2 ]
    s3  = [ s2[-1] + s/((2.0*np.pi)/lat_a) for s in s3 ]

    plt.axhline(y=0, color="gray", linestyle="dashed")
    for k1 in s1:
        plt.axvline(x=k1, color="gray", linestyle="dashed")
    for k2 in s2:
        plt.axvline(x=k2, color="gray", linestyle="dashed")
    for k3 in s3:
        plt.axvline(x=k3, color="gray", linestyle="dashed")

    plt.axvline(x=s1[0],  lw=3, color="dimgray", linestyle="solid")
    plt.axvline(x=s1[-1], lw=3, color="dimgray", linestyle="solid")
    plt.axvline(x=s2[-1], lw=3, color="dimgray", linestyle="solid")
    plt.axvline(x=s3[-1], lw=3, color="dimgray", linestyle="solid")

    dl = 0.08
    xticks = s1[:-1] + [s1[-1]-dl] + [s2[0]+dl] + s2[1:-1] + [s2[-1]-dl] + [s3[0]+dl] + s3[1:]
    plt.xticks(xticks, (r"$\Gamma$", "M", "K", r"$\Gamma$", "A", "L", "H", "A", "L", "M", "H", "K", r"H$_2$"))
        
    bands = read_bands(bandfile)
    for bi, band in enumerate(bands):
        if bi==elec_num-1:
            plt.plot(band[:, 0], band[:, 1]-fermi, color="red")
        else:
            plt.plot(band[:, 0], band[:, 1]-fermi, color="black")

    plt.xlim([s1[0], s3[-1]])
    plt.ylim([-2, 4])
    #plt.ylim([-0.1, 0.1])

    plt.ylabel("Energy (eV)")
    
    plt.show()
    #plt.savefig("./bands_wide.png")

if __name__=="__main__":
    main()

これは、今回のために私がてきとうに作ったものなので、 あまり汎用性がないことに注意してください。

plot.pyなどの名前をつけて保存し、実行してください。

% pwd
/home/username/lk99/04plotbands
% python plot.py

つぎのようなバンド図が得られます。

QEで計算したLK-99のバンド図

フェルミ面近傍のバンドを赤で表示しています。参考にした論文同様、平坦なバンドが得られていることが確認できますね。


ほいよ
これ、フリーのソフトで作ったLK-99のフラットバンド図ね
追加でHubbard U足してもいいし、VASPに30万出すことはないんだよ



あとがき

現在はまだ、世界各国のグループが検証を進めている段階ですので、 このLK-99が本当に常温常圧超伝導なのかは定かではありません。 この辺の最新の情報については、日本語版Wikipediaでも活発に 情報が更新されているので、気になる方は随時チェックすると良いかも知れません。

ja.wikipedia.org

こうした世間の盛り上がりとは裏腹に、専門家の方々からはどこか落ち着いている印象を受けます。 それもそのはずで、この超伝導に関しては、過去多くの科学者が振り回されてきた歴史があるのです。 有名なものとしては、当時ベル研にいたヘンドリック・シェーンの起こした、通称シェーン事件というものがあります。 こちらについては、下記の動画に大変わかり易く、そして面白くまとまっているので、是非見てみてください。

www.nicovideo.jp

LK-99の現状については、下記のサイトに非常に分かり易くまとまっていました。

science-log.com



[参考文献]

・ チャールズ キッテル (著)「キッテル 固体物理学入門 第8版〈上〉」
・ 小口 多美夫 (著)「バンド理論―物質科学の基礎として (材料学シリーズ)」
・ 前園 涼 (著), 市場 友宏 (著)「動かして理解する 第一原理電子状態計算:DFTパッケージによるチュートリアル」