Monday, 26 March 2018

Computational chemistry allow properties of molecules to be determined with incredible accuracy. This includes how much energy is stored in a molecule (and how much is released when you combust it), at which frequencies does it vibrate (which can be used to identify chemicals), how quickly a reaction will occur etc. . One of the most widely used method to calculate these properties as a first approximation is the Hartree-Fock (HF) method. Many very successful widely used density functional methods (B3LYP, PBE0, B97) will add in exact exhange calculated using HF to improve the calculation. One of the classic texts in computational chemistry is Modern Quantum Chemistry by Szabo and Ostlund. The book includes a simple calculation of the molecule $HeH^+$ which has the same electronic configuration as the hydrogen molecule but due to the +2 charge on the helium nucleus it is not symmetrical and must be calculated iteratively by converging on a consistent set of parameters which is why the method is often called self consistent field theory (SCF). I have rewritten this example in the python language based on the original Fortran code and try to explain how the calculations are performed and why the operations are done. I have also added in figures to show what the results look like.

We begin with the time-independent Schrodinger equation for a molecule is

$$H\Psi = E\Psi$$

The Hamiltonian $H$ performs some mathematical operations on the function $\Psi$ and returns a number the energy and the $\Psi$ function unchanged. This equation is what is called a eigenvalue equation and turns out to be fairly straightforward to solve particularly if you make use of matrix algebra which is the method used in with the Roothmaan-Hall equations.

You will notice that there are no units because we are making use of atomic units (a.u.) which simplifies the equations by equating these constants to one $m_e=e=h/2\pi=1/(4\pi\epsilon_{0})=1$.

The sums $\Sigma$ have electron $i,j$ and nuclei $A,B$ indices. With the electron positions $\boldsymbol{r}$ and nuclei position $\boldsymbol{R}$ being used to calculate distances between nuclei-electrons $r_{Ai}$, nuclei-nuclei $r_{AB}$ and electron-electron $r_{ij}$. Atomic charges are needed $Z_{A/B}$ in our case $Z_{He}=+2$ and $Z_{H}=+1$.

Using classical mechanics we can determine the kinetic energy from the velocity $\boldsymbol{v}$ and mass $m$ of the particle (multiplied together to give the momentum $\boldsymbol{p}=m\boldsymbol{v}$) $$E_k=T=\frac{1}{2}m\boldsymbol{v}^2=\frac{\boldsymbol{p}^2}{2m}$$.

This classical definition cannot be applied in quantum mechanics as you cannot define a definite position $\boldsymbol{r}$ due to the uncertainity principle which states that the momentum and position cannot be known to any degree of accuracy.

We must make use of the quantum mechanical momentum operator $\hat{p} =-i\hbar\nabla$ which comes from the wave description of a free particle. $$T_e = \frac{\boldsymbol{p}^2}{2m} = \frac{(-i\hbar\nabla)^2}{2m} = -\frac{\hbar^2}{2m}\nabla^2$$ ($\nabla$ is the symbol used to specify derivatives in multivariable calculus. i.e $\partial /\partial x$)

Think of this as a kinetic energy pressure which due to quantum uncertainity holds the electron away from the nuclei allowing it not to collapse to the positive nuclei.

This fluctuation is dependent inversely on the mass so the fluctuation in the light electron is much larger than the heavy nuclei which means it is only neccessary (in most cases) to consider the qunatum mechanical fluctutations of the electrons around fixed nuclei.

this is an attractive energy between nuclei and electrons and is just coulombs law which is usually written $E_{Coulomb} = \frac{q_A q_B}{4\pi \epsilon_0R_{AB}}$ with two charge $q_A$ and $q_B$ a distance apart $R_{AB}$

remember atomic units removes the $4\pi\epsilon_0$ constants, where $\epsilon_0$ is the electric permittivity.

electron-electron repulsion similar to the electrostatic term between the nuclei and electrons.

