Skip to content

FHI-aims with Pymatgen

Pymatgen version

This tutorial is known to work when using Pymatgen version 2024.09.10. Pymatgen version is stored in the pymatgen.core.__version__ attribute, so you can use this code snippet to check it:

>>> import pymatgen.core
>>> pymatgen.core.__version__
'2024.9.10'

This tutorial is devoted to getting acquainted with the Pymatgen interface to FHI-aims. It is loosely based on the Jupyter Notebook in the Pymatgen examples Git repository. The tutorial consists of four parts:

  • Pymatgen configuration;
  • creating FHI-aims inputs;
  • using FHI-aims input sets;
  • analyzing FHI-aims outputs.

Configuration

Before any work can be done with Pymatgen, it has to be configured.
Pymatgen configuration is by default stored in YAML format in ~/.config/.pmgrc.yaml file. The configuration can be done in two ways: either manually editing the file, or use pmg command line tool to insert the necessary config lines. The variable that needs to be set is AIMS_SPECIES_DIR, and it should point to the root species_defaults directory inside FHI-aims folder: 

AIMS_SPECIES_DIR: ~/FHIaims/species_defaults

The path should be either absolute or using tilde as home directory expansion; bash variables like ${HOME} do not work here. To use pmg to add this line to the config, run in the activated environment: 

(atomate2-env) user@host:~$ pmg config -a AIMS_SPECIES_DIR ~/FHIaims/species_defaults

This command will back the existing Pymatgen config up and automatically expand the home directory.

Pymatgen configuration file

Please note that the given location of the Pymatgen configuration file is default for all Pymatgen installations in all existing Python environments under a given user. Normally this should not be an issue; if, however, you would like to change that, please refer to the Pymatgen documentation for the guidance.

Creating FHI-aims inputs

All Aims-related Pymatgen objects reside in pymatgen.io.aims package. To create FHI-aims inputs, we'll use AimsControlIn and AimsGeometryIn objects from pymatgen.aims.io module. To use the objects, we import them to our namespace (importing also the necessary auxiliary packages): 

import numpy as np
from pymatgen.core import Structure
from pymatgen.io.aims.inputs import AimsControlIn, AimsGeometryIn

Let's first create the atomic structure object in Pymatgen. Such objects are described using Structure or Molecule class instances, depending on the periodicity of the structure. To create a Structure, one needs to know at least lattice vectors, atomic species and coordinates. Also, site-specific properties, such as initial charges (charge), spin moments (magmom) and velocities (velocity), can be added to the structure as the keys to site_properties dictionary. The next line creates the Structure object for Si primitive cell with the initial charge of 1 \(e^-\) on the first atom: 

structure = Structure(
    lattice=np.array([[0, 2.715, 2.715], [2.715, 0, 2.715], [2.715, 2.715, 0]]),
    species=["Si", "Si"],
    coords=np.array([[0., 0., 0.], [0.25, 0.25, 0.25]]),
    site_properties={"charge": [1, 0]}
)

We can create an AimsGeometryIn object from this structure: 

geo_in = AimsGeometryIn.from_structure(structure)
and check whether the object has been generated successfully by examining its content: 
>>> print(geo_in.content)
lattice_vector  0.000000000000e+00  2.715000000000e+00  2.715000000000e+00
lattice_vector  2.715000000000e+00  0.000000000000e+00  2.715000000000e+00
lattice_vector  2.715000000000e+00  2.715000000000e+00  0.000000000000e+00
atom  0.000000000000e+00  0.000000000000e+00  0.000000000000e+00 Si
     initial_charge 1.000000000000e+00
atom  1.357500000000e+00  1.357500000000e+00  1.357500000000e+00 Si

After having created the geometry, we will create the control.in file containing calculation parameters. For that we will instantiate AimsControlIn object with a dictionary of control keywords as input: 

cont_in = AimsControlIn(
    {
        "xc": "pw-lda",
        "relax_geometry": "trm 0.01",
        "relax_unit_cell": "full",
        "k_grid": [8, 8, 8],
        "species_dir": "light"
    }
)
This will create an object that contains FHI-aims keywords describing relaxation run. species_dir keyword can be set either to be a string (a path to the species defaults relative to ${AIMS_SPECIES_DIR} or ${AIMS_SPECIES_DIR}/defaults_2020), or a dictionary mapping species label to such a string.

control.in keywords

As of now, this dictionary can contain any key-value pairs and is not checked against any kind of schema. A special care might be necessary from the user's side to confirm that the input dictionary really leads to a meaningful control.in file.

AimsControlIn object behaves like a dictionary; that means that one can add the keywords to it and alter the existing keywords by addressing to the different keys: 

cont_in["k_grid"] = [4, 4, 4]

To facilitate the creation of output cube blocks, we recommend using AimsCube objects, which translate the entities defining the 3D space in cube files (like origin, edges, points) to the properly formatted control.in block. One can assign a list of AimsCube objects to cont_in["cubes"] key, and this will result in a several output cube blocks being appended to control.in.

Other output keywords

