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

# はじめに  
物性研スパコンのohtakaでCrySPYとQuantum ESPRESSOを使って結晶構造探索を行う．
動作確認程度のチュートリアル．詳細は[ドキュメント](https://tomoki-yamashita.github.io/CrySPY_doc/index.html)を参照してほしい．

- 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からインストール
``` bash
pip3 install csp-cryspy
```
必要なライブラリや実行スクリプトがインストールされる．下記コマンドで実行スクリプトの場所を確認可能
``` bash
$ which cryspy
~/local/miniforge3.10/bin/cryspy
```

# Si8の構造探索
## 入力ファイルの準備

/work以下の自分のディレクトリに移動
``` bash
cd /work/kxxxx/kxxxxyy
```

[ここ](https://tomoki-yamashita.github.io/CrySPY_doc/cryspy_utility/examples/qe_si8_rs/index.html)にexampleがおいてあるのでダウンロードする

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

展開して中に移動
```bash
$ tar zxf qe_Si8_RS.tar.gz
$ cd qe_Si8_RS
$ ls
calc_in/  cryspy.in
```

cryspyのインプットcryspy.inをohataka向けに下記のように編集．
（`jobcmd = sbatch`に変えた）
``` txt
[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/に入っている．
```bash
$ 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`
```txt
 &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`
```txt
 &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`
```txt
#!/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直下に戻る
```bash
$ cd ..
$ ls
calc_in/  cryspy.in
```

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

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

```bash
$ ls
calc_in/  cryspy.in  cryspy.stat  data/  err_cryspy  log_cryspy
```
`log_cryspy`にログが出力される

```bash
$ 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にしているので）投げる．

```bash
cryspy &
```
ログチェック．

```bash
$ 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
```

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

下記はQEの計算が終わっている場合の例．
```bash
$ cat stat_job
0         # Structure ID
1         # Stage
done
```

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

```bash
cryspy &
```
```bash
$ 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の構造最適化が自動的に始まる．

```bash
cryspy &
```
```bash
$ 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の値を増やす
```txt
[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]
```
```bash
cryspy &
```
```bash
$ 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を増やす
```txt
[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]
```
```bash
cryspy &
```
tot_strucを増やした後の実行は構造生成モードに入る．この時QEの計算結果の回収やジョブのサブミットは行われない
```bash
$ 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以下にある