Patch v1.3.1 (#68)

* Adding box highlighting functionality to MatrixPlotter and argument to disable energies in ContourPlotter

* Refactoring TEI code and adding benchmark

* Expanding matrix visualization and updating manual

* Fixing bug in BlenderRender

Author: ifilot

docs/_static/img/overlap-co-boxes.png

@@ -6,43 +6,100 @@ Basis functions
 .. contents:: Table of Contents
     :depth: 3
 
-Gaussian type orbitals (GTO)
-============================
+Gaussian Type Orbitals (GTOs)
+=============================
 
-:program:`PyQInt` uses cartesian Gaussian type orbitals as given by
+:program:`PyQInt` employs *Cartesian Gaussian Type Orbitals* (GTOs), defined as
 
 .. math::
 
-    \Phi(\alpha,l,m,n,\vec{R}) = N (x - X)^{l} (y - Y)^{m} (z - Z)^{n} \exp \left(- \alpha |\vec{r} - \vec{R}|^{2} \right)
+    \Phi(\alpha, l, m, n, \vec{R}) =
+    N (x - X)^{l} (y - Y)^{m} (z - Z)^{n}
+    \exp \left(- \alpha \lvert \vec{r} - \vec{R} \rvert^{2} \right)
 
-wherein :math:`\alpha` is the exponent, :math:`\vec{R} = \left(X,Y,Z\right)` the
-position of the orbital, :math:`(l,m,n)` the orders of the pre-exponential
-polynomial, and :math:`N` a normalization constant such that
+Here, :math:`\alpha` denotes the orbital exponent,
+:math:`\vec{R} = (X, Y, Z)` the position of the orbital center,
+:math:`(l, m, n)` the orders of the Cartesian polynomial prefactor,
+and :math:`N` a normalization constant chosen such that
 
 .. math::
 
-    \left< \Phi | \Phi \right> = 1
+    \langle \Phi \mid \Phi \rangle = 1
 
-GTOs are a fundamental building block of CGF (see below) and typically a user
-would not directly work with them (a notable exception is provided below).
-Nevertheless, GTO objects can be constructed as follows::
+GTO Construction
+****************
 
-    from pyqint import PyQInt, CGF, GTO
+Gaussian Type Orbitals form the fundamental building blocks of *Contracted
+Gaussian Functions* (CGFs; see below). In typical workflows, users interact
+with CGFs rather than individual GTOs. Nevertheless, GTOs can be constructed
+explicitly when needed via::
 
-    coeff = 1.0    # coefficients only have meaning for GTOs within a CGF
-    alpha = 0.5
-    l,m,n = 0,0,0
-    p = (0,0,0)
-    G = GTO(coeff, p, alpha, l, m, n)
+    from pyqint import GTO
+
+    coeff = 1.0      # linear expansion coefficient
+    alpha = 0.5      # exponent
+    l, m, n = 0, 0, 0
+    position = (0.0, 0.0, 0.0)
+
+    gto = GTO(coeff, position, alpha, l, m, n)
 
 .. note::
-    * The normalization constant is automatically calculated by `PyQInt` based
-      on the value of :math:`\alpha` and :math:`(l,m,n)` and does not have
-      to be supplied by the user.
-    * If you work with individual GTOs, the first parameter to construct the GTO
-      should have a value of 1.0. This first parameter corresponds to the linear
-      expansion coefficient used in the formulation of Contracted Gaussian Functions
-      (see below).
+
+    - The normalization constant :math:`N` is computed automatically by
+      :program:`PyQInt` based on the values of :math:`\alpha` and
+      :math:`(l, m, n)` and must not be supplied explicitly.
+    - When working with individual GTOs (i.e. not as part of a CGF), the
+      coefficient should be set to ``1.0``. This parameter represents the
+      linear expansion coefficient used internally by CGFs.
+
+Retrieving GTO Data Members
+***************************
+
+The :code:`GTO` class is intentionally designed as a *flat and accessible*
+Python object. All defining parameters of the orbital are stored as public
+data members and can be accessed directly.
+
+The following attributes are available:
+
+- ``c`` — linear expansion coefficient
+- ``p`` — orbital center :math:`(X, Y, Z)`
+- ``alpha`` — Gaussian exponent
+- ``l``, ``m``, ``n`` — Cartesian polynomial orders
+
+For example::
+
+    coeff = gto.c
+    position = gto.p
+    alpha = gto.alpha
+    l = gto.l
+    m = gto.m
+    n = gto.n
+
+    print("Coefficient:", coeff)
+    print("Position:", position)
+    print("Exponent:", alpha)
+    print("Orders:", (l, m, n))
+
+These attributes fully define the mathematical form of the primitive Gaussian
+orbital and may be inspected or reused when constructing higher-level objects,
+such as CGFs.
+
+Evaluating a GTO
+****************
+
+In addition to direct data access, a GTO provides methods for numerical
+evaluation. For example, the orbital amplitude at a point
+:math:`(x, y, z)` can be obtained using::
+
+    value = gto.get_amp(x, y, z)
+
+The normalization constant used internally may be retrieved via::
+
+    norm = gto.get_norm()
+
+This separation between *data members* and *evaluation routines* allows users
+to both inspect the orbital parameters and efficiently compute values when
+needed.
 
 Contracted Gaussian Functions (CGF)
 ===================================

docs/basis_functions.rst

@@ -108,4 +108,66 @@ In this example:
 
 .. figure:: _static/img/coefficient-co.png
 
-    Coefficient matrix of CO, visualized using the :code:`MatrixPlotter` class.
+    Coefficient matrix of CO, visualized using the :code:`MatrixPlotter` class.
+
+Highlighting elements
+---------------------
+
+In some cases, it is useful not only to visualize a matrix, but also to
+emphasize specific elements or regions within it. This can be achieved using
+the ``boxes`` argument.
+
+The ``boxes`` argument expects a list of 5-tuples. Each tuple defines a
+rectangular region and consists of the following elements, in order:
+
+- the starting x-coordinate
+- the starting y-coordinate
+- the width of the rectangle
+- the height of the rectangle
+- the color of the rectangle
+
+The example below demonstrates how to draw red, green, and blue rectangles on
+top of a matrix visualization.
+
+.. code-block:: python
+
+    from pyqint import MoleculeBuilder, HF, MatrixPlotter
+
+    def main():
+        mol = MoleculeBuilder().from_name('CO')
+        res = HF(mol, 'sto3g').rhf(verbose=False)
+
+        labels = MatrixPlotter.generate_ao_labels([
+            ('O', ('1s', '2s', '2p')),
+            ('C', ('1s', '2s', '2p')),
+        ])
+
+        # highlighting boxes
+        boxes = [
+            (0,0,1,1,'#BB0000'),
+            (1,2,3,1,'#00BB00'),
+            (3,4,1,3,'#0000BB'),
+        ]
+
+        MatrixPlotter.plot_matrix(
+            mat=res['overlap'],
+            filename='overlap-co.png',
+            xlabels=labels,
+            ylabels=labels,
+            xlabelrot=0,
+            title='Overlap',
+            boxes=boxes
+        )
+
+    if __name__ == '__main__':
+        main()
+
+.. figure:: _static/img/overlap-co-boxes.png
+
+    Coefficient matrix of CO, visualized using the :code:`MatrixPlotter` class,
+    where a number of elements have been highlighted.
+
+.. note::
+
+    The ``(x, y)`` coordinates use **zero-based indexing**. Consequently, the
+    top-left corner of the first box is specified as ``(0, 0)``.

docs/matrix_visualization.rst

@@ -5,7 +5,7 @@
 # for this basis set
 integrator = PyQInt()
 st = time.perf_counter()
-mol = MoleculeBuilder.from_name('cubane')
+mol = MoleculeBuilder.from_name('benzene')
 cgfs, nuclei = mol.build_basis('sto3g')
 S, T, V, tetensor = integrator.build_integrals_openmp(cgfs, nuclei)
 end = time.perf_counter()

examples/benzene_tei.py

@@ -0,0 +1,12 @@
+from pyqint import MoleculeBuilder, PyQInt
+import time
+
+# simple test to see how long it takes to build the two-electron integrals
+# for this basis set
+integrator = PyQInt()
+st = time.perf_counter()
+mol = MoleculeBuilder.from_name('cubane')
+cgfs, nuclei = mol.build_basis('sto3g')
+S, T, V, tetensor = integrator.build_integrals_openmp(cgfs, nuclei)
+end = time.perf_counter()
+print('Elapsed time: %.2f s' % (end - st))

examples/cubane.py

@@ -0,0 +1,30 @@
+from pyqint import MoleculeBuilder, HF, MatrixPlotter
+
+def main():
+    mol = MoleculeBuilder().from_name('CO')
+    res = HF(mol, 'sto3g').rhf(verbose=False)
+
+    labels = MatrixPlotter.generate_ao_labels([
+        ('O', ('1s', '2s', '2p')),
+        ('C', ('1s', '2s', '2p')),
+    ])
+
+    # highlighting boxes
+    boxes = [
+        (0,0,1,1,'#BB0000'),
+        (1,2,3,1,'#00BB00'),
+        (3,4,1,3,'#0000BB'),
+    ]
+
+    MatrixPlotter.plot_matrix(
+        mat=res['overlap'],
+        filename='overlap-co.png',
+        xlabels=labels,
+        ylabels=labels,
+        xlabelrot=0,
+        title='Overlap',
+        boxes=boxes
+    )
+
+if __name__ == '__main__':
+    main()

examples/matrixplotter.py

@@ -43,7 +43,7 @@ module-path = "src"
 # ---------------------------------------------------------------------------
 [project]
 name = "pyqint"
-version = "1.3.0"
+version = "1.3.1"
 description = "Python package for evaluating integrals of Gaussian type orbitals"
 readme = "README.md"
 license = { text = "GPL-3.0-only" }

pyproject.toml

@@ -0,0 +1,20 @@
+from pyqint import PyQInt, MoleculeBuilder
+import timeit
+import statistics
+
+mol = MoleculeBuilder.from_name('benzene')
+cgfs, nuclei = mol.build_basis('sto3g')
+integrator = PyQInt()
+
+def tei():
+    tetensor = integrator.build_tei_openmp(cgfs)
+
+runs = timeit.repeat(
+    tei,
+    repeat=5,     # number of independent runs
+    number=1      # iterations per run
+)
+
+print(f"Mean: {statistics.mean(runs):.6f}s")
+print(f"Std dev: {statistics.stdev(runs):.6f}s")
+print(f"Min: {min(runs):.6f}s")

scripts/bench.py

@@ -237,11 +237,13 @@ def __build_isosurface(self, filename, scalarfield, isovalue, sz=5):
         isovalue = np.abs(isovalue)
         unitcell = np.diag(np.ones(3) * 2 * sz)
 
-        # try to import PyTessel but do not throw an error if it cannot be loaded
         try:
             from pytessel import PyTessel
-        except ModuleNotFoundError:
-            print('Cannot find module PyTessel')
+        except ModuleNotFoundError as exc:
+            raise ImportError(
+                "Isosurface rendering requires the optional dependency 'pytessel'. "
+                "Please install it (e.g. `pip install pytessel`)."
+            ) from exc
 
         pytessel = PyTessel()
         vertices, normals, indices = pytessel.marching_cubes(scalarfield.flatten(), scalarfield.shape, unitcell.flatten(), isovalue)

src/pyqint/blenderrender.py

@@ -24,6 +24,7 @@ def build_contourplot(
         dpi: int = 144,
         ngrid: int = 5,
         labels=None,
+        plot_energies = True,
     ):
         """
         Generate a grid of contour plots for molecular orbitals.
@@ -41,6 +42,7 @@ def build_contourplot(
         plane : str or tuple
             Plane specification:
             - str: one of {'xy', 'xz', 'yz'}
+            - list with str
             - tuple: (atom_i, atom_j, atom_k, up_direction)
         sz : float
             Half-width of the plotted plane.
@@ -80,6 +82,10 @@ def build_contourplot(
                     grid = ContourPlotter.__create_cartesian_grid(sz, npts, plane)
                     xlabel = '$%s$ [a.u.]' % plane[0]
                     ylabel = '$%s$ [a.u.]' % plane[1]
+                elif isinstance(plane, list):
+                    grid = ContourPlotter.__create_cartesian_grid(sz, npts, plane[orb_idx])
+                    xlabel = '$%s$ [a.u.]' % plane[orb_idx][0]
+                    ylabel = '$%s$ [a.u.]' % plane[orb_idx][1]
                 else:
                     # Arbitrary plane defined by three atoms
                     grid = ContourPlotter.__build_plane_from_atoms(
@@ -140,14 +146,15 @@ def build_contourplot(
                 if ylabel is not None:
                     ax[i,j].set_ylabel(ylabel)
 
-                # Title handling
                 if labels is not None:
-                    ax[i, j].set_title(
-                        r"%s (%.4f Ht)"
-                        % (labels[orb_idx], res["orbe"][orb_idx])
-                    )
+                    orblabel = labels[orb_idx]
                 else:
-                    ax[i, j].set_title(r"$\psi_{%i}$ (%.4f Ht)" % (orb_idx + 1, res["orbe"][orb_idx]))
+                    orblabel = r'$\psi_{%i}$' % (orb_idx + 1)
+
+                if plot_energies:
+                    orblabel += r' (%.4f Ht)' % res["orbe"][orb_idx]
+
+                ax[i, j].set_title(orblabel)
 
         plt.tight_layout()
         plt.savefig(filename)