Arbitrary-order multipole integrals (and gradients) with McMurchie-Davidson #2496

maxscheurer · 2022-03-23T22:49:18Z

Description

This is the second PR of the M-D series, this time replacing the OS code for arbitrary-order multipole integrals.
Some notable points:

I added a new routine ao_multipoles to the Py API to conveniently call those integrals.
There's also a new multipole_grad which computes first derivatives of arbitrary-order multipole integrals (new feature!).
The dipole derivative code in DipoleInt (which uses l2 for dipole ints, used OS86 for derivs) is removed.
To keep the API stable, the dipole_grad is hard-wired to multipole_grad with the appropriate arguments.
➡️ more features, less code 💡

ToDos

Implement arbitrary-order multipole integrals using M-D
first derivatives
finite-difference tests
more equation numbers
small benchmark against OS86 implementation

Questions

Question1

Checklist

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

Status

Ready for review
Ready for SQUASH-merge (wait for @andysim and @loriab approval)

JonathonMisiewicz · 2022-03-24T14:27:01Z

Why do we have DipoleInt as a class unrelated to MultipoleInt?

maxscheurer · 2022-03-24T15:14:59Z

Why do we have DipoleInt as a class unrelated to MultipoleInt?

I don't know. There's also QuadrupoleInt and TracelessQuadrupoleInt which both use L2. Some more senior developer would have to comment on whether a) we should keep all of the multipole classes as-is or b) hard-wire some to the general MultipoleInt or c) remove the specialized integral classes entirely. ping @loriab @andysim @jturney

jturney · 2022-03-24T15:26:59Z

Back in the early days, we didn't have a generic multipole integral code. Things were just coded up as they were needed. It was quick to code up DipoleInt to get the job done.

Going with @maxscheurer option (c) sounds good to me as long as it can provide both the regular and traceless variants.

maxscheurer · 2022-03-24T15:31:25Z

That's what I figured, thanks @jturney!

(c) sounds good to me as long as it can provide both the regular and traceless variants.

Generalizing the traceless version for higher-order moments is really tedious IMHO, and I don't see an immediate use case in Psi4 currently. Not even the TracelessQuadrupoleInt class is used anywhere currently. Am I missing something? If the traceless multipole moments are required though, I will look into it again.

andyj10224 · 2022-03-24T15:31:25Z

This PR will be useful for my CFMM code :)

maxscheurer · 2022-03-24T15:32:15Z

This PR will be useful for my CFMM code :)

I know 😄

andyj10224 · 2022-03-24T15:35:21Z

This PR will be useful for my CFMM code :)

I know 😄

Also looks like I do not need to change any of my multipole integral calls in my CFMM code too. Double bonus :)

zachglick

This is fantastic Max! I really like the generalized MultipoleInt class. I'm good with options (b) and (c). I think that DipoleInt and QuadrupoleInt classes which function as light wrappers around MultipoleInt are more user-friendly than requiring users to pluck the appropriate integrals out of the MultipoleInt return.

Have you done any performance comparisons between the new MD code and the old OS code? I don't know if one is expected to be faster than the other. It would be good to do some simple timings (maybe one low angmom system and one high angmom system?) before completely ditching the OS code.

psi4/src/psi4/libmints/dipole.h

zachglick · 2022-03-24T15:19:16Z

psi4/src/psi4/libmints/mintshelper.h

    /// Vector AO Traceless Quadrupole Integrals
    std::vector<SharedMatrix> ao_traceless_quadrupole();
+    /// Vector AO Multipole Integrals up to given order
+    std::vector<SharedMatrix> ao_multipoles(const std::vector<double>& origin = {0., 0., 0.}, int order = 1);


I like this general ao_multipoles function a lot! Could we get a comment here explaining the order that the multipoles are returned in? Would it be [x, y, z, xx, xy, xz, yy, yz, zz, xxx, ...]? You might've already written that comment somewhere else.

I like the general function, too. However, order shouldn't have a default value; the caller should explicitly specify it.

