Category: OpenSource

DFT geometry optimizers

These are undeservedly less attention to optimizers than density functionals (concerning Jacob’s ladder). It is not even covered in the recent review: Best-Practice DFT Protocols for Basic Molecular Computational Chemistry. At the same time, in my current projects, the most resource-demanding was geometry optimization – the time spent on optimizing structures was much longer than a single-point calculation. Papers that introduce new (AI-based) optimizers promise significant speed-up. However, there are always some problems:

  1. The tested systems are different from my electrochemical interfaces.
  2. The code is not available or difficult to install.
  3. The code is outdated and contains bugs.
  4. Optimizers perform worse than the common ones, like QuasiNewton in ASE.

ASE wiki lists all internal and some external optimizers and provides their comparison. I have checked the most promising on a high-entropy alloy slab.

Observation 1. QuasiNewton outperforms all other optimizers. Point. I have run a standard GPAW/DFT/PBE/PW optimization with various optimizers:

Observation 2. Pre-optimizing the slab with a cheaper method does not reduce the number of optimization steps. I have preoptimized the geometry with TBLITE/DFTB/GFN1-xTB to continue with GPAW/DFT/PBE/PW. Preoptimization takes just some minutes and the obtained geometry looks similar to the DFT one but that does not reduce the number of DFT optimization steps.

OptimizerN steps*Time$N steps*#Total time#
BFGS1602:44:271703:01:26
LBFGS1502:30:351602:55:04
BondMin1202:46:271302:45:07
GPMin1205:26:233108:14:22
MLMin38verylong2812:31:29
FIRE3805:06:564405:56:54
QuasiNewton801:36:23902:00:10

Note * – the printed number of steps might different from the actuall number of calculations because each calculator has a different way of reporting that number.

Note $ – the time between the end of the first and last steps.

Note # – started from the TBLITE/DFTB/GFN1-xTB preoptimized geometry.

N.B! I have done my test only once in two runs: starting with slab.xyz and preoptized geometry. Runs were on similar nodes and all optimizations were done on the same node.

Conclusion. Do not believe in claims in articles advertizing new optimizers – Run your tests before using them.

A practical finding. The usual problem with calculations that require many optimization steps is that they need to fit into HPC time limits. On the restart, ASE usually rewrites the trajectory. Some optimizers (GPMin and AI-based) could benefit from reading the full trajectory. So, I started writing two trajectories and a restart file like this.

# Restarting
if os.path.exists(f'{name}_last.gpw') == True and os.stat(f'{name}_last.gpw').st_size > 0:
    atoms,calc = restart(f'{name}_last.gpw', txt=None)
    parprint(f'Restart from the gpw geometry.')
elif os.path.exists(f'{name}_full.traj') == True and os.stat(f'{name}_full.traj').st_size > 0:
    atoms = read(f'{name}_full.traj',-1)
    parprint(f'Restart with the traj geometry.')
else:
    atoms = read(f'{name}_init.xyz')
    parprint(f'Start with the initial xyz geometry.')

# Optimizing
opt = QuasiNewton(atoms, trajectory=f'{name}.traj', logfile=f'{name}.log')
traj= Trajectory(f'{name}_full.traj', 'a', atoms)
opt.attach(traj.write, interval=1)
def writegpw():
    calc.write(f'{name}_last.gpw')
opt.attach(writegpw, interval=1)
opt.run(fmax=0.05, steps=42)

Here are some details on the tests.

My gpaw_opt.py for DFT calculations on 24 cores:

# Load modules
from ase import Atom, Atoms
from ase.build import add_adsorbate, fcc100, fcc110, fcc111, fcc211, molecule
from ase.calculators.mixing import SumCalculator
from ase.constraints import FixAtoms, FixedPlane, FixInternals
from ase.data.vdw_alvarez import vdw_radii
from ase.db import connect
from ase.io import write, read
from ase.optimize import BFGS, GPMin, LBFGS, FIRE, QuasiNewton
from ase.parallel import parprint
from ase.units import Bohr
from bondmin import BondMin
from catlearn.optimize.mlmin import MLMin
from dftd4.ase import DFTD4
from gpaw import GPAW, PW, FermiDirac, PoissonSolver, Mixer, restart
from gpaw.dipole_correction import DipoleCorrection
from gpaw.external import ConstantElectricField
from gpaw.utilities import h2gpts
import numpy as np
import os

atoms = read('slab.xyz')
atoms.set_constraint([FixAtoms(indices=[atom.index for atom in atoms if atom.tag in [1,2]])])

# Set calculator
kwargs = dict(poissonsolver={'dipolelayer':'xy'},
              xc='RPBE',
              kpts=(4,4,1),
              gpts=h2gpts(0.18, atoms.get_cell(), idiv=4),
              mode=PW(400),
              basis='dzp',
              parallel={'augment_grids':True,'sl_auto':True,'use_elpa':True},
             )
