Fix Hessian bug

[andysim]

Description

A typo in the overlap integral second derivatives caused errors in the analytic hessians. The error seems to be confined to one of the three contributions to matrix elements where the angular momentum in the bra and ket differ, and only when the derivatives both refer to the same perturbation; which is why the code made it through the initial tests. I'm still trying to wrap my head around exactly why those tests work, so I'd like a day or two before this is ready to merge. My sincere apologies to all whom this bug inconvenienced. Fixes #791 and #901.

Todos

Notable points that this PR has either accomplished or will accomplish.

• User-Facing for Release Notes
  • Fixed a bug in the analytic Hessian code.

Questions

• Which extra tests do you have in mind, @loriab?

Status

• Ready to go

[dgasmith]

Awesome.

[dgasmith]

Hmm, I think you just never tested anything beyond P orbitals. This if statement will not be hit without D or higher.

I wish we could put weights on commits so that we could crank this one up. Three lines changed, but I am sure quite a bit of time looking into it.

[loriab]

Three characters changed, but painful ones to track down, I'm sure.

I'll apply these changes to pyvib2 and see if tests/python/vibanalysis/ magically heals.

A thousand thanks for investigating this.

[andysim]

It was tested with higher angular momentum, which is why I’m confused. Lori and I sat down at the last hackathon to try and find problem cases, but couldn’t reproduce the error. What I’m really struggling with is the origin of the importance of atom order in the test provided in #901. It’s probably related to the fact that I only take derivatives on center A and use translational invariance relations to get the rest, but my flu addled brain is failing me right now. For the same reason, yes, this took me a while to track down. Oh well, hopefully this is the last bug we ever have in the hessians. Or anywhere...

…

On Jan 27, 2018, at 8:16 PM, Daniel Smith ***@***.***> wrote: @dgasmith commented on this pull request. Awesome. In psi4/src/psi4/libmints/overlap.cc: > @@ -336,7 +336,7 @@ void OverlapInt::compute_pair_deriv2(const GaussianShell& s1, const GaussianShel buffer_[(0*size)+ao12] += (4.0*a1*a1*x[l1+2][l2]*y[m1][m2]*z[n1][n2] - 2.0*a1*(2*l1+1)*x[l1][l2]*y[m1][m2]*z[n1][n2]) * over_pf; if (l1 > 1) Hmm, I think you just never tested anything beyond P orbitals. This if statement will not be hit without D or higher. I wish we could put weights on commits so that we could crank this one up. Three lines changed, but I am sure quite a bit of time looking into it. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, or mute the thread.

[loriab]

Looking great, @andysim! Gone from 1 pass, 1 fail w/normco errors in the tenths, and 3 radically fail to 4 pass and 1 fail with normco errors in the ten thousandths. Systems are formaldehyde, ammonia, methane, ethene, carbon dioxide. I think HOOH was the original suspicious one, so I'll dig that up.

[loriab]

Looks wonderful to me

[loriab]

Here's some of the errors in frequencies pre-fix; max is for HOOH TS. I've cherry-picked the commit into pyvib2, and the tests, if handy, are here.

• CO2

	exp:  [ 761.1526  761.1526  1513.3122  2580.1495]
	obs:  [  765.97192061+0.j   954.84604992+0.j  1511.98600598+0.j
  2577.29764966+0.j]
	dif:  [   4.81932061+0.j  193.69344992+0.j   -1.32619402+0.j   -2.85185034+0.j]

• Ethene

	exp:  [ 885.4386  1073.4306  1080.4104  1135.9390  1328.7166  1467.8428
  1565.5328  1831.6445  3287.3173  3312.1889  3371.7974  3399.4232]
	obs:  [  885.43769105+0.j  1073.42980093+0.j  1080.48361108+0.j
  1135.93848940+0.j  1341.57874983+0.j  1467.71701985+0.j
  1565.53257843+0.j  1831.21804884+0.j  3287.31942795+0.j
  3312.16817472+0.j  3373.24482206+0.j  3399.42544494+0.j]
	dif:  [ -9.08947590e-04+0.j  -7.99067003e-04+0.j   7.32110847e-02+0.j
  -5.10595379e-04+0.j   1.28621498e+01+0.j  -1.25780152e-01+0.j
  -2.21567931e-04+0.j  -4.26451164e-01+0.j   2.12794906e-03+0.j
  -2.07252822e-02+0.j   1.44742206e+00+0.j   2.24494358e-03+0.j]

• CH3 & NH3 ok
• Formaldehyde

	exp:  [ 1325.3286  1359.7579  1637.4774  2013.4255  3108.9786  3183.3975]
	obs:  [ 1325.82947096+0.j  1376.06130208+0.j  1637.37151242+0.j
  2010.83034420+0.j  3108.93196717+0.j  3185.06227338+0.j]
	dif:  [  0.50087096+0.j  16.30340208+0.j  -0.10588758+0.j  -2.59515580+0.j
  -0.04663283+0.j   1.66477338+0.j]

• HOOH TS

	exp:  [    0.0000+278.6136j  1128.8155  +0.j      1364.3495  +0.j
  1698.2924  +0.j      4140.0257  +0.j      4146.4749  +0.j    ]
	obs:  [    0.00000000+273.62786912j  1122.20829758  +0.j          1385.41184315  +0.j
  1715.97750485  +0.j          4143.48648522  +0.j          4306.47193962  +0.j        ]
	dif:  [   0.00000000-4.98573088j   -6.60720242+0.j           21.06234315+0.j
   17.68510485+0.j            3.46078522+0.j          159.99703962+0.j        ]

[andysim]

Alrighty, this one should be ready to go now. I tried to cook up a simple test using f orbitals, but they're all too costly; in the end I just added a permutation to the atom ordering in our existing cc-pVDZ water test, because this is enough to reveal the bug. The fact that we've always had a working cc-pVDZ test case shows how subtle the bug is; the affected d components in water must be zero by symmetry. I did have a distorted water in my test suite to check that case, but didn't detect problems. Oh well, live and learn I guess.

This should be a trivial review, and it clearly helps to address a couple of high priority tickets, so please have at it whenever you get a chance. Sorry again for the error. Next time you see me, I'll be at a chalkboard, writing "I will not cut and paste", à la Bart in the Simpsons opening credits. Except, instead of writing it, I'll be cutting and pasting it.

[jturney]

Oh man, I remember debugging this type of typos in other integrals. Kudos to you for finding them.