Enable whole testsuite

.github/workflows/build_docs.yml

@@ -26,6 +26,7 @@ jobs:
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
+        pip install .
         pip install --upgrade --upgrade-strategy eager -r docs/requirements.txt
         echo ${{ github.workspace }} > $(python -c "from distutils.sysconfig import get_python_lib; print(get_python_lib())")/quemb.pth
     - name: Test building docs & Archive the results

.github/workflows/quemb_unittest.yml

@@ -37,12 +37,17 @@ jobs:
         python -m pip install --upgrade pip
         pip install --upgrade --upgrade-strategy eager pytest
         if [ -f requirements.txt ]; then pip install --upgrade --upgrade-strategy eager -r requirements.txt; fi
+        pip install git+https://github.com/pyscf/dmrgscf
+        PYSCFHOME=$(pip show pyscf-dmrgscf | grep 'Location' | tr ' ' '\n' | tail -n 1)
+        wget https://raw.githubusercontent.com/pyscf/dmrgscf/master/pyscf/dmrgscf/settings.py.example
+        mv settings.py.example ${PYSCFHOME}/pyscf/dmrgscf/settings.py
+        pip install .
         echo ${{ github.workspace }} > $(python -c "from distutils.sysconfig import get_python_lib; print(get_python_lib())")/quemb.pth
 
     - name: Test with pytest
       run: |
         cd tests
-        pytest --doctest-modules --junitxml=junit/quemb-test-results_${{ matrix.python-version }}.xml
+        QUEMB_SKIP_EXPENSIVE_TESTS=true pytest --doctest-modules --junitxml=junit/quemb-test-results_${{ matrix.python-version }}.xml
       if: always()
 
     - name: Upload pytest junit results

.ruff.toml

@@ -0,0 +1,14 @@
+[lint]
+# see the rules [here](https://docs.astral.sh/ruff/rules/)
+select = ["E", "F", "I", "NPY", "ARG"]
+# select = ["E501"]
+ignore = [
+    "S101",
+        # https://docs.astral.sh/ruff/rules/assert/
+        # We want to use assert especially because it is disabled for
+        # optimized builds.
+    "E741",
+        # https://docs.astral.sh/ruff/rules/ambiguous-variable-name/
+        # Prevents the use of the characters 'l', 'O', or 'I' as variable names.
+        # Overly restrictive, in particular when implementing mathematical expressions.
+]

docs/requirements.txt

@@ -1,5 +1,2 @@
+# the additional (!) requirements for building the doc
 sphinx_rtd_theme
-numpy
-scipy
-h5py
-pyscf

docs/source/conf.py

@@ -8,24 +8,24 @@
 
 import os
 import sys
-sys.path.insert(0, os.path.abspath('../..'))
 
-project = 'QuEmb'
-copyright = '2024, Oinam Romesh Meitei'
-author = 'Oinam Romesh Meitei'
+sys.path.insert(0, os.path.abspath("../.."))
+
+project = "QuEmb"
+copyright = "2024, Oinam Romesh Meitei"
+author = "Oinam Romesh Meitei"
 
 # -- General configuration ---------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration
 
-extensions = ['sphinx.ext.autodoc', 'sphinx.ext.napoleon','sphinx_rtd_theme']
+extensions = ["sphinx.ext.autodoc", "sphinx.ext.napoleon", "sphinx_rtd_theme"]
 napoleon_google_docstring = False
 napoleon_include_init_with_doc = True
 napoleon_numpy_docstring = True
 exclude_patterns = []
 
 
-
 # -- Options for HTML output -------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output
 
-html_theme = 'sphinx_rtd_theme'
+html_theme = "sphinx_rtd_theme"

example/kbe_polyacetylene.py

@@ -3,24 +3,25 @@
 # pyscf's periodic HF calculation
 
 import numpy
-from pyscf.pbc import gto, scf, df
-from kbe import fragpart, BE
+from pyscf.pbc import df, gto, scf
+
+from kbe import BE, fragpart
 
 kpt = [1, 1, 3]
 cell = gto.Cell()
 
-a = 8.
-b = 8.
-c = 2.455 *2.
+a = 8.0
+b = 8.0
+c = 2.455 * 2.0
 
 lat = numpy.eye(3)
-lat[0,0] = a
-lat[1,1] = b
-lat[2,2] = c
+lat[0, 0] = a
+lat[1, 1] = b
+lat[2, 2] = c
 
 cell.a = lat
 
