Search for adsorption structure using Minima Hopping method

Using the Minima Hopping method implemented in ASE (Atomic Simulation Environment), we will explore the adsorption structure of H _{2 </ sub> molecules on the Pt surface. This article was created with reference to the Official Documentation Tutorial.}

The reason why EMT is used is because of the calculation speed. Especially for organic molecules (C, H, O, N), as mentioned in the Official Documents, It's just a play model. It is not reliable. Select an appropriate theory when using it in practice.

Create a slab model

Create a model in which H _{2 </ sub> molecules are adsorbed on the surface of Pt (211).}

from ase.build import surface, add_adsorbate
from ase.data import atomic_numbers, covalent_radii
from ase.visualize import view

#Pt bulk
Pt_bulk= bulk('Pt', 'fcc', a=3.9, cubic=True) 

# Pt(211)Slab
Pt_211 = surface(Pt_bulk, (2,1,1), 3, vacuum=10, periodic=True) 
d_PtH = (covalent_radii[atomic_numbers['H']] + covalent_radii[atomic_numbers['Pt']])
z_max =  np.max(Pt_211.positions[:,2]) #Maximum z-coordinate value of Pt atom
xy_center = ((Pt_211.cell[0]+Pt_211.cell[1])/2)[0:2]

#Add H2 molecule
add_adsorbate(Pt_211, molecule('H2'), d_PtH, xy_center)

#Check the structure
view(Pt_211)

Set constraint conditions

Use the FixAtoms and Hookean classes to set constraints. Under the constraint condition of Hookean, a spring-like force is applied between atoms.

--Pt slab is completely fixed --H _{2 </ sub> Molecules should not break and be not too far from the surface}

from ase.calculators.emt import EMT
from ase.constraints import FixAtoms, Hookean
from ase.optimize.minimahopping import MinimaHopping

#Restraint condition
const_Pt = [FixAtoms(indices=[atom.index for atom in Pt_211 if atom.symbol=='Pt'])]
d_HH = covalent_radii[atomic_numbers['H']] * 2
H_indices = [atom.index for atom in Pt_211 if atom.symbol=='H']
const_H2 = [Hookean(a1=H_indices[0], a2=H_indices[1], rt=d_HH*1.3, k=15.), #rt is the distance threshold[A],k is the spring constant
            Hookean(a1=H_indices[0], a2=(0., 0., 1., -(z_max+d_PtH*5)), k=5.),
            Hookean(a1=H_indices[1], a2=(0., 0., 1., -(z_max+d_PtH*5)), k=5.)]
surf.set_constraint(const_Pt + const_H2)
surf.set_calculator(EMT())

Run Minima Hopping

Execute Minima Hopping. By default, the log is output to the hop.log file.

hop = MinimaHopping(surf, Ediff0=0.5, T0=2000)
hop(totalsteps=10)

Visualization of results

The ʻase.optimize.minimahoppingmodule implements theMHPlot` function that visualizes the result of MinimaHopping. It is a wrapper for matplotlib. A black check mark indicates the newly found tiny structure.

import matplotlib.pyplot as plt
from ase.optimize.minimahopping import MHPlot

plt.rcParams["font.size"] = 18 #Adjust font size
MHPlot()

Let's also check the tiny structure found by Minima Hopping.

reference

Constrained minima hopping (global optimization) https://wiki.fysik.dtu.dk/ase/tutorials/minimahopping/minimahopping.html