Using Skala with the Atomic Simulation Environment (ASE)

Using Skala with the Atomic Simulation Environment (ASE)#

This tutorial provides a comprehensive overview of how to use the Skala neural network-based exchange-correlation functional with the Atomic Simulation Environment (ASE). The Skala functional is available as an ASE calculator, enabling accurate and scalable density functional theory calculations on molecular systems.

import numpy as np
from ase.build import molecule
from ase.optimize import LBFGSLineSearch as Opt
from ase.units import Hartree
from skala.ase import Skala

Basic Calculator Setup#

We can create the Skala calculator and set options like the basis set and density fitting when constructing the calculator. The calculator supports several key parameters:

xc: Exchange-correlation functional (default: “skala”)
basis: Basis set for the calculation (e.g., “def2-svp”, “def2-tzvp”)
with_density_fit: Enable density fitting for faster calculations
with_dftd3: Include DFT-D3 dispersion correction (default: True)
verbose: Control output verbosity (0-4)
with_retry: Enable retry mechanism for SCF convergence (default: True)
ks_config: Additional configuration options for the KS solver (e.g., convergence criteria)

# Create a water molecule
atoms = molecule("H2O")

# Set up the Skala calculator with specific parameters
atoms.calc = Skala(xc="skala-1.1", basis="def2-svp", verbose=4)

# Display the calculator parameters
print("Calculator parameters:")
for key, value in atoms.calc.parameters.items():
    print(f"  {key}: {value}")

Calculator parameters:
  xc: skala-1.1
  basis: def2-svp
  with_density_fit: False
  auxbasis: None
  with_newton: False
  with_dftd3: True
  with_retry: True
  charge: None
  multiplicity: None
  verbose: 4
  ks_config: None
  device: cpu

Energy Calculations#

Calculations can be performed by requesting properties from the atoms object, which is connected to the calculator. For example, we can calculate the total energy of a water molecule. The energy is returned in eV (ASE’s default energy unit):

# Calculate the total energy
energy_eV = atoms.get_potential_energy()  # eV
energy_Hartree = energy_eV / Hartree  # Convert to Hartree

print(f"Total energy: {energy_eV:.6f} eV")
print(f"Total energy: {energy_Hartree:.6f} Hartree")

System: uname_result(system='Linux', node='runnervm8k3b6', release='6.17.0-1022-azure', version='#22-Ubuntu SMP Mon Jul 27 17:24:03 UTC 2026', machine='x86_64')  Threads 2
Python 3.12.14 (main, Sep  2 2026, 23:27:36) [GCC 15.3.0]
numpy 2.4.6  scipy 1.18.1  h5py 3.16.0
Date: Thu Oct 15 18:17:17 2026
PySCF version 2.14.0
PySCF path  /home/runner/work/skala/skala/.pixi/envs/docs/lib/python3.12/site-packages/pyscf

[CONFIG] conf_file None
[INPUT] verbose = 4
[INPUT] num. atoms = 3
[INPUT] num. electrons = 10
[INPUT] charge = 0
[INPUT] spin (= nelec alpha-beta = 2S) = 0
[INPUT] symmetry False subgroup None
[INPUT] Mole.unit = Angstrom
[INPUT] Symbol           X                Y                Z      unit          X                Y                Z       unit  Magmom
[INPUT]  1 O      0.000000000000   0.000000000000   0.119262000000 AA    0.000000000000   0.000000000000   0.225372517068 Bohr   0.0
[INPUT]  2 H      0.000000000000   0.763239000000  -0.477047000000 AA    0.000000000000   1.442312677587  -0.901488178545 Bohr   0.0
[INPUT]  3 H      0.000000000000  -0.763239000000  -0.477047000000 AA    0.000000000000  -1.442312677587  -0.901488178545 Bohr   0.0

nuclear repulsion = 9.08829376913928
number of shells = 12
number of NR pGTOs = 38
number of NR cGTOs = 24
basis = def2-svp
ecp = {}
CPU time:         4.40
******** <class 'skala.pyscf.dft.SkalaRKS'> ********
method = SkalaRKS
initial guess = minao
damping factor = 0
level_shift factor = 0
DIIS = <class 'pyscf.scf.diis.CDIIS'>
diis_start_cycle = 1
diis_space = 8
diis_damp = 0
SCF conv_tol = 1e-09
SCF conv_tol_grad = None
SCF max_cycles = 50
direct_scf = True
direct_scf_tol = 1e-13
chkfile to save SCF result = /tmp/tmpuo_wvpp1
max_memory 4000 MB (current use 505 MB)
XC library skala.pyscf.libxc version None
    None
XC functionals = custom
radial grids: 
    Treutler-Ahlrichs [JCP 102, 346 (1995); DOI:10.1063/1.469408] (M4) radial grids
becke partition: Becke, JCP 88, 2547 (1988); DOI:10.1063/1.454033
pruning grids: <function nwchem_prune at 0x7fac6e547920>
grids dens level: 3
symmetrized grids: False
atomic radii adjust function: <function treutler_atomic_radii_adjust at 0x7fac6e547740>
small_rho_cutoff = 0
Set gradient conv threshold to 3.16228e-05
Initial guess from minao.
tot grids = 33698
init E= -76.2374175220662
  HOMO = -0.432688651289054  LUMO = -0.000586541850088107  gap/eV = 11.75810
