Structure optimization: H₂O

Let’s calculate the structure of the H₂O molecule.

Exercise

Create an Atoms object representing an H₂O molecule by providing chemical symbols and a guess for the positions. Visualize it, making sure the molecule is V shaped.

Exercise

Run a self-consistent calculation of the approximate H2O molecule using GPAW.

Optimizers

We will next want to optimize the geometry. ASE provides several optimization algorithms that can run on top of Atoms equipped with a calculator:

from ase.optimize import BFGS
opt = BFGS(atoms, trajectory='opt.traj', logfile='opt.log')
opt.run(fmax=0.05)

Exercise

Run a structure optimization, thus calculating the equilibrium geometry of H₂O.

The trajectory keyword above ensures that the trajectory of intermediate geometries is written to opt.traj.

Exercise

Visualize the output trajectory and play it as an animation. Use the mouse to drag a box around and select the three atoms — this will display the angles between them. What is H–O–H angle of H₂O?

As always in ASE, we can do things programmatically, too, if we know the right incantations:

from ase.io import read
atoms = read('opt.traj')
print(atoms.get_angle(0, 1, 2))
print(atoms.get_angle(2, 0, 1))
print(atoms.get_angle(1, 2, 0))

The documentation on the Atoms object provides a long list of methods.

G2 molecule dataset

ASE knows many common molecules, so we did not really need to type in all the molecular coordinates ourselves. As luck would have it, the ase.build.molecule() function does exactly what we need:

from ase.build import molecule
atoms = molecule('H2O', vacuum=3.0)

This function returns a molecule from the G2 test set, which is nice if we remember the exact name of that molecule, in this case \('H2O'\). In case we don’t have all the molecule names memorized, we can work with the G2 test set using the more general ase.collections.g2 module:

from ase.collections import g2

print(g2.names)  # These are the molecule names
atoms = g2['CH3CH2OH']
view(atoms)
view(g2)  # View all 162 systems

Use another calculator

We could equally well substitute another calculator, often accessed through imports like from ase.calculators.emt import EMT or from ase.calculators.aims import Aims. For a list, see ase.calculators or run:

$ ase info --calculators

For these tutorials we also have FHI-Aims installed. Let’s run the same relaxation with FHI-Aims then. But in the list above, Aims (probably) wasn’t listed as available. We first need to tell ASE how to run Aims – a typical configuration step for many ASE calculators. This means specifying 1) the command used to run Aims, and 2) where to find information about chemical species. We can do this by setting environment variables in the shell:

$ export ASE_AIMS_COMMAND=aims.x
$ export AIMS_SPECIES_DIR=/home/alumne/software/FHIaims/species_defaults/light

Now ase info --calculators should tell us that it thinks Aims is installed as aims.x.

However, if we open a new shell it will forget this. And we don’t want to modify .bashrc on these computers. Let’s instead set these variables in our Python script:

import os
os.environ['ASE_AIMS_COMMAND'] = 'aims.x'
os.environ['AIMS_SPECIES_DIR'] = '/home/alumne/software/FHIaims/species_defaults/light'

Exercise

Run a structure optimization of H₂O using the FHI-Aims calculator.

To enable the calculation of forces, you will need compute_forces=True. Aims will want an explicitly given XC functional, so we put xc='LDA'. The xc keyword is supported by several ASE calculators to make it easier to specify common XC functionals.

After running the calculation, some new files will be present. ASE has generated control.in and geometry.in, then ran FHI-aims on them, producing aims.out. Be sure to briefly inspect the files. Being perfectionist and/or paranoid, we of course want to be sure that the ASE interface set the parameters the way we wanted them.

Most ASE calculators can be made to generate a file without triggering a calculation using calc.write_input_file(atoms). This is useful, say, if you want to generate the files now but run them later, with or without ASE.

ASE knows many file formats. ase.io.read() can read both the input file and the output file, returning Atoms. These files can also be opened directly with the ASE GUI.

