Configuration interaction (CISD and FCI)

Modules: pyscf.ci, pyscf.fci, pyscf.pbc.ci

PySCF has separate molecular implementations for configuration interaction singles and doubles (CISD) (ci) and full configuration interaction (FCI) (fci). The functionalities of the CISD implementation are similar to the functionalities of MP2 (Second-order Møller–Plesset perturbation theory (MP2)) and CCSD (Coupled-cluster theory) in PySCF.

Introduction

Configuration interaction (see e.g. [31, 32, 33]) is a post-Hartree–Fock method, diagonalising the many-electron Hamiltonian matrix. FCI includes all Slater determinants of appropriate symmetry in the eigenvector basis. CISD includes only those that differ from the Hartree–Fock determinant by at most two occupied spinorbitals.

Configuration interaction singles and doubles (CISD)

A simple example (see examples/ci/00-simple_cisd.py) of running a restricted and an unrestricted CISD calculation is

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

'''
A simple example to run CISD calculation.
'''

import pyscf

mol = pyscf.M(
    atom = 'H 0 0 0; F 0 0 1.1',
    basis = 'ccpvdz')

mf = mol.HF().run()
mycc = mf.CISD().run()
print('RCISD correlation energy', mycc.e_corr)

mf = mol.UHF().run()
mycc = mf.CISD().run()
print('UCISD correlation energy', mycc.e_corr)

which outputs

converged SCF energy = -99.9873974403488
E(RCISD) = -100.196182975762  E_corr = -0.2087855354132202
RCISD correlation energy -0.2087855354132202

namely, the (restricted) Hartree–Fock energy, the RCISD energy and their difference, the RCISD correlation energy, as well as

converged SCF energy = -99.9873974403486  <S^2> = 5.3734794e-13  2S+1 = 1
E(UCISD) = -100.196182975774  E_corr = -0.2087855354254317
UCISD correlation energy -0.20878553542543174

namely, the (unrestricted) Hartree–Fock energy, the UCISD energy and their difference, the UCISD correlation energy.

Note

mycc = mf.CISD().run()

could have been replaced by

myci = pyscf.ci.CISD(mf)
myci.kernel()

or

myci = pyscf.ci.CISD(mf).run()

for the same result.

Note

Density fitting is not implemented for CISD.

Spin symmetry

The CISD module in PySCF supports a number of reference wavefunctions with broken spin symmetry. In particular, CISD can be performed with a spin-restricted, spin-unrestricted, and general (spin-mixed) Hartree-Fock solution, leading to the RCISD, UCISD, and GCISD methods. The small example above shows this for RCISD and UCISD.

The module-level ci.CISD(mf) constructor can infer the correct method based on the level of symmetry-breaking in the mean-field argument. For more explicit control or inspection, the respective classes and functions can be found in cisd.py (restricted), ucisd.py (unrestricted), and gcisd.py (general).

For example, a spin-unrestricted calculation on triplet oxygen can be performed as follows:

from pyscf import gto, scf, ci
mol = gto.M(
    atom = 'O 0 0 0; O 0 0 1.2',  # in Angstrom
    basis = 'ccpvdz',
    spin = 2
)
mf = scf.HF(mol).run() # this is UHF
myci = ci.CISD(mf).run() # this is UCISD
print('UCISD total energy = ', myci.e_tot)

Properties

A number of properties are available at the CISD level.

Unrelaxed 1- and 2-electron reduced density matrices can be calculated. They are returned in the MO basis:

dm1 = myci.make_rdm1()
dm2 = myci.make_rdm2()

Analytical nuclear gradients can be calculated (see [34, 35] and references therein):

mygrad = myci.nuc_grad_method().run()

Frozen orbitals

By default, CISD calculations in PySCF correlate all electrons in all available orbitals. To freeze the lowest-energy core orbitals, use the frozen keyword argument:

myci = ci.CISD(mf, frozen=2).run()

To freeze occupied and/or unoccupied orbitals with finer control, a list of 0-based orbital indices can be provided as the frozen keyword argument:

