● 多電子原子の全電子軌道

水素原子の1つの電子の波動関数が計算できたので、 複数の電子を持つ多電子原子の各電子の波動関数を求めよう。 電子間のCoulomb相互作用のため、方程式を代数的に解くことは 不可能であり、まさにコンピュータに頼るところである。

■ 多電子原子系の支配方程式

N 電子系の波動関数 Ψ とそれを支配する方程式の 導出の導出の過程を歴史にそって簡単に紹介しよう。

▲ Hartree方程式

1928年に D.R.Hartree は電子 i が位置 r_1 で感じる他の電子の potential Vh_i(r_1) を次の平均値で扱うことを提案した。

(1)

(2)

(3)

(4)

(5)

(6)

電子 i のスピンが上向きの場合
(7)

電子 i のスピンが下向きの場合
(8)

(9)

複矢印同順ということである。

このHartree-Fock方程式によりPauli原理の要請を満たす軌道波動関数と エネルギーが求まる。

---

▲ Hartree-Fock-Slater方程式

Hartree-Fock 方程式の Vx の演算子の作用の計算はその定義式を見ての 通り複雑である。1951年、Slaterは自由電子ガスを想定してこの項の計算を 劇的に簡単にする近似式を導き、さらにその方法を束縛系であるこの 多電子原子系にも応用することを提案した。この方法はその単純化の割には 良好な計算結果を出してくれる。

Slaterによると Vx^{↑↓}(r_1) は次のように Vsl^{↑↓}(r_1) で近似される。

(10)

(11)

(12)

(13)

(14)

方程式(14)式での 固有値 E_i は軌道エネルギーと呼ばれるエネルギーの次元を持った 量であるが、これは電子 i のエネルギーでもイオン化エネルギーでもない。 従って軌道エネルギーの 総和は N 電子系の全エネルギーではない。 全エネルギーは N 電子系の波動関数 Ψ を Schrodinger方程式に入れ直して求める。その結果、 全エネルギー E は次式で表される。

(15)
(16)

(17)

このことを利用して、イオン化エネルギー I_i は次のように表される。 (18)

基礎方程式
(19)

電子間Coulomb相互作用項
(20)

電子間交換相互作用項
(21)

全エネルギー
(23)

(24)

基礎方程式
(25)
(26)

電子間Coulomb相互作用項
(27)

電子間交換相互作用項
(28)

全エネルギー
(29)

プログラムを作る準備が整った。ソースファイルは atom.ccである。 1電子方程式を解く関数群は水素原子の所で作成したものを利用するので、 多電子版として新たに付け加える部分のみを載せることにしよう。

まずは基本定数、外部変数の定義である。

const int ELE_NUM      = 6;                // number of electrons
const int ORB_NUM      = 3;                // number of orbitals
const double Z         = (double)ELE_NUM;  // atomic number

//                  1S 2S 2P 3S 3P 
const double g[] = { 2, 2, 2, 0, 0 };       // degeneracy of each orbitals
const int    n[] = { 1, 2, 2, 3, 3 };       // quantum number of each orbitals
const int    l[] = { 0, 0, 1, 0, 1 };

Orbital classのpropertyに軌道の縮重度を表す double g を 追加しよう。縮重度が整数でなく実数なのはイオン化エネルギーを計算する ためである。 今回はPotentialはSCFループで調整するものなので、その変化分に合わせて 軌道エネルギーを修正する次のmethodを Orbital classに追加しよう。

//---- modifies energy by first order perturbation Hamiltonian
void Orbital::Perturb( double dU[] )
{
  for( int i=0 ; i<=R_FAR ; i++ ) E += sqr(G[i])*dU[i]*DX;
}

//---- represents density at given grid
double Orbital::Density( int i ){
  return( g*sqr(G[i]) );
}
//---- represents energy of this orbital
double Orbital::Energy( void ){
  return( g*E );
}
//---- represents energy of one electron in this orbital
double Orbital::OneEnergy( void ){
  return( E );
}

Potentialを管理するclass Potential の設計が本プログラムでの 主題とも言えよう。いくつかのmethodが加わる。

Form() methodは密度汎関数法の1電子potentialをその時点での電子 波動関数から構築して U[] に設定するのだが、単純にそのものの値を U[] に代入するのではなく、前回に計算した U[] からの変化が 緩くなるな値を U[] に代入する。これは potentialの急進的な変更は SCFループの収束を不安定にしてしまうからである。 さらに、potentialの変化分を各軌道のobjectに通知して、 各軌道エネルギーを摂動で調整してもらう。

