球面ガウシアン(s関数)のみの基底関数を使った第一原理計算
(Ab initio)

第一原理計算を理解するために基底関数としてs関数(球面ガウシアン)だけを使って

SCF (Self Consistency Field) ループ計算を行う Maple プログラムを紹介します。

[Plot]

単位は全て原子単位 (Å: オングストローム) とします。

最も簡単な水素分子の計算例です。

原子の数および座標を設定すれば他の分子でも計算できます。

> restart:

ちなみにガウシアン s 関数は以下のものです。

> g:=x->exp(-x^2);

g := proc (x) options operator, arrow; exp(-x^2) end proc

> with(plots):with(plottools):
plot3d(g(sqrt(x^2+y^2)),x=-2..2,y=-2..2);

Warning, the name changecoords has been redefined

Warning, the assigned name arrow now has a global binding

[Plot]

指数のx^2に係数が入り異なる基底関数を表現します。

exp(-c*x^2)

また基底関数は確立密度関数なので積分すると1になるように正規化されている必要がありますのでこの関数全体に係数がかかります。

基本データ入力(原子核データ、基底関数セット、占有行列)

線形代数パッケージ(linalg)を呼び出します。固有値など線形計算を行う関数群があります。

> with(linalg):

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

Mapleの計算精度(桁数)を設定します。これはDigitsという変数に桁数を与えることで設定できます。

> Digits:=10:                   #     精度指定

NumberofNucleiという変数に原子核の数を与えます。ここでは水素分子なので2個です。

> NumberofNuclei := 2:          #   例 : 水素分子(2原子(核)

次元の指定を行います。全体の基底関数の数です。

> Dim := 4:                     #    核に2個の基底関数すなわち4基底関数(4次元)

原子核データの定義

原子番号(atomic number)

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

> print(%);

vector([1, 1])

1は水素です。

次に各原子の座標を定義します。原子核の座標(X Y Z)を与えます。

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

二つの原子核の座標です。 一個は原点で間隔 1.4 Å(x軸方向)あけてもう一個の原子があります。

> print(%);

matrix([[0, 0, 0], [1.4, 0, 0]])

水素分子 H2 をプロットしてみます。  原子の半径を1Åでプロットしています。

> display3d({sphere([0,0,0],1),sphere([1.4,0,0],1)},scaling=constrained,lightmodel=light4,style=patchnogrid);

[Plot]

基底関数セットの定義:

ガウス関数の指数(Exponent)で基底を表現します。

> exp(-c*x^2);

exp(-c*x^2)

各基底関数の指数を与え基底関数のセットを定義します。

> Exponent := array([1.2,0.3,1.2,0.3]):
# 4つの(ガウス関数の)指数 (各原子核に二つの基底関数)

> print(%);

vector([1.2, .3, 1.2, .3])

次に各基底関数の座標を与えます。 X Y Z

> BasisCoordinate := array([[0,0,0],[0,0,0],[1.4,0,0],[1.4,0,0]]):
# 各原子核毎に 2 個の基底関数 (2 basis functions per nucleus)

> print(%);

matrix([[0, 0, 0], [0, 0, 0], [1.4, 0, 0], [1.4, 0, 0]])

分子軌道表現(MO-基底)における占有行列(Occupation-Matrix) を定義します。

> Occupation := array([[2,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]):

最下位のMO(軌道)が2電子で占有されています。

> print(%);

matrix([[2, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]])

パラメータとその他の関数

> X:=1:
Y:=2:
Z:=3:
# プログラムが見やすいようにx、yおよびz座標として配列に保存されます。

計算に使われる関数群

sqrDistance - 基底関数間の距離の平方根

二つの基底関数間の距離の二乗を計算します。

> 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 proc:

ProductGaussian - ガウシアン関数の積

二つのガウシアン関数の積を計算

二つのガウス関数の積は新しい指数cnewでまたガウス関数となります。

これは元の指数の和です。

そして新しい座標はRxnew, Rynew, Rznew,で古い指数に重み付けられたものです。

> 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 proc:

AuxInt - 積分で使われる式

ポテンシャルおよび2電子積分の計算で使われる積分関数です。erfは誤差関数(ガウシアン関数の積分)です。

> AuxInt := proc(X)
 if X=0 then
   1
 else
   evalf(1/2*sqrt(Pi/X)*erf(sqrt(X)))
 end if;
end proc:

Diagonalisation - 対角化

Jacobi法による対角化を行う関数です。固有値とその固有ベクトルを計算します。

a: 対角化を行う対象の行列

d: 固有値

v: 固有ベクトル

最初に対称化を行います。また、終了時に固有値と対応する固有ベクトルをソートします。

> 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;
   end do:
   v[ip,ip]:=1;
 end do:
 for ip from 1 to Dim do
   b[ip]:=a[ip,ip];
   d[ip]:=b[ip];
   z[ip]:=0;
 end do:
 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]);
     end do:
   end do:
   if(sm<10**(-Digits-2)) then break end if;
   if(i<4)then
     tresh:=0.2*sm/Dim**2;
   else
     tresh:=0;
   end if;
   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;
           end if;
         end if;
         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);
         end do;
         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);
         end do:
         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):
         end do:
         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):
         end do:
       end if:
     end do:
   end do:
   for ip from 1 to Dim do
     b[ip]:=b[ip]+z[ip]:
     d[ip]:=b[ip]:
     z[ip]:=0:
   end do:
 end do:
 Sort(d,v);
