Skip to content

Update docs for revised SGX algorithm #188

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
kylebystrom wants to merge 1 commit into
base: master
Choose a base branch
from
+45 −24
Open

Update docs for revised SGX algorithm #188

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 10 additions & 0 deletions source/user/ref.bib
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -840,6 +840,16 @@ @article{Neese2009
pages = {98},
year = {2009}
}
@article{Helmich2021,
title = {An Improved Chain of Spheres for Exchange Algorithm},
author = {{Helmich-Paris}, Benjamin and {de Souza}, Bernardo and Neese, Frank and Izs{\'a}k, R{\'o}bert},
year = 2021,
journal = {J. Chem. Phys.},
volume = {155},
number = {10},
pages = {104109},
doi = {10.1063/5.0058766},
}
@article{Cances2013,
author = {Cances, Eric and Maday, Yvon and Stamm, Benjamin },
title = {Domain decomposition for implicit solvation models},
Expand Down
59 changes: 35 additions & 24 deletions source/user/sgx.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,59 +5,70 @@ Seminumerical exchange (SGX)

*Modules*: :py:mod:`pyscf.sgx`

The SGX module implements seminumerical computation of the exchange matrix.
The SGX module implements seminumerical computation of the exchange matrix (and also the Coulomb matrix).

Introduction
============

Direct computation of the Hartree-Fock exchange matrix in the atomic orbital basis scales poorly with system size.
To achieve better scaling, one three-dimensional integral in the 6-dimensional two-electron integrals can be computed analytically, while the other can be evaluated on a real-space grid, as proposed by Friesner :cite:`Friesner1985`.
The PySCF implementation resembles the chain-of-spheres (COSX) algorithm of Neese et al. :cite:`Neese2009`, but uses more conservative grids and a slightly different P-junction screening function.
The PySCF implementation resembles the chain-of-spheres (COSX) algorithm of Neese et al. :cite:`Neese2009`, but uses slightly different default grids and a modified P-junction screening algorithm.
Overlap fitting is used to reduce aliasing errors :cite:`Izsak2011`.
SGX scales as :math:`O(N^2)` with system size, as opposed to the :math:`O(N^4)` scaling of analytical exchange :cite:`Neese2009`.
Before screening of negligible integrals, SGX scales as :math:`O(N^3)` with system size, as opposed to the :math:`O(N^4)` scaling of analytical exchange :cite:`Neese2009`. This makes its performance favorable for large basis sets and large systems.
In principle, both SGX and analytical exchange scale asymptotically as :math:`O(N^2)` if negligible integrals are screened out, and as :math:`O(N)` if the system has a gap and negligible density matrix components are screened.
However, the linear scaling is difficult to achieve in practice and only seen for some very large systems.

Usage and Example
=================

Any :attr:`scf.hf.SCF` object :attr:`mf` can be converted to an equivalent object that computes the Coulomb and exchange matrices with SGX instead of analytical integration by calling :code:`sgx.sgx_fit(mf)`.
Any :mod:`scf.hf.SCF` object :code:`mf` can be converted to an equivalent object that computes the Coulomb and exchange matrices with SGX instead of analytical integration by calling :code:`sgx.sgx_fit(mf)`.

* :source:`examples/sgx/00-simple_sgx.py`

.. literalinclude:: ../../examples/sgx/00-simple_sgx.py

In this case, the error of DFJ+SGX compared to analytical exchange is about 0.03 mEh. The line
In this case, the error of DF-J+SGX-K compared to analytical exchange is about 0.03 mEh. The line

.. code::

mf.with_df.dfj = True

specifies to compute the Coulomb matrix using density fitting (DF-J) while using SGX for the exchange matrix.

The :attr:`pjs` option turns on P-junction screening, which means that the three-center integrals are screened by the density matrix.
If a three-center integral gets contracted only with negligibly small density matrix elements, it is not computed.
This feature can only be used with :attr:`dfj=True`.
If :attr:`dfj=False` and :attr:`pjs=True`, :attr:`dfj` is set to :attr:`True`, and a warning is issued.
The P-junction screening threshold is determined by :attr:`mf.direct_scf_tol`.
Passing :code:`pjs=True` to :code:`sgx_fit` sets :code:`dfj=True` and also turns on P-junction screening (:code:`optk=True`), which means that the three-center integrals are screened by the density matrix.
If a three-center integral gets contracted only with negligibly small density matrix elements, it is not computed. Because the J-matrix cannot be screened in this manner, :code:`optk=True` only activates P-junction screening when used with :code:`dfj=True`.
If :code:`dfj=False`, :code:`optk` is ignored.

Direct evaluation of the J matrix can be used by setting :attr:`mf.with_df.direct_j = True`, but this is much slower than SGX-J or DF-J and only intended for testing purposes.
When :code:`optk=True`, the P-junction screening threshold is determined by two parameters: :code:`mf.with_df.sgx_tol_energy` and :code:`mf.with_df.sgx_tol_potential`.
These parameters correspond to approximate upper bounds on the error of the exchange energy and K-matrix due to P-junction screening, respectively.
By default, :code:`sgx_tol_energy=sgx_tol_potential="auto"`, in which case :code:`sgx_tol_energy` is set to :code:`mf.direct_scf_tol` and :code:`sgx_tol_potential=sgx_tol_energy**0.5`.
For more aggressive screening that is still numerically stable and sufficiently precise for most cases, we recommend :code:`sgx_tol_energy=mf.conv_tol` or :code:`0.1*mf.conv_tol` and :code:`sgx_tol_potential="auto"`.
The first thing to try if you have numerical issues is setting :code:`mf.rebuild_nsteps=1`.
This turns off incremental building of the K-matrix, which is faster but can accumulate errors between SCF steps.

Direct evaluation of the J matrix can be used by setting :code:`mf.with_df.direct_j = True`, but this is much slower than SGX-J or DF-J and only intended for testing purposes.

Adjustable Parameters
=====================

Calling the :attr:`sgx_fit` function on an :attr:`scf.hf.SCF` object returns an equivalent :attr:`SGXHF` object with a :attr:`with_df` attribute that handles SGX integration. To use a non-default auxiliary basis (for :attr:`dfj=True`), :attr:`auxbasis` can be specified in the call to :attr:`sgx_fit`. In addition, there are various adjustable parameters for SGX:
Calling the :code:`sgx_fit` function on an :mod:`scf.hf.SCF` object returns an equivalent :code:`SGXHF` object with a :code:`with_df` attribute that handles SGX integration. To use a non-default auxiliary basis (for :code:`dfj=True`), :code:`auxbasis` can be specified in the call to :code:`sgx_fit`. The other input option to :code:`sgx_fit` is :code:`pjs`. If :code:`pjs=False` (default), :code:`dfj=optk=False`, otherwise :code:`dfj=optk=True`. In addition, there are various adjustable parameters for SGX that can be set after initialization.

:code:`with_df` attributes:

:attr:`with_df` attributes:
* :code:`grids_level_i`: The grid level to use for initial SCF iterations.
* :code:`grids_level_f`: The grid level to use for final SCF iterations.
* :code:`dfj`: If :code:`True`, density fitting is used for the J-matrix. If :code:`False`, SGX is used.
* :code:`optk`: If :code:`True`, P-junction screening is used. Threshold determined by :code:`sgx_tol_energy` and :code:`sgx_tol_potential`. Only applies if :code:`dfj=True`.
* :code:`sgx_tol_energy` and :code:`sgx_tol_potential`: Sets the approximate upper bound on the energy and K-matrix errors due to density matrix screening. Ignored if :code:`optk=False` or :code:`dfj=False`. See the Usage and Example section above for instructions on how to set these, as well as the docstrings for more details.
* :code:`use_opt_grids`: Whether to use optimized grids for SGX based on COSX in ORCA :cite:`Helmich2021`. Default is True.
* :code:`fit_ovlp`: Whether to numerically fit the overlap matrix to improve numerical precision. Default is True.
* :code:`grids_thrd`: When the values of all atomic orbitals, multiplied by the qudrature weight, are less than this threshold for a given grid point, it is removed from the SGX numerical grid.
* :code:`grids_switch_thrd`: The threshold for the magnitude of the change in the density matrix that triggers the switch from the initial grid specified by :code:`grids_level_i` to the final grid specified by :code:`grids_level_f`.
* :code:`bound_aglo`: Determines how to estimate the maximum value of the three-center SGX integrals. Options are :code:`"ovlp"`, :code:`"sample"`, and :code:`"sample_pos"` (default). :code:`"ovlp"` assumes the maximum value of each integral is roughly equal to the maximum overlap between orbitals in a shell (with a value of 1 for shells on the same atom to account for orthogonality). :code:`"sample"` constructs a coarse grid for each atom pair and computes the maximum value of the integrals over these grids. :code:`"sample_pos"` does the same but also contructs an approximate estimate for the integral inside each grid block (position-dependent estimate), which can produce a small speedup for global exchange and a larger speedup for short-range exchange. Usually the default value :code:`"sample_pos"` should be fine.
* :code:`direct_j`: If :code:`True`, direct evaluation of the J-matrix is used (slow, for testing only).

* :attr:`grids_level_i`: The grid level to use for initial SCF iterations.
* :attr:`grids_level_f`: The grid level to use for final SCF iterations.
* :attr:`grids_thrd`: The grid points where no atomic orbital is significant (has a value greater than this threshold) are removed from consideration.
* :attr:`grids_switch_thrd`: The threshold for the magnitude of the change in the density matrix that triggers the switch from the initial grid specified by :attr:`grids_level_i` to the final grid specified by :attr:`grids_level_f`.
* :attr:`blockdim`: The maximum number of grid points to loop over at once. The number of grid points per batch is the minimum of :attr:`blockdim` and the maximum number of points allowed by the memory available for the calculation. The maximum memory can be adjusted by setting the :attr:`max_memory` attribute, which is initially set to :attr:`mol.max_memory`, the max memory setting for the Mole object.
* :attr:`dfj`: If :attr:`True`, density fitting is used for the J-matrix. If :attr:`False`, SGX is used.
* :attr:`direct_j`: If :attr:`True`, direct evaluation of the J-matrix is used (slow, for testing only).
* :attr:`pjs`: If :attr:`True`, P-junction screening is used. Threshold determined by :attr:`direct_scf_tol`.

:attr:`SGXHF` attribute:
:code:`SGXHF` attributes:

* :attr:`direct_scf_sgx`: Whether to use direct SCF within the SGX module, meaning that the J and K matrices are evaluated from the difference in the density matrix from the previous iteration.
* :attr:`rebuild_nsteps`: Rebuild the SGX JK matrix from scratch every :attr:`rebuild_nsteps` steps (default 5). Set to 0 to turn off rebuilding the JK matrix (Warning: This can cause creeping numerical error).
* :code:`auxbasis`: The auxiliary basis for J-matrix density fitting. Used if :code:`dfj=True`.
* :code:`rebuild_nsteps`: By default, when :code:`mf.direct_scf=True`, the K-matrix is constructed incrementally from the change in the density matrix between SCF steps (as is the case for other exchange algorithms in PySCF). The SGX J and K matrices are then rebuilt from scratch every :code:`rebuild_nsteps` steps (default 5). Set to 1 to rebuild from scratch at every step, and set to greater than the max number of cycles to turn off rebuilding the JK matrix (Warning: This can cause creeping numerical error). The SGX matrix is built from scratch when the numerical grid changes, regardless of the value of :code:`rebuild_nsteps`, because the same density matrix will have different energies for different grid sizes.