Geometry optimization

Modules: geomopt

Basics

PySCF implements geometry optimization via interfaces to geomeTRIC and PyBerny (see Installation 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 with the format described here. One needs to pass the name of this file 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)

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()