Atomic Simulation Environment (ASE)とは

第一原理計算/量子化学計算を行うpythonのライブラリ
- https://wiki.fysik.dtu.dk/ase/
第一原理計算の主要部分はcalculatorとして条件を設定し、実際の計算のインプットを編集することなく計算を実行できる
さまざまな計算パッケージをcalculatorに用いることができる

環境設定

ASEのインストール: pip install ase
ASEに組み込まれているポテンシャル(EMT)を用いる場合はこれ以外の設定は不要

VASPを用いる場合

第一原理計算としてメジャーなソフトウェアであるVASPを用いる場合は設定が必要
VASPがインストールされていることが前提(実行ファイルは~/vasp.X.X.X/bin/vasp_stdとする)
ASEから参照するVASP実行用のPythonスクリプト(run_vasp.pyとする)を用意する

import os
exitcode = os.system('~/vasp.X.X.X/bin/vasp_std')
# exitcode = os.system('mpirun -np 2 ~/vasp.X.X.X/bin/vasp_std') # MPI並列の場合はこちらを使う

~/.bashrcに上記のVASP_SCRIPTの場所と擬ポテンシャルの場所(VASP_PP_PATH)を設定する

$ export VASP_SCRIPT=$HOME/somewhere/run_vasp.py
$ export VASP_PP_PATH=$HOME/vasp/potentials/

run_vasp.pyや擬ポテンシャルの場所は適宜設定すること

エネルギー計算

ASEを用いた計算の最も基本的な使い方としてエネルギー計算を行う
Pt結晶のself-consistent field (SCF)計算を行う
ここでは以下のcalculatorを用いる。それぞれ異なる計算手法となっており、計算コスト・計算精度が違ってくる

EMT (多体ポテンシャル)
DFT (VASP)
M3GNet (ニューラルネットワーク・ポテンシャル)
CHGNet (ニューラルネットワーク・ポテンシャル)

1. EMT

Effective medium theory (EMT): 多体ポテンシャルの一種
対応している元素: Al, Cu, Ag, Au, Ni, Pd, Pt, H, C, N, O

from ase.build import bulk
from ase.calculators.emt import EMT

solid = bulk(name="Pt", crystalstructure="fcc", a=3.9)

solid.calc = EMT()

energy = solid.get_potential_energy()
print(f"energy = {energy:5.3} eV")

2. DFT

Vienna ab initio simulation package (VASP)によるDFT計算からエネルギーを求める

from ase.build import bulk
from ase.calculators.vasp import Vasp

solid = bulk(name="Pt", crystalstructure="fcc", a=3.9)

solid.calc = Vasp(xc="pbe", pp="potpaw_PBE.64")

energy = solid.get_potential_energy()
print(f"energy = {energy:5.3} eV")

Vaspのxcは汎函数, ppは擬ポテンシャルの種類を指定する
RuntimeError: Looking for PP for potpaw_PBE/Pt/POTCARが出力される場合はVASP_PP_PATHもしくはppの指定がうまくいっていない
上記だと$VASP_PP_PATH/potpaw_PBE.64以下が参照される

3. M3GNet

ニューラルネットワーク・ポテンシャルの1種
DFTでフィットされたポテンシャルで、DFTと同等の精度であるが計算が早い
M3GNetのインストール: pip install matgl

from ase.build import bulk
import matgl
from matgl.ext.ase import PESCalculator

solid = bulk(name="Pt", crystalstructure="fcc", a=3.9)

potential = matgl.load_model("M3GNet-MP-2021.2.8-PES")
solid.calc = PESCalculator(potential=potential)

energy = solid.get_potential_energy()
print(f"energy = {energy:5.3} eV")

4. CHGNet

M3GNet同様でニューラルネットワーク・ポテンシャルの1種
CHGNetのインストール: pip install chgnet
M3GNetより精度が高いきがすr

from ase.build import bulk
from chgnet.model.dynamics import CHGNetCalculator
from chgnet.model.model import CHGNet
    
solid = bulk(name="Pt", crystalstructure="fcc", a=3.9)

chgnet = CHGNet.load()
potential = CHGNetCalculator(potential=chgnet, properties="energy")
solid.calc = potential
        
energy = solid.get_potential_energy()
print(f"energy = {energy:5.3} eV")

構造最適化

Ptのfcc(111)面を構築して構造最適化計算を行う
原子層は4層とし、上2層のみを最適化(下2層は固定)

1. EMT, M3GNet, CHGNet

from ase.calculators.emt import EMT
from ase.build import fcc111
from ase.optimize import BFGS
from ase.constraints import FixAtoms

# M3Gnet
import matgl
from matgl.ext.ase import PESCalculator

# CHGNet
from chgnet.model.dynamics import CHGNetCalculator
from chgnet.model.model import CHGNet
from pymatgen.io.ase import AseAtomsAdaptor

surf = fcc111(symbol="Pt", size=[1, 1, 4], a=3.9, vacuum=12.0)
surf.pbc = True
c = FixAtoms(indices=[atom.index for atom in surf if atom.tag >= 3])
surf.constraints = c

# --- EMT
surf.calc = EMT()

# --- M3GNet
# potential = matgl.load_model("M3GNet-MP-2021.2.8-PES")
# surf.calc = PESCalculator(potential=potential)

# --- CHGNet
# chgnet = CHGNet.load()
# potential = CHGNetCalculator(potential=chgnet, properties="energy")
# surf.calc = potential
        
opt = BFGS(surf, trajectory="pt-relax.traj")
opt.run(fmax=0.05, steps=100)

fcc111はlayer numberに応じてAtomのtagを自動的に決める
- 第一層: tag=1, 第二層: tag=2, 吸着子: tag=0, などとなる。このtagを用いて原子層を選択する。tagを確認するにはprint(atoms.get_tags())
.trajファイルの可視化: ase gui pt-relax.traj

2. VASP

from ase.calculators.vasp import Vasp
from ase.build import fcc111
from ase.constraints import FixAtoms

surf = fcc111(symbol="Pt", size=[1, 1, 4], a=3.9, vacuum=12.0)
surf.pbc = True
c = FixAtoms(indices=[atom.index for atom in surf if atom.tag >= 3])
surf.constraints = c

surf.calc = Vasp(xc="pbe", pp="potpaw_PBE.64", ibrion=2, nsw=10,potim=0.1, directory="pt_relax")
surf.get_potential_energy()

EMTのようにASEのoptimizerを使うのでも良いが、VASPの最適化を用いた方が収束が早い(と思う)
Vaspのibrion, nsw, potimはVASPのINCARで指定するものと同じ
directoryは必須ではないが、設定した場合そのディレクトリで(ない場合は作成される)VASPの計算が行われるので設定しておくとよい。以下のどちらかで設定できる
- Atoms.calc.directory = "some_directory"
- Vasp(direcoty="some_directory")

おわりに

ここではASEの基本的な使い方としてエネルギー計算と構造最適化の例を示した
CalculatorにはEMTとVASPを用いた。EMTはASEに含まれておりそのまま利用でき計算も早いが、計算精度が相当悪いのと、利用できる元素の種類が非常に少ない(H, C, N, O, Al, Cu, Ag, Au, Ni, Pd, Ptのみ)が欠点であり、コードの確認だけに用いた方がよい
他の機能は別ページで説明する予定