# freeze 2 core orbitals
myci = ci.CISD(mf, frozen=[0,1]).run()
# freeze 2 core orbitals and 3 unoccupied orbitals
myci = ci.CISD(mf, frozen=[0,1,16,17,18]).run()

Wavefunction overlap

The following example shows how to evaluate the overlap of two different CISD wavefunctions.

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

'''
Overlap of two CISD wave functions (they can be obtained from different
geometries).
'''

from functools import reduce
import numpy
from pyscf import gto, scf, ci

#
# RCISD wavefunction overlap
#
myhf1 = gto.M(atom='H 0 0 0; F 0 0 1.1', basis='6-31g', verbose=0).apply(scf.RHF).run()
ci1 = ci.CISD(myhf1).run()
print('CISD energy of mol1', ci1.e_tot)

myhf2 = gto.M(atom='H 0 0 0; F 0 0 1.2', basis='6-31g', verbose=0).apply(scf.RHF).run()
ci2 = ci.CISD(myhf2).run()
print('CISD energy of mol2', ci2.e_tot)

s12 = gto.intor_cross('cint1e_ovlp_sph', myhf1.mol, myhf2.mol)
s12 = reduce(numpy.dot, (myhf1.mo_coeff.T, s12, myhf2.mo_coeff))
nmo = myhf2.mo_energy.size
nocc = myhf2.mol.nelectron // 2
print('<CISD-mol1|CISD-mol2> = ', ci.cisd.overlap(ci1.ci, ci2.ci, nmo, nocc, s12))

#
# UCISD wavefunction overlap
#
myhf1 = gto.M(atom='H 0 0 0; F 0 0 1.1', basis='6-31g', spin=2, verbose=0).apply(scf.UHF).run()
ci1 = ci.UCISD(myhf1).run()
print('CISD energy of mol1', ci1.e_tot)

myhf2 = gto.M(atom='H 0 0 0; F 0 0 1.2', basis='6-31g', spin=2, verbose=0).apply(scf.UHF).run()
ci2 = ci.UCISD(myhf2).run()
print('CISD energy of mol2', ci2.e_tot)

s12 = gto.intor_cross('cint1e_ovlp_sph', myhf1.mol, myhf2.mol)
mo1a, mo1b = myhf1.mo_coeff
mo2a, mo2b = myhf2.mo_coeff
s12 = (reduce(numpy.dot, (mo1a.T, s12, mo2a)),
       reduce(numpy.dot, (mo1b.T, s12, mo2b)))
nmo = (myhf2.mo_energy[0].size, myhf2.mo_energy[1].size)
nocc = myhf2.nelec
print('<CISD-mol1|CISD-mol2> = ', ci.ucisd.overlap(ci1.ci, ci2.ci, nmo, nocc, s12))

Full configuration interaction (FCI)

A simple example (see examples/fci/00-simple-fci.py) of running a restricted and an unrestricted FCI calculation is:

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

'''
A simple example to run FCI
'''

import pyscf

mol = pyscf.M(
    atom = 'H 0 0 0; F 0 0 1.1',  # in Angstrom
    basis = '6-31g',
    symmetry = True,
)
myhf = mol.RHF().run()

#
# create an FCI solver based on the SCF object
#
cisolver = pyscf.fci.FCI(myhf)
print('E(FCI) = %.12f' % cisolver.kernel()[0])

#
# create an FCI solver based on the SCF object
#
myuhf = mol.UHF().run()
cisolver = pyscf.fci.FCI(myuhf)
print('E(UHF-FCI) = %.12f' % cisolver.kernel()[0])

#
# create an FCI solver based on the given orbitals and the num. electrons and
# spin of the mol object
#
cisolver = pyscf.fci.FCI(mol, myhf.mo_coeff)
print('E(FCI) = %.12f' % cisolver.kernel()[0])

The second and third FCI calculation in the example show two different ways of running unrestricted FCI calculations.

Setting Spin and Symmetry of Wavefunction

The following example discusses how to set the spin for an FCI calculations:

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

