DFT (density functional theory) program for atoms utilizing

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

copyright: P. Vogt, H. Huber, Jan. 2000

All properties in atomic units!

In this program we kept the names and structure as close as possible to the corresponding MAPLE-sheet for SCF-calculations, to show the analogies and differences

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:

> NumberofNuclei := 1: Example: helium atom i.e. 1 nucleus in the origin

> Dim := 3: 3 basis functions on the nucleus

define nuclear data:

atomic number

X Y Z

> AtomicNumber := array([2]);

> NuclearCoordinate := array([[0,0,0]]);

define basis-set data:

Exponent

X Y Z

> Exponent := array([3.6,1.2,0.4]);

> BasisCoordinate := array([[0,0,0],
[0,0,0],
[0,0,0]]);

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

> Occupation := array([[2,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- and Vx-Operator)

SCF like functions

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

as in SCF, but without exchange integral!

> 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);

> od;

> od;

> Sum;

> end:

Exchange functional

rho - calculates the density (called by ValueOfVx)

> rho := proc(r)

> local row, column, sum, mu, nu, factor, prefactor;

prefactor = (2/Pi)**(3/2)

> prefactor := .50794908747392775826;

> sum := 0;

> for row from 1 to Dim do

> for column from 1 to row do

> mu := Exponent[row];

> nu := Exponent[column];

> factor := 2*(mu * nu)**0.75;

> if row=column then

> factor:=0.5*factor;

> fi;

> sum := sum + factor * DensityMatrix[row,column] * exp(-(mu + nu) * r**2);

> od:od:

> prefactor*sum:

> end:

ValueOfVx - calculates the value of the exchange functional

(uses rho and is called by the numerical integration trapzd, which in turn is called by DensityIntegral)

> ValueOfVx := proc(i,j,x)

> local rhotmp,result;

> rhotmp := rho(x);

> if rhotmp = 0 then

> result := 0;

> else

> result := x**2 * exp(-(Exponent[i]+Exponent[j]) * x**2+ln(rhotmp) / 3);

> fi;

> result;

> end:

DensityIntegral - numerical integration over exchange density (for a spherical density only!)

This is the trapezoidal integration routine used in DensityIntegral

(this routine from Numerical Recipe uses ValueofX (see above))

> trapzd := proc(i,j,A,B,st::evaln,K,it::evaln)

> local result, tnm, loop, del, x, sum;

> if eval(K) = 1 then

> result := 0.5 * (B-A) * (ValueOfVx(i,j,A) + ValueOfVx(i,j,B));

> it := 1;

> else

> tnm := eval(it);

> del := (B-A) / tnm;

> x := A + 0.5 * del;

> sum := 0;

> for loop from 1 to eval(it) do

> sum := sum + ValueOfVx(i,j,x);

> x := x + del;

> od;

> result := 0.5 * (eval(st) + (B-A) * sum / tnm);

> it := eval(it) * 2;

> fi;

> result;

> end:

DensityIntegral - calculates the integral used in VxOperator

(this is procedure QSIMP from Numerical Recipes; it calles trapzd)

> DensityIntegral := proc(i,j)

> local A,B,EPS, s, st, ost, os, K, it, loop;

> A:=0;B:=10;

> EPS:=1e-6;

> K:=1;it := 0;loop:=true;

> st:=-1e30;s:=-1e30;

> os:=s;ost:=st;

> while (loop) do

> os := s;

> ost := st;

> st := evalf(trapzd(i,j,A,B,ost,K,it));

> K := K + 1;

> s := (4 * st - ost) / 3;

> loop := evalb(abs(s-os) > EPS * abs(os));

> od;

> s;

> end:

VxOperator - Exchange Operator

> VxOperator := proc(row,column)

> local prefactor;

prefactor = -(2/Pi)**(3/2) * 4*Pi * (3/Pi)**(1/3)

> prefactor := -6.2857027940461438505;

> evalf(prefactor * sqrt(sqrt((Exponent[row] * Exponent[column])**3))*DensityIntegral(row,column));

> end:

Construction of integral matrices

The same as in the SCF program

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 Dim do

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

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

> od;

> od;

> print(VMatrix);

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

The same as in the SCF program

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 DFT-Steps

Two-Electron-Coulomb-Operator Matrix

M mn = Element mn of the M-Matrix = Element mn of the 2J-Matrix

= S l S r P lr ( mn / lr )

> 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);

Two-Electron-Exchange-Operator Matrix

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

> for row to Dim do

> for column to row do

> VxMatrix[row,column] := VxOperator(row,column);

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

> od;

> od;

> print(VxMatrix);

Fock Operator:

F = Hcore + M + Vx

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

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 &+ VxMatrix &+ FMatrix &+ (-.5 * VxMatrix) ));

> 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 DFT-Step

>