diff --git a/.ruff.toml b/.ruff.toml
index b94483a7..708d0329 100644
--- a/.ruff.toml
+++ b/.ruff.toml
@@ -11,3 +11,4 @@ ignore = [
         # Prevents the use of the characters 'l', 'O', or 'I' as variable names.
         # Overly restrictive, in particular when implementing mathematical expressions.
 ]
+preview = true
diff --git a/src/kbe/__init__.py b/src/kbe/__init__.py
index a532531e..186920fd 100644
--- a/src/kbe/__init__.py
+++ b/src/kbe/__init__.py
@@ -1,2 +1,4 @@
 from kbe.fragment import fragpart
 from kbe.pbe import BE
+
+__all__ = ["fragpart", "BE"]
diff --git a/src/kbe/_opt.py b/src/kbe/_opt.py
index b996e769..2b622d94 100644
--- a/src/kbe/_opt.py
+++ b/src/kbe/_opt.py
@@ -2,6 +2,7 @@
 
 import sys
 
+from kbe.misc import print_energy
 from molbe._opt import BEOPT
 
 
@@ -53,8 +54,6 @@ def optimize(
     J0 : list of list of float
        Initial Jacobian.
     """
-    from .misc import print_energy
-
     # Check if only chemical potential optimization is required
     if not only_chem:
         pot = self.pot
diff --git a/src/kbe/autofrag.py b/src/kbe/autofrag.py
index e90beeca..8c20a98a 100644
--- a/src/kbe/autofrag.py
+++ b/src/kbe/autofrag.py
@@ -3,7 +3,9 @@
 import sys
 
 import numpy
+from pyscf import lib
 
+from kbe.misc import sgeom
 from molbe.helper import get_core
 
 
@@ -324,10 +326,6 @@ def autogen(
         Weights for each fragment. Each entry contains a weight and a list of LO
         indices.
     """
-    from pyscf import lib
-
-    from .misc import sgeom
-
     if not float(unitcell).is_integer():
         print("Fractional unitcell is not supported!")
         sys.exit()
diff --git a/src/kbe/chain.py b/src/kbe/chain.py
index c84986bc..650ce988 100644
--- a/src/kbe/chain.py
+++ b/src/kbe/chain.py
@@ -3,10 +3,10 @@
 
 # Old fragmentation code - do not use other than debugging!
 
+import math
 
-def findH(mol, nh, tmphlist=[]):
-    import math
 
+def findH(mol, nh, tmphlist=[]):
     coord_ = mol.atom_coords()
 
     hidx = []
diff --git a/src/kbe/fragment.py b/src/kbe/fragment.py
index 80c7342e..2f659a7d 100644
--- a/src/kbe/fragment.py
+++ b/src/kbe/fragment.py
@@ -2,10 +2,9 @@
 
 import sys
 
+from kbe.autofrag import autogen
 from molbe.helper import get_core
 
-from .autofrag import autogen
-
 
 def print_mol_missing():
     print("Provide pyscf gto.Cell object in fragpart() and restart!", flush=True)
@@ -104,7 +103,7 @@ def __init__(
         self.allcen = allcen
         self.valence_basis = valence_basis
         self.kpt = kpt
-        self.molecule = False  ### remove this
+        self.molecule = False  # remove this
 
         # Check for frozen core approximation
         if frozen_core:
@@ -156,4 +155,6 @@ def __init__(
             print("exiting", flush=True)
             sys.exit()
 
-    from .chain import polychain
+    # This import makes polychain a method of the class and
+    #  cannot be moved to the top of the file
+    from kbe.chain import polychain  # noqa: PLC0415
diff --git a/src/kbe/helper.py b/src/kbe/helper.py
index 09e52097..a6e1b0f8 100644
--- a/src/kbe/helper.py
+++ b/src/kbe/helper.py
@@ -3,11 +3,10 @@
 import functools
 
 import numpy
+from pyscf import scf
 
 
 def get_veff(eri_, dm, S, TA, hf_veff, return_veff0=False):
-    from pyscf import scf
-
     """
     Calculate the effective HF potential (Veff) for a given density matrix
     and electron repulsion integrals.
diff --git a/src/kbe/lo.py b/src/kbe/lo.py
index 192b0777..6fedb920 100644
--- a/src/kbe/lo.py
+++ b/src/kbe/lo.py
@@ -1,12 +1,24 @@
 # Author(s): Oinam Meitei
 #            Henry Tran
 #
+import functools
+import os
 import sys
 from functools import reduce
 
 import numpy
-
+import scipy.linalg
+from libdmet.lo import pywannier90
+
+from kbe.lo_k import (
+    get_iao_k,
+    get_pao_native_k,
+    get_xovlp_k,
+    remove_core_mo_k,
+    symm_orth_k,
+)
 from molbe.external.lo_helper import get_aoind_by_atom, reorder_by_atom_
+from molbe.helper import ncore_
 
 # iao tmp
 
@@ -51,17 +63,11 @@ def localize(
     iao_wannier : bool
        Whether to perform Wannier localization in the IAO space
     """
-    import functools
-
-    import scipy.linalg
-
-    from molbe.helper import ncore_
-
     if lo_method == "iao":
         if valence_basis == "sto-3g":
-            from .basis_sto3g_core_val import core_basis, val_basis
+            from .basis_sto3g_core_val import core_basis, val_basis  # noqa: PLC0415
         elif valence_basis == "minao":
-            from .basis_minao_core_val import core_basis, val_basis
+            from .basis_minao_core_val import core_basis, val_basis  # noqa: PLC0415
         elif iao_val_core:
             sys.exit(
                 "valence_basis="
@@ -128,18 +134,6 @@ def localize(
         self.cinv = cinv_
 
     elif lo_method == "iao":
-        import os
-
-        from libdmet.lo import pywannier90
-
-        from .lo_k import (
-            get_iao_k,
-            get_pao_native_k,
-            get_xovlp_k,
-            remove_core_mo_k,
-            symm_orth_k,
-        )
-
         if not iao_val_core or not self.frozen_core:
             Co = self.C[:, :, : self.Nocc].copy()
             S12, S2 = get_xovlp_k(self.cell, self.kpts, basis=valence_basis)
@@ -395,11 +389,6 @@ def localize(
         self.cinv = cinv_
 
     elif lo_method == "wannier":
-        # from pyscf.pbc.tools import pywannier90
-        from libdmet.lo import pywannier90
-
-        from .lo_k import remove_core_mo_k
-
         nk, nao, nmo = self.C.shape
         lorb = numpy.zeros((nk, nao, nmo), dtype=numpy.complex128)
         lorb_nocore = numpy.zeros((nk, nao, nmo - self.ncore), dtype=numpy.complex128)
diff --git a/src/kbe/lo_k.py b/src/kbe/lo_k.py
index aba3bd6e..ea0da688 100644
--- a/src/kbe/lo_k.py
+++ b/src/kbe/lo_k.py
@@ -6,6 +6,7 @@
 
 import numpy
 import scipy
+from pyscf.pbc import gto as pgto
 
 
 def dot_gen(A, B, ovlp):
@@ -67,9 +68,6 @@ def get_xovlp_k(cell, kpts, basis="sto-3g"):
         S12 - Overlap of two basis sets
         S22 - Overlap in new basis set
     """
-
-    from pyscf.pbc import gto as pgto  # intor_cross
-
     cell_alt = cell.copy()
     cell_alt.basis = basis
     cell_alt.build()
diff --git a/src/kbe/misc.py b/src/kbe/misc.py
index ff50c2b2..53efa3b7 100644
--- a/src/kbe/misc.py
+++ b/src/kbe/misc.py
@@ -2,6 +2,7 @@
 
 import numpy
 from pyscf import lib
+from pyscf.pbc import tools
 
 
 def sgeom(cell, kmesh=None):
@@ -15,11 +16,7 @@ def sgeom(cell, kmesh=None):
     kmesh : list of int
        Number of k-points in each lattice vector dimension
     """
-    from pyscf.pbc import tools
-
-    scell = tools.super_cell(cell, kmesh)
-
-    return scell
+    return tools.super_cell(cell, kmesh)
 
 
 def get_phase(cell, kpts, kmesh):
diff --git a/src/kbe/pbe.py b/src/kbe/pbe.py
index 9e76c0ee..4db4a694 100644
--- a/src/kbe/pbe.py
+++ b/src/kbe/pbe.py
@@ -3,14 +3,18 @@
 import os
 import pickle
 import sys
+from multiprocessing import Pool
 
 import h5py
 import numpy
+from libdmet.basis_transform.eri_transform import get_emb_eri_fast_gdf
+from pyscf import ao2mo
+from pyscf.pbc import df, gto
+from pyscf.pbc.df.df_jk import _ewald_exxdiv_for_G0
 
 import molbe.be_var as be_var
-
-from .misc import storePBE
-from .pfrag import Frags
+from kbe.misc import storePBE
+from kbe.pfrag import Frags
 
 
 class BE:
@@ -322,11 +326,11 @@ def __init__(
         if not restart:
             self.initialize(mf._eri, compute_hf)
 
-    # this is a molbe method not BEOPT
-    from molbe.external.optqn import get_be_error_jacobian
-
-    from ._opt import optimize
-    from .lo import localize
+    # The following import of these functions turns them into
+    # proper methods of the class.
+    from kbe._opt import optimize  # noqa: PLC0415
+    from kbe.lo import localize  # noqa: PLC0415
+    from molbe.external.optqn import get_be_error_jacobian  # noqa: PLC0415
 
     def print_ini(self):
         """
@@ -349,8 +353,6 @@ def print_ini(self):
         print(flush=True)
 
     def ewald_sum(self, kpts=None):
-        from pyscf.pbc.df.df_jk import _ewald_exxdiv_for_G0
-
         dm_ = self.mf.make_rdm1()
         nk, nao, nao = dm_.shape
 
@@ -380,13 +382,6 @@ def initialize(self, eri_, compute_hf, restart=False):
         restart : bool, optional
             Whether to restart from a previous calculation, by default False.
         """
-        import os
-        from multiprocessing import Pool
-
-        import h5py
-        from libdmet.basis_transform.eri_transform import get_emb_eri_fast_gdf
-        from pyscf import ao2mo
-
         if compute_hf:
             E_hf = 0.0
         EH1 = 0.0
@@ -609,8 +604,9 @@ def oneshot(
         clean_eri : bool, optional
             Whether to clean up ERI files after calculation, by default False.
         """
-        from .be_parallel import be_func_parallel
-        from .solver import be_func
+        # has to be here because of circular dependency
+        from .be_parallel import be_func_parallel  # noqa: PLC0415
+        from .solver import be_func  # noqa: PLC0415
 
         print("Calculating Energy by Fragment? ", calc_frag_energy)
         if nproc == 1:
@@ -768,10 +764,6 @@ def eritransform_parallel(a, atom, basis, kpts, C_ao_emb, cderi):
     """
     Wrapper for parallel eri transformation
     """
-    from pyscf.pbc import df, gto
-
-    from molbe.external.eri_transform import get_emb_eri_fast_gdf
-
     cell = gto.Cell()
     cell.a = a
     cell.atom = atom
diff --git a/src/kbe/pfrag.py b/src/kbe/pfrag.py
index eea12d18..bebc018d 100644
--- a/src/kbe/pfrag.py
+++ b/src/kbe/pfrag.py
@@ -6,12 +6,11 @@
 import h5py
 import numpy
 
+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 .helper import *
-from .misc import *
-from .solver import schmidt_decomp_svd
-
 
 class Frags:
     """
@@ -141,9 +140,6 @@ def sd(
         kmesh : list of int
             Number of k-points in each lattice vector direction
         """
-
-        from .misc import get_phase
-
         nk, nao, nlo = lao.shape
         rdm1_lo_k = numpy.zeros((nk, nlo, nlo), dtype=numpy.result_type(lmo, lmo))
         for k in range(nk):
diff --git a/src/kbe/solver.py b/src/kbe/solver.py
index 2cdc5f3d..a30dd068 100644
--- a/src/kbe/solver.py
+++ b/src/kbe/solver.py
@@ -1,6 +1,7 @@
 # Author(s): Oinam Romesh Meitei
 
 import numpy
+import scipy.linalg
 
 
 def schmidt_decomp_svd(rdm, Frag_sites):
@@ -23,8 +24,6 @@ def schmidt_decomp_svd(rdm, Frag_sites):
     numpy.ndarray
         Transformation matrix (TA) including both fragment and entangled bath orbitals.
     """
-    import scipy.linalg
-
     thres = 1.0e-10
     Tot_sites = rdm.shape[0]
 
diff --git a/src/molbe/__init__.py b/src/molbe/__init__.py
index a8df64ff..b5daa904 100644
--- a/src/molbe/__init__.py
+++ b/src/molbe/__init__.py
@@ -1,3 +1,5 @@
 from molbe.fragment import fragpart
 from molbe.mbe import BE
 from molbe.ube import UBE
+
+__all__ = ["fragpart", "BE", "UBE"]
diff --git a/src/molbe/_opt.py b/src/molbe/_opt.py
index d350bbb5..16308be3 100644
--- a/src/molbe/_opt.py
+++ b/src/molbe/_opt.py
@@ -4,8 +4,10 @@
 
 import numpy
 
-from .be_parallel import be_func_parallel
-from .solver import be_func
+from molbe.be_parallel import be_func_parallel
+from molbe.external.optqn import FrankQN
+from molbe.misc import print_energy
+from molbe.solver import be_func
 
 
 class BEOPT:
@@ -182,10 +184,6 @@ def optimize(self, method, J0=None, trust_region=False):
         trust_region : bool, optional
            Use trust-region based QN optimization, by default False
         """
-        import sys
-
-        from molbe.external.optqn import FrankQN
-
         print("-----------------------------------------------------", flush=True)
         print("             Starting BE optimization ", flush=True)
         print("             Solver : ", self.solver, flush=True)
@@ -289,8 +287,6 @@ def optimize(
     trust_region : bool, optional
        Use trust-region based QN optimization, by default False
     """
-    from .misc import print_energy
-
     # Check if only chemical potential optimization is required
     if not only_chem:
         pot = self.pot
diff --git a/src/molbe/autofrag.py b/src/molbe/autofrag.py
index 6bdb1993..622e3499 100644
--- a/src/molbe/autofrag.py
+++ b/src/molbe/autofrag.py
@@ -4,7 +4,7 @@
 
 import numpy
 
-from .helper import get_core
+from molbe.helper import get_core
 
 
 def autogen(
diff --git a/src/molbe/be_parallel.py b/src/molbe/be_parallel.py
index f8a76f98..3c95003d 100644
--- a/src/molbe/be_parallel.py
+++ b/src/molbe/be_parallel.py
@@ -1,16 +1,17 @@
 # Author(s): Oinam Romesh Meitei, Leah Weisburn
 
 import functools
+import os
 import sys
+from multiprocessing import Pool
 
 import numpy
+from pyscf import ao2mo, fci, mcscf
 
 from molbe.external.ccsd_rdm import make_rdm1_uccsd, make_rdm2_uccsd
 from molbe.external.unrestricted_utils import make_uhf_obj
-
-from .helper import *
-from .helper import get_frag_energy
-from .solver import (
+from molbe.helper import get_eri, get_frag_energy, get_frag_energy_u, get_scfObj
+from molbe.solver import (
     make_rdm1_ccsd_t1,
     make_rdm2_urlx,
     solve_ccsd,
@@ -127,13 +128,11 @@ def run_solver(
                 relax=True,
             )
     elif solver == "FCI":
-        from pyscf import fci
-
         mc_ = fci.FCI(mf_, mf_.mo_coeff)
         efci, civec = mc_.kernel()
         rdm1_tmp = mc_.make_rdm1(civec, mc_.norb, mc_.nelec)
     elif solver == "HCI":
-        from pyscf import ao2mo, hci
+        from pyscf import hci  # noqa: PLC0415  # hci is an optional module
 
         nao, nmo = mf_.mo_coeff.shape
         eri = ao2mo.kernel(mf_._eri, mf_.mo_coeff, aosym="s4", compact=False).reshape(
@@ -159,7 +158,7 @@ def run_solver(
         rdm2s = rdm2aa + rdm2ab + rdm2ab.transpose(2, 3, 0, 1) + rdm2bb
 
     elif solver == "SHCI":
-        from pyscf.shciscf import shci
+        from pyscf.shciscf import shci  # noqa: PLC0415    # shci is an optional module
 
         nao, nmo = mf_.mo_coeff.shape
         nelec = (nocc, nocc)
@@ -177,7 +176,9 @@ def run_solver(
         rdm1_tmp, rdm2s = mch.fcisolver.make_rdm12(0, nmo, nelec)
 
     elif solver == "SCI":
-        from pyscf import ao2mo, cornell_shci, mcscf
+        from pyscf import (  # noqa: PLC0415  # cornell_shci is an optional module
+            cornell_shci,
+        )
 
         nao, nmo = mf_.mo_coeff.shape
         nelec = (nocc, nocc)
@@ -469,9 +470,6 @@ def be_func_parallel(
         Depending on the parameters, returns the error norm or a tuple containing
         the error norm, error vector, and the computed energy.
     """
-    import os
-    from multiprocessing import Pool
-
     nfrag = len(Fobjs)
     # Create directories for fragments if required
     if writeh1 and solver == "SCI":
@@ -639,9 +637,6 @@ def be_func_parallel_u(
     float
         Returns the computed energy
     """
-    import os
-    from multiprocessing import Pool
-
     # Set the number of OpenMP threads
     os.system("export OMP_NUM_THREADS=" + str(ompnum))
     nprocs = int(nproc / ompnum)
diff --git a/src/molbe/eri_onthefly.py b/src/molbe/eri_onthefly.py
index c3b6ac71..c5cfa8f5 100644
--- a/src/molbe/eri_onthefly.py
+++ b/src/molbe/eri_onthefly.py
@@ -8,7 +8,7 @@
 from pyscf.gto.moleintor import getints3c, make_cintopt, make_loc
 from scipy.linalg import cholesky, solve_triangular
 
-from . import be_var
+from molbe import be_var
 
 
 def integral_direct_DF(mf, Fobjs, file_eri, auxbasis=None):
diff --git a/src/molbe/external/ccsd_rdm.py b/src/molbe/external/ccsd_rdm.py
index 49ecd2c8..6fed74d1 100644
--- a/src/molbe/external/ccsd_rdm.py
+++ b/src/molbe/external/ccsd_rdm.py
@@ -4,6 +4,7 @@
 #
 
 import numpy
+from pyscf.cc.uccsd_rdm import make_rdm1, make_rdm2
 
 
 def make_rdm1_ccsd_t1(t1):
@@ -53,8 +54,6 @@ def make_rdm2_urlx(t1, t2, with_dm1=True):
 
 
 def make_rdm1_uccsd(ucc, relax=False):
-    from pyscf.cc.uccsd_rdm import make_rdm1
-
     if relax == True:
         rdm1 = make_rdm1(ucc, ucc.t1, ucc.t2, ucc.l1, ucc.l2)
     else:
@@ -65,8 +64,6 @@ def make_rdm1_uccsd(ucc, relax=False):
 
 
 def make_rdm2_uccsd(ucc, relax=False, with_dm1=True):
-    from pyscf.cc.uccsd_rdm import make_rdm2
-
     if relax == True:
         rdm2 = make_rdm2(ucc, ucc.t1, ucc.t2, ucc.l1, ucc.l2, with_dm1=with_dm1)
     else:
diff --git a/src/molbe/external/cpmp2_utils.py b/src/molbe/external/cpmp2_utils.py
index d43a454a..65915554 100644
--- a/src/molbe/external/cpmp2_utils.py
+++ b/src/molbe/external/cpmp2_utils.py
@@ -5,7 +5,10 @@
 
 import numpy as np
 import scipy.linalg as slg
-from pyscf import ao2mo
+from pyscf import ao2mo, scf
+
+from molbe.external.cphf_utils import cphf_kernel_batch as cphf_kernel
+from molbe.external.cphf_utils import get_cpuhf_u_batch as cpuhf_kernel
 
 """ RMP2 implementation
 """
@@ -17,14 +20,10 @@ def get_Diajb_r(moe, no):
     eo = moe[:no]
     ev = moe[no:]
     Dia = (eo.reshape(-1, 1) - ev).ravel()
-    Diajb = (Dia.reshape(-1, 1) + Dia).reshape(no, nv, no, nv)
-
-    return Diajb
+    return (Dia.reshape(-1, 1) + Dia).reshape(no, nv, no, nv)
 
 
 def get_dF_r(no, V, C, Q, u):
-    from pyscf import scf
-
     n = C.shape[0]
     nv = n - no
     Co = C[:, :no]
@@ -33,15 +32,11 @@ def get_dF_r(no, V, C, Q, u):
     dP = -Co @ uov @ Cv.T
     dP += dP.T
     vj, vk = scf.hf.dot_eri_dm(V, dP * 2.0, hermi=1)
-    dF = Q + vj - 0.5 * vk
-
-    return dF
+    return Q + vj - 0.5 * vk
 
 
 def get_dmoe_F_r(C, dF):
-    de = np.einsum("pi,qi,pq->i", C, C, dF)
-
-    return de
+    return np.einsum("pi,qi,pq->i", C, C, dF)
 
 
 def get_full_u_F_r(no, C, moe, dF, u):
@@ -108,8 +103,6 @@ def get_dPmp2_batch_r(C, moe, V, no, Qs, aorep=True):
     Diajb = get_Diajb_r(moe, no)
     t2 = Vovov / Diajb
 
-    from .cphf_utils import cphf_kernel_batch as cphf_kernel
-
     us = cphf_kernel(C, moe, V, no, Qs)
     nQ = len(Qs)
     dPs = [None] * nQ
@@ -311,9 +304,7 @@ def get_dPmp2_batch_u(C, moe, V, no, Qs, aorep=True):
     Diajb = get_Diajb_u(moe, no)
     t2 = [Vovov[s] / Diajb[s] for s in range(3)]
 
-    from .cphf_utils import get_cpuhf_u_batch as cphf_kernel
-
-    us = cphf_kernel(C, moe, V, no, Qs)
+    us = cpuhf_kernel(C, moe, V, no, Qs)
     nQ = len(Qs)
     dPs = [None] * nQ
     Phf = [np.diag([1 if i < no[s] else 0 for i in range(n[s])]) for s in [0, 1]]
diff --git a/src/molbe/external/jac_utils.py b/src/molbe/external/jac_utils.py
index 20cee9c6..1cf8b18a 100644
--- a/src/molbe/external/jac_utils.py
+++ b/src/molbe/external/jac_utils.py
@@ -6,6 +6,9 @@
 import numpy as np
 from pyscf import ao2mo
 
+from molbe.external.cphf_utils import cphf_kernel_batch
+from molbe.external.cpmp2_utils import get_dF_r, get_Diajb_r, get_dmoe_F_r
+
 """ Derivative of approximate t1 amplitudes
 t_ia = ((2*t2-t2)_ibjc g_cjba - g_ikbj (2*t2-t2)_jbka) / (e_i - e_a)
 This approximate t1 is by substituting the MP2 t2 amplitudes into the CCSD equation and
@@ -99,8 +102,6 @@ def get_Dia_r(moe_, no_):
         return moe_[:no_].reshape(-1, 1) - moe_[no_:]
 
     # prepare integrals
-    from .cpmp2_utils import get_Diajb_r
-
     Vovov = get_Vmogen_r(no, V, C, "ovov")
     Vvovv = get_Vmogen_r(no, V, C, "vovv")
     Voovo = get_Vmogen_r(no, V, C, "oovo")
@@ -111,16 +112,12 @@ def get_Dia_r(moe_, no_):
 
     # solve CPHF get u
     if us is None:
-        from .cphf_utils import cphf_kernel_batch
-
         us = cphf_kernel_batch(C, moe, V, no, Qs)
 
     dt1s = [None] * len(Qs)
     iu = 0
     for u, Q in zip(us, Qs):
         # get A
-        from .cpmp2_utils import get_dF_r
-
         A = -get_dF_r(no, V, C, Q, u)
         Aoo = Co.T @ A @ Co
         Avv = Cv.T @ A @ Cv
@@ -137,8 +134,6 @@ def get_Dia_r(moe_, no_):
         dVoovo = get_dVmogen_r(no, V, C, u, "oovo")
 
         # get dmoe and deov
-        from .cpmp2_utils import get_dmoe_F_r
-
         dmoe = get_dmoe_F_r(C, -A)
         deov = get_Dia_r(dmoe, no)
 
@@ -170,8 +165,6 @@ def get_Dia_r(moe_, no_):
 
 def get_dPccsdurlx_batch_u(C, moe, eri, no, vpots):
     # cphf
-    from .cphf_utils import cphf_kernel_batch
-
     us = cphf_kernel_batch(C, moe, eri, no, vpots)
     # get mp2 part
     dPt1aos = get_dt1ao_an(no, eri, C, moe, vpots, us=us)
diff --git a/src/molbe/external/lo_helper.py b/src/molbe/external/lo_helper.py
index 4d5dad11..43888eb3 100644
--- a/src/molbe/external/lo_helper.py
+++ b/src/molbe/external/lo_helper.py
@@ -48,8 +48,6 @@ def get_aoind_by_atom(mol, atomind_by_motif=None):
 
 
 def reorder_by_atom_(Clo, aoind_by_atom, S, thr=0.5):
-    import numpy as np
-
     natom = len(aoind_by_atom)
     nlo = Clo.shape[1]
     X = get_symm_mat_pow(S, 0.5)
@@ -61,8 +59,8 @@ def reorder_by_atom_(Clo, aoind_by_atom, S, thr=0.5):
     loshift = 0
     for ia in range(natom):
         ra = aoind_by_atom[ia]
-        poplo_by_atom = np.sum(Clo_soao[ra] ** 2.0, axis=0)
-        loind_a = np.where(poplo_by_atom > thr)[0].tolist()
+        poplo_by_atom = numpy.sum(Clo_soao[ra] ** 2.0, axis=0)
+        loind_a = numpy.where(poplo_by_atom > thr)[0].tolist()
         loind_reorder += loind_a
         nlo_a = len(loind_a)
         loind_by_atom[ia] = list(range(loshift, loshift + nlo_a))
diff --git a/src/molbe/external/optqn.py b/src/molbe/external/optqn.py
index 143dff70..91273c13 100644
--- a/src/molbe/external/optqn.py
+++ b/src/molbe/external/optqn.py
@@ -8,7 +8,11 @@
 
 import numpy
 
-from .. import be_var
+from molbe import be_var
+from molbe.external.cphf_utils import cphf_kernel_batch, get_rhf_dP_from_u
+from molbe.external.cpmp2_utils import get_dPmp2_batch_r
+from molbe.external.jac_utils import get_dPccsdurlx_batch_u
+from molbe.helper import get_eri, get_scfObj
 
 
 def line_search_LF(func, xold, fold, dx, iter_):
@@ -307,8 +311,6 @@ def get_be_error_jacobian(self, jac_solver="HF"):
 
 
 def get_atbe_Jblock_frag(fobj, res_func):
-    from molbe.helper import get_eri, get_scfObj
-
     vpots = get_vpots_frag(fobj.nao, fobj.edge_idx, fobj.fsites)
     eri_ = get_eri(fobj.dname, fobj.nao, eri_file=fobj.eri_file)
     dm0 = (
@@ -416,8 +418,6 @@ def get_be_error_jacobian_selffrag(self, jac_solver="HF"):
 
 
 def hfres_func(mf, vpots, eri, nsocc):
-    from molbe.external.cphf_utils import cphf_kernel_batch, get_rhf_dP_from_u
-
     C = mf.mo_coeff
     moe = mf.mo_energy
     eri = mf._eri
@@ -431,8 +431,6 @@ def hfres_func(mf, vpots, eri, nsocc):
 
 
 def mp2res_func(mf, vpots, eri, nsocc):
-    from molbe.external.cpmp2_utils import get_dPmp2_batch_r
-
     C = mf.mo_coeff
     moe = mf.mo_energy
     eri = mf._eri
@@ -446,8 +444,6 @@ def mp2res_func(mf, vpots, eri, nsocc):
 
 
 def ccsdres_func(mf, vpots, eri, nsocc):
-    from molbe.external.jac_utils import get_dPccsdurlx_batch_u
-
     C = mf.mo_coeff
     moe = mf.mo_energy
     eri = mf._eri
diff --git a/src/molbe/external/unrestricted_utils.py b/src/molbe/external/unrestricted_utils.py
index 17284c37..11b85635 100644
--- a/src/molbe/external/unrestricted_utils.py
+++ b/src/molbe/external/unrestricted_utils.py
@@ -1,14 +1,14 @@
 # Authors: Leah Weisburn, Hongzhou Ye, Henry Tran
 
+import ctypes
 import functools
 
 import h5py
 import numpy
+from pyscf import ao2mo, lib, scf
 
 
 def make_uhf_obj(fobj_a, fobj_b, frozen=False):
-    from pyscf import scf
-
     """
     Constructs UHF object from the alpha and beta components
     """
@@ -57,8 +57,6 @@ def make_uhf_obj(fobj_a, fobj_b, frozen=False):
 
 
 def uccsd_restore_eris(symm, fobj_a, fobj_b, pad0=True, skip_Vab=False):
-    from pyscf import ao2mo
-
     """
     restore ERIs in the correct spin spaces
     """
@@ -98,53 +96,50 @@ def restore_eri_gen(targetsym, eri, norb1, norb2):
 
 
 def _convert_eri_gen(origsym, targetsym, eri, norb1, norb2):
-    import ctypes
-
-    from pyscf import lib
-
-    libao2mo = lib.load_library("libao2mo")
     """
-    #NOTE: IF YOU GET AN ERROR ABOUT THIS ATTRIBUTE:
-    This requires a custom PySCF compilation
-
-    Add to your PySCF and recompile: /pyscf/lib/ao2mo/restore_eri.c
-
-void AO2MOrestore_nr4to1_gen(double *eri4, double *eri1, int norb1, int norb2)
-{
-        size_t npair1 = norb1*(norb1+1)/2;
-        size_t npair2 = norb2*(norb2+1)/2;
-        size_t i, j, ij;
-        size_t d2 = norb2 * norb2;
-        size_t d3 = norb1 * d2;
-
-        for (ij = 0, i = 0; i < norb1; i++) {
-        for (j = 0; j <= i; j++, ij++) {
-                NPdunpack_tril(norb2, eri4+ij*npair2, eri1+i*d3+j*d2, HERMITIAN);
-                if (i > j) {
-                        memcpy(eri1+j*d3+i*d2, eri1+i*d3+j*d2,
-                               sizeof(double)*d2);
-                }
-        } }
-}
-
-void AO2MOrestore_nr1to4_gen(double *eri1, double *eri4, int norb1, int norb2)
-{
-        size_t npair1 = norb1*(norb1+1)/2;
-        size_t npair2 = norb2*(norb2+1)/2;
-        size_t i, j, k, l, ij, kl;
-        size_t d1 = norb2;
-        size_t d2 = norb2 * norb2;
-        size_t d3 = norb1 * d2;
-
-        for (ij = 0, i = 0; i < norb1; i++) {
-        for (j = 0; j <= i; j++, ij++) {
-                for (kl = 0, k = 0; k < norb2; k++) {
-                for (l = 0; l <= k; l++, kl++) {
-                        eri4[ij*npair2+kl] = eri1[i*d3+j*d2+k*d1+l];
-                } }
-        } }
-}
+        #NOTE: IF YOU GET AN ERROR ABOUT THIS ATTRIBUTE:
+        This requires a custom PySCF compilation
+
+        Add to your PySCF and recompile: /pyscf/lib/ao2mo/restore_eri.c
+
+    void AO2MOrestore_nr4to1_gen(double *eri4, double *eri1, int norb1, int norb2)
+    {
+            size_t npair1 = norb1*(norb1+1)/2;
+            size_t npair2 = norb2*(norb2+1)/2;
+            size_t i, j, ij;
+            size_t d2 = norb2 * norb2;
+            size_t d3 = norb1 * d2;
+
+            for (ij = 0, i = 0; i < norb1; i++) {
+            for (j = 0; j <= i; j++, ij++) {
+                    NPdunpack_tril(norb2, eri4+ij*npair2, eri1+i*d3+j*d2, HERMITIAN);
+                    if (i > j) {
+                            memcpy(eri1+j*d3+i*d2, eri1+i*d3+j*d2,
+                                   sizeof(double)*d2);
+                    }
+            } }
+    }
+
+    void AO2MOrestore_nr1to4_gen(double *eri1, double *eri4, int norb1, int norb2)
+    {
+            size_t npair1 = norb1*(norb1+1)/2;
+            size_t npair2 = norb2*(norb2+1)/2;
+            size_t i, j, k, l, ij, kl;
+            size_t d1 = norb2;
+            size_t d2 = norb2 * norb2;
+            size_t d3 = norb1 * d2;
+
+            for (ij = 0, i = 0; i < norb1; i++) {
+            for (j = 0; j <= i; j++, ij++) {
+                    for (kl = 0, k = 0; k < norb2; k++) {
+                    for (l = 0; l <= k; l++, kl++) {
+                            eri4[ij*npair2+kl] = eri1[i*d3+j*d2+k*d1+l];
+                    } }
+            } }
+    }
     """
+    libao2mo = lib.load_library("libao2mo")
+
     fn = getattr(libao2mo, "AO2MOrestore_nr%sto%s_gen" % (origsym, targetsym))
 
     if targetsym == 1:
diff --git a/src/molbe/fragment.py b/src/molbe/fragment.py
index 1eee2c07..1cd74dd2 100644
--- a/src/molbe/fragment.py
+++ b/src/molbe/fragment.py
@@ -2,8 +2,8 @@
 
 import sys
 
-from .autofrag import autogen
-from .helper import get_core
+from molbe.autofrag import autogen
+from molbe.helper import get_core
 
 
 class fragpart:
@@ -133,7 +133,8 @@ def __init__(
             print("exiting", flush=True)
             sys.exit()
 
-    from .lchain import chain
+    # importing the function here turns it into a proper method
+    from molbe.lchain import chain  # noqa: PLC0415
 
     def hchain_simple(self):
         """Hard coded fragmentation feature"""
diff --git a/src/molbe/helper.py b/src/molbe/helper.py
index 30448a0a..cdb93a43 100644
--- a/src/molbe/helper.py
+++ b/src/molbe/helper.py
@@ -6,7 +6,7 @@
 
 import h5py
 import numpy
-from pyscf import ao2mo
+from pyscf import ao2mo, gto, lib, scf
 
 
 def get_veff(eri_, dm, S, TA, hf_veff):
@@ -36,10 +36,6 @@ def get_veff(eri_, dm, S, TA, hf_veff):
         Effective HF potential in the embedding basis.
     """
 
-    import functools
-
-    from pyscf import scf
-
     # Transform the density matrix
     ST = numpy.dot(S, TA)
     P_ = functools.reduce(numpy.dot, (ST.T, dm, ST))
@@ -96,8 +92,6 @@ def get_scfObj(
         The SCF object after running the Hartree-Fock calculation.
     """
     # from 40-customizing_hamiltonian.py in pyscf examples
-    from pyscf import gto, scf
-
     nao = h1.shape[0]
 
     # Initialize a dummy molecule with the required number of electrons
@@ -175,9 +169,6 @@ def get_eri(i_frag, Nao, symm=8, ignore_symm=False, eri_file="eri_file.h5"):
     numpy.ndarray
         Electron repulsion integrals, possibly restored with symmetry.
     """
-    import h5py
-    from pyscf import lib
-
     # Open the HDF5 file and read the ERI for the specified fragment
     r = h5py.File(eri_file, "r")
     eri__ = numpy.array(r.get(i_frag))
diff --git a/src/molbe/lo.py b/src/molbe/lo.py
index 1cb66945..197fc0e5 100644
--- a/src/molbe/lo.py
+++ b/src/molbe/lo.py
@@ -1,13 +1,17 @@
 # Author(s): Henry Tran, Oinam Meitei, Shaun Weatherly
 #
+import functools
 import sys
 
 import numpy
+from numpy.linalg import eigh
+from pyscf.gto import intor_cross
 
 from molbe.external.lo_helper import (
     get_aoind_by_atom,
     reorder_by_atom_,
 )
+from molbe.helper import ncore_
 
 
 def dot_gen(A, B, ovlp):
@@ -92,8 +96,6 @@ def get_xovlp(mol, basis="sto-3g"):
         S12 - Overlap of two basis sets
         S22 - Overlap in new basis set
     """
-    from pyscf.gto import intor_cross
-
     mol_alt = mol.copy()
     mol_alt.basis = basis
     mol_alt.build()
@@ -214,11 +216,11 @@ def get_pao_native(Ciao, S, mol, valence_basis):
 
 def get_loc(mol, C, method, pop_method=None, init_guess=None):
     if method.upper() == "ER":
-        from pyscf.lo import ER as Localizer
+        from pyscf.lo import ER as Localizer  # noqa: PLC0415
     elif method.upper() == "PM":
-        from pyscf.lo import PM as Localizer
+        from pyscf.lo import PM as Localizer  # noqa: PLC0415
     elif method.upper() == "FB" or method.upper() == "BOYS":
-        from pyscf.lo import Boys as Localizer
+        from pyscf.lo import Boys as Localizer  # noqa: PLC0415
     else:
         raise NotImplementedError("Localization scheme not understood")
 
@@ -262,12 +264,6 @@ def localize(
        basis in the IAO partitioning.
        This is an experimental feature.
     """
-    import functools
-
-    from numpy.linalg import eigh
-
-    from .helper import ncore_
-
     if lo_method == "lowdin":
         es_, vs_ = eigh(self.S)
         edx = es_ > 1.0e-15
diff --git a/src/molbe/mbe.py b/src/molbe/mbe.py
index e9d2c2e3..fbcb2242 100644
--- a/src/molbe/mbe.py
+++ b/src/molbe/mbe.py
@@ -5,10 +5,13 @@
 
 import h5py
 import numpy
+from pyscf import ao2mo
 
 import molbe.be_var as be_var
-
-from .pfrag import Frags
+from molbe.be_parallel import be_func_parallel
+from molbe.eri_onthefly import integral_direct_DF
+from molbe.pfrag import Frags
+from molbe.solver import be_func
 
 
 class storeBE:
@@ -300,11 +303,12 @@ def __init__(
         else:
             self.initialize(None, compute_hf, restart=True)
 
-    from molbe.external.optqn import get_be_error_jacobian
-
-    from ._opt import optimize
-    from .lo import localize
-    from .rdm import compute_energy_full, rdm1_fullbasis
+    # The following imports turn the imported functions into proper methods
+    #  cannot be moved to head of file.
+    from molbe._opt import optimize  # noqa: PLC0415
+    from molbe.external.optqn import get_be_error_jacobian  # noqa: PLC0415
+    from molbe.lo import localize  # noqa: PLC0415
+    from molbe.rdm import compute_energy_full, rdm1_fullbasis  # noqa: PLC0415
 
     def print_ini(self):
         """
@@ -339,9 +343,6 @@ def initialize(self, eri_, compute_hf, restart=False):
         restart : bool, optional
             Whether to restart from a previous calculation, by default False.
         """
-        import h5py
-        from pyscf import ao2mo
-
         if compute_hf:
             E_hf = 0.0
 
@@ -410,8 +411,6 @@ def initialize(self, eri_, compute_hf, restart=False):
                 if (
                     self.integral_direct_DF
                 ):  # Use density fitting to generate fragment ERIs on-the-fly
-                    from .eri_onthefly import integral_direct_DF
-
                     integral_direct_DF(
                         self.mf, self.Fobjs, file_eri, auxbasis=self.auxbasis
                     )
@@ -497,9 +496,6 @@ def oneshot(
         clean_eri : bool, optional
             Whether to clean up ERI files after calculation, by default False.
         """
-        from .be_parallel import be_func_parallel
-        from .solver import be_func
-
         self.scratch_dir = scratch_dir
         self.solver_kwargs = solver_kwargs
 
diff --git a/src/molbe/misc.py b/src/molbe/misc.py
index 69260ea6..48b14c3a 100644
--- a/src/molbe/misc.py
+++ b/src/molbe/misc.py
@@ -4,9 +4,13 @@
 import sys
 import time
 
+import h5py
 import numpy
-from pyscf import gto, scf
+from pyscf import ao2mo, df, gto, qmmm, scf
 from pyscf.lib import chkfile
+from pyscf.tools import fcidump
+
+from molbe.fragment import fragpart
 
 
 def libint2pyscf(
@@ -90,8 +94,6 @@ def libint2pyscf(
     mol.incore_anyway = True
     if use_df:
         mf = scf.UHF(mol).density_fit() if unrestricted else scf.RHF(mol).density_fit()
-        from pyscf import df
-
         mydf = df.DF(mol).build()
         mf.with_df = mydf
     else:
@@ -116,10 +118,6 @@ def be2fcidump(be_obj, fcidump_prefix, basis):
         'embedding' to get the integrals in the embedding basis
         'fragment_mo' to get the integrals in the fragment MO basis
     """
-    import h5py
-    from pyscf import ao2mo
-    from pyscf.tools import fcidump
-
     for fidx, frag in enumerate(be_obj.Fobjs):
         # Read in eri
         read = h5py.File(frag.eri_file, "r")
@@ -171,10 +169,6 @@ def ube2fcidump(be_obj, fcidump_prefix, basis):
         'embedding' to get the integrals in the embedding basis
         'fragment_mo' to get the integrals in the fragment MO basis
     """
-    import h5py
-    from pyscf import ao2mo
-    from pyscf.tools import fcidump
-
     for fidx, frag in enumerate(be_obj.Fobjs_a):
         # Read in eri
         read = h5py.File(frag.eri_file, "r")
@@ -324,9 +318,10 @@ def be2puffin(
         syntax: {'Atom_X': 'ECP_for_X'; 'Atom_Y': 'ECP_for_Y'}
         By default None
     """
-    from .fragment import fragpart
-    from .mbe import BE
-    from .ube import UBE
+    # The following imports have to happen here to avoid
+    # circular dependencies.
+    from molbe.mbe import BE  # noqa: PLC0415
+    from molbe.ube import UBE  # noqa: PLC0415
 
     # Check input validity
     assert os.path.exists(xyzfile), "Input xyz file does not exist"
@@ -373,8 +368,6 @@ def be2puffin(
                 use_df = False
             if hcore is None:
                 if pts_and_charges:
-                    from pyscf import qmmm
-
                     print(
                         "Using QM/MM Point Charges: Assuming QM structure in Angstrom "
                         "and MM Coordinates in Bohr !!!"
@@ -396,8 +389,6 @@ def be2puffin(
                 mf = scf.UHF(mol).set(max_cycle=200).newton()
         else:  # restricted
             if pts_and_charges:  # running QM/MM
-                from pyscf import qmmm
-
                 print(
                     "Using QM/MM Point Charges: Assuming QM structure in Angstrom and "
                     "MM Coordinates in Bohr !!!"
diff --git a/src/molbe/pfrag.py b/src/molbe/pfrag.py
index c7469fa9..408566ed 100644
--- a/src/molbe/pfrag.py
+++ b/src/molbe/pfrag.py
@@ -4,9 +4,10 @@
 
 import h5py
 import numpy
+import scipy.linalg
 
-from .helper import *
-from .solver import schmidt_decomposition
+from molbe.helper import get_eri, get_scfObj, get_veff
+from molbe.solver import schmidt_decomposition
 
 
 class Frags:
@@ -201,8 +202,6 @@ def get_nsocc(self, S, C, nocc, ncore=0):
         numpy.ndarray
             Projected density matrix.
         """
-        import scipy.linalg
-
         C_ = functools.reduce(numpy.dot, (self.TA.T, S, C[:, ncore : ncore + nocc]))
         P_ = numpy.dot(C_, C_.T)
         nsocc_ = numpy.trace(P_)
diff --git a/src/molbe/rdm.py b/src/molbe/rdm.py
index 57d6b451..ac1ed285 100644
--- a/src/molbe/rdm.py
+++ b/src/molbe/rdm.py
@@ -1,6 +1,7 @@
 # Author(s): Oinam Romesh Meitei
 
 import numpy
+from pyscf import ao2mo, scf
 
 
 def rdm1_fullbasis(
@@ -48,8 +49,6 @@ def rdm1_fullbasis(
         The two-particle RDM in the MO basis
         (if return_ao is False and return_RDM2 is True).
     """
-    from pyscf import ao2mo
-
     # Copy the molecular orbital coefficients
     C_mo = self.C.copy()
     nao, nmo = C_mo.shape
@@ -236,9 +235,6 @@ def compute_energy_full(
     approximate or true cumulants, and to return the reduced density matrices (RDMs).
     The energy components are printed as part of the function's output.
     """
-
-    from pyscf import ao2mo, scf
-
     # Compute the one-particle reduced density matrix (RDM1) and the cumulant (Kumul)
     #   in the full basis
     rdm1f, Kumul, rdm1_lo, rdm2_lo = self.rdm1_fullbasis(
diff --git a/src/molbe/solver.py b/src/molbe/solver.py
index b3015ae5..e8f1a9e9 100644
--- a/src/molbe/solver.py
+++ b/src/molbe/solver.py
@@ -5,6 +5,8 @@
 import sys
 
 import numpy
+from pyscf import ao2mo, cc, fci, mcscf, mp
+from pyscf.cc.ccsd_rdm import make_rdm2
 
 from molbe import be_var
 from molbe.external.ccsd_rdm import (
@@ -13,6 +15,9 @@
     make_rdm2_uccsd,
     make_rdm2_urlx,
 )
+from molbe.external.uccsd_eri import make_eris_incore
+from molbe.external.unrestricted_utils import make_uhf_obj
+from molbe.helper import get_frag_energy, get_frag_energy_u
 
 
 def be_func(
@@ -89,12 +94,6 @@ def be_func(
         Depending on the options, it returns the norm of the error vector, the energy,
         or a combination of these values.
     """
-    import os
-
-    from pyscf import ao2mo, fci
-
-    from .helper import get_frag_energy
-
     rdm_return = False
     if relax_density:
         rdm_return = True
@@ -140,8 +139,7 @@ def be_func(
             rdm1_tmp = mc.make_rdm1(civec, mc.norb, mc.nelec)
 
         elif solver == "HCI":
-            from pyscf import hci
-            # pilot pyscf.hci only in old versions
+            from pyscf import hci  # noqa: PLC0415    # optional module
 
             nao, nmo = fobj._mf.mo_coeff.shape
 
@@ -174,7 +172,7 @@ def be_func(
             rdm2s = rdm2aa + rdm2ab + rdm2ab.transpose(2, 3, 0, 1) + rdm2bb
 
         elif solver == "SHCI":
-            from pyscf.shciscf import shci
+            from pyscf.shciscf import shci  # noqa: PLC0415    # shci is optional
 
             if scratch_dir is None and be_var.CREATE_SCRATCH_DIR:
                 tmp = os.path.join(be_var.SCRATCH, str(os.getpid()), str(fobj.dname))
@@ -205,7 +203,7 @@ def be_func(
             rdm1_tmp, rdm2s = mch.fcisolver.make_rdm12(0, nmo, nelec)
 
         elif solver == "SCI":
-            from pyscf import ao2mo, cornell_shci, mcscf
+            from pyscf import cornell_shci  # noqa: PLC0415  # optional module
 
             nao, nmo = fobj._mf.mo_coeff.shape
             nelec = (fobj.nsocc, fobj.nsocc)
@@ -390,10 +388,6 @@ def be_func_u(
         Depending on the options, it returns the norm of the error vector, the energy,
         or a combination of these values.
     """
-    from molbe.external.unrestricted_utils import make_uhf_obj
-
-    from .helper import get_frag_energy_u
-
     rdm_return = False
     if relax_density:
         rdm_return = True
@@ -591,8 +585,6 @@ def solve_mp2(mf, frozen=None, mo_coeff=None, mo_occ=None, mo_energy=None):
     pyscf.mp.mp2.MP2
         The MP2 object after running the calculation.
     """
-    from pyscf import mp
-
     # Set default values for optional parameters
     if mo_coeff is None:
         mo_coeff = mf.mo_coeff
@@ -670,9 +662,6 @@ def solve_ccsd(
         - cc__ (pyscf.cc.ccsd.CCSD, optional): CCSD object
             (if rdm_return is True and rdm2_return is False).
     """
-    from pyscf import cc
-    from pyscf.cc.ccsd_rdm import make_rdm2
-
     # Set default values for optional parameters
     if mo_coeff is None:
         mo_coeff = mf.mo_coeff
@@ -793,7 +782,7 @@ def solve_block2(mf: object, nocc: int, frag_scratch: str = None, **solver_kwarg
 
 
     """
-    from pyscf import dmrgscf, mcscf
+    from pyscf import dmrgscf  # noqa: PLC0415   # optional module
 
     use_cumulant = solver_kwargs.pop("use_cumulant", True)
     norb = solver_kwargs.pop("norb", mf.mo_coeff.shape[1])
@@ -950,10 +939,6 @@ def solve_uccsd(
         - rdm2 (tuple, numpy.ndarray, optional): Two-particle density matrix
             (if rdm2_return is True and rdm_return is True).
     """
-    from pyscf import ao2mo, cc
-
-    from molbe.external.uccsd_eri import make_eris_incore
-
     C = mf.mo_coeff
     nao = [C[s].shape[0] for s in [0, 1]]
 
@@ -1073,9 +1058,6 @@ def schmidt_decomposition(
         returns TA (above), number of orbitals in the fragment space, and number of
         orbitals in bath space
     """
-
-    import functools
-
     # Threshold for eigenvalue significance
     thres = 1.0e-10
 
diff --git a/src/molbe/ube.py b/src/molbe/ube.py
index 1743c957..3341925e 100644
--- a/src/molbe/ube.py
+++ b/src/molbe/ube.py
@@ -14,12 +14,15 @@
 
 import os
 
+import h5py
 import numpy
+from pyscf import ao2mo
 
 import molbe.be_var as be_var
-
-from .mbe import BE
-from .pfrag import Frags
+from molbe.be_parallel import be_func_parallel_u
+from molbe.mbe import BE
+from molbe.pfrag import Frags
+from molbe.solver import be_func_u
 
 
 class UBE(BE):  # 🍠
@@ -177,9 +180,6 @@ def __init__(
         self.initialize(mf._eri, compute_hf)
 
     def initialize(self, eri_, compute_hf):
-        import h5py
-        from pyscf import ao2mo
-
         if compute_hf:
             E_hf = 0.0
         EH1 = 0.0
@@ -418,9 +418,6 @@ def initialize(self, eri_, compute_hf):
     def oneshot(
         self, solver="UCCSD", nproc=1, ompnum=4, calc_frag_energy=False, clean_eri=False
     ):
-        from .be_parallel import be_func_parallel_u
-        from .solver import be_func_u
-
         if nproc == 1:
             E, E_comp = be_func_u(
                 None,
diff --git a/tests/chem_dm_kBE_test.py b/tests/chem_dm_kBE_test.py
index 892678eb..b26dab7e 100644
--- a/tests/chem_dm_kBE_test.py
+++ b/tests/chem_dm_kBE_test.py
@@ -13,13 +13,13 @@
 
 from kbe import BE, fragpart
 
+try:
+    import libdmet
+except ImportError:
+    libdmet = None
 
-class Test_kBE_Full(unittest.TestCase):
-    try:
-        import libdmet
-    except ImportError:
-        libdmet = None
 
+class Test_kBE_Full(unittest.TestCase):
     @unittest.skipIf(libdmet is None, "Module `libdmet` not imported correctly.")
     def test_kc2_sto3g_be1_chempot(self):
         kpt = [1, 1, 1]
diff --git a/tests/dmrg_molBE_test.py b/tests/dmrg_molBE_test.py
index cb151a9e..038a344b 100644
--- a/tests/dmrg_molBE_test.py
+++ b/tests/dmrg_molBE_test.py
@@ -11,13 +11,13 @@
 
 from molbe import BE, fragpart
 
+try:
+    from pyscf import dmrgscf
+except ImportError:
+    dmrgscf = None
 
-class TestBE_DMRG(unittest.TestCase):
-    try:
-        from pyscf import dmrgscf
-    except ImportError:
-        dmrgscf = None
 
+class TestBE_DMRG(unittest.TestCase):
     @unittest.skipIf(
         not os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true",
         reason="This test is known to fail.",