BAND Fragment Analysis Recipe

BAND Bonding Analysis EDA

Use the PLAMS BAND fragment recipe to run two fragment calculations plus a full-system calculation and extract energy decomposition terms.

Downloads: Notebook | Script ?

Related documentation
../_images/band_frag_4_0_e5d62dee.png

Initialization

from scm.plams import Settings, fromASE, view, Units
from scm.plams.recipes.bandfragment import BANDFragmentJob

# build the surface
from ase import Atoms
from ase.build import fcc111, add_adsorbate

Build Surface & Fragments

We first build a gold surface and add the hydrogen adsorbate.

mol = fcc111("Au", size=(2, 2, 3))
add_adsorbate(mol, "H", 1.5, "ontop")
mol.center(vacuum=10.0, axis=2)
view(fromASE(mol), padding=-4, width=400)
image generated from notebook

The system is then split into two fragments based on the species.

surface = mol.copy()
symbols = surface.get_chemical_symbols()
del surface[[i for i in range(len(symbols)) if symbols[i] != "Au"]]
adsorbate = mol.copy()
del adsorbate[[i for i in range(len(symbols)) if symbols[i] == "Au"]]

If available, optimized fragments can also be loaded.

# from ase import io
# surface_opt = io.read("surface_opt.xyz")
# adsorbate_opt = io.read("adsorbate_opt.xyz")
# assert len(surface_opt) == len(surface)
# assert len(adsorbate_opt) == len(adsorbate)

Set Up & Run Job

For efficiency in this example, we use a minimal basis and reduced computational details to run the job.

base_settings = Settings()
base_settings.input.ams.task = "SinglePoint"
base_settings.input.band.basis.type = "SZ"
base_settings.input.band.basis.core = "Large"
base_settings.input.band.dos.calcdos = "No"
base_settings.input.band.kspace.regular.numberofpoints = "3 3"
base_settings.input.band.beckegrid.quality = "Basic"
base_settings.input.band.zlmfit.quality = "Basic"
base_settings.input.band.usesymmetry = "No"
base_settings.input.band.xc.gga = "PBE"
base_settings.input.band.xc.dispersion = "Grimme4"

eda_settings = Settings()
eda_settings.input.band.peda = ""

eda_job = BANDFragmentJob(
    fragment1=fromASE(surface),
    fragment2=fromASE(adsorbate),
    settings=base_settings,
    full_settings=eda_settings,
    #    fragment1_opt=fromASE(surface_opt),
    #    fragment2_opt=fromASE(adsorbate_opt),
)

eda_job.run();

[23.03|17:02:01] JOB plamsjob STARTED
[23.03|17:02:01] JOB plamsjob RUNNING
[23.03|17:02:01] JOB plamsjob/frag1 STARTED
[23.03|17:02:01] JOB plamsjob/frag1 RUNNING
[23.03|17:02:53] JOB plamsjob/frag1 FINISHED
[23.03|17:02:53] JOB plamsjob/frag1 SUCCESSFUL
[23.03|17:02:53] JOB plamsjob/frag2 STARTED
[23.03|17:02:53] JOB plamsjob/frag2 RUNNING
[23.03|17:02:56] JOB plamsjob/frag2 FINISHED
[23.03|17:02:56] JOB plamsjob/frag2 SUCCESSFUL
[23.03|17:02:56] JOB plamsjob/full STARTED
[23.03|17:02:56] JOB plamsjob/full RUNNING
[23.03|17:03:36] JOB plamsjob/full FINISHED
[23.03|17:03:36] JOB plamsjob/full SUCCESSFUL
[23.03|17:03:36] JOB plamsjob FINISHED
[23.03|17:03:36] JOB plamsjob SUCCESSFUL

Energy Decomposition Results

Finally, we extract the results of the energy decomposition:

import pandas as pd

results = eda_job.results
eda_res = eda_job.results.get_energy_decomposition(unit="eV")

df = pd.DataFrame([{"Term": t, "Energy [eV]": e} for t, e in eda_res.items()])

df

Term

Energy [eV]

0

E_int

-1.886021

1

E_int_disp

-0.094696

2

E_Pauli

11.635861

3

E_elstat

-5.387854

4

E_orb

-8.039060

5

E_1

-19.288663

6

E_2

-0.018088

See also

Python Script

#!/usr/bin/env python
# coding: utf-8

# ## Initialization

from scm.plams import Settings, fromASE, view, Units
from scm.plams.recipes.bandfragment import BANDFragmentJob

# build the surface
from ase import Atoms
from ase.build import fcc111, add_adsorbate


# ## Build Surface & Fragments

# We first build a gold surface and add the hydrogen adsorbate.

mol = fcc111("Au", size=(2, 2, 3))
add_adsorbate(mol, "H", 1.5, "ontop")
mol.center(vacuum=10.0, axis=2)
view(fromASE(mol), padding=-4, width=400, picture_path="picture1.png")


# The system is then split into two fragments based on the species.

surface = mol.copy()
symbols = surface.get_chemical_symbols()
del surface[[i for i in range(len(symbols)) if symbols[i] != "Au"]]
adsorbate = mol.copy()
del adsorbate[[i for i in range(len(symbols)) if symbols[i] == "Au"]]


# If available, optimized fragments can also be loaded.

# from ase import io
# surface_opt = io.read("surface_opt.xyz")
# adsorbate_opt = io.read("adsorbate_opt.xyz")
# assert len(surface_opt) == len(surface)
# assert len(adsorbate_opt) == len(adsorbate)


# ## Set Up & Run Job

# For efficiency in this example, we use a minimal basis and reduced computational details to run the job.

base_settings = Settings()
base_settings.input.ams.task = "SinglePoint"
base_settings.input.band.basis.type = "SZ"
base_settings.input.band.basis.core = "Large"
base_settings.input.band.dos.calcdos = "No"
base_settings.input.band.kspace.regular.numberofpoints = "3 3"
base_settings.input.band.beckegrid.quality = "Basic"
base_settings.input.band.zlmfit.quality = "Basic"
base_settings.input.band.usesymmetry = "No"
base_settings.input.band.xc.gga = "PBE"
base_settings.input.band.xc.dispersion = "Grimme4"

eda_settings = Settings()
eda_settings.input.band.peda = ""

eda_job = BANDFragmentJob(
    fragment1=fromASE(surface),
    fragment2=fromASE(adsorbate),
    settings=base_settings,
    full_settings=eda_settings,
    #    fragment1_opt=fromASE(surface_opt),
    #    fragment2_opt=fromASE(adsorbate_opt),
)


eda_job.run()
# ## Energy Decomposition Results

# Finally, we extract the results of the energy decomposition:

import pandas as pd

results = eda_job.results
eda_res = eda_job.results.get_energy_decomposition(unit="eV")

df = pd.DataFrame([{"Term": t, "Energy [eV]": e} for t, e in eda_res.items()])

print(df)