Brueckner self-consistency broken for open-shell systems

[obackhouse]

Brueckner orbital self-consistency has not been written in a way that supports open shell systems where the number of occupied alpha and beta electrons is inequal, as per #41 (comment).

[vvp-nsk]

Hej!

It seems to work now. However, I am puzzled a bit why DF-CC does not support single precision. Specifically, an input given below doesn't work:

import numpy as np
from pyscf import gto, scf

from ebcc import UEBCC
from ebcc.precision import single_precision


# Define the molecule using PySCF
mol = gto.Mole()
mol.atom = "H 0 0 0; F 0 0 1.1"
mol.basis = "cc-pvdz"
mol.charge = 1
mol.spin = 1
mol.build()

# Run a UHF calculation using PySCF
mf = scf.UHF(mol)
mf = mf.density_fit()
mf.kernel()

# Run a CCSD calculation
with single_precision():
     ccsd = UEBCC(mf, ansatz="CC2")
     ccsd.kernel()

It fails with the following error message:

  File "ebcc/ebcc/uebcc.py", line 152, in init_amps
    tn[comb] = eris[comb_t][key_t].swapaxes(1, 2) / self.energy_sum(key, comb)
  File "ebcc/ebcc/cderis.py", line 186, in __getattr__
    e1 = getattr(v1, key[:2])
  File "ebcc/ebcc/cderis.py", line 94, in __getattr__
    block = ao2mo._ao2mo.nr_e2(
  File "lib/python3.10/site-packages/pyscf/ao2mo/_ao2mo.py", line 160, in nr_e2
    assert (mo_coeff.dtype == numpy.double)

In my eye, it is natural to combine DF and single precision to speedup CC calculations further.

With best regards,
Victor

[obackhouse]

Hej!

It seems to work now. However, I am puzzled a bit why DF-CC does not support single precision. Specifically, an input given below doesn't work:

import numpy as np
from pyscf import gto, scf

from ebcc import UEBCC
from ebcc.precision import single_precision


# Define the molecule using PySCF
mol = gto.Mole()
mol.atom = "H 0 0 0; F 0 0 1.1"
mol.basis = "cc-pvdz"
mol.charge = 1
mol.spin = 1
mol.build()

# Run a UHF calculation using PySCF
mf = scf.UHF(mol)
mf = mf.density_fit()
mf.kernel()

# Run a CCSD calculation
with single_precision():
     ccsd = UEBCC(mf, ansatz="CC2")
     ccsd.kernel()

It fails with the following error message:

  File "ebcc/ebcc/uebcc.py", line 152, in init_amps
    tn[comb] = eris[comb_t][key_t].swapaxes(1, 2) / self.energy_sum(key, comb)
  File "ebcc/ebcc/cderis.py", line 186, in __getattr__
    e1 = getattr(v1, key[:2])
  File "ebcc/ebcc/cderis.py", line 94, in __getattr__
    block = ao2mo._ao2mo.nr_e2(
  File "lib/python3.10/site-packages/pyscf/ao2mo/_ao2mo.py", line 160, in nr_e2
    assert (mo_coeff.dtype == numpy.double)

In my eye, it is natural to combine DF and single precision to speedup CC calculations further.

With best regards, Victor

yep that's another bug

[obackhouse]

Should be fixed by #47. Thanks for identifying so many edge cases @vvp-nsk, it's much appreciated.

A couple notes:

• Density fitting has another bug if your PySCF calculation uses out-of-core integrals (i.e. for large systems). If you encounter this then the workaround is to set the max_memory higher so that the integrals are incore, as per Density fitting fails when not incore #48.
• I'm not sure of the validity of Brueckner self-consistency for single precision, the matrix logarithm seems to be poorly conditioned -- I've cast the output to a real valued matrix for now, again I'm not an expert on this so I'm not sure the best approach here.

Again, let me know if you catch anything else 😄

[vvp-nsk]

Hej!

Thanks for the prompt bug fixing, indeed. I do have a challenging example when Bruckner CC implementation doesn't work in EBCC while PySCF works correctly. Let me narrow the testbed to a tiny problem size first.

I also have a general question. For CC methods of practical interest, such as CC2 and CCSD, both right and left CC eigenvectors are implemented in ECBB. Since 1pdm and 2pdm are available, would it be possible then to implement an expectation value of <S^2> as it can be done in PySCF via call to the pyscf.fci.spin_op.spin_square_general function?

Thank you very much!

With best regards,
Victor

[obackhouse]

I also have a general question. For CC methods of practical interest, such as CC2 and CCSD, both right and left CC eigenvectors are implemented in ECBB. Since 1pdm and 2pdm are available, would it be possible then to implement an expectation value of <S^2> as it can be done in PySCF via call to the pyscf.fci.spin_op.spin_square_general function?

There is some support of the EOM eigenvectors, yes -- but they're not needed to calculate the 1DM and 2DM, which are also separately available for both CC2 and CCSD. Be sure to call solve_lambda first to solve the lambda equations. I think it should be as easy as this to use that PySCF function:

import numpy as np
from pyscf import gto, scf
from ebcc import UEBCC

mol = gto.M(
    atom="Be 0 0 0; H 0 0 2",
    basis="cc-pvdz",
    spin=1,
    verbose=0,
)

mf = scf.UHF(mol)
mf.kernel()

uebcc = UEBCC(mf, ansatz="CC2")
uebcc.kernel()
uebcc.solve_lambda()

dm1 = uebcc.make_rdm1_f()
dm2 = uebcc.make_rdm2_f()

from pyscf.fci.spin_op import spin_square_general
ovlp = np.eye(mol.nao)
mo_coeff = (np.eye(mol.nao),) * 2
print(spin_square_general(dm1.aa, dm1.bb, dm2.aaaa, dm2.aabb, dm2.bbbb, mo_coeff, ovlp=ovlp))

where the overlap and mo_coeff are set to identity, since these DMs are calculated in the (orthogonal) MO basis.

I should note that I have an open issue about a possible hiccup in my CC2 implementation (#9), which I haven't been able to resolve, so you should verify any results that use the CC2 lambda equations or DMs. This (possible) issue doesn't extend to CCSD, though.

[vvp-nsk]

Fantastic!

Could you please do me a favor and rework the given above code on computing <S^2> for the frozen-core case as well:

import numpy as np
from pyscf import gto, scf

from ebcc import UEBCC, Space

# Define the molecule using PySCF
mol = gto.Mole()
mol.atom = "H 0 0 0; F 0 0 1.1"
mol.basis = "cc-pvdz"
mol.charge = 1
mol.spin = 1
mol.build()

# Run a UHF calculation using PySCF
mf = scf.UHF(mol)
#mf = mf.density_fit()
mf.kernel()

frozen_a = np.zeros_like(mf.mo_occ[0])
frozen_a[:1] = 1
frozen_a = frozen_a.astype(bool)
frozen_b = np.zeros_like(mf.mo_occ[1])
frozen_b[:1] = 1 ;
frozen_b = frozen_b.astype(bool)

space_a = Space(
        mf.mo_occ[0] > 0,
        frozen_a,
        np.zeros_like(mf.mo_occ[0]),
)
space_b = Space(
        mf.mo_occ[1] > 0,
        frozen_b,
        np.zeros_like(mf.mo_occ[1]),
)

# Run a CCSD calculation
ccsd = UEBCC(mf, ansatz="CC2", space=(space_a,space_b))
ccsd.kernel()
...

?

Thank you in advance!

With best regards,
Victor

[obackhouse]

This should do the trick:

import numpy as np
from pyscf import gto, scf
from ebcc import UEBCC, Space
from pyscf.fci.spin_op import spin_square_general

mol = gto.M(
    atom="Be 0 0 0; H 0 0 2",
    basis="cc-pvdz",
    spin=1,
    verbose=0,
)

mf = scf.UHF(mol)
mf.kernel()

frozen_a = np.zeros_like(mf.mo_occ[0])
frozen_a[:1] = 1
frozen_a = frozen_a.astype(bool)
frozen_b = np.zeros_like(mf.mo_occ[1])
frozen_b[:1] = 1
frozen_b = frozen_b.astype(bool)

space_a = Space(
        mf.mo_occ[0] > 0,
        frozen_a,
        np.zeros_like(mf.mo_occ[0]),
)
space_b = Space(
        mf.mo_occ[1] > 0,
        frozen_b,
        np.zeros_like(mf.mo_occ[1]),
)
space = (space_a, space_b)

uebcc = UEBCC(mf, ansatz="CC2", space=space)
uebcc.kernel()
uebcc.solve_lambda()

dm1_act = uebcc.make_rdm1_f()
dm2_act = uebcc.make_rdm2_f()

ovlp = np.eye(mol.nao)
mo_coeff = (np.eye(mol.nao),) * 2

dm1_a, dm1_b = mf.make_rdm1(mo_coeff, mf.mo_occ)
dm2_aa, dm2_ab, dm2_bb = mf.make_rdm2(mo_coeff, mf.mo_occ)

dm1_a[np.ix_(space_a.correlated, space_a.correlated)] = dm1_act.aa
dm1_b[np.ix_(space_b.correlated, space_b.correlated)] = dm1_act.bb

dm2_aa[np.ix_(space_a.correlated, space_a.correlated, space_a.correlated, space_a.correlated)] = dm2_act.aaaa
dm2_ab[np.ix_(space_a.correlated, space_a.correlated, space_b.correlated, space_b.correlated)] = dm2_act.aabb
dm2_bb[np.ix_(space_b.correlated, space_b.correlated, space_b.correlated, space_b.correlated)] = dm2_act.bbbb

print(spin_square_general(dm1_a, dm1_b, dm2_aa, dm2_ab, dm2_bb, mo_coeff, ovlp=ovlp))

It finds the bare (uncorrelated) 1DM and 2DM using the mean-field function in PySCF, and then embeds the correlated (non-frozen) part using the boolean masks in the Space classes. Not the prettiest -- I'll aim to add a with_frozen flag to the make_rdm1_f and make_rdm2_f functions, because this is definitely of general interest.

[vvp-nsk]

Hej!

Thank you very much for your assistance and support!

Attached please find a problem reproducer. The UBCCD malfunctions at the 2nd iteration:

Converged.
E(corr) = -0.6577317650
E(tot)  = -2519.1924758110
Brueckner options:
 > e_tol:  1e-06
 > t_tol:  1e-05
 > max_iter:  20
 > diis_space:  12
Solving for Brueckner orbitals.
Iter   Energy (corr.)  Converged        Δ(Energy)             |T1|
   0    -0.6577317650       True
   1     0.0000000000       True          0.65773                0

The problem is not related to certain CC ansatz and remains with CC2 as well.
bccd_incorrect.zip

With best regards,
Victor

[obackhouse]

I'm struggling to get your chkfile to work on my end -- can you instead send the script including the full Mole setup and UHF calculations?

[vvp-nsk]

Hej!

Sorry for the late response. Indeed, there was one file missing. I just re-uploaded testbed and now it should work out.
bccd_incorrect.zip

With best regards,
Victor

[obackhouse]

I think this should be fixed on #49. Sorry this one took a bit longer!

CC2 seems to converge well for this system but CCSD is a little trickier. I think doing BCCD in single precision will be ambitious, numerically speaking, if that was your plan.

[vvp-nsk]

Hej!

Everything works like a charm now. Thank you!

My guess that often SP would be sufficient to get Brueckner orbitals converged. Anyway, among open-source CC packages, UBCC2 is available only in EBCC 👍

With best regards,
Victor

[obackhouse]

Great -- thanks for finding all these bugs.