cycle= 1 E= -76.2119597533397  delta_E= 0.0255  |g|= 0.754  |ddm|= 1.39
  HOMO = -0.0801446470724816  LUMO = 0.073851221320731  gap/eV = 4.19044
cycle= 2 E= -76.052793398178  delta_E= 0.159  |g|= 1.18  |ddm|= 1.34
  HOMO = -0.273083856722571  LUMO = 0.0310070487513516  gap/eV = 8.27474
cycle= 3 E= -76.3222362299332  delta_E= -0.269  |g|= 0.0274  |ddm|= 0.892
  HOMO = -0.268311627019257  LUMO = 0.0426599492366839  gap/eV = 8.46197
cycle= 4 E= -76.3223972537899  delta_E= -0.000161  |g|= 0.00519  |ddm|= 0.0171
  HOMO = -0.269282488248618  LUMO = 0.0426208997239579  gap/eV = 8.48732
cycle= 5 E= -76.3224012355806  delta_E= -3.98e-06  |g|= 0.000995  |ddm|= 0.0029
  HOMO = -0.269349032424208  LUMO = 0.0426066142257211  gap/eV = 8.48875
cycle= 6 E= -76.3224014295731  delta_E= -1.94e-07  |g|= 2.74e-05  |ddm|= 0.000635
  HOMO = -0.269358285800721  LUMO = 0.0426082051977845  gap/eV = 8.48904
cycle= 7 E= -76.322401459116  delta_E= -2.95e-08  |g|= 1.87e-06  |ddm|= 3.58e-05
  HOMO = -0.269357661744508  LUMO = 0.0426082680384963  gap/eV = 8.48903
cycle= 8 E= -76.3224014529258  delta_E= 6.19e-09  |g|= 1.59e-07  |ddm|= 1.68e-06
  HOMO = -0.269357640040278  LUMO = 0.0426082500414267  gap/eV = 8.48902
cycle= 9 E= -76.3224014496424  delta_E= 3.28e-09  |g|= 1.15e-08  |ddm|= 1.26e-07
  HOMO = -0.269357639430068  LUMO = 0.0426082489118706  gap/eV = 8.48902
cycle= 10 E= -76.3224014533022  delta_E= -3.66e-09  |g|= 9.09e-09  |ddm|= 6.7e-09
  HOMO = -0.269357641154076  LUMO = 0.0426082516987845  gap/eV = 8.48902
cycle= 11 E= -76.3224014537603  delta_E= -4.58e-10  |g|= 4.01e-09  |ddm|= 7.31e-09
  HOMO = -0.269357642247033  LUMO = 0.0426082515652162  gap/eV = 8.48902
Extra cycle  E= -76.3224014501361  delta_E= 3.62e-09  |g|= 6.16e-09  |ddm|= 7.08e-09
converged SCF energy = -76.3224014501361
Dipole moment(X, Y, Z, Debye): -0.00000, -0.00000, -1.98636
Total energy: -2076.838328 eV
Total energy: -76.322401 Hartree

Overwritten attributes  pre_kernel post_kernel  of <class 'skala.pyscf.dft.SkalaRKS'>

Updating Calculator Parameters#

To update the calculator parameters, we can use the set method. For example, we can activate density fitting to speed up the calculation and make the output less verbose. Note that changing certain parameters will reset the calculator and clear previous results.

# Update calculator settings
changed_params = atoms.calc.set(
    with_density_fit=True,
    verbose=0,
    auxbasis="def2-universal-jkfit",
    ks_config={"conv_tol": 1e-6},
)
print(changed_params)
print(f"Changed parameters: {changed_params}")

# Show current parameters
print("\nCurrent calculator parameters:")
for key, value in atoms.calc.parameters.items():
    print(f"  {key}: {value}")

print(atoms.calc)

{'with_density_fit': True, 'verbose': 0, 'auxbasis': 'def2-universal-jkfit', 'ks_config': {'conv_tol': 1e-06}}
Changed parameters: {'with_density_fit': True, 'verbose': 0, 'auxbasis': 'def2-universal-jkfit', 'ks_config': {'conv_tol': 1e-06}}

Current calculator parameters:
  xc: skala-1.1
  basis: def2-svp
  with_density_fit: True
  auxbasis: def2-universal-jkfit
  with_newton: False
  with_dftd3: True
  with_retry: True
  charge: None
  multiplicity: None
  verbose: 0
  ks_config: {'conv_tol': 1e-06}
  device: cpu
<skala.ase.calculator.Skala object at 0x7fac66175f10>

Force Calculations#

We can continue to use the calculator for further calculations, such as computing the forces on the atoms. Since we changed the settings above, the calculator was reset and all previous results have been cleared. Requesting the forces will trigger a new calculation with the updated parameters.

# Calculate forces
forces = atoms.get_forces()  # eV/Å

