在本教程中我們將C、N、O三種原子分別吸附在具有1層,2層和3層的7種不同FCC金屬(111)晶面上,然后使用數據庫文件儲存結果。
大塊晶體(Bulk)
首先,我們使用EMT勢,分別計算7種元素的體相FCC晶格常數。
from ase.build import bulk
from ase.calculators.emt import EMT
from ase.eos import calculate_eos
from ase.db import connect
db = connect('bulk.db')
for symb in ['Al', 'Ni', 'Cu', 'Pd', 'Ag', 'Pt', 'Au']:
atoms = bulk(symb, 'fcc')
# 尋找能量最小時的晶格常數
atoms.calc = EMT()
eos = calculate_eos(atoms)
v, e, B = eos.fit()
# 在最優體系盡心計算,寫入數據庫中
atoms.cell *= (v / atoms.get_volume()) ** (1 / 3)
atoms.get_potential_energy()
db.write(atoms, bm=B)
運行bulk.py腳本,得到如下結果:
$ python3 bulk.py
$ ase db bulk.db -c +bm # show also the bulk-modulus column
id|age|formula|calculator|energy| fmax|pbc|volume|charge| mass| bm
1|10s|Al |emt |-0.005|0.000|TTT|15.932| 0.000| 26.982|0.249
2| 9s|Ni |emt |-0.013|0.000|TTT|10.601| 0.000| 58.693|1.105
3| 9s|Cu |emt |-0.007|0.000|TTT|11.565| 0.000| 63.546|0.839
4| 9s|Pd |emt |-0.000|0.000|TTT|14.588| 0.000|106.420|1.118
5| 9s|Ag |emt |-0.000|0.000|TTT|16.775| 0.000|107.868|0.625
6| 9s|Pt |emt |-0.000|0.000|TTT|15.080| 0.000|195.084|1.736
7| 9s|Au |emt |-0.000|0.000|TTT|16.684| 0.000|196.967|1.085
Rows: 7
Keys: bm
$ ase gui bulk.db
bulk.db是一個保存在一個簡單文件中的SQLite3數據庫
$ file bulk.db
bulk.db: SQLite 3.x database
如果你想要看到文件中的內容,你可以將數據庫文件轉換為json文件,然后使用文本編輯器將其打開:
$ ase db bulk.db --insert-into bulk.json
Added 0 key-value pairs (0 pairs updated)
Inserted 7 rows
或者,你可以像下面一樣查看某一行的內容:
$ ase db bulk.db Cu -j
{"1": {
"calculator": "emt",
"energy": -0.007036492048371201,
"forces": [[0.0, 0.0, 0.0]],
"key_value_pairs": {"bm": 0.8392875566787444},
...
...
}
json文件格式是人類可以直接閱讀的格式,但其與SQLite3文件相比,處理效率要低得多;
吸附體系(adsorbates)
現在我們進行吸附計算(運行ads.py腳本)
from ase.calculators.emt import EMT
from ase.db import connect
from ase.build import fcc111, add_adsorbate
from ase.constraints import FixAtoms
from ase.optimize import BFG5
db1 = connect('bulk.db')
db2 = connect('ads.db')
def run(symb, a, n, ads):
"""
symb:吸附表面元素種類(字符)
a :以晶胞為單位指定系統的大小(元組)
n :晶面常數(取1,2,3)
ads :吸附質元素種類(字符)
"""
atoms = fcc111(symb, (1, 1, n), a=a)
add_adsorbate(atoms, ads, height=1.0, position='fcc')
# 將除吸附質以外的所有原子約束住
fixed = list(range(len(atoms) - 1))
atoms.constraints = [FixAtoms(indices=fixed)]
atoms.calc = EMT()
opt = BFG5(atoms, logfile=None)
opt.run(fmax=0.01)
return atoms
for row in db1.select():
a = row.cell[0, 1] * 2
symb = row.symbols[0]
for n in [1, 2, 3]:
for ads in 'CNO':
atoms = run(symb, a, n, ads)
db2.write(atoms, layers=n, surf=symb, ads=ads)
我們現在擁有了一個新的數據庫文件,包括63行:
$ ase db ads.db -n
63 rows
如果使用EMT方法,這63次計算只需要很小一段時間。假設你想要在超算上使用DFT方法進行計算,在這種情況下你也許希望在計算機上的不同jobs上運行多個計算。另外,有些工作可能會因為超時而無法完成。在這種情況下,最好對腳本進行一些修改,我們在內循環中添加一行:
for row in db1.select():
a = row.cell[0, 1] * 2
symb = row.symbols[0]
for n in [1, 2, 3]:
for ads in 'CNO':
id = db2.reserve(layers=n, surf=symb, ads=ads)
if id is not None:
atoms = run(symb, a, n, ads)
db2.write(atoms, layers=n, surf=symb, ads=ads)
del db2[id]
上面提到的reserve()方法會檢測一行里面是否包含關鍵詞layer=n,surf=symb和ads=ads。如果包含,計算會被跳過;如果不包含,然后包含這些關鍵詞的空行會被寫入到文件中,計算開始。計算完成后,包含計算結果的內容寫入,空行被刪除。修改后的腳本可以在多個並行運行的作業中運行,並且不會進行兩次計算。
萬一計算崩潰,你可以看到數據庫中被保留下來的空行:
$ ase db ads.db natoms=0 -c ++
id|age|user |formula|pbc|charge| mass|ads|layers|surf
17|31s|jensj| |FFF| 0.000|0.000| N| 1| Cu
Rows: 1
Keys: ads, layers, surf
刪除這些內容,解決問題,然后再次運行腳本:
$ ase db ads.db natoms=0 --delete
Delete 1 row? (yes/No): yes
Deleted 1 row
$ python ads.py # or sbatch ...
$ ase db ads.db natoms=0
Rows: 0
參考能量
接着,讓我們計算清潔吸附表面和孤立吸附質的能量:
from ase import Atoms
from ase.calculators.emt import EMT
from ase.db import connect
from ase.build import fcc111
db1 = connect('bulk.db')
db2 = connect('ads.db')
def run(symb, a, n):
atoms = fcc111(symb, (1, 1, n), a=a)
atoms.calc = EMT()
atoms.get_forces()
return atoms
# 清潔表面體系能量計算
for row in db1.select():
a = row.cell[0, 1] * 2
symb = row.symbols[0]
for n in [1, 2, 3]:
id = db2.reserve(layers=n, surf=symb, ads='clean')
if id is not None:
atoms = run(symb, a, n)
db2.write(atoms, id=id, layers=n, surf=symb, ads='clean')
# 孤立原子
for ads in 'CNO':
a = Atoms(ads)
a.calc = EMT()
a.get_potential_energy()
db2.write(a)
假設我們想要上述24種參比能量(包括清潔表面X21和孤立吸附質X3)保存在refs.db文件中,而不是較大的ads.db文件。我們可以改變上面的refs.db腳本,重新計算。我們也可以使用ase db工具操作得到的文件,首先,移動清潔表面體系:
$ ase db ads.db ads=clean --insert-into refs.db
Added 0 key-value pairs (0 pairs updated)
Inserted 21 rows
$ ase db ads.db ads=clean --delete --yes
Deleted 21 rows
然后,是三個孤立原子(pdc=FFF,無周期性)
$ ase db ads.db pbc=FFF --insert-into refs.db
Added 0 key-value pairs (0 pairs updated)
Inserted 3 rows
$ ase db ads.db pbc=FFF --delete --yes
Deleted 3 rows
$ ase db ads.db -n
63 rows
$ ase db refs.db -n
24 rows
分析(Analysis)
現在我們擁有了計算吸附能和吸附高度所需要的數據:
for ase.db import connect
refs = connect('refs.db')
db = connect('ads.db')
for row in db.select():
ea = (row.energy -
refs.get(formula=row.ads).energy -
refs.get(layers=row.layers, surf=row.surf).energy)
h = row.position[-1, 2] - row.position[-2, 2]
db.update(row.id, height=h, ea=ea)
在這里我們展示三層Pt金屬FCC(111)晶面上吸附能的數據:
$ python3 ea.py
$ ase db ads.db Pt,layers=3 -c formula,ea,height
formula| ea|height
Pt3C |-3.715| 1.504
Pt3N |-5.419| 1.534
Pt3O |-4.724| 1.706
Rows: 3
Keys: ads, ea, height, layers, surf
Note:
EMT方法可以很好地描述對Ni,Cu,Pd,Ag,Pt,Au和Al純金屬晶體,但不適用於C、N或O的參雜體系,因此上述代碼的計算結果並不具備實際意義。