calc = GPAW(**kwargs)

#atoms.calc = SumCalculator([DFTD4(method='RPBE'), calc])
#atoms.calc = calc

# Optimization paramters
maxf = 0.05

# Run optimization
###############################################################################

# 2.A. Optimize structure using MLMin (CatLearn).
initial_mlmin = atoms.copy()
initial_mlmin.set_calculator(calc)
mlmin_opt = MLMin(initial_mlmin, trajectory='results_mlmin.traj')
mlmin_opt.run(fmax=maxf, kernel='SQE', full_output=True)

# 2.B Optimize using GPMin.
initial_gpmin = atoms.copy()
initial_gpmin.set_calculator(calc)
gpmin_opt = GPMin(initial_gpmin, trajectory='results_gpmin.traj', logfile='results_gpmin.log', update_hyperparams=True)
gpmin_opt.run(fmax=maxf)

# 2.C Optimize using LBFGS.
initial_lbfgs = atoms.copy()
initial_lbfgs.set_calculator(calc)
lbfgs_opt = LBFGS(initial_lbfgs, trajectory='results_lbfgs.traj', logfile='results_lbfgs.log')
lbfgs_opt.run(fmax=maxf)

# 2.D Optimize using FIRE.
initial_fire = atoms.copy()
initial_fire.set_calculator(calc)
fire_opt = FIRE(initial_fire, trajectory='results_fire.traj', logfile='results_fire.log')
fire_opt.run(fmax=maxf)

# 2.E Optimize using QuasiNewton.
initial_qn = atoms.copy()
initial_qn.set_calculator(calc)
qn_opt = QuasiNewton(initial_qn, trajectory='results_qn.traj', logfile='results_qn.log')
qn_opt.run(fmax=maxf)

# 2.F Optimize using BFGS.
initial_bfgs = atoms.copy()
initial_bfgs.set_calculator(calc)
bfgs_opt = LBFGS(initial_bfgs, trajectory='results_bfgs.traj', logfile='results_bfgs.log')
bfgs_opt.run(fmax=maxf)

# 2.G. Optimize structure using BondMin.
initial_bondmin = atoms.copy()
initial_bondmin.set_calculator(calc)
bondmin_opt = BondMin(initial_bondmin, trajectory='results_bondmin.traj',logfile='results_bondmin.log')
bondmin_opt.run(fmax=maxf)

# Summary of the results
###############################################################################

fire_results = read('results_fire.traj', ':')
parprint('Number of function evaluations using FIRE:',
         len(fire_results))

lbfgs_results = read('results_lbfgs.traj', ':')
parprint('Number of function evaluations using LBFGS:',
         len(lbfgs_results))

gpmin_results = read('results_gpmin.traj', ':')
parprint('Number of function evaluations using GPMin:',
         gpmin_opt.function_calls)

bfgs_results = read('results_bfgs.traj', ':')
parprint('Number of function evaluations using BFGS:',
         len(bfgs_results))

qn_results = read('results_qn.traj', ':')
parprint('Number of function evaluations using QN:',
         len(qn_results))

catlearn_results = read('results_mlmin.traj', ':')
parprint('Number of function evaluations using MLMin:',
         len(catlearn_results))

bondmin_results = read('results_bondmin.traj', ':')
parprint('Number of function evaluations using BondMin:',
         len(bondmin_results))

Initial slab.xyz file: 

