index | modules | gitlab | page source

# Structure optimization: H2O

Let’s calculate the structure of the H2O molecule.

Exercise

Create an `Atoms` object representing an H2O 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 H2O.

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 H2O?

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

Documentation v: 3.19.0

Versions
  3.20.0b1 (development)
  3.19.0 (latest stable)