Atomic Simulation Environmentの使い方: エネルギー計算と構造最適化

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

## ASEを使う利点
* 第一原理計算や分子動力学計算を行う計算プログラムは非常に数が多く、それぞれが独自のインプット・アプトプットを持つ
* これらのプログラムを用いて計算を行うためには、インプットの作成方法やアウトプットの解析方法に習熟していることが望ましく、学習コストが高い
* ASEを用いることで、`calculator`の設定を通して計算プログラムのインプット作成処理ができる --> 計算プログラムのインプットの詳細を知っていなくても計算ができる(もちろん知っていた方が利点は多い)
* Pythonをベースとして用いるため、for文によるループ処理やif文による分岐が簡単に導入できる

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

## VASPを用いる場合
* 第一原理計算としてメジャーなソフトウェアであるVASPを用いる場合は設定が必要
* VASPがインストールされていることが前提(実行ファイルは`~/vasp.X.X.X/bin/vasp_std`とする)
* ASEから参照するVASP実行用のPythonスクリプト(`run_vasp.py`とする)を用意する
```python
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`)を設定する
```bash
$ 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を用いる。それぞれ異なる計算手法となっており、計算コスト・計算精度が違ってくる
1. EMT (多体ポテンシャル)
2. DFT (VASP)
3. M3GNet (ニューラルネットワーク・ポテンシャル)
4. CHGNet (ニューラルネットワーク・ポテンシャル)

## 1. EMT
* Effective medium theory (EMT): 多体ポテンシャルの一種
* 対応している元素: Al, Cu, Ag, Au, Ni, Pd, Pt, H, C, N, O
```python
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計算からエネルギーを求める
```python
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`
    * matglはdgl(deep graph library)を使うが、PyTorchのバージョンが新しいと対応していないことがある。`torch==2.3.0`, `dgl==2.2.1`なら動いた(2025/2/22)
```python
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より精度が高い気がする
```python
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
```python
from ase.calculators.emt import EMT
from ase.build import fcc111
from ase.optimize import FIRE
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

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
calc = EMT()

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

# --- CHGNet
# potential = CHGNet.load()
# calc = CHGNetCalculator(potential=potential, properties="energy")

surf.calc = calc
        
opt = FIRE(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())`
* ASEでよく使われるoptimizerはBFGSだが、FIREやFIRE2のほうが素性がよいかもしれない
* `.traj`ファイルの可視化: `ase gui pt-relax.traj`

## 2. VASP
```python
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のみ)が欠点であり、コードの確認だけに用いた方がよい
* NNPを使うとより多くの元素を扱うことができ
* 他の機能は別ページで説明する予定