end proc:

Sort - ソート

固有値をソート(小さい順に)します。

固有値の大きさに従って固有値と対応する固有ベクトルを並べ替えます。

> 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];
     end if;
   end do;
   if MaxIndex < Dim+1 then
     Eigenvect:=swapcol(Eigenvect,column1,MaxIndex);
     Eigenvalues[MaxIndex]:=Eigenvalues[column1];
     Eigenvalues[column1]:=MaxValue;
   end if;
 end do;
end proc:

Symmetrisation - 対称化

行列対称化 行列 A を対称化します。

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

積分関数 - S, T, V, 2e-integral, M

積分を行う関数群です。

Overlapintegral - S オーバラップ積分

オーバーラップ積分 オーバーラップ積分 <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 proc:

KineticIntegral - T 運動積分

運動積分 <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 proc:

PotentialIntegral - V ポテンシャル積分

ポテンシャル積分  <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);
 end do;
 evalf(-Overlapintegral(func1,func2) * 2/sqrt(Pi) * sqrt(c) * V);
end proc:

TwoElectronIntegral - 2電子積分 (2-e integral)

2電子積分 2-e積分  <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 proc:

MOperator - M 閉殻2電子オペレータ

閉殻の2電子オペレータです。

> 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));
   end do;
 end do;
 Sum;
end proc:

積分行列の構成

積分行列の構成

最初に積分行列 S を構成します。

基底関数間のオーバーラップ積分 (内積) で行列要素を決定します。

オーバーラップ行列 S

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

> for irow to Dim do
 for column to irow do
   SMatrix[irow,column] := Overlapintegral(irow,column);  # i行j列のオーバーラップ積分
   SMatrix[column,irow] := SMatrix[irow,column];        # 対称
 end do:
end do:

> print (SMatrix):