Most often, however, we only consider the electronic problem for a given set of nuclear positions $\boldsymbol{R}$ as the electrons are significantly lighter and therefore can be considered to instantanously adjust for any nuclear position. (However, the energy here does not include the nuclear repulsions which need to be added at the end to get the total energy.)

One of the challenges in solving this equation is that you do not know where the electron resides only where it is statisically likely to be if you took a measurement - called the many-body effect. This makes the terms involving the electron repulsions ($+\sum_{i>j}\frac{1}{r_{ij}}$) difficult to calculate.

In the Hartree-Fock approach the first approximation used is to treat each electron separately interacting with an averaged distribution of the other electrons (mean field approach).

The functions chosen to represent each electron is based on the hydrogen atomic solutions, which are exact, for an electron around an atom of hydrogen the atomic orbitals.

So we consider a product of one-electron wavefunctions $\psi_i(\boldsymbol{r})$ interacting with a mean field of all the other electrons. $$\Psi(\boldsymbol{r}) = \psi_1(\boldsymbol{r}_1)\psi_2(\boldsymbol{r}_2)\psi_3(\boldsymbol{r}_3)...=\prod_i \psi_i(\boldsymbol{r}_i)$$

However this does not have the correct symmetry requirements which means if you swap electrons the sign of the wavefunction inverts. This requirement means that the Fermions (electrons) cannot occupy the same quantum states and do not collapse towards the nucleus. We can construct something called a slater determinant which does have the correct symmetry properties.

What then is $v_{HF}(i)$ we can consider two interactions a coulomb interaction $\mathscr{J}$ (repulsion between the "averaged out" one-electron orbitals) and an exchange interactions $\mathscr{K}$ which is not classical and is due to electrons of the same spin being indistinguishable which lowers the energy of the system.

$$ v_{HF}(i) = \sum_j \left[ 2\mathscr{J}_j - \mathscr{K}_j\right] $$

Where the sum is over all of the other $j$ electrons. There is a double contribution from the coulomb operator as we are considering closed shell systems where there are two electrons per orbital but the exhange operator is only considered once as this interaction can only be between electrons with like-spins. This means for our HeH+ system it will not contribute as we have

The first two terms are often considered as the core terms as they are the kinetic energy and potential energy for an isolated system such as the hydrogen atom and therefore they are refered to as $\mathscr{H}^{CORE}$

The most successful method to construct these one-electron wavefunctions $\psi_i$ is to consider them delocalised over the whole molecule. Therefore there will be a set of orbitals, one for each electron, which are spread over the whole molecule. We call this type of treatment molecular orbital theory. The question then is what functions do we use to treat these one-electron delocalised orbitals. We know the solutions exactly for a hydrogen like atom so we can use atomic orbitals $\phi_i(\boldsymbol{r}_i)$ centred at each atom $\boldsymbol{r}_i$ to be the basis for our delocalised molecular orbitals. We need to be able to optimise the amounts of each of these atomic orbitals in each of our one-electron delocalised orbitals so we consider a linear combination of atomic orbitals (LCAO).

where the constants $c_i$ are going to be optimised to minimise the energy.

The method relies on the vartionial principle which is true for the Hartree-Fock equations which is that any approximation you make (i.e. describing the electrons as a LCAO) will raise the energy above the real energy of the system. The advantage of a variational method is that you can systematically improve the wavefunction and improve your result knowing that you will be able to converge upon a value (which will be above the actual energy).

The atomic orbitals are called the slater type orbitals for the 1s orbital this can be written as

# Need to import some libraries to have the maths functions and plottingimportnumpyasnpimportscipy.specialasspimportmatplotlib.pyplotasplt%matplotlib notebook # Allows plotting in the notebook
x=np.linspace(-5,5,num=1000)r=abs(x)zeta=1.0psi_STO=(zeta**3/np.pi)**(0.5)*np.exp(-zeta*r)plt.figure(figsize=(4,3))plt.plot(x,psi_STO)

Out[1]:

[<matplotlib.lines.Line2D at 0x1065af0b8>]

