Ab initio program for small molecules utilizing

basis-sets consisting of s-functions (spherical Gaussians) only

copyright: P. Vogt, H. Huber, Dec. 1999

All properties in atomic units!

Input (nuclear data, basis-set data, occupation matrix)

> restart:

> with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

> Digits:=10: numerical accuracy

> NumberofNuclei := 2: Example: hydrogen molecule i.e. 2 atoms (nuclei)

> Dim := 4: 2 basis functions on each nucleus -> 4 basis functions

define nuclear data:

atomic number

X Y Z

> AtomicNumber := array([1,1]);

> NuclearCoordinate := array([[0,0,0],
[1.4,0,0]]); 2. nucleus 1.4 a o along x-coordinate

define basis-set data:

Exponent

X Y Z

> Exponent := array([1.2,0.3,1.2,0.3]); 4 exponents (2 basis functions on each nucleus)

> BasisCoordinate := array([[0,0,0],
[0,0,0],
[1.4,0,0],
[1.4,0,0]]); 2 basis functions per nucleus

Occupation-Matrix in the molecular orbital representation (MO-basis)

> Occupation := array([[2,0,0,0],
[0,0,0,0],
[0,0,0,0],
[0,0,0,0]]); lowest MO occupied with 2 electrons

Parameters and auxiliary functions

> X:=1:Y:=2:Z:=3: makes the program more readable, as x, y ans z-ccordinates are stored in arrays

sqrDistance - Function to calculate the square of the distance between two basis functions

> sqrDistance := proc (func1, func2)

> (BasisCoordinate[func1,X]-BasisCoordinate[func2,X])^2 +
(BasisCoordinate[func1,Y]-BasisCoordinate[func2,Y])^2 +
(BasisCoordinate[func1,Z]-BasisCoordinate[func2,Z])^2;

> end:

ProductGaussian - Calculates the product of two Gaussian functions

The product of two Gaussian fucntions is again a Gaussian with a new exponent cnew,

which is the sum of the original exponents, and new coordinates Rxnew, Rynew, Rznew,

which are found from the old ones by weighting with the exponents;

> ProductGaussian := proc(func1,func2,cnew::evaln,Rxnew::evaln,Rynew::evaln,Rznew::evaln)

> local a,b;

> cnew := Exponent[func1] + Exponent[func2];

> a := Exponent[func1] / eval(cnew);

> b := Exponent[func2] / eval(cnew);

> Rxnew := a * BasisCoordinate[func1,X] + b * BasisCoordinate[func2,X];

> Rynew := a * BasisCoordinate[func1,Y] + b * BasisCoordinate[func2,Y];

> Rznew := a * BasisCoordinate[func1,Z] + b * BasisCoordinate[func2,Z];

> end:

AuxInt - Auxiliary function for the calculation of the potential and the two-electron-integrals

> AuxInt := proc(X)

> if X=0 then 1 else evalf(1/2*sqrt(Pi/X)*erf(sqrt(X))) fi;

> end:

Diagonalisation of a matrix by the Jacobi method

a: Matrix to be diagonalised

d: Eigenvalues

v: Eigenvectors

includes symmetrisation of the matrix at the beginning and

sorting of the eigenvalues and -vectors at the end

