After talking to @gabrielasd I think it might be a good idea to remove a dependency on the specific optimizer (e.g cma). We can still have predefined optimization in fanci, but I propose to add support of arbitrary optimizer of choice. Potentially, it should make code more modular and flexible.
This problem arose from a collaborative work with Ramon's group, where they have a serious problem of using a layers of wrappers around the FanCI and FanPy just to perform a regular minimization. I think such splitting should help them
Also, we can remove dependency on cma package which doesn't seem very necessary to me either.
@PaulWAyers @msricher do you think we should do it?

My tests for the wave function data are fine, but for the energies, CI coefficients, and RDMs, they are lacking. We need to pull a bunch of results from PySCF (or other sources) and compare them against PyCI's results. See the test files pyci/test/test_*.py.

Semidefinite Constraints

In the Semidefinite Constraints the first step is to map the the density matrix, $\Gamma$ , to a positive semidefinite constraint. Do we have the definition of $\mathcal{L}(\Gamma)$ ?

Equality Constraints

In the Equality Constraints part if the equality is not satisfied do we need to check if $Tr[\mathbf{B}_m \Gamma] > 0$. If it's not, should we then replace $\mathbf{B}_m$ with $-\mathbf{B}_m$? After this adjustment, is the projection based on the constraint $Tr[\mathbf{B}_m \Gamma] \leq 0$?"

@PaulWAyers

I'm trying to run seniority 0+2+4 CI for H6. There should be determinants of seniority-4, but it raises:
ValueError: invalid seniority number in seniorities = (0, 2, 4)

Minimal code using H6.FCIDUMP (attached as a .txt)

H6_2.00.txt

import pyci

filename = "H6_2.00.FCIDUMP"
occs = (3,3)
ham = pyci.hamiltonian(filename)

seniorities = [0,2,4]
wfn = pyci.fullci_wfn(ham.nbasis, *occs)
pyci.add_seniorities(wfn, *seniorities)

Some useful formulas

After getting from @RichRick1 a FCIDUMP file of a Hubbard model I am getting the values of the cost function, the residuals, the optimality and active_mask equal to zero when optimizing the AP1roG wavefunction. The solution found is:

x: [ 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 3.000e+01]

Optimizing the pCCD+S wave function give me residuals near zero ($10^{-15}$ order) and a solution x:

x: [-1.334e+01 1.334e+01 ... -5.454e-14 3.000e+01]

demo pyci.ipynb.txt
wavefunctions.ipynb.txt

This function can go into the rdm.constraints module which can be found at rdm/constraints.py file

The task is to select Slater Determinants based on a given cost function.
There are methods for selecting SDs based on excitation level and seniority are already implemented.
Also, user is able to select SDs based on number of orbital nodes. This is also implemented with GKCI algorithm
We need to implement selection of SD given energy of the orbitals (similar to the number of orbital nodes)
ToDo:

• #33
• #32
• #37
  1. Create almost the same version of odometer_one_spin and odometer_two_spin functions in the cost.py file, but changing the selection criteria by removing smoothness factor.
  2. (Recommended) Implement general version of odometer_one_spin and odometer_two_spin functions in utility.py file that takes selection criteria (probably in the form of another function) as input.

UPD:

• #36

We are implementing FanPT. The basic idea is in the recent paper.

Recall that the FanCI equations are

$$ \braket{\mathbf{m}|\hat{H}|\Psi(\mathbf{c})} = E \braket{\mathbf{m}|\Psi(\mathbf{c})} \qquad \qquad \forall \mathbf{m} \in \mathcal{P} $$

where $\mathcal{P}$ is the projection space (user-specified), $\mathbf{c}$ are the wavefunction parameters, and the nonlinear FanCI equations are usually overdetermined and solved in a least-squares sense. The set of determinants $\mathbf{n} \in \mathcal{S}$ is defined as the space of determinants that are linked to determinants in $\mathcal{P}$ by the Hamiltonian, typically just single- and double-replacements of the occupied orbitals.

We write the Hamiltonian as a zeroth-order part (for which the FANCI equations have been solved) and a first-order part,

$$ \hat{H}(\lambda) = \hat{H}(\lambda_0) + (\lambda - \lambda_0) \hat{H}'(\lambda_0) + \tfrac{1}{2} (\lambda - \lambda_0)^2 \hat{H}''(\lambda_0) +\cdots $$

Going through some math, we get the equation

