E.4 cffi — Python 2/3 interface to FHI-aims

The extension can be linked to aims using both Make. Note that when you compile with Python 2 or 3, you then need to use that particular version in your user scripts.

• •

  To compile with Make, define USE_CFFI=yes. This also requires specifying a Python interpreter with PYTHON=<path to python>. The Python include files and libraries are detected automatically and printed at the top of the Make output.

The interface is written using the iso_c_binding module, which is a Fortran 2003 feature. It is supported by all major compilers.

The build extension links FHI-aims to a dynamic Python library. Under ideal circumstances, everything should work out without any intervention, but manual intervention may be needed in non-standard environments. The Python installation needs to have packages cffi ($\ge$1.5.2), Numpy/Scipy. The mpi4py package is not required, but its absence severely limits the functionality.

In Linux, dynamic libraries are usually recorded only by their names during linking and found dynamically at runtime at several standard locations. If you link against your own Python installation, such as Anaconda, you need to make sure that FHI-aims can find the Python library at runtime. This can be done by adding anaconda/lib to $LD_LIBRARY_PATH. Since Anaconda now packages includes its own version of the Intel Math Kernel Library, which could clash with any system MKL, it is recommended to uninstall MKL from Anaconda if you include it in $LD_LIBRARY_PATH.¹¹ 1 A more user-friendly behaviour could be setup using RPATH in the future.

On OS X, dynamic libraries are supposed to contain their absolute installation path, in which case these paths are recorded during linking and used at runtime. This is the case with the system Python and Homebrew Python, but not with Anaconda Python. To use this extension with Anaconda Python on OS X, either set the variable $DYLD_LIBRARY_PATH to the Anaconda library directory anaconda/lib or (better) change the install path of the Python library in the compiled FHI-aims binary with the system install_name_tool. Furthermore, some version of Anaconda Python on OS X do not seem to properly initialize Python paths when embedded, in which case $PYTHONHOME needs to be set when running FHI-aims.

The Python interface can be used either in an interactive or non-interactive mode. The interactive mode is activated by the option REPL (see python_hook). In this mode, AIMS launches an interactive Python console at a given location. This is available only in a serial run. In the console, connection to the running FHI-aims instance is provided via a local variable ctx. Various quantities are accessible as attributes of this context object. For example,

  Self-consistency cycle converged.
  [...]
------------------------------------------------------------
  Executing Python hook: post_scf

There is a local variable ‘ctx‘.
See ‘help(ctx)‘ for details.
Press CTRL+D to
continue the aims run.
Type ‘exit(1)‘ to abort aims.
+>>> ctx
<AimsContext ’rho, batches, coords, partition_tab’>
+>>> ctx.coords
array([[-3.77945226,  3.77945226],
       [ 0.,  0.],
       [ 0.,  0.]])
+>>> ctx.rho
array([[  5.81201802e-02,   3.72428091e-01,   3.72426419e-01, ...,
          1.54065044e-05,   4.86078607e-12,   1.17640761e-29]])
+>>> ctx.rho[:] = 0
+>>>

The context object provides also several convenience functions which are documented in help(ctx). When the console is quit normally (CTRL+D), the FHI-aims run continues. When the return code is non-zero (for example with exit(1)), the FHI-aims run is aborted.

In the non-interactive mode, which is available also in parallel runs, the context object is passed to the user-defined function run which is loaded from the specified file. Two locations are currently available at which the run function is executed, as specified by the argument to the python_hook keyword:

post_scf

immediately after the end of the self-consistent loop, before any post-processing

post_hirshfeld

immediately after the end of the Hirshfeld analysis, that is, before any van der Waals routines

A minimal user script (can be both Python 2 and 3) could look for example like this:

import json

def run(ctx):
    with open(’coords.json’, ’w’) as f:
        json.dump(ctx.coords.tolist(), f)