-cell.atom='''
+cell.atom = """
 H      1.4285621630072645    0.0    -0.586173422487319
 C      0.3415633681566205    0.0    -0.5879921146011252
 H     -1.4285621630072645    0.0     0.586173422487319
@@ -29,11 +30,11 @@
 C      0.3415633681566205    0.0     1.867007885398875
 H     -1.4285621630072645    0.0     3.041173422487319
 C     -0.3415633681566205    0.0     3.0429921146011254
-'''
+"""
 
-cell.unit='Angstrom'
-cell.basis = 'sto-3g'
-cell.verbose=0
+cell.unit = "Angstrom"
+cell.basis = "sto-3g"
+cell.verbose = 0
 cell.build()
 
 kpts = cell.make_kpts(kpt, wrap_around=True)
@@ -47,11 +48,9 @@
 kpoint_energy = kmf.kernel()
 
 # Define fragment in the supercell
-kfrag = fragpart(be_type='be2', mol=cell,
-                 kpt=kpt, frozen_core=True)
+kfrag = fragpart(be_type="be2", mol=cell, kpt=kpt, frozen_core=True)
 # Initialize BE
-mykbe = BE(kmf, fobj,
-           kpts=kpts)
+mykbe = BE(kmf, kfrag, kpts=kpts)
 
 # Perform BE density matching
-mykbe.optimize(solver='CCSD')
+mykbe.optimize(solver="CCSD")

example/molbe_dmrg_block2.py

@@ -2,9 +2,12 @@
 # `block2` is a DMRG and sparse tensor network library developed by the
 # Garnet-Chan group at Caltech: https://block2.readthedocs.io/en/latest/index.html
 
-import os, numpy, sys
+import os
+
 import matplotlib.pyplot as plt
-from pyscf import gto, scf, fci, cc
+import numpy
+from pyscf import cc, fci, gto, scf
+
 from molbe import BE, fragpart
 
 # We'll consider the dissociation curve for a 1D chain of 8 H-atoms:
@@ -18,8 +21,8 @@
 for a in seps:
     # Hartree-Fock serves as the starting point for all BE calculations:
     mol = gto.M()
-    mol.atom = [['H', (0.,0.,i*a)] for i in range(8)]
-    mol.basis = 'sto-3g'
+    mol.atom = [["H", (0.0, 0.0, i * a)] for i in range(8)]
+    mol.basis = "sto-3g"
     mol.charge = 0
     mol.spin = 0
     mol.build()
@@ -43,18 +46,18 @@
     # any clear advantage to using any one scheme over another,
     # the Pipek-Mezey scheme continues to be the most popular. With
     # BE-DMRG, localization takes place prior to fragmentation:
-    fobj = fragpart(be_type='be1', mol=mol)
+    fobj = fragpart(be_type="be1", mol=mol)
     mybe = BE(
         mf,
         fobj,
-        lo_method='pipek-mezey',                        # or 'lowdin', 'iao', 'boys'
-        pop_method='lowdin'                             # or 'meta-lowdin', 'mulliken', 'iao', 'becke'
-        )
+        lo_method="pipek-mezey",  # or 'lowdin', 'iao', 'boys'
+        pop_method="lowdin",  # or 'meta-lowdin', 'mulliken', 'iao', 'becke'
+    )
 
     # Next, run BE-DMRG with default parameters and maxM=100.
     mybe.oneshot(
         solver='block2',                                # or 'DMRG', 'DMRGSCF', 'DMRGCI'
-        scratch=scratch,                                # Scratch dir for fragment DMRG
+        scratch_dir=scratch,                            # Scratch dir for fragment DMRG
         maxM=100,                                       # Max fragment bond dimension
         force_cleanup=True,                             # Remove all fragment DMRG tmpfiles
         )
@@ -70,66 +73,61 @@
 # Finally, plot the resulting potential energy curves:
 fig, ax = plt.subplots()
 
-ax.plot(seps, fci_ecorr, 'o-', linewidth=1, label='FCI')
-ax.plot(seps, ccsd_ecorr, 'o-', linewidth=1, label='CCSD')
-ax.plot(seps, ccsdt_ecorr, 'o-', linewidth=1, label='CCSD(T)')
-ax.plot(seps, bedmrg_ecorr, 'o-', linewidth=1, label='BE1-DMRG')
+ax.plot(seps, fci_ecorr, "o-", linewidth=1, label="FCI")
+ax.plot(seps, ccsd_ecorr, "o-", linewidth=1, label="CCSD")
+ax.plot(seps, ccsdt_ecorr, "o-", linewidth=1, label="CCSD(T)")
+ax.plot(seps, bedmrg_ecorr, "o-", linewidth=1, label="BE1-DMRG")
 ax.legend()
 
-plt.savefig(os.path.join(scratch, f'BEDMRG_H8_PES{num_points}.png'))
+plt.savefig(os.path.join(scratch, f"BEDMRG_H8_PES{num_points}.png"))
 
 # (See ../quemb/example/figures/BEDMRG_H8_PES20.png for an example.)
 
 # For larger fragments, you'll want greater control over the fragment
 # DMRG calculations. Using the same setup as above for a single geometry:
 mol = gto.M()
-mol.atom = [['H', (0.,0.,i*1.2)] for i in range(8)]
-mol.basis = 'sto-3g'
+mol.atom = [["H", (0.0, 0.0, i * 1.2)] for i in range(8)]
+mol.basis = "sto-3g"
 mol.charge = 0
 mol.spin = 0
 mol.build()
-fobj = fragpart(be_type='be2', mol=mol)
-mybe = BE(
-    mf,
-    fobj,
-    lo_method='pipek-mezey',
-    pop_method='lowdin'
-    )
+fobj = fragpart(be_type="be2", mol=mol)
+mybe = BE(mf, fobj, lo_method="pipek-mezey", pop_method="lowdin")
 
 # We automatically construct the fragment DMRG schedules based on user keywords. The following
 # input, for example, yields a 60 sweep schedule which uses the two-dot algorithm from sweeps 0-49,
 # and the one-dot algo from 50-60. The fragment orbitals are also reordered according the Fiedler
 # vector procedure, along with a few other tweaks:
 
 mybe.optimize(
-    solver='block2',                                # or 'DMRG', 'DMRGSCF', 'DMRGCI'
-    scratch=scratch,                                # Scratch dir for fragment DMRG
-    startM=20,                                      # Initial fragment bond dimension (1st sweep)
-    maxM=200,                                       # Maximum fragment bond dimension
-    max_iter=60,                                    # Max number of sweeps
-    twodot_to_onedot=50,                            # Sweep num to switch from two- to one-dot algo.
-    max_mem=40,                                     # Max memory (in GB) allotted to fragment DMRG
-    max_noise=1e-3,                                 # Max MPS noise introduced per sweep
-    min_tol=1e-8,                                   # Tighest Davidson tolerance per sweep
-    block_extra_keyword=['fiedler'],                # Specify orbital reordering algorithm
-    force_cleanup=True,                             # Remove all fragment DMRG tmpfiles
+    solver="block2",  # or 'DMRG', 'DMRGSCF', 'DMRGCI'
+    scratch=scratch,  # Scratch dir for fragment DMRG
+    startM=20,  # Initial fragment bond dimension (1st sweep)
+    maxM=200,  # Maximum fragment bond dimension
+    max_iter=60,  # Max number of sweeps
+    twodot_to_onedot=50,  # Sweep num to switch from two- to one-dot algo.
+    max_mem=40,  # Max memory (in GB) allotted to fragment DMRG
+    max_noise=1e-3,  # Max MPS noise introduced per sweep
+    min_tol=1e-8,  # Tighest Davidson tolerance per sweep
+    block_extra_keyword=["fiedler"],  # Specify orbital reordering algorithm
+    force_cleanup=True,  # Remove all fragment DMRG tmpfiles
     only_chem=True,
 )
 
 # Or, alternatively, we can construct a full schedule by hand:
-schedule={
-    'scheduleSweeps': [0, 10, 20, 30, 40, 50],                  # Sweep indices
-    'scheduleMaxMs': [25, 50, 100, 200, 500, 500],              # Sweep maxMs
-    'scheduleTols': [1e-5,1e-5, 1e-6, 1e-6, 1e-8, 1e-8],        # Sweep Davidson tolerances
-    'scheduleNoises': [0.01, 0.01, 0.001, 0.001, 1e-4, 0.0],    # Sweep MPS noise
+schedule = {
+    "scheduleSweeps": [0, 10, 20, 30, 40, 50],  # Sweep indices
+    "scheduleMaxMs": [25, 50, 100, 200, 500, 500],  # Sweep maxMs
+    "scheduleTols": [1e-5, 1e-5, 1e-6, 1e-6, 1e-8, 1e-8],  # Sweep Davidson tolerances
+    "scheduleNoises": [0.01, 0.01, 0.001, 0.001, 1e-4, 0.0],  # Sweep MPS noise
 }
 
 # and pass it to the fragment solver through `schedule_kwargs`:
 mybe.optimize(
-    solver='block2',
+    solver="block2",
     scratch=scratch,
     schedule_kwargs=schedule,
-    block_extra_keyword=['fiedler'],
+    block_extra_keyword=["fiedler"],
     force_cleanup=True,
     only_chem=True,
 )
@@ -139,4 +137,4 @@
 # by `block2`.
 
 #NOTE: Parameters in `schedule_kwargs` will overwrite any other DMRG kwargs.
-#NOTE: The DMRG schedule kwargs and related syntax follows the standard notation used in block2.
+#NOTE: The DMRG schedule kwargs and related syntax follows the standard notation used in block2.

example/molbe_h8_chemical_potential.py

@@ -1,11 +1,13 @@
 # Illustrates a simple molecular BE calculation with chemical
 # potential matching
 
-from pyscf import gto,scf, fci
+from pyscf import fci, gto, scf
+
 from molbe import BE, fragpart
 
 # PySCF HF generated mol & mf (molecular desciption & HF object)
-mol = gto.M(atom='''
+mol = gto.M(
+    atom="""
 H 0. 0. 0.
 H 0. 0. 1.
 H 0. 0. 2.
