BoothGroup / Vayesta

A Python package for wave function-based quantum embedding
Apache License 2.0
33 stars 8 forks source link

Fix afm single site #80

Closed cjcscott closed 1 year ago

cjcscott commented 1 year ago

For a system which has strong AFM character running the default single-site EWF calculation in Vayesta with a spin-symmetry broken reference results in an error due to CISD coefficient generation on master. Here's is a minimal example using a dissociated hydrogen molecule, but I originally came across it in H-ring dissociations and I'd expect the same for the high-U Hubbard model.

from pyscf import gto, scf
from vayesta.misc.molecules import ring

mol = gto.M()
mol.atom = "H  0  0  0;  H  0  0  10.0"
mol.basis="sto-3g"
mol.build()

umf = scf.UHF(mol).density_fit()
umf.kernel()
umf.kernel(dm0=umf.make_rdm1(mo_coeff=umf.stability()[0]))

uewf = vayesta.ewf.EWF(umf, solver="FCI", bath_options={"bathtype":"dmet"})
uewf.kernel()

from which I get the error

 Traceback (most recent call last):
  File "/home/charlie/calc/bugfix/cisd_amps/minex.py", line 18, in <module>
    uewf.kernel()
  File "/home/charlie/miniconda3/envs/vayesta/lib/python3.9/site-packages/vayesta/ewf/ewf.py", line 179, in kernel
    x.kernel()
  File "/home/charlie/miniconda3/envs/vayesta/lib/python3.9/site-packages/vayesta/ewf/fragment.py", line 293, in kernel
    cluster_solver.kernel(eris=eris, **init_guess)
  File "/home/charlie/miniconda3/envs/vayesta/lib/python3.9/site-packages/vayesta/solver/fci.py", line 164, in kernel
    self.c0, self.c1, self.c2 = self.get_cisd_amps(self.civec)
  File "/home/charlie/miniconda3/envs/vayesta/lib/python3.9/site-packages/vayesta/solver/fci.py", line 290, in get_cisd_amps
    c2ab = einsum('i,j,ij->ij', t1signa, t1signb, civec[t1addra[:,None],t1addrb])
TypeError: list indices must be integers or slices, not tuple

This arises because in these systems each cluster consists of a single occupied orbital of one spin, and a single virtual orbital of the other spin. This means that there are no valid excitations within the cluster, and so runs into an issue where the call to pyscf to obtain indexing and sign changes t1addra, t1signa = pyscf.ci.cisd.tn_addrs_signs(norba, nocca, 1) sets both variables to a zero-length list when there are no possible excitations, rather than a zero-length array which might be expected given the usual return value is an array. This then leads to a later error due to python lists not supporting zero-length indexing, unlike numpy arrays.

This has already come up in the to_cisdtq branch (where obtaining an empty list for the t3 and t4 indexing is a more common occurrence), so I've just borrowed the fix already implemented there to ensure these are always array-valued. Given to_cisdtq affects many of the same classes and is expected to be merged soon it seems more straightforward to add it here rather than to master.

This branch also adds a test for this case to avoid future regression. Doing this required a way to perform tests on a UHF solution which is stable under spin stability analysis, so I've added the TestMolecule.uhf_stable() functionality in vayesta/tests/testsystems.py for this purpose without breaking lots of existing UHF tests.

maxnus commented 1 year ago

Thanks for this Charlie. Is the return type inconsistent already for the PySCF function? Then it should be fixed on their side already. Although I agree we should keep the conversion on our side too, no matter what.

Are the tests not running because the PR is not for master? If this fix is independent of the to_cisdtq branch, maybe we can apply it to master instead.

@obackhouse Should we change to run tests for PRs to all branches?

cjcscott commented 1 year ago

The return type is inconsistent in the PySCF function- specifically this function. Having this be consistent would obviously be preferable, but ensuring what we get is consistent seems sensible for now. While this is independent of to_cisdtq in direct functionality, given that branch refactors a decent amount of connected code basing this fix on it seemed more practical than essentially implementing twice (once with and once without the refactor). Assuming it's merging to master soon it shouldn't lead to a major delay, and we've not shied away from having branches with multiple functional changes in them before.

Any unexpected changes should be caught before getting into master, either in to_cisdtq if this PR is accepted first or in this PR if to_cisdtq is merged first, so there's not to much risk. But I agree that we should have tests on all PR's- I'll do a PR changing ci.yaml as such. Wonder if that'll start tests running on this one...