CCSDt with UHF results in ValueError

[abhishekkhedkar09]

Here's the error for MWE for N2 with UHF and atomic fragmentation, full bath and DMET active space. (inp attached below)

          |      ebcc:
          |       > Version:  1.4.0
          |       > Git hash: ff52401
          |      OMP_NUM_THREADS = 24
          |      
          |      UCCSDt
          |      ******
          |      Options:
          |       > e_tol:  1e-08
          |       > t_tol:  1e-08
          |       > max_iter:  200
          |       > diis_space:  12
          |      Ansatz: CCSDt
          |      Space: ((7o, 53v) [(5o, 3v) active], (7o, 53v) [(5o, 3v) active])
          |      Iter   Energy (corr.)        Δ(Energy)    Δ(Amplitudes)
          |         0    -0.8248985092
          |      Time for CCSDt solver: 3.3 s
Traceback (most recent call last):
  File "/users/k2262473/Calculations/testing/CCSDt/Solver_CCSDt/input.py", line 26, in <module>
    emb.kernel()
  File "/users/k2262473/dev/Vayesta/vayesta/ewf/ewf.py", line 169, in kernel
    x.kernel()
  File "/users/k2262473/dev/Vayesta/vayesta/ewf/fragment.py", line 236, in kernel
    cluster_solver.kernel()
  File "/users/k2262473/dev/Vayesta/vayesta/solver/ebcc.py", line 33, in kernel
    mycc.kernel()
  File "/users/k2262473/dev/ebcc/ebcc/rebcc.py", line 380, in kernel
    amplitudes = self.update_amps(amplitudes=amplitudes, eris=eris)
  File "/users/k2262473/dev/ebcc/ebcc/uebcc.py", line 238, in update_amps
    res = func(**kwargs)
  File "/users/k2262473/dev/ebcc/ebcc/codegen/UCCSDt.py", line 120, in update_amps
    t3new_aaaaaa += einsum(v.aabb.OVov, (0, 1, 2, 3), t3.abaaba[np.ix_(soa,sob,sOfa,sva,svb,sVfa)], (4, 2, 5, 6, 3, 7), (4, 5, 0, 6, 7, 1)) * 2.0
ValueError: operands could not be broadcast together with shapes (7,7,5,53,53,3) (7,5,5,53,3,3) (7,7,5,53,53,3) 

input.py.txt output.txt

The same with RHF runs without an error and produces correct results.

[abhishekkhedkar09]

If I use density fitting I get the following error in EBCC

Traceback (most recent call last):
  File "/users/k2262473/Calculations/testing/ebcc_CCSDt/input.py", line 15, in <module>
    eb = EBCC(mf, ansatz='CCSDt')
  File "/users/k2262473/dev/ebcc/ebcc/__init__.py", line 145, in EBCC
    return UEBCC(mf, *args, **kwargs)
  File "/users/k2262473/dev/ebcc/ebcc/rebcc.py", line 279, in __init__
    self._eqns = self.ansatz._get_eqns(self.spin_type)
  File "/users/k2262473/dev/ebcc/ebcc/ansatz.py", line 139, in _get_eqns
    eqns = importlib.import_module("ebcc.codegen.%s" % name)
  File "/scratch/grp/morty/k2262473/miniconda/envs/pygnme/lib/python3.10/importlib/__init__.py", line 126, in import_module
    return _bootstrap._gcd_import(name[level:], package, level)
  File "<frozen importlib._bootstrap>", line 1050, in _gcd_import
  File "<frozen importlib._bootstrap>", line 1027, in _find_and_load
  File "<frozen importlib._bootstrap>", line 1004, in _find_and_load_unlocked
ModuleNotFoundError: No module named 'ebcc.codegen.UDFCCSDt'

here's a better MWE without Vayesta:

import pyscf
import pyscf.gto
import pyscf.scf
import ebcc
from ebcc import EBCC 
mol=pyscf.gto.Mole()
mol.atom = ["N 0 0 0", "N 0 0 2"]
mol.basis = "CCpvtz"
mol.output = "pyscf.out"
mol.build()
mf = pyscf.scf.UHF(mol)
#mf = mf.density_fit()
mf.kernel()
eb = EBCC(mf, ansatz='CCSDt')
eb.kernel()

the error:

numpy:
 > Version:  1.25.2
 > Git hash: N/A
pyscf:
 > Version:  2.4.0
 > Git hash: b27e4f9fb
ebcc:
 > Version:  1.4.0
 > Git hash: ff52401
OMP_NUM_THREADS = 

UCCSDt
******
Options:
 > e_tol:  1e-08
 > t_tol:  1e-08
 > max_iter:  200
 > diis_space:  12
Ansatz: CCSDt
Space: ((7o, 53v), (7o, 53v))
Solving for excitation amplitudes.
Iter   Energy (corr.)        Δ(Energy)    Δ(Amplitudes)
   0    -0.9081471029
Traceback (most recent call last):
  File "/users/k2262473/Calculations/testing/ebcc_CCSDt/input.py", line 16, in <module>
    eb.kernel()
  File "/users/k2262473/dev/ebcc/ebcc/rebcc.py", line 380, in kernel
    amplitudes = self.update_amps(amplitudes=amplitudes, eris=eris)
  File "/users/k2262473/dev/ebcc/ebcc/uebcc.py", line 238, in update_amps
    res = func(**kwargs)
  File "/users/k2262473/dev/ebcc/ebcc/codegen/UCCSDt.py", line 122, in update_amps
    t3new_aaaaaa += einsum(v.aabb.OVov, (0, 1, 2, 3), t3.abaaba[np.ix_(soa,sob,sOfa,sva,svb,sVfa)], (4, 2, 5, 6, 3, 7), (4, 5, 0, 6, 7, 1)) * 2.0
