Frozen Core / Natural Orbitals

[PaulWAyers]

@msricher I remember you were looking at frozen core. These notes seem useful
http://vergil.chemistry.gatech.edu/notes/ci.pdf

[PaulWAyers]

This is probably not for the release

[msricher]

@PaulWAyers It really should be.... it would let us compute much, much larger systems and it's always been a big limiter for us when we try to use PyCI for computational studies.

Also, we need the natural orbital transformation. That should really be there. It's an easy way to do pseudo orbital-optimization.

[PaulWAyers]

Natural orbital transformation is pretty easy. Frozen core should be doable too, but we've failed in the past....

Maybe start with n.o. transformation then take another crack at frozen core after that?

[msricher]

@PaulWAyers Sure, I can do this.

I'm trying to transform the restricted MO integral arrays using $\gamma^{\alpha\alpha}$. See this code:

import numpy as np

import pyci

from pyci.test import datafile


# Make CISD wfn
ham = pyci.secondquant_op(datafile("be_ccpvdz.fcidump"))
wfn = pyci.fullci_wfn(ham.nbasis, 2, 2)
pyci.add_excitations(wfn, 0, 1, 2)

# Solve wfn, compute RDMs
op = pyci.sparse_op(ham, wfn)
es, cs = op.solve(n=1, ncv=30, tol=1.0e-6)
d1, d2 = pyci.compute_rdms(wfn, cs[0])

# Use eigenvectors of \alpha\alpha 1-RDM to transform 1- and 2- electron MO integrals
_, u = np.linalg.eig(d1[0])
ecore = ham.ecore
one_mo = np.einsum("pi,qj,ij->pq", u, u, ham.one_mo)
two_mo = np.einsum("pi,qj,rk,sl,ijkl->pqrs", u, u, u, u, ham.two_mo)

# Compute energy of 1-determinant (Hartree-Fock) wfn using transformed MO integrals
no_ham = pyci.secondquant_op(ecore, one_mo, two_mo)
no_wfn = pyci.fullci_wfn(no_ham.nbasis, 2, 2)
pyci.add_excitations(no_wfn, 0)
no_op = pyci.sparse_op(no_ham, no_wfn)
no_es, no_cs = no_op.solve(n=1, ncv=1, tol=1.0e-6)

# Compare
print(f"  CISD E:      {es[0]:16.9e}")
print(f"  Nat. Orb. E: {no_es[0]:16.9e}")

Which gives an incorrect result:

[master]-[michelle@t14s-g4]-[~/Git/pyci]-% python no.py
  CISD E:      -1.461735579e+01
  Nat. Orb. E: -1.457209016e+01

This might be the wrong way to go about this?

[msricher]

@PaulWAyers
I can also try using $\gamma^{\alpha\alpha} + \gamma^{\beta\beta}$, d1[0] + d1[1]:

  CISD E:      -1.461735579e+01
  Nat. Orb. E: -1.457209030e+01

Still wrong.

[msricher]

Regarding the frozen core approximation, from Sherill's notes... I see that it uses the C matrix from Hartree-Fock. I think you have to manipulate the AOs. This is actually outside of the scope I envisioned for PyCI... maybe it doesn't belong here. It should, then, be a post-processing routine in the mean-field code. Do you agree @PaulWAyers

[PaulWAyers]

It could be done with MOs, but I can see that it belongs in meanfield or maybe even better, gbasis, as it is basically computing the integrals for PyCI.

[msricher]

We've decided this is outside the scope of this repo.