- Email: [email protected]
Use of contour integral method in molecular orbital theory
CHEMICI\LPHYSICSLETTERS
Volume 5. number 3
USE
OF
CONTOUR
INTEGRAL
METHOD
15 March 1970
IN MOLECULAR
ORBITAL
THEORY
JAN LINDERBERG Department
of Chemistry,
Aauhus
University,
Denmark
Received 12 January 1970
Caulsan’s contour integration method has been implemented in matrix form. It is found to be effective for the calculation of the charge and bond order matrix and total energies. Orbital energies may alsobc
dctcrmincd.
Coulson [l] introduced an elegant integration technique whereby he could evaluate the results of molecular orbital calculations without direct reference to the eigenvalues and eigenvectors of the secular problem. His method has not enjoyed the popularity it deserves. High speed computing devices have become available with matrix diagonalization routines that may seem to make Coulson’s method obsolete. It is the purpose of this note to point out that his technique is particularly useful in conjunction with computers and large scale calculations. The amount of informationthat is desired from a molecular orbital analysis is very limited in relation to the number of operations that are required to obtain eigenvectors and eigenvalues in any large molecule, particular when overlap is included. We will
while matrix elements of the effective ian in the atomic basis are k?-s = s UjA,ffU, .
PYS +-t&S
=4-s
erical illustrations. Molecular orbital theory, as developed in the approximation of linear combinations of atomic orbitals, leads to a secular problem that will be written as
of the complex
Eere we have denoted by cyk the coefficient the rth atomic orbital in the kth molecular tal, *k =F
Overlap integrals
a matrix
134
R,
zR,&)
is included_
the elements
of which
are
functions
variable 2, through the relations = 6,-, +F
&&s(s)
-
(6)
There are no solutions to eq. (6) whenz equals one of the values ek - LY, but we can derive the spectral representation RY&)
(7)
=.+kZ+Of_ek4k~
dsk = Csk + ’ t
(3)
%tCtk
*
(8)
and find simple poles ae the singularities. The orbitals are normalized and we have that
molecular
- &s,
(5)
with
are defined as
%-s =sv*,QhJ
- @4-Y, - @SYS.
The further development will start with the resolvent related to the problem (1). We define
of orbi-
(2)
%%-k-
(4)
The parameter (Y is to be chosen such that it is larger than any of the orbital energies of occupied molecular orbitals and less than the orbital energy of any virtual orbital. It is easily seen that for the Hffckel method applied to hydrocarbons these notations conform to the conventional ones. The p-matrix defined in (5) is not symmetric when overlap
(1)
dv .
They define together the matrix elements S,, through the system of equations
show here that Coulson’s method only requires the solution of a relatively small number of linear equation systems. First we present the method in matrix form and proceed later to num-
&k- 0)c)crk = F &sCsk-
hamilton-
Volume 5, number 3
CHEiNICALPHYSICSLETTERS
1s &farch 1970
By means of the substitution
9 %kd,k =‘YS-
(9)
We find in the case of no overlap that the charge and bond order matrix is simply related to the residues at the poles of the resolvent:
Y =BY#,
8>0
we obtain the a!ternative
(19)
representation
occ qr = formal charge in orbital r = 2 $
IQ/~,
(10)
occ prs = bond order between Y and s = 2 ck c rk c*sk *
(211
(11)
(22)
Similarly we find in the case with overlap that the population [2] of the orbital Y is occ !?r = a F
c&;k,
!W
$
[c&&
+ d3_kc:k
1.
(13)
By means of the Cauchy residue theorem Coulson and Longuet-Higgins [4] derived the formulae qr = 1 - (l/i;) s-“, R,(iY)
dy ,
Pys = -U/274 _fIw [R&9
+ $&91
(14)
dy. (15)
These equations are fundamental for the method. The imaginary part of the resolvent R(iy) will not contribute to the integrals in eqs. (14) and (15). The real part, Rl(iy), is obtained from the equations
Ek = ff + a/.%,
. (10
It follows that only positive y-values need to be considered. We carried out the integration by a GaussMehler quadrature. In order to find a suitable number of points we investigated in detail one term of eq. (7). It holds that (l&j-_:
(a - Ek) [Y2 + (0!-Ek)2]-1 dy = sgn(Q -Ek) , (17)
where sgn(x) =
1
x >o;
I -1
X co.
(18)
(25)
where we find from the Gauss-MehIer sgn(x) = :El
procedure
[a+ b cos(2j-l)ii/2tt]-1.
(261
This approximative representation of the sign function is plotted in fig. 1 and it can be seen that for 0.2 < Ixl ,C 5 and n 2 8 the error is less than lb. The corresponding interval For n = 4 is 0.33 < x < 3. Good accuracy in the numerical quadrature can be expected when (Y and p are chosen such that the orbital energy falls in the intervais (c+- B/xn, a - &J and (oz t@,, e +p,‘x,J, where the particular value xl1depends upon the curacy.
3, = P rs +~[~B_BVt]&!W
(241
and when
number of integration
-y2R :,(i
(23) between the cases when
Ek=CY+xp
OCC =
quadrature oc-
Ek = a, & p. There is symmetry
and the bond order [3] is firs
The ideal case for a numerical curs for b = 0 or when
points and the desired
ac-
The previous discussion can now be directly continued for the evaluation of charges and bond orders. Most applications of the HSckel method can be carried out for or and p being the standard carbon values. The more general case requires only a rough idea of the energy spectrum. A computer program will generally have these as input parameters that might be changed. ;‘he choice of Q! is also a clue to the orbital energies. We have that the total number of electrons is given by N= Tqr.
(27)
It is a function of 01 which has a step of 2 units at each point (Y = E,$. The step is not quite sharp in the approximate integral form but has the shape of the functions in fig. 1. We have calculated this 135
Volume 5. number 3
CHEMICAL PHYSICS LETTERS
15 March 1970
I
N 6
s jnhd
4
0
-. 4
-.2
0
.2
+zpe
-4
Fig. 2. Approximation to the eignfunction from eq. (26). function
in the HBckel approximation can
be deternined
with an accuracy
of
p/l000 from the figure. Eigenvectors are found from the change in the charge and bond order matrix. Knowledge of the charge and bond order rratrix gives the possibility of calculating the total energy and several other ground state properties. Other quantities, for instance mutual polarizabilities, can be calculated from the resolvent as well [4]. Our calculations
were performed
nal to a GE 230/35 system
on a termi-
in the BASIC programming language which permits algebraic operations with matrices. The computing time for one run on a 6 x 6 matrix is about 4 sec.
136
4
-L-f%
Fig. 2. Number of electrons as a function of the “ohemical potential” (Y. Example from the HUckel problem for furan.
for furan
and the result is displayed in fig. 2. The orbital energies
d2&c
Coulson’s
integration procedure
offers a conven-
ient alternative to the standard diagonalisation methods used in molecular orbital theory and is numerically direct without the iterative character that otherwise may cause problems
REFERENCES [l] C. A. Coulson, Proc. Cambridge Phil. Sot. 36 (1940) 201; J. Chem. Sot. (1954) 3111. [ZJ R. S. Mulliken. J. Chem. Phys. 23 (1955) 1833. [3] B.H. Chirgwin and C. A. Coulson. Proc. Roy. Sot. A201 (1950) 197. [4] C. A. Coulson and H. C. Longuet:Higgins. Proc. Roy. Sot. A191 0947) 39.
Volume 5. number 3
USE
OF
CONTOUR
INTEGRAL
METHOD
15 March 1970
IN MOLECULAR
ORBITAL
THEORY
JAN LINDERBERG Department
of Chemistry,
Aauhus
University,
Denmark
Received 12 January 1970
Caulsan’s contour integration method has been implemented in matrix form. It is found to be effective for the calculation of the charge and bond order matrix and total energies. Orbital energies may alsobc
dctcrmincd.
Coulson [l] introduced an elegant integration technique whereby he could evaluate the results of molecular orbital calculations without direct reference to the eigenvalues and eigenvectors of the secular problem. His method has not enjoyed the popularity it deserves. High speed computing devices have become available with matrix diagonalization routines that may seem to make Coulson’s method obsolete. It is the purpose of this note to point out that his technique is particularly useful in conjunction with computers and large scale calculations. The amount of informationthat is desired from a molecular orbital analysis is very limited in relation to the number of operations that are required to obtain eigenvectors and eigenvalues in any large molecule, particular when overlap is included. We will
while matrix elements of the effective ian in the atomic basis are k?-s = s UjA,ffU, .
PYS +-t&S
=4-s
erical illustrations. Molecular orbital theory, as developed in the approximation of linear combinations of atomic orbitals, leads to a secular problem that will be written as
of the complex
Eere we have denoted by cyk the coefficient the rth atomic orbital in the kth molecular tal, *k =F
Overlap integrals
a matrix
134
R,
zR,&)
is included_
the elements
of which
are
functions
variable 2, through the relations = 6,-, +F
&&s(s)
-
(6)
There are no solutions to eq. (6) whenz equals one of the values ek - LY, but we can derive the spectral representation RY&)
(7)
=.+kZ+Of_ek4k~
dsk = Csk + ’ t
(3)
%tCtk
*
(8)
and find simple poles ae the singularities. The orbitals are normalized and we have that
molecular
- &s,
(5)
with
are defined as
%-s =sv*,QhJ
- @4-Y, - @SYS.
The further development will start with the resolvent related to the problem (1). We define
of orbi-
(2)
%%-k-
(4)
The parameter (Y is to be chosen such that it is larger than any of the orbital energies of occupied molecular orbitals and less than the orbital energy of any virtual orbital. It is easily seen that for the Hffckel method applied to hydrocarbons these notations conform to the conventional ones. The p-matrix defined in (5) is not symmetric when overlap
(1)
dv .
They define together the matrix elements S,, through the system of equations
show here that Coulson’s method only requires the solution of a relatively small number of linear equation systems. First we present the method in matrix form and proceed later to num-
&k- 0)c)crk = F &sCsk-
hamilton-
Volume 5, number 3
CHEiNICALPHYSICSLETTERS
1s &farch 1970
By means of the substitution
9 %kd,k =‘YS-
(9)
We find in the case of no overlap that the charge and bond order matrix is simply related to the residues at the poles of the resolvent:
Y =BY#,
8>0
we obtain the a!ternative
(19)
representation
occ qr = formal charge in orbital r = 2 $
IQ/~,
(10)
occ prs = bond order between Y and s = 2 ck c rk c*sk *
(211
(11)
(22)
Similarly we find in the case with overlap that the population [2] of the orbital Y is occ !?r = a F
c&;k,
!W
$
[c&&
+ d3_kc:k
1.
(13)
By means of the Cauchy residue theorem Coulson and Longuet-Higgins [4] derived the formulae qr = 1 - (l/i;) s-“, R,(iY)
dy ,
Pys = -U/274 _fIw [R&9
+ $&91
(14)
dy. (15)
These equations are fundamental for the method. The imaginary part of the resolvent R(iy) will not contribute to the integrals in eqs. (14) and (15). The real part, Rl(iy), is obtained from the equations
Ek = ff + a/.%,
. (10
It follows that only positive y-values need to be considered. We carried out the integration by a GaussMehler quadrature. In order to find a suitable number of points we investigated in detail one term of eq. (7). It holds that (l&j-_:
(a - Ek) [Y2 + (0!-Ek)2]-1 dy = sgn(Q -Ek) , (17)
where sgn(x) =
1
x >o;
I -1
X co.
(18)
(25)
where we find from the Gauss-MehIer sgn(x) = :El
procedure
[a+ b cos(2j-l)ii/2tt]-1.
(261
This approximative representation of the sign function is plotted in fig. 1 and it can be seen that for 0.2 < Ixl ,C 5 and n 2 8 the error is less than lb. The corresponding interval For n = 4 is 0.33 < x < 3. Good accuracy in the numerical quadrature can be expected when (Y and p are chosen such that the orbital energy falls in the intervais (c+- B/xn, a - &J and (oz t@,, e +p,‘x,J, where the particular value xl1depends upon the curacy.
3, = P rs +~[~B_BVt]&!W
(241
and when
number of integration
-y2R :,(i
(23) between the cases when
Ek=CY+xp
OCC =
quadrature oc-
Ek = a, & p. There is symmetry
and the bond order [3] is firs
The ideal case for a numerical curs for b = 0 or when
points and the desired
ac-
The previous discussion can now be directly continued for the evaluation of charges and bond orders. Most applications of the HSckel method can be carried out for or and p being the standard carbon values. The more general case requires only a rough idea of the energy spectrum. A computer program will generally have these as input parameters that might be changed. ;‘he choice of Q! is also a clue to the orbital energies. We have that the total number of electrons is given by N= Tqr.
(27)
It is a function of 01 which has a step of 2 units at each point (Y = E,$. The step is not quite sharp in the approximate integral form but has the shape of the functions in fig. 1. We have calculated this 135
Volume 5. number 3
CHEMICAL PHYSICS LETTERS
15 March 1970
I
N 6
s jnhd
4
0
-. 4
-.2
0
.2
+zpe
-4
Fig. 2. Approximation to the eignfunction from eq. (26). function
in the HBckel approximation can
be deternined
with an accuracy
of
p/l000 from the figure. Eigenvectors are found from the change in the charge and bond order matrix. Knowledge of the charge and bond order rratrix gives the possibility of calculating the total energy and several other ground state properties. Other quantities, for instance mutual polarizabilities, can be calculated from the resolvent as well [4]. Our calculations
were performed
nal to a GE 230/35 system
on a termi-
in the BASIC programming language which permits algebraic operations with matrices. The computing time for one run on a 6 x 6 matrix is about 4 sec.
136
4
-L-f%
Fig. 2. Number of electrons as a function of the “ohemical potential” (Y. Example from the HUckel problem for furan.
for furan
and the result is displayed in fig. 2. The orbital energies
d2&c
Coulson’s
integration procedure
offers a conven-
ient alternative to the standard diagonalisation methods used in molecular orbital theory and is numerically direct without the iterative character that otherwise may cause problems
REFERENCES [l] C. A. Coulson, Proc. Cambridge Phil. Sot. 36 (1940) 201; J. Chem. Sot. (1954) 3111. [ZJ R. S. Mulliken. J. Chem. Phys. 23 (1955) 1833. [3] B.H. Chirgwin and C. A. Coulson. Proc. Roy. Sot. A201 (1950) 197. [4] C. A. Coulson and H. C. Longuet:Higgins. Proc. Roy. Sot. A191 0947) 39.