Tutorial 2: ASE ONETEP interface

Author:: Tom Demeyere
Date:: August 2023, updated June 2024

This tutorial will guide you through the use of the ASE interface to ONETEP. The Atomic Simulation Environment (ASE) is a set of tools and Python modules for setting up, manipulating, running, visualizing and analyzing atomistic simulations. The ASE interface to ONETEP allows you to set up and run ONETEP calculations from Python scripts.

Some tutorials (1, 4, 10) have Jupyter notebooks that you can use to run examples interactively. You can download them below, along with the needed pseudopotentials, input files and ONETEP launching script. Alternatively, you can clone the entire repository, and have them almost ready to run (you will still need to tweak the onetep_launch.sh for your machine.)

To run these tutorials, please have a look inside the launch_onetep.sh script you downloaded in order to load correctly all modules and environment variables needed on your machine. Similarly, you might need to change the paths of the OnetepProfile object inside the Jupyter Notebooks to match your paths. The OnetepProfile is an ASE object that must be created and passed to the ASE ONETEP calculator in order to specify both ONETEP command and pseudopotential path. More information is available just below.

This tutorial will mainly focus on running ONETEP calculations with ASE, if you are not familiar with ASE as a whole, feel free to consult the mini-tutorial here:

learn_ase_in_y_minutes.py

ASE Configuration File

To run ONETEP with ASE, you need to have a configuration file that specifies the path to the ONETEP binary and the location of the pseudopotentials. There is no default configuration file, so you need to create one.

ASE looks for a configuration file named config.ini in the default location $HOME/.config/ase/. The configuration file should follow the pattern (do not put quotes around the values):

[onetep]
command = mpirun -np 10 -v /path/to/onetep/binary
pseudo_path = /path/to/pseudos

Replace /path/to/onetep/binary with the actual path to your ONETEP binary. If you want to use a different location for the configuration file, you can set the ASE_CONFIG_PATH environment variable. For example:

$ export ASE_CONFIG_PATH="/path/to/custom/config.ini"

Alternatively, if you don’t want to use the configuration file, you can create a OnetepProfile directly in your script. For example:

from ase.calculators.onetep import Onetep, OnetepProfile

profile = OnetepProfile(
  command="mpirun -np 10 -v /path/to/onetep/binary",
  pseudo_path="/path/to/pseudos"
)
calc = Onetep(profile=profile)

This will override the configuration file and use the OnetepProfile object instead.

Pseudopotentials

If no pseudopotentials are passed ASE will try to guess the files based on the element used and the pseudo_path variable. Otherwise you can pass a dictionary with the element as key and the path to the pseudopotential as value. The path can be either absolute or relative to the pseudo_path.

# Explicitly providing each path:
calc = Onetep(pseudopotentials = {'H': '/path/to/pseudos/H.usp', 'O': '/path/to/pseudos/O.usp'})
# Using relative paths:
calc = Onetep(pseudopotentials = {'H': 'H.usp', 'O': 'O.usp'})
# ASE will try to guess them if you don't provide them:
calc = Onetep()

For ASE to correctly guess the pseudopotentials, it is best to use a pseudo_path that contains only one pseudopotential file for each element. If there are multiple files for the same element, ASE will not be able to guess which one to use.

ONETEP Calculator

Simple calculations can be setup calling the Onetep calculator without any parameters, in this case ONETEP’s default parameters will be used. For more complex cases using the keywords parameters is necessary. The keywords parameters is a dictionary, in which each of the keys is a string that should be a ONETEP keyword, and the corresponding value is what you want to set that keyword to in the input.

from ase.calculators.onetep import Onetep

# Default parameters
calc = Onetep()

# Custom parameters
keywords = {
    'xc' : 'PBE',
    'do_properties' : True,
    'cutoff_energy' : 35,
    'output_detail': 'verbose',
    'elec_energy_tol': 1.0e-5,
}

calc = Onetep(keywords=keywords)

Alternatively you can read an already existing input file with the function read_onetep_keywords

from ase.io.onetep import read_onetep_keywords

keywords = read_onetep_keywords('input_file.dat')

# Let's change one specific keyword
keywords['xc'] = 'PBE0'

calc = Onetep(keywords=keywords)

Examples

