adressed Leah's and Minsik's comments

introduced explict `unused` function

src/kbe/helper.py

@@ -5,6 +5,8 @@
 import numpy
 from pyscf import scf
 
+from molbe.helper import unused
+
 
 def get_veff(eri_, dm, S, TA, hf_veff, return_veff0=False):
     """
@@ -30,7 +32,8 @@ def get_veff(eri_, dm, S, TA, hf_veff, return_veff0=False):
     """
 
     # construct rdm
-    nk, _, neo = TA.shape
+    nk, nao, neo = TA.shape
+    unused(nao)
     P_ = numpy.zeros((neo, neo), dtype=numpy.complex128)
     for k in range(nk):
         Cinv = numpy.dot(TA[k].conj().T, S[k])

src/kbe/lo.py

@@ -7,8 +7,8 @@
 from functools import reduce
 
 import numpy
-import scipy.linalg
 from libdmet.lo import pywannier90
+from numpy.linalg import eigh, svd
 
 from kbe.lo_k import (
     get_iao_k,
@@ -18,7 +18,7 @@
     symm_orth_k,
 )
 from molbe.external.lo_helper import get_aoind_by_atom, reorder_by_atom_
-from molbe.helper import ncore_
+from molbe.helper import ncore_, unused
 
 # iao tmp
 