Note that by default, subsequent calculations will overwrite each other. Hence the Aims input and output files correspond to the final step of the structure relaxation. The documentation on ase.optimize will tell us that we can override this behaviour by adding an observer, or using the even more flexible ase.optimize.Dynamics.irun() method to force different steps into different directories.

Appendix: Communication between calculators and codes

What follows is not necessary knowledge for normal usage of ASE. Unless you are interested in how to optimize the communication between ASE and external calculators you may skip ahead.

Different calculators communicate with computational codes in different ways. GPAW is written in Python, so ASE and GPAW run within the same process. However FHI-aims is a separate program. What the Aims calculator does for us is to generate an input file, run FHI-aims, read the output, and return the results.

We just ran a relaxation which involved multiple geometry steps. Each step, a new Aims process is started and later stopped. This is inefficient because the ground-state density and wavefunctions of one step would be an excellent initial guess for the next step, lowering the number of steps necessary to converge. But these quantities are lost when the program terminates. To get the best performance in structure optimisations and dynamics, we need to avoid this loss of efficiency.

Many ASE calculators support more advanced ways of communicating. These calculators can communicate with persistent external processes over pipes (Lammpsrun, CP2K) or sockets (Siesta, Aims, Espresso), or they can work within the same process through direct library calls (Lammpslib, GPAW).

ASE can communicate with FHI-aims over sockets using the i-PI protocol (http://ipi-code.org/). This is done by wrapping the calculator in a ase.calculators.socketio.SocketIOCalculator. The socket calculator will use the calculator it wraps to launch a calculation, then run it.

The documentation on the socket I/O calculator already provides full examples, so we only need minor adjustments to run them on our local machine.

Optional exercise

Based on our previous relaxation with FHI-aims, write a script which runs the same calculation using the ase.calculators.socketio.SocketIOCalculator.

You can run time python3 myscript.py to see how long time the calculation takes in total. How much of a speedup do you get from running the relaxation over a socket? INET sockets often have high latency. If you don’t see much of a speedup, this is probably why. In that case, try switching to a UNIX socket.

The socket I/O calculator automatically generated an input file and also immediately launched the calculation. Since it only launches the process once, subsequent steps don’t overwrite each other and we can find all the intermediate steps in aims.out.

Solutions

GPAW optimisation:

from ase import Atoms
from ase.optimize import BFGS
from gpaw import GPAW

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

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

opt = BFGS(atoms, trajectory='opt.traj')
opt.run(fmax=0.05)

FHI-aims optimisation:

import os
from ase import Atoms
from ase.optimize import BFGS
from ase.calculators.aims import Aims

os.environ['ASE_AIMS_COMMAND'] = 'aims.x'
os.environ['AIMS_SPECIES_DIR'] = '/home/alumne/software/FHIaims/species_defaults/light'

atoms = Atoms('HOH',
              positions=[[0, 0, -1], [0, 1, 0], [0, 0, 1]])

calc = Aims(xc='LDA', compute_forces=True)
atoms.calc = calc

opt = BFGS(atoms, trajectory='opt-aims.traj')
opt.run(fmax=0.05)

FHI-aims/socket-io optimisation:

import os
from ase import Atoms
from ase.optimize import BFGS
from ase.calculators.aims import Aims
from ase.calculators.socketio import SocketIOCalculator


os.environ['ASE_AIMS_COMMAND'] = 'aims.x'
os.environ['AIMS_SPECIES_DIR'] = '/home/alumne/software/FHIaims/species_defaults/light'

atoms = Atoms('HOH',
              positions=[[0, 0, -1], [0, 1, 0], [0, 0, 1]])
opt = BFGS(atoms, trajectory='opt-aims-socketio.traj')

aims = Aims(xc='LDA',
            compute_forces=True,
            use_pimd_wrapper=('UNIX:mysocket', 31415))

with SocketIOCalculator(aims, unixsocket='mysocket') as calc:
    atoms.calc = calc
    opt.run(fmax=0.05)