As you may know, FHI-aims supports multiple strings starting with output keyword. As of now, this is supported by Pymatgen if one assigns multiple strings to cont_in["output], one string at a time. So, the following Python snippet: 

cont_in["output"] = "hirshfeld"
cont_in["output"] = ["eigenvectors"]
will result in the following control.in lines: 
output                             hirshfeld
output                             eigenvectors
We understand that this is a questionable design choice and will be changing it in the future. When it happens, this tutorial will be changed as well.

Both created objects have a write_file method, which writes a corresponding file into a given directory. One can specify if overwriting is permitted; also AimsControlIn.write_file needs information about the structure to append species defaults. The following code snippet will make a directory named workdir and write both files to it: 

work_dir = Path.cwd() / "workdir"
work_dir.mkdir(exist_ok=True)

geo_in.write_file(work_dir, overwrite=True)
cont_in.write_file(structure, work_dir, overwrite=True)

The complete API for the described objects one can find in Pymatgen documentation.

Working with Pymatgen input sets for FHI-aims

Input sets are the Pymatgen take on the code-agnostic calculations. The basic concept behind an input set is to specify a scheme to get a consistent set of code inputs from a structure without further user intervention. Input sets translate the inputs in general language (like timestep in MD simulation) into specific code input keywords.

Pymatgen has several input sets for different types of FHI-aims calculations, which are located in pymatgen.io.aims.sets:

Let's generate the input set of the relaxation calculation. To do that, we first import the InputSetGenerator class from pymatgen.io.aims.sets.core and instantiate it with the default values. The only parameter we need to provide is the species defaults' directory: 

from pymatgen.io.aims.sets.core import RelaxSetGenerator
set_generator = RelaxSetGenerator(user_params={"species_dir": "light"})

After that, we can generate the input set and write it to the provided directory. The input set needs structure at generating time: 

input_set = set_generator.get_input_set(structure)
input_set.write_input(work_dir)

This populates work_dir with three files: control.in, geometry.in and parameters.json, where all the parameters used for the input set are written out in JSON format. The control.in file has the following content: 

xc                                 pbe
relativistic                       atomic_zora scalar
relax_geometry                     trm 1.000000e-03
relax_unit_cell                    full
k_grid                             12 12 12

These are the defaults for this type of calculation. These defaults can be altered in two ways: the keywords pertaining to the calculation type (like relaxation method, force threshold or a relax_unit_cell flag) can be set as keywords to RelaxSetGenerator at the instantiating time (method, max_force and relax_cell bool flag, respectively). Other keywords can be set as the key-value pairs of the user_params dictionary, in the same way as we have set species_dir.

Default k-grid

If k_grid keyword is not set, then it is chosen automatically, based on the k-grid density of 20Å.

The keywords for other input sets can also be found in Pymatgen documentation.

Parsing FHI-aims outputs

Pymatgen also includes the parser for FHI-aims standard output file, located in AimsOutput class. To parse FHI-aims standard output, we import AimsOutput class and instantiate it with the path to the FHI-aims output file.

Let's run the FHI-aims calculation in the work_dir. If you're using Jupyter Notebook to follow the tutorial, you can add exclamation mark in front of the command to run it in the shell. Assuming you have aims.x in your PATH, run the following command in the notebook cell: 

!cd workdir; mpirun -np 4 aims.x | tee aims.out
This changes the directory to workdir (the value of work_dir variable) and runs FHI-aims executable on 4 cores, copying the standard output also to aims.out file. After the calculation is finished, we can create AimsOutput class instance from the resulting file: 
output = AimsOutput.from_outfile(work_dir / 'aims.out')

After that, the calculation results are available as the attributes of the object: 

>>> output.final_energy
-15802.6578366142

>>> output.stress
array([[-0.02356802,  0.        ,  0.        ],
       [ 0.        , -0.02355199, -0.        ],
       [ 0.        , -0.        , -0.02355199]])

>>> output.structure_summary
{'initial_structure': Structure Summary
 Lattice
     abc : 3.839589821842953 3.839589821842953 3.839589821842953
  angles : 60.00000000000001 60.00000000000001 60.00000000000001
  volume : 40.02575174999999
       A : np.float64(0.0) np.float64(2.715) np.float64(2.715)
       B : np.float64(2.715) np.float64(0.0) np.float64(2.715)
       C : np.float64(2.715) np.float64(2.715) np.float64(0.0)
     pbc : True True True
 PeriodicSite: Si (0.0, 0.0, 0.0) [0.0, 0.0, 0.0]
 PeriodicSite: Si (1.357, 1.357, 1.357) [0.25, 0.25, 0.25],
 'initial_lattice': Lattice
     abc : 3.839589821842953 3.839589821842953 3.839589821842953
  angles : 60.00000000000001 60.00000000000001 60.00000000000001
  volume : 40.02575174999999
       A : np.float64(0.0) np.float64(2.715) np.float64(2.715)
       B : np.float64(2.715) np.float64(0.0) np.float64(2.715)
       C : np.float64(2.715) np.float64(2.715) np.float64(0.0)
     pbc : True True True,
 'is_relaxation': True,
 'is_md': False,
 'n_atoms': 2,
 'n_bands': 20,
 'n_electrons': 28,
 'n_spins': 1,
 'electronic_temperature': 0.0,
 'n_k_points': 868,
 'k_points': None,
 'k_point_weights': None}