物性研スパコンohtakaで結晶構造探索を行う

はじめに

物性研スパコンのohtakaでCrySPYとQuantum ESPRESSOを使って結晶構造探索を行う．
動作確認程度のチュートリアル．詳細はドキュメントを参照してほしい．

• Python: miniforge3を入れる
• CrySPY: version 1.0.0．PyPIからインストールする
• QE: ohtakaにデフォルトでインストールされているものを使う

準備

Python

ohtakaにはpython3のいくつかのバージョンがプリインストールされているが，どれもpymatgenが動かないので使わない．自分でminiforgeを入れる
https://github.com/conda-forge/miniforge
Miniforge3-Linux-x86_64は問題なくインストールできる．

CrySPY

PyPIからインストール

pip3 install csp-cryspy

必要なライブラリや実行スクリプトがインストールされる．下記コマンドで実行スクリプトの場所を確認可能

$ which cryspy
~/local/miniforge3.10/bin/cryspy

Si8の構造探索

入力ファイルの準備

/work以下の自分のディレクトリに移動

cd /work/kxxxx/kxxxxyy

ここにexampleがおいてあるのでダウンロードする

wget https://tomoki-yamashita.github.io/CrySPY_doc/examples/qe_Si8_RS.tar.gz

展開して中に移動

$ tar zxf qe_Si8_RS.tar.gz
$ cd qe_Si8_RS
$ ls
calc_in/  cryspy.in

cryspyのインプットcryspy.inをohataka向けに下記のように編集．
（jobcmd = sbatchに変えた）

[basic]
algo = RS
calc_code = QE
tot_struc = 5
nstage = 2
njob = 1
jobcmd = sbatch
jobfile = job_cryspy

[structure]
natot = 8
atype = Si
nat = 8

[QE]
qe_infile = pwscf.in
qe_outfile = pwscf.out
kppvol =  40  80

[option]

QE用のファイルがcalc_in/に入っている．

$ cd calc_in
$ ls
job_cryspy/  pwscf.in_1  pwscf.in_2

CrySPYではステージ分の入力ファイルをxxxx_1, xxxx_2の用に準備する．

• pwscf.in_1: ステージ1の入力ファイル（自由に設定できるがここではセルを固定した構造最適化）
• pwscf.in_2: ステージ2の入力ファイル（セルも含めたフル緩和）
• job_cryspy: ジョブスクリプト