Here is an example python script which sets up a calculation on a water molecule:

from ase.build import molecule
from ase.calculators.onetep import Onetep

water = molecule('H2O', vacuum=10)

calc = Onetep(xc='PBE', paw=True)
water.calc = calc

water.get_potential_energy()

Here is a more complex example, this time setting up a $\mathrm{Pt}_{13}$ cluster and running a geometry optimisation, note that here as far as ONETEP is concerned we are running singlepoint calculations, the geometry optimisation is done by ASE’s BFGS optimiser:

import numpy as np

from ase.build import molecule
from ase.calculators.onetep import Onetep
from ase.cluster import Octahedron
from ase.optimize import BFGSLineSearch

# Pt13 from ase.cluster
nano = Octahedron('Pt', 3, 1)
nano.center(vacuum=10)

# ONETEP default are atomic units, one can specify 'cutoff_energy' : '600 eV' if needed.
keywords = {
    'xc' : 'rpbe',
    'do_properties' : True,
    'cutoff_energy' : 35,
    'output_detail': 'verbose',
    'elec_energy_tol': 1.0e-5/len(atoms),
    'edft': True,
}

# append = True will not overwrite file at each step
calc = Onetep(
    append = True,
    keywords = keywords)

nanoparticle.calc = calc

opt = BFGSLineSearch(atoms = nano)
opt.run(fmax=0.1)

Here is an example of setting up an EELS and LDOS calculation on an N-substituted graphene sheet, demonstrating several more advanced functionalities (tags, species groups, and overrides to pseudopotentials and atomic solver strings)

import numpy as np

from ase.build import graphene_nanoribbon
from ase.calculators.onetep import Onetep
from ase.io import write

sheet = graphene_nanoribbon(10, 10, type='zigzag', vacuum = 10)

# Get all distances to center of mass
com = sheet.get_center_of_mass()
distances_to_com = np.linalg.norm(sheet.positions - com, axis = 1)

# Find atoms close to com and change one randomly to N
p, = np.where(distances_to_com < 5)
to_nitro = choice(p)
sheet[to_nitro].symbol = 'N'

shell_rad = np.array([1.5, 2.5, 3.0, 4.0, 4.5])

tags = np.zeros(len(sheet), dtype=np.int32)

# We want to tag atoms that are close to the introduced nitrogen
for idx, rad in enumerate(reversed(shell_rad)):
    # All distances N-C
    dist = norm(sheet[to_nitro].position - sheet.get_positions(), axis = 1)
    # Which ones are closest to rad?
    p, = np.where(dist < rad)
    # Cannot be the nitrogen itself
    p = p[p != to_nitro]
    # Tags them
    tags[p] = len(shell_rad) - idx

sheet.set_tags(tags)

tags = ['' if i == 0 else i for i in tags]

species = np.unique(np.char.add(sheet.get_chemical_symbols(), tags))

keywords = {
    'species_core_wf' : ['N /path/to/pseudo/corehole.abinit'],
    'species_solver' : ['N SOLVE conf=1s1 2p4'],
    'pseudo_path': '/Users/tomdm/PseudoPotentials/SSSP_1.2.1',
    'xc' : 'PBE',
    'paw': True,
    'do_properties': True,
    'cutoff_energy' : '500 eV',
    'species_ldos_groups': species,
    'task' : 'GeometryOptimization'
}

calc = Onetep(
    keywords = keywords
)

# Checking the input before running the calculation
write('to_check.dat', sheet, format='onetep-in', keywords = keywords)

sheet.calc = calc
# Will actually run the geometry optimisation
# using ONETEP internal BFGS
sheet.get_potential_energy()

Quickly restart with solvation effect using the soft sphere solvation model:

from ase.io import read
from ase.io.onetep import get_onetep_keywords

# Read from the previous run...
optimized_sheet = read("onetep.out")

# Function to retrieve keywords dict from input file...
keywords = get_onetep_keywords('onetep.dat')

# We add solvation keywords
keywords.update(
    {
    'is_implicit_solvent': True,
    'is_include_apolar': True,
    'is_smeared_ion_rep': True,
    'is_dielectric_model': 'fix_cavity',
    'is_dielectric_function' : 'soft_sphere'
    }
)