ValueError: operands could not be broadcast together with shapes (7,7,0,53,53,0) (7,0,0,53,0,0) (7,7,0,53,53,0) 

[obackhouse]

Sorry I haven't got round to this yet -- the DF issue is intended behaviour. If there is a with_df option, ebcc will try to use a DF-optimised ansatz. If you want to use a DF mean-field with the non-DF-optimised ansatzes, then either make a new mean-field with the MOs, or do something like

mf._eri = np.dot(with_df._cderi.T, with_df._cderi)
del mf.with_df

[abhishekkhedkar09]

Sorry I haven't got round to this yet -- the DF issue is intended behaviour. If there is a with_df option, ebcc will try to use a DF-optimised ansatz. If you want to use a DF mean-field with the non-DF-optimised ansatzes, then either make a new mean-field with the MOs, or do something like

mf._eri = np.dot(with_df._cderi.T, with_df._cderi)
del mf.with_df

Thanks!

I did the following and it worked.

with_df = mf.with_df
_eri = dot(with_df._cderi.T,with_df._cderi)

newmf = pyscf.scf.UHF(mol)
newmf.mo_coeff = mf.mo_coeff
newmf.mo_occ = mf.mo_occ
newmf._eri = _eri

cc_  = EBCC(newmf, ... )

Just to confirm, occ, coeff, and eri suffice.. right?

[obackhouse]

you'll probably need e_tot for energies - i don't think mo_energy is actually needed but not sure. I'd just be complete and assigned mo_energy, mo_coeff, mo_occ, e_tot, converged.

the easiest way is probably just

from pyscf import lib

with lib.temporary_env(mf, with_df=None, _eri=np.dot(with_df._cderi.T, with_df._cderi):
    cc = EBCC(mf)
    cc.kernel()

[obackhouse]

@abhishekkhedkar09 can you see if everything is working on 16998b4fd4da8471318f58a6352231ff038f18f0?

[abhishekkhedkar09]

@abhishekkhedkar09 can you see if everything is working on 16998b4?

I still get the following error

UCCSDt
******
Options:
 > e_tol:  1e-08
 > t_tol:  1e-08
 > max_iter:  200
 > diis_space:  12
Ansatz: CCSDt
Space: ((5o, 14v) [(1o, 1v) active], (5o, 14v) [(1o, 1v) active])
Solving for excitation amplitudes.
Iter   Energy (corr.)        Δ(Energy)    Δ(Amplitudes)
   0    -0.2113674627
Traceback (most recent call last):
  File "/users/k2262473/Calculations/testing/ebcc_CCSDt/ebcc-examples-folder/input.py", line 30, in <module>
    ccsdt.kernel()
  File "/users/k2262473/dev/ebcc/ebcc/rebcc.py", line 380, in kernel
    amplitudes = self.update_amps(amplitudes=amplitudes, eris=eris)
  File "/users/k2262473/dev/ebcc/ebcc/uebcc.py", line 238, in update_amps
    res = func(**kwargs)
  File "/users/k2262473/dev/ebcc/ebcc/codegen/UCCSDt.py", line 118, in update_amps
    t3new_aaaaaa[np.ix_(soa,soa,sOa,sva,sva,sVa)] += einsum(f.aa.OO, (0, 1), t3.aaaaaa[np.ix_(soa,soa,sOfa,sva,sva,sVfa)], (2, 3, 1, 4, 5, 6), (2, 3, 0, 4, 5, 6)) * -6.0
IndexError: index 4 is out of bounds for axis 2 with size 1

here's my input based on the RCCSDt example in ebcc, perhaps to check for correctness as well.

import numpy as np
from pyscf import gto, scf

from ebcc import UEBCC, Space

mol = gto.Mole()
mol.atom = "H 0 0 0; F 0 0 1.1"
mol.basis = "cc-pvdz"
mol.build()
mf = scf.UHF(mol)
mf.kernel()

occupied = mf.mo_occ > 0
active = np.zeros_like(occupied)

active[0][np.sum(mol.nelectron) // 2 - 1] = True  # HOMO
active[1][np.sum(mol.nelectron) // 2 - 1] = True  # HOMO
active[0][np.sum(mol.nelectron) // 2] = True      # LUMO
active[1][np.sum(mol.nelectron) // 2] = True      # LUMO

space = tuple(Space(o > 0, np.zeros_like(a), a) for o, a in zip(mf.mo_occ, active))

# Run a CCSDt calculation
ccsdt = UEBCC(mf, ansatz="CCSDt", space=space)
ccsdt.kernel()

[obackhouse]

it's running on cec2afd420d879b1203ac6c67ad8d5fc3c24202f - can you check correctness and the vayesta interface?

[obackhouse]

it's running on cec2afd - can you check correctness and the vayesta interface?

@abhishekkhedkar09 did you see this?

[abhishekkhedkar09]

Yes!

I'm sorry for not getting back to you sooner. Thank you for the fix! I did test it today, also with the Vayesta interface. I confirm it all works fine. I will add an example and test in Vayesta.

it's running on cec2afd - can you check correctness and the vayesta interface?

@abhishekkhedkar09 did you see this?

[obackhouse]

cool no rush -- just checking. I'll make a PR.

[obackhouse]

Fixed in #39