QEのポテンシャルは今回はGBRV( https://www.physics.rutgers.edu/gbrv/ )のPBEを使用．ポテンシャルは自分で用意する．
ポテンシャルの置き場所やファイル名，カットオフなどを下記ように編集．
（自分のアカウントに合わせる）
CrySPYでは結晶構造とk点の情報は自動で追加されるので，書かなくて良い．

pwscf.in_1

 &control
    title = 'Si8'
    calculation = 'relax'
    nstep = 100
    restart_mode = 'from_scratch',
    pseudo_dir = '/home/kxxxx/kxxxxyy/local/gbrv/all_pbe_UPF_v1.5/'
    outdir='./out.d/'
 /

 &system
    ibrav = 0
    nat = 8
    ntyp = 1
    ecutwfc = 40
    ecutrho = 200
    occupations = 'smearing'
    degauss = 0.01
 /

 &electrons
 /

 &ions
 /

 &cell
 /

ATOMIC_SPECIES
  Si  28.086  si_pbe_v1.uspp.F.UPF

pwscf.in_2

 &control
    title = 'Si8'
    calculation = 'vc-relax'
    nstep = 200
    restart_mode = 'from_scratch',
    pseudo_dir = '/home/k0362/k036200/local/gbrv/all_pbe_UPF_v1.5/'
    outdir='./out.d/'
 /

 &system
    ibrav = 0
    nat = 8
    ntyp = 1
    ecutwfc = 40
    ecutrho = 200
    occupations = 'smearing'
    degauss = 0.01
 /

 &electrons
 /

 &ions
 /

 &cell
 /

ATOMIC_SPECIES
  Si  28.086  si_pbe_v1.uspp.F.UPF

QEのジョブスクリプトは物性研ウェブサイトの情報を参考にして編集する
https://mdcl.issp.u-tokyo.ac.jp/scc/system/quantum-espresso_guide

job_cryspy

#!/bin/sh

#SBATCH -p i8cpu
#SBATCH -N 1
#SBATCH -n 4
#SBATCH -c 32
#SBATCH -t 00:10:00

set -e

source /home/issp/materiapps/intel/espresso/espressovars.sh
module list


srun pw.x -nk 4  < pwscf.in > pwscf.out


# ---------- for CrySPY
if [ -e "CRASH" ]; then
    sed -i -e '3 s/^.*$/skip/' stat_job
    exit 1
fi

sed -i -e '3 s/^.*$/done/' stat_job

i8cpuはデバッグ用なので利用には注意

CrySPYではジョブスクリプトの最後に
sed -i -e '3 s/^.*$/done/' stat_jobと書くことが決まり事．

実行

qe_Si8_RS直下に戻る

$ cd ..
$ ls
calc_in/  cryspy.in

初回実行（pythonライブラリのインストール直後の初回起動はコンパイルだとかキャッシュとかで時間がかかる．CrySPYが遅い訳ではない）

cryspy &

CrySPYを実行すると最初は構造生成モードに入り，cryspy.inのtot_strucで指定した構造数のランダム構造が作られる．今回は5構造．

$ ls
calc_in/  cryspy.in  cryspy.stat  data/  err_cryspy  log_cryspy

log_cryspyにログが出力される

$ cat log_cryspy
2023/04/05 17:58:59
CrySPY 1.0.0
Start cryspy.py

Read input file, cryspy.in
Save input data in cryspy.stat

# --------- Generate initial structures
# ------ mindist
Si - Si 1.11
Structure ID      0 was generated. Space group:  96 -->  96 P4_32_12
Structure ID      1 was generated. Space group: 132 --> 132 P4_2/mcm
Structure ID      2 was generated. Space group: 225 --> 225 Fm-3m
Structure ID      3 was generated. Space group:  27 -->  27 Pcc2
Structure ID      4 was generated. Space group: 216 --> 216 F-43m

もう一度cryspyを実行する．CrySPYはcryspy.statファイルがあるかどうかで最初から開始するか途中から再開するか判断している．
構造生成が終わったので，次は構造最適化のジョブを一つ（njob = 1にしているので）投げる．

cryspy &

ログチェック．

$ cat log_cryspy
2023/04/05 17:58:59
CrySPY 1.0.0
Start cryspy.py

Read input file, cryspy.in
Save input data in cryspy.stat

# --------- Generate initial structures
# ------ mindist
Si - Si 1.11
Structure ID      0 was generated. Space group:  96 -->  96 P4_32_12
Structure ID      1 was generated. Space group: 132 --> 132 P4_2/mcm
Structure ID      2 was generated. Space group: 225 --> 225 Fm-3m
Structure ID      3 was generated. Space group:  27 -->  27 Pcc2
Structure ID      4 was generated. Space group: 216 --> 216 F-43m



2023/04/05 18:05:38
CrySPY 1.0.0
Restart cryspy.py



# ---------- job status
ID      0: submit job, Stage 1

計算がちゃんとできているか確認する．

cd work/000000

ID 0，stage 1の計算がここで行われている．job_cryspyやpwscf.inがコピーされて来ているのも確認できるはず．
stat_jobというファイルはCrySPYが作っており，ジョブが終わるとここの3行目がdoneになるようにしている．初期はsubmitted．CrySPYはこのファイルがdoneになっているかどうかでジョブが終わったかをチェックする．

下記はQEの計算が終わっている場合の例．

$ cat stat_job
0         # Structure ID
1         # Stage
done

もう一度cryspyを実行すると，stage 1が終わっているのでstage 2に移る．

cryspy &

$ tail log_cryspy
2023/04/05 18:12:27
CrySPY 1.0.0
Restart cryspy.py



# ---------- job status
ID      0: Stage 1 Done!
    submitted job, ID      0 Stage 2

stage 2の計算が終わった後にもう一度cryspyを実行すると，ID 0の最適化が完了するのでエネルギーと最適化後の構造データを収集して，次のIDの構造最適化が自動的に始まる．

cryspy &

$ tail log_cryspy
# ---------- job status
ID      0: Stage 2 Done!
    collect results: E = -130.04242681082383 eV/atom


 recheck 1


# ---------- job status
ID      1: submit job, Stage 1

あとは繰り返すだけでOK．構造最適化ジョブが一つずつじゃめんどくさいので同時にいっぱいサブミットしたい時はcryspy.inのnjobの値を増やす

[basic]
algo = RS
calc_code = QE
tot_struc = 5
nstage = 2
njob = 5
jobcmd = sbatch
jobfile = job_cryspy

[structure]
natot = 8
atype = Si
nat = 8

[QE]
qe_infile = pwscf.in
qe_outfile = pwscf.out
kppvol =  40  80

[option]

cryspy &

$ tail -n 15 log_cryspy                                                   


2023/04/05 18:18:02
CrySPY 1.0.0
Restart cryspy.py


Changed njob from 1 to 5

# ---------- job status
ID      1: Stage 1 Done!
    submitted job, ID      1 Stage 2
ID      2: submit job, Stage 1
ID      3: submit job, Stage 1
ID      4: submit job, Stage 1

構造数が5じゃ少なければ追加可能．cryspy.inのtot_strucを増やす

[basic]
algo = RS
calc_code = QE
tot_struc = 20
nstage = 2
njob = 5
jobcmd = sbatch
jobfile = job_cryspy

[structure]
natot = 8
atype = Si
nat = 8

[QE]
qe_infile = pwscf.in
qe_outfile = pwscf.out
kppvol =  40  80

[option]

cryspy &

tot_strucを増やした後の実行は構造生成モードに入る．この時QEの計算結果の回収やジョブのサブミットは行われない

$ tail -n 30 log_cryspy                                                                                                                                                            [/work/k0362/k036200/qe_Si8_RS]


2023/04/05 18:20:47
CrySPY 1.0.0
Restart cryspy.py


Changed tot_struc from 5 to 20

Backup data

# ---------- Append structures
# ------ mindist
Si - Si 1.11
Structure ID      5 was generated. Space group:   3 -->   3 P2
Structure ID      6 was generated. Space group:  54 -->  67 Cmme
Structure ID      7 was generated. Space group:  42 -->  42 Fmm2
Structure ID      8 was generated. Space group:  39 -->  39 Aem2
Structure ID      9 was generated. Space group: 126 --> 123 P4/mmm
Structure ID     10 was generated. Space group: 156 --> 156 P3m1
Structure ID     11 was generated. Space group:  33 -->  33 Pna2_1
Structure ID     12 was generated. Space group:  34 -->  34 Pnn2
Structure ID     13 was generated. Space group:   7 -->   7 Pc
Structure ID     14 was generated. Space group:  54 -->  54 Pcca
Structure ID     15 was generated. Space group: 163 --> 194 P6_3/mmc
Structure ID     16 was generated. Space group:  95 -->  95 P4_322
Structure ID     17 was generated. Space group:  61 -->  47 Pmmm
Structure ID     18 was generated. Space group:  21 -->  21 C222
Structure ID     19 was generated. Space group:  89 --> 123 P4/mmm

cryspyの実行を繰り返すと最適化が終わりデータが収集できる

結果

• data/以下にたまる
• 基本的にはcryspy_rslt_energy_ascを見る
• init_POSCARS（初期構造）やopt_POSCARS, opt_CIFS（最適化構造）はVESTAで開くことができて，構造を可視化できる
• pythonのpklデータがdata/pkl_data以下にある