CI 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 from file (basis-set data, occupation matrix, MO-representation of Hcore)

> restart:

> with(linalg):

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

> Digits:=10:

> read("SCF.save"); Example: hydrogen molecule i.e. 2 atoms (nuclei)

Parameters and auxiliary functions

> X:=1:Y:=2:Z:=3:

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

> 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, 2e-integral and its MO-representation)

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:

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:

Transformation of 2-e-integral s from the AO-Basis ( ) to the MO-Basis (ijkl)

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

> local Sum,mu,nu,lambda,sigma;

> Sum:=0;

> for mu from 1 to Dim do

> for nu from 1 to Dim do

> for lambda from 1 to Dim do

> for sigma from 1 to Dim do

> Sum:=Sum+MOEigenVectors[mu,i]*MOEigenVectors[nu,j]*MOEigenVectors[lambda,k]*MOEigenVectors[sigma,l]*TwoElectronIntegral(mu,nu,lambda,sigma);

> od:od:od:od:

> Sum;

> end:

Construct CI-Matrix

> Element11:=2*Hij[1,1]+ijklIntegral(1,1,1,1);

> Element22:=2*Hij[2,2]+ijklIntegral(2,2,2,2);

> Element33:=2*Hij[3,3]+ijklIntegral(3,3,3,3);

> Element44:=Hij[2,2]+Hij[3,3]+ijklIntegral(2,2,3,3);

> Element55:=Element44;

> Element12:=ijklIntegral(2,1,2,1);

> Element13:=ijklIntegral(3,1,3,1);

> Element14:=ijklIntegral(3,1,2,1);

> Element15:=Element14;

> Element23:=ijklIntegral(3,2,3,2);

> Element24:=Hij[2,3]+ijklIntegral(3,2,2,2);

> Element25:=Element24;

> Element34:=Hij[2,3]+ijklIntegral(3,3,3,2);

> Element35:=Element34;

> Element45:=Element23;

> CIMatrix:=array([[Element11,Element12,Element13,Element14,Element15],
[Element12,Element22,Element23,Element24,Element25], [Element13,Element23,Element33,Element34,Element35], [Element14,Element24,Element34,Element44,Element45], [Element15,Element25,Element35,Element45,Element55]]);

> evalf(Diagonalisation(CIMatrix,CIEigenValues,CIEigenVec));

> print(CIEigenVec);

>

>