ASE: 窒素分子のCu(111)の二原子層における吸着構造の最適化

Python module のASE を使って、吸着構造を最適化。
ここでは簡単な有効媒質理論を使うが、calculatorを選ぶことで、AbinitとかVASPとかも使えるらしい。
Introduction: Nitrogen on copper — ASE documentation
元サイトのチュートリアルに少し追加して、吸着した時に結合距離がどう変わるかを見る。
吸着エネルギーを出しているが、ここではあまり関係ない。

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
from ase.visualize import view

h = 1.85 # initial adsorption distance
d = 1.10 # initial bond distance

# make Cu structure
slab = fcc111('Cu', size=(4, 4, 2), vacuum=10.0)
view( slab )

# calculate Cu energy
slab.set_calculator(EMT())
e_slab = slab.get_potential_energy()

# make N2 structure
molecule = Atoms('2N', positions=[(0., 0., 0.), (0., 0., d)])
molecule.set_calculator(EMT())
##e_N2 = molecule.get_potential_energy()

# optimise N2 structure
dyn = QuasiNewton(molecule, trajectory='N2.traj')
dyn.run(fmax=0.05)
print( molecule.get_all_distances() )

# calculate N2 energy
e_N2 = molecule.get_potential_energy()

# make a whole structure
add_adsorbate(slab, molecule, h, 'ontop')
view( slab )
print( [ dis for dis in slab.get_all_distances()[-2] if dis < 3. ] ) #  distance Cu-N1, N1-N1, N1-N2

# optimize adsorption
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)

print('Adsorption energy:', e_slab + e_N2 - slab.get_potential_energy())
view( slab )
print( [ dis for dis in slab.get_all_distances()[-2] if dis < 3. ] ) #  distance Cu-N1, N1-N1, N1-N2

add_adsorbateで足された原子群はリストの後ろから入る。
距離を出力する部分は、ちょっとサボって、3オングストローム以下だけ出力するようにした。
本当は、短い距離順にソートして、そこから三つ持ってくるのが良いと思う。

結果：
[1.8499999999999996, 0.0, 0.99809885325797332]
[2.0610952626456882, 0.0, 1.0480140019810147]

1.85は適当に決めた吸着距離だからどうでもよい。
N-Nの距離に注目すると、吸着前は 0.998 だったが、吸着後は 1.048 に伸びることが予言された。
「Cuに近い方の窒素原子の方がより引っ張られるから、結合が伸びる」というイメージを支持するものと思われる。

今はN2分子がCu表面に対して垂直だが、横にした時どうかとか、そういう比較をし始めるときには、Adsorption energyが必要になってくるだろう。