24. Miscellaneous

24.1. Decoration pipe

24.1.1. SCF

There are three decoration function for Hartree-Fock class density_fit(), sfx2c(), newton() to apply density fitting, scalar relativistic correction and second order SCF. The different ordering of the three decoration operations have different effects. For example

#!/usr/bin/env python
#
# Author: Qiming Sun <osirpt.sun@gmail.com>
#

import numpy
from pyscf import gto
from pyscf import scf

'''
Mixing decoration, for density fitting, scalar relativistic effects, and
second order (Newton-Raphson) SCF.

Density fitting and scalar relativistic effects can be applied together,
regardless to the order you apply the decoration.

NOTE the second order SCF (New in version 1.1) decorating operation are not
commutable with scf.density_fit operation
        [scf.density_fit, scf.sfx2c1e    ] == 0
        [scf.newton     , scf.sfx2c1e    ] != 0
        [scf.newton     , scf.density_fit] != 0
* scf.density_fit(scf.newton(scf.RHF(mol))) is the SOSCF for regular 2e
  integrals, but with density fitting integrals for the Hessian.  It's an
  approximate SOSCF optimization method;
* scf.newton(scf.density_fit(scf.RHF(mol))) is the exact second order
  optimization for the given scf object which is a density-fitted-scf method.
  The SOSCF is not an approximate scheme.
* scf.density_fit(scf.newton(scf.density_fit(scf.RHF(mol))), auxbasis='ahlrichs')
  is an approximate SOSCF scheme for the given density-fitted-scf method.
  Here we use small density fitting basis (ahlrichs cfit basis) to approximate
  the Hessian for the large-basis-density-fitted-scf scheme.
'''

mol = gto.Mole()
mol.build(
    verbose = 0,
    atom = '''8  0  0.     0
              1  0  -0.757 0.587
              1  0  0.757  0.587''',
    basis = 'ccpvdz',
)

#
# 1. spin-free X2C-HF with density fitting approximation on 2E integrals
#
mf = scf.density_fit(scf.sfx2c1e(scf.RHF(mol)))
mf = scf.RHF(mol).x2c().density_fit()  # Stream style
energy = mf.kernel()
print('E = %.12f, ref = -76.075408156180' % energy)

#
# 2. spin-free X2C correction for density-fitting HF.  Since X2C correction is
# commutable with density fitting operation, it is fully equivalent to case 1.
#
mf = scf.sfx2c1e(scf.density_fit(scf.RHF(mol)))
mf = scf.RHF(mol).density_fit().x2c()  # Stream style
energy = mf.kernel()
print('E = %.12f, ref = -76.075408156180' % energy)

#
# 3. The order to apply X2C or newton method matters.  If relativistic effects
# need to be included in the calculation, you should first call x2c then
# newton method.
#
e1 = scf.RHF(mol).kernel()
e2 = scf.RHF(mol).x2c().kernel()
print('E(NR)         = %.12f  E(X2C)        = %.12f' % (e1, e2))
e1 = scf.RHF(mol).newton().x2c().kernel()
e2 = scf.RHF(mol).x2c().newton().kernel()
print('E(Newton,X2C) = %.12f  E(X2C,Newton) = %.12f' % (e1, e2))

#
# 4. Newton method for non-relativistic HF
#
mf = scf.newton(scf.RHF(mol))
mf = scf.RHF(mol).newton()  # Stream style
energy = mf.kernel()
print('E = %.12f, ref = -76.026765673120' % energy)

#
# 5. Newton method for non-relativistic HF with density fitting for orbital
# hessian of newton solver.  Note the answer is equal to case 3, but the
# solver "mf" is different.
#
mf = scf.density_fit(scf.newton(scf.RHF(mol)))
mf = scf.RHF(mol).newton().density_fit()
energy = mf.kernel()
print('E = %.12f, ref = -76.026765673120' % energy)

#
# 6. Newton method to solve the density-fitting approximated HF object.  There
# is no approximation for newton method (orbital hessian).  Note the density
# fitting is applied on HF object only.  It does not affect the Newton solver.
#
mf = scf.newton(scf.density_fit(scf.RHF(mol)))
mf = scf.RHF(mol).density_fit().newton()
energy = mf.kernel()
print('E = %.12f, ref = -76.026744737357' % energy)

#
# 7. Newton method for density-fitting HF, and the hessian of Newton solver is
# also approximated with density fitting.  Note the anwser is equivalent to
# case 6, but the solver "mf" is different.  Here the fitting basis for HF and
# Newton solver are different.  HF is approximated with the default density
# fitting basis (Weigend cfit basis).  Newton solver is approximated with
# Ahlrichs cfit basis.
#
mf = scf.density_fit(scf.newton(scf.density_fit(scf.RHF(mol))), 'ahlrichs')
mf = scf.RHF(mol).density_fit().newton().density_fit(auxbasis='ahlrichs')
energy = mf.kernel()
print('E = %.12f, ref = -76.026744737357' % energy)

24.1.2. FCI

Direct FCI solver cannot guarantee the CI wave function to be the spin eigenfunction. Decoration function fci.addons.fix_spin_() can fix this issue.

24.1.3. CASSCF

mcscf.density_fit(), and scf.sfx2c() can be used to decorate CASSCF/CASCI class. Like the ordering problem in SCF decoration operation, the density fitting for CASSCF solver only affect the CASSCF optimization procedure. It does not change the 2e integrals for CASSCF Hamiltonian. For example

#!/usr/bin/env python
#
# Author: Qiming Sun <osirpt.sun@gmail.com>
#

from pyscf import gto, scf, mcscf

