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が必要になってくるだろう。