45
Lattice="8.529357696932532 0.0 0.0 4.264678848466266 7.386640443507905 0.0 0.0 0.0 29.190908217261956" Properties=species:S:1:pos:R:3:tags:I:1 pbc="T T F"
Ir       0.00000000       1.62473838      10.00000000        5
Ru       2.81412943       1.62473838      10.00000000        5
Pt       5.62825885       1.62473838      10.00000000        5
Pd       1.40706471       4.06184595      10.00000000        5
Ag       4.22119414       4.06184595      10.00000000        5
Ag       7.03532356       4.06184595      10.00000000        5
Ag       2.81412943       6.49895353      10.00000000        5
Ru       5.62825885       6.49895353      10.00000000        5
Pt       8.44238828       6.49895353      10.00000000        5
Pt       0.00000000       0.00000000      12.29772705        4
Ag       2.81412943       0.00000000      12.29772705        4
Ru       5.62825885       0.00000000      12.29772705        4
Ru       1.40706471       2.43710757      12.29772705        4
Ir       4.22119414       2.43710757      12.29772705        4
Ag       7.03532356       2.43710757      12.29772705        4
Ag       2.81412943       4.87421514      12.29772705        4
Ir       5.62825885       4.87421514      12.29772705        4
Pd       8.44238828       4.87421514      12.29772705        4
Pd       1.40706471       0.81236919      14.59545411        3
Ir       4.22119414       0.81236919      14.59545411        3
Pt       7.03532356       0.81236919      14.59545411        3
Ag       2.81412943       3.24947676      14.59545411        3
Ir       5.62825885       3.24947676      14.59545411        3
Ir       8.44238828       3.24947676      14.59545411        3
Pd       4.22119414       5.68658433      14.59545411        3
Pt       7.03532356       5.68658433      14.59545411        3
Ag       9.84945299       5.68658433      14.59545411        3
Pd       0.00000000       1.62473838      16.89318116        2
Pd       2.81412943       1.62473838      16.89318116        2
Ag       5.62825885       1.62473838      16.89318116        2
Pt       1.40706471       4.06184595      16.89318116        2
Ag       4.22119414       4.06184595      16.89318116        2
Ag       7.03532356       4.06184595      16.89318116        2
Ru       2.81412943       6.49895353      16.89318116        2
Ru       5.62825885       6.49895353      16.89318116        2
Ru       8.44238828       6.49895353      16.89318116        2
Ir       0.00000000       0.00000000      19.19090822        1
Ag       2.81412943       0.00000000      19.19090822        1
Pt       5.62825885       0.00000000      19.19090822        1
Pd       1.40706471       2.43710757      19.19090822        1
Ag       4.22119414       2.43710757      19.19090822        1
Pd       7.03532356       2.43710757      19.19090822        1
Ag       2.81412943       4.87421514      19.19090822        1
Ru       5.62825885       4.87421514      19.19090822        1
Ir       8.44238828       4.87421514      19.19090822        1

My tblite_opt.py for DFTB calcualation with just one core. It takes some minutes but eventually crashes 🙁

# Load modules
from ase import Atom, Atoms
from ase.build import add_adsorbate, fcc100, fcc110, fcc111, fcc211, molecule
from ase.calculators.mixing import SumCalculator
from ase.constraints import FixAtoms, FixedPlane, FixInternals
from ase.data.vdw_alvarez import vdw_radii
from ase.db import connect
from ase.io import write, read
from ase.optimize import BFGS, GPMin, LBFGS, FIRE, QuasiNewton
from ase.parallel import parprint
from ase.units import Bohr
from tblite.ase import TBLite
import numpy as np
import os

# https://tblite.readthedocs.io/en/latest/users/ase.html

atoms = read('slab.xyz')
atoms.set_constraint([FixAtoms(indices=[atom.index for atom in atoms if atom.tag in [1,2]])])

# Set calculator
calc = TBLite(method="GFN1-xTB",accuracy=1000,electronic_temperature=300,max_iterations=300)
atoms.set_calculator(calc)
qn_opt = QuasiNewton(atoms, trajectory='results_qn.traj', logfile='results_qn.log', maxstep=0.1)
qn_opt.run(fmax=0.1)

To compare structures I have used MDanalysis, which unfortunately does not work with ASE traj, so I prepared xyz-files with “ase convert -n -1 file.traj file.xyz”

import MDAnalysis as mda
from MDAnalysis.analysis.rms import rmsd
import sys

def coord(file_name):
    file  = mda.Universe(f"{file_name}.xyz")
    atoms = file.select_atoms("index 1:9")
    return  atoms.positions.copy()

print(rmsd(coord(sys.argv[1]),coord(sys.argv[2])))

An instruction on installation of GPAW. TBLITE can be installed as “conda install -c conda-forge tblite”.

GPAW installation with pip

Between installation with conda and compilation of libraries, an intermediate path – installation of GPAW with pip – is a compromise for those who wish to text specific GPAW branches or packages.

For example, I wish to text self-interaction error correction (SIC) and evaluate Bader charges with pybader. Neither SIC nor pybader is compatible with the recent GPAW. Here is not to get a workable version.

# numba in pybader is not compatible with python 3.11, so create a conda environment with python 3.10
conda create -n gpaw-pip python=3.10 
conda activate gpaw-pip

conda install -c conda-forge libxc libvdwxc
conda install -c conda-forge ase
# ensure that you install the right openmpi (not external)
conda install -c conda-forge openmpi ucx
conda install -c conda-forge compilers
conda install -c conda-forge openblas scalapack
conda install -c conda-forge pytest
pip install pybader

# Get a developer version of GPAW with SIC
git clone -b dm_sic_mom_update https://gitlab.com/alxvov/gpaw.git
cd gpaw
cp siteconfig_example.py siteconfig.py

# In the siteconfig.py rewrite
'''
fftw = True
scalapack = True
if scalapack:
    libraries += ['scalapack']
'''

