Quantum Espressoによる構造最適化

https://qiita.com/hiro949/items/2ea86ef94df2f5d7ab41

の補足資料です。

第一原理計算では実際の原子が作るポテンシャルとは似せてはいますが異なる擬ポテンシャルを使用しているので、擬ポテンシャルのもとで形成されるべき結晶構造が実際とは異なる可能性があります。そのため、実際の結晶構造で計算すると、エネルギー的に不安定だったり圧力やひずみがかかった状態になっていたりこともあります。
方針の１つとしてエネルギー的に不安定でも実験的にはその構造が正しいのだからそのまま計算するという方針(こちらの方が多いかも)と、擬ポテンシャルの下で安定な構造に微調整してから計算するという方針があります。構造最適化は後者の方針を取りたい場合に安定な構造を探索するためのツールです。

また、圧力やひずみを加えた結晶の計算を行う場合にもこのプロセスが必要となります。

# 入力ファイル

基本的には電荷密度分布計算の入力ファイルと同じですが、いくつか異なる部分があります。

```Bi2Se3.ecut80.relax.in
 &CONTROL
                 calculation = 'vc-relax' ,
                  wf_collect = .true. ,
                      outdir = './output' ,
                  pseudo_dir = '~/q-e-qe-6.5/pseudo/' ,
                      prefix = 'Bi2Se3' ,
                     disk_io = 'low' ,
                   verbosity = 'high' ,
                     tstress = .true. ,
                     tprnfor = .true. ,
  max_seconds = 86000,
  restart_mode = from_scratch,
 /
 &SYSTEM
  ibrav = 0,
  celldm(1) = 7.81213,
  nat = 5,
  ntyp = 2,
                       ecutwfc = 80 ,
                     ecutrho = 800 ,
                 occupations = 'smearing' ,
                     degauss = 0.00001 ,
		     noncolin = .true. ,
		     lspinorb = .true. ,
  /
 &ELECTRONS
            electron_maxstep = 10000000,
                    conv_thr = 0.0000000001 ,
                 mixing_beta = 0.5 ,
/
&IONS
/
&CELL
            press = 0,
/
CELL_PARAMETERS {alat}
  0.577350269189626   0.000000000000000   2.306079664570230 
 -0.288675134594813   0.500000000000000   2.306079664570230 
 -0.288675134594813  -0.500000000000000   2.306079664570230 
ATOMIC_SPECIES
  Bi  208.98000  Bi.rel-pbe-dn-kjpaw_psl.1.0.0.UPF
  Se   78.96000  Se.pbe-dn-kjpaw_psl.1.0.0.UPF
  
ATOMIC_POSITIONS {crystal}
Bi   0.400460000000000   0.400460000000000   0.400460000000000 0 0 0
Bi   0.599540000000000   0.599540000000000   0.599540000000000 0 0 0
Se   0.209700000000000   0.209700000000000   0.209700000000000 0 0 0
Se   0.790300000000000   0.790300000000000   0.790300000000000 0 0 0
Se   0.000000000000000   0.000000000000000   0.000000000000000 0 0 0

K_POINTS automatic 
  4 4 4   0 0 0 
```

### パラメーター

よくいじるパラメーターとしては

- `&CELL`の`press`：結晶にかかる圧力
- `ATOMIC_POSITIONS` の　内部座標のあとの`0 0 0`：内部座標の最適化を行わないことを指定するモードです。このときは基本格子ベクトル$\vec{a}_1,\vec{a}_2,\vec{a}_3$が最適化パラメーターになります。これを外すと内部座標も最適化されます。ただし、外した場合はパラメーター数が多くなって計算がさらに重くなります。結晶の対称性はちゃんと保つように計算されます(https://www.mail-archive.com/users@lists.quantum-espresso.org/msg34576.html)。


# 実行コマンド

ファイルの構成は

https://qiita.com/hiro949/items/2ea86ef94df2f5d7ab41

と同じだとします。MPI並列で60並列して計算するとします。

```
mpirun -np 60 ~/q-e-qe-6.5/biin/pw.x <Bi2Se3.ecut80.scf.in >Bi2Se3.ecut80.scf.out
```

# 出力ファイル

最適化したときのエンタルピー、結晶にかかるひずみ(圧力)、新しい基本並進ベクトル、内部座標(今回の場合は最適化する前と同じ)などが出力されます。

# カットオフエネルギーの決定

先の入力では

```
          ecutwfc = 80 ,
          ecutrho = 800 ,
```

としましたが、`ecutwfc = 80, 100, 150, 200, 250, 300`と変えていった入力ファイルを作り、カットオフ依存性を見ていきます。(`ecutrho`は`ecutwfc`の10倍とします。)
これで新しい構造の格子長$a$のカットオフ$E_{cut}$に対する依存性を調べると、ある$E_{cut}^*$が存在して、$E_{cut}>E_{cut}^*$では格子長がほとんど変わらなくなります。
そのような$E_{cut}^*$を大体見積もることで最適なカットオフエネルギーを見積もることができます。

## コメント
あくまでこれは決め方の一つの流派みたいなものです。
必ずしもこのようにする必要はなく、擬ポテンシャルのソースファイルのコメントにある推奨カットオフエネルギーをそのまま使う方がよいという方々もいます。

## ジョブシーケンス

次のコードで`ecutwfc`を変えていった入力を順次実行させられます。(ファイル名も`ecut80`の部分を変えてください。)

```sub.sh
### input parameters
dir_exe='../..'
prefix='Bi2Se3'
ecutArr=(80 100 150 200 250 300)

Exe_mpi(){
    #var1: executable file
    #var2: input file
    #var3: output file
    mpirun -np 60 $1 <$2 >$3
}

cd ${PBS_O_WORKDIR}
for cutoff in ${ecutArr[@]}; do
    Sub_qe_mpi "$dir_exe/bin/pw.x" "input/$prefix.ecut$cutoff.relax.in" "log/$prefix.ecut$cutoff.relax.out"
done
```