Quickstart
==========


Coulomb matrix computation
--------------------------

This example demonstrates the complete workflow: building a molecule with PySCF, converting bases to `libintx` structures, constructing the DF metric `(P|Q)`, and computing the Coulomb matrix `J`.

.. code-block:: python

   import numpy as np
   from pyscf import gto, scf, df
   import libintx
   from libintx.tools import pyscf_to_libintx_basis

   # Build a PySCF molecule
   mol = gto.Mole()
   mol.atom = '''
   O  0.000000  0.000000  0.000000
   H  0.000000 -0.757160  0.586260
   H  0.000000  0.757160  0.586260
   '''
   mol.basis = 'def2-svp'
   # Spherical harmonic GTOs are not yet supported:
   mol.cart = True
   mol.build()

   # Build an auxiliary basis and its two-center metric (P|Q)
   auxmol = df.addons.make_auxmol(mol)  # set up the default auxiliary basis
   V = auxmol.intor('int2c2e', hermi=1)  # 2-center 2-electron integrals matrix

   # Convert PySCF bases to libintx Basis
   basis = pyscf_to_libintx_basis(mol)
   aux_basis = pyscf_to_libintx_basis(auxmol)

   # Create the J-engine and compute J
   jeng = libintx.Jengine(basis, aux_basis, V, vinv_type='cholesky', vinv_device='cuda')

   D = np.eye(mol.nao)

   J = jeng.compute_j(D, scale=0.5)

   # Check the results
   mf = scf.RHF(mol).density_fit()
   J_ref = mf.get_j(dm=D)

   print(f'|J_ref - J| = {np.linalg.norm(J_ref-J)}')


Schwarz Screening
-----------------

Schwarz screening can dramatically reduce compute time, especially for large systems. You can enable screening using per-shell two-center two-electron :math:`(P|P)` and four-center two-electron :math:`(ij|ij)` estimates.

**Computation:**

* :math:`(P|P)` can be computed from the diagonal AO blocks of `(P|Q)` aggregated per auxiliary shell.
  It should be a vector of length `nsh_aux`, where `nsh_aux` is the number of shells in the auxiliary basis (`len(aux_basis)`).
* :math:`(ij|ij)` can be initialized via PySCF's RHF helper.
  It should be a matrix of shape `(nsh, nsh)`, where `nsh` is the number of shells in the primary basis (`len(basis)`).

.. note::
   Some PySCF bases are general-contracted; convert them to segmented-contraction variants before computing screening quantities.

An example code:

.. code-block:: python

    import numpy as np
    from pyscf import scf
    import libintx

    # Compute screening quantities
    pp = np.sqrt(np.array([np.max(np.abs(V[i:j,i:j])) for i,j in aux_basis.aos_by_shell()]))
    ijij = scf.RHF(mol).init_direct_scf().q_cond

    # Enable screening
    jeng.set_screening(pp, ijij, cutoff=1e-12)

Solver Options
--------------

To use CUDA-based solver, pass `vinv_device='cuda'`. Optionally adjust `vinv_mem_fraction` (default `0.5`) to limit the maximum fraction of GPU memory for the solver.

.. code-block:: python

    jeng = libintx.Jengine(
        basis, aux_basis, V,
        vinv_device='cuda',
        vinv_mem_fraction=0.5
    )

If memory limits are exceeded, the Cholesky/LU factors will be stored on the host and loaded to the GPU when needed.
Consider distributed Cholesky solver (`cholesky_mp`) for very large molecules (requires MPI).