BoothGroup / ebcc

Coupled cluster models for both purely electron and coupled electron-boson models, with a focus on generality and model extensibility.
https://boothgroup.github.io/ebcc/
MIT License
11 stars 2 forks source link

Brueckner self-consistency broken for open-shell systems #44

Closed obackhouse closed 8 months ago

obackhouse commented 8 months ago

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 https://github.com/BoothGroup/ebcc/issues/41#issuecomment-1929989963.

vvp-nsk commented 8 months ago

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 commented 8 months ago

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 commented 8 months ago

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

A couple notes:

Again, let me know if you catch anything else :smile:

vvp-nsk commented 8 months ago

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 commented 8 months ago

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 commented 8 months ago

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 commented 8 months ago

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 commented 8 months ago

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 commented 8 months ago

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 commented 8 months ago

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 commented 8 months ago

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 commented 8 months ago

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 :+1:

With best regards, Victor

obackhouse commented 8 months ago

Great -- thanks for finding all these bugs.