Structure optimization: H2O¶
Let’s calculate the structure of the H2O molecule.
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.
Run a self-consistent calculation of the approximate H2O molecule using GPAW.
We will next want to optimize the geometry.
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)
Run a structure optimization, thus calculating the equilibrium geometry of H2O.
trajectory keyword above ensures that the trajectory of intermediate
geometries is written to
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
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
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
ase info --calculators should tell us that it thinks
Aims is installed as
However, if we open a new shell it will forget this. And we don’t want to
.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'
Run a structure optimization of H2O
using the FHI-
To enable the calculation of forces, you will need
Aims will want an explicitly given XC functional, so we put
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
ran FHI-aims on them, producing
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
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
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
will tell us that we can override this behaviour by adding an observer,
or using the even more flexible
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 (
CP2K) or sockets (
or they can work within the same process
through direct library calls
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.
Based on our previous relaxation with FHI-aims, write a script
which runs the same calculation using the
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
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)
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/myname/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)
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/myname/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)