【Atomic Simulation Environment文檔】使用ASE庫進行表面吸附研究


在本教程中我們將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=symbads=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的參雜體系,因此上述代碼的計算結果並不具備實際意義。


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM