Geometry optimization

Modules: pyscf.geomopt


PySCF implements geometry optimization via interfaces to geomeTRIC and PyBerny, and through PySCF extension qsdopt <>`_(see :ref:`installing for installation instructions).

There are two ways to invoke geometry optimization. The first is to import the optimize() function from the respective modules, i.e., pyscf.geomopt.geometric_solver and pyscf.geomopt.berny_solver:

from pyscf import gto, scf
mol = gto.M(atom='N 0 0 0; N 0 0 1.2', basis='ccpvdz')
mf = scf.RHF(mol)

# geometric
from pyscf.geomopt.geometric_solver import optimize
mol_eq = optimize(mf, maxsteps=100)

# pyberny
from pyscf.geomopt.berny_solver import optimize
mol_eq = optimize(mf, maxsteps=100)

The second way is to create an optimizer() from the Gradients class:

# geometric
mol_eq = mf.Gradients().optimizer(solver='geomeTRIC').kernel()

# pyberny
mol_eq = mf.Gradients().optimizer(solver='berny').kernel()

For the geomeTRIC backend, the convergence criteria are controlled by the following parameters:

conv_params = { # These are the default settings
    'convergence_energy': 1e-6,  # Eh
    'convergence_grms': 3e-4,    # Eh/Bohr
    'convergence_gmax': 4.5e-4,  # Eh/Bohr
    'convergence_drms': 1.2e-3,  # Angstrom
    'convergence_dmax': 1.8e-3,  # Angstrom
mol_eq = optimize(mf, **conv_params)
mol_eq = mf.Gradients().optimizer(solver='geomeTRIC').kernel(conv_params)

For the PyBerny backend, the convergence criteria are controlled by the following parameters:

conv_params = {  # These are the default settings
    'gradientmax': 0.45e-3,  # Eh/[Bohr|rad]
    'gradientrms': 0.15e-3,  # Eh/[Bohr|rad]
    'stepmax': 1.8e-3,       # [Bohr|rad]
    'steprms': 1.2e-3,       # [Bohr|rad]
mol_eq = optimize(mf, **conv_params)
mol_eq = mf.Gradients().optimizer(solver='berny').kernel(conv_params)


geomeTRIC supports user defined constraints. The constraints can be specified in a text file. Its format can be found in the online documentation or in the template file. Then the name of the constraints file needs to be passed to PySCF:

params = {"constraints": "constraints.txt",}
mol_eq = optimize(mf, **params)
mol_eq = mf.Gradients().optimizer(solver='geomeTRIC').kernel(params)

The geometry optimization can also be carried out based on custom energy and gradient functions:

#!/usr/bin/env python

For the customized energy and gradients function (e.g. adding DFT-D3
correction), a fake pyscf method need to be created before passing to

import numpy as np
from pyscf import gto, scf
from pyscf.geomopt import berny_solver, geometric_solver, as_pyscf_method

mol = gto.M(atom='N 0 0 0; N 0 0 1.8', unit='Bohr', basis='ccpvdz')
mf = scf.RHF(mol)

grad_scan = scf.RHF(mol).nuc_grad_method().as_scanner()
def f(mol):
    e, g = grad_scan(mol)
    r = mol.atom_coords()
    penalty = np.linalg.norm(r[0] - r[1])**2 * 0.1
    e += penalty
    g[0] += (r[0] - r[1]) * 2 * 0.1
    g[1] -= (r[0] - r[1]) * 2 * 0.1
    print('Customized |g|', np.linalg.norm(g))
    return e, g

# Function as_pyscf_method is a wrapper that convert the "energy-gradients"
# function to berny_solver.  The "energy-gradients" function takes the Mole
# object as geometry input, and returns the energy and gradients of that
# geometry.
fake_method = as_pyscf_method(mol, f)
new_mol = berny_solver.optimize(fake_method)

print('Old geometry (Bohr)')

print('New geometry (Bohr)')

# Geometry can be also optimized with geomeTRIC library
new_mol = geometric_solver.optimize(fake_method)

Transition state optimization

Transition state optimization are available in geomeTRIC and qsdopt.

The PySCF extension qsdopt performs transition state optimizations through quadratic steepest descent method This is a second order method that requires computation of the hessian at some steps during the optimization process. The following is a minimal usage example:

from pyscf imort gto, scf
from pyscf.qsdopt.qsd_optimizer import QSD

mol = gto.M(atom='''
O 0 0 0
H 0 0 1.2
H 0, 0.5, -1.2''',
basis='minao', verbose=0, unit="Bohr")
mf = scf.RHF(mol)

optimizer = QSD(mf, stationary_point="TS")

Several keyword arguments can be passed to kernel:

  • hess_update_freq: Frequency for numerical reevaluation of hessian. = 0 evaluates the numerical hessian in the first iteration and is updated with an BFGS rule, unless approaching a trap region, where it is reevaluated. (Default: 0)

  • numhess_method: Method for evaluating numerical hessian. Forward and central differences are available with “forward” and “central”, respectively. (Default: “forward”)

  • max_iter: Maximum number of optimization steps. (Default: 100)

  • step: Maximum norm between two optimization steps. (Default: 0.1)

  • hmin: Minimum distance between current geometry and stationary point of the quadratic form to consider convergence reached. (Default: 1e-6)

  • gthres: Gradient norm to consider convergence. (Default: 1e-5)

from pyscf.qsdopt.qsd_optimizer import QSD
from pyscf import gto, scf

mol = gto.M(atom='''
C  -2.17177  -0.46068  -0.00000
C  -0.02194  -0.50455  0.00000
H  0.35133   0.30506   -0.59169
H  0.33677   -0.40796  1.00344
H  0.31586   -1.43259  -0.41175
H  -2.68285  -0.25342  0.93134
H  -2.40018  0.15510   -0.90287
H   -2.53568   -1.52338   -0.17790''', basis='minao', verbose=0)
mf = scf.RHF(mol)

optimizer = QSD(mf, stationary_point="TS")
optimizer.kernel(hess_update_freq=0, step=0.5, hmin=1e-3)

When using geomeTRIC to search transition state, you can simply add the keyword transition in geomeTRIC input configuration to trigger the TS search module:

params = {'transition': True}
mol_ts = mf.Gradients().optimizer(solver='geomeTRIC').kernel(params)

params = {'transition': True, 'trust': 0.02, 'tmax': 0.06}
mol_ts = mf.Gradients().optimizer(solver='geomeTRIC').kernel(params)

Keywords trust and tmax control the trust region in TS search. They are not mandatory input settings. The default settings in geomeTRIC library are chosen to maximize job success rates. Tuning trust region may bring small improvements to the TS optimization performance. For more detailed discussions of available options, we refer to the online documentation of geomeTRIC transition states.

Except the initial step, geomeTRIC library does not support the ability to read analytical Hessian from quantum chemistry calculations at runtime. The initial Hessian can be fed to the optimizer after enabling the keyword hessian:

params = {'transition': True, 'hessian': True}
mol_ts = mf.Gradients().optimizer(solver='geomeTRIC').kernel(params)

The interface implemented in PySCF can generate a temporary file to save the initial Hessian and pass to geomeTRIC. This keyword will be ignored if the analytical Hessian of the underlying method was not implemented.

Excited states

For excited-state geometry optimizations, the state to be optimized needs to be specified in the respective Gradients objects:

#!/usr/bin/env python

Optimize the geometry of the excited states

Note when optiming the excited states, states may flip and this may cause
convergence issue in geometry optimizer.

from pyscf import gto
from pyscf import scf
from pyscf import ci, tdscf, mcscf
from pyscf import geomopt

mol = gto.Mole()
mol.atom="N; N 1, 1.1"
mol.basis= "6-31g"
mol1 = mol.copy()

mf = scf.RHF(mol).run()

mc = mcscf.CASCI(mf, 4,4)
mc.fcisolver.nstates = 3
excited_grad = mc.nuc_grad_method().as_scanner(state=2)
mol1 = excited_grad.optimizer().kernel()
# or
#mol1 = geomopt.optimize(excited_grad)

td = tdscf.TDHF(mf)
td.nstates = 5
excited_grad = td.nuc_grad_method().as_scanner(state=4)
mol1 = excited_grad.optimizer().kernel()
# or
#mol1 = geomopt.optimize(excited_grad)

myci = ci.CISD(mf)
myci.nstates = 2
excited_grad = myci.nuc_grad_method().as_scanner(state=1)
mol1 = excited_grad.optimizer().kernel()
# or

For examples of state-specific and state-averaged CASSCF geometry optimizations, see examples/geomopt/


Callback functions can be invoked at each optimization step. The following example shows how to add charge analysis for each geometry during the optimization.

#!/usr/bin/env python

Optimize molecular geometry within the environment of QM/MM charges.

from pyscf import gto, scf
from pyscf.geomopt import berny_solver
from pyscf.geomopt import geometric_solver

mol = gto.M(atom='''
C        0.000000    0.000000             -0.542500
O        0.000000    0.000000              0.677500
H        0.000000    0.9353074360871938   -1.082500
H        0.000000   -0.9353074360871938   -1.082500

mf = scf.RHF(mol)

# Run analyze function in callback
def cb(envs):
    mf = envs['g_scanner'].base

# Method 1: Pass callback to optimize function
geometric_solver.optimize(mf, callback=cb)

berny_solver.optimize(mf, callback=cb)

# Method 2: Add callback to geometry optimizer
opt = mf.nuc_grad_method().as_scanner().optimizer()
opt.callback = cb