$$ \sum_{\mathbf{n} \in \mathcal{S}} f_\mathbf{n}(\lambda_0,\mathbf{c}) E'(\lambda_0) + \sum_{\mathbf{n} \in \mathcal{S}} \left(\braket{\mathbf{m}|\hat{H}(\lambda_0)|\mathbf{n}} - E(\lambda_0) \delta_{\mathbf{m} \mathbf{n}}\right) \frac{d f_{\mathbf{n}}(\mathbf{c})}{d \mathbf{c}(\lambda_0)} \cdot \mathbf{c}'(\lambda_0) = \sum_{\mathbf{n} \in \mathcal{S}} \braket{\mathbf{m}|\hat{H}'(\lambda_0) |\mathbf{n}} f_\mathbf{n}(\lambda_0,\mathbf{c}) $$

The key ingredients are the Hamiltonian matrix elements. The right hand side is easy to evaluate in terms of transition-2DMs, which appear on the right-hand-side,

$$
\text{RHS} = \text{Tr}[\hat{H}'(\lambda_{0}){pqrs} \Gamma{pqrs}(\lambda_{0})]
$$

My tendency would be to implement this as some functions that build the linear equations to solve, and then solve them. I realize that the (old) implementation is somewhat more general, but we may also want to focus on a specific one or two flavors if it makes the code a lot simpler/faster. Ramon, what do you think the best flavors are? (It will depend on how much computational complexity and code complexity we add by writing a general treatment.)

There are inequality constraints and (wrapper of inequality constraints) equality constraints.

Implement the projection onto an inequality constraints.
It is helpful to rewrite the inequality constraint so that it is traceless, since the trace of the density matrix is fixed. Consider the cosntraint

$$ Tr\left[\mathbf{A}_l \Gamma\right] \geq 0 $$

and define

$$ \tilde{\mathbf{A}}_l=\mathbf{A}_l-Tr\left[\mathbf{A}_l\right] \mathbf{I} $$

The inequality can be rewritten as

If the inequality is satisfied, then there is no constraint. If the inequality is not satisfied, then the projection is

Given $\Gamma$ and number of electrons, function should return projection $P(\Gamma)$

UPD: all of these can go into the rdm.constraints module which can be found at rdm/constraints.py file

It would be useful to be able to perform orbital optimization. With the recent work of @ShuoyangW , we have the gradient (and soon, the Hessian) with respect to orbitals. Being able to wrap PyCI in an orbital optimizer is a good thing.

This basically means that one needs to be able to consider PyCI a function where 1- and 2-electron integrals are read in, and the energy and 1- and 2-electron density matrices come out. This then needs to be wrapped in another function that takes 1- and 2-electron integrals and performs a unitary transformation of them based on the step accepted from an optimizer, as the optimizer expects to receive the energy, gradient, and (perhaps) Hessian as a function of an argument (probably the antiHermitian matrix that generates the unitary transformation).

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

There are 3 algorithms:

• #51
• #52
• #53

UPD: all of these can go into the rdm.algorithms module which can be found at rdm/algorithms.py file

Some useful formulas

There are 6 semidefinite constraints, 3 easy and 3 hard. One could use:

• #47
• #48
• #49
• #50

For the constraints, the function should be, for example:

def calc_P(gamma):
      ...
      return P

UPD: all of these can go into the rdm.conditions module which can be found at rdm/conditions.py file

To proceed with implementing FanPT capabilities in PyCI, we need to calculate the derivatives of the overlap matrices.
This feature is already implemented for AP1roG, but is not currently available in the code for pCCD+S.

@alnaba1, to test the odometer algorithm we can use H2 and H4 systems in the STO-3G basis set. By providing different cost function we can see which determinants are selected.
We can use pyscf package to generate .fcidump file and test the results on it:

from pyscf import gto
mol = gto.Mole()
mol.build( atom = "H 0 0 0; H  0 1 0", basis = 'sto-3g')

# Import modules
import pyci

# System information
filename = datafile("lih_sto6g.fcidump")
ham = pyci.hamiltonian(filename)

@PaulWAyers do you think it's good testing systems?

Todo

• Clean up the contributed code from the course. I've just been doing black -l 100 on my own Python files so far, but I will add some proper pre-commit hooks to take care of this once the course is done.
• Improve C++ code comments.
• Bump version number to 0.7.0.
• Include in the README and website a section like the following with a bibtex entry once the paper is published:
  Citing PyCI:
  -----------------
  Please cite PyCI in accordance with the following if you use it for published research:

  @article{2024pyci
  ...
  }

I would also appreciate others' thoughts on what should be done during this time.

For the odometer algorithm we need to

• Move odometer_one_spin, odometer_two_spin into utils.py file
• Add input (user) parameter q_max
• Change name of parameter nodes into parameter cost

NotaBene:
The current selection rule isn't changing:
$$f(\tau, \mathbf{c})=\sum_{p \in \text { occupied orbitals }} c_p+\tau \max _{p \in \text { occupied orbitals }} c_p$$

I made the domain pyci.qcdevs.org

It would be nice to have 3- and 4-electron reduced density matrices, especially (perhaps even exclusively) for the seniority-zero case, where the number of indices in these matrices is very small. This might go together with a special implementation for the DMs of seniority-zero wavefunctions.

I'm tentatively tagging @DanielCalero1 on this, as it may be useful for that project. It isn't a super high priority (yet) though.

1. Import odometers from utility.py
2. Leave q_max calculation as it is, but move it outside of odometers and then pass it as a parameter