\maketitle
There are a few different functions built into ase that let’s you create simple bulk systems. The one we will use is ase.build.bulk. Let us consider how to construct an fcc structure. You can either create the primitive cell (default), a orthorhombic cell, or a cubic cell. Note these will have different unit cell sizes and different number of atoms, which will affect the computational cost of your simulation.
from ase.build import bulk
from vasp import *
import matplotlib.pyplot as plt
from ase.visualize import view
# a is the lattice constant in Angstroms
Pd_cubic = bulk('Pd', 'fcc', a = 3.89, cubic = True)
view(Pd_cubic)
print(Pd_cubic)
Pd_ortho = bulk('Pd', 'fcc', a = 3.89, orthorhombic = True)
view(Pd_ortho)
print(Pd_ortho)
Pd_primitive = bulk('Pd', 'fcc') # if a is not specified, default experimental value is used.
view(Pd_primitive)
print(Pd_primitive)Calculations can be done in a similar way to what we performed in the last lab and homework, with a few additional parameters. The most important of these is the k-point grid.
For modeling bulk systems/surfaces, we need to specify a k-point mesh. A few things to note:
- The accuracy (and cost) of your simulation usually goes up with increasing the k-point grid used.
- In this lab (for demonstrative purposes only) a k-point grid of (8x8x8) has been used. We typically need to do a convergence study for k-points, similar to the one we did for the planewave cutoff.
- The number of k-points required scales inversely with the size of the unit cell.
- We will usually specify a parameter
kpts=(k1,k2,k3)invasp, which creates aKPOINTSfile that tellsVASPto automatically create a Monkhorst-Pack k-mesh. - Metals (conductors) generally require more k-points than insulators to reach the same level of convergence.
ase by default uses the experimental lattice constants (if known). Our computationally calculated lattice constants will usually be slightly different, depending on the exchange-correlation functional used. Let us consider a series of lattice constants, to find out which one minimizes the total energy.
import numpy as np
from ase.utils.eos import EquationOfState
from vasp.vasprc import VASPRC
VASPRC['queue.q'] = '*long'
VASPRC['queue.nprocs'] = 24
VASPRC['queue.pe'] = 'mpi-24'
# Generate an array of 7 points 10% around 3.89 angstroms
A = np.linspace(3.89 * 0.9, 3.89 * 1.1, 7)
energies = []
volumes = []
ready = True
for a in A:
# We will use the cubic cell for simplicity
Pd_cubic = bulk('Pd', 'fcc', a = a, cubic = True)
calc = Vasp('EOS/Pd-a-{0:1.2f}'.format(a),
xc = 'PBE',
encut = 400,
ismear = 1, # Use MP smearing for metlas
kpts = [8,8,8],
atoms = Pd_cubic)
calc.calculate()
energies.append(Pd_cubic.get_potential_energy())
volumes.append(Pd_cubic.get_volume())
if None in energies:
ready = False
if ready == True:
plt.plot(A, energies, 'bo-')
plt.xlabel('Lattice constant ($\AA$)')
plt.ylabel('Total energy (eV)')
plt.savefig('images/Pd-fcc-lattice.png')
plt.show()
To find the ‘optimal’ lattice constant we need to fit our data to an equation of state, which describes the energy as a function of volume. The Murnaghan or Birch-Murnaghan EOS is commonly used. Let us use ase.eos.EquationOfState to fit the data we calculated above to the Birch-Murnaghan EOS.
# Lets fit this to an equation of state
eos = EquationOfState(volumes, energies, eos='birchmurnaghan')
v0, e0, b = eos.fit()
a0 = v0 **(1/3)
eos.plot(filename='images/Pd-EOS.png', show=True)
print('Minimum Energy = {0:1.3f} eV'.format(e0))
print('Optimal Volume = {0:1.3f} cubic angstroms'.format(v0))
print('Optimal lattice constant = {0:1.3f} angstroms'.format(a0))
ase provides functions to create surfaces too. Surfaces are layers of atoms formed by cleaving the bulk structure in a given direction. In our models, we add vacuum in the direction perpendicular to the surface. Thus, the atoms are finite in the direction perpendicular to the surface, but infinite in the other two directions. Here is an example of how to make a surface.
from ase.build.surface import fcc111
from ase.constraints import FixAtoms
from ase.io import write
a = 3.932 # Optimal lattice constant from EOS
# Create a surface with 3 unit cells in x and y
# 3 layers deep
atoms = fcc111('Pd', size=(2,2,3), vacuum=10.0, a=a)
view(atoms)
for atom in atoms:
print(atom.symbol, atom.tag)
write('images/Pd-slab.png', atoms, rotation='90x', show_unit_cell=2)
The tag on the atom indicates which layer of the surface it is in.
We can see that there are actually two surfaces, one in the top layer and one at the bottom layer. Surface atoms will tend to contract toward the bulk due to decreased coordination.
To simulate bulk like behavior in regions away from the surface, we can do two things:
- increase the number of layers in the slab (requires many atoms, large cost)
- Constrain(freeze) the the bottom layer(s) in their bulk positions (common, lower cost). The bottom layer is now representative of bulk behavior.
Let us now optimize the geometry of our surface. Note that only one /k/-point is required in the direction perpendicular to the surface.
from ase.lattice.surface import fcc111
a = 3.9322 # Optimal lattice constat from EOS
atoms = fcc111('Pd', size = (2,2,3), vacuum = 10.0, a=a)
constraint = FixAtoms(mask = [atom.tag >= 3 for atom in atoms])
atoms.set_constraint(constraint)
calc = Vasp('surfaces/Pd-slab-relaxed',
xc = 'PBE',
ismear = 1,
kpts = [8,8,1],
encut = 400,
ibrion = 2,
nsw = 20,
atoms = atoms)
calc.calculate()
print(calc): -----------------------------
VASP calculation from /afs/crc.nd.edu/user/p/pmehta1/computational-chemistry/Lab4/surfaces/Pd-slab-relaxed
converged: True
Energy = -58.019294 eV
Unit cell vectors (angstroms)
x y z length
a0 [ 5.561 0.000 0.000] 5.561
a1 [ 2.780 4.816 0.000] 5.561
a2 [ 0.000 0.000 24.540] 24.540
a,b,c,alpha,beta,gamma (deg):5.561 5.561 24.540 90.0 90.0 90.0
Unit cell volume = 657.154 Ang^3
Stress (GPa):xx, yy, zz, yz, xz, xy
0.012 0.012 0.000-0.000 -0.000 -0.000
Atom# sym position [x,y,z]tag rmsForce constraints
0 Pd [1.390 0.803 10.000] 3 0.00 F F F
1 Pd [4.171 0.803 10.000] 3 0.00 F F F
2 Pd [2.780 3.210 10.000] 3 0.00 F F F
3 Pd [5.561 3.210 10.000] 3 0.00 F F F
4 Pd [5.561 1.605 12.270] 2 0.04 T T T
5 Pd [2.780 1.605 12.270] 2 0.04 T T T
6 Pd [6.951 4.013 12.270] 2 0.04 T T T
7 Pd [4.171 4.013 12.270] 2 0.04 T T T
8 Pd [0.000 0.000 14.548] 1 0.02 T T T
9 Pd [2.780 0.000 14.548] 1 0.02 T T T
10 Pd [1.390 2.408 14.548] 1 0.02 T T T
11 Pd [4.171 2.408 14.548] 1 0.02 T T T
--------------------------------------------------
INCAR Parameters:
-----------------
nbands: 72
ismear: 1
nsw: 20
ibrion: 2
encut: 400.0
magmom: None
kpts: (8, 8, 1)
reciprocal: False
xc: PBE
txt: -
gamma: False
Pseudopotentials used:
----------------------
Pd: potpaw_PBE/Pd/POTCAR (git-hash: 04426435b178dfad58ed91b470847d50ff70b858)
Note that Vasp is a little unintuitive. The constraint ‘F’ means frozen.
We can go back to the calculation directory and see how our surface relaxed with jaspsum -t.
Now let us add an adsorbate on our surface. There are multiple places where it could adsorb. Here is a picture of a fcc(111) gold surface, showing the possible adsorption sites.
Let’s go back to our Pd surface and perform a calculation with an Oxygen adsorbate at the fcc site.
from ase.build.surface import fcc111, add_adsorbate
a = 3.932 # Optimal lattice constant from EOS
atoms = fcc111('Pd', size=(2,2,3), vacuum=10.0)
add_adsorbate(atoms, 'O', height=1.2, position='fcc')
# Note that constraints only work after adding the adsorbate
constraint = FixAtoms(mask=[atom.tag >= 3 for atom in atoms])
atoms.set_constraint(constraint)
view(atoms)
calc = Vasp('surfaces/O-on-Pd-fcc',
xc='PBE',
ismear=1,
kpts=[8, 8, 1],
encut=400,
ibrion=2, # Conjugate Gradient
nsw=20, # relaxation steps
atoms=atoms)
calc.calculate()
print(calc)
write('images/Pd-slab-O-fcc.png', atoms, show_unit_cell=2) : -----------------------------
VASP calculation from /afs/crc.nd.edu/user/p/pmehta1/computational-chemistry/Lab4/surfaces/O-on-Pd-fcc
converged: True
Energy = -64.436725 eV
Unit cell vectors (angstroms)
x y z length
a0 [ 5.501 0.000 0.000] 5.501
a1 [ 2.751 4.764 0.000] 5.501
a2 [ 0.000 0.000 24.492] 24.492
a,b,c,alpha,beta,gamma (deg):5.501 5.501 24.492 90.0 90.0 90.0
Unit cell volume = 641.919 Ang^3
Stress (GPa):xx, yy, zz, yz, xz, xy
0.004 0.004 0.002-0.000 -0.000 -0.000
Atom# sym position [x,y,z]tag rmsForce constraints
0 Pd [1.375 0.794 10.000] 3 0.00 F F F
1 Pd [4.126 0.794 10.000] 3 0.00 F F F
2 Pd [2.751 3.176 10.000] 3 0.00 F F F
3 Pd [5.501 3.176 10.000] 3 0.00 F F F
4 Pd [5.505 1.586 12.288] 2 0.02 T T T
5 Pd [2.747 1.586 12.288] 2 0.02 T T T
6 Pd [6.877 3.970 12.405] 2 0.03 T T T
7 Pd [4.126 3.975 12.288] 2 0.02 T T T
8 Pd [-0.031 -0.018 14.660] 1 0.01 T T T
9 Pd [2.782 -0.018 14.660] 1 0.01 T T T
10 Pd [1.375 2.418 14.660] 1 0.01 T T T
11 Pd [4.126 2.382 14.541] 1 0.01 T T T
12 O [1.375 0.794 15.820] 0 0.02 T T T
--------------------------------------------------
INCAR Parameters:
-----------------
nbands: 80
ismear: 1
nsw: 20
ibrion: 2
encut: 400.0
magmom: None
kpts: (8, 8, 1)
reciprocal: False
xc: PBE
txt: -
gamma: False
Pseudopotentials used:
----------------------
O: potpaw_PBE/O/POTCAR (git-hash: 592f34096943a6f30db8749d13efca516d75ec55)
Pd: potpaw_PBE/Pd/POTCAR (git-hash: 04426435b178dfad58ed91b470847d50ff70b858)
The adsorption energy is given by $Eads = Esurface+O - Esurface - 0.5 EO_{2}$. This can easily be calculated from the two calculations we performed and the O2 calculation from the last homework.
It is possible to plot out the density of states (DOS) from VASP calculations. The density of states describes the number of states per interval of energy at each energy level that are available to be occupied (Wikipedia).
We can get the total density of states from an old DFT calculation without having to run a new calculation (Though the VASP manual recommends an additional run at ismear=-5). The DOS depends on the k-point grid you choose.
Let’s read in our calculation from our bulk Pd lattice constant studies.
from ase.dft.dos import DOS
calc = Vasp('./EOS/Pd-a-3.89')
dos = DOS(calc, width = 0.2)
energies = dos.get_energies()
dos = dos.get_dos()
plt.plot(energies, dos, linewidth = 2)
# Add a vertical line at the fermi level
plt.axvline(x=0, color = 'r', linestyle='--', linewidth = 2)
plt.ylim(0,12)
plt.xlim(-10,10)
plt.xlabel('energy (eV)')
plt.ylabel('DOS (arb. units)')
plt.savefig('images/Pd-bulk-dos.png')
plt.show()
States to the left of the fermi level (indicated by the red line) are the occupied states.
To figure out which density of states belong to which atoms in a molecule, we need to perform an additional calculation. We can compute the atom-projected density of states (ADOS), which is done by projecting the wave function onto localized atomic orbitals. These are only a qualitative representation of the orbitals, because the atoms will often form molecular orbitals, hybridize, etc.
In VASP we can specify an RWIGS parameter, which is the radius around the atom at which to cutoff the projection. The choice of RWIGS is somewhat arbitrary, one can choose the ionic radius of an atom, or a value that minimizes overlap of neighboring spheres. Another way to calculate the ADOS is by specifying the LORBIT parameter to be 10 or 11, but this only works for PAW potentials (this is what we will use).
In transition metals, the s and p states are dispersed, and the only states that matter in terms of bonding are the d-states. Here is an example to plot the DOS projected on to the d states for clean Pd surface atoms.
# get the geometry the previous calculation
calc = Vasp('surfaces/Pd-slab-relaxed')
atoms = calc.get_atoms()
#Now submit a calculation for the ados
calc = Vasp('surfaces/Pd-ados',
xc='PBE',
ismear=1,
kpts=[8, 8, 1],
encut=400,
lorbit=10,
atoms=atoms)
calc.calculate()
# Atom index 10 is a surface atom (tag=1)
print(atoms[10])
energies, d_dos = calc.get_ados(10, 'd')
plt.plot(energies, d_dos, lw=2)
plt.axvline(lw=2, ls='--', color='r')
plt.xlabel('energy (eV)')
plt.ylabel('DOS (arb. units)')
plt.savefig('images/Pd-ados.png')
plt.show()
Now let us plot the density of states for the adsorbed O atom.
calc = Vasp('surfaces/O-on-Pd-fcc')
atoms = calc.get_atoms()
calc = Vasp('surfaces/O-on-Pd-fcc-ados',
xc = 'PBE',
ismear = 1,
kpts = [8,8,1],
encut = 400,
lorbit = 10,
atoms = atoms)
calc.calculate()
s_energies, s_dos = calc.get_ados(12,'s')
p_energies, p_dos = calc.get_ados(12,'p')
plt.plot(s_energies, s_dos, label='O$_{s}$', lw=2)
plt.plot(p_energies, p_dos, label='O$_{p}$', lw=2)
# Now plot the clean surface ados for comparison
calc = Vasp('surfaces/Pd-ados')
clean_energies, d_dos = calc.get_ados(11,'d')
plt.plot(clean_energies, d_dos, label='Pd$_{d}$', lw=2)
plt.xlim(-20,8)
plt.axvline(ls='-.', color='k', lw=2)
plt.xlabel('energy (eV)')
plt.ylabel('DOS (arb.units)')
plt.legend()
plt.savefig('images/adsorbate-dos.png')
plt.show()
The blue line indicates the Oxygen s-states. The two peaks of the green line left and right of the Pd d-band are the bonding and antibonding Oxygen p-states. Note that the antibonding peak is to the right of the fermi level, meaning that the antibonding states are unoccupied.