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 and periodic systems.
import numpy as np
from ase.build import molecule
from ase.optimize import LBFGSLineSearch as Opt
from ase.units import Bohr, Hartree
from skala.ase import Skala
/home/giulialuise/miniconda3/envs/skala/lib/python3.13/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
from .autonotebook import tqdm as notebook_tqdm
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 calculationswith_dftd3: Include DFT-D3 dispersion correction (default: True)verbose: Control output verbosity (0-4)
# Create a water molecule
atoms = molecule("H2O")
# Set up the Skala calculator with specific parameters
atoms.calc = Skala(xc="skala", 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
basis: def2-svp
with_density_fit: False
with_newton: False
with_dftd3: True
charge: None
multiplicity: None
verbose: 4
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='GCRAZGDL1498', release='6.8.0-1034-azure', version='#39~22.04.1-Ubuntu SMP Wed Aug 13 22:25:47 UTC 2025', machine='x86_64') Threads 24
Python 3.13.7 | packaged by conda-forge | (main, Sep 3 2025, 14:30:35) [GCC 14.3.0]
numpy 2.3.3 scipy 1.16.2 h5py 3.14.0
Date: Tue Oct 7 15:43:37 2025
PySCF version 2.9.0
PySCF path /home/giulialuise/miniconda3/envs/skala/lib/python3.13/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: 1.87
/home/giulialuise/.cache/huggingface/hub/models--microsoft--skala/snapshots/9eedaca41deeaebe197bec77d1898f777770d3ef/skala-1.0.fun
Downloaded in 0.16 seconds
******** <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/tmpcw83w2ua
max_memory 4000 MB (current use 736 MB)
XC library 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 0x764481e2b560>
grids dens level: 3
symmetrized grids: False
atomic radii adjust function: <function treutler_atomic_radii_adjust at 0x764476f1aca0>
small_rho_cutoff = 1e-07
Set gradient conv threshold to 3.16228e-05
Initial guess from minao.
tot grids = 33704
init E= -76.2302482032812
HOMO = -0.430502394424911 LUMO = -0.00557962018270088
cycle= 1 E= -76.1844346685567 delta_E= 0.0458 |g|= 0.818 |ddm|= 1.41
HOMO = -0.0571497729893296 LUMO = 0.092111968345332
cycle= 2 E= -75.9921484773828 delta_E= 0.192 |g|= 1.28 |ddm|= 1.4
HOMO = -0.266575634970147 LUMO = 0.0399392427556757
cycle= 3 E= -76.3081751739264 delta_E= -0.316 |g|= 0.0337 |ddm|= 0.956
HOMO = -0.25652282325852 LUMO = 0.0420512491571007
cycle= 4 E= -76.30841420647 delta_E= -0.000239 |g|= 0.0064 |ddm|= 0.0228
HOMO = -0.257181682273231 LUMO = 0.0406822667985362
cycle= 5 E= -76.3084203471352 delta_E= -6.14e-06 |g|= 0.00122 |ddm|= 0.00337
HOMO = -0.257376647610882 LUMO = 0.0405720409497939
cycle= 6 E= -76.3084205799498 delta_E= -2.33e-07 |g|= 3.1e-05 |ddm|= 0.000715
HOMO = -0.257386810929883 LUMO = 0.0405644826339652
cycle= 7 E= -76.3084206025217 delta_E= -2.26e-08 |g|= 3.89e-06 |ddm|= 3.08e-05
HOMO = -0.257387527784119 LUMO = 0.0405648211980371
cycle= 8 E= -76.3084205960069 delta_E= 6.51e-09 |g|= 3.05e-07 |ddm|= 3.24e-06
HOMO = -0.25738749260452 LUMO = 0.0405648168862294
cycle= 9 E= -76.3084206012037 delta_E= -5.2e-09 |g|= 3.17e-08 |ddm|= 2.17e-07
HOMO = -0.257387492441066 LUMO = 0.0405648161651567
cycle= 10 E= -76.3084206024703 delta_E= -1.27e-09 |g|= 2.04e-08 |ddm|= 1.73e-08
HOMO = -0.257387491334681 LUMO = 0.0405648167913778
cycle= 11 E= -76.3084206001009 delta_E= 2.37e-09 |g|= 1.44e-08 |ddm|= 7.84e-09
HOMO = -0.257387491019448 LUMO = 0.0405648174191108
cycle= 12 E= -76.3084205990203 delta_E= 1.08e-09 |g|= 1.38e-08 |ddm|= 4.07e-09
HOMO = -0.257387490637195 LUMO = 0.0405648176151493
cycle= 13 E= -76.3084206006191 delta_E= -1.6e-09 |g|= 2.06e-08 |ddm|= 7.24e-09
HOMO = -0.257387491245932 LUMO = 0.0405648173875793
cycle= 14 E= -76.3084206022113 delta_E= -1.59e-09 |g|= 1.76e-08 |ddm|= 4.61e-09
HOMO = -0.257387491753683 LUMO = 0.0405648174842908
cycle= 15 E= -76.3084205999321 delta_E= 2.28e-09 |g|= 2.05e-08 |ddm|= 4.88e-09
HOMO = -0.25738749162096 LUMO = 0.0405648180109603
cycle= 16 E= -76.3084206010219 delta_E= -1.09e-09 |g|= 1.75e-08 |ddm|= 3.15e-09
HOMO = -0.257387492013501 LUMO = 0.0405648177311237
cycle= 17 E= -76.3084206001081 delta_E= 9.14e-10 |g|= 1.46e-08 |ddm|= 3.2e-09
HOMO = -0.257387493925225 LUMO = 0.0405648195794535
Extra cycle E= -76.3084206016336 delta_E= -1.53e-09 |g|= 1.48e-08 |ddm|= 1.74e-08
converged SCF energy = -76.3084206016336
Dipole moment(X, Y, Z, Debye): -0.00000, 0.00000, -2.00440
Total energy: -2076.457890 eV
Total energy: -76.308421 Hartree
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)
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}
Changed parameters: {'with_density_fit': True, 'verbose': 0}
Current calculator parameters:
xc: skala
basis: def2-svp
with_density_fit: True
with_newton: False
with_dftd3: True
charge: None
multiplicity: None
verbose: 0
<skala.ase.calculator.Skala object at 0x7643f91ae900>
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)):
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/Å")
/home/giulialuise/.cache/huggingface/hub/models--microsoft--skala/snapshots/9eedaca41deeaebe197bec77d1898f777770d3ef/skala-1.0.fun
Downloaded in 0.10 seconds
Forces on atoms (eV/Å):
Atom 1 (O): [ 0.0000, 0.0000, -0.5867]
Atom 2 (H): [ -0.0000, -0.6148, 0.2933]
Atom 3 (H): [ 0.0000, 0.6148, 0.2933]
Maximum force component: 0.6148 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 15:43:56 -2076.458638 0.681230
LBFGSLineSearch: 1 15:44:02 -2076.466242 0.308524
LBFGSLineSearch: 2 15:44:13 -2076.467743 0.153213
LBFGSLineSearch: 3 15:44:26 -2076.468485 0.002282
Optimization completed in 3 steps
Initial energy: -2076.458638 eV
Final energy: -2076.468485 eV
Energy change: -0.009847 eV
Maximum final force: 0.0023 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.4200]
Dipole moment magnitude: 0.4200 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",
basis="def2-tzvp", # Triple-zeta basis set
with_density_fit=True, # Enable density fitting for efficiency
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(f"\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
/home/giulialuise/.cache/huggingface/hub/models--microsoft--skala/snapshots/9eedaca41deeaebe197bec77d1898f777770d3ef/skala-1.0.fun
Downloaded in 0.11 seconds
Calculation results:
Total energy: -1102.070344 eV
Maximum force: 0.351115 eV/Å
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