psi4 / psi4numpy

Combining Psi4 and Numpy for education and development.
BSD 3-Clause "New" or "Revised" License
337 stars 157 forks source link

AO basis Overlap matrix is printed incorrectly via psi4numpy RHF_Hessian script. #106

Closed AlexB67 closed 3 years ago

AlexB67 commented 3 years ago

Hello,

As stated above, it appears fine in the PSI4 output file however. When the overlap is transformed to the MO basis it appears correctly also (unit matrix).

To reproduce use the script here

https://github.com/psi4/psi4numpy/blob/master/Self-Consistent-Field/RHF_Hessian.py and Insert the lines near the top

S_ao = np.asarray(mints.ao_overlap()) print(S_ao)

Not all 1s along the diagonal as expected, whereas psi4 ouput does as well as my own code while testing Usiing the latest release build from the psi4 donwload page(linux).

Thank you for all the great software.I am learning a lot from it :)

ashutoshvt commented 3 years ago

Hi Alex, when I run S_ao = np.asarray(mints.ao_overlap()) print(S_ao), I do get the proper ao_overlap matrix with 1s on the diagonal.

[[ 1. 0.236703936510848 0. 0. 0. 0.038405592106858 0.038405592106858] [ 0.236703936510848 1. -0. 0. 0. 0.386138791494766 0.386138791494766] [ 0. -0. 1. 0. 0. 0.209726921614688 0.209726921614688] [ 0. 0. 0. 1. 0. 0. 0. ] [ 0. 0. 0. 0. 1. -0.268438218366334 0.268438218366334] [ 0.038405592106858 0.386138791494766 0.209726921614688 0. -0.268438218366334 1. 0.181759849558527] [ 0.038405592106858 0.386138791494766 0.209726921614688 0. 0.268438218366334 0.181759849558527 1. ]]

Could you post the matrix you are getting?

AlexB67 commented 3 years ago

Hi thanks for the prompt reply. Try with 6-31+G** i didn't mention that., ... my bad

