Google Colaboratoryを使って無料で第一原理計算

最近、新型コロナで実験ができない人のために「物性実験家のための無料でできる第一原理計算入門」
https://cometscome.github.io/DFT/build/
という記事を書いているのですが、その過程でGoogle Colabだけで第一原理計算をする方法がわかりましたので、こちらでもシェアしたいと思います。

## 必要なもの
- ブラウザ
- Googleのアカウント

## 使うソフト
- Quantum Espresso:第一原理計算ソフト
- Atomic Simulation Environment (ASE):原子分子シミュレーションのためのPythonライブラリ
https://qiita.com/cometscome_phys/items/4af134de6d959a7718b9



## Quantum Espressoのインストール
Google ColabでQuantum Espressoを使うためには、Qunatum Espressoをソースからコンパイルする必要がありますが、以下のような方法で可能であることがわかりました。
ポイントは、FFTW3をインストールしておくことです。

```
!git clone https://github.com/QEF/q-e.git
!apt-get install -y libfftw3-3 libfftw3-dev libfftw3-doc
%cd q-e
!DFLAGS='-D__OPENMP -D__FFTW3 -D__MPI -D__SCALAPACK' FFT_LIBS='-lfftw3'  ./configure --enable-openmp
```
ソースをダウンロードし、FFTW3をインストール、そしてコンパイル準備をしています。
次に、

```
!make pw
```
でQuantum Espressoのコアのコードであるpw.xを実行できます。
あとはポストプロセスツールである、pp:

```
make pp
```
もインストールしておきましょう。

さて、[Google colaboratoryでOSSをビルドする](https://qiita.com/kuronekodaisuki/items/d7910259120ef805f686)
にありますように、ここでコンパイルしたバイナリは12時間後に消滅してしまいます。
ですので、

```
from google.colab import drive
drive.mount('/content/drive')
```
を実行してご自分のGoogle Driveにバイナリを保存しておきましょう。次に使うときはこれを解凍して使います。

```
%cd /content/
!zip -r /content/drive/'My Drive'/q-e.zip q-e 
```
としてq-e.zipを保存しておきます。

Quantum Espressoを実行できるように、環境変数を設定します。つまり、

```
import os
os.environ['PATH'] = "/content/q-e/bin:"+os.environ['PATH']
```
とします。

## ASEのインストール
ASEのインストールは簡単で、

```
!pip install ase
```
で入ります。

## 第一原理計算のテスト
では、実際に第一原理計算をやってみましょう。

### NaClの構造最適化
NaClのディレクトリを作成し、移動します。

```
%cd /content
!mkdir NaCl
%cd NaCl
```
そして、擬ポテンシャルをダウンロードします。

```
!wget https://www.quantum-espresso.org/upf_files/Na.pbesol-spn-kjpaw_psl.1.0.0.UPF
!wget https://www.quantum-espresso.org/upf_files/Cl.pbesol-n-kjpaw_psl.1.0.0.UPF
```
これはNaClディレクトリに入りました。

次に、

```python
from ase.build import bulk
from ase.calculators.espresso import Espresso
from ase.constraints import UnitCellFilter
from ase.optimize import LBFGS
import ase.io 

pseudopotentials = {'Na': 'Na.pbesol-spn-kjpaw_psl.1.0.0.UPF',
                    'Cl': 'Cl.pbesol-n-kjpaw_psl.1.0.0.UPF'}  
rocksalt = bulk('NaCl', crystalstructure='rocksalt', a=6.0)
calc = Espresso(pseudopotentials=pseudopotentials,pseudo_dir = './',
                tstress=True, tprnfor=True, kpts=(3, 3, 3))

rocksalt.set_calculator(calc)

ucf = UnitCellFilter(rocksalt)
opt = LBFGS(ucf)
opt.run(fmax=0.005)

# cubic lattic constant
print((8*rocksalt.get_volume()/len(rocksalt))**(1.0/3.0))
```
を実行すれば、NaClの構造最適化ができます。このコードでは擬ポテンシャルの場所を```pseudo_dir```で指定しています。今回は今いるディレクトリですね。

### Cuのバンド図
次に、Cuのバンド図を計算してみます。

Cuのディレクトリを作成します。

```
%cd /content
!mkdir Cu
%cd Cu
```

擬ポテンシャルをダウンロードします。

```
!wget https://www.quantum-espresso.org/upf_files/Cu.pz-d-rrkjus.UPF
```

自己無撞着計算をして、電子密度を決定します。電子密度が決定されれば、各k点での計算をすることでバンド図を 描くことができます。

```python
from ase import Atoms
from ase.build import bulk
from ase.calculators.espresso import Espresso
atoms = bulk("Cu")
pseudopotentials = {'Cu':'Cu.pz-d-rrkjus.UPF'}

input_data = {
    'system': {
        'ecutwfc': 30,
        'ecutrho': 240,
        'nbnd' : 35,
    'occupations' : 'smearing',
        'smearing':'gauss',
        'degauss' : 0.01},
    'disk_io': 'low'}  # automatically put into 'control'

calc = Espresso(pseudopotentials=pseudopotentials,kpts=(4, 4, 4),input_data=input_data,pseudo_dir = './')
atoms.set_calculator(calc)

atoms.get_potential_energy()
fermi_level = calc.get_fermi_level()
print(fermi_level)
```

そして、バンド図のための計算をします。

```python
input_data.update({'calculation':'bands',
                              'restart_mode':'restart',
                               'verbosity':'high'})
calc.set(kpts={'path':'GXWLGK', 'npoints':100},
          input_data=input_data)
calc.calculate(atoms)
```
ここで、kptsのpathに好きなブリルアンゾーンの点を入れることで、簡単にバンド図を描くことができるのがASEの面白い点です。

最後にバンド図を計算します。

```python
import matplotlib.pyplot as plt

bs = calc.band_structure()
bs.reference = fermi_level
bs.plot(emax=40,emin=5)
```

<img width="411" alt="Cu.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/246113/8c9ad68d-a9d1-e466-dde6-0095d5e3f9b8.png">


これで、ブラウザだけで第一原理計算を実行できるようになりました。

## 計算の再開
別に日に計算を再開したい場合は、以下のようにやるとできると思います。

```python
from google.colab import drive
drive.mount('/content/drive')
!cp /content/drive/'My Drive'/q-e.zip ./
!unzip q-e.zip
import os
os.environ['PATH'] = "/content/q-e/bin:"+os.environ['PATH']
!pip install ase
```
あとは、ディレクトリを作って、そこに擬ポテンシャルをダウンロードしたりすればよいでしょう。

```python
!mkdir NaCl
%cd NaCl
!wget https://www.quantum-espresso.org/upf_files/Na.pbe-spn-kjpaw_psl.1.0.0.UPF
!wget https://www.quantum-espresso.org/upf_files/Cl.pbe-n-rrkjus_psl.1.0.0.UPF
```