> Diagonalisation := proc(a,d,v)
local i,ip,iq,j,c,g,h,s,sm,t,tau,theta,tresh,b,z,Maxiter;
Symmetrisation(a);
Maxiter:=10**Digits;
d:=array(1..Dim);
v:=array(1..Dim,1..Dim);
b:=array(1..Dim);
z:=array(1..Dim);
for ip from 1 to Dim do
for iq from 1 to Dim do
v[ip,iq]:=0;
od:
v[ip,ip]:=1;
od:
for ip from 1 to Dim do
b[ip]:=a[ip,ip];
d[ip]:=b[ip];
z[ip]:=0;
od:
for i from 1 to Maxiter do
sm:=0;
for ip from 1 to Dim-1 do
for iq from ip+1 to Dim do
sm:=sm+abs(a[ip,iq]);
od:
od:
if(sm<10**(-Digits-2)) then break fi;
if(i<4)then
tresh:=0.2*sm/Dim**2;
else
tresh:=0;
fi;
for ip from 1 to Dim-1 do
for iq from ip+1 to Dim do
g:=100*abs(a[ip,iq]);
if ((i>4) and (abs(d[ip])+g=abs(d[ip]))and(abs(d[iq])+g=abs(d[iq]))) then
a[ip,iq]:=0;
elif(abs(a[ip,iq])>tresh) then
h:=evalm(d[iq]-d[ip]);
if(abs(h)+g = abs(h))then
t:=a[ip,iq]/h;
else
theta:=0.5*h/a[ip,iq];
t:=1/(abs(theta)+sqrt(1.+theta**2));
if(theta<0) then t:=-t; fi;
fi;
c:=1/sqrt(1+t**2);
s:=t*c;
tau:=s/(1+c);
h:=t*a[ip,iq];
z[ip]:=z[ip]-h;
z[iq]:=z[iq]+h;
d[ip]:=d[ip]-h;
d[iq]:=d[iq]+h;
a[ip,iq]:=0;
for j from 1 to ip-1 do
g:=a[j,ip];
h:=a[j,iq];
a[j,ip]:=g-s*(h+g*tau);
a[j,iq]:=h+s*(g-h*tau);
od;
for j from ip+1 to iq-1 do
g:=a[ip,j];
h:=a[j,iq];
a[ip,j]:=g-s*(h+g*tau);
a[j,iq]:=h+s*(g-h*tau);
od:
for j from iq+1 to Dim do
g:=a[ip,j]:
h:=a[iq,j]:
a[ip,j]:=g-s*(h+g*tau):
a[iq,j]:=h+s*(g-h*tau):
od:
for j from 1 to Dim do
g:=v[j,ip]:
h:=v[j,iq]:
v[j,ip]:=g-s*(h+g*tau):
v[j,iq]:=h+s*(g-h*tau):
od:
fi:
od:
od:
for ip from 1 to Dim do
b[ip]:=b[ip]+z[ip]:
d[ip]:=b[ip]:
z[ip]:=0:
od:od:
Sort(d,v);

> print(d);

> end:

Sort - Sort the eigenvalues und eigenvectors

To sort the eigenvalues and corresponding eigenvectors after a diagonalisation according to the

size of the eigenvalues.

> Sort:=proc(Eigenvalues,Eigenvect)

> local MaxValue,MaxIndex,column1,column2, buffer;

> for column1 to Dim-1 do column1;

> MaxValue:=Eigenvalues[column1];

> MaxIndex:=Dim+1;

> for column2 from column1+1 to Dim do

> if MaxValue > Eigenvalues[column2] then

> MaxIndex:=column2;

> MaxValue:=Eigenvalues[MaxIndex];

> fi;

> od;

> if MaxIndex < Dim+1 then

> Eigenvect:=swapcol(Eigenvect,column1,MaxIndex);

> Eigenvalues[MaxIndex]:=Eigenvalues[column1];

> Eigenvalues[column1]:=MaxValue;

> fi;

> od;

> evalm(Eigenvalues);

> end:

Symmetrisation - Symmetrizes a matrix A

> Symmetrisation:=proc(A)

> A:=scalarmul(matadd(A,transpose(A)),0.5);

> end:

Integral functions (S, T, V, 2e-integral, M)

Overlapintegral - Function to calculate the overlapinteg rals <j|j>

> Overlapintegral := proc (func1, func2)

> local alpha, beta, cinv,Q ,aux;

> alpha := Exponent[func1];

> beta := Exponent[func2];

> cinv := 1/(alpha+beta);

> Q := exp(-alpha * beta * cinv * sqrDistance(func1,func2));

> aux := (4*alpha*beta*cinv^2);

> Q*sqrt(sqrt(aux^3));

> end:

