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”.

Erasmus

A group of students and tutors from Paris 13 visited Chemicum. Aurélie and Chris on the left and David with sunglasses in the middle. Georgi (on the right) made a nice excursion. Next year two Erastus students from Paris 13 will be studying in Tartu!

Visiting Tartu in November

A list of things to bring with you:

  • Gloves, hat and scarf (the average temperature is −1.5°C)
  • Waterproof boots or trekking boots (with a good grip in case there is ice)
  • Layered clothing (like pullovers and cardigans, so that you can remove or add layers according to the weather and how fast you are moving)
  • Swimming equipment (for SPA and sauna or why not doing some winter swimming?)
  • Napkins for a runny nose
  • A postcard to pin in the office 5072 where Vladislav works

A valuable message

My friend Olga Jasnovidova has recently finished her PhD studies in Brno in Czech republic. In her acknowledgement speech she concluded the 7 years-long work and gave a message for young ones. Here it is.

“””
“Science is one of the most creative fields to work in”

Through my PhD studies I have learned that science has two dimensions: one scientific, and one human… and that, of the two, the human dimension is the more difficult to grasp.

I understood that, in order to achieve your goals, you must not only care about your work or yourself. You have got to care about the people around you: senior and junior students, technicians, facility managers, senior colleagues, your supervisors. As you step into PhD studies, you should not expect them to support and motivate you, but rather you yourself should start by supporting and motivating them.

I have also learned that one of the most effective ways to grow and develop is to ask for feedback, then learn to accept it and moreover learn to give constructive feedback yourself. This process can be painful for our egos, but it is the only way to grow.

Most importantly, I have learned that there are always several ways to achieve the same goal. There is no one correct way to reach one’s target. Therefore, you should always stay open to new ways to achieve your goals.

Based on my experience, I would like to say to younger students that they should not fear anything new or unfamiliar. Please don’t create any mental barriers for yourselves. Academic research seems like a very conventional and strict field after Bachelor’s or Master’s studies or peer review, but it is not so. Science is one of the most creative fields to work in, providing endless opportunities to grow and discover. Just learn the rules and then use them to create. You can do something truly unique for the first time in human history, something that will lay a path for many to follow. Believe in your own abilities. Every single day you have an opportunity to do something amazing — don’t waste it.
“””

P.S. Olga, congratulations!

P.P.S. Meanwhile Samual Coles defended his PhD thesis in Oxford. Grande Sam! Sam, congratulations!

The second visit of Iuliia

By Iuliia:

“””

I have been collaborating with the Electrical Double Layer group from the University of Tartu since the beginning of 2016. I had been to Tartu once, in March’2016, and this August I have visited the group again. During this visit, I was accompanied by Dr. Marco Preto, Researcher in Novelmar Project from Interdisciplinary Centre of Marine and Environmental Research of the University of Porto.
The host institution received us very warmly. There was no need to settle any bureaucracy procedures – Estonian efficiency does not cease to amaze me. Everything was taken care of in advance, and we immediately got out a working spaces, keys or anything we could need for work. I think such attitude is very important for these short visits.
In Estonia we spent two wonderful weeks with work and leisure interconnected. Most of the time in Tartu we worked closely with Dr. Vladislav Ivaništšev and his team, where very productive work was carried out, with social activity interludes that recharged us with a relaxed exchange of ideas. During this visit, the work on developing of an approach to an analysis of electrical double layer in ionic liquids systems was conducted, and an article on our previously done work was prepared for submission.
Among all the Master and PhD students, that are being trained at the group, Meeri Lembinen must be acknowledged especially. Meeri, besides being a brilliant student, is a perfect manager. I suspect, due to her care and attention we have not got a single problem at the university and during the whole stay were accompanied by her and felt like at home.
I hope, our fruitful collaboration is to be continued!

“””

Iuliia (left), Meeri (right)

Predictions of Physicochemical Properties of Ionic Liquids with DFT

We have published an article in Computation, where we discuss predicting both static and dynamic physicochemical properties of ionic liquids with density functional theory. We prepared a workflow using NaRIBaS and conducted the calculations with 48 common ionic pairs. Thence, we estimated relevant properties for practical electrochemical applications, such as electrochemical stability and viscosity.

This work is an example how a simplistic computational model can be used in combination with informatics techniques to obtain relevant information about the ionic liquids. It can be found at the MDPI Journal website (Computation 2016, 4(3), 25; doi:10.3390/computation4030025).


Figure 6 from the article.
Figure 6 from the article. The relation between the electrochemical stability window (EW) and the estimated activation energy of viscosity (Eaest).

Sergey Sosnin’s visit to Tartu

I was in a three-days academic trip in December 2016 within a collaboration between Skoltech research group of Prof. Maxim Fedorov and Tartu theoretical electrochemistry group. My primary goal was to study 3D-RISM method for calculation of the solvatation of molecules and applying this method in my QSAR researches. During my visit I’ve worked closely with Maksim Mišin. He has trained me the 3D-RISM methodology, workflow, and theoretical background. After that we have done a small project related with application of RISM in deep learning for chemoinformatics tasks. I was inspired by good theoretical and practical skills of researchers in the theoretical electrochemistry group. We discussed a lot about current science challenges, and I hope that mutual exchange of ideas was helpful for both.

Because it was my first visit to Tartu I had have intended to familiarize myself with this town. Maksim Mišin and Dr. Vladislav Ivaništšev have represented me this city and have told me a lot of interesting things about Estonia. I was really enjoyed in this trip, and I hope that I will visit it again.

Sergey Sosnin

Soaked to the skin: tuning ionic liquids for electrochemical devices

A post in JPhys+ about our recent article: Researchers at the Universities of Santiago de Compostela, A Coruña, Tartu, Stratchclyde and Cambridge, shed light on the structure of the electrified interface in mixtures of contaminated ionic liquids in their recently published JPCM letter.

https://jphysplus.iop.org/2016/10/28/soaked-to-the-skin-tuning-ionic-liquids-for-electrochemical-devices/

The article is published here: http://iopscience.iop.org/article/10.1088/0953-8984/28/46/464001

P.S. This publication became possible thanks to the COST CM1206.

P.P.S. Finally they have corrected the authors names!