'''
Assign spin state for FCI wavefunction.

By default, the FCI solver will take Mole attribute spin for the spin state.
It can be overwritten by passing kwarg ``nelec`` to the kernel function of FCI
solver.  The nelec argument is a two-element tuple.  The first is the number
of alpha electrons; the second is the number of beta electrons.

If spin-contamination is observed on FCI wavefunction, we can use the
decoration function :func:`fci.addons.fix_spin_` to level shift the energy of
states which do not have the target spin.
'''

import numpy
from pyscf import gto, scf, fci

mol = gto.M(atom='Ne 0 0 0', basis='631g', spin=2)
m = scf.RHF(mol)
m.kernel()
norb = m.mo_energy.size

fs = fci.FCI(mol, m.mo_coeff)
e, c = fs.kernel()
print('E = %.12f  2S+1 = %.7f' %
      (e, fs.spin_square(c, norb, (6,4))[1]))

e, c = fs.kernel(nelec=(5,5))
print('E = %.12f  2S+1 = %.7f' %
      (e, fs.spin_square(c, norb, (5,5))[1]))


fs = fci.addons.fix_spin_(fci.FCI(mol, m.mo_coeff), shift=.5)
e, c = fs.kernel()
print('E = %.12f  2S+1 = %.7f' %
      (e, fs.spin_square(c, norb, (6,4))[1]))


#
# Example 2:  Oxygen molecule singlet state
#

nelec = (8,8)
mol = gto.M(atom='O 0 0 0; O 0 0 1.2', spin=2, basis='sto3g',
            symmetry=1, verbose=0)
mf = scf.RHF(mol).run()
mci = fci.FCI(mol, mf.mo_coeff)
mci.wfnsym = 'A1g'
mci = fci.addons.fix_spin_(mci, ss=0)
# Use keyword argument nelec to explicitly control the spin. Otherwise
# mol.spin is applied.
e, civec = mci.kernel(nelec=nelec)
print('A1g singlet E = %.12f  2S+1 = %.7f' %
      (e, mci.spin_square(civec, mf.mo_coeff.shape[1], nelec)[1]))

mci.wfnsym = 'A2g'
mci = fci.addons.fix_spin_(mci, ss=0)
e, civec = mci.kernel(nelec=nelec)
print('A2g singlet E = %.12f  2S+1 = %.7f' %
      (e, mci.spin_square(civec, mf.mo_coeff.shape[1], nelec)[1]))

mol = gto.M(atom='O 0 0 0; O 0 0 1.2', spin=2, basis='sto3g',
            verbose=0)
mf = scf.RHF(mol).run()
mci = fci.FCI(mol, mf.mo_coeff)
mci = fci.addons.fix_spin_(mci, ss=0)
e, civec = mci.kernel(nelec=nelec)
print('Singlet E = %.12f  2S+1 = %.7f' %
      (e, mci.spin_square(civec, mf.mo_coeff.shape[1], nelec)[1]))


#
# Example 3:  Triplet and Sz=0
#
# It's common that the default initial guess has no overlap with the ground
# state. In this example, the default initial guess is singlet. Computing
# multiple roots to overcome the initial guess issue.
#
mol = gto.M(verbose=0,
            atom = '''
                  H     1  -1.      0
                  H     0  -1.     -1
                  H     0  -0.5    -0
                  H     0  -0.     -1
                  H     1  -0.5     0
                  H     0   1.      1''',
            basis='sto-3g')
mf = scf.RHF(mol).run()
mci = fci.FCI(mol, mf.mo_coeff, singlet=False)
mci = fci.addons.fix_spin_(mci, ss=2)

e, civec = mci.kernel(nroots=2)
nelec = (3,3)
print('Triplet E = %.12f  2S+1 = %.7f' %
      (e[0], mci.spin_square(civec[0], mf.mo_coeff.shape[1], nelec)[1]))


#
# Be careful with the trick of energy penalty.  Numerical problems may be
# observed when function fix_spin_ was applied.
#

mol = gto.M(atom='O 0 0 0', basis='6-31G')
m = scf.RHF(mol).run()
norb = m.mo_coeff.shape[1]
nelec = mol.nelec
fs = fci.addons.fix_spin_(fci.FCI(mol, m.mo_coeff), .5)
fs.nroots = 15
e, fcivec = fs.kernel(verbose=5)
# The first 5 states should be degenerated. The degeneracy may be broken.
for i, c in enumerate(fcivec):
    print('state = %d, E = %.9f, S^2=%.4f' %
          (i, e[i], fci.spin_op.spin_square(c, norb, nelec)[0]))

