Add Libcint bindings #137
Add Libcint bindings #137
Conversation
|
Some progress notes: For non-segmented contractions, I can evaluate all basic integrals (overlap, nuclear-attraction, kinetic-energy, electron-repulsion) correctly.
Finally, segmented contractions don't work and cause libcint to segfault, so I will fix that once I've got the orderings all figured out. @PaulWAyers also mentioned other functionality from libcint he wanted me to bind. Can get I get those written down? (See the manual.) |
|
One issue I still have (I think...) is the ordering of the spherical integrals. Libcint gives component spherical ordering as in page 5 of its Programmer's Manual. For example, the ordering of the d orbitals in Libcint is xy, yz, z^2, xz, x^2-y^2. GBasis gives the spherical component ordering in the Notes in terms of m and l, with strings of "c{n}" and "s{n}". I'm not sure how to translate between the two notations here. I've also read the IOData page, but I don't understand well enough what to do. I've left room for putting in this reordering [1], [2]. |
|
The direct Cartesian notation is sometimes used instead of the more instructive notation in terms of solid harmonics. One has to use Eq. (6) in the notes to construct it. One can work through the math, as done in Eq. (2.56) of my notes on 1-electron atoms The easiest thing to do, probably, is to run the hydrogen atom in a huge basis set (ccpv8z should suffice; Ne with ccpv9z would go even further) with a program whose conventions are known and a program whose conventions are not known, then compute the overlap. This can be used to deduce the sign convention for the "unknown" program. This would be a useful utility to have, though it's not failsafe (one may need to put a constellation of point charges around the atom to break the symmetry and resolve the degeneracy satisfactorially.) The general convention, though, is:
This isn't quite enough to specify the sign convention for the orbital basis, however. According to the PySCF documentation they use the Condon-Shortley convention, which is the same thing we use in |
|
The latest commit fixes ordering conventions for spherical and cartesian orbitals. I did, however, add this note to the code: This can easily be verified and fixed up if needed, once we fix problem (1.) below. What remains is:
|
|
I fixed the ordering conventions now. Libcint uses (S{l,l-1} S{l,1} ... C_{l,0} ... C_{l,l}), EXCEPT for l=1, in which case it goes (X Y Z), or (C_{1,1} S_{1,1} C_{1,0}). I forced the Libcint build to use (Y Z X). |
|
I fixed the crash with generalized contractions! I also added proper normalization for spherical basis sets, and I believe those are correct now. We can then compare numbers: things look good, but there are some differences between GBasis' answer and Libcint's answers... I will add in the more complicated Cartesian normalization constants soon. I'm wondering if I can use Paul and Farnaz, if you would like to add more integrals and gradients, please see this section of the Libcint README (and the script, for examples), and list the functions you would like me to add. |
|
We use different recursions from If you've linked to the integrals that are currently supported in
We would like to be able to support analytical gradients too. This means the integrals that appear in PySCF for this would be needed. To me, it seems like the key terms are
In the longer run, even analytic Hessians would be interesting, requiring the additional terms here. Not all of these integrals would be required (certainly not for now, as some of them are relevant for PySCF functionality that I doubt we'll ever support). So I'd suggest making a new issue for the analytic gradient/Hessian integrals @msricher |
|
@PaulWAyers OK, I might need some help with identifying which integrals are which, though. Using https://github.com/sunqm/libcint/tree/master#generating-integrals you can see what the autocode symbols mean. Here's the autocode, which gives the names with the symbolic expressions: https://github.com/sunqm/libcint/blob/master/scripts/auto_intor.cl And the header file which lists all the included integrals: The integrals are not given English names, and so I have a hard time knowing which is, for example, the angular momentum integral... |
|
@PaulWAyers I just added the functionality for Gradients and Hessians. I'm not sure what I'll test them against... probably PySCF. I still need some help with interpreting which Libcint functions we want, so could we schedule a meeting soon maybe, or something else? |
|
My interpretation is that the first 4 moments (powers of x,y,z) are generated by The point-charge integral is, I think The derivative integrals we support in gbasis are only supported in combinations (e.g., dot and cross products of momentum operators) in The angular momentum and momentum integral (for order = 1) is I think: |
OK, great! I will add this one.
We'll need to discuss how I will name these, but I have the machinery for this implemented.
It's the notation for complex numbers in Common Lisp. #C(X, Y) = X + iY. |
|
My bindings are complete now. I just have a few issues.
Apart from this, everything is in good shape. I just added in all of the libcint optimizers to make the functions faster, it's easy to add functions on the Python side, and my libcint build script allows for transparent generation of arbitrary integrals using their auto-code feature. |
|
We should figure out the Cartesian normalization. I think it is in the (PDF) notes; admittedly the code is not (yet) documented the way it should be. Where you notice lapses in documentation, please fix it. Where is a complex number appearing? I was not expecting complex results anywhere. |
|
The GBasis momentum, angular_momentum, and moment functions all return complex arrays. Same for Libcint. In GBasis, there's also already If I try to do it myself, then I get lost in the notes. Going from page 2:
I'm not even sure how to get the dimensions right. Where does |
|
Oh, I think I figured it out. The ordering is just wrong now... |
|
Some functions with Cartesian basis still fail. It's the cc-pVDZ/ANO-RCC ones, i.e. the ones where there are D orbitals. I'm not sure why, since it appears it's not just an ordering issue there. |
|
Is it true for all of the basis sets with d orbitals? Does 6-31G(d,p) fail too, for example? |
|
6-31G(d,p) also fails. So in Libcint, it says Cartesian GTOs are un-normalized. I'm only applying the radial part of the normalization to get the simpler things working. But it seems this isn't enough, going to these other basis sets. |
|
I think I need to apply this second part: n = np.power(
np.prod(factorial2(2 * shell.angmom_components_cart - 1), axis=1) / \
factorial2(2 * shell.angmom - 1),
0.5,
)But it has shape (L,), and I'm not sure how to apply it to the coefficients of shape (K, M) |
|
I have no idea what's happening here. tester.py: import numpy as np
from pyscf import gto
from gbasis.parsers import make_contractions, parse_nwchem
from gbasis.integrals.libcint import CBasis
from gbasis.integrals.overlap import overlap_integral
# Disable basis normalization otherwise the basis contraction
# coefficients are scaled to meet basis normalization
gto.mole.NORMALIZE_GTO = False
ps_basis = gto.M(atom="C 0 0 0", unit="Bohr", basis="ccpvdz")
ps_int = ps_basis.intor("int1e_ovlp_cart")
atcoords = np.asarray([[0., 0., 0.]])
atsyms = ["C"]
gb_basis = make_contractions(parse_nwchem("tests/data_ccpvdz.nwchem"), atsyms, atcoords)
gb_int = overlap_integral(gb_basis, coord_type="cartesian")
gc_basis = CBasis(gb_basis, atsyms, atcoords, coord_type="cartesian")
gc_int = gc_basis.olp()
print("GBASIS")
for row in gb_int:
print(", ".join(map("{:12.5e}".format, row)))
print("CBASIS")
for row in gc_int:
print(", ".join(map("{:12.5e}".format, row)))
print("PYSCF")
for row in ps_int:
print(", ".join(map("{:12.5e}".format, row)))python ./tester.py: (G|C)Basis both use PY PZ PX convention, while PySCF uses PX PY PZ, but otherwise the conventions are the same. These matrices are just not the same. |
|
So I am a big confused.
|
|
It works for S and P, not for D and F. |
|
D and F not working makes sense because it's where N_2(a) (from equation 2 in the notes) starts to differ from 1. I think I need to apply this second part: n = np.power(
np.prod(factorial2(2 * shell.angmom_components_cart - 1), axis=1),
-0.5,
)But it has shape (L,), and I'm not sure how to apply it to the coefficients of shape (K, M) Originally posted by @msricher in #137 (comment) |
|
Does Eq. (4) help? I would have assumed we used the normalization from this document. |
|
I'll try that. I didn't have this document. |
|
So, (2 a_x - 1)!! (2 a_y - 1)!! (2 a_z - 1)!! = (2L - 1)!!. I looked at Libcint, and it says to leave the CGBFs unnormalized, but there's a flag to normalize them too. I tried a few contraction normalizations to no success. I also can't figure out what is happening here.... none of these are the same. The only different convention is we have PY PX PZ while PySCF has PX PY PZ. |
|
I looked around a bit more, and saw this issue (linked here to a specific comment). This is what it means for the Cartesian integrals to not be normalized in Libcint... not that you don't need to normalize (just) the radial part, but that the result itself is unnormalized due to how Libcint works. This is fine though, since the SCF will still work and you can normalize the coefficients later using the S matrix. This is how PySCF operates. This weirdness is due to Libcint going from sphericals to Cartesian, almost the opposite of how GBasis works, I guess. Is this fine to leave it be? I would recommend testing the Cartesian integrals with the mean-field repo, once that's ready. |
|
If you talk to Shuoyang and Marco you can use the mean-field repository as we have it right now. I think it works. I'll add you to it. Right now it only has QC-SCF I think. |
Ensure that we get a specific version such that the patch is guaranteed to be applied correctly (and for API/result consistency)
for more information, see https://pre-commit.ci
leila-pujal
left a comment
leila-pujal
left a comment
There was a problem hiding this comment.
@msricher and I have rebased the PR into master including all the changes that were introduced to some functions when #140 and #111 were merged. I have reviewed the code and checked that all the tests pass (except for the angular momentum ones). @FarnazH, on my side, I approve the changes. Let me know if you would like to wait for @marco-2023 and @Ali-Tehrani reviews or if you would be good to merge this into master with only mine.
|
@leila-pujal One more concern... The Github actions do not compile the C library before running the tests (and it seems it fails before even getting to Pytest). Do we want the C library and my tests to be run by Github Actions? Or is it OK as-is? |
|
you are right @msricher. I think we can still merge the PR because you added the scripts to build the libcint library so the functionality can be used once you clone the repo, install it and execute the scripts. We should include github actions to build the library but maybe that would be easier when all the PR related to github actions are merged, like #151. @Ali-Tehrani, feel free to comment on this. |
|
It's better for this to be merged before and then I can add testing it with GIthub actions within that pull-request #151. It also makes it easier for me to fix the merge conflicts. I'm assuming that libcint is optional and if so, I'll have to figure out how to exclude it's test in that case. |
|
Yes, building Libcint is only required for |
|
btw @Ali-Tehrani @leila-pujal, the pyproject.toml (which doesn't seem to work yet?) must ensure that the |
|
For the angular momentum, the right values are: @msricher can you say which states are which in the data you are printing? I presume that the three columns are Lx, Ly, Lz, but I'm not sure the order of elements in the rows. |
|
Yes, it's Lx, Ly, Lz.
…On Wed, Jan 17, 2024 at 3:29 PM Paul W. Ayers ***@***.***> wrote:
For the angular momentum, the right values are:
$$ \begin{align} -i \braket{y|x \partial_y - y \partial_x|x} &= +i \\ -i
\braket{x|x \partial_y - y \partial_x|y} &= -i \\ +i \braket{x|z \partial_x
- x \partial_z|z} &= -i \\ +i \braket{z|z \partial_x - x \partial_z|x} &=
+i \\ -i \braket{z|y \partial_z - z \partial_y|y} &= +i \\ -i \braket{y|y
\partial_z - z \partial_y|z} &= -i \end{align} $$
@msricher <https://github.com/msricher> can you say which states are
which in the data you are printing? I presume that the three columns are
Lx, Ly, Lz, but I'm not sure the order of elements in the rows.
—
Reply to this email directly, view it on GitHub
<#137 (comment)>, or
unsubscribe
<https://github.com/notifications/unsubscribe-auth/ADLODOSGAN2MT32JOMMJB43YPAYCFAVCNFSM6AAAAAA7ADYGUCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQOJWGY2DAMBTGY>
.
You are receiving this because you were mentioned.Message ID:
***@***.***>
|
|
As per Paul's comment in #149, maybe we should just leave out the angular momentum function for now? I will push a commit adding a NotImplementedError onto that function, and then it should be ready to merge if the reviews are ok. |
|
@leila-pujal @Ali-Tehrani @marco-2023 can this be merged? I could start a new PR for the nuclear gradients once it's been merged. |
|
I am good with this. If no one has an objection to it I will merge this at the end of today. @Ali-Tehrani, @marco-2023, @FarnazH, @PaulWAyers. |
5 participants
Steps
Description
This is a draft pull request for adding a C interface (to Libcint) to GBasis.
Libcint (and most other C(++) codes) work by having a list of atoms, and a list of basis shells, where the atom to which a shell belongs is specified with an index (see Libcint programmer's manual). In GBasis (and IOData), I noticed that this concept is not followed, and is combined into one List[GeneralizedContractionShell] structure. This is fine, but there are some problems:
1
The
Attributessection of the GeneralizedContractionShell documentation indicates that it should have achargeattribute, but this seems to be unused (it doesn't exist).2
The IOData
Shellclass has this attribute:@leila-pujal and I have worked on exposing
icenterin the GeneralizedContractionShell object, but I might need more help with that since I'm not sure where it's supposed to have been assigned, when running the test I wrote.(Note: this has been fixed. There was confusion about whether certain quantities (atomic {numbers,coordinates}) were accessible through GBasis objects or passed as arguments to the integral functions.)
Remaining tasks
Figure out how to give to my bindings theicenterandchargeattributes of each atom/shell it needs.Can anyone give me some advice on this? (Also it's my first time making a PR, so I hope I followed the guidelines well enough.)
Type of Changes
Related Issue