原子分子シミュレーションPythonパッケージAtomic Simulation Environment (ASE)を使ってみる その１

Atomic Simulation Environment (ASE)というパッケージがあります。これは、第一原理計算や分子動力学法などの様々なパッケージをPythonによって統一的に操作できるフリーのソフトウエアです。
https://wiki.fysik.dtu.dk/ase/index.html

面白いのは、ある原子配置が与えられた時に、その時のエネルギーと原子間に働く力さえ返せる外部パッケージであれば、様々なものをPythonで扱うことができるようになります。有名どころでは、第一原理計算パッケージのVASPやQuantumEspressoなどがあります。
（このページ
https://wiki.fysik.dtu.dk/ase/ase/calculators/calculators.html
には、「ABINIT, AMBER, CP2K, CASTEP, deMon2k, DFTB+, ELK, EXCITING, FHI-aims, FLEUR, GAUSSIAN, Gromacs, Jacapo, LAMMPS, MOPAC, NWChem, Octopus, ONETEP, Quantum ESPRESSO, SIESTA, TURBOMOLE and VASP」と書いてあります）。
そして、エネルギーと力を計算できるのであれば、それらを使った分子動力学法ができてしまいます。

インストール

インストールにはPythonが必要です。Macの場合はcondaあるいはbrewなどでPython 2.7以上（Python 3系ももちろんOKです）を入れてください。
ASEはpipでインストールできて、

pip install --upgrade --user ase

で入れることができます。
https://wiki.fysik.dtu.dk/ase/install.html

チュートリアル　

チュートリアルを順番にやってみましょう。

銅への窒素の吸着エネルギー

https://wiki.fysik.dtu.dk/ase/tutorials/surface.html
窒素分子が銅のスラブへ吸着する際のエネルギーを見積もってみましょう。このエネルギーは、「独立した窒素分子のエネルギー＋銅原子スラブのエネルギー」と「銅原子スラブに窒素分子が吸着しているエネルギー」の差となります。

コードは以下に示します。適宜コメントを入れておきました。

from ase import Atoms
from ase.calculators.emt import EMT
from ase.constraints import FixAtoms
from ase.optimize import QuasiNewton
from ase.build import fcc111, add_adsorbate

h = 1.85
d = 1.10

slab = fcc111('Cu', size=(4, 4, 2), vacuum=10.0) #銅原子スラブのセット

slab.set_calculator(EMT()) #銅原子スラブの計算にはEMTを使用
e_slab = slab.get_potential_energy() #スラブのポテンシャルエネルギーを計算

molecule = Atoms('2N', positions=[(0., 0., 0.), (0., 0., d)]) #窒素分子のセット。(0,0,0)が一つ目のNの位置、(0,0,d)が二つ目のNの位置。
molecule.set_calculator(EMT()) #窒素分子の計算にはEMTを使用
e_N2 = molecule.get_potential_energy() #窒素分子のポテンシャルエネルギーの計算

add_adsorbate(slab, molecule, h, 'ontop') #窒素分子を上にのせる
constraint = FixAtoms(mask=[a.symbol != 'N' for a in slab]) #拘束条件としては、計算を高速化するため、銅原子の位置を緩和させずに固定
slab.set_constraint(constraint) #拘束条件をセット
dyn = QuasiNewton(slab, trajectory='N2Cu.traj') #準ニュートン法を設定
dyn.run(fmax=0.05) #構造緩和スタート。全ての原子に働く力がfmax以下になるまで。

print('Adsorption energy:', e_slab + e_N2 - slab.get_potential_energy()) #ポテンシャルエネルギーを計算し、先ほどの二つとの差を取る

となります。

Atomsオブジェクト

ASEではAtomsオブジェクトというものを作成し、そこに操作を加える形になります。上の例のように、分子やスラブの設定ができます。このAtomsオブジェクトに対して、計算方法を設定したり、エネルギーを得るように指示したりできます。

Calculators

Calculatorsは計算方法を指定します。

atom.set_calculator(calc)

とすると、atomというAtomsオブジェクトに対して、calcという計算方法をセットします。ここで、calcというのは様々なものを使うことができて、VASPやQuantum Espressoなどを選ぶ事ができます。

保存

作った原子の位置座標は

from ase.io import write
write('slab.xyz', slab)

で保存できます。

読み出しは

from ase.io import read
slab_from_file = read('slab.xyz')

です。

Julia実装

PyCall.jlの練習を兼ねて、Julia言語でも実装してみましょう。
PyCall.jlについては下の記事
「Juliaで機械学習：PyCall.jlを使ってTensorFlowのKerasを使ってみる」
https://qiita.com/cometscome_phys/items/e09e0300e303ba184bcc

using PyCall
const ase = pyimport("ase")
const Atoms = pyimport("ase")["Atoms"]
const EMT = pyimport("ase.calculators.emt")["EMT"]
const FixAtoms = pyimport("ase.constraints")["FixAtoms"]
const QuasiNewton = pyimport("ase.optimize")["QuasiNewton"]
const fcc111 = pyimport("ase.build")["fcc111"]
const add_adsorbate = pyimport("ase.build")["add_adsorbate"]

h = 1.85
d = 1.10

slab = fcc111("Cu", size=(4, 4, 2), vacuum=10.0)
slab[:set_calculator](EMT())
e_slab = slab[:get_potential_energy]()

molecule = Atoms("2N", positions=[(0., 0., 0.), (0., 0., d)])
molecule[:set_calculator](EMT())
e_N2 = molecule[:get_potential_energy]()

add_adsorbate(slab, molecule, h, "ontop")
constraint = FixAtoms(mask=[a[:symbol] != "N" for a in slab])
slab[:set_constraint](constraint)
dyn = QuasiNewton(slab, trajectory="N2Cu.traj")
dyn[:run](fmax=0.05)

println("Adsorption energy:", e_slab + e_N2 - slab[:get_potential_energy]())

JuliaでPythonのfrom a import bをやるときは、pyimport("a")["b"]となるようです。

分子動力学

チュートリアルで作った窒素分子に対して分子動力学法を実行することもできます。
その場合は、

from ase.md.verlet import VelocityVerlet
from ase import units
dyn = VelocityVerlet(molecule, dt=1.0 * units.fs)
for i in range(10):
     pot = molecule.get_potential_energy()
     kin = molecule.get_kinetic_energy()
     print('%2d: %.5f eV, %.5f eV, %.5f eV' % (i, pot + kin, pot, kin))
     dyn.run(steps=20)

となります。

Juliaバージョンは

using Printf
const VelocityVerlet = pyimport("ase.md.verlet")["VelocityVerlet"]
const units = pyimport("ase")["units"]
dyn = VelocityVerlet(molecule, dt=1.0 * units[:fs])
for i in 0:10-1
     pot = molecule[:get_potential_energy]()
     kin = molecule[:get_kinetic_energy]()
     @printf("%2d: %.5f eV, %.5f eV, %.5f eV \n",i, pot + kin, pot, kin)     
     dyn[:run](steps=20)
end

です。