options = {'BASIS':'6-31+G**', 'PUREAM':'FALSE', 'GUESS':"CORE", 'E_CONVERGENCE':1e-8, 'D_CONVERGENCE':1e-8 }` and coordinates

mol = psi4.geometry(""" O -0.000000000003 0.000000000000 -0.112929827494 H -1.483747300978 0.000000000000 0.896139156594 H 1.483747301023 0.000000000000 0.896139156604 unit bohr symmetry c1 """)

Best regards.

AlexB67 commented 3 years ago
[afb@localhost ~/psi4data/scripts/psi4/Self-Consistent-Field]$  python RHF_Hessian.py
[[ 1.       0.23369 -0.       0.      -0.       0.16728  0.       0.      -0.       0.03353  0.       0.       0.03353  0.       0.03353  0.07232 -0.       0.      -0.       0.03617  0.0696   0.05003  0.      -0.03403  0.03617  0.0696  -0.05003  0.      -0.03403]
 [ 0.23369  1.       0.       0.       0.       0.76364 -0.       0.       0.       0.54707  0.      -0.       0.54707  0.       0.54707  0.40953  0.       0.       0.       0.25458  0.38326  0.28602  0.      -0.19452  0.25458  0.38326 -0.28602  0.      -0.19452]
 [-0.       0.       1.       0.       0.      -0.       0.50152  0.       0.      -0.       0.      -0.      -0.       0.      -0.       0.       0.16806  0.       0.      -0.2419  -0.12211 -0.21751  0.       0.24863  0.2419   0.12211 -0.21751  0.      -0.24863]
 [ 0.       0.       0.       1.       0.       0.       0.       0.50152  0.       0.      -0.       0.       0.      -0.       0.       0.       0.       0.16806  0.       0.       0.       0.       0.14808  0.       0.       0.       0.       0.14808  0.     ]
 [-0.       0.       0.       0.       1.       0.       0.       0.       0.50152 -0.       0.      -0.      -0.       0.      -0.      -0.       0.       0.       0.16806  0.16451  0.08304  0.24863  0.      -0.02101  0.16451  0.08304 -0.24863  0.      -0.02101]
 [ 0.16728  0.76364 -0.       0.       0.       1.       0.       0.       0.       0.69901  0.       0.       0.69901  0.       0.69901  0.78665  0.       0.       0.       0.43218  0.68774  0.21659  0.      -0.1473   0.43218  0.68774 -0.21659  0.      -0.1473 ]
 [ 0.      -0.       0.50152  0.       0.       0.       1.       0.       0.       0.       0.       0.       0.       0.       0.       0.       0.67035  0.       0.      -0.48396 -0.39657  0.01277  0.       0.18236  0.48396  0.39657  0.01277  0.      -0.18236]
 [ 0.       0.       0.       0.50152  0.       0.       0.       1.       0.       0.       0.       0.       0.       0.       0.       0.       0.       0.67035  0.       0.       0.       0.       0.28092  0.       0.       0.       0.       0.28092  0.     ]
 [-0.       0.       0.       0.       0.50152  0.       0.       0.       1.       0.       0.       0.       0.       0.       0.       0.       0.       0.       0.67035  0.32913  0.2697   0.18236  0.       0.1569   0.32913  0.2697  -0.18236  0.       0.1569 ]
 [ 0.03353  0.54707 -0.       0.      -0.       0.69901  0.       0.       0.       1.       0.       0.       0.33333  0.       0.33333  0.47078  0.       0.       0.       0.43836  0.45101  0.14833  0.      -0.36406  0.43836  0.45101 -0.14833  0.      -0.36406]
 [ 0.       0.       0.      -0.       0.       0.       0.       0.       0.       0.       0.33333  0.       0.       0.       0.       0.       0.       0.       0.       0.       0.       0.      -0.19349  0.       0.       0.       0.       0.19349  0.     ]
 [ 0.      -0.      -0.       0.      -0.       0.       0.       0.       0.       0.       0.       0.33333  0.       0.       0.       0.       0.       0.       0.      -0.18204 -0.03265 -0.13676  0.      -0.01099  0.18204  0.03265 -0.13676  0.       0.01099]
 [ 0.03353  0.54707 -0.       0.      -0.       0.69901  0.       0.       0.       0.33333  0.       0.       1.    

Not all rows and cols shown ,see the 0.33333 entries on diagonal

AlexB67 commented 3 years ago

I see sorry. I wasn't aware even if PUREAM is true ?

Best ...

On Wed, 9 Dec 2020 at 21:15, Justin Turney notifications@github.com wrote:

Hi Alex,

If your basis set has d or higher orbitals then you may not have 1s on the diagonal. For example, the Cartesian dx2, dy2, and dz2 will be 1 but dxy, dxz, and dyz will not. A similar pattern is found the higher angular momentum orbitals. This is simply down to the normalization that psi4 uses.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/psi4/psi4numpy/issues/106#issuecomment-742061537, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMQS3FQOATNY5US3XEH42Q3ST7SGFANCNFSM4UT7T5LQ .

AlexB67 commented 3 years ago

I meant PUREAM false.

On Wed, 9 Dec 2020 at 22:43, Alexander Borro alexander.borro@gmail.com wrote:

I see sorry. I wasn't aware even if PUREAM is true ?

Best ...

On Wed, 9 Dec 2020 at 21:15, Justin Turney notifications@github.com wrote:

Hi Alex,

If your basis set has d or higher orbitals then you may not have 1s on the diagonal. For example, the Cartesian dx2, dy2, and dz2 will be 1 but dxy, dxz, and dyz will not. A similar pattern is found the higher angular momentum orbitals. This is simply down to the normalization that psi4 uses.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/psi4/psi4numpy/issues/106#issuecomment-742061537, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMQS3FQOATNY5US3XEH42Q3ST7SGFANCNFSM4UT7T5LQ .

ashutoshvt commented 3 years ago

Great! I am putting @jturney's email here for visibility:

If your basis set has d or higher orbitals then you may not have 1s on the diagonal. For example, the Cartesian dx2, dy2, and dz2 will be 1 but dxy, dxz, and dyz will not. A similar pattern is found the higher angular momentum orbitals. This is simply down to the normalization that psi4 uses.

Please close the issue if this is resolved.

AlexB67 commented 3 years ago

I am veeery Sorry for the inconvenience.

On further investigation, so much output I apologise. PSI output does agree with psi4numpy. I should have been more careful. My mistake in the input file, so case closed.

Can anyone shed light (point to a resource) how the normalisation works for D orbitals in psi4 ? I can look at the source

In my code they are all 1 on the diagonal. Is this wrong or difference in convention ?

To come full circle I wrote code for evaluating the Hessian when i was testing against psi4. When I get D orbitals my answers do indeed begin to depart from psi4 in some cases (Linear molecules still okay). numerical from analytic gradient agrees however to within a wavenumber for those same test cases.

I always set puream to false since my code works with Cartesians at this time i.e 6d 10f etc. No transformation to spherical.

andysim commented 3 years ago

Psi4 follows the CCA integral standard; the normalization of Cartesians is described at the end of the paper. Good luck!

AlexB67 commented 3 years ago

Thanks for the paper, but alas I do not have access to that resource (Not in university) , unless I buy it. It also made my project extra challenging. :D Still google gives a lot for free that got me this far.

Once upon a time (in the nineties) I was doing all that sexy rotation vibration Hamiltonian stuff with Prof Mills... so thought it would fun to learn a bit of electronic structure theory 30 or more years later ... very rusty ... but getting into the groove again. :)

On Wed, 9 Dec 2020 at 23:31, Andy Simmonett notifications@github.com wrote:

Psi4 follows the CCA integral standard https://onlinelibrary.wiley.com/doi/epdf/10.1002/jcc.20815; the normalization of Cartesians is described at the end of the paper. Good luck!

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/psi4/psi4numpy/issues/106#issuecomment-742129913, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMQS3FTTBIGVQ7YJDGQVXGTSUACGRANCNFSM4UT7T5LQ .

Found in a pdf :D

andysim commented 3 years ago

Kudos to you! That's not something that many people can just pick up without siginificant guidance. I don't think I can post the whole paper in a public place because of copyright issues, but here's the pertinent section: image

Only the basis functions with all quanta in x,y or z will have unit norm as a result.

AlexB67 commented 3 years ago

Thank you kindly sir. My d orbital overlaps now do agree with PS4 after adopting the new normalisation.

Luckily it didn't affect my code wither way. In terms of answers. Everything still comes out the same, energies gradients etc ... as far as I can tell so far anyway.

Now just to fix those pesky second derivative potential ints and all should be set for my fully analytic Hessian.

Cheers :)