Biased Gradient Squared Descent

Biased Gradient Squared Descent is a saddle point finding method that does not require knowledge of a product state. The method converts a potential energy surface (PES) into the square of the gradient which converts all critical points into local minima. A bias term is added to stabilize critical points around a specified energy, beta, and destablize all other critical points leading to the following family of objective functions:

\(H(r;\alpha,\beta) = |\nabla V(r)| ^ {2} + \alpha (V(r) - \beta)^{2}\) (1)

Implementation details of BGSD are below with citation. 1

BGSD Potential

This class creates an atoms calculator (See ASE for details), which converts a different atoms object describing the PES, pot, into the H landscape. Below defines the various options for this class:

class tsase.bgsd. BGSD_potential (pot, alpha=1.0, beta=1.0, dr=1e-5):

pot : atoms object defining the PES and the reactant state

alpha : parameter \(\alpha\) in equation (1)

beta : parameter \(\beta\) in equation (1)

dr : finite difference size for gradient of \(|\nabla V(r)|^{2}\)

NOTE: The initial configuration of pot needs to be at the local minimum of the reactant state. The parameter, \(\beta\), is a relative energy to the local minimum.

BGSD

This class follows the procedure described in 1 to find a first order saddle point connected to a reactant state using BGSD. First the isosurface corresponding to the value of \(\beta\) used is sampled. There options for various local displacement schemes while sampling this isosurface. Next, \(H\) and the \(|\nabla V(r)|^{2}\) are minimized using the SDLBFGS optimizer. Finally the class tests to see if the critical point found is connected to the reactant state and a 1st order saddle point. Below defines the various options for this class:

class tsase.bgsd. BGSD (p, alpha=5.0, beta=0.0, dr=1e-5, numstep=1000, stepsize=0.05, kT_MC=0.01, k=5.0, displace_atomlist=[], displace_radius=3.3, displace_all_listed=True, CC1=0.01, CC2=0.001*0.001, CC3 = 0.0001, maxstepsize=0.2, maxnumstep=1000, memory=100, eigcheck = ‘dimer’)

BGSD:

p, alpha, beta, dr : these parameters are exactly the same as the parameters defined in the BGSD Potential class

NOTE : The initial configuration of pot needs to be at the local minimum of the reactant state. The parameter, \(\beta\), is a relative energy to the local minimum.

Monte Carlo (MC) Sampling of Isosurface:

numstep : number of steps MC steps

stepsize : MC max per atom stepsize

kT_MC : temperature in kT for MC sampling

k : spring constant of isosurface

Biasing Initial Configuration:

displacement_atomlist : index number of atoms to be moved in sampling of isocontour surface

displacement_radius : displace atoms within this specified radius from atoms in atomlist

displace_all_listedIf True then atomlist uses all atoms listed. If False, the one atom is randomly choosen one form the atomlist.

(Note: once the one atom is selected, all atoms within the displacement_radius of this atom are moved)

SDLBFGS optimizer options:

CC1 : magnitude of the L2 norm used as the convergence criteria

CC2 : magnitude of \(|\nabla V(r)|^{2}\) used as convergence criteria for the \(|\nabla V(r)|^{2}\) landscape

CC3 : magnitude of the L2 norm used as the convergence criteria of \(|\nabla V(r)|^{2}\); this additional convergence critieria is used to look for inflection points

maxstepsize : maximum step size for the SDLBFGS optimizer

maxnumstep : maximum number of steps in optimization

memory : number of memory terms in SDLBFGS

Options for Determining if a Saddle Point is Connected and a First Order:

eigcheck : option to use hessian or dimer to determine the minimum mode

Example Script

Below is an example python script of how the two classes described above can be used.

#!/usr/bin/env python

import tsase

######## load V #########################
al = tsase.calculators.al()
p = tsase.io.read_con('al.con')
p.set_calculator(al)

#### minimize p so the system is at a local minimum ##########
bmin = tsase.optimize.SDLBFGS(p)
bmin.run()

##### set calculator to be H landscape using the PES of p and the BGSD_potential #######
p1 = tsase.io.read_con('al.con')
Hpotential = tsase.bgsd.BGSD_potential(pot = p,alpha = 5.0,beta = 0.2 ,dr=1.e-5)
p1.set_calculator(Hpotential)

####### run BGSD ##########################
bgsd = tsase.bgsd.BGSD(p,alpha=5.0,beta=0.27,displace_atomlist=[25,26,27],displace_radius=0.0)

"""

The function 'find saddle' samples the beta isosurface, minimizes H and the gradient squared landscape, 
and checks to see if a 1st order saddle point connected to the reactant state is found.  
It returns two parameters. The first parameter is True if a first order connected saddle point is found and 
False otherwise.  The second parameters is the number force calls required. 
The script produces a file of called 'SP.con' which is the critical point found in the optimization.  

"""

converged, ForceCalls = bgsd.find_saddle()

###########################################
print('Force Calls:',ForceCalls)
print('Was first order saddle point found',converged)

References

1(1,2)

10. Duncan, Q. Wu, K. Promislow, and G. Henkelman, “Biased gradient squared descent saddle point finding method”, J. Chem. Phys. 140, 194102 (2014). DOI