Optimising even simple intermediate on metal surfaces might be tricky. They tend change the adsorption site. OOH also very flexible in changing geometry. Thus, some optimizers get stuck (like BFGSLineSearch) while others just crush (like GPMin).
Here is my test with OCPCalculator and EquiformerV2-31M-S2EF-OC20-All+MD (on 1 CPU in Google Colab). Optimisation of OH on Pt(111) is simple.
Here SciPyFminBFGS failed and other optimizers ended up with the same structure, shown below.
Optimisation of OOH on Pt(111) is challenging. GoodOldQuasiNewton, FIRE, FIRE2, MDMin, and BFGSLineSearch (QuasiNewton) do not converge in 95 cycles (which I set as a maximum number of cycles). GPMin, SciPyFminBFGS, SciPyFminCG, and LBFGSLineSearch failed.
Geometry obtained with FIRE2 (similar to BFGS and FIRE)
Geometry obtained with LBFGS. It is clearly far from the optimum.
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:
The tested systems are different from my electrochemical interfaces.
The code is not available or difficult to install.
The code is outdated and contains bugs.
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.
Optimizer
N steps*
Time$
N steps*#
Total time#
BFGS
16
02:44:27
17
03:01:26
LBFGS
15
02:30:35
16
02:55:04
BondMin
12
02:46:27
13
02:45:07
GPMin
12
05:26:23
31
08:14:22
MLMin
38
verylong
28
12:31:29
FIRE
38
05:06:56
44
05:56:54
QuasiNewton
8
01:36:23
9
02: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))
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”
TS09 and D4 are atomic-charge dependent dispersion corrections (see TS09 PRL paper and D4 homepage for the refs). The D4 code is available at github. According to GPAW documentation, TS09 and D4 show for the S26 test set smaller mean deviation than vdW-DF. Herewith, D4 correction does not depend on the actual calculation as it is added to the calculated energy.
Here is how D4 correction can be added with ASE (see Readme) after installing it (for example, as conda install -c conda-forge dftd4 dftd4-python):
from ase.build import molecule
from ase.calculators.mixing import SumCalculator
from ase.optimize import BFGS
from dftd4.ase import DFTD4
from gpaw import GPAW
atoms = molecule('H2O')
atoms.center(vacuum=4)
gpaw = GPAW(mode='fd',txt='H2O_D4.txt',xc='PBE')
atoms.calc = SumCalculator([DFTD4(method='PBE'), gpaw])
#atoms.get_potential_energy()
opt = BFGS(atoms,trajectory='H2O_D4.traj', logfile='H2O_D4.log')
opt.run(fmax=0.05)
Let me stress that before choosing TS09 or D4 one should consider all pro and contra. TS09 method used Hirshfeld charges while D4 uses the electronegativity equilibration method to obtain charges. The former naturally accounts for the interfacial charge transfer while the latter does not. The TS09 correction requires vdW radii and is implemented for a limited set on functionals (see ASE code), like PBE, RPBE, and BLYP. The D4 correction supports much more functionals (see parameters). Regarding the vdW radii values for TS09 bare in mind that there are four data sources – one in GPAW, two in ASE and one more in ASE.
Here is how TS09 correction can be added with ASE and GPAW:
from ase.build import molecule
from ase.calculators.vdwcorrection import vdWTkatchenko09prl
from ase.data.vdw_alvarez import vdw_radii
from ase.optimize import BFGS
from gpaw.analyse.hirshfeld import HirshfeldPartitioning
from gpaw.analyse.vdwradii import vdWradii
from gpaw import GPAW
atoms = molecule('H2O')
atoms.center(vacuum=4)
gpaw = GPAW(mode='fd',txt='H2O_TS.txt',xc='PBE')
atoms.calc = vdWTkatchenko09prl(HirshfeldPartitioning(gpaw), vdWradii(atoms.get_chemical_symbols(), 'PBE'))
#atoms.get_potential_energy()
opt = BFGS(atoms,trajectory='H2O_TS.traj', logfile='H2O_TS.log')
opt.run(fmax=0.05)
N.B! Note that the TS09 and D4 energies are no outputted to the H2O.txt. They are written to the log-file.
Bader analysis is a fast and simple way of getting atomic charges. It is especially useful for periodic calculations. The analysis can be done on the fly using pybader tool from pypi.org/project/pybader. I recommend using it within conda environments.
The installation is straightforward:
pip install pybader
The usage is less obvious. Here is a function for obtaining xmol-type xyz that can obtained with GPAW and visualized with ASE:
def xyzb(atoms, filename, nCPU):
from pybader.io import gpaw
from pybader.interface import Bader
import os
bader = Bader(*gpaw.read_obj(atoms.calc)) # read ASE object 'atoms'
bader(threads=nCPU) # specify the number of CPUs
f = open('{0}.xyz'.format(filename), 'w') # set xmol format
b = bader.atoms_charge # get number of electrons per atom
n = atoms.get_atomic_numbers() # get atomic numbers
a = atoms.get_chemical_symbols() # get chemical symbols
p = atoms.get_positions() # get positions of the atoms
f.write('{0}\n'.format(len(a)))
f.write('Properties=species:S:1:pos:R:3:charge:R:1\n') # ensure compatibility with ASE
for i in range(len(a)): # print symbol, positions, and charge
s = '{0}'.format(a[i])
x = '{0:.6f}'.format(round(p[i][0],6))
y = '{0:.6f}'.format(round(p[i][1],6))
z = '{0:.6f}'.format(round(p[i][2],6))
c = '{0:.3f}'.format(round(float(n[i]) - float(b[i]),3))
f.write('{0:<4}{1:>16}{2:>16}{3:>16}{4:>10}\n'.format(s,x,y,z,c))
f.close()
del bader
os.remove('bader.p')
Similarly one can obtain xmol-type xyz file with Hirshfeld charges:
def xyzh(atoms, filename):
from gpaw.analyse.hirshfeld import HirshfeldPartitioning
f = open('{0}.xyz'.format(filename), 'w') # set xmol format
a = atoms.get_chemical_symbols() # get chemical symbols
p = atoms.get_positions() # get positions of the atoms
h = HirshfeldPartitioning(atoms.calc).get_charges()
f.write('{0}\n'.format(len(a)))
f.write('Properties=species:S:1:pos:R:3:charge:R:1\n') # ensure compatibility with ASE
for i in range(len(a)): # print symbol, positions, and charge
s = '{0}'.format(a[i])
x = '{0:.6f}'.format(round(p[i][0],6))
y = '{0:.6f}'.format(round(p[i][1],6))
z = '{0:.6f}'.format(round(p[i][2],6))
c = '{0:.3f}'.format(round(h[i],3))
f.write('{0:<4}{1:>16}{2:>16}{3:>16}{4:>10}\n'.format(s,x,y,z,c))
f.close()
In the LCAO mode of GPAW one can also get the Mulliken charges. Test before using:
def xyzm(atoms, filename):
from gpaw.lcao.tools import get_mulliken
f = open('{0}.xyz'.format(filename), 'w') # set xmol format
a = atoms.get_chemical_symbols() # get chemical symbols
p = atoms.get_positions() # get positions of the atoms
m = get_mulliken(atoms.calc, range(len(a)))
f.write('{0}\n'.format(len(a)))
f.write('Properties=species:S:1:pos:R:3:charge:R:1\n') # ensure compatibility with ASE
for i in range(len(a)): # print symbol, positions, and charge
s = '{0}'.format(a[i])
x = '{0:.6f}'.format(round(p[i][0],6))
y = '{0:.6f}'.format(round(p[i][1],6))
z = '{0:.6f}'.format(round(p[i][2],6))
c = '{0:.3f}'.format(round(m[i],3))
f.write('{0:<4}{1:>16}{2:>16}{3:>16}{4:>10}\n'.format(s,x,y,z,c))
f.close()
Archives
Categories
My work was supported by the Estonian Research Council under grants PUT1107, PRG259 and STP52. My research was supported by the from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101031656. All related posts are tagged with MSCA.