第７回、乱数発生 --- 可視化 ---

目次

• 小テスト
• 乱数発生 --- 可視化 ---
• 今日の課題
• レポート提出用form

乱数発生 --- 可視化 ---

• ２x２行列、ベクトルの入力
• 最小２乗法をかける解析データの入力

水素原子の波動関数の表示

例題１

水素原子1s軌道の波動関数Ψ、Ψ²を電子と原子核との距離に ついて計算し出力する。

解説

水素類似原子の波動関数はアトキンス物理化学（上）P.513を参照せよ。 s軌道は電子と原子核との距離のみの関数なので、角度部分はすでに積分した 動径波動関数を計算する。適当な間隔(下のプログラムでは0.01 A 間隔)で Ψ、Ψ²を計算し表として出力する。

解答

      implicit real*8(a-h,o-z)
      real*8 x, y, y2
      integer i
c
      do i = 0 , 500
         x = i / 100.0d0
         y = f1s(x)
         y2 = y**2
         write(6,'(3e20.10)') x , y, y2 
      end do
      end
c==================================================================
      real*8 function f1s(r)
c
c     1s function for hydrogen atom
c
      implicit real*8(a-h,o-z)
      parameter(r0 = 0.529d0)
c
      f1s = ( 1/r0 )**(1.5d0) * exp(-r/r0) 
c
      return
      end

さて、上のプログラムを"Intro"を使ってコンパイル、実行してみましょう。 ただし計算の出力結果は１００行にもわたるので、出力ファイルを指定して おきましょう。Introの設定を変更していなければ出力ファイル名を指定 できないのですが、これをできるようにする方法を 山本君に教えてもらいました。もちろん私が口頭で以前説明したshellから 利用する方法でもよいのですが、こちらの方がわかりやすいと思います。ここでは その出力ファイル名をprog7-1.outと しておきます。
次にグラフを作成するソフトを起動します。現在みなさんが使っている NextStepの中にグラフ作成ツールはいくつかありますが、ここでは"GNUPLOT"を 使用します。
"GNUPLOT"はファイルビューアを使って起動します。場所は

電子雲の表示

例題２

例題１と同じ1s波動関数の可視化を行なう。今度はある平面内に電子が どのように分布しているか、を表示する。アトキンス物理化学（上）P.522 の 図13.8をある平面（例えばxy平面）で輪切りにしたものを想像して欲しい。

解説

電子は原子核の回りをある軌道にそって運動しているのではなく、任意の 領域にみいだす確率のみ意味を持つ。これを計算機内でシミュレーションする ために以下のようなテクニックを使用する。

1. 一様乱数を１つ発生させる
2. それに対する波動関数の値とまた別の乱数とを比較し、波動関数の値の 方が大きければそこに電子がいるとし、小さければ棄却する。これで 原子核との距離を決める
3. 偏角をもう一度乱数を発生させて決定する

解答

      implicit real*8(a-h,o-z)
      real*8 pi, x0, dm, s, x, y, a, rr, xx, yy
      integer i
c
      pi = atan(1.0d0)*4.0d0
      x0 = 2.0d1
      dm = 0.1d0
      s = 100.0d0 / x0
c
      do i = 1 , 50000
         x = x0 * irnd() * 1.0d-8
         y = dm * irnd() * 1.0d-8
         if( f1s(x) .gt. y ) then
            rr = s * x
            a = 2.0d0 * pi * irnd() * 1.0d-8
            xx = rr*cos(a)
            yy = rr*sin(a)
            write(6,'(2e20.10)') xx , yy 
         end if 
      end do
      end
c==================================================================
      real*8 function f1s(r)
c
c     1s function for hydrogen atom
c
      implicit real*8(a-h,o-z)
      parameter(r0 = 0.529d0)
c
      f1s = ( 1/r0 )**3 * r*r * exp(-2.0d0*r/r0) 
c
      return
      end
c==================================================================
      integer function irnd()
      integer m1, m2, m3, i
      save m1, m2, m3
c
      data m1, m2, m3/92658393, 76438567, 65428733/
 1    continue
      i = m1 + m2 + m3
      if( m2 .lt. 50000000) i = i + 1357
      if( i .ge. 100000000) i = i - 100000000
      if( i .ge. 100000000) i = i - 100000000
      m1 = m2
      m2 = m3
      m3 = i
      irnd = i
c
      return
      end

• ５ー８行目
  必要なパラメータを与えます。piはπ、x0は動径方向の最大値、dmは採用する 乱数の最大値をそれぞれ表します。またx0,dmとあわせて次のdo文の繰り返し 回数である"50000"を適当に選んで、うまく表示させるように調整します。 これらをうまく選ばないと絵がつぶれてしまったり、採用される電子の数が 極端に少なくなったりします。
• do文の中
  "irnd()"によって乱数を発生させますが、得られるのは８桁の整数なので、 10^-8をかけて[0,1]の乱数とします。xは動径部分の値、yは 比較する乱数です。yとf1s(x)の値を比較し、f1s(x)の方が大きければ以下の 処理を行ない電子の座標を計算します。

今日の課題

またコンパイル回数、実行回数の報告は特にしなくて結構です。むしろ試行錯誤を いろいろやらないとうまくいかないものもあると思います。

1. 小テスト
2. [7-1]例題１の応用である。水素原子の2s,3s軌道を計算し、1s軌道と合わせて その違いを議論せよ。なおプログラムの提出はどちらか１つでよい。もちろん １つのプログラムで両方とも出力するようにしていてもかまわない。
3. [7-2]例題１の応用である。水素原子の2s,2p軌道を計算し、その違いを 議論せよ。なおプログラムの提出はどちらか１つでよい。もちろん １つのプログラムで両方とも出力するようにしていてもかまわない。
4. [7-3]例題１の応用である。水素分子の1σ軌道(1s+1s)を図示せよ。 原子間距離は平衡核間距離でよい。原子の場合と比較して共有結合の 本質を語ってくれるとなおよい。
5. [7-4]例題２の応用である。水素原子の2s,3s軌道を計算し、1s軌道と合わせて その違いを議論せよ。
6. [7-5]例題２の応用である。水素分子の1σ軌道、２σ軌道を図示せよ。 原子間距離は平衡核間距離でよい。
7. [7-6]例題２を応用して、水素原子の2p軌道を図示せよ。