Good point! But then, origin won't have a default value either if I leave the order of function arguments unchanged. Should we go for non-default values in both cases? I think this would be the most clean solution.

Then I'd switch the order. More than likely the caller will want to run it at that default origin.

I'm not sure about that, I wanted to keep the order as in similar mintshelper functions.

Yeah, that sounds good @JonathonMisiewicz .

Ok, so we change the order for the multipole function already here, rest comes later?

Yes. I just created an issue to change the other.

Does having a default origin make sense, for arbitrary order multipoles? The moments are origin dependent from n+1, upwards, where n is the order of the first nonzero moment. So, e.g., the dipole is center-dependent in charged molecules but not in neutral molecules. If we're talking about allowing high orders in this function, the majority of them are going to be origin-dependent. Max already mentioned external potentials, but for things like CFMM one must also specify the expansion center (before then translating that using fancy operators). That's all a verbose way of me trying to say that I don't think either the origin or the expansion order should be given default arguments here.

I'm all for no default values for the generalized functions (should be changed to ao_multipole_potential, too, I can do that when I add M-D for these integrals). If one is not aware that the moments are indeed origin-dependent, better throw an error than returning a probably unintended result IMHO.
For the specialized functions that are exposed via MintsHelper there is not even the possibility to change the origin, but this is fine for now I guess.

psi4/src/psi4/libmints/multipoles.h

jturney · 2022-03-24T15:38:37Z

Hmm, I thought I had added the traceless quadrupole for @lothian years ago. If I'm completely wrong, he can certainly correct me. 😄

maxscheurer · 2022-03-24T15:42:24Z

Hmm, I thought I had added the traceless quadrupole for @lothian years ago. If I'm completely wrong, he can certainly correct me. 😄

You're right, I overlooked ao_traceless_quadrupole on the py side -- I can just code up a hard-wired version of the traceless quadrupole integrals as @zachglick mentioned above. I think that'd be a good compromise.

lothian · 2022-03-24T15:44:14Z

Yes, we use the traceless quadrupole integrals for calculations of ROA spectra in ccresponse.

andyj10224 · 2022-03-24T15:52:21Z

This is fantastic Max! I really like the generalized MultipoleInt class. I'm good with options (b) and (c). I think that DipoleInt and QuadrupoleInt classes which function as light wrappers around MultipoleInt are more user-friendly than requiring users to pluck the appropriate integrals out of the MultipoleInt return.

Have you done any performance comparisons between the new MD code and the old OS code? I don't know if one is expected to be faster than the other. It would be good to do some simple timings (maybe one low angmom system and one high angmom system?) before completely ditching the OS code.

@maxscheurer If you want, I can test CFMM with your new code to see if the multipole calculations are indeed faster. Just let me know.

maxscheurer · 2022-03-24T16:22:06Z

If you want, I can test CFMM with your new code to see if the multipole calculations are indeed faster. Just let me know.

It's okay, I already have a script to benchmark the old code vs. the new one -- will report some of the benchmarks soon 👍

andysim · 2022-03-24T16:23:02Z

Regarding the general case MultipoleInt vs. special routines like DipoleInt. Keep in mind that asking for MultipoleInts with order=2 will give overlap, dipole, quadrupole. Asking for QuadrupoleInts will only give quadrupoles. Computing the extra integrals isn't really a big deal in terms of efficiency, but it might be a little surprising for the user to find that the indexing doesn't start from zero. The current quadrupole integral implementation just calls Libint2 and picks out only the quadrupole components. Going with only MultipoleInts is better for maintenance, but changes the API and could lead to some surprises. However, there isn't really any efficiency penalty for doing that, so I don't really have a strong opinion either way.

andysim · 2022-03-24T16:25:19Z

The conversion to traceless form is very similar to the Cartesian->Spherical transform process and should be applied after the fact, in my opinion. We only really care about the quadrupoles for now so we can stick with that. I have some code somewhere that makes coefficients for the general case, but it's a little bit overkill to implement that, IMHO.

