Frozen orbitals¶
DIRAC allows the freezing (or elimination) of orbitals during Hartree-Fock or Kohn-Sham calculations. We shall illustrate this functionality using bonding in the water molecule as an example and essentially reproducing a calculation reported by Jarvie and co-workers in 1973 [Jarvie1973].
The water molecule has a bond angle of 104.5 \(^{\circ}\). This is sometimes rationalized by invoking \(sp^3\) hybridized oxygen. Perfect hybridization would, however, lead to a bond angle of 109.5 \(^{\circ}\) \(=\mathrm{arccos}(-1/3)\), and so the deviation is rationalized by invoking valence shell electron pair repulsion (VSEPR) theory and in particular the stronger repulsion between lone pairs.
Jarvie and co-workers therefore proposed to freeze the oxygen \(2s\) orbital in water and see how this affected the bond angle. Let us do this calculation the DIRAC way:
Atomic calculation¶
First we need an oxygen \(2s\) orbital. We therefore set up an atomic calculation using the input files O.inp
**DIRAC
.TITLE
Oxygen atom
.WAVE F
.ANALYZE
**WAVE FUNCTIONS
.SCF
*SCF
.CLOSED SHELL
4
.OPEN SHELL
1
4/6
**ANALYZE
.MULPOP
*END OF
and O.mol
INTGRL
Water
C 1 2 X Y
8. 1
O 0.0D0 0.0D0 0.0D0
LARGE BASIS cc-pVDZ
FINISH
Our plan now is to import this orbital into a subsequent molecular calculation. There is, however, one important restriction: The atomic calculation has to be carried out in the same symmetry as the molecular one, which in this case will be \(C_{2v}\), hence the explicit symmetry specification in the O.mol input file (for more details, see How to specify the molecule and basis set using xyz files) . We run the calculation using the command:
pam --inp=O --mol=O --get "DFCOEF=cf.O"
Geometry optimization with frozen orbitals: how it works¶
For the molecular calculation we use the input files H2O.inp
**DIRAC
.TITLE
Water molecule.
.OPTIMI
*OPTIMI
.NUMGRA
**WAVE FUNCTIONS
.SCF
*SCF
.CLOSED SHELL
10
.OWNBAS
.PROJEC
1
DFOXXX
1
2
.FROZEN
2
*END OF
and H2O.mol
INTGRL
Water
C 2 2 X Y
8. 1
O .0000000000 0.0000000000 -.2249058930
LARGE BASIS cc-pVDZ
1. 1
H 1.4523499293 .0000000000 .8996235720
LARGE BASIS cc-pVDZ
FINISH
The latter file specifies the experimental geometry of water with bond angle 104.5 \(^{\circ}\) and bond distance 0.972 Å. The former file includes the .NUMGRA keyword to activate geometry optimization using a numerical gradient. This will be needed here since we have included the constraint of a frozen orbital.
Now let us look at the specification of the frozen orbital: We have to tell DIRAC that the orbital was not calculated in the full molecular basis, rather using only the basis of the oxygen atom. This information is indicated by the keyword .OWNBAS. Next, we have to tell where to find the orbital. This is specified using the keyword .PROJECT and its arguments. The oxygen \(2s\) orbital is specified by the coefficients found on the file DFCOEF that we recovered in the atomic calculation and renamed cf.O. The oxygen atom is considered a fragment of the water molecule, and the coefficient file therefore constitute a fragment file. The first argument of the .PROJECT keyword is the specification of the number of fragment files to read, which in this case is just one.
For each fragment we now have to specify the name of the fragment file. As seen in the FORTRAN snippet in the manual under the .PROJECT keyword, the fragment file should be given a six-character name, so we in this case choose DFOXXX. When the fragment has been calculated in its own basis, as specified by the .OWNBAS keyword, we specify it by indicating the number of symmetry-independent nuclei it contains. Let us see how this works: Going back to the molecular input file H2O.mol file we see that it gives the coordinate of one oxygen atom and one hydrogen atom. We do not need to specify the coordinates of the second hydrogen atom because it is then known from symmetry. The two hydrogen atoms are symmetry-dependent and have to be treated together as a fragment. We have, however, selected the oxygen atom as our fragment and specify that our fragment consists of a singe symmetry-independent center. Oxygen will be chosen since it comes first in the molecular input file.
We now have to specify what orbital(s) to select on the fragment file. This is done through an orbital string (see the section Specification of orbital strings). The oxygen atom has five occupied orbitals and the \(2s\) orbital is the second one, as we can for instance infer from looking at the Mulliken population analysis in the output of the atomic calculation.
Now we have to understand the .FROZEN keyword. This will take a little bit more theory: At each iteration in the SCF-cycle we have to solve the generalized eigenvalue equation
where \(F\) is the Fock (or Kohn-Sham) matrix in the AO-basis (symmetry-adapted or not) and \(\mathbf{c}\) and \(\varepsilon\) are the coefficients and eigenvalues we seek. The equation also features the overlap matrix \(S\) since the AO-basis is generally not orthogonal. A first step towards the solution of this equation is to introduce a non-unitary transformation \(V\) that will take us to an orthonormal basis, that is
where \(I\) is the identity matrix. We insert \(VV^{-1}=I\) in front of the coefficients and pre-multiply with \(V^{\dagger}\) to give
Since the matrix \(V\) transforms the overlap matrix to the identity matrix, the transformed equation simplifies to
This is a normal eigenvalue problem and it suffices now to give the transformed Fock matrix \(F^{\prime}\) to a standard routine for matrix diagonalization in order to obtain eigenvalues \(\varepsilon\) and the transformed eigenvectors \(\mathbf{c}^{\prime}\). In a final step we get the MO-coefficients \(\mathbf{c}\) by the backtransformation
This diagonalization procedure is repeated in each iteration of the SCF-cycle and provides the optimization of the molecular orbitals.
The freezing of orbitals is done in the following manner: The imported fragment orbitals are projected out of the transformation matrix \(V\), which means that they are not present in the orthonormal basis of the transformed Fock matrix \(F^{\prime}\). They have in other words been removed from the variational space. If this was our goal we could in fact have stopped here. This is for example what was done in [Fossgaard2003b] where the effect of the lanthanide contraction on bonding in the CsAu molecule was investigated by projecting out the gold \(4f\) orbitals from the variational space (hence the name of the .PROJECT keyword).
In the present case we do not want to delete the oxygen \(2s\) orbital, rather freeze it. Therefore we have to put the oxygen \(2s\) orbital back when backtransforming coefficients. The .FROZEN keyword takes as argument an orbital string which indicates in what position one would like to put the frozen orbitals in the total set of orbitals.
Geometry optimization of the water molecule with a frozen oxygen 2s orbital¶
We run the geometry optmization with the command:
pam --inp=H2O --mol=H2O --put "cf.O=DFOXXX"
After optimization we find that the bond distance is 0.972 Å and the bond angle 95.9 \(^{\circ}\). This is perhaps a surprising result since one would expect that breaking hybdridization would lead to a bond angle of 90.0 \(^{\circ}\). Our result suggests that there must be some additional repulsion between the hydrogens of water. The nature of this repulsion (electrostatic vs. Pauli repulsion) has been a subject of considerable discussion in the literature. A full account is found in [Dubillard2006].