13
$\begingroup$

I implemented a restricted Hartree-Fock (RHF) calculation in the STO-3G basis set, as described in Szabo and Ostlund's book [1]. I managed to reproduce the energies of all their calculations ($\ce{H_2}$, $\ce{HeH^+}$, $\ce{H_2O}$, $\ce{CO}$, $\ce{CH_4}$ and $\ce{FH}$) with precision. For example, for $\ce{N2}$ I found $$ E_\text{tot}^\text{HM} = -107.49583784393724 \text{ a.u.} $$ while in the book they report the following (p. 192) $$ E_\text{tot}^\text{Szabo} = -107.496 \text{ a.u.} $$ The situation is analogous for all other molecules I tested.

Now I am trying to implement the calculation of forces, in order to perform Born-Oppenheimer molecular dynamics. I will start to compute the gradient of the PES by finite differences (computationally costly but easy to implement). In order to check my implementation I wanted to use Gaussian09 (for the first time!) as a benchmark because we don't find many force calculations on standar books. In order to check that my input file is correct (i.e. the options can reproduce my "dumb" implementation) I performed a single point calculation for the $\ce{N2}$ dimer. Surprisingly I found $$ E_\text{tot}^\text{G09} = -106.765838750 \text{ a.u.} $$ which is quite different from my result. The input file I used is the following:

%NProcShared=2
%Mem=1GB
#P RHF/STO-3G SP Symmetry=None Units=Bohr Pop=None Guess=Core

RHF

0 1
N 2.074 0.000 0.000
N 0.000 0.000 0.000

I tested the same structure on the input file for $\ce{H2O}$:

%NProcShared=2
%Mem=1GB
#P RHF/STO-3G SP Symmetry=None Units=Bohr Pop=None Guess=Core

RHF

0 1
H +1.4305507125e+00 +0.0000000000e+00 +0.0000000000e+00
H -1.4305507125e+00 +0.0000000000e+00 +0.0000000000e+00
O +0.0000000000e+00 +1.1072513982e+00 +0.0000000000e+00

In this case I find $$ E_\text{tot}^\text{G09} = -74.9629400524 \text{ a.u.} $$ wigh Gaussian09, while with my program I have $$ E_\text{tot}^\text{HM} = -74.962937313769473 \text{ a.u.} = E_\text{tot}^\text{Szabo} $$ For $\ce{H2O}$ everything seems to work fine. In a similar way, I get matching results also for $\ce{CH_4}$, $\ce{FH}$ and $\ce{CO}$.

There is a problem with the $\ce{N2}$ molecule (or simply with the input file) I am not considering?