and the next example how to set symmetry:

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

'''
Assign FCI wavefunction symmetry
'''

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

mol = gto.M(atom='H 0 0 0; F 0 0 1.1', basis='631g', symmetry=True)
m = scf.RHF(mol).run()
norb = m.mo_coeff.shape[1]
nelec = mol.nelec

fs = fci.FCI(mol, m.mo_coeff)
fs.wfnsym = 'E1x'
e, c = fs.kernel(verbose=5)
print('Energy of %s state %.12f' % (fs.wfnsym, e))

#
# Symmetry adapted FCI solver can be called with arbitrary Hamiltonian.  In the
# following example, you need to prepare 1-electron and 2-electron integrals,
# core energy shift, and the symmetry of orbitals.
#
cas_idx = numpy.arange(2,10)
core_idx = numpy.arange(2)
mo_cas = m.mo_coeff[:,cas_idx]
mo_core = m.mo_coeff[:,core_idx]
dm_core = mo_core.dot(mo_core.T) * 2
vhf_core = m.get_veff(mol, dm_core)
h1 = mo_cas.T.dot( m.get_hcore() + vhf_core ).dot(mo_cas)
h2 = ao2mo.kernel(mol, mo_cas)
ecore = (numpy.einsum('ij,ji', dm_core, m.get_hcore())
         + .5 * numpy.einsum('ij,ji', dm_core, vhf_core) + m.energy_nuc())
orbsym = symm.label_orb_symm(mol, mol.irrep_id, mol.symm_orb, mo_cas)

norb = mo_cas.shape[1]
ncore = mo_core.shape[1]
nelec = mol.nelectron - ncore * 2
fs = fci.direct_spin0_symm.FCI(mol)
fs.wfnsym = 'A2'
e, c = fs.kernel(h1, h2, norb, nelec, ecore=ecore, orbsym=orbsym, verbose=5)
print('Energy of %s state %.12f' % (fs.wfnsym, e))

#
# mcscf module has a function h1e_for_cas to generate 1-electron Hamiltonian and
# core energy.
#
cas_idx = numpy.arange(2,10)
core_idx = numpy.arange(2)
mo_cas = m.mo_coeff[:,cas_idx]
mo_core = m.mo_coeff[:,core_idx]
norb = mo_cas.shape[1]
ncore = mo_core.shape[1]
nelec = mol.nelectron - ncore * 2
h1, ecore = mcscf.casci.h1e_for_cas(m, m.mo_coeff, norb, ncore)
h2 = ao2mo.kernel(mol, mo_cas)
orbsym = scf.hf_symm.get_orbsym(mol, m.mo_coeff)[cas_idx]
fs = fci.direct_spin0_symm.FCI(mol)
fs.wfnsym = 'A2'
e, c = fs.kernel(h1, h2, norb, nelec, ecore=ecore, orbsym=orbsym)
print('Energy of %s state %.12f' % (fs.wfnsym, e))

#
# Using the CASCI object, initialization of 1-electron and 2-electron integrals
# can be further simplified.
#
cas_idx = numpy.arange(2,10)
norb = len(cas_idx)
ncore = 2
nelec = mol.nelectron - ncore * 2
mc = mcscf.CASCI(m, norb, nelec)
mo = mc.sort_mo(cas_idx, base=0)
h1, ecore = mc.get_h1eff(mo)
h2 = mc.get_h2eff(mo)
orbsym = scf.hf_symm.get_orbsym(mol, m.mo_coeff)[cas_idx]
fs = fci.direct_spin0_symm.FCI(mol)
fs.wfnsym = 'A2'
e, c = fs.kernel(h1, h2, norb, nelec, ecore=ecore, orbsym=orbsym)
print('Energy of %s state %.12f' % (fs.wfnsym, e))