KineticIntegral - Function to calculate the kinetic integral <j| T |j>

> KineticIntegral := proc (func1, func2)

> local alpha, beta, E;

> alpha := Exponent[func1];

> beta := Exponent[func2];

> E := alpha * beta /(alpha+beta);

> Overlapintegral(func1,func2) * E *(3 - 2*E * sqrDistance(func1,func2));

> end:

PotentialIntegral - Function to calculate the potential integral <j| V |j>

> PotentialIntegral := proc(func1, func2)

> local c,Rx,Ry,Rz,V,argument, nuc;

> ProductGaussian(func1, func2,c,Rx,Ry,Rz);

> V := 0;

> for nuc to NumberofNuclei do

> argument := c*((Rx-NuclearCoordinate[nuc,X])^2 + (Ry-NuclearCoordinate[nuc,Y])^2 + (Rz-NuclearCoordinate[nuc,Z])^2);

> V := V + AtomicNumber[nuc] * AuxInt(argument);

> od;

> evalf(-Overlapintegral(func1,func2) * 2/sqrt(Pi) * sqrt(c) * V);

> end:

TwoElectronIntegral - Function to calculate the 2-e-integral <jj|1/ r |jj>

> TwoElectronIntegral := proc(i,j,k,l)

> local Argument, cnew, c1, c2, Rx1, Rx2, Ry1, Ry2, Rz1, Rz2;

> ProductGaussian(i,j,c1,Rx1,Ry1,Rz1);

> ProductGaussian(k,l,c2,Rx2,Ry2,Rz2);

> cnew := c1 * c2 / (c1 + c2);

> Argument := cnew * ((Rx1 - Rx2)^2 + (Ry1 - Ry2)^2 + (Rz1 - Rz2)^2);

> evalf(2/sqrt(Pi)) * Overlapintegral(i, j) * Overlapintegral(k, l) * sqrt(cnew) * AuxInt(Argument);

> end:

MOperator - 2-e-Operator for closed shells

> MOperator := proc(i,j)

> local Sum, k, l;

> Sum := 0;

> for k to Dim do

> for l to Dim do

> Sum := Sum + DensityMatrix[k, l] * (TwoElectronIntegral(i, j, k, l) - 0.5 * TwoElectronIntegral(i, k, j, l));

> od;

> od;

> Sum;

> end:

Construction of integral matrices

Overlap Matrix S

> SMatrix := array(1..Dim,1..Dim):

> for row to Dim do

> for column to row do

> SMatrix[row,column] := Overlapintegral(row,column);

> SMatrix[column,row] := SMatrix[row,column];

> od:

> od:

> print (SMatrix):

Kinetic Energy Matrix T

> TMatrix := array(1..Dim,1..Dim):

> for row to Dim do

> for column to row do

> TMatrix[row,column] := KineticIntegral(row,column);

> TMatrix[column,row] := TMatrix[row,column];

> od:

> od:

> print(TMatrix);

Potential Energy Matrix V

> VMatrix := array(1..Dim,1..Dim):

> for row to Dim do

> for column to row do

> VMatrix[row,column] := PotentialIntegral(row,column);

> VMatrix[column,row] := VMatrix[row,column];

> od;

> od;

> print(VMatrix);

1. SCF-Step (Solution of the general eigenvalue problem)

Orthogonalisation of the basis-set

The general eigenvalue problem can be simplified to the special eigenvalue problem by a transformation

to an orthonormal basis set, This is solved by diagonalisation. Othonormalisation means that the

overlap-matrix has to be transformed to the unity-matrix, i.e. a transformation matrix A has to

be found such that A T SA = I. To this purpose we first diagonalise : B T SB = D, where

B = "EigenVectors" is the transformation matrix, and D = "SSpur" is diagonal.

> SSpur := diag(eigenvals(SMatrix));

> evalf(Diagonalisation(SMatrix,EigenValues,EigenVectors));

> print(EigenVectors);

Then we form the square-root of the diagonal matrix = "RootSMatrix"

