Skip to main content
`Ctrl`+`K`

# 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 gpaw import GPAW

    from ase import Atoms
    from ase.io import Trajectory

    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 gpaw import GPAW

    from ase import Atoms

    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()