Note that the provided variables are local to a given MPI process. For synchronization, one can use the mpi4py Python package. Some convenience synchronisers are defined in ctx. For exampeles of how to use the mpi4py package, see the commented implementation of ctx.gather_all_grids() in cffi/python_interface.py. The user script can optionally define also function parse, which is then called without any arguments during parsing of control.in. This can be used to verify or precompute certain data before launching the full calculation. Any uncaught exceptions in the both parse() and run() cause aims to abort immediately. To recapitulate, a user script has one or two entry points: (i) the mandatory run is executed at the location specified in control.in and takes the single context argument and (ii) the optional parse function, which does not take any arguments, and is always called at the end of parsing control.in.

At the time of writing, two example hooks are provided in aimsfiles/utilities/python_hooks/. The first one plots the electron density in the $xy$ plane. Below is a shortened glimpse to get a general idea, see the original file for the complete implementation.

# ~~ imports ~~

bohr = 0.52917721067

def run(ctx):
    points, rho = ctx.gather_all_grids([’rho’])
    if ctx.rank == 0:  # the grids are gathered only on the root process
        points = points.T  # (3, npts) -> (npts, 3)
        rho = rho.sum(0)  # sum over spin
        in_plane = abs(points[:, 2]) < 1e-10  # points in xy plane
        points = bohr*points[in_plane, 0:2]  # filter & take x, y & scale
        rho = np.log10(1+rho[in_plane])  # filter & log scale
        X, Y = np.mgrid[-4:4:400j, -2:2:200j]  # get rectangular grid
        # interpolate density to rectangular grid
        rho = griddata(points, rho, (X, Y), method=’cubic’)
        # ~~ matplotlib plotting ~~

Saving this file into plot-density.py and putting python_hook plot-density.py parallel into control.in produces Figure E.1 for an argon dimer.

The other example hook loads Hirshfeld volumes from the output of a previous FHI-aims calculation. The parsing part of the hook is executed when FHI-aims parses control.in and searches for the Hirshfeld analysis section in hirshfeld.out. The actual modification in all processes happens after the normal Hirshfeld analysis is performed by specifying python_hook load-hirshfeld.py parallel.

import numpy as np

volumes = None

def parse():
    global volumes  # otherwise volumes would be local to function
    with open(’hirshfeld.out’) as f:
        while ’Performing Hirshfeld analysis’ not in next(f):
            pass  # find the Hirshfeld section
        volumes = np.array(list(get_volumes(f)))

def get_volumes(f):
    for line in f:
        if not line.strip():
            return
        if ’Free atom volume’ in line:
            free = float(line.split()[-1])
        elif ’Hirshfeld volume’ in line:
            hirsh = float(line.split()[-1])
            yield hirsh/free

def run(ctx):
    ctx.hirshfeld_volume[:] = volumes  # replace values in-place

The two-step operation saves computation time if the parsing ended with an error, since it would abort FHI-aims right after start. By specifying also

sc_iter_limit             0
postprocess_anyway        .true.

in control.in, one can skip the SCF cycle completely and evaluate only the van der Waals part of the energy.

To extend the range of objects provided by the context object, several simple steps need to be followed.

python_interface.f90

The Fortran type AimsContext_t needs to be extended and the added component needs to be initialized in get_aims_context() and potentially deallocated in destroy_aims_context() if it required any allocation. The latter is of concern only for mapping Fortran types, not simple arrays.

cffi/python_interface.h

The corresponding C struct AimsContext_t needs to be updated accordingly. When these two steps are done, the object is already available in its raw form in ctx._c_ctx.

cffi/python_interface.py

The initialization of the Python object AimsContext in call_python_cffi_inner() needs to be updated accordingly. For simple arrays, this amounts to simply wrapping the raw C pointers with Numpy arrays. For wrapping Fortran types, see how BatchOfPoints and GridPoints are wrapped.

To add new hook labels, add the following lines to the desired location:

use python_interface, only: run_python_hook, python_hooks

% ...

if (python_hooks%<hookname>%registered) then
    call run_python_hook(python_hooks%<hookname>)
end if

and extend HookRegister_t and register_python_hook() in python_interface.f90 by adding <hookname>.