jjgoings / McMurchie-Davidson

do a simple closed shell Hartree-Fock using McMurchie-Davidson to compute integrals
BSD 3-Clause "New" or "Revised" License
78 stars 17 forks source link

Wrong orbital energy of He with cc-pVTZ #2

Closed erikkjellgren closed 7 years ago

erikkjellgren commented 7 years ago

I have run a calculation of the orbital energies of He by running:

from mmd.molecule import * 
from mmd.scf import * 
from mmd.postscf import * 
import numpy as np

# read in geometry
geometry = './geoms/he.dat'

# init molecule and build integrals
mol = Molecule(filename=geometry,basis='cc-pvtz')
mol.build()

# do the SCF
scf = SCF(mol)
scf.RHF()

A = np.dot(np.transpose(mol.C), np.dot(mol.F, mol.C))
for i in range(len(A)):
    print(A[i,i])

This gives: (-0.91763022293+0j) (0.590153664669+0j) (1.50050938826+0j) (1.50050938826+0j) (1.50050938826+0j) (3.44176528807+0j) (6.35759720056+0j) (6.35759720056+0j) (6.35759720056+0j) (6.35759720056+0j) (6.35759720056+0j) (7.84834677038+0j) (7.84834677038+0j) (7.84834677038+0j) (8.33659446818+0j)

When the same calculation is performed with the Dalton program the following energies are found:

Symmetry 1A . Hartree-Fock orbital energies: -0.91762549 0.63664267 1.50051327 1.50051327 1.50051327 4.60856244 6.35760376 6.35760376 6.35760376 6.35760376 6.35760376 7.84835826 7.84835826 7.84835826

It can be seen that the second orbital energy is quite different (0.590153664669 and 0.63664267). I have only been able to observe differences larger than the SCF convergence when D orbitals are included in the basisset. The orbital energies might be correct when only S and P orbitals are included, maybe. I have no idea right now what the problem is. Here is data.py with cc-pVTZ for He included (rename it from .pdf to .py, github cannot attach .py in issues). data.pdf

jjgoings commented 7 years ago

This is not a bug. DALTON (as well as Gaussian, for cc-pVTZ) is doing the calculation using spherical Gaussians. The mmd code does the calculations with Cartesian Gaussians. Look at the orbital number: mmd had 15 MOs, whereas DALTON has 14 MOs.

In Gaussian, using spherical gaussians: Input:

#p rhf/cc-pvtz

helium

0 1
He

Output:

 Alpha  occ. eigenvalues --   -0.91763
 Alpha virt. eigenvalues --    0.63664   1.50051   1.50051   1.50051   4.60856
 Alpha virt. eigenvalues --    6.35760   6.35760   6.35760   6.35760   6.35760
 Alpha virt. eigenvalues --    7.84836   7.84836   7.84836

Whereas in Gaussian, with Cartesian gaussians (note 6D keyword) gives: Input:

#p rhf/cc-pvtz 6D

helium

0 1
He

Output:

 Alpha  occ. eigenvalues --   -0.91763
 Alpha virt. eigenvalues --    0.59015   1.50051   1.50051   1.50051   3.44177
 Alpha virt. eigenvalues --    6.35760   6.35760   6.35760   6.35760   6.35760
 Alpha virt. eigenvalues --    7.84835   7.84835   7.84835   8.33659

Thus Gaussian agrees with both sets of orbital energies.

erikkjellgren commented 7 years ago

Oh I feel stupid now. I thought that the orbital energies had some physical significans, and thus should have been independent of coordinate system.


Fra: Joshua Goings notifications@github.com Sendt: 11. juli 2017 02:02:00 Til: jjgoings/McMurchie-Davidson Cc: Erik Kjellgren; Author Emne: Re: [jjgoings/McMurchie-Davidson] Wrong orbital energy of He with cc-pVTZ (#2)

This is not a bug. DALTON (as well as Gaussian, for cc-pVTZ) is doing the calculation using spherical Gaussians. The mmd code does the calculations with Cartesian Gaussians. Look at the orbital number: mmd had 15 MOs, whereas DALTON has 14 MOs.

In Gaussian, using spherical gaussians: Input:

p rhf/cc-pvtz

helium

0 1 He

Output:

Alpha occ. eigenvalues -- -0.91763 Alpha virt. eigenvalues -- 0.63664 1.50051 1.50051 1.50051 4.60856 Alpha virt. eigenvalues -- 6.35760 6.35760 6.35760 6.35760 6.35760 Alpha virt. eigenvalues -- 7.84836 7.84836 7.84836

Whereas in Gaussian, with Cartesian gaussians (note 6D keyword) gives: Input:

p rhf/cc-pvtz 6D

helium

0 1 He

Output:

Alpha occ. eigenvalues -- -0.91763 Alpha virt. eigenvalues -- 0.59015 1.50051 1.50051 1.50051 3.44177 Alpha virt. eigenvalues -- 6.35760 6.35760 6.35760 6.35760 6.35760 Alpha virt. eigenvalues -- 7.84835 7.84835 7.84835 8.33659

Thus Gaussian agrees with both sets of orbital energies.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/jjgoings/McMurchie-Davidson/issues/2#issuecomment-314283569, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ALa952qN9aj2TL7m2zJMZ8np2iN-MlKCks5sMrt4gaJpZM4OTjXm.