print("Forces on atoms (eV/Å):")
for i, (symbol, force) in enumerate(
    zip(atoms.get_chemical_symbols(), forces, strict=True)
):
    print(
        f"  Atom {i + 1} ({symbol}): [{force[0]:8.4f}, {force[1]:8.4f}, {force[2]:8.4f}]"
    )

# Calculate maximum force component
max_force = np.max(np.abs(forces))
print(f"\nMaximum force component: {max_force:.4f} eV/Å")

Forces on atoms (eV/Å):
  Atom 1 (O): [ -0.0000,  -0.0000,  -0.0771]
  Atom 2 (H): [  0.0000,  -0.3913,   0.0385]
  Atom 3 (H): [ -0.0000,   0.3913,   0.0385]

Maximum force component: 0.3913 eV/Å

Geometry Optimization#

To use the calculator for geometry optimization or molecular dynamics simulations, you can use the optimize or dynamics modules from ASE. Here we’ll perform a geometry optimization to find the minimum energy structure:

# Store initial geometry for comparison
initial_positions = atoms.positions.copy()
initial_energy = atoms.get_potential_energy()

# Set up and run geometry optimization
opt = Opt(atoms, trajectory="water_opt.traj")
opt.run(fmax=0.01)  # Optimize until forces are below 0.01 eV/Å

# Show optimization results
final_energy = atoms.get_potential_energy()
final_forces = atoms.get_forces()
max_final_force = np.max(np.abs(final_forces))

print(f"Optimization completed in {opt.get_number_of_steps()} steps")
print(f"Initial energy: {initial_energy:.6f} eV")
print(f"Final energy:   {final_energy:.6f} eV")
print(f"Energy change:  {final_energy - initial_energy:.6f} eV")
print(f"Maximum final force: {max_final_force:.4f} eV/Å")

                 Step     Time          Energy          fmax
LBFGSLineSearch:    0 18:17:46    -2076.839069        0.393157
LBFGSLineSearch:    1 18:17:54    -2076.841931        0.220730
LBFGSLineSearch:    2 18:18:09    -2076.844634        0.155762
LBFGSLineSearch:    3 18:18:15    -2076.844883        0.020198
LBFGSLineSearch:    4 18:18:22    -2076.844888        0.000616
Optimization completed in 4 steps
Initial energy: -2076.839069 eV
Final energy:   -2076.844888 eV
Energy change:  -0.005819 eV
Maximum final force: 0.0006 eV/Å

Dipole Moment Calculations#

The Skala calculator also supports dipole moment calculations, which are useful for understanding molecular polarity:

# Calculate dipole moment
dipole = atoms.get_dipole_moment()  # Debye

print(
    f"Dipole moment vector (Debye): [{dipole[0]:8.4f}, {dipole[1]:8.4f}, {dipole[2]:8.4f}]"
)
print(f"Dipole moment magnitude: {np.linalg.norm(dipole):.4f} Debye")

Dipole moment vector (Debye): [  0.0000,  -0.0000,  -0.4190]
Dipole moment magnitude: 0.4190 Debye

Working with Different Molecules and Basis Sets#

Let’s explore how to work with different molecular systems and basis sets. Here we’ll calculate properties for methane (CH₄) using a larger basis set:

# Create a methane molecule
atoms = molecule("CH4")

# Set up calculator with a larger basis set and density fitting
atoms.calc = Skala(
    xc="skala-1.1",
    basis="def2-tzvp",  # Triple-zeta basis set
    with_density_fit=True,  # Enable density fitting for efficiency
    auxbasis="def2-universal-jkfit",
    with_dftd3=True,  # Include dispersion correction
    verbose=1,
)

print("Methane molecule:")
print(f"Chemical formula: {atoms.get_chemical_formula()}")
print(f"Number of atoms: {len(atoms)}")

# Calculate properties
energy = atoms.get_potential_energy()
forces = atoms.get_forces()
max_force = np.max(np.abs(forces))

print("\nCalculation results:")
print(f"Total energy: {energy:.6f} eV")
print(f"Maximum force: {max_force:.6f} eV/Å")

Methane molecule:
Chemical formula: CH4
Number of atoms: 5

Calculation results:
Total energy: -1102.336883 eV
Maximum force: 0.060181 eV/Å

Using Skala on GPU#

In case you have access to a compatible GPU, you can enable GPU acceleration by setting the device parameter to "cuda" when initializing the calculator. This can significantly speed up calculations for larger systems.

calc = Skala(
    xc="skala-1.1",
    basis="def2-tzvp",
    with_density_fit=True,
    auxbasis="def2-universal-jkfit",
    device="cuda",  # Enable GPU acceleration
)

The device can be changed at any time, causing the calculator to reset and use the new device for subsequent calculations.

calc.set(device="cpu")

{'device': 'cpu'}

Summary#

This tutorial covered the essential features of the Skala ASE calculator:

Setting up the calculator with various parameters
Calculating energies, forces, and dipole moments
Updating calculator parameters dynamically
Performing geometry optimizations
Working with different molecules and basis sets