Move to LibECPInt

[andysim]

Description

This will upgrade our ECP engine from the native implementation to the superior standalone library, also written by Rob Shaw. That library also provides gradient and Hessian terms, which will greatly improve performance.

Todos

• Hook up external build
• Get energies working
• Get gradients working
• Get Hessians working
• Remove old native implementation
• Make all tests / references to ECPs conditional depending on libECP activation. LAB: conditional build marked by USING_ecpint compile definition. This turned out uglier than should strictly be needed because as soon as you ifdef the ao_ecp declarations in integrals.cc, all calcs, w/ or w/o ecp segfault. Possibly I'm missing something, but present ifdef pattern is working stably.
• Linux and Mac conda packages are available off psi4/label/dev. Windows untested. These packages have stylistic problems, but they'll do for now. There are also Linux and Mac conda packages off conda-forge. At least the Linux works fine. Installation is mildly involved with the mixed base channels, so the psi4 packages exist for simplicity of availability.
• NYI message on stability code added to close Cannot use stability analysis #2577
• build docs will be in a separate PR

Questions

• Can anyone think of how one might access code with ecp ifdefs without building a psi4 BasisSet (besides zerobasis)? Particularly though some non-driver mintshelper API call? This is important because the nearly sole whoa-your-mol+basis-needs-ecp-but-libecpint-not-built warning happens at the BasisSet build in export_mints. If users get around that, they'll be surprised by missing electrons b/c ecp code is deactivated.

Checklist

• Tests added for any new features
• All or relevant fraction of full tests run

Status

• Ready for review
• Ready for merge

[andysim]

I just pushed this up so that others can play with it - there has been a lot of interest in getting it implemented. The current system is very janky from a build point of view. I currently just conda install pugixml -c conda-forge first, and build LibECPInt as a standalone library. From there, you can hardwire the location to those libraries by something awful like cmake .. -DCMAKE_CXX_FLAGS="-L/Users/simmonettac/programming/libecpint/installdir/lib -lFaddeeva -lecpint /Users/simmonettac/opt/anaconda3/envs/psi4dev/lib/libpugixml.a -I/Users/simmonettac/programming/libecpint/installdir/include -I/Users/simmonettac/opt/anaconda3/envs/psi4dev/include". As things stand, the implementation looks very straightforward but the answers don't appear to be correct.

[andysim]

Thanks for putting the library together, @robashaw. The API you wrote looks very clean and straightforward to use. Within a shell I appear to be getting the correct nonzero structure, but if I compare your hand-rolled version in Psi4 vs. the results from the new library, I seem to be off by a non constant factor. Some combinations of l and n also seem to be giving nonzero results in the new code, but not in the old version. It's almost certainly something silly that I'm doing - if you get a few minutes to skim the code I added to ecpint.cc and see if anything looks obviously wrong, I'd greatly appreciate it. No problem if you have too much going on. I tried a couple of normalization schemes for the coefficients, but that doesn't seem to be the issue.

[loriab]

neat, thanks Andy! I started looking at the build/packaging for libecpint last week, so hopefully I can keep up with your interface progress.

[andysim]

Thanks for helping with that, Lori. I think the CMake and libmints integration tasks are pretty well separated at this point. The stuff I added to external was mostly placeholder so please overwrite at will - you have write access to my branch. It looks like some upstream changes to the CMake config file naming scheme might be needed, so I just decided to build a standalone lib and worry about integration later - therefore you won't be interfering with anything if you feel like changing the CMake stuff.

[robashaw]

The libecpint results have been very thoroughly tested against near-exact numerical integration routines, so I would be surprised if they are wrong. I will have to take a closer look at what you've done so far to know whether the problem is bugs in the original code or a change in the way certain things are handled. I can clarify though that libecpint does no normalization of the basis functions.

I do know the reason there are now some non-zero quantities though - I was screening them assuming everything would eventually be transformed to spherical gaussians! So it won't have affected results as long as cartesians weren't being used.

@PeterKraus has also noted a problem with linking against pugixml, and I'm not sure what the problem is, but hopefully we'll work out how to fix it soon.

That said, it is worth noting that I have almost finished wrapping a python package of libecpint, so that might prove an easier root to install the library?

[andysim]

Thanks for the rapid response, @robashaw. I have no doubts about the correctness of the code - it's something wrong with my implementation for sure. The biggest mystery was really the nonzero entries, so I'm glad that has a simple explanation. I'll double check the normalization.

[robashaw]

If you need some help working out what's going wrong, please send over the test output you're getting. I've just looked through your changes and I'm not seeing anything that looks out of place.

I remember there being one bug that I fixed in libecpint that will have been in the Psi4 code, but it will have only been noticeable with a specific class of ECPs (those that have an n=1 term (or n=-1 depending on your convention). The line is

psi4/psi4/src/psi4/libmints/ecpint.cc

Line 473 in d16d02e

tooSmall = intValues(l, i) < tolerance;

and should be changed to

tooSmall = tooSmall || ( intValues(l, i) < tolerance );

[andysim]

I'm surprised you're still awake at this time! I applied the fix that you suggested to the old code, but nothing changed. I have attached an output.txt that matches the debug printing in the latest commit (sorry for the mess). The overall results are in the lines starting with "native" and "libecp" in the first column, with some intermediate quantities in the indented lines above those. I'm still adding more debug info to get to the source of the problem. The structure of the nonzero values looks good, so I think it's just normalization or something similar at this point. No obligation to go through this mess of course, but any insights you might have are very greatly appreciated :) The old code values are the numbers being dumped into the buffer right now, so the test passes. This is just the scf-ecp test case for now.