//---- modifies potential
void  Potential::Form( void )
{
  static double dU[R_NUM];
  for( int i=0 ; i<=R_FAR ; i++ ){
    dU[i] = 0.60*( -Z + UCoulomb(i) + UExchange(i) - U[i] );
    U[i] += dU[i];
  }
  extern Orbital O[ORB_NUM];
  for( int e=0 ; e<ORB_NUM ; e++ ) O[e].Perturb( dU );
}

//---- calculates coulomb term at given index
double Potential::UCoulomb( int i )
{
  extern Orbital O[ORB_NUM];
  int j, e;
  double Uh, sigma;

  for( j=0, Uh=0.0 ; j<=R_FAR ; j++ ){
    for( e=0, sigma=0.0 ; e<ORB_NUM ; e++ ) sigma += O[e].Density(j);
    if( j<i ) Uh += sigma * R[j]*DX;
    else      Uh += sigma * R[i]*DX;
  }
  return( Uh );
}

//---- calculates exchange term at given index
double Potential::UExchange( int i )
{
  extern Orbital O[ORB_NUM];
  int e;
  double Uxc, sigma;

  for( e=0, sigma=0.0 ; e<ORB_NUM ; e++ ) sigma += O[e].Density(i);
  const double rs = cbrt(6.0*R[i]*R[i]/sigma);
  const double be = 1.0 + 0.0368*rs*log(1.0+21.0/rs);
  Uxc = -2*be*cbrt(3.0/(32*M_PI*M_PI)*R[i]*sigma);
  return( Uxc );
}

//---- calculates total energy
double Potential::TotalEnergy( void )
{
  extern Orbital O[ORB_NUM];
  int e, i;
  double Etot, sigma;

  for( e=0, Etot=0.0 ; e<ORB_NUM ; e++ ) Etot += O[e].Energy();
  for( i=0 ; i<=R_FAR ; i++ ){
    for( e=0, sigma=0.0 ; e<ORB_NUM ; e++ ) sigma += O[e].Density(i);
    Etot -= 0.25*sigma*(2*UCoulomb(i)+UExchange(i)) * DX;
  }
  return( Etot );
}

最後に main() 関数である。 各軌道を水素原子として初期化して、 SCFループでpotentailを構成して、各軌道を調整して、全エネルギーを 求める。その変化が十分に小さくなるまでSCFループを繰り返す。

//---- definition of fundermental
static Potential U;
static Orbital O[ORB_NUM];

//---- main function
int main( void )
{
  double Etot, Epre;
  int e;

  for( e=0 ; e<ORB_NUM ; e++ ){
    O[e].Assign( g[e], n[e], l[e] );
    O[e].Adjust();
  }
  Etot = U.TotalEnergy();

  do{  // SCF loop
    U.Form();
    for( e=0 ; e<ORB_NUM ; e++ ) O[e].Adjust();
    Epre = Etot;
    Etot = U.TotalEnergy();
  }while( fabs(Etot-Epre)>1e-6 );

  fprintf( stderr, "Total energy is %lf [a.u]\n\n", Etot );
  for( e=0 ; e<ORB_NUM ; e++ ) O[e].Output();

  return(0);
}

このプログラムが算出する種々の原子の電子のイオン化エネルギーと 全エネルギーを載せる。ただしイオン化エネルギーを求めるには軌道の 占有数を 0.5 だけ減らして再計算する若干の改造が必要である。

より多電子の系を計算する場合には、より慎重に解かなくてはならない。 Form() method でのpotentialの更新をより緩やかに 更新するとうまくいく。セルの数や空間の範囲にも改良が必要である。 単純に広げるだけでは指数関数の効果により波動関数の値が 0 もしくは ∞ となり計算ができなくなるので注意が必要である。

多電子系の電子の準位で興味深い現象は 20番目の原子 Caの最外殻電子2つは3D軌道を空にして 4S軌道に入ることである。 この原因は、このような占有数の状態と、3D軌道を先に 詰めた占有数の状態の全エネルギーを計算して比較すれば明らかになる。 つまり、電子間の相互作用により、4S軌道に先に入る方が 全エネルギーがわずかに下がるのである。