> RootSMatrix := SSpur:

> for i to Dim do

> RootSMatrix[i,i] := sqrt(EigenValues[i]);

> od:

> print(RootSMatrix);

Multiplying the "inverse RootSMatrix" W with the EigenVectors B from the left yields the

desired transformation matrix: A = BW, as (BW) T SBW = W(B T SB)W =

WDW = I

> TransMatrix := evalm(EigenVectors &* inverse(RootSMatrix));

Fock Operator, transformation to the orthogonal basis-set and diagonalisation

In the first SCF-step, the density matrix is put to zero, i.e. the problem is treated as one-electron problem. The Fock operator in a one electron problem is F = Hcore = T + V = TVMatrix.

Transformation to the orthonormal basis-set yields F' = A T FA = A T TVMatrix A

> TVMatrix := evalm( TMatrix &+ VMatrix); Hcore = T + V

> FockOp := evalm(transpose(TransMatrix) &* TVMatrix &* TransMatrix);

The Fock operator F' is diagonalised and yields as eigenvalues the orbital energies and as

eigenvectors the coefficients c' representing the molecular orbitals in the orthonormal basis-set

> evalf(Diagonalisation(FockOp,EigenValues,OBEigenVec));

> print(OBEigenVec);

Backtransformation of the coefficients and formation of the density matrix

Backtransformation of the coefficients to the original AO-Basis: c = Ac'

> MOEigenVectors := evalm(TransMatrix &* OBEigenVec);

Transformation der Besetzungs-Matrix in der MO-Basis (s. Eingabe) in die ursprngliche

AO-Basis -> Dichtematrix P

> DensityMatrix := evalm(MOEigenVectors &* Occupation &* transpose(MOEigenVectors));

Further SCF-Steps

Two-Electron-Operator for closed shells:

M mn = Element mn of the M-Matrix = Element mn of the (2J - K)-Matrix

= S l S r P lr [( mn / lr )-1/2( ml / nr )]

> MMatrix := array(1..Dim,1..Dim):

> for row to Dim do

> for column to row do

> MMatrix[row,column] := MOperator(row,column);

> MMatrix[column,row] := MMatrix[row,column];

> od;

> od:

> print(MMatrix);

Fock Operator:

F = Hcore + M

> FMatrix := evalm(TVMatrix &+ MMatrix);

F' = A T FA, Transformation of the Fock Operator to the orthonormal basis-set

> FockOp := evalm(transpose(TransMatrix) &* FMatrix &* TransMatrix);

> evalf(Diagonalisation(FockOp,EigenValues,OBEigenVec));

> print(OBEigenVec);

The Fock operator F' is diagonalised and yields as eigenvalues the orbital energies and as

eigenvectors the coefficients c' representing the molecular orbitals in the orthonormal basis-set

Calculate density matrix

Backtransformation of the coefficients to the original AO-Basis: c = Ac'

> MOEigenVectors := evalm(TransMatrix &* OBEigenVec);

Transformation of the Occupation-Matrix from the MO-Basis (see input) to the original

AO-Basis -> Densitymatrix P

> DensityMatrix := evalm(MOEigenVectors &* Occupation &* transpose(MOEigenVectors));

Convergence

Energy = trace of "EMatrix" = trace [P ( H + 1/2 M)]

> EMatrix := evalm(.5 * DensityMatrix &* (TVMatrix &+ FMatrix));

> Etot_old := Etot; The energy is saved in Etot_old , for a comparison in the next step

> Etot := trace(EMatrix); form trace

> Dif := Etot - Etot_old; Energy lowering in the last SCF-Step -> Convergence ?

>

Next SCF-Step

Transform Hcore, i.e. TVMatrix, to a representation in the MO's for the CI calculation

> Hij:= evalm(transpose(MOEigenVectors) &* TVMatrix &* MOEigenVectors);

Save important data for correlation calculation

> save(Dim,Exponent,BasisCoordinate,Occupation,EigenValues,MOEigenVectors,Hij,"SCF.save"):