@@ -14,7 +16,10 @@
 H 0. 0. 5.
 H 0. 0. 6.
 H 0. 0. 7.
-''',basis='sto-3g', charge=0)
+""",
+    basis="sto-3g",
+    charge=0,
+)
 
 mf = scf.RHF(mol)
 mf.conv_tol = 1e-12
@@ -23,39 +28,38 @@
 # Perform PySCF FCI to get reference energy
 mc = fci.FCI(mf)
 fci_ecorr = mc.kernel()[0] - mf.e_tot
-print(f'*** FCI Correlation Energy: {fci_ecorr:>14.8f} Ha', flush=True)
+print(f"*** FCI Correlation Energy: {fci_ecorr:>14.8f} Ha", flush=True)
 
 # Perform BE calculations with different fragment schemes:
 
 # Define BE1 fragments
-fobj = fragpart(be_type='be1', mol=mol)
+fobj = fragpart(be_type="be1", mol=mol)
 # Initialize BE
 mybe = BE(mf, fobj)
 # Perform chemical potential optimization
-mybe.optimize(solver='FCI', only_chem=True)
+mybe.optimize(solver="FCI", only_chem=True)
 
 # Compute BE error
 be_ecorr = mybe.ebe_tot - mybe.ebe_hf
-err_ = (fci_ecorr - be_ecorr)*100./fci_ecorr
-print(f'*** BE1 Correlation Energy Error (%) : {err_:>8.4f} %')
+err_ = (fci_ecorr - be_ecorr) * 100.0 / fci_ecorr
+print(f"*** BE1 Correlation Energy Error (%) : {err_:>8.4f} %")
 
 # Define BE2 fragments
-fobj = fragpart(be_type='be2', mol=mol)
+fobj = fragpart(be_type="be2", mol=mol)
 mybe = BE(mf, fobj)
-mybe.optimize(solver='FCI', only_chem=True)
+mybe.optimize(solver="FCI", only_chem=True)
 
 # Compute BE error
 be_ecorr = mybe.ebe_tot - mybe.ebe_hf
-err_ = (fci_ecorr - be_ecorr)*100./fci_ecorr
-print(f'*** BE2 Correlation Energy Error (%) : {err_:>8.4f} %')
+err_ = (fci_ecorr - be_ecorr) * 100.0 / fci_ecorr
+print(f"*** BE2 Correlation Energy Error (%) : {err_:>8.4f} %")
 
 # Define BE3 fragments
-fobj = fragpart(be_type='be3', mol=mol)
+fobj = fragpart(be_type="be3", mol=mol)
 mybe = BE(mf, fobj)
-mybe.optimize(solver='FCI', only_chem=True)
+mybe.optimize(solver="FCI", only_chem=True)
 
 # Compute BE error
 be_ecorr = mybe.ebe_tot - mybe.ebe_hf
-err_ = (fci_ecorr - be_ecorr)*100./fci_ecorr
-print(f'*** BE3 Correlation Energy Error (%) : {err_:>8.4f} %')
-
+err_ = (fci_ecorr - be_ecorr) * 100.0 / fci_ecorr
+print(f"*** BE3 Correlation Energy Error (%) : {err_:>8.4f} %")