optimized_sheet.calc = Onetep(keywords=keywords)
optimized_sheet.get_potential_energy()

Important note

If you are not keen about using ASE to run ONETEP calculations, it is always possible to use ASE to write ONETEP input files and run them manually. This should be done by using the general ASE IO modules ase.io.write and ase.io.read to write and read ONETEP input files. In every example above, all you need to do is to replace the get_potential_energy() call by a write call to write the input file, such as write('input_file.dat', atoms, format='onetep-in', keywords=keywords). You can then run the ONETEP binary manually as you always do.

How to use ASE on HPCs

If the HPC you are using has a module system, you can load the conda module and create an environment with the required packages. If you don’t have access to a module system, you can install miniforge in your home directory and create an environment there. A tutorial to do so is available at the end of this document.

How does python launch ONETEP under the hood?

When you run a python script with ASE and ONETEP, ASE will both construst the command to be launched and the input file. The command will be constructed based on the command key in the ASE configuration file. Or based on the command key in the OnetepProfile object if you send the profile manually. The command will be executed with the subprocess module using the check_call function. The inner working of the check_call function is to run the command in a subprocess and wait for it to finish. If the command fails, an exception will be raised. To run the command no new shell is created, and all the environment variables are inherited from the parent process. All stdout and stderr will be redirected to the onetep.out and onetep.err files.

The input file will be created by the IO functions of ASE, namely ase.io.onetep.write_onetep_input. This function will write the input file in the format expected by ONETEP. This will be automatically done if a calculation is launched via atoms.get_potential_energy() or else.

General case

There are two ways to submit job using ASE on HPC, you can directly sbatch the python script by putting the correct shebang at the top of the script, or you can use an additional bash script to submit the job. The bash script will have to activate the environment and run the python ASE script. Here is an example of such a script:

#!/bin/bash
#SBATCH --job-name=ASE_ONETEP
...

conda activate myenv

module load ... # Load all the modules needed by ONETEP
export ... # Set all the environment variables needed by ONETEP

export ASE_CONFIG_PATH="/path/to/scratch/.ase_config.ini"

python my_ase_script.py

# Your python script can look like this

from ase.build import molecule

from ase.calculators.onetep import Onetep

water = molecule('H2O')

keywords = {
    'xc' : 'PBE',
    'do_properties' : True,
    'cutoff_energy' : 35,
    'output_detail': 'verbose',
    'elec_energy_tol': 1.0e-5/len(water),
}

calc = Onetep(keywords=keywords)

water.calc = calc
water.get_potential_energy()

Make sure that the ONETEP command being used contains srun for example: command = srun /path/to/onetep/binary. Otherwise the job will not dispatch correctly on the compute nodes. This is no different from launching a normal job, with the expection that ASE takes care of the input file and the command to be launched.

Archer2

Archer2 is a Cray system, and the conda module is not available. You should install it by having a look at the instruction at the end of this document. One of the Archer2’s particularity to keep in mind is that compute nodes only have access to the scratch space and not to the home directory. You should make sure that every file which will be used during the calculation is accessible from the scratch space, most likely this will be: the input files, the pseudopotentials, the executable and conda. This also means that if you are using the ase config file, you should make sure to change its location with the ASE_CONFIG_PATH environment variable to the scratch space. Once this is done you should have a working environment to run ASE on Archer2.

Iridis5

Iridis5 is an Intel based HPC, with conda available as a module. You can alternavely install your own Conda, following the instruction at the end of this document if you want it. There is no particularity to keep in mind when running ASE on Iridis5, you can use the conda module to create an environment with the required packages. You can then submit a job with the python script directly or with a bash script as shown above. Make sure to use srun in the command to dispatch the job on the compute nodes.

Young

The only particularity of Young is that srun is not available, instead a home-made wrapper around mpirun is made avaible (gerun). This will not cause limitations as long as you keep each job to serial execution. For example, if you use the ASE NEB module with threading, i.e. launching multiple ONETEP in parallel in the same PBS job, gerun will most likely not distribute the job correctly, and the calculation will either fail, or be very slow. The only way around this is to make use of the mpirun command directly and specifying the node to use for each job. Which will not be detailed here, you should probably use another HPC for this kind of calculation.