#
# In the following example, the default initial guess of the regular FCI
# solver is different to the symmetry of HF wfn.  In the symmetry adapted FCI
# solver, the initial guess of required wfnsym is generated.  The FCI solution
# has the same symmetry as the HF wfn. 
#
mol = gto.M(atom='C 0. 0. 0.; C 0. 0. 1.24253', basis='6-31g', symmetry=True)
mf = scf.RHF(mol).run()
cas_idx = [2, 3, 4, 5, 6]
norb = len(cas_idx)
ncore = 2
nelec = mol.nelectron - ncore * 2
mc = mcscf.CASCI(mf, norb, nelec)
h1, ecore = mc.get_h1eff(mf.mo_coeff)
h2 = mc.get_h2eff(mf.mo_coeff)
orbsym = scf.hf_symm.get_orbsym(mol, mf.mo_coeff)[cas_idx]

fs = fci.direct_spin0.FCI(mol)  # Regular FCI solver
fs.davidson_only = True
e, c = fs.kernel(h1, h2, norb, nelec, ecore=ecore)
print('''Using regular FCI solver, FCI converges to
  wfnsym = %s  E(FCI) = %.12f''' % (fci.addons.guess_wfnsym(c, len(cas_idx), nelec, orbsym), e))
hf_as_civec = numpy.zeros_like(c)
hf_as_civec[0,0] = 1
e, c = fs.kernel(h1, h2, norb, nelec, ecore=ecore, ci0=hf_as_civec)
print('''Using HF determinant as initial guess, FCI converges to
  wfnsym = %s  E(FCI) = %.12f''' % (fci.addons.guess_wfnsym(c, len(cas_idx), nelec, orbsym), e))

fs = fci.direct_spin0_symm.FCI(mol)  # Symmetry adapted FCI solver
fs.wfnsym = 0
e, c = fs.kernel(h1, h2, norb, nelec, ecore=ecore, orbsym=orbsym)
print('''Using symmetry adapted FCI solver, FCI converges to
  wfnsym = %s  E(FCI) = %.12f''' % (fci.addons.guess_wfnsym(c, len(cas_idx), nelec, orbsym), e))

Properties

Reduced density matrices can be evaluated, see the following example:

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

'''
Compute FCI 1,2,3,4-particle density matrices
'''

#
# Note: Environment variable LD_PRELOAD=...libmkl_def.so may cause this script
# crashing
#

import numpy
from pyscf import gto, scf, fci

mol = gto.Mole()
mol.build(
    atom = 'H 0 0 0; F 0 0 1.1',  # in Angstrom
    basis = '6-31g',
    spin = 2,
)
myhf = scf.RHF(mol)
myhf.kernel()

#
# First create FCI solver with function fci.FCI and solve the FCI problem
#
cisolver = fci.FCI(mol, myhf.mo_coeff)
e, fcivec = cisolver.kernel()

#
# Spin-traced 1-particle density matrix
#
norb = myhf.mo_coeff.shape[1]
# 6 alpha electrons, 4 beta electrons because spin = nelec_a-nelec_b = 2
nelec_a = 6
nelec_b = 4
dm1 = cisolver.make_rdm1(fcivec, norb, (nelec_a,nelec_b))

#
# alpha and beta 1-particle density matrices
#
dm1a, dm1b = cisolver.make_rdm1s(fcivec, norb, (nelec_a,nelec_b))
assert(numpy.allclose(dm1a+dm1b, dm1))

#
# Spin-traced 1 and 2-particle density matrices
#
dm1, dm2 = cisolver.make_rdm12(fcivec, norb, (nelec_a,nelec_b))

#
# alpha and beta 1-particle density matrices
# For 2-particle density matrix,
# dm2aa corresponds to alpha spin for both 1st electron and 2nd electron
# dm2ab corresponds to alpha spin for 1st electron and beta spin for 2nd electron
# dm2bb corresponds to beta spin for both 1st electron and 2nd electron
#
(dm1a, dm1b), (dm2aa,dm2ab,dm2bb) = cisolver.make_rdm12s(fcivec, norb, (nelec_a,nelec_b))
assert(numpy.allclose(dm2aa+dm2ab+dm2ab.transpose(2,3,0,1)+dm2bb, dm2))