Which is the best set of options for Gaussian09 that can emulate a home-made Hartree-Fock program (written following closely Szabo's book [1])?

[1] A. Szabo and N. Ostlund, Modern Quantum Chemistry, Dover, 1996.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
12
$\begingroup$

Diagonalization of the core Hamiltonian provides usually not the best guess for the SCF procedure to say the least, and thus, by default Gaussian uses a more sophisticated guess obtained by diagonalizing the Harris functional (Guess=Harris). With this default guess one get the same energy as OP reported:

 SCF Done:  E(RHF) =  -107.495842181     A.U. after    5 cycles

Now, Pop=None was a terrible idea, of course. It is always advisable to look at the final orbitals. Everything looks fine for the case of the default guess

 Alpha  occ. eigenvalues --  -15.51806 -15.51612  -1.44283  -0.72249  -0.57311
 Alpha  occ. eigenvalues --   -0.57311  -0.53949
 Alpha virt. eigenvalues --    0.28131   0.28131   1.12344

but results with the core Hamiltonian guess look suspicios

 Alpha  occ. eigenvalues --  -15.44896 -15.44853  -1.34440  -0.64250  -0.60670
 Alpha  occ. eigenvalues --   -0.52980  -0.16921
 Alpha virt. eigenvalues --   -0.07243   0.31047   1.20826

The energy of the first virtual orbital is negative which is odd. It is difficult to say what went wrong with the core Hamiltonian guess, but some clue can be provided if turning the symmetry back on (this is another very sensible default turning which off was a bad idea, I think).

The default guess yields the electronic state with a particular symmetry ($\sigma_{\mathrm{g}}$)

 Initial guess orbital symmetries:
       Occupied  (SGU) (SGG) (SGG) (SGU) (PIU) (PIU) (SGG)
       Virtual   (PIG) (PIG) (SGU)
 The electronic state of the initial guess is 1-SGG.

and converges within this symmetry

 Orbital symmetries:
       Occupied  (SGG) (SGU) (SGG) (SGU) (PIU) (PIU) (SGG)
       Virtual   (PIG) (PIG) (SGU)
 The electronic state is 1-SGG.

But the core Hamiltonian guess yields the electronic state with no symmetry

 Initial guess orbital symmetries:
       Occupied  (SGG) (SGU) (SGG) (SGU) (PIU) (PIU) (PIG)
       Virtual   (PIG) (SGG) (SGU)

and converges within this state with no symmetry

 Orbital symmetries:
       Occupied  (SGG) (SGU) (SGG) (SGU) (PIU) (PIU) (PIG)
       Virtual   (SGG) (PIG) (SGU)
 Unable to determine electronic state:  partially filled degenerate orbitals.

which is again suspicious.

So, there is something wrong with the core Hamiltonian guess for this system (I suspect the basis is too small for the core Hamiltonian guess to yield sensible results), thus, SCF spits out some "strange" numbers. Strictly, speaking there is nothing strange with the numbers. Remember that SCF usually converges to a local minimum, and which particular depends on which part of the wave function space the initial guess placed the system in. And besides, you're not even guaranteed to have a minimum, SCF might converge to a stationary point of any kind.

Improve this answer
$\endgroup$
5
  • $\begingroup$ I know the options I choose are bad for an actual calculation. However, in this case the goal is to have Gaussian09 behave as close as possible to my implementation. In particular I don't use symmetry and I don't perform population analysis (but I agree that looking at least at the orbitals would have been useful!). I thought that the option Guess=Core corresponds to setting the density matrix elements to zero (approximate Fock's operator with the core Hamiltonian), which is what I do. This is not the case? $\endgroup$
    – user23061
    Dec 6, 2015 at 13:56
  • $\begingroup$ @R.M., yes, you're right, the core Hamiltonian guess ignores all the electron-electron interactions, and is thus physically terrible. For that reason it can lead to strange results. $\endgroup$
    – Wildcat
    Dec 6, 2015 at 13:58
  • $\begingroup$ It seemed strange to me that this guess works in my implementation (and apparently also in Szabo's one) and not in Gaussian09. In any case I will try to implement a better initial guess! Thank you. $\endgroup$
    – user23061
    Dec 6, 2015 at 14:01
  • $\begingroup$ @R.M., in your particular case the core Hamiltonian yields the electronic state that is too far away from the global minimum, so that SCF couldn't find it is way to the global minimum. Just increasing the basis set to 6-31G "fixes" the core Hamiltonian guess, so that it yields a much better starting point. $\endgroup$
    – Wildcat
    Dec 6, 2015 at 14:01
  • $\begingroup$ @R.M., that is a bit strange, I agree. But different implementations of the same algorithm could work differently. The devil, as usual, is in the details. $\endgroup$
    – Wildcat
    Dec 6, 2015 at 14:03
0
$\begingroup$

@R.M. I'm guessing that your HF uses direct SCF, while gaussian09 by default uses an iterative (DIIS and EDIIS) procedure to find the minimum. As previously noted, SCF is not guaranteed to find a global minimum, and different minimization procedures can end up in different minima, even starting from the same spot (Guess=Core). I Recall this being an issue when I wrote an HF program last year in matlab and tried to benchmark results against gaussian. You can definitely turn it off, I believe that SCF=(NoDIIS, DM) (I can never remember the exact gaussian synax) should make it more like your home written program

Improve this answer
$\endgroup$
1
  • $\begingroup$ Direct SCF and DIIS/EDIIS are different things. Direct SCF refers to the way in which integrals are generated and stored. DIIS/EDIIS refers to how convergence of the SCF equations is achieved. They are completely independent of each other. $\endgroup$
    – pentavalentcarbon
    Sep 21, 2016 at 13:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.