unset CC
python -m pip install -e .
gpaw info

Visualizing ASE atoms in Jupyter notebooks

For a long time I wanted to see ASE atoms in my Jupyter notebook. My previous attempts were usually unsuccessful. Today I decided to try again. First ASE wiki suggests x3d and webngl:

view(atoms, viewer='x3d')
view(atoms, viewer='ngl')

Łucasz Mentel gives some useful tips in his blogpost from 2017.

In my case x3d works and webngl fails. The x3d picture is not enought, and I do not want to spend much time on fixing webngl.

ASE-notebook is what works for me.

conda create -n "jupyter"
conda activate jupyter
conda install -c conda-forge ase-notebook
conda install -c conda-forge jupyterlab

By the way, the model is from my “Surface Curvature Effect on Dual-Atom Site Oxygen Electrocatalysis” paper, which you can read at chemRxiv until it turns Gold Open Access.

Bayesian Error Estimation (BEE) for RPBE

The Beyesian Error Estimation (BEE) is implemented in GPAW only for PBE, BEEF-vdW, and mBEEF-vdW.

Here is a trick for making the BEE with the RPBE functional. Just edit the lines in ASE and GPAW codes by adding RPBE as an exception.

To find the needed files, run

find ./ -name "bee.py"

In ase/dft/bee.py change one line:

class BEEFEnsemble:



            if self.xc in ['BEEF-vdW', 'BEEF', 'PBE', 'RPBE']: # add RPBE
                self.beef_type = 'beefvdw'

In gpaw/xc/bee.py add two lines:

class BEEFEnsemble:
    """BEEF ensemble error estimation."""
    def __init__(self, calc):



        # determine functional and read parameters
        self.xc = self.calc.get_xc_functional()
        if self.xc == 'BEEF-vdW':
            self.bee_type = 1
        elif self.xc == 'RPBE': # catch the RPBE exchange functional
            self.bee_type = 1   # assign BEEF coefficients the RBPE

Below we use BEEF-vdW, RPBE, and PBE dimensionless density (n) with gradient (s) and apply BEEF coefficients (E₀, ΔEᵢ) to evaluate the BEE as the standard deviation for the ensemble total energies with the variable enhancement factor (F(s,θᵢ)).


from ase import Atoms
from ase.dft.bee import BEEFEnsemble
from ase.parallel import parprint
from gpaw import GPAW
import time

for xc in ['BEEF-vdW','RPBE','PBE']:
    start_time = time.time()

    h2 = Atoms('H2',[[0.,0.,0.],[0.,0.,0.741]]) #exp. bond length
    h2.center(vacuum=3)
    cell = h2.get_cell()

    calc = GPAW(xc=xc,txt='H2_{0}.txt'.format(xc))
    h2.calc = calc
    e_h2 = h2.get_potential_energy()
    ens = BEEFEnsemble(calc)
    de_h2 = ens.get_ensemble_energies()
    del h2, calc, ens

    h = Atoms('H')
    h.set_cell(cell)
    h.center()
    calc = GPAW(xc=xc,txt='H_{0}.txt'.format(xc), hund=True)
    h.calc = calc
    e_h = h.get_potential_energy()
    ens = BEEFEnsemble(calc)
    de_h = ens.get_ensemble_energies()
    del h, calc, ens

    E_bind = 2*e_h - e_h2
    dE_bind = 2*de_h[:] - de_h2[:]
    dE_bind = dE_bind.std()
    
    parpting('{0} functional'.format(xc))
    parprint('Time: {0} s'.format(round(time.time()-start_time,0)))
    parprint('E_bind: {0} eV'.format(round(E_bind,4)))
    parprint('Error bar {0} eV'.format(round(dE_bind,4)))

Synchronizing calendars

This post describes ways of pushing MS outlook and Google calendars to Nextcloud.

My main working calendar is the Nextcloud app because I can easily sync it is my Sailfish phone. I also use Google calendar (for sharing family events) and MS outlook calendar (for work). Today I decided to merge all these calendars into a single one that I can sync on all my devices. Here is how.

MS to Nextcloud

Use outlookcaldavsynchronizer as recommended in the Nextcloud blog.

Google to Nextcould

  • Get the iCal link from Google Calendar as follows:
  • In the left calendar list menu of Google Calendar, go to the ⋮ menu of the calendar to be shared
  • Click on “Settings and sharing”
  • On the Calendar settings page, scroll down to “Secret address in iCal format”
  • In Nextcloud Calendar’s left menu, click on “+New Calendar” > “New subscription from link (read-only)”
  • Insert the “Secret address in iCal format”
  • Your new calendar subscription will appear in the list; you can change its name or color in the menu of your calendar