最近、新型コロナで実験ができない人のために「物性実験家のための無料でできる第一原理計算入門」
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をビルドする
にありますように、ここでコンパイルしたバイナリは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ディレクトリに入りました。
次に、
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点での計算をすることでバンド図を 描くことができます。
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)
そして、バンド図のための計算をします。
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の面白い点です。
最後にバンド図を計算します。
import matplotlib.pyplot as plt
bs = calc.band_structure()
bs.reference = fermi_level
bs.plot(emax=40,emin=5)
これで、ブラウザだけで第一原理計算を実行できるようになりました。
計算の再開
別に日に計算を再開したい場合は、以下のようにやるとできると思います。
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
あとは、ディレクトリを作って、そこに擬ポテンシャルをダウンロードしたりすればよいでしょう。
!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