matrix([[1.000000000, .7155417528, .3085103151, .4470363683], [.7155417528, 1.000000000, .4470363683, .7452764914], [.3085103151, .4470363683, 1.000000000, .7155417528], [.4470363683, .7452764914, .71...

運動エネルギー 行列 T の構成

運動積分により運動エネルギー行列 T を構成します。

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

> for irow to Dim do
 for column to irow do
   TMatrix[irow,column] := KineticIntegral(irow,column);
   TMatrix[column,irow] := TMatrix[irow,column];
 end do:
end do:

> print(TMatrix);

matrix([[1.800000000, .5151900621, .1199488105, .2209289495], [.5151900621, .4500000000, .2209289495, .2696410346], [.1199488105, .2209289495, 1.800000000, .5151900621], [.2209289495, .2696410346, .51...

ポテンシャル・エネルギー行列 V

ポテンシャル積分によりポテンシャル行列 V を構成します。

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

> for irow to Dim do
 for column to irow do
   VMatrix[irow,column] := PotentialIntegral(irow,column);
   VMatrix[column,irow] := VMatrix[irow,column];
 end do;
end do;

> print(VMatrix);

matrix([[-2.460820054, -1.492136295, -.7711679084, -.9726350940], [-1.492136295, -1.498951325, -.9726350940, -1.185642164], [-.7711679084, -.9726350940, -2.460820054, -1.492136295], [-.9726350940, -1....

SCFループ

SCFループと呼ばれる繰り返しによりセルフ・コンシステント(自己整合)なポテンシャル場を決定します。

基底セットの対角化

一般化固有値問題は直交正規基底セットに変換することにより特殊固有値問題に単純化できます。

これは対角化によって解くことができます。

直交正規化はオーバーラップ行列が単位行列に変換されなければなりません。

すなわち変換行列 A は A^t*S*A = I を満たさなければなりません。

この目的のために最初にBTSB = Dと対角化します。

ここで B = "固有ベクトル" は変換行列で D = "SSpur" は対角行列です。

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

> print(%);

matrix([[.1270050008, 0, 0, 0], [0, .3439821383, 0, 0], [0, 0, .8192081927, 0], [0, 0, 0, 2.709804668]])

S 行列の対角化を行います。

固有値および固有ベクトルを解きます。

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

上が固有値で下が同時に得られた固有ベクトル(座標変換行列)です。

> print(EigenVectors);

matrix([[.3037351969, -.5442000863, -.6385490822, .4514933730], [-.6385490820, .4514933731, -.3037351971, .5442000864], [-.3037351973, -.5442000865, .6385490817, .4514933732], [.6385490821, .451493373...

対角行列の平方根の行列”RootS行列"を形成します。

> RootSMatrix := SSpur:

> for i to Dim do
 RootSMatrix[i,i] := sqrt(EigenValues[i]);
end do:

> print(RootSMatrix);

matrix([[.3563776093, 0, 0, 0], [0, .5864999047, 0, 0], [0, 0, .9051012058, 0], [0, 0, 0, 1.646148434]])

この行列を平方すると、先ほどの対角行列の行列 S になります。

固有ベクトル B を左から RoosS 行列の逆行列 W にかけることにより目的の変換行列 A が求まります。

A = B*W

(B*W)^T*S*B*W = W*B^T*S*B*W = W*D*W = I

となります。

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

> print(%);

matrix([[.8522847367, -.9278775357, -.7055002006, .2742725768], [-1.791776659, .7698097979, -.3355814744, .3305899243], [-.8522847378, -.9278775361, .7055002001, .2742725769], [1.791776659, .769809797...

Fock オペレータ:

直交基底セットを変換し対角化を行います。

最初の SCF ステップにおいて密度行列はゼロにされます。すなわり問題は1電子問題として扱われます。

1電子問題の Fock オペレータは F = Hcors = T + V = TV 行列です。

直交基底への変換は F^`'` = A^T*F*A = A^T*TVMatrix*A で得られ、

直交基底への変換は F^`'` = A^T*TV*A です。

Hcore = T + V となります。

> TVMatrix := evalm( TMatrix &+ VMatrix):

> print(%);

matrix([[-.660820054, -.9769462329, -.6512190979, -.7517061445], [-.9769462329, -1.048951325, -.7517061445, -.9160011294], [-.6512190979, -.7517061445, -.660820054, -.9769462329], [-.7517061445, -.916...

> Etot:=0;

Etot := 0

Fock オペレータは TV 行列に両側から変換行列を求めることができます。

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

> print(%);

matrix([[.5082484784, 0.85e-9, -.5889457768, 0.85e-9], [0.2e-9, .3509235191, 0.16e-9, -0.177128680e-2], [-.5889457779, 0.6e-9, -.2528066151, 0.50e-9], [0.1e-8, -0.17712863e-2, 0.2e-9, -1.253855330]])

Fock オペレータF'は対角化され、固有値として軌道エネルギーを導きます。

そして固有ベクトルとして係数 c' は直交基底関数セットでの分子軌道を表します。

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

> print(OBEigenVec);

matrix([[-0.7994702426e-9, .4781771539, -0.2099661022e-9, .8782634055], [0.1103755396e-2, -0.6305032602e-9, .9999993909, 0.5833564796e-9], [-0.8204018660e-9, .8782634055, 0.8336017124e-9, -.4781771539...

係数の逆変換と密度行列の構成

係数を元のAO基底に逆変換します。

c = Ac'

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

> print(%);

matrix([[.2732482598, -.2120719182, -.9281797011, 1.085884573], [.3314394063, -1.151515592, .7694444387, -1.413184476], [.2732482601, .2120719190, -.9281796999, -1.085884575], [.3314394030, 1.15151559...

MO基底に関する占有(Besetzungs)行列を元のAO基底に変換すると密度行列 P が導かれます。 

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

> print(%);

matrix([[.1493292230, .1811304820, .1493292231, .1811304802], [.1811304820, .2197041601, .1811304822, .2197041579], [.1493292231, .1811304822, .1493292233, .1811304804], [.1811304802, .2197041579, .18...

さらなるSCFステップ

さらなる SCF ステップはここから繰り返されます。(ループの先頭)

求まった密度行列 P を元に閉じた殻 (shell) に対する二つの電子オペレータの計算を行い M 行列を計算します。

ここで電子間の相互作用が考慮されます。

Mmn = M 行列の要素 mn = (2J-K)-行列の要素 mn

   = SlSr Plr[(mn/lr)-1/2(ml/nr)]

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

> for irow to Dim do
 for column to irow do
   MMatrix[irow,column] := MOperator(irow,column);
   MMatrix[column,irow] := MMatrix[irow,column];
 end do;
end do:

> print(MMatrix);

matrix([[1.093540332, .5807486174, .1319176057, .2548980657], [.5807486174, .7281860308, .2548980655, .4815359172], [.1319176057, .2548980655, 1.093540331, .5807486178], [.2548980657, .4815359172, .58...

Fock オペレータは

F = Hcore + M です。

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

> print(%);

matrix([[.432720278, -.3961976155, -.5193014922, -.4968080788], [-.3961976155, -.3207652942, -.4968080790, -.4344652122], [-.5193014922, -.4968080790, .432720277, -.3961976151], [-.4968080788, -.43446...

F^`'` = A^T*F*A

Fock オペレータの直交基底セットへの変換

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

> print(%);

matrix([[1.498565174, -0.73e-9, -.8113324122, 0.85e-9], [-0.6e-9, 1.507262042, 0.68e-9, -.1695717342], [-.8113324121, 0.92e-9, 1.068588456, 0.16e-9], [0.10e-8, -.1695717342, 0.1e-9, -.5019856068]])

Fock オペレータを対角化します。

次は固有値です。

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

次はその固有ベクトルです。

> print(OBEigenVec);

matrix([[-0.6044510743e-9, .6098598854, 0.1501947462e-8, .7925092557], [0.8351020165e-1, -0.9959021925e-10, .9965069223, -0.1748228330e-8], [-0.4333489726e-9, .7925092557, -0.9543928666e-9, -.60985988...

Fock オペレータ F' は対角化され、固有値として軌道エネルギーを導き、固有ベクトルとして直交基底セットに関する分子軌道を表す係数 c' を導きます。

密度行列の計算

係数の元の AO 基底への変換: c = Ac'

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

> print(%);

matrix([[.1958272811, -0.3934116670e-1, -.9475409436, 1.105699815], [.3937221207, -1.358684132, .7395131589, -1.215341908], [.1958272816, 0.3934116630e-1, -.9475409480, -1.105699813], [.3937221183, 1....

占有行列をMO基底 (入力を参照) から元の AO 基底に変換すると密度行列 P が求められます。

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

> print(%);

matrix([[0.7669664805e-1, .1542030648, 0.7669664824e-1, .1542030639], [.1542030648, .3100342167, .1542030652, .3100342148], [0.7669664824e-1, .1542030652, 0.7669664844e-1, .1542030643], [.1542030639, ...

次に収束性の判定を行います。

エネルギーは行列 E のトレースで求められます。まずエネルギー行列 E を求めます。

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

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

> print(%);

matrix([[-.2557686012, -.3102664674, -.2557686012, -.3102664674], [-.5142376260, -.6238087505, -.5142376255, -.6238087505], [-.2557686019, -.3102664682, -.2557686018, -.3102664682], [-.5142376225, -.6...

前回までの全エネルギーを保存しておきます。

> Etot_old := evalf(Etot);   # 次のループのために前の全エネルギーを Etot_old に保存します。

Etot_old := 0.

トレースで全エネルギーが求まります。

> Etot := trace(EMatrix);  # トレースすなわち全エネルギーを計算

Etot := -1.759154701

> Dif := Etot - Etot_old;  #

Dif := -1.759154701

差をとり最後の SCF ステップにおけるエネルギー低下を求めます。

      

最後の SCF ステップでのエネルギー低下が小さければ収束と判定します。

全エネルギーの前回のものと差をとり、充分小さければ (0 に近ければ) 収束したものと判定します。

収束していなければ
SCF ループの先頭へ戻り収束するまで繰り返し計算を行います。

これが SCF ループです。

SCFループ

これは次のように for ループで繰り返すことができます。nn は繰り返しの回数。

> nn:=3;

nn := 3

> for kk to nn do
 print("kk=",kk);
 for irow to Dim do
   for column to irow do
     MMatrix[irow,column] := MOperator(irow,column);
     MMatrix[column,irow] := MMatrix[irow,column];
   end do:
 end do:
 FMatrix := evalm(TVMatrix &+ MMatrix);
 FockOp := evalm(transpose(TransMatrix) &* FMatrix &* TransMatrix);
 evalf(Diagonalisation(FockOp,EigenValues,OBEigenVec));
 MOEigenVectors := evalm(TransMatrix &* OBEigenVec);
 DensityMatrix := evalm(MOEigenVectors &* Occupation &* transpose(MOEigenVectors));
 EMatrix := evalm(.5 * DensityMatrix &* (TVMatrix &+ FMatrix)):
 print("E=",%);
 Etot_old := evalf(Etot);  # 次のループのために前の全エネルギーを Etot_old に保存します。
 Etot := trace(EMatrix);  # トレースすなわち全エネルギーを計算
 Dif := Etot - Etot_old;
 print("Etot=",Etot,"Dif=",Dif);
end do:

kk=, 1

E=, matrix([[-.2442245148, -.2968078420, -.2442245146, -.2968078420], [-.5333862150, -.6482281750, -.5333862150, -.6482281750], [-.2442245152, -.2968078426, -.2442245152, -.2968078426], [-.533386212...

Etot=, -1.784905377, Dif=, -0.25750676e-1

kk=, 2

E=, matrix([[-.2424313434, -.2947088574, -.2424313432, -.2947088574], [-.5360622745, -.6516579015, -.5360622745, -.6516579015], [-.2424313436, -.2947088578, -.2424313436, -.2947088578], [-.536062272...

Etot=, -1.788178487, Dif=, -0.3273110e-2

kk=, 3

E=, matrix([[-.2421670548, -.2943993046, -.2421670548, -.2943993046], [-.5364505855, -.6521559240, -.5364505855, -.6521559240], [-.2421670552, -.2943993050, -.2421670550, -.2943993050], [-.536450583...

Etot=, -1.788645956, Dif=, -0.467469e-3

Eはエネルギー行列

Etotは全エネルギーです。

配置相関(CI)計算

Hcore を変換すなわり配置相関 (Cofiguration Interaction: CI 計算) の為に、TV 行列を分子軌道 MO に関する表現に変換します。

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

> print(%);

matrix([[-1.238922025, -0.1e-8, .1540934355, 0.], [-0.94e-9, -.5431346098, 0., .2039882007], [.1540934356, -0.3e-9, .3359902135, -0.5e-9], [0.102e-8, .2039881993, 0.129e-8, .7985764710]])

データの保存

相関計算を行うためにデータ(固有値および固有ベクトル)を保存します。

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

フローチャートの作成

> hako:=proc(x,y,l,h)
 return(rectangle([x-l,y-h],[x+l,y+h])):
end proc:

> p:=textplot([0,0,'SCFループのフローチャート']):

> yy:=-2:

>

> yy:=yy-3:

> p:=p,textplot([0,yy,' ']):

> yy:=yy-1:

> p:=p,textplot([0,yy,'オーバーラップ積分行列の生成']):

> p:=p,hako(0,yy+0.2,10,1.5):

> p:=p,line([0,yy-1.2],[0,yy-2.2]):

>

> yy:=yy-3:

> p:=p,textplot([0,yy,' ']):

> yy:=yy-1:

> p:=p,textplot([0,yy,'運動およびポテンシャル・エネルギー行列の生成']):

> p:=p,hako(0,yy+0.2,10,1.5):

> p:=p,line([0,yy-1.2],[0,yy-2.2]):

>

> yy:=yy-3:

> p:=p,textplot([0,yy,' ']):

> yy:=yy-1:

> p:=p,textplot([0,yy,'固有値問題を解き対角化']):

> p:=p,hako(0,yy+0.2,10,1.5):

> p:=p,line([0,yy-1.2],[0,yy-2.2]):

>

> yy:=yy-3:

> p:=p,textplot([0,yy,' ']):

> yy:=yy-1:

> p:=p,textplot([0,yy,'対角化し逆行列を求める']):

> p:=p,hako(0,yy+0.2,10,1.5):

> p:=p,line([0,yy-1.2],[0,yy-2.2]):

>

> yy:=yy-3:

> p:=p,textplot([0,yy,' ']):

> yy:=yy-1:

> p:=p,textplot([0,yy,'変換行列を求める']):

> p:=p,hako(0,yy+0.2,10,1.5):

> p:=p,line([0,yy-1.2],[0,yy-2.2]):

>

> yy:=yy-3:

> p:=p,textplot([0,yy,' ']):

> yy:=yy-1:

> p:=p,textplot([0,yy,'密度行列をゼロとしFockオペレータを計算']):

> p:=p,hako(0,yy+0.2,10,1.5):

> p:=p,line([0,yy-1.2],[0,yy-2.2]):

> yy:=yy-3:

> p:=p,textplot([0,yy,' ']):

> yy:=yy-1:

> p:=p,textplot([0,yy-0.1,`Fockオペレータを直交基底に座標変換`]):

> p:=p,hako(0,yy+0.2,10,1.5):

> p:=p,line([0,yy-1.2],[0,yy-2.2]):

> yy:=yy-3:

> p:=p,textplot([0,yy,' ']):

> yy:=yy-1:

> p:=p,textplot([0,yy-0.1,`Fockオペレータを対角化`]):

> p:=p,hako(0,yy+0.2,10,1.5):

> p:=p,line([0,yy-1.2],[0,yy-2.2]):

>

> yy:=yy-3:

> p:=p,textplot([0,yy,' ']):

> yy:=yy-1:

> p:=p,textplot([0,yy,'密度行列を計算']):

> p:=p,hako(0,yy+0.2,10,1.5):

> p:=p,line([0,yy-1.2],[0,yy-2.2]):

> p:=p,line([0,yy-2],[15,yy-2]):                 yyy:=yy-2:

> yy:=yy-3:

> p:=p,textplot([0,yy,' ']):

> yy:=yy-1:

> p:=p,textplot([0,yy,'閉じた殻に対する2電子オペレータ計算']):

> p:=p,hako(0,yy+0.2,10,1.5):

> p:=p,line([0,yy-1.2],[0,yy-2.2]):

> yy:=yy-3:

> p:=p,textplot([0,yy,' ']):

> yy:=yy-1:

> p:=p,textplot([0,yy,'Fockオペレータを計算']):

> p:=p,hako(0,yy+0.2,10,1.5):

> p:=p,line([0,yy-1.2],[0,yy-2.2]):

> yy:=yy-3:

> p:=p,textplot([0,yy,' ']):

> yy:=yy-1:

> p:=p,textplot([0,yy-0.1,`Fockオペレータを直交基底に座標変換`]):

> p:=p,hako(0,yy+0.2,10,1.5):

> p:=p,line([0,yy-1.2],[0,yy-2.2]):

> yy:=yy-3:

> p:=p,textplot([0,yy,' ']):

> yy:=yy-1:

> p:=p,textplot([0,yy-0.1,`Fockオペレータを対角化`]):

> p:=p,hako(0,yy+0.2,10,1.5):

> p:=p,line([0,yy-1.2],[0,yy-2.2]):

> yy:=yy-3:

> p:=p,textplot([0,yy,' ']):

> yy:=yy-1:

> p:=p,textplot([0,yy-0.1,`密度行列の計算`]):

> p:=p,hako(0,yy+0.2,10,1.5):

> p:=p,line([0,yy-1.2],[0,yy-2.2]):

> yy:=yy-3:

> p:=p,textplot([0,yy,' ']):

> yy:=yy-1:

> p:=p,textplot([0,yy-0.1,`エネルギー行列の計算`]):

> p:=p,hako(0,yy+0.2,10,1.5):

> p:=p,line([0,yy-1.2],[0,yy-2.2]):

> yy:=yy-3:

> p:=p,line([10,yy-1],[15,yy-1]):

> p:=p,line([15,yy-1],[15,yyy]):

> p:=p,textplot([0,yy,' ']):

> yy:=yy-1:

> p:=p,textplot([0,yy-0.1,`収束の判定`]):

> p:=p,hako(0,yy+0.2,10,1.5):

> p:=p,line([0,yy-1.2],[0,yy-2.2]):

> yy:=yy-3:

> p:=p,textplot([0,yy,' ']):

> yy:=yy-1:

> p:=p,textplot([0,yy-0.1,`SCFループの収束`]):

> p:=p,hako(0,yy+0.2,10,1.5):

>

> display([p],axes=none);

[Plot]