Geometry optimization¶
Modules: pyscf.geomopt
Basics¶
PySCF implements geometry optimization via interfaces to geomeTRIC and PyBerny, and through PySCF extension qsdopt <https://github.com/pyscf/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)
print(mol_eq.atom_coords())
# pyberny
from pyscf.geomopt.berny_solver import optimize
mol_eq = optimize(mf, maxsteps=100)
print(mol_eq.atom_coords())
The second way is to create an optimizer() from the Gradients class:
# geometric
mol_eq = mf.Gradients().optimizer(solver='geomeTRIC').kernel()
print(mol_eq.atom_coords())
# pyberny
mol_eq = mf.Gradients().optimizer(solver='berny').kernel()
print(mol_eq.atom_coords())
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)
Constraints¶
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
berny_solver.
'''
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(mol.atom_coords())
print('New geometry (Bohr)')
print(new_mol.atom_coords())
#
# 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")
optimizer.kernel()
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"
mol.build()
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
#geomopt.optimize(excited_grad)
For examples of state-specific and state-averaged CASSCF geometry optimizations, see examples/geomopt/12-mcscf_excited_states.py.
Callback¶
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
''',
basis='3-21g')
mf = scf.RHF(mol)
# Run analyze function in callback
def cb(envs):
mf = envs['g_scanner'].base
mf.analyze(verbose=4)
#
# 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
opt.kernel()