Forte Tutorial 3: Forte's determinant class

In this tutorial we are going to explore how to interact with forte in Jupyter notebooks using the Python API.

Import modules

The first step necessary to interact with forte is to import psi4 and forte



In [ ]:

    
import psi4
import forte
import forte.utils

First we will run psi4 using the function forte.utils.psi4_scf



In [ ]:

    
# setup xyz geometry
geom = """
O
H 1 1.0
H 1 1.0 2 180.0
"""
(E_scf, wfn) = forte.utils.psi4_scf(geom,basis='sto-3g',reference='rhf')
print(f'SCF Energy = {E_scf}')

Starting Forte and reading options



In [ ]:

    
forte.startup()

from forte import forte_options

options = psi4.core.get_options() # options = psi4 option object
options.set_current_module('FORTE') # read options labeled 'FORTE'
forte_options.get_options_from_psi4(options)

Setting the molecular orbital spaces



In [ ]:

    
# Setup forte and prepare the active space integral class
mos_spaces = {'FROZEN_DOCC' :     [1,0,0,0,0,0,0,0], # freeze the O 1s orbital
              'RESTRICTED_DOCC' : [1,0,0,0,0,1,0,0]}
mo_space_info = forte.make_mo_space_info_from_map(wfn,mos_spaces,[])



In [ ]:

    
mo_space_info.size('ACTIVE')

Building a `ForteIntegral` object to read integrals from psi4

In Forte there are two classes responsible for handling integrals:

ForteIntegral: reads the integrals from psi4 and stores them in varios formats (conventional, density fitting, Cholesky, ...).
ActiveSpaceIntegrals: stores a copy of all integrals and it is used by active space methods. This class only stores a subset of the integrals and includes an effective potential due to non-active doubly occupied orbitals.

We will first build the ForteIntegral object via the function make_forte_integrals



In [ ]:

    
ints = forte.make_forte_integrals(wfn, options, mo_space_info)
print(f'Number of molecular orbitals: {ints.nmo()}')
print(f'Number of correlated molecular orbitals: {ints.ncmo()}')

# the space that defines the active orbitals. We select only the 'ACTIVE' part
active_space = 'ACTIVE'
# the space(s) with non-active doubly occupied orbitals
core_spaces = ['RESTRICTED_DOCC']

as_ints = forte.make_active_space_ints(mo_space_info, ints, active_space, core_spaces)

print(f'Frozen-core energy = {as_ints.frozen_core_energy()}')
print(f'Nuclear repulsion energy = {as_ints.nuclear_repulsion_energy()}')
print(f'Scalar energy = {as_ints.scalar_energy()}')

Creating determinants



In [ ]:

    
d = forte.Determinant()



In [ ]:

    
nact = mo_space_info.size('ACTIVE')
d.str(nact)

d.create_alfa_bit(1)



In [ ]:

    
d.str(nact)
d.create_alfa_bit(2)



In [ ]:

    
d.str(nact)



In [ ]:

    
d.destroy_alfa_bit(2)



In [ ]:

    
d.str(nact)

Creating the HF determinant



In [ ]:

    
nirrep = mo_space_info.nirrep()
nactpi = mo_space_info.dimension('ACTIVE').to_tuple()

# compute the number of alpha electrons per irrep
nact_aelpi = wfn.nalphapi() - mo_space_info.dimension('FROZEN_DOCC') - mo_space_info.dimension('RESTRICTED_DOCC')
nact_aelpi = nact_aelpi.to_tuple()      
# compute the number of beta electrons per irrep
nact_belpi = wfn.nbetapi() - mo_space_info.dimension('FROZEN_DOCC') - mo_space_info.dimension('RESTRICTED_DOCC')
nact_belpi = nact_belpi.to_tuple()           

print(f'Number of alpha electrons per irrep: {nact_aelpi}')
print(f'Number of beta electrons per irrep:  {nact_belpi}')
print(f'Number of active orbtials per irrep: {nactpi}')

ref = forte.Determinant()

# we loop over each irrep and fill the occupied orbitals 
irrep_start = [sum(nactpi[:h]) for h in range(nirrep)]
for h in range(nirrep):
    for i in range(nact_aelpi[h]): ref.set_alfa_bit(irrep_start[h] + i, True)
    for i in range(nact_belpi[h]): ref.set_beta_bit(irrep_start[h] + i, True)        
    
print(f'Reference determinant: {ref.str(nact)}')

We can now compute the energy of the determinant as $\langle \Phi | \hat{H} | \Phi \rangle$ using the slater_rules function in the ActiveSpaceIntegrals class



In [ ]:

    
as_ints.slater_rules(ref,ref) + as_ints.scalar_energy() + as_ints.nuclear_repulsion_energy()

We can also print all the integrals at once (see the output.dat file)



In [ ]:

    
import itertools
import functools

dets = []
orbs = range(nact)

# get the symmetry of each active orbital
act_sym = mo_space_info.symmetry('ACTIVE')

nact_ael = sum(nact_aelpi)
nact_bel = sum(nact_belpi)
print(f'Number of alpha electrons: {nact_ael}')
print(f'Number of beta electrons:  {nact_bel}')

# specify the target symmetry
sym = 0

# generate all the alpha strings
for astr in itertools.combinations(orbs, nact_ael):
    # compute the symmetry of the alpha string
    asym = functools.reduce(lambda i, j:  act_sym[i] ^ act_sym[j], astr)
    # generate all the beta strings
    for bstr in itertools.combinations(orbs, nact_bel):
        # compute the symmetry of the beta string
        bsym = functools.reduce(lambda i, j:  act_sym[i] ^ act_sym[j], bstr)    
        # if the determinant has the correct symmetry save it
        if (asym ^ bsym) == sym:
            d = forte.Determinant()
            for i in astr: d.set_alfa_bit(i, True)
            for i in bstr: d.set_beta_bit(i, True)                
            dets.append(d)

print(f'==> List of FCI determinants <==')
for d in dets:
    print(f'{d.str(4)}')



In [ ]:

    
import numpy as np

ndets = len(dets)
H = np.ndarray((ndets,ndets))
for I, detI in enumerate(dets):
    for J, detJ in enumerate(dets):
        H[I][J] = as_ints.slater_rules(detI,detJ)

# or we could use the more fancy looping below that avoid computing half of the matrix elements
# for I, detI in enumerate(dets):
#     H[I][I] = as_ints.slater_rules(detI,detI) # diagonal term
#     for J, detJ in enumerate(dets[:I]):
#         HIJ = as_ints.slater_rules(detI,detJ) # off-diagonal term (only upper half)
#         H[I][J] = H[J][I] = HIJ
   
print(H)
evals, evecs = np.linalg.eigh(H)

psi4_fci = -74.846380133240530
print(f'FCI Energy = {evals[0] + as_ints.scalar_energy() + as_ints.nuclear_repulsion_energy()}')
print(f'FCI Energy Error = {evals[0] + as_ints.scalar_energy() + as_ints.nuclear_repulsion_energy()- psi4_fci}')

index_hf = dets.index(ref)
print(f'Index of the HF determinant in the FCI vector {index_hf}')



In [ ]:

    
forte.cleanup()