wavefunction91 / integratorxx Goto Github PK
View Code? Open in Web Editor NEWReusable DFT Grids for the Masses
License: BSD 3-Clause "New" or "Revised" License
Reusable DFT Grids for the Masses
License: BSD 3-Clause "New" or "Revised" License
According to the Catch2 documentation
We strongly recommend against using Approx when writing new code. You should be using floating point matchers instead.
The tests should accordingly use
WithinAbs(double target, double margin)
,WithinRel(FloatingPoint target, FloatingPoint eps)
andWithinULP(FloatingPoint target, uint64_t maxUlpDiff)
An improvement over #87 is to use minimal basis non-integer Slater-type orbital wave functions. Such wave functions have been reported by Koga et al in Int. J. Quantum Chem. 62, 1 (1997).
Unfortunately, this paper does not specify what the employed orbital occupations are. However, the occupations are likely the same as in other works reporting (differences to) numerical Hartree-Fock energies; occupations have been given in e.g. Atomic Data and Nuclear Data Tables 95 (2009) 836
I am not sure if these are already implemented, but...
Murray et al describe a scheme which is exact for all spherical harmonics up to angular momentum
Treutler and Ahlrichs describe a more compact scheme with Lobatto polynomials with
Baker, Andzelm, Scheiner, and Delley describe a radial grid in eq (12) of J. Chem. Phys. 101, 8894 (1994), which is not currently available.
Shizgal, Ho, and Yang propose a Gauss-Maxwell quadrature in J. Math. Chem. 55, 413 (2017).
It is not immediately obvious what this scheme actually is.
addendum: it looks like a much better reference for the actual nodes and weights is J. Comput. Phys. 41, 309 (1981)
For some weight functions, (e.g.
It requires the solution of a Hankel-type tridiagonal eigenvalue problem, which tends to be numerically unstable for large orders, but having it as an option would lower the barrier to entry for testing out new quadratures and an important validation tool moving forward.
Requires:
LAPACK++
to solve the EVPThe Gauss-Jacobi quadrature integrates functions of the form
where
The other quadratures in quadrature/
have consistent file names, but the Gauss-Chebyshev rules have shorter names. They should be renamed, e.g. gausscheby1.hpp
-> gausschebyshev1.hpp
.
Slater's rules offer a simple way to compute somewhat realistic atomic electron densities. They are also used in the formation of some quadrature grids. The rules used in Molpro have been summarized in J. Chem. Phys. 149, 044103 (2018); in addition, one needs to assign formal occupations to the orbitals.
The angular grids are precomputed and stored statically.
I wonder whether it would make sense to also precompute the primitive quadratures, since this could be done in arbitrary precision with e.g. sympy
.
Storing the weights and nodes for
Storage cost | |
---|---|
100 | 80 kB |
250 | 52 kB |
500 | 2 MB |
2500 | 50 MB |
It might be a good idea to to have this functionality anyway, since one should test the accuracy of the computed nodes and weights by comparison to accurate reference values.
#22 fixed the GaussChebyshev2 implementation, but it still needs a unit test.
Lindh, Malmqvist and Gagliardi describe a radial quadrature rule for Gaussian basis sets in Theor. Chem. Acc. 106, 178 (2001). This would be nice to have as well.
The Gauss-Legendre routine only computes half of the nodes and weights, since they are symmetric.
However, I think the Gauss-Chebyshev rules 1 and 2 (and its modification) also generate symmetric nodes and weights, allowing this symmetry to be taken advantage of.
This is a requirement for integration with legacy codes
At
, gausscheby1 has a comment // TODO: Tranform on bound
points[i] = pt;
weights[i] = wgt;
What do you mean by this?
Add the following weight partition schemes:
#46 includes tests for low orders of nodes.
I think it would be good to also include tests for higher numbers of nodes, since quadrature grids often use even several hundreds of points.
However, not every number of points needs to be tested. I think that going from 100, it might be reasonable to test rules every 20 points, (100, 120, 140, 160, 180, 200), then maybe every 25 points (225, 250, 275, 300, 325, 350, 375, 400), then every 50 points (450, 500, 550, 600).
Reducing the spacing of the tests is a big reduction storage space, and the cost to generate the tests; even the 100-point rules take some seconds to generate on my laptop with SymPy...
The Gauss-Chebyshev rules 1 and 2 are found in SciPy; one just has to modify the weights, which can be done in the generator.
Delley's paper on numerical atomic orbital calculations, J. Chem. Phys. 92, 508 (1990), states employing a partition function (eq. 3)
where "preferred choices are the functions"
where
Delley writes
These partition functions automatically depend on atom sizes via
$p_A$ and have given consistently good results for all sorts of compounds including heteronuclear compounds containing heavy transition metals and light elements including hydrogen. The last function leads to continuous vanishing of all partition functions and all their derivatives on all nuclei where they are not centered.
These schemes are clearly similar to Becke's, and should be quite easy to implement. The only thing necessary are just the atomic electron densities. The adaptive Molpro grid wavefunction91/GauXC#51 implementation of J. Chem. Phys. 157, 234106 (2022) employed Slater atomic densities:
In Ref. 79, these are defined for elements up to
$n = 6$ with effective principal numbers of$n^\star = 3.7$ , 4.0, 4.2 for$n = 4, 5, 6$ . The implementation in Molpro uses$n^\star = 4.4$ for seventh row elements based on a simple linear extrapolation of the$n^\star$ values given by Leach.79 The advantage of using Slater’s densities over tabulated numerical densities is that they are given analytically and can be easily computed at any given grid point.
In Int. J. Quantum Chem. 34, 55 (1988), Savin employs a partitioning where the atomic weight is given by
where
This partitioning is already available in CRYSTAL.
The source files don't mention any license at all.
The license should be included at the start of every source file.
NWChem implements an error function based partitioning described in J. Chem. Phys. 152, 184102 (2020) on page 8
where
and
Weber, Daul and Baltensperger discuss adaptive radial grids using sinc quadrature in Computer Physics Communications 163, 133 (2004). They also present equations for STOs for using Lindh's scheme from Theor Chem Acc 106, 178 (2001).
I don't know if this will be useful, but I am opening the ticket as a note.
Per discussion in #27 (comment), the "canonical" implementation of the MK quadrature uses a Uniform Trapezoid base quadrature. However Lehtola, Marques J. Chem. Phys. 157, 174114 (2022) points out that accuracy of this scheme can be improved via using Gauss-Chebyshev (second kind) as the base quadrature.
Ahrens and Beylkin have described rotationally invariant quadratures on the sphere in Proc. R. Soc. A (2009) 465, 3103–3125.
All quadratures seem to be on the same footing, separate into (1) priimtive quadratures (GC, GL, etc), (2) radial quadratures and (3) angular quadratures.
At present, the implementations of the radial rules hardcode the primitive quadrature rule, and the radial size adjustment. This leads to code duplication, and limits the reusability of the code.
The radial rules should be constructed from an array of nodes and weights. The size adjustment, which consists of scaling the radial nodes by R
and the weights by R^3
should likewise happen in a dedicated routine.
In particular, each of the primitive quadratures are guaranteed to produce exact integrals for polynomials of a particular order subject to a particular weight function. This can be tested directly by randomly generating polnomials for the unit tests.
For the radial quadratures, each ansatz was designed to integrate function of a specific form (e.g. radial integrands). These tests should be updated to test these classes of functions.
#35 pointed out that the bounds transformations from the "natural" bounds of a quadrature should be decoupled from the quadrature implementation to reduce code replication. This will require two steps:
Womersley has published t designs ranging up to 16 k points.
The weight used in gausscheby2 are incorrect.
There are two problems here. Compared to the equations given on Wikipedia
Self explanatory
The Gauss-Legendre quadrature routine hardcodes an absolute tolerance eps = 3.0e-11
. The origin of this value is not documented, and it also seems a bit too large to me. Since eps
should be of the order of machine epsilon of the employed data type.
Another thing that needs to be fixed is to employ a for
loop instead of the current while
loop which risks hanging. Instead, one should take e.g. a maximum of 100 Newton iterations, and then throw an error if convergence was not achieved.
The code should check that the point and weight arrays have the same size. This check should be configurable and possible to turn off for production level.
Hayle and Trefethen discuss quadrature formulas from conformal maps in SIAM J. Numer. Anal. 46, 930 (2008), which may prove useful.
#23 added formatting files, the corresponding GH action should be enabled
https://github.com/wavefunction91/MACIS/blob/master/.github/workflows/clang_format_lint.yaml
Mieke Peels' PhD thesis states that they have obtained new angular integration grids
The new integration grids have been generated to integrate polynomials, and therefore spherical harmonics, up through order 𝐿 = 105. Many of the new grids contain fewer points than the Lebedev grids. None of the grids contain negative weights, and the weight spreads are generally smaller than for the Lebedev grids.
The thesis also states that
In the future, grid generators containing optimized seed points will be released open source as part of the angular integration grid generator (aigg) program, which would allow for the new grids to be implemented in existing quantum chemistry software packages.
I have asked Gerald Knizia about these grids in October 2022, and in March 2023, but have not received a reply. I did contact Peels as well, who said I should be in touch with Knizia about the code.
My assumption is that Gerald is simply too busy, so if someone runs into him at a conference, please ask about this.
Add Becke's original radial transform
as a quadrature method
The tests are currently using Catch2. However, many Linux distributions have already switched to Catch3.
The migration guide is available at
https://github.com/catchorg/Catch2/blob/devel/docs/migrate-v2-to-v3.md
(Sorry, hit enter too soon.) Anyway, to be brief, I think you're mentioning IntegratorXX::integratorxx
but exporting IntegratorXX::IntegratorXX
.
https://github.com/wavefunction91/IntegratorXX/blob/master/cmake/IntegratorXXConfig.cmake.in#L48 is the lowercase. Otherwise, easy enough to add at psi4/psi4#3041
Gauss-Lobatto quadrature should be added as well
The double-exponential radial quadrature of Takahasi and Mori, also known as tanh-sinc quadrature, was suggested by Mitani in Theor. Chem. Acc. 130, 645 (2011) and Theor. Chem. Acc. 131, 1169 (2012) for density functional quadrature.
This quadrature is used in the SG-2 and SG-3 standard grids, for instance.
Franchini et al proposed an atomic partitioning in Journal of Computational Chemistry 2013, 34, 1819–1827
where the atomic partition function has a very simple form:
The same partitioning was also used in their paper on Coulomb fitting J. Chem. Theory Comput. 2014, 10, 1994.
Hirshfeld partitioning could be easily implemented in IntegratorXX. Here, the weight function is defined by the ratio of electron densities
Messages such as the following are less-than-meaningful:
/home/dbwy/Development/IntegratorXX/test/gausslegendre.cxx:42824: FAILED:
REQUIRE_THAT( wgt[i], Catch::Matchers::WithinAbs(ref_wgt[i],w_tolerance) )
with expansion:
0.0000203866 is within 0.0 of 0.1000203866
The following information should be printed in a sane manner
Delley published high order integration schemes on the unit sphere in J. Comput. Chem. 17, 1152 (1996).
These schemes are similar to Lebedev's schemes, but are determined to higher numerical precision. I have also heard that Delley's grids do not suffer from the problem of negative weights, which affects some Lebedev grids.
I have asked Delley by email whether he would be willing to contribute his grids.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.