rough draft of some docs to make sure things are working

ciderpress/dft/settings.py

@@ -226,39 +226,47 @@ def __init__(self, slist, nk0, nk1, l1_dots, nd1=0, ld_dots=None, ndd=0):
         on the indexing of the features, fl_rho is the part of the density
         ingredient vector AFTER the semilocal part.
 
+        Note on the l=1 feature indexing. If l1_dots[i] = (j, k), then::
+
+            fl_feat[i+nk0] = einsum(
+                'xg,xg->g',
+                fl_rho[3*j+nk0 : 3*j+3+nk0],
+                fl_rho[3*k+nk0 : 3*k+3+nk0],
+            )
+
+        If ld_dots[i] = (j, k), then::
+
+            nstart = nk0 + 3 * nk1
+            fl_feat[i + nk0 + len(l1_dots)] = einsum(
+                'xg,xg->g',
+                fl_rho[3*j+nstart : 3*j+3+nstart],
+                fl_rho[3*k+nstart : 3*k+3+nstart],
+            )
+            nstart = nk0 + 3 * nk1 + 3 * nd1
+            fl_feat[i + nk0 + len(l1_dots) + len(ld_dots)] = fl_rho[i + nstart]
+
         Args:
             slist (list of float): s parameters for the fractional
-                Laplacian. For each s, (-Laplacian)**s phi_i is computed
+                Laplacian. For each s, :math:`(-\\Delta)^s \\phi_i` is computed
                 for each single-particle orbital.
             nk0 (int): Number of scalar (l=0) features. Must be <= len(slist).
-                fl_feat[i] = fl_rho[i] = (-Laplacian')**slist[i] rho(r, r') |_r'=r
+                fl_feat[i] = fl_rho[i] = :math:`(-\\Delta')^s \\rho(r, r') |_{r'=r}`,
+                with s=slist[i].
             nk1 (int): Number of vector (l=1) features F_s^1. Must be <= len(slist).
-                fl_rho[3*i+nk0 : 3*i+3+nk0] = nabla'' (-Laplacian')**slist[i]
-                                              rho(r'', r') |_r'=r,r''=r=
+                fl_rho[3*i+nk0 : 3*i+3+nk0] =
+                :math:`\\nabla'' (-\\Delta')^s \\rho(r'', r') |_{r'=r,r''=r}`,
+                with s=slist[i].
             l1_dots (list of (int, int)): List of index tuples for contracting
                 the l=1 F_s^1 features. -1 indicates to use the semilocal
-                density gradient. If l1_dots[i] = (j, k), then
-                    fl_feat[i+nk0] = einsum(
-                        'xg,xg->g',
-                        fl_rho[3*j+nk0 : 3*j+3+nk0],
-                        fl_rho[3*k+nk0 : 3*k+3+nk0],
-                    )
-            nd1 (int): Number of vector features F_s^d. Must be <= len(slist)
-                nstart = nk0 + 3 * nk1
-                fl_rho[i * 3 + nstart : (i+1) * 3 + nstart] =
-                    nabla' (-Laplacian')**slist[i] rho(r, r') |_r'=r
+                density gradient.
+            nd1 (int): Number of vector features F_s^d. Must be ``<= len(slist)``
+                With ``nstart = nk0 + 3 * nk1``,
+                ``fl_rho[i * 3 + nstart : (i+1) * 3 + nstart] =``
+                :math:`\\nabla' (-\\Delta')^s \\rho(r, r') |_{r'=r}`,
+                with s=slist[i].
             ld_dots (list of (int, int)): Same as l1_dots but for the F_s^d
                 features. -1 indicates to use the semilocal density gradient.
-                If ld_dots[i] = (j, k), then
-                    nstart = nk0 + 3 * nk1
-                    fl_feat[i + nk0 + len(l1_dots)] = einsum(
-                        'xg,xg->g',
-                        fl_rho[3*j+nstart : 3*j+3+nstart],
-                        fl_rho[3*k+nstart : 3*k+3+nstart],
-                    )
-            ndd (int): Numer of dot product features F_s^{dd}. Must be <= nd1
-                nstart = nk0 + 3 * nk1 + 3 * nd1
-                fl_feat[i + nk0 + len(l1_dots) + len(ld_dots)] = fl_rho[i + nstart]
+            ndd (int): Numer of dot product features :math:`F_s^{dd}`. Must be <= nd1
         """
         self.slist = slist
         assert nk0 <= self.npow
@@ -449,7 +457,7 @@ class SADMSettings(SDMXBaseSettings):
     the self-repulsion of the exchange hole constructed only from
     the spherically averaged density matrix at a point.
 
-    WARNING DEPRECATED, especially 'exact' which is not numerically
+    WARNING: DEPRECATED, especially 'exact' which is not numerically
     accurate. Use SDMXSettings and its variants instead.
     """
 
@@ -546,7 +554,7 @@ def __init__(self, pows, ndt):
         Initialize SDMXGSettings
 
         Args:
-            pows (list of int): list or 0, 1, 2, see SDMXSettings docstring.
+            pows (list of int): list of 0, 1, 2, see SDMXSettings docstring.
             ndt (int): Number of gradient features H_n^d. Will compute features
                 for the first ndt values of n in pows.
         """
@@ -603,7 +611,7 @@ def __init__(self, pows, n1):
         Initialize SDMX1Settings
 
         Args:
-            pows (list of int): list or 0, 1, 2, see SDMXSettings docstring.
+            pows (list of int): list of 0, 1, 2, see SDMXSettings docstring.
             n1 (int): Number of gradient features H_n^1. Will compute features
                 for the first n1 values of n in pows.
         """
@@ -648,7 +656,7 @@ def __init__(self, pows, nd, n1):
         Initialize SDMXG1Settings
 
         Args:
-            pows (list of int): list or 0, 1, 2, see SDMXSettings docstring.
+            pows (list of int): list of 0, 1, 2, see SDMXSettings docstring.
             nd (int): Number of gradient features H_n^d. Will compute features
                 for the first nd values of n in pows.
             n1 (int): Number of gradient features H_n^1. Will compute features
@@ -893,50 +901,61 @@ def get_reasonable_normalizer(self):
         raise NotImplementedError
 
 
+ALLOWED_I_SPECS_L0 = ["se", "se_r2", "se_apr2", "se_ap", "se_ap2r2", "se_lapl"]
 """
+Allowed spec strings for version i l=0 features.
+
 se: squared-exponential
+
 se_ap: squared-exponential times the exponent of the
-       r' coordinate.
+r' coordinate.
+
 se_apr2: squared-exponential times the exponent of the
-         r' coordinate times r^2
+r' coordinate times r^2
+
 se_ap2r2: squared-exponential times exponent^2 times r^2
+
 se_lapl: laplacian of squared-exponential
 """
-ALLOWED_I_SPECS_L0 = ["se", "se_r2", "se_apr2", "se_ap", "se_ap2r2", "se_lapl"]
-
 
+ALLOWED_I_SPECS_L1 = ["se_grad", "se_rvec"]
 """
+Allowed spec strings for version i l=1 features
+
 se_grad: gradient of squared-exponential
+
 se_rvec: squared-exponential times vector (r'-r)
 """
-ALLOWED_I_SPECS_L1 = ["se_grad", "se_rvec"]
 
+ALLOWED_J_SPECS = ["se", "se_ar2", "se_a2r4", "se_erf_rinv"]
 """
+Allowed spec strings for version j features.
+Version k features have the same allowed spec strings
+
 se: squared-exponential
+
 se_ar2: squared-exponential * a * r^2
+
 se_a2r4: squared-exponential * a^2 * r^4
+
 se_erf_rinv: squared-exponential * 1/r with short-range
-             erf damping
+erf damping
 """
-ALLOWED_J_SPECS = ["se", "se_ar2", "se_a2r4", "se_erf_rinv"]
 ALLOWED_K_SPECS = ALLOWED_J_SPECS
-
 """
-TODO docs here
+See ``ALLOWED_J_SPECS``
 """
-ALLOWED_RHO_MULTS = ["one", "taumix", "dampmix", "expnt"]
 
+ALLOWED_RHO_MULTS = ["one", "taumix", "dampmix", "expnt"]
 """
 TODO docs here
 """
-ALLOWED_RHO_DAMPS = ["none", "exponential", "asymptotic_const"]
 
+ALLOWED_RHO_DAMPS = ["none", "exponential", "asymptotic_const"]
 """
-Uniform-scaling powers (USPs) descibe how features scale as the
-density is scaled by n_lambda(r) = lambda^3 n(lambda r). If the
-USP of a functional F is u, then
-F[n_lambda](r) = lambda^u F[n](lambda r)
+TODO docs here
 """
+
 SPEC_USPS = {
     "se": 0,
     "se_r2": -2,
@@ -951,6 +970,12 @@ def get_reasonable_normalizer(self):
     "se_rvec": -1,
     "grad_rho": 4,
 }
+"""
+Uniform-scaling powers (USPs) descibe how features scale as the
+density is scaled by n_lambda(r) = lambda^3 n(lambda r). If the
+USP of a functional F is u, then
+F[n_lambda](r) = lambda^u F[n](lambda r)
+"""
 RHO_MULT_USPS = {
     "one": 0,
     "taumix": 0,

docs/ciderpress/dft/dft.rst

@@ -3,3 +3,5 @@ DFT Module
 
 .. toctree::
    :maxdepth: 2
+
+   settings

docs/ciderpress/dft/settings.rst

@@ -0,0 +1,7 @@
+.. _settings_module:
+
+The Settings Module
+===================
+
+.. automodule:: ciderpress.dft.settings
+   :members:

docs/conf.py

@@ -22,13 +22,16 @@
 extensions = [
     'sphinx.ext.autodoc',
     'sphinx_rtd_theme',
+    'sphinx.ext.mathjax',
     'sphinx.ext.napoleon',
     'breathe',
 ]
 
 templates_path = ['_templates']
 exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store']
 
+autoclass_content = 'both'
+
 breathe_projects = {"CiderPress": "xml"}
 breathe_default_project = "CiderPress"
 breathe_default_members = ("members", "undoc-members")

docs/features/features.rst

@@ -0,0 +1,40 @@
+Density and Orbital Features in CiderPress
+==========================================
+
+To predict the exchange-correlation energy :math:`E_{\text{xc}}`, we need
+to train a machine learning model :math:`e_{\text{xc}}(\mathbf{x})` such that
+
+.. math:: E_{\text{xc}} = \int \text{d}^3\mathbf{r}\, e_\text{xc}(\mathbf{x}[n_1](\mathbf{r}))
+
+where :math:`\mathbf{x}[n_1](\mathbf{r})` is a position-dependent feature vector that is
+a functional of the density :math:`n(\mathbf{r})` in "pure" Kohn-Sham DFT and a 
+functional of the density matrix :math:`n_1(\mathbf{r}, \mathbf{r}')` in the case of
+"generalized" Kohn-Sham DFT. Cider provides several types of feature that can
+be used as input to the ML model. These features
+can be divided into four categories:
+
+* :ref:`Semilocal Features (SL) <sl_feat>`: Same features as in conventional GGA/meta-GGA functionals (i.e. :math:`n`, :math:`\nabla n`, :math:`\tau`).
+  NOTE: All Cider models must include semilocal features.
+  They need not be used explicitly in the model, but evaluating
+  them is required to evalute baseline functionals and other quantities in the code.
+* :ref:`Nonlocal Density Features (NLDF) <nldf_feat>`: These features are constructed by integrating the density
+  over a real-space kernel function to characterize the shape of the density around a point :math:`\mathbf{r}`.
+* :ref:`Nonlocal Orbital Features (NLOF) <nlof_feat>`: EXPERIMENTAL, NOT TESTED, NOT RECOMMENDED FOR USE.
+  One coordinate of the density matrix is operated on by a fractional Laplacian.
+* :ref:`Smooth Density Matrix Exchange (SDMX) <sdmx_feat>`:
+  One coordinate of the density matrix is convolved at different length scales. Then these convolutions
+  are contracted to approximately characterize the shape of the density matrix around a point :math:`\mathbf{r}`.
+
+The set of features to be used in a model is specified using the :class:`FeatureSettings` object. To see
+the code API for setting up feature settings, see the :ref:`Settings module <settings_module>` section. To see
+mathematical descriptions and physical intuition for the different types of features, see
+the subsections below.
+
+.. toctree::
+   :maxdepth: 1
+
+   sl
+   nldf
+   nlof
+   sdmx
+

docs/features/nldf.rst

@@ -0,0 +1,7 @@
+.. _nldf_feat:
+
+Nonlocal Density Features (NLDF)
+================================
+
+TODO
+

docs/features/nlof.rst

@@ -0,0 +1,10 @@
+.. _nlof_feat:
+
+Nonlocal Orbital Features (NLOF)
+================================
+
+**The nonlocal orbital features are still highly experimental. We note here
+that they are implemented in the code so that their presence does not cause
+confusion, but their use is not encouraged, and they might be removed at a
+later date.**
+

docs/features/sdmx.rst

@@ -0,0 +1,7 @@
+.. _sdmx_feat:
+
+Smooth Density Matrix Exchange (SDMX)
+=====================================
+
+TODO
+

docs/features/sl.rst

@@ -0,0 +1,66 @@
+.. _sl_feat:
+
+Semilocal Features (SL)
+=======================
+
+Semilocal features are the ingredient commonly used in semilocal DFT.
+They most basic ingredient is the electron density :math:`n(\mathbf{r})`,
+which is defined as
+
+.. math:: n(\mathbf{r}) = \sum_i f_i |\phi_i(\mathbf{r})|^2
+
+with :math:`\phi_i(\mathbf{r})` being the Kohn-Sham orbitals. A functional
+constructed from :math:`n(\mathbf{r})` only is called a local density
+approximation (LDA). The gradient of the density :math:`\nabla n(\mathbf{r})`
+is also commonly used and results in a generalized-gradient approximation (GGA).
+Lastly, the kinetic energy density
+
+.. math:: \tau(\mathbf{r}) = \frac{1}{2} \sum_i f_i |\nabla\phi_i(\mathbf{r})|^2
+
+can be used, resulting in a meta-generalized gradient approximation (meta-GGA or MGGA).
+To help enforce physical constraints such as uniform scaling, regularized
+features are often introduced, including the reduced gradient
+
+.. math:: p = \frac{|\nabla n|^2}{2(3\pi^2)^{1/3}n^{4/3}}
+
+and the iso-orbital indicator
+
+.. math:: \alpha = \frac{\tau - \tau_W}{\tau_0}
+
+where :math:`\tau_W=|\nabla n|^2/8n` is the single-electron kinetic energy,
+and :math:`\tau_0=\frac{3}{10}(3\pi^2)^{2/3}n^{5/3}` is the kinetic energy
+density of the uniform electron gas.
+
+Note that for simplicity, all of these definitions are provided for the
+non-spin-polarized case, where the orbital occupations :math:`f_i\in [0,2]`
+
+In Cider, the choice of semilocal features is specified using the class:`SemilocalSettings`
+class in :py:mod:`ciderpress.dft.settings`. There is only one argument to this
+class, :py:obj:`mode`, which specifies which features to compute. There are four choices:
+
+* ``ns``: :math:`\mathbf{x}_\text{sl} = [n, |\nabla n|^2]`
+* ``np``: :math:`\mathbf{x}_\text{sl} = [n, p]`
+* ``nst``: :math:`\mathbf{x}_\text{sl} = [n, |\nabla n|^2, \tau]`
+* ``npa``: :math:`\mathbf{x}_\text{sl} = [n, p, \alpha]`
+
+**IMPORTANT NOTE 1**: If the mode is ``ns`` or ``np``, the kinetic energy density :math:`\tau`
+is not computed, so it cannot be used at any point in the calculation. This means
+that :ref:`Nonlocal Density Features <nldf_feat>` must be instantiated with
+``sl_level="GGA"``. If desired, orbital-dependent features can still be incorporated
+via :ref:`Smooth Density Matrix Exchange <sdmx_feat>`. If the mode is ``nst`` or ``npa``,
+:math:`\tau` is always computed.
+
+**IMPORTANT NOTE 2**: The regularized features :math:`p` and :math:`\alpha` are
+scale-invariant, meaning that
+an exchange functional trained with these features and
+a proper exchange functional baseline will obey the
+uniform scaling rule (see :ref:`Uniform Scaling Constraints <unif_scaling>`).
+The raw features :math:`n`, :math:`|\nabla n|^2`, and :math:`\tau` are not
+scale-invariant, and CiderPress does not automatically regularize these features
+as it does with nonlocal features. Therefore, these features must be regularized in
+the :py:mod:`ciderpress.dft.transform_data` module to enforce the uniform scaling
+constraint in trained models. (TODO need to elaborate on this).
+
+See the :class:`SemilocalSettings` class in the :ref:`Feature Settings <settings_module>`
+documentation for more details on the API for setting up these features.
+

docs/index.rst

@@ -14,10 +14,12 @@ kylebystrom@gmail.com.
 
 
 .. toctree::
-   :maxdepth: 4
+   :maxdepth: 2
    :caption: Contents:
 
    installation/installation
+   theory/theory
+   features/features
    ciderpress/dft/dft
    ciderpress/pyscf/pyscf
    ciderpress/gpaw/gpaw

docs/theory/theory.rst

@@ -0,0 +1,12 @@
+Theory Overview
+===============
+
+The theory section provides an overview of some of the important concepts
+involved in designing exchange-correlation functionals for DFT,
+particularly machine learning-based functionals.
+
+.. toctree::
+   :maxdepth: 1
+
+   uniform_scaling
+