Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Problem with diatomic molecules where both atoms have orbitals different than s #2

Closed
aromanro opened this issue Oct 18, 2016 · 11 comments
Closed

Comments

@aromanro
Copy link

@aromanro aromanro commented Oct 18, 2016

Hi,

I think I found a problem with such setup. For example for N2 molecule, with the following configuration:

basis = sto-3g
charge = 0
atoms = 2
units = bohr
N 0.0000 0.0000 0.0000
N 2.147932 0.0000 0.0000

the given result is:

Total energy [Hartree]: -144.1995252

This is with the fix of the bug I previously signaled.

According to some paper I was looking into, the Hartree-Fock limit for N2 is -108.994.

I am comparing your results with a program I've written. It has different basis sets than yours but the values are quite similar. I too have issues with diatomic molecules for this kind of configuration. It works nicely for example for H2O or any diatomic molecule I tried that has a hydrogen in it, but fails once I put two atoms with both more than s orbitals.

I'm using recurrence relations (I think they are called Hamilton & Schaefer) to solve the integrals, for this described case I could spot different terms in both the nuclear matrix and in electron-electron integrals than given by your program. Most of the values are basically identical, but some come out differently.

It doesn't mean yours are wrong since I am also getting bad results for such molecules, but it seems that I'm getting results closer to the expected value, for example for the above molecule I get -113.677. Still wrong :(

Thanks,
Adrian

@ifilot ifilot mentioned this issue Oct 18, 2016
@aromanro
Copy link
Author

@aromanro aromanro commented Oct 18, 2016

As a help to diagnose this, I just compared with the values my program gives for N2, the overlap and kinetic matrices are basically identical, despite the different methods of calculating them, so I suppose the problems originate elsewhere. I'll let you know if I found more about this.

@aromanro
Copy link
Author

@aromanro aromanro commented Oct 18, 2016

For the core matrix, I found only one difference - it comes out from the nuclear matrix, since the kinetic matrix is ok. I don't know which one is correct (they could be both wrong), the matrix element corresponding to a px orbital from one N and the px orbital from the other N comes out different, all other matrix values seem to match. There are also electron-electron integral differences :(

@ifilot
Copy link
Owner

@ifilot ifilot commented Oct 18, 2016

I have calculated N2 using the geometry as you described above in Gaussian; the result with a STO-3G basis set is -107.50 a.u.; which is within the HF limit as you mentioned above. HFCXX indeed gives -144.20 a.u., which is clearly wrong. Interestingly, also PyQuante is off with a value of -149.46 a.u.

I am still in the dark where the problem originates from.

@aromanro
Copy link
Author

@aromanro aromanro commented Oct 18, 2016

Could you please look at the value PyQuante gives for the above matrix element (line 3, col 8, numbering starting with 1 or the symmetric correspondent)? I'm curious if I have the wrong value...

@aromanro
Copy link
Author

@aromanro aromanro commented Oct 18, 2016

Your program has -3.44603 there, mine gives 3.487057. If I have it wrong, probably e-e integrals are also wrong in a similar manner, since the calculations have similarities.

@ifilot
Copy link
Owner

@ifilot ifilot commented Oct 18, 2016

I make a mistake with PyQuante; I now get -106.82 a.u. for the total energy (yay!). That number is correct.

Furthermore, PyQuante gives

<2px | hcore | 2px> = 3.487

So I probably have some error in my code or in my definition for the basis set.

For reproduction of the PyQuante calculation, see below:

from PyQuante.Molecule import Molecule
from PyQuante.hartree_fock import rhf

n2 = Molecule('n2',[(7,(0,0,0)),(7,(2.147932,0,0))])

en,orbe,orbs = rhf(n2, basis="sto-3g")
print en
@aromanro
Copy link
Author

@aromanro aromanro commented Oct 18, 2016

Thank you! That should help me a lot!

@aromanro
Copy link
Author

@aromanro aromanro commented Oct 18, 2016

By the way, by looking at your code it seems that you use some code similar with the one discussed here: http://chemistry.stackexchange.com/questions/40958/computing-two-electron-integrals-with-an-sto-3g-basis-set There is an errata here: http://spider.shef.ac.uk/ It's unlikely that you have that mistake, but it's worth a try...

@aromanro
Copy link
Author

@aromanro aromanro commented Oct 19, 2016

Here is more info about the bug: I simply forced H[2][7] = H[7][2] = 3.487; With this change, the code converges to -107.5005469, which seems to be correct. So you have an issue somewhere with the nuclear matrix, I have an issue with the e-e integrals :(

I tried replacing the Boys functions calculations with mine with the same result, so it's not originating from Boys/Fgamma implementation. Anyway, I think I can rely on your code for e-e integrals values to check against mine.

@aromanro
Copy link
Author

@aromanro aromanro commented Oct 21, 2016

Hi,

Now I have more info. I used your program to check the electron-electron integrals results from my program. Now I get basically identical results for the N2 molecule. I had a nasty bug in the electron transfer relation implementation combined with a mistake in the way I stored integral computation results :(

So the issue with the N2 molecule in your program is certain to originate only from the nuclear matrix computation.

ifilot pushed a commit that referenced this issue Dec 14, 2016
ifilot
@ifilot
Copy link
Owner

@ifilot ifilot commented Dec 14, 2016

OK; it turned out to be a really stupid mistake. See:
https://github.com/ifilot/hfcxx/blob/master/src/nuclear.cpp#L65

That line had arrA[iI] = A_term(... instead of arrA[iI] += A_term(...

@aromanro : I am very grateful that you found this mistake. Many thanks again!

@ifilot ifilot closed this Dec 14, 2016
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Linked pull requests

Successfully merging a pull request may close this issue.

None yet
2 participants
You can’t perform that action at this time.