Ctrl+K

Intro to ASE Calculators

2.3. Intro to ASE Calculators#

In this lesson, we will cover the basics of ASE calculators and how to use them. ASE has the ability to attach a “calculator” to an atoms object. At a fundamental level, these are just classes that will return energies and forces. That means it can be a DFT program, a quantum chemistry program, or even just a classical potential like “Leonard-Jones” There are a ton of calculators implmented. Here we will just use a simple one, the EMT calculator. A calculator acts as a black box that takes in atomic numbers and positions to give out energy, forces and stresses.

The EMT calculator is just a simple pair potential calculator for a few metals. In practice it is a toy calculator used for testing.

You start by making an instance of your calculator, then using set_calculator to attach it to an atoms object. Once that is done you can call get_potential_energy and get_forces to calculate the energy and forces.

import matplotlib.pyplot as plt

from ase.calculators.emt import EMT

calc = EMT()
water.set_calculator(calc)
energy = water.get_potential_energy()
forces = water.get_forces()

print(energy)
print(forces)

In our group we use quatum espresso, VASP, and Abinit to do DFT in general.

A primary use of calculators is to perform structural optimizations. This allows us to find the “lowest energy configuration” of a given structure. ASE has tools to do this. Below we’re using the BFGSLineSearch optimization method

from ase.optimize import BFGSLineSearch

relax = BFGSLineSearch(atoms = water)
relax.run(fmax = 0.05) # relax the structure until the maximum force is 0.05 eV/A

view(water)

2.3.1. Example:#

Use the EMT calculator to plot the potential energy as a function of distance between two H atoms.

#Solution to the exercise:

import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
rs = np.linspace(0.4, 3, 50)
engs = []
## calculate energy as a function of r
for r in rs:
    atoms = molecule('H2')
    atoms.set_distance(0,1,r)
    atoms.set_calculator(EMT())
    engs.append(atoms.get_potential_energy())
    
plt.plot(rs, engs)

2.3.2. Calculating Reaction Energies#

Reaction energies are the core of computational catalysis and surface science because they provide fundamental information about the energy for a chemical reaction. This section is meant to cover the basics of how to calculate a reaction energy using DFT.

The adsorption energy is defined as the energy difference between the combined system and the separate systems:

\(E_{rxn} = \Sigma E_{products} - \Sigma E_{reactants}\)

here we’re going to calculate the reaction energy of forming water for O\(_2\) and H\(_2\)

#We'll start by making an O$_2$ molecule
from ase.build import molecule

O2 = molecule('O2')

#let's visualize it next
from ase.visualize import view

view(O2)

#To calculate the energy of this molecule, we're going to use the EMT calculator
from ase.calculators.emt import EMT

O2.set_calculator(EMT())

#Next we need to optimize the structure, because chances are it's not in its lowest energy configuration
from ase.optimize import QuasiNewton

dyn = QuasiNewton(O2)
dyn.run(fmax=0.05)


E_O2 = O2.get_potential_energy()

Next, let’s build and optimize H\(_2\) and H\(_2\)O. ASE has some nice tools for doing this.

H2 = molecule('H2')
view(H2, viewer='nglview')

H2.set_calculator(EMT())
dyn = QuasiNewton(H2)
dyn.run(fmax=0.05)

E_H2 = H2.get_potential_energy()

H2O = molecule('H2O')
view(H2O, viewer='nglview')

H2O.set_calculator(EMT())
dyn = QuasiNewton(H2O)
dyn.run(fmax=0.05)

E_H2O = H2O.get_potential_energy()

E_rxn = E_H2O - E_H2 - 0.5*E_O2
print(E_rxn) # answer is in eV

Note that this is not very accurate! The real value is approximately -285.8261 kJ/mol. The reason is that we used the “EMT” calculator, which is not appropriate for the physics of molecules. However, it is very fast, so it allows us to see the principles behind how reaction energies can be computed. If we want accurate numbers, we need to use quantum mechanics, approximated through Density Functional Theory (DFT). This will be much more computationally expensive, and will require a supercomputer. You will learn more about this in future lessons.

from ase.units import kJ, mol
true_value = -285.8261 * kJ / mol
print(true_value)

error = (E_rxn-true_value)/true_value
print(error)

2.3.3. More about ASE Calculators#

2.3.3.1. Supported calculators in ASE#

ASE supports a variety of calculators. List of calculators either available through native ASE interfaces or python implementations in ASE package can be found here. The Calculator takes in keywords which can be parsed in by user based on the type of calculations they want to perform. Some example of keywords for a DFT calculator are: xc (exchange-correlation functional), kpts (k-points), smearing, nbands. Each calculator will have different native keywords specific to that calculator. These arguments can also be set or changed later by using the set() method.

2.3.3.2. Calculator interface in ASE#

An ASE calculator will store a copy of atoms object that is used for calculation. To print out the properties after a calculation has been performed, one can use get_potential_energy(), get_forces(), get_stress() methods.

Contents