Other python packages

Other packages that can be used with Onetep + ASE are numerous, here we do mini-tutorials for some of them.

DFTD3/DFTD4

DFTD3 and DFTD4 are dispersion correction methods that can be used with ONETEP. These packages also interface with ASE, which is why they can be used in conjunction with ONETEP. To install DFTD3 or DFTD4, you can use the conda package manager. Here is how to install them:

conda install -c conda-forge dftd3-python
conda install -c conda-forge dftd4-python

If you really care about the performance you should probably compile them yourself, although the performance gain should probably be minimal. After installation they can be used in the ASE calculator as follows:

from ase.build import molecule
from ase.calculators.mixing import SumCalculator
from ase.calculators.onetep import Onetep
from dftd4.ase import DFTD4

atoms = molecule('H2O')

calc = SumCalculator([DFTD4(method="PBE"), Onetep(xc="PBE")])
atoms.calc = calc

atoms.get_potential_energy()

For DFTD3 the code is pretty much the same, just replace DFTD4 by DFTD3. The DFTD3 version requires to have method and damping parameters set at all times. With both versions you can pass an additional parameter params_tweaks where you can manually override the internal D3 parameters, see the documentation for more information.

Alloy Catalysis Automated Toolkit (ACAT)

ACAT (https://gitlab.com/asm-dtu/acat) is a python package that can be used to automate the setup of ONETEP calculations for (alloy) catalysis. ACAT can be used in conjunction with ASE, and can be installed using pip:

pip install acat

The package allows many operations on both surfaces and nanoclusters, the two main classes are the ClusterAdsorptionSites and the SlabAdsorptionSites. Which are used to detect all possible binding sites of your systems. Here is a complete example to create ONETEP input files for an alloyed nanocluster:

from pathlib import Path

from acat.adsorption_sites import ClusterAdsorptionSites
from acat.build.action import add_adsorbate_to_site
from ase.cluster import Octahedron
from ase.io import write

calc_dir = Path("alloy_project_tutorial")
calc_dir.mkdir(exist_ok=True)

atoms = Octahedron("Ni", length=7, cutoff=2)

# Let's create our alloy
for atom in atoms:
    if atom.index % 2 == 0:
        atom.symbol = "Pt"

atoms.center(vacuum=5.0)

# We create the ACAT object with our parameters,
# Many more are available, check the documentation
cas = ClusterAdsorptionSites(
    atoms,
    composition_effect=True,
    label_sites=True,
    surrogate_metal="Ni",
)

# Only unique sites, we don't want to duplicate calculations
sites = cas.get_unique_sites(unique_composition=True)

for site in sites:
    # add_adsorbate_to_site is modifies the object in place
    # so we copy it to avoid modifying the original object
    tmp = atoms.copy()

    add_adsorbate_to_site(tmp, "O", site)

    # We create a unique custom label based on the information
    label = (f"{tmp.get_chemical_formula(mode='metal').lower()}"
            f"_{site['surface']}_{site['site']}_{site['label']}")

    # The directory for this specific calculation
    current_dir = calc_dir / label
    current_dir.mkdir(exist_ok=True)

    # ASE can of course, write onetep input files
    # In practice you would have to specify keywords and pseudopotentials
    write(current_dir / "onetep.dat", tmp, format="onetep-in")

You will have a directory called alloy_project_tutorial with a subdirectory for each adsorption site, each containing an input file for ONETEP. You can then run these input files manually or with ASE as shown in the previous examples. Alternatively you can visualise them using the ase gui tool.

Phonopy

Phonopy (https://github.com/phonopy/phonopy) is a python package that can be used to calculate phonon properties of materials. and can be installed using pip or conda:

pip install phonopy

Phonopy can be used to calculate the phonon band structure of a material. Usually everything is done using the CLI but I personnaly prefer to use the API directly, here is an example for a water molecule:

from ase.build import molecule
from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms

from ase.calculators.onetep import Onetep

water = molecule('H2O', vacuum=10)

calc = Onetep()

phonopy_atoms = PhonopyAtoms(symbols=water.get_chemical_symbols(),
                             positions=water.get_positions(),
                             cell=water.get_cell())

phonopy = Phonopy(phonopy_atoms, supercell_matrix=[[1, 0, 0], [0, 1, 0], [0, 0, 1]])

phonopy.generate_displacements(distance=0.01)

displacements = phonopy.supercells_with_displacements

forces = []

for i, disp in enumerate(displacements):

    disp_dir = Path(f"displacement_{i}")
    disp_dir.mkdir(exist_ok=True)

    atoms = Atoms(disp.get_chemical_symbols(),
                  disp.get_positions(),
                  cell=disp.get_cell()
    )

    calc.directory = str(disp_dir)

    atoms.calc = calc

    forces.append(atoms.get_forces())

phonopy.forces = forces
phonopy.produce_force_constants()

phonon.save("ifc.yaml", settings={'force_constants': True})

print(phonon.get_frequencies_with_eigenvectors((0, 0, 0))[0]*33.356)

With the annoying fact that the Atoms` object has to be manually transfered to PhonopyAtoms back and forth. The phonon frequencies are in THz, to convert them to cm-1 you have to multiply by 33.356. The ifc.yaml file can be used for further processing. See the phonopy documentation for more information.

Many more

There are many more packages that can be used with ONETEP and ASE. Some of them are listed below:

pymatgen: A python package for materials analysis, which can be used to generate structures, calculate band structures, and much more. (https://github.com/materialsproject/pymatgen)
phono3py: A python package for calculating phonon lifetime and thermal conductivity. (https://github.com/phonopy/phono3py)
HiPhive: A python package to compute higher order force constants without using a specific set of configurations. (https://hiphive.materialsmodeling.org/index.html)
Sella: Sella is a utility primarily intended for refining approximate saddle point geometries. Interfaces well with ASE. (https://github.com/zadorlab/sella)
QuAcc: The Quantum Accelerator (QuAcc) is a python package that can be used to create automated workflows and run them concurrently with workflow managers like Parsl, Dask or Covalent. ONETEP has an interface and a few recipes. (https://github.com/Quantum-Accelerators/quacc)

Conda for the Impatient

Why Conda?

Do not pollute your system-wide python, you might regret it: Conda creates isolated environments, keeping your system Python clean and preventing conflicts between different projects.
Stop compiling your tools, use binaries by Conda: Conda can manage packages for various languages, including R, C++, and Fortran, making it a versatile tool for scientific computing.
Complement Conda with pip: While Conda handles most python package installations, you might occasionally need pip for packages not available in Conda repositories.
Conda is self-contained: Install it everywere, no need for root access. Even HPC systems encourage the use of Conda. Conda will not break your system, and you can remove it easily.

Installing Mambaforge on Linux

Download the Mambaforge installer (Linux x86_64) from the Conda Forge repository:

wget https://github.com/conda-forge/miniforge/releases/latest/download/Mambaforge-Linux-x86_64.sh

Run the installer:

bash Mambaforge-Linux-x86_64.sh

Follow the prompts, agreeing to the license and choosing the installation location.
Initialize Mambaforge by running:

conda init

Close and reopen your terminal for the changes to take effect.

Installing Conda on Windows

To install Conda on Windows, follow these steps:

Visit the official Anaconda website (https://www.anaconda.com) and download the Anaconda Navigator.
Run the installer and follow the installation prompts. Make sure to select the option to add Conda to your system’s PATH environment variable.
Once the installation is complete, open the Anaconda Navigator application to manage packages and environments. You can create environments, install packages, and launch Jupyter notebooks directly from the Navigator interface.
If you want to install python packages that are only available through pip you can launch a terminal from the navigator inside the environment you want to install the package and run pip install package_name

Creating and Managing Environments

Create a new environment:: conda create --name myenv
Activate the environment:: conda activate myenv
Deactivate the environment:: conda deactivate

Installing Packages

Install packages in the active environment:: conda install numpy pandas
For packages not available in Conda repositories, use pip:: pip install somepackage

Updating and Removing Packages

Update a package:: conda update somepackage
Remove a package:: conda remove somepackage
Update all packages in the current environment:: conda update --all

Managing Environments

List all environments:: conda env list
Remove an environment:: conda env remove --name myenv
Export an environment to a YAML file:: conda env export > environment.yml
Create an environment from a YAML file:: conda env create -f environment.yml