maxscheurer · 2022-03-24T16:41:38Z

Keep in mind that asking for MultipoleInts with order=2 will give overlap, dipole, quadrupole.

That's not correct, the multipole moments start at dipole (to keep the behavior consistent with the previous OS impl).

andysim · 2022-03-24T16:42:51Z

That's not correct, the multipole moments start at dipole (to keep the behavior consistent with the previous OS impl).

Ooops. I should probably have known that, given that I wrote that code in the first place!

maxscheurer · 2022-03-24T16:57:46Z

I think now that we've gathered some opinions on the MultipoleInt/DipoleInt/... splitting, I suggest to move this reconciliation to another PR. I want to move the M-D transition forward as fast as possible, and ripping out L2 out if the existing integral classes seems a bit beyond the scope. I'll open an issue so we don't forget about it in the future.

JonathonMisiewicz · 2022-03-24T17:01:21Z

I'll approve this when Zach's comments are addressed. Otherwise, this looks great!

Thanks a lot for the tests.

maxscheurer · 2022-03-24T23:23:16Z

Small benchmark

The following timings were obtained for H2O from 3 "takes", where in each take, the
integral is repeatedly evaluated 15x. The result shown here is that of the fastest take, divided by the number of repeats (i.e., 15 here). The performance of M-D is actually quite nice -- of course one could tweak it a little bit here and there, but that's not needed right now IMHO.

The code snippet from my benchmark script looks sth like:

takes = 3
repeat = 15
for b in bas:
    for m in [1, 4, 8, 16, 24, 36]:
        basis = psi4.core.BasisSet.build(mol, 'orbital', b)
        mints = psi4.core.MintsHelper(basis)
        key = f"{b}/m = {m}"
        keys.append(key)
        print(key)
        for _ in tqdm(range(takes)):
            with timer.record(key):
                for i in range(repeat):
                    M = mints.ao_multipoles(order=m, origin=[1.0, 2.0, 3.0])
for k in keys:
    best = timer.best(k) / repeat
    print(f"{k:<15} |   {best * 1000:8.2f} ms")

Results

basis set	order	M-D	OS
cc-pvdz	1	0.11 ms	0.20 ms
	4	0.46 ms	0.73 ms
	8	1.87 ms	2.19 ms
	16	11.25 ms	10.88 ms
	24	47.11 ms	43.30 ms
	36	193.20 ms	177.40 ms
cc-pvtz	1	0.27 ms	0.52 ms
	4	1.34 ms	2.19 ms
	8	5.97 ms	7.36 ms
	16	46.39 ms	47.05 ms
	24	172.76 ms	170.22 ms
	36	742.18 ms	725.64 ms
cc-pvqz	1	0.81 ms	1.53 ms
	4	4.16 ms	7.12 ms
	8	20.47 ms	26.56 ms
	16	153.60 ms	168.69 ms
	24	543.86 ms	579.47 ms
	36	2327.06 ms	2436.51 ms
cc-pv5z	1	2.09 ms	4.37 ms
	4	12.59 ms	22.29 ms
	8	63.34 ms	84.20 ms
	16	457.24 ms	522.29 ms
	24	1634.66 ms	1807.03 ms
	36	6654.10 ms	7198.05 ms

andyj10224 · 2022-03-24T23:34:48Z

Seems to me like there is no reason to keep the old OS code. The only case where it is faster is smaller basis sets in high orders, where it is unstable anyway :)

zachglick

Thanks for the very thorough benchmark Max, LGTM!

andysim

Outstanding effort! Thanks for fixing my bug in the M recursion limit. It's really great that this PR deletes code, while adding more stable, faster integrals with derivatives and better testing. LGTM!

andysim · 2022-03-25T14:30:06Z

Heads-up, @edeprince3. You will need these derivative integrals.

