Atoms and calculators

ASE allows atomistic calculations to be scripted with different computational codes. In this introductory exercise, we go through the basic concepts and workflow of ASE and will eventually calculate the binding curve of N2.

These tutorials often use the electronic structure code GPAW. They can be completed just as well using other supported codes, subject to minor adjustments.

Python

In ASE, calculations are performed by writing and running Python scripts. A very short primer on Python can be found in the ASE documentation. If you are new to Python it would be wise to look through this to understand the basic syntax, datatypes, and things like imports. Or you can just wing it — we won’t judge.

Atoms

Let’s set up a molecule and run a DFT calculation. We can create simple molecules by manually typing the chemical symbols and a guess for the atomic positions in Ångström. For example N2:

from ase import Atoms
atoms = Atoms('N2', positions=[[0, 0, -1], [0, 0, 1]])

Just in case we made a mistake, we should visualize our molecule using the ASE GUI:

from ase.visualize import view
view(atoms)

Equivalently we can save the atoms in some format, often ASE’s own trajectory format:

from ase.io import write
write('myatoms.traj', atoms)

Then run the GUI from a terminal:

$ ase gui myatoms.traj

ASE supports quite a few different formats. For the full list, run:

$ ase info --formats

Although we won’t be using all the ASE commands any time soon, feel free to get an overview:

$ ase --help

Exercise

Write a script which sets up and saves an N2 molecule, then visualize it.

Calculators

Next let us perform an electronic structure calculation. ASE uses calculators to perform calculations. Calculators are abstract interfaces to different backends which do the actual computation. Normally, calculators work by calling an external electronic structure code or force field code. To run a calculation, we must first create a calculator and then attach it to the Atoms object. Here we use GPAW and set a few calculation parameters as well:

from gpaw import GPAW

calc = GPAW(mode='lcao', basis='dzp', txt='gpaw.txt', xc='LDA')
atoms.calc = calc

Different electronic structure codes have different input parameters. GPAW can use real-space grids (mode='fd'), planewaves (mode='pw'), or localized atomic orbitals (mode='lcao') to represent the wavefunctions. Here we have asked for the faster but less accurate LCAO mode, together with the standard double-zeta polarized basis set ('dzp'). GPAW and many other codes require a unit cell (or simulation box) as well. Hence we center the atoms within a box, leaving 3 Å of empty space around each atom:

atoms.center(vacuum=3.0)
print(atoms)

The printout will show the simulation box (or cell) coordinates, and the box can also be seen in the GUI.

Once the Atoms have a calculator with appropriate parameters, we can do things like calculating energies and forces:

e = atoms.get_potential_energy()
print('Energy', e)
f = atoms.get_forces()
print('Forces')
print(f)

This will give us the energy in eV and the forces in eV/Å. (ASE also provides atoms.get_kinetic_energy(), referring to the kinetic energy of the nuclei if they are moving. In DFT calculations, we normally just want the Kohn–Sham ground state energy which is the “potential” energy as provided by the calculator.)

Calling get_potential_energy() or get_forces() triggers a selfconsistent calculation and gives us a lot of output text. Inspect the gpaw.txt file. You can review the text file to see what computational parameters were chosen. Note how the get_forces() call did not actually trigger a new calculation — the forces were evaluated from the ground state that was already calculated, so we only ran one calculation.

Binding curve

The strong point of ASE is that things are scriptable. atoms.positions is a numpy array with the atomic positions:

print(atoms.positions)

We can move the atoms by adding or assigning other values into some of the array elements. Then we can trigger a new calculation by calling atoms.get_potential_energy() or atoms.get_forces() again.

Exercise

Move one of the atoms so you trigger two calculations in one script.

This way we can implement any series of calculations. When running multiple calculations, we often want to write them into a file. We can use the standard trajectory format to write multiple calculations (atoms and energy) like this:

from ase.io.trajectory import Trajectory
traj = Trajectory('mytrajectory.traj', 'w')
...
traj.write(atoms)

Exercise

Write a loop, displacing one of the atoms in small steps to trace out a binding energy curve \(E(d)\) around the equilibrium distance. Save each step to a trajectory and visualize. What is the equilibrium distance?

Be careful that the atoms don’t move too close to the edge of the simulation box (or the electrons will squeeze against the box, increasing energy and/or crashing the calculation).

Note

The binding will be much too strong because our calculations are spin-paired, and the atoms would polarise as they move apart.

In case we forgot to write the trajectory, we can also run ASE GUI on the gpaw.txt file although its printed precision is limited.

Although the GUI will plot the energy curve for us, publication quality plots usually require some manual tinkering. ASE provides two functions to read trajectories or other files:

  • ase.io.read() reads and returns the last image, or possibly a list of images if the index keyword is also specified.

  • ase.io.iread() reads multiple images, one at a time.

Use ase.io.iread() to read the images back in, e.g.:

for atoms in iread('mytrajectory.traj'):
    print(atoms)

Exercise

Plot the binding curve (energy as a function of distance) with matplotlib. You will need to collect the energies and the distances when looping over the trajectory. The atoms already have the energy. Hence, calling atoms.get_potential_energy() will simply retrieve the energy without calculating anything.

Optional exercise

To get a more correct binding energy, set up an isolated N atom and calculate its energy. Then calculate the molecular atomisation energy \(E_{\mathrm{atomisation}} = E[\mathrm N_2] - 2 E[\mathrm N]\) of the N2 molecule.

You can use atoms.set_initial_magnetic_moments([3]) before triggering the calculation to tell GPAW that your atom is spin polarized.

Solutions

from ase import Atoms
from ase.io import Trajectory
from gpaw import GPAW

atoms = Atoms('N2', positions=[[0, 0, -1], [0, 0, 1]])
atoms.center(vacuum=3.0)

calc = GPAW(mode='lcao', basis='dzp', txt='gpaw.txt')
atoms.calc = calc

traj = Trajectory('binding_curve.traj', 'w')

step = 0.05
nsteps = int(3 / step)

for i in range(nsteps):
    d = 0.5 + i * step
    atoms.positions[1, 2] = atoms.positions[0, 2] + d
    atoms.center(vacuum=3.0)
    e = atoms.get_potential_energy()
    f = atoms.get_forces()
    print('distance, energy', d, e)
    print('force', f)
    traj.write(atoms)

import matplotlib.pyplot as plt
from ase.io import iread

energies = []
distances = []

for atoms in iread('binding_curve.traj'):
    energies.append(atoms.get_potential_energy())
    distances.append(atoms.positions[1, 2] - atoms.positions[0, 2])

ax = plt.gca()
ax.plot(distances, energies)
ax.set_xlabel('Distance [Å]')
ax.set_ylabel('Total energy [eV]')
plt.show()

from ase import Atoms
from gpaw import GPAW

atoms = Atoms('N')
atoms.center(vacuum=3.0)
atoms.set_initial_magnetic_moments([3])

calc = GPAW(mode='lcao', basis='dzp')
atoms.calc = calc
atoms.get_potential_energy()