@@ -88,7 +88,7 @@ def localize(
             cinv_ = numpy.zeros((nk, nmo, nao), dtype=numpy.complex128)
 
         for k in range(self.nkpt):
-            es_, vs_ = scipy.linalg.eigh(self.S[k])
+            es_, vs_ = eigh(self.S[k])
             edx = es_ > 1.0e-14
 
             W[k] = numpy.dot(vs_[:, edx] / numpy.sqrt(es_[edx]), vs_[:, edx].conj().T)
@@ -110,7 +110,7 @@ def localize(
 
                 S_ = functools.reduce(numpy.dot, (C_.conj().T, self.S[k], C_))
 
-                es_, vs_ = scipy.linalg.eigh(S_)
+                es_, vs_ = eigh(S_)
                 edx = es_ > 1.0e-14
                 W_ = numpy.dot(vs_[:, edx] / numpy.sqrt(es_[edx]), vs_[:, edx].conj().T)
                 W_nocore[k] = numpy.dot(C_, W_)
@@ -268,7 +268,7 @@ def localize(
                     numpy.dot, (Ciao_[k].conj().T, self.FOCK[k], Ciao_[k])
                 )
                 S_iao = reduce(numpy.dot, (Ciao_[k].conj().T, self.S[k], Ciao_[k]))
-                e_iao, _ = scipy.linalg.eigh(fock_iao, S_iao) # e_iao, v_iao
+                e_iao, _ = eigh(fock_iao, S_iao)  # e_iao, v_iao
                 mo_energy_.append(e_iao)
             iaomf = KMF(self.mol, kpts=self.kpts, mo_coeff=Ciao_, mo_energy=mo_energy_)
 
@@ -363,9 +363,10 @@ def localize(
                     self.Nocc - self.ncore,
                 )
                 # Find virtual orbitals that lie in the span of LOs
-                _, l, vt = numpy.linalg.svd(
+                u, l, vt = svd(
                     self.W[k].conj().T @ self.S[k] @ Cv[k], full_matrices=False
                 )
+                unused(u)
                 nvlo = nlo - self.Nocc - self.ncore
                 assert numpy.allclose(numpy.sum(l[:nvlo]), nvlo)
                 C_ = numpy.hstack([Co_nocore[k], Cv[k] @ vt[:nvlo].conj().T])
@@ -393,7 +394,7 @@ def localize(
         lorb = numpy.zeros((nk, nao, nmo), dtype=numpy.complex128)
         lorb_nocore = numpy.zeros((nk, nao, nmo - self.ncore), dtype=numpy.complex128)
         for k in range(nk):
-            es_, vs_ = scipy.linalg.eigh(self.S[k])
+            es_, vs_ = eigh(self.S[k])
             edx = es_ > 1.0e-14
             lorb[k] = numpy.dot(
                 vs_[:, edx] / numpy.sqrt(es_[edx]), vs_[:, edx].conj().T
@@ -414,7 +415,8 @@ def localize(
                 S_lnc = reduce(
                     numpy.dot, (lorb_nocore[k].conj().T, self.S[k], lorb_nocore[k])
                 )
-                e__, _ = scipy.linalg.eigh(fock_lnc, S_lnc)
+                e__, v__ = eigh(fock_lnc, S_lnc)
+                unused(v__)
                 mo_energy_nc.append(e__)
             lmf = KMF(
                 self.mol, kpts=self.kpts, mo_coeff=lorb_nocore, mo_energy=mo_energy_nc

src/kbe/lo_k.py

@@ -8,6 +8,8 @@
 import scipy
 from pyscf.pbc import gto as pgto
 
+from molbe.helper import unused
+
 
 def dot_gen(A, B, ovlp):
     if ovlp is None:
@@ -111,7 +113,8 @@ def get_iao_k(Co, S12, S1, S2=None, ortho=True):
         S2: valence AO ovlp
     """
 
-    nk, nao, _ = S12.shape
+    nk, nao, nmo = S12.shape
+    unused(nmo)
     P1 = numpy.zeros_like(S1, dtype=numpy.complex128)
     P2 = numpy.zeros_like(S2, dtype=numpy.complex128)
 
@@ -160,7 +163,8 @@ def get_pao_k(Ciao, S, S12, S2):
         Cpao (orthogonalized)
     """
 
-    nk, nao, _ = Ciao.shape
+    nk, nao, niao = Ciao.shape
+    unused(niao)
     Cpao = []
     for k in range(nk):
         s12 = scipy.linalg.inv(S[k]) @ S12[k]
@@ -172,9 +176,7 @@ def get_pao_k(Ciao, S, S12, S2):
         numpy.o0 = cpao_.shape[-1]
         Cpao.append(cano_orth(cpao_, ovlp=S[k]))
         numpy.o1 = Cpao[k].shape[-1]
-    Cpao = numpy.asarray(Cpao)
-
-    return Cpao
+    return numpy.asarray(Cpao)
 
 
 def get_pao_native_k(Ciao, S, mol, valence_basis, kpts, ortho=True):

src/kbe/pbe.py

@@ -230,9 +230,7 @@ def __init__(
             if not restart:
                 self.Nocc -= self.ncore
 
-                nk, nao, nao_2nd_dim = self.hf_dm.shape
-                assert nao == nao_2nd_dim
-                del nao_2nd_dim
+                nk, nao = self.hf_dm.shape[:2]
 
                 dm_nocore = numpy.zeros(
                     (nk, nao, nao), dtype=numpy.result_type(self.C, self.C)
@@ -353,9 +351,7 @@ def print_ini(self):
 
     def ewald_sum(self, kpts=None):
         dm_ = self.mf.make_rdm1()
-        nk, nao, nao_2nd_dim = dm_.shape
-        assert nao == nao_2nd_dim
-        del nao_2nd_dim
+        nk, nao = dm_.shape[:2]
 
         vk_kpts = numpy.zeros(dm_.shape) * 1j
         _ewald_exxdiv_for_G0(

src/kbe/pfrag.py

@@ -9,7 +9,7 @@
 from kbe.helper import get_veff
 from kbe.misc import get_phase, get_phase1
 from kbe.solver import schmidt_decomp_svd
-from molbe.helper import get_eri, get_scfObj
+from molbe.helper import get_eri, get_scfObj, unused
 
 
 class Frags:
@@ -207,7 +207,8 @@ def cons_h1(self, h1):
             One-electron Hamiltonian matrix.
         """
 
-        nk, _, teo = self.TA.shape
+        nk, nao, teo = self.TA.shape
+        unused(nao)
         h1_eo = numpy.zeros((teo, teo), dtype=numpy.complex128)
         for k in range(nk):
             h1_eo += functools.reduce(

src/kbe/solver.py

@@ -3,6 +3,8 @@
 import numpy
 import scipy.linalg
 
+from molbe.helper import unused
+
 
 def schmidt_decomp_svd(rdm, Frag_sites):
     """
@@ -33,7 +35,8 @@ def schmidt_decomp_svd(rdm, Frag_sites):
     nfs = len(Frag_sites)
 
     Denv = rdm[Env_sites1][:, Fragsites]
-    U, sigma, _ = scipy.linalg.svd(Denv, full_matrices=False, lapack_driver="gesvd")
+    U, sigma, V = scipy.linalg.svd(Denv, full_matrices=False, lapack_driver="gesvd")
+    unused(V)
     nbath = (sigma >= thres).sum()
     TA = numpy.zeros((Tot_sites, nfs + nbath), dtype=numpy.complex128)
     TA[Fragsites, :nfs] = numpy.eye(nfs)

src/molbe/autofrag.py

@@ -4,7 +4,7 @@
 
 import numpy
 
-from molbe.helper import get_core
+from molbe.helper import get_core, unused
 
 
 def autogen(
@@ -88,7 +88,8 @@ def autogen(
         cell.basis = valence_basis
         cell.build()
 
-    _, _, core_list = get_core(cell)
+    ncore, no_core_idx, core_list = get_core(cell)
+    unused(ncore, no_core_idx)
     coord = cell.atom_coords()
     ang2bohr = 1.88973
     normdist = 3.5 * ang2bohr

src/molbe/be_parallel.py

@@ -10,7 +10,7 @@
 
 from molbe.external.ccsd_rdm import make_rdm1_uccsd, make_rdm2_uccsd
 from molbe.external.unrestricted_utils import make_uhf_obj
-from molbe.helper import get_eri, get_frag_energy, get_frag_energy_u, get_scfObj
+from molbe.helper import get_eri, get_frag_energy, get_frag_energy_u, get_scfObj, unused
 from molbe.solver import (
     make_rdm1_ccsd_t1,
     make_rdm2_urlx,
@@ -129,7 +129,8 @@ def run_solver(
             )
     elif solver == "FCI":
         mc_ = fci.FCI(mf_, mf_.mo_coeff)
-        _, civec = mc_.kernel()
+        efci, civec = mc_.kernel()
+        unused(efci)
         rdm1_tmp = mc_.make_rdm1(civec, mc_.norb, mc_.nelec)
     elif solver == "HCI":
         # pylint: disable-next=E0611
@@ -151,7 +152,8 @@ def run_solver(
 
         nelec = (nocc, nocc)
         h1_ = functools.reduce(numpy.dot, (mf_.mo_coeff.T, h1, mf_.mo_coeff))
-        _, civec = ci_.kernel(h1_, eri, nmo, nelec)
+        eci, civec = ci_.kernel(h1_, eri, nmo, nelec)
+        unused(eci)
         civec = numpy.asarray(civec)
 
         (rdm1a_, rdm1b_), (rdm2aa, rdm2ab, rdm2bb) = ci_.make_rdm12s(civec, nmo, nelec)
@@ -184,7 +186,8 @@ def run_solver(
         nao, nmo = mf_.mo_coeff.shape
         nelec = (nocc, nocc)
         cas = mcscf.CASCI(mf_, nmo, nelec)
-        h1, _ = cas.get_h1eff(mo_coeff=mf_.mo_coeff)
+        h1, ecore = cas.get_h1eff(mo_coeff=mf_.mo_coeff)
+        unused(ecore)
         eri = ao2mo.kernel(mf_._eri, mf_.mo_coeff, aosym="s4", compact=False).reshape(
             4 * ((nmo),)
         )

src/molbe/helper.py

@@ -517,3 +517,8 @@ def contract_2e(jmaxs, rdm2_, V_, s, sym):
             ec_tmp += efac[s][0] * ec[s][i]
 
     return [e1_tmp, e2_tmp, ec_tmp]
+
+
+def unused(*args) -> None:
+    for arg in args:
+        del arg

src/molbe/ube.py

@@ -20,6 +20,7 @@
 
 import molbe.be_var as be_var
 from molbe.be_parallel import be_func_parallel_u
+from molbe.helper import unused
 from molbe.mbe import BE
 from molbe.pfrag import Frags
 from molbe.solver import be_func_u
@@ -323,9 +324,10 @@ def initialize(self, eri_, compute_hf):
             )
 
             if compute_hf:
-                eh1_a, ecoul_a, _ = fobj_a.energy_hf(
+                eh1_a, ecoul_a, ef_a = fobj_a.energy_hf(
                     return_e1=True, unrestricted=True, spin_ind=0
                 )
+                unused(ef_a)
                 EH1 += eh1_a
                 ECOUL += ecoul_a
                 E_hf += fobj_a.ebe_hf
@@ -345,9 +347,10 @@ def initialize(self, eri_, compute_hf):
             )
 
             if compute_hf:
-                eh1_b, ecoul_b, _ = fobj_b.energy_hf(
+                eh1_b, ecoul_b, ef_b = fobj_b.energy_hf(
                     return_e1=True, unrestricted=True, spin_ind=1
                 )
+                unused(ef_b)
                 EH1 += eh1_b
                 ECOUL += ecoul_b
                 E_hf += fobj_b.ebe_hf
@@ -414,7 +417,7 @@ def oneshot(
         self, solver="UCCSD", nproc=1, ompnum=4, calc_frag_energy=False, clean_eri=False
     ):
         if nproc == 1:
-            E, _ = be_func_u(
+            E, E_comp = be_func_u(
                 None,
                 zip(self.Fobjs_a, self.Fobjs_b),
                 solver,
@@ -427,7 +430,7 @@ def oneshot(
                 frozen=self.frozen_core,
             )
         else:
-            E, _ = be_func_parallel_u(
+            E, E_comp = be_func_parallel_u(
                 None,
                 zip(self.Fobjs_a, self.Fobjs_b),
                 solver,
@@ -441,6 +444,7 @@ def oneshot(
                 nproc=nproc,
                 ompnum=ompnum,
             )
+        unused(E_comp)
 
         print("-----------------------------------------------------", flush=True)
         print("             One Shot BE ", flush=True)