…vidson (psi4#2496) * multipole integrals with M-D * fix M matrix * multipole derivatives * fdiff test and some cleanup * more eq numbers * reviews * cleanup, no default origin + full tests gotcha

maxscheurer mentioned this pull request Mar 23, 2022

Libint2 2022 Strategy #2442

Open

maxscheurer force-pushed the multipole_md branch from 1c897f3 to 36b2171 Compare March 24, 2022 12:56

zachglick reviewed Mar 24, 2022

View reviewed changes

maxscheurer mentioned this pull request Mar 24, 2022

Unify MultipoleInt/DipoleInt/QuadrupoleInt #2498

Open

maxscheurer added 6 commits March 24, 2022 23:42

multipole integrals with M-D

98be17a

fix M matrix

8309174

multipole derivatives

caf22fa

fdiff test and some cleanup

c30933d

more eq numbers

a5baa88

reviews

6d4a111

maxscheurer force-pushed the multipole_md branch from 865bbab to 6d4a111 Compare March 24, 2022 22:57

zachglick approved these changes Mar 25, 2022

View reviewed changes

maxscheurer added libmints For things that are wrong with libmints (but not wavefunction). cleanup For issues where the goal is to make Psi4 a little cleaner. labels Mar 25, 2022

cleanup, no default origin + full tests gotcha

5581a43

maxscheurer requested review from JonathonMisiewicz, andysim and loriab March 25, 2022 10:10

JonathonMisiewicz approved these changes Mar 25, 2022

View reviewed changes

jturney approved these changes Mar 25, 2022

View reviewed changes

andysim approved these changes Mar 25, 2022

View reviewed changes

loriab approved these changes Mar 25, 2022

View reviewed changes

maxscheurer merged commit 69bc9ce into psi4:master Mar 25, 2022

maxscheurer deleted the multipole_md branch March 25, 2022 15:47

loriab added this to the Psi4 1.6 milestone Mar 25, 2022

maxscheurer mentioned this pull request Mar 29, 2022

Multipole Potential Integrals with McMurchie-Davidson #2504

Merged

10 tasks

maxscheurer mentioned this pull request Apr 5, 2022

Implement McMurchie Davidson Integrals #2414

Closed

Arbitrary-order multipole integrals (and gradients) with McMurchie-Davidson #2496

Arbitrary-order multipole integrals (and gradients) with McMurchie-Davidson #2496

Uh oh!

Conversation

maxscheurer commented Mar 23, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Description

ToDos

Questions

Checklist

Status

Uh oh!

JonathonMisiewicz commented Mar 24, 2022

Uh oh!

maxscheurer commented Mar 24, 2022

Uh oh!

jturney commented Mar 24, 2022

Uh oh!

maxscheurer commented Mar 24, 2022

Uh oh!

andyj10224 commented Mar 24, 2022

Uh oh!

maxscheurer commented Mar 24, 2022

Uh oh!

andyj10224 commented Mar 24, 2022

Uh oh!

zachglick left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

maxscheurer Mar 24, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

jturney commented Mar 24, 2022

Uh oh!

maxscheurer commented Mar 24, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

lothian commented Mar 24, 2022

Uh oh!

andyj10224 commented Mar 24, 2022

Uh oh!

maxscheurer commented Mar 24, 2022

Uh oh!

andysim commented Mar 24, 2022

Uh oh!

andysim commented Mar 24, 2022

Uh oh!

maxscheurer commented Mar 24, 2022

Uh oh!

andysim commented Mar 24, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

maxscheurer commented Mar 24, 2022

Uh oh!

JonathonMisiewicz commented Mar 24, 2022

Uh oh!

maxscheurer commented Mar 24, 2022

maxscheurer commented Mar 23, 2022 •

edited

Loading

maxscheurer Mar 24, 2022 •

edited

Loading

maxscheurer commented Mar 24, 2022 •

edited

Loading

andysim commented Mar 24, 2022 •

edited

Loading

andyj10224 commented Mar 24, 2022 •

edited by maxscheurer

Loading