'''
Density fitting for orbital optimzation.

Note mcscf.density_fit function follows the same convention of decoration
ordering which is applied in the SCF decoration.  See pyscf/mcscf/df.py for
more details and pyscf/example/scf/23-decorate_scf.py as an exmple.
'''

mol = gto.Mole()
mol.build(
    atom = [
    ["C", (-0.65830719,  0.61123287, -0.00800148)],
    ["C", ( 0.73685281,  0.61123287, -0.00800148)],
    ["C", ( 1.43439081,  1.81898387, -0.00800148)],
    ["C", ( 0.73673681,  3.02749287, -0.00920048)],
    ["C", (-0.65808819,  3.02741487, -0.00967948)],
    ["C", (-1.35568919,  1.81920887, -0.00868348)],
    ["H", (-1.20806619, -0.34108413, -0.00755148)],
    ["H", ( 1.28636081, -0.34128013, -0.00668648)],
    ["H", ( 2.53407081,  1.81906387, -0.00736748)],
    ["H", ( 1.28693681,  3.97963587, -0.00925948)],
    ["H", (-1.20821019,  3.97969587, -0.01063248)],
    ["H", (-2.45529319,  1.81939187, -0.00886348)],],
    basis = 'ccpvtz'
)

mf = scf.RHF(mol)
mf.conv_tol = 1e-8
e = mf.kernel()

#
# DFCASSCF uses density-fitting 2e integrals overall, regardless the
# underlying mean-filed object
#
mc = mcscf.DFCASSCF(mf, 6, 6)
mo = mc.sort_mo([17,20,21,22,23,30])
mc.kernel(mo)
print('E(CAS) = %.12f, ref = -230.845892901370' % mc.e_tot)

#
# Assign DF basis
#
mc = mcscf.DFCASSCF(mf, 6, 6, auxbasis='ccpvtzfit')
mo = mc.sort_mo([17,20,21,22,23,30])
mc.kernel(mo)
print('E(CAS) = %.12f, ref = -230.845892901370' % mc.e_tot)

24.2. Customizing Hamiltonian

PySCF supports user-defined Hamiltonian for many modules. To customize Hamiltonian for Hartree-Fock, CASSCF, MP2, CCSD, etc, one need to replace the methods get_hcore(), get_ovlp() and attribute _eri of SCF class for new Hamiltonian. E.g. the user-defined Hamiltonian for Hartree-Fock

#!/usr/bin/env python
#
# Author: Qiming Sun <osirpt.sun@gmail.com>
#

import numpy
from pyscf import gto, scf, ao2mo

'''
Customizing Hamiltonian for SCF module.

Three steps to define Hamiltonian for SCF:
1. Specify the number of electrons. (Note mole object must be "built" before doing this step)
2. Overwrite three attributes of scf object
    .get_hcore
    .get_ovlp
    ._eri
3. Specify initial guess (to overwrite the default atomic density initial guess)

Note you will see warning message on the screen:

        Overwritten attributes  get_ovlp get_hcore  of <class 'pyscf.scf.hf.RHF'>

'''

mol = gto.M()
n = 10
mol.nelectron = n

mf = scf.RHF(mol)
h1 = numpy.zeros((n,n))
for i in range(n-1):
    h1[i,i+1] = h1[i+1,i] = -1.0
h1[n-1,0] = h1[0,n-1] = -1.0  # PBC
eri = numpy.zeros((n,n,n,n))
for i in range(n):
    eri[i,i,i,i] = 4.0

mf.get_hcore = lambda *args: h1
mf.get_ovlp = lambda *args: numpy.eye(n)
# ao2mo.restore(8, eri, n) to get 8-fold permutation symmetry of the integrals
# ._eri only supports the two-electron integrals in 4-fold or 8-fold symmetry.
mf._eri = ao2mo.restore(8, eri, n)

mf.kernel()

# If you need to run post-HF calculations based on the customized Hamiltonian,
# setting incore_anyway=True to ensure the customized Hamiltonian (the _eri
# attribute) to be used.  Without this parameter, some post-HF method
# (particularly in the MO integral transformation) may ignore the customized
# Hamiltonian if memory is not enough.
mol.incore_anyway = True

and the user-defined Hamiltonian for CASSCF

#!/usr/bin/env python
#
# Author: Qiming Sun <osirpt.sun@gmail.com>
#

import numpy
from pyscf import gto, scf, ao2mo, mcscf

'''
User-defined Hamiltonian for CASSCF module.

Defining Hamiltonian once for SCF object, the derivate post-HF method get the
Hamiltonian automatically.
'''

mol = gto.M()
mol.nelectron = 6

# incore_anyway=True ensures the customized Hamiltonian (the _eri attribute)
# to be used.  Without this parameter, the MO integral transformation may
# ignore the customized Hamiltonian if memory is not enough.
mol.incore_anyway = True

#
# 1D anti-PBC Hubbard model at half filling
#
n = 12

h1 = numpy.zeros((n,n))
for i in range(n-1):
    h1[i,i+1] = h1[i+1,i] = -1.0
h1[n-1,0] = h1[0,n-1] = -1.0
eri = numpy.zeros((n,n,n,n))
for i in range(n):
    eri[i,i,i,i] = 2.0

mf = scf.RHF(mol)
mf.get_hcore = lambda *args: h1
mf.get_ovlp = lambda *args: numpy.eye(n)
mf._eri = ao2mo.restore(8, eri, n)
mf.init_guess = '1e'
mf.kernel()

mycas = mcscf.CASSCF(mf, 4, 4)
mycas.kernel()