MP2-program for small molecules utilizing basis-sets consisting of s-functions (spherical Gaussians) only. Reads data from file, which were generated by the ab initio SCF-program.

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

All properties in atomic units!

Input from file

> restart:

> with(linalg):

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

> Digits:=10:

> read("SCF.save");

Parameter und Auxiliary Functions

> X:=1:Y:=2:Z:=3:rootPi:=evalf(sqrt(Pi)):

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:

Occupied orbitals - Function to calculate the number of occupied MOs from the occupation matrix

> OccupiedOrbitals := proc(OccupationMatrix)

> local counter,i;

> counter:=0;

> for i to Dim do

> counter := counter + OccupationMatrix[i,i];

> od;

> counter/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 two-lectron-integrals

> AuxInt := proc(X)

> local rootX;

> if X=0 then 1 else

> rootX:=sqrt(X);

> evalf(0.5*rootPi/rootX*erf(rootX)) fi;

> end:

Integral-functions and -transformation

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

> Q*aux^0.75;

> 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/rootPi) * 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:

MP2-Energy

> MP2Energy := proc(occupied)

> local Sum,iajb,ibja,i,j,a,b;

> Sum:=0;

> for i from 1 to occupied do

> for j from 1 to occupied do

> for a from occupied+1 to Dim do

> for b from occupied+1 to Dim do

> iajb:=IJKLIntegral(i,a,j,b);

> ibja:=IJKLIntegral(i,b,j,a);

> Sum := Sum + iajb*(2*iajb-ibja)/(EigenValues[i]+EigenValues[j]-EigenValues[a]-EigenValues[b]);

> od:od:od:od:

> Sum;

> end:

> MP2Energy(OccupiedOrbitals(Occupation));

>