[robashaw]

This is the output I get from running the same system directly through libecpint's main API (with normalised contraction coefficients in the orbital basis). All the results seem to agree with the output you have labelled as `native'. This suggests to me that something is going wrong when the libecpint containers are being created, but I'm struggling to see what exactly it is. It might be worth comparing what your ECP object for argon looks like with one created by libecpint by loading from file. output.txt

[andysim]

That's a great suggestion - I'll give it a shot. It's possible that the orbital basis normalisation isn't correct, but I'll try the high level API as a first step. Thanks so much for your help.

[andysim]

Ok, so progress! It appears that at some point in the past I mixed up the angular momentum and n in the data structures - at least it seems that way when I compare to the test cases. If I switch them, I can get exact agreement between the old and new codes. The pesky nonzero values where the old code gave zeros seem to be problematic, regardless of pure/Cartesian basis.

[robashaw]

I've finally managed to get your branch to compile on my machine (Libint was throwing a fit like it always does whenever I want to do anything), so hopefully I can get to the bottom of the non-zero values soon. I get the correct SCF energy with libecpint through my (non-Psi4) code, so it has to be a wrapping issue. Edit: see below - it was indeed a wrapping issue.

[robashaw]

There are two problems now that I can see your code more closely. The minor one is that when you switched n and l, you forgot to move the +2. But, you don't actually need that +2 anyway (libecpint and Psi4 use the same convention for values of n nowadays). Offending line is here:

https://github.com/andysim/psi4/blob/a6bc397b3ce31130bbef862fc993d11bf8c23afb/psi4/src/psi4/libmints/ecpint.cc#L678

The much bigger problem is that all the things you're passing to libecpint are getting handled as Type1 ECP integrals when they should be Type2. This is because the internal psi4 format is clearly to have the ECPs as separate angular momentum shells, but libecpint assumes that all the shells are in a single ECP object and the highest angular momentum one (as per convention) is the Type1 integral.

There are two possible solutions: make the ECP objects in the way that libecpint is expecting them (I've hacked this in to your code by checking whether the center changes), or I can add a Type1/Type2 flag to the libecpint::ECP object. I'd have to think about whether that would affect anything, but I don't think it would, and might be the more generally useful answer.

Anyway the good news is once I fixed these two things, everything agrees with the `native' results and those I get from the libecpint main API :)

[andysim]

Thanks so much for taking the time to look at this, @robashaw, and for the Libint2 struggles - it's quite a beast that we're only just learning to tame on this end. We have conda packages for some platforms and I think the plan it so have them for all platforms in the long run, but thanks for persevering and getting it compiled locally. The +2 that's present in my version is the result of my desperate hack to try and get things to match, so I'll make sure I change that. I wasn't aware of the the way things are grouped in the new code, so I'll make that change also. I think that using the centre to identify common terms like a very good solution, and you don't need to worry about writing any new type1/type2 classes as far as I'm concerned. Thanks again for your time - I would surely be floundering for the next few weeks without your help.

[PeterKraus]

@robashaw One thought I had while looking at this was: How hard would it be to put pugixml, the associated libecpint internal database reads, and the ecp files themselves behind a CMAKE keyword within libecpint? That way we might be able to keep the basis sets in one place instead of having a duplicated copy within libecpint, and also can ditch the pugixml dependency. I'm not sure this is a helpful functionality for other QM codes, or how much work it'd involve.

[andysim]

That's a good idea, Peter. As far as I know, the XML stuff is only for the high level API, which we are not using. We still parse the basis sets with the same parser and use the same internal data structures - there's just a simple conversion to the library's data structures and then they are handed off. So your concern about duplicate basis set libraries is not a concern, but it would be nice to be able to disable the high level API at compile time, to remove the XML dependency. It's not a huge problem if not - the XML lib is easily obtained from conda forge.

[robashaw]

This is definitely possible and I am currently looking into doing it. CMake gives me headaches though 😭 Also I need to be careful to make sure I don't break anyone else's build involving the library.

[robashaw]

Just to let you know you can now turn off pugixml in libecpint via the option LIBECPINT_USE_PUGIXML

[susilehtola]

This would be really nice to get in Psi4. I don't think there is any free quantum chemistry program that can do geometry optimizations with ECPs. I think in principle PySCF can compute gradients, but it doesn't have a geometry optimizer; I just tried out pyberny and geometric and neither worked with modern Python 💩

[maxscheurer]

LGTM. Thanks for the detailed comments on ECPint component ordering 👍

[fevangelista]

Maybe in the future these could return unique_ptrs directly.

[JonathonMisiewicz]

There's already an issue about this: #2493

[fevangelista]

It's so cool we get both gradients and hessians. Thanks Andy for adding this and Rob for providing the library!