#########################################
#
# Transition density matrices
#
#########################################

#
# First generate two CI vectors, of which the transition density matrices will
# be computed
#
cisolver.nroots = 2
(e0,e1), (fcivec0,fcivec1) = cisolver.kernel()

#
# Spin-traced 1-particle transition density matrix
# <0| p^+ q |1>
#
norb = myhf.mo_coeff.shape[1]
nelec_a = 6
nelec_b = 4
dm1 = cisolver.trans_rdm1(fcivec0, fcivec1, norb, (nelec_a,nelec_b))

#
# alpha and beta 1-particle transition density matrices
#
dm1a, dm1b = cisolver.trans_rdm1s(fcivec0, fcivec1, norb, (nelec_a,nelec_b))
assert(numpy.allclose(dm1a+dm1b, dm1))

#
# Spin-traced 1 and 2-particle transition density matrices
#
dm1, dm2 = cisolver.trans_rdm12(fcivec0, fcivec1, norb, (nelec_a,nelec_b))

#
# alpha and beta 1-particle transition density matrices
# For 2-particle density matrix,
# dm2aa corresponds to alpha spin for both 1st electron and 2nd electron
# dm2ab corresponds to alpha spin for 1st electron and beta spin for 2nd electron
# dm2ba corresponds to beta spin for 1st electron and alpha spin for 2nd electron
# dm2bb corresponds to beta spin for both 1st electron and 2nd electron
#
(dm1a, dm1b), (dm2aa, dm2ab, dm2ba, dm2bb) = \
        cisolver.trans_rdm12s(fcivec0, fcivec1, norb, (nelec_a,nelec_b))
assert(numpy.allclose(dm2aa+dm2ab+dm2ba+dm2bb, dm2))


#########################################
#
# 3 and 4-particle density matrices
#
#########################################

#
# Spin-traced 3-particle density matrix
# 1,2,3-pdm can be computed together.
# Note make_dm123 computes  dm3[p,q,r,s,t,u] = <p^+ q r^+ s t^+ u>  which is
# NOT the 3-particle density matrices.  Funciton reorder_dm123 transforms it
# to the true 3-particle DM  dm3[p,q,r,s,t,u] = <p^+ r^+ t^+ u s q> (as well
# as the 2-particle DM)
#
dm1, dm2, dm3 = fci.rdm.make_dm123('FCI3pdm_kern_sf', fcivec0, fcivec0, norb,
                                   (nelec_a,nelec_b))
dm1, dm2, dm3 = fci.rdm.reorder_dm123(dm1, dm2, dm3)

#
# Spin-traced 3-particle transition density matrix
#
dm1, dm2, dm3 = fci.rdm.make_dm123('FCI3pdm_kern_sf', fcivec0, fcivec1, norb,
                                   (nelec_a,nelec_b))
dm1, dm2, dm3 = fci.rdm.reorder_dm123(dm1, dm2, dm3)

#
# NOTE computing 4-pdm is very slow
#
#
# Spin-traced 4-particle density matrix
# Note make_dm1234 computes  dm4[p,q,r,s,t,u,v,w] = <p^+ q r^+ s t^+ u v^+ w>  which is
# NOT the 4-particle density matrices.  Funciton reorder_dm1234 transforms it
# to the true 4-particle DM  dm4[p,q,r,s,t,u,v,w] = <p^+ r^+ t^+ v^+ w u s q>
# (as well as the 2-particle and 3-particle DMs)
#
dm1, dm2, dm3, dm4 = fci.rdm.make_dm1234('FCI4pdm_kern_sf', fcivec0, fcivec0, norb,
                                         (nelec_a,nelec_b))
dm1, dm2, dm3, dm4 = fci.rdm.reorder_dm1234(dm1, dm2, dm3, dm4)

#
# Spin-traced 4-particle transition density matrix
#
dm1, dm2, dm3, dm4 = fci.rdm.make_dm1234('FCI4pdm_kern_sf', fcivec0, fcivec1, norb,
                                         (nelec_a,nelec_b))
dm1, dm2, dm3, dm4 = fci.rdm.reorder_dm1234(dm1, dm2, dm3, dm4)