Slater type orbitals (STO) are the exact solutions for the hydrogen atom and provide an accurate basis set for many electron molecules however the calculations of the integrals are expensive as their is no simple exact solution for the integrals. One way around this is to approximate the Slater type orbitals using a sum of contracted Gaussian functions (CGF). There are simple analytical expressions for the integrals between two Gaussians so this can save a lot of computing time. Let's look at this for the case of the 1s orbital

# Coeff is the d_n variable in the equation aboveCoeff=np.array([[1.00000,0.0000000,0.000000],[0.678914,0.430129,0.000000],[0.444635,0.535328,0.154329]])# Expon is the alpha variable in the equation aboveExpon=np.array([[0.270950,0.000000,0.000000],[0.151623,0.851819,0.000000],[0.109818,0.405771,2.227660]])psi_CGF_STO1G=Coeff[0,0]*(2*Expon[0,0]/np.pi)**(0.75)*np.exp(-Expon[0,0]*r**2)psi_CGF_STO2G=Coeff[1,0]*(2*Expon[1,0]/np.pi)**(0.75)*np.exp(-Expon[1,0]*r**2) \
+Coeff[1,1]*(2*Expon[1,1]/np.pi)**(0.75)*np.exp(-Expon[1,1]*r**2) \
+Coeff[1,2]*(2*Expon[1,2]/np.pi)**(0.75)*np.exp(-Expon[1,2]*r**2)psi_CGF_STO3G=Coeff[2,0]*(2*Expon[2,0]/np.pi)**(0.75)*np.exp(-Expon[2,0]*r**2) \
+Coeff[2,1]*(2*Expon[2,1]/np.pi)**(0.75)*np.exp(-Expon[2,1]*r**2) \
+Coeff[2,2]*(2*Expon[2,2]/np.pi)**(0.75)*np.exp(-Expon[2,2]*r**2)# Plot the three functionsplt.figure(figsize=(5,3))plt.title("Approximations to a STO with CGF")plt.plot(x,psi_STO,label="STO")plt.plot(x,psi_CGF_STO1G,label="STO-1G")plt.legend()plt.figure(figsize=(5,3))plt.plot(x,psi_STO,label="STO")plt.plot(x,psi_CGF_STO2G,label="STO-2G")plt.legend()plt.figure(figsize=(5,3))plt.plot(x,psi_STO,label="STO")plt.plot(x,psi_CGF_STO3G,label="STO-3G")plt.legend()

Out[19]:

<matplotlib.legend.Legend at 0x106ddda58>

As the number of Gaussians is increased the function more closely describes the slater type orbitals. You will also see that nearest the centre (x=0) the approximation is poorest. This region is called the cusp.

We will use this very poor basis set for our HeH+ molecule for ease of calculation but for accurate calculations at the very least a 6-311G(d,p) basis set will used which would have 6 contracted Gaussian functions to describe hydrogen and helium 1s orbital.

The Hartree-Fock equations need to be converted into a tabular form (matrix form) using the basis set of atomic orbitals to allow for a solution to be determined using the computer. We insert our basis set of orbitals $\sum_i c_i \phi^{CGF}_i(\boldsymbol{r}_i)$ (we will drop the superscript and just have the orbitals denoted as $\phi_i$) We are also going to be building up a table by combining different basis set indices so will use the subscripts $\nu$ and $\mu$ to denotes the indices as we build up the table/matrix. We will also replace $(\boldsymbol{r}_\nu)$ with $(1)$ for ease of reading.

This last term $\int d\nu_1 \phi_{\mu}(1)\phi_{\nu}(1)$ is called the overlap integral written as $S_{\mu \nu}$. This term will help to illustrate the population of a table if we want to know how much the 1s orbital on the He interacts with the 1s orbital on the H atom we would calculate the values on the diagonal of the overlap matrix (for our calculation below this is 0.45077041 in the code the variable is called S12). However, the overlap of the basis function with itself is one by definition.

We have another matrix called the Fock matrix $F_{\mu \nu}$ which is made up of the core contributions $\mathscr{H}^{CORE}$ and the effective potential $v_{HF}$ lets see how the matrices are constructed from these.

Where $H^{CORE}_{\mu \nu}$ is another matrix that has to be calculated containing a kinetic energy and potential energy contribution. These are called one electron integrals (as they are the interactions of a single electron in a molecular orbital). In the code the core contributions are broken up into kinetic terms (T11, T12 and T22, note we do not need T21 as the matrix is symmetrical) and potential terms (V11A, V12A, V22A for atom one and V11B, V12B, V22B for the second atom, note the sum is also over the nuclei for this term)

Where $N$ is the total number of electrons and for closed shell molecules we want half of these electrons, due to two electrons residing in each occupied molecular orbital. This integrate is tricky as it involves the impact of exchange which requires that electrons are indistingishable and therefore instead of a matrix we have 3D matrix or a cube of values (often called a tensor in mathematics). We therefore have to introduce two new indices $\lambda$ and $\sigma$ to denote these indistinguishable versions of the electrons. Which creates quite a construction.

We are trying to find out the values of $\boldsymbol{C}$ which minimise the energy however $\boldsymbol{C}$ is on both sides of the equation so we have to iterative solve the equations until the $\boldsymbol{C}$ on both sides are equal and the energy is the lowest possible this is called the self consistent field method or SCF for short.

Computers are best at solving equations in the form of $\boldsymbol{FC}=\boldsymbol{CE}$, however, that overlap integral is a bit of a problem ($\boldsymbol{F} \boldsymbol{C} = \boldsymbol{S} \boldsymbol{C} \boldsymbol{E}$) so we have to do some matrix algebra to rearrange the equation into an equivalent equation $\boldsymbol{F'C'}=\boldsymbol{C'E}$ (where $\boldsymbol{F'}=\boldsymbol{S}^{-1/2}\boldsymbol{F}\boldsymbol{S}^{-1/2}$ and $\boldsymbol{C'}=\boldsymbol{S}^{-1/2}\boldsymbol{C}$ (In the code $\boldsymbol{S}^{-1/2}$ and $\boldsymbol{S}^{1/2}$ are denoted as X and X.T )

defIntgrl(N,R,Zeta1,Zeta2,Za,Zb):""" Declares the variables and compiles the integrals. """globalS12,T11,T12,T22,V11A,V12A,V22A,V11B,V12B,V22B,V1111,V2111,V2121,V2211,V2221,V2222S12=0.0T11=0.0T12=0.0T22=0.0V11A=0.0V12A=0.0V22A=0.0V11B=0.0V12B=0.0V22B=0.0V1111=0.0V2111=0.0V2121=0.0V2211=0.0V2221=0.0V2222=0.0R2=R*R# The coefficients for the contracted Gaussian functions are belowCoeff=np.array([[1.00000,0.0000000,0.000000],[0.678914,0.430129,0.000000],[0.444635,0.535328,0.154329]])Expon=np.array([[0.270950,0.000000,0.000000],[0.151623,0.851819,0.000000],[0.109818,0.405771,2.227660]])D1=np.zeros([3])A1=np.zeros([3])D2=np.zeros([3])A2=np.zeros([3])# This loop constructs the contracted Gaussian functionsforiinrange(N):A1[i]=Expon[N-1,i]*(Zeta1**2)D1[i]=Coeff[N-1,i]*((2.0*A1[i]/np.pi)**0.75)A2[i]=Expon[N-1,i]*(Zeta2**2)D2[i]=Coeff[N-1,i]*((2.0*A2[i]/np.pi)**0.75)# Calculate one electron integrals # Centre A is first atom centre B is second atom# Origin is on second atom# V12A - off diagonal nuclear attraction to centre A etc.foriinrange(N):forjinrange(N):# Rap2 - squared distance between centre A and centre PRap=A2[j]*R/(A1[i]+A2[j])Rap2=Rap**2Rbp2=(R-Rap)**2S12=S12+S_int(A1[i],A2[j],R2)*D1[i]*D2[j]T11=T11+T_int(A1[i],A1[j],0.0)*D1[i]*D1[j]T12=T12+T_int(A1[i],A2[j],R2)*D1[i]*D2[j]T22=T22+T_int(A2[i],A2[j],0.0)*D2[i]*D2[j]V11A=V11A+V_int(A1[i],A1[j],0.0,0.0,Za)*D1[i]*D1[j]V12A=V12A+V_int(A1[i],A2[j],R2,Rap2,Za)*D1[i]*D2[j]V22A=V22A+V_int(A2[i],A2[j],0.0,R2,Za)*D2[i]*D2[j]V11B=V11B+V_int(A1[i],A1[j],0.0,R2,Zb)*D1[i]*D1[j]V12B=V12B+V_int(A1[i],A2[j],R2,Rbp2,Zb)*D1[i]*D2[j]V22B=V22B+V_int(A2[i],A2[j],0.0,0.0,Zb)*D2[i]*D2[j]# Calculate two electron integralsforiinrange(N):forjinrange(N):forkinrange(N):forlinrange(N):Rap=A2[i]*R/(A2[i]+A1[j])Rbp=R-RapRaq=A2[k]*R/(A2[k]+A1[l])Rbq=R-RaqRpq=Rap-RaqRap2=Rap*RapRbp2=Rbp*RbpRaq2=Raq*RaqRbq2=Rbq*RbqRpq2=Rpq*RpqV1111=V1111+TwoE(A1[i],A1[j],A1[k],A1[l],0.0,0.0,0.0)*D1[i]*D1[j]*D1[k]*D1[l]V2111=V2111+TwoE(A2[i],A1[j],A1[k],A1[l],R2,0.0,Rap2)*D2[i]*D1[j]*D1[k]*D1[l]V2121=V2121+TwoE(A2[i],A1[j],A2[k],A1[l],R2,R2,Rpq2)*D2[i]*D1[j]*D2[k]*D1[l]V2211=V2211+TwoE(A2[i],A2[j],A1[k],A1[l],0.0,0.0,R2)*D2[i]*D2[j]*D1[k]*D1[l]V2221=V2221+TwoE(A2[i],A2[j],A2[k],A1[l],0.0,R2,Rbq2)*D2[i]*D2[j]*D2[k]*D1[l]V2222=V2222+TwoE(A2[i],A2[j],A2[k],A2[l],0.0,0.0,0.0)*D2[i]*D2[j]*D2[k]*D2[l]return

In [8]:

defColect(N,R,Zeta1,Zeta2,Za,Zb):""" Takes the basic integrals and assembles the relevant matrices, that are S,H,X,XT and Two electron integrals """# Form core hamiltonianH[0,0]=T11+V11A+V11BH[0,1]=T12+V12A+V12BH[1,0]=H[0,1]H[1,1]=T22+V22A+V22B# Form overlap matrixS[0,0]=1.0S[0,1]=S12S[1,0]=S12S[1,1]=1.0# This is S^-1/2X[0,0]=1.0/np.sqrt(2.0*(1.0+S12))X[1,0]=X[0,0]X[0,1]=1.0/np.sqrt(2.0*(1.0-S12))X[1,1]=-X[0,1]# This is the coulomb and exchange term (aa|bb) and (ab|ba)TT[0,0,0,0]=V1111TT[1,0,0,0]=V2111TT[0,1,0,0]=V2111TT[0,0,1,0]=V2111TT[0,0,0,1]=V2111TT[1,0,1,0]=V2121TT[0,1,1,0]=V2121TT[1,0,0,1]=V2121TT[0,1,0,1]=V2121TT[1,1,0,0]=V2211TT[0,0,1,1]=V2211TT[1,1,1,0]=V2221TT[1,1,0,1]=V2221TT[1,0,1,1]=V2221TT[0,1,1,1]=V2221TT[1,1,1,1]=V2222

Solve the eigenvalue problem using the secular equation $|\boldsymbol{F'}-\boldsymbol{E}\boldsymbol{I}|=0$ giving the eigenvalues $\boldsymbol{E}$ and eigenvectors $\boldsymbol{C'}$ which can be solved by diagonalising $\boldsymbol{F'}$

Calculate the new density matrix $\boldsymbol{P}$ from the matrix $\boldsymbol{C}$

Check to see if the energy has converged if not then head back to step 6 and repeat.

I have labelled the steps in the program below using the tag ######## STEP 1. ########

In [143]:

defSCF(N,R,Zeta1,Zeta2,Za,Zb,G):""" Performs the SCF iterations """Crit=1e-11# Convergence criteraMaxit=250# Maximum number of iterationsIter=0######## STEP 1. Guess an initial density matrix ######### Use core hamiltonian for initial guess of F, I.E. (P=0)P=np.zeros([2,2])Energy=0.0while(Iter<Maxit):Iter+=1print(Iter)######## STEP 2. calculate the Fock matrix ######### Form two electron part of Fock matrix from PG=np.zeros([2,2])# This is the two electron contribution in the equations aboveforiinrange(2):forjinrange(2):forkinrange(2):forlinrange(2):G[i,j]=G[i,j]+P[k,l]*(TT[i,j,k,l]-0.5*TT[i,j,k,l])# Add core hamiltonian H^CORE to get fock matrixF=H+G# Calculate the electronic energyEnergy=np.sum(0.5*P*(H+F))print('Electronic energy = ',Energy)######## STEP 3. Calculate F' (remember S^-1/2 is X and S^1/2 is X.T) ########G=np.matmul(F,X)Fprime=np.matmul(X.T,G)######## STEP 4. Solve the eigenvalue problem ######### Diagonalise transformed Fock matrixDiag(Fprime,Cprime,E)######## STEP 5. Calculate the molecular orbitals coefficients ######### Transform eigen vectors to get matrix CC=np.matmul(X,Cprime)######## STEP 6. Calculate the new density matrix from the old P ########Oldp=np.array(P)P=np.zeros([2,2])# Form new density matrixforiinrange(2):forjinrange(2):#Save present density matrix before creating a new oneforkinrange(1):P[i,j]+=2.0*C[i,k]*C[j,k]######## STEP 7. Check to see if the energy has converged ########Delta=0.0# Calculate delta the difference between the old density matrix Old P and the new PDelta=(P-Oldp)Delta=np.sqrt(np.sum(Delta**2)/4.0)print("Delta",Delta)#Check for convergenceif(Delta<Crit):# Add nuclear repulsion to get the total energyEnergytot=Energy+Za*Zb/Rprint("Calculation converged with electronic energy:",Energy)print("Calculation converged with total energy:",Energytot)print("Density matrix",P)print("Mulliken populations",np.matmul(P,S))print("Coeffients",C)break

In [144]:

defFormG():""" Calculate the G matrix from the density matrix and two electron integals """foriinrange(2):forjinrange(2):G[i,j]=0.0forkinrange(2):forlinrange(2):G[i,j]=G[i,j]+P[k,l]*(TT[i,j,k,l]-0.5*TT[i,j,k,l])defMult(A,B,C_,IM,M):""" Multiples two square matrices A and B to get C """foriinrange(M):forjinrange(M):forkinrange(M):C_[i,j]=A[i,j]*B[i,j]defDiag(Fprime,Cprime,E):""" Diagonalises F to give eigenvectors in C and eigen values in E, theta is the angle describing the solution """# importmath# Angle for heteronuclear diatonicTheta=0.5*math.atan(2.0*Fprime[0,1]/(Fprime[0,0]-Fprime[1,1]))#print('Theta', Theta)Cprime[0,0]=np.cos(Theta)Cprime[1,0]=np.sin(Theta)Cprime[0,1]=np.sin(Theta)Cprime[1,1]=-np.cos(Theta)E[0,0]=Fprime[0,0]*np.cos(Theta)**2+Fprime[1,1]*np.sin(Theta)**2+Fprime[0,1]*np.sin(2.0*Theta)E[1,1]=Fprime[1,1]*np.cos(Theta)**2+Fprime[0,0]*np.sin(Theta)**2-Fprime[0,1]*np.sin(2.0*Theta)if(E[1,1]<=E[0,0]):Temp=E[1,1]E[1,1]=E[0,0]E[0,0]=TempTemp=Cprime[0,1]Cprime[0,1]=Cprime[0,0]Cprime[0,0]=TempTemp=Cprime[1,1]Cprime[1,1]=Cprime[1,0]Cprime[1,0]=Tempreturn

In [145]:

defHFCALC(N,R,Zeta1,Zeta2,Za,Zb,G):""" Calculates the integrals constructs the matrices and then runs the SCF calculation """# Calculate one and two electron integralsIntgrl(N,R,Zeta1,Zeta2,Za,Zb)# Put all integals into arrayColect(N,R,Zeta1,Zeta2,Za,Zb)# Perform the SCF calculationSCF(N,R,Zeta1,Zeta2,Za,Zb,G)return

In [146]:

"""Let's set up the variables and perform the calculations"""globalH,S,X,XT,TT,G,C,P,Oldp,F,Fprime,Cprime,E,ZbH=np.zeros([2,2])S=np.zeros([2,2])X=np.zeros([2,2])XT=np.zeros([2,2])TT=np.zeros([2,2,2,2])G=np.zeros([2,2])C=np.zeros([2,2])P=np.zeros([2,2])Oldp=np.zeros([2,2])F=np.zeros([2,2])Fprime=np.zeros([2,2])Cprime=np.zeros([2,2])E=np.zeros([2,2])Energy=0.0Delta=0.0N=3R=1.4632Zeta1=2.0925Zeta2=1.24Za=2.0Zb=1.0HFCALC(N,R,Zeta1,Zeta2,Za,Zb,G)

The energy calculated is getting close to the Hartree-Fock limited which is not the real energy of the molecule as it can get even lower if the correlation energy is included. These are the many-body effects which are being ignored in the mean-field approximation to speed up the calculations. The most accurate energy for HeH+ is __. So the mean-field approximation is getting very close to the full energy.

The larger basis set provides more electron density on the Helium atom. You can also see how well the cusp is modelled with many basis sets added together.

If we grab the fchk file we can plot the molecular orbitals in three-dimensions using Avogadro software. Here is the occupied molecular orbital HOMO the surface is the isosurface at the 0.02 density value.

Occupied orbital HOMO

And the unoccupied orbital LUMO

We can also plot the electron density at the 0.002 $e$Å$^3$ value which is roughly where the exchange energy is very strong and you have repuslion between the molecules.

You may be considering that these calculations overly complicated for a simple diatomic such as HHe+ and attempting anything larger would be a nightmare. Fortunately many free publically accessible codes are available for these calculations. If you want to get your hands on something right now there is the website called http://molcalc.org that allows you to do calculations on simple molecules right from your web browser. If you want to do some proper calculations I would recommend getting a copy of the software packages NWChem or GAMESS. To construct your own geometries you will need a program like Avogadro which can be used to build molecules and the input text files you need as input to these quantum chemistry codes.

About Me

Jacob is a PhD student in the CoMo group and a member of Churchill College. He completed a Bachelor of Science with First Class Honours in Chemistry and Physics followed by Masters in Chemistry at the University of Auckland (New Zealand). This included research in the areas of ultrafast spectroscopy, Raman spectroscopy, Bayesian data analysis, computational chemistry and microfluidics. Jacob has strong interests in renewable energy, pollution reduction and carbon nanomaterials. He uses physical models and simulations to describe the chemical world and is developing instruments to measure chemical properties. Jacob is currently studying the formation of soot in engines using molecular dynamics and quantum chemistry to look at gas-soot interactions and self-assembly processes within carbon materials. Jacob is also a founding member of New Zealand Christians in Science (www.nzcis.org) which aims to bring together the two cultures of science and faith through academic study.