Add IntegrationGrid settings

[EBB2675]

The settings of numerical integration grids for molecular systems should be added to somewhere in NumericalSettings.

The angular grid scheme, the radial grid scheme, atom partitioning, and pruning method are the most important quantities that I can think of right now.

[ndaelman-hu]

angular grid scheme, the radial grid scheme

These should already be covered. We already have a Mesh class for describing the geometry.
I think that we also have a section to denote integration precision. At least in the old schema, maybe this wasn't ported.

[JosePizarro3]

At least in the old schema, maybe this wasn't ported.

Can you point exactly to what you meant here?

The angular grid scheme, the radial grid scheme

@EBB2675 what do you mean here? Are these grids in some spherical coordinates system?

[EBB2675]

Indeed. A product between a radial and an angular grid is built for each atom, then transformed into a molecular grid.

I was thinking of (roughly) something like this:

class MolecularGrid(Mesh):
    """Settings for integration grids.
    TODO: connect to AtomsState
    """

    type = Quantity(
        type=str,
        description= """ type of the grid, e.g. exchange-correlation or COSX""",
    ) 

    angular_scheme = Quantity(
        type=str,  
        description=""" 
        the angular quadrature scheme.
        e.g. Lebedev: A. D. Becke, J. Chem. Phys. 88, 2547 (1988)
        """,
    )

    angular_pruning_method = Quantity(
        type=str,
        description=""" 
        Angular grid pruning method, usually adaptive.
        """,
    )

    radial_scheme = Quantity(
        type=str,
        description=""" 
        the radial quadrature scheme, e.g., Gauss-Chebyschev 
        """,
    )

    atom_partitioning = Quantity(
        type=str,
        description=""" 
        Weight generation scheme, typically C. W. Murray, N. C. Handy and G. J. Laming, Mol. Phys. 78, 997 (1993) 
        """,

    )

[JosePizarro3]

Ok, so I think it is time for restructuring better Mesh, and work out on what we are actually doing here.

Most importantly, let's do some Math. These schema is about handling the following N-dimensional integrations:

$$\int_{\vec{r}_a}^{\vec{r}_b} d^3 \vec{r} F(\vec{r})$$

where this integration is typically numerically obtained by discrete summations (https://en.wikipedia.org/wiki/Numerical_integration, https://en.wikipedia.org/wiki/Gaussian_quadrature):

$$\approx \sum_{n=a}^{b} w(\vec{r}_n) F(\vec{r}_n)$$

Note:

1. We need to define the space for our integration, i.e., our $\vec{r}_n$ which is equivalent to define Mesh.
2. Function $F$ changes depending on the problem at hand, but a well-placed string could be enough to characterize the type of function which is being integrated and which rule (e.g., which Gaussian quadrature rule) is being applied.
3. Given that we tend to treat also orbital wavefunctions (so that $F(\vec{r}) \rightarrow R(r) Y(\theta, \phi)$), this means that we might have two distinct integration limits, i.e., two different meshes.

So, I propose the following schema (let me know what you think):

class Mesh(NumericalSettings):
    """
    A base section used to define the mesh or space partitioning over which a discrete numerical integration is performed.
    """

    dimensionality = Quantity(
        type=np.int32,
        default=3,
        description="""
        Dimensionality of the mesh: 1, 2, or 3. Defaults to 3.
        """,
    )

    kind = Quantity(
        type=MEnum('equidistant', 'logarithmic', 'tan'),
        shape=['dimensionality'],
        description="""
        Kind of mesh identifying the spacing in each of the dimensions specified by `dimensionality`. It can take the values:

        | Name      | Description                      |
        | --------- | -------------------------------- |
        | `'equidistant'`  | Equidistant grid (also known as 'Newton-Cotes') |
        | `'logarithmic'`  | log distance grid |
        | `'Tan'`  | Non-uniform tan mesh for grids. More dense at low abs values of the points, while less dense for higher values |
        """,
    )

    grid = Quantity(
        type=np.int32,
        shape=['dimensionality'],
        description="""
        Amount of mesh point sampling along each axis.
        """,
    ) 

    n_points = Quantity(
        type=np.int32,
        description="""
        Number of points in the mesh.
        """,
    )

    points = Quantity(
        type=np.complex128,
        shape=['n_points', 'dimensionality'],
        description="""
        List of all the points in the mesh.
        """,
    )

    multiplicities = Quantity(
        type=np.float64,
        shape=['n_points'],
        description="""
        The amount of times the same point reappears. A value larger than 1, typically indicates
        a symmetry operation that was applied to the `Mesh`.
        """,
    )

    pruning = Quantity(
        type=str,
        description="""
        Pruning method applied for reducing the amount of points in the Mesh. This is typically
        used for numerical integration near the core levels in atoms, and it takes the value
        `adaptive`.
        """
    )

    def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
        super().normalize(archive, logger)

        if self.dimensionality not in [1, 2, 3]:
            logger.error('`dimensionality` meshes different than 1, 2, or 3 are not supported.')


class NumericalIntegration(NumericalSettings):
    """
    Numerical integration settings used to resolve the following type of integrals by discrete
    numerical integration:
    
    ```math
    \int_{\vec{r}_a}^{\vec{r}_b} d^3 \vec{r} F(\vec{r}) \approx \sum_{n=a}^{b} w(\vec{r}_n) F(\vec{r}_n)   
    ```

    Here, $F$ can be any type of function which would define the type of rules that can be applied
    to solve such integral (e.g., 1D Gaussian quadrature rule or multi-dimensional `angular` rules like the 
    Lebedev quadrature rule).
    
    These multidimensional integral has a `Mesh` defined over which the integration is performed, i.e., the 
    $\vec{r}_n$ points. 
    """

    coordinate = Quantity(
        type=MEnum('all', 'radial', 'angular'),
        description="""
        Coordinate over which the integration is performed. `all` means the integration is performed in 
        all the space. `radial` and `angular` describe cases where the integration is performed for
        functions which can be splitted into radial and angular distributions (e.g., orbital wavefunctions). 
        """,
    )

    integration_rule = Quantity(
        type=str,  # ? extend to MEnum?
        description="""
        Integration rule used. This can be any 1D Gaussian quadrature rule or multi-dimensional `angular` rules, 
        e.g., Lebedev quadrature rule (see e.g., Becke, Chem. Phys. 88, 2547 (1988)).
        """
    )

    weight_partitioning = Quantity(
        type=str,
        description="""
        Approximation applied to the weight when doing the numerical integration. See e.g., C. W. Murray, N. C. Handy 
        and G. J. Laming, Mol. Phys. 78, 997 (1993).
        """
    )

    mesh = SubSection(sub_section=Mesh.m_def)

    def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
        super().normalize(archive, logger)

If you want radial and angular integrations, you need to instantiate twice the NumericalIntegration class. An alternative design would be to have multiple classes, each for each numerical integration.

Now, I was reading a bit the ORCA docu page (https://www.faccts.de/docs/orca/6.0/manual/contents/detailed/numint.html), and I understand you wanted to add xc and cosx. However, I cannot understand whether this is very ORCA-specific, or something which could be shown in these 2 base classes.

But! This shouldn't prevent us to inherit from Mesh and redefine stuff if you need it. E.g., if the formula in the ORCA docu page describes another kind of Mesh, you can inherit from it and define it so that it is understandable.

[EBB2675]

thanks a lot! @JosePizarro3

• Since the meshes can be atom dependent, e.g. better grids for heavier atoms, would it be a good idea to add a reference for AtomsState inside NumericalIntegration?

• In addition, the accuracy of integration is controlled by two further cutoff values, basis function cutoff and integration weight cutoff, which I think should be added to the NumericalIntegration.

• Pure DFT exchange has one-type of grid, and HF exchange is another grid “COSX” (often smaller) for hybrid functionals. This applies to RIJCOSX approximation which is now default for hybrids. Turbomole also supports COSX and but it is called SENEX (seminumeric exchange) in there. So it would be great to add a quantity to signal for what kind of grid is NumericalIntegration for.

• Following the previous point, in RIJCOSX approximation for instance, the SCF iterations are done with a small grid, and gradients and final energies are evaluated on a larger, more accurate grid.

So at the end, a typical calculation can look like:

DFT GRID
…settings….

COSX GRID X1
…settings….

COSX GRID X2
…settings….

COSX GRID X3
…settings….

and thats why I wanted to distinguish between XC and COSX grids. What would be the best way to accommodate these four grids in the schema?

[JosePizarro3]

• Since the meshes can be atom dependent, e.g. better grids for heavier atoms, would it be a good idea to add a reference for AtomsState inside NumericalIntegration?

Question: does the atom-dependent integration sums up to give a global property? If so, I would define something called AtomicNumericalIntegration(NumericalIntegration) and add atoms_state_ref = Quantity(type=AtomsState) there. Otherwise, yes, go ahead and add atoms_state_ref directly there.

• In addition, the accuracy of integration is controlled by two further cutoff values, basis function cutoff and integration weight cutoff, which I think should be added to the NumericalIntegration.

Yeah, I was writing a long answer but I see the problem now. I think we could keep the definition of the function to be integrated, $F(\vec{r})$ as part of another numerical settings section called BasisSet. In that case, the cutoff for that BasisSet is defined there, while for the integration weights could be defined here.

I wonder now whether we maintain this structure, we should add basis_set_ref = Quantity(type=BasisSet) in the NumericalIntegration class.

• Pure DFT exchange has one-type of grid, and HF exchange is another grid “COSX” (often smaller) for hybrid functionals. This applies to RIJCOSX approximation which is now default for hybrids. Turbomole also supports COSX and but it is called SENEX (seminumeric exchange) in there. So it would be great to add a quantity to signal for what kind of grid is NumericalIntegration for.
• Following the previous point, in RIJCOSX approximation for instance, the SCF iterations are done with a small grid, and gradients and final energies are evaluated on a larger, more accurate grid.

Should this be added to Mesh.kind? What are these grids, i.e., how do they partition the space?

[ndaelman-hu]

Starting from @JosePizarro3's proposal (#137 (comment)). I largely agree with the suggestion, but have several remarks. I'll just list them here, but can make an example schema if needed.

1. Mesh.dimensionality: agreed. Note that the positions (which are contingent on this) are set to 3D in other places. We should apply a dimensionality generalization in a separate MR.
2. Mesh.kinds: should we add numerical parameters like "spacing"?
3. Mesh.points and Mesh.multiplicities: the interdependence has to be clearer. When dealing with analysis methods, e.g. k_line_density, points is useless, as it is unclear up to which degree they were symmetry-reduced. I'd have multiplicities default to np.ones((n_points, 1)) (unless already set) via normalization. I also need an optional quantity to denote the symmetry reduction.
4. Mesh.pruning: if it takes the value "adaptive", how about using that as a default? Moreover, pls clarify what "adaptive" stands for exactly.
5. Mesh: change the name grid for grid_shape. Maybe also change kind (does it have special GUI visualization rights)?
6. NumericalIntegration: I would make this more flexible than what coordinate suggests. There are many kinds of meshes used for similar ends. The partition between radial and angular aside, there's e.g. rectangular grids (whose boundary adapts to the system's geometry). I would move all coordinate system (boundary) responsibility to Mesh (Geometry) and make it so more complex Meshes may be built up from simpler ones. Whether a 3D rectangular mesh should conventionally be deconstructed into 3 1D meshes, is an open question. I'm curious how NeXus solved this.
7. NumericalIntegration.coordinate: in case you forego the previous remark, change all to full and correct the phrasing in the description from "all" to "whole".
8. NumericalIntegration.weight_partitioning: this name may lead to confusion given Mesh.multiplicities.

[ndaelman-hu]

Since the meshes can be atom dependent, e.g. better grids for heavier atoms, would it be a good idea to add a reference for AtomsState inside NumericalIntegration?

I'd reverse the reference direction.

In addition, the accuracy of integration is controlled by two further cutoff values, basis function cutoff and integration weight cutoff, which I think should be added to the NumericalIntegration.

Could you elaborate on both? I suppose that the basis function cutoff that you have in mind is specific to atom-centered basis sets. If so, I agree with Chema's answer, that is should be implemented by the child class of BasisSet.