Kalkulasi Energi Atom Helium: Aplikasi Hartree-Fock (ditulis dengan Fortran 90)

November 16, 2011

Menyelesaikan persamaan Schrödinger atom Helium pada keadaan dasar (l=0) dengan metoda Hartree-Fock dengan basis yang digunakan adalah gaussian.

Persamaan Schrödinger atom Helium dengan pendekatan Born-Oppenheimer bebas dari spin dengan mengikuti aturan antisimetrik fermion:

$\left[-\frac{1}{2}\bigtriangledown_{r1}^2 - \frac{1}{2}\bigtriangledown_{r2}^2 - \frac{Z}{r_1} - \frac{Z}{r_2} + \frac{1}{r_{12}} \right]\phi(r_1)\phi(r_2)=E\phi(r_1)\phi(r_2)$

Untuk mendekati energi keadaan dasar dan fungsi gelombang dari atom Helium (l=0), digunakan persamaan gelombang trivial yang merupakan kombinasi linear dari basis set fungsi Gaussian:

$\phi(r)=\sum_{p=1}^{N}C_{p}\chi_{p}$

$\chi_{p}=\exp{(-\alpha_{p}r^2)}$

Nilai α merupakan konstanta variasional yang diperoleh dari tabel. Pada atom Helium nilai dari
α=[0.297104, 1.236745, 5.74982, 38.216677]

Setelah disubtitusikan, diperoleh persamaan :

$\sum_{pq}\left(h_{pq}+\sum_{rs}C_rC_sQ_{pqrs}\right)=E\sum_{pq}S_{pq}C_{q}$

dengan elemen-elemen matriks:

$S_{pq}=\left(\frac{\pi}{\alpha_p+\alpha_q} \right)^{\frac{3}{2}}$

$h_{pq}=3\frac{\alpha_p\alpha_q\pi^{\frac{3}{2}}}{(\alpha_p+\alpha_q)^{\frac{5}{2}}} - 4\frac{\pi}{\alpha_p+\alpha_q}$

$Q_{pqrs}=2\frac{\pi^{\frac{5}{2}}}{(\alpha_p+\alpha_q)(\alpha_r+\alpha_s)\sqrt{\alpha_p+\alpha_q+\alpha_r+\alpha_s}}$

Persamaan di atas dapat diselesaikan dengan langkah-langkah di bawah ini:

Hitung komponen matriks NxN seperti h_pq dan S_pq, serta tensor NxNxNxN Q_pqrs.
Memberi nilai inisial sembarang dari vektor C, dalam contoh ini Cinisial =(1 0 0 0), dengan tujuan telah ternormalisasi
Nilai-nilai yang telah di dapat tersebut, digunakan untuk membangun fock matriks F_pq, yaitu :

$F_{pq}=h_{pq}+\sum_{rs}{C_rC_sQ_{pqrs}}$

Persamaan Fock di atas merupakan permasalahan generalzed eigen problem yang dapat diselesaikan. Untuk memperoleh nilai terkecil, fungsi eigen dari nilai eigen yang terkecil di ambil untuk digunakan sebagai input Fock di atas. Bagian ini merupakan bagian dari self consistency. Sebagai catatan, fungsi eigen C haruslah normalized, yaitu dengan syarat :

$\sum_{pq}C_pS_{pq}C_q = 1$

Dengan menentukan akurasi kekonvergenan, dalam contoh ini adalah 10^-8, setelah perbedaan tiap elemen C lama dengan elemen C baru memenuhi bilangan paling sedikit 10^-8, maka iterasi selesai. Bagian ini merupakan bagian akhir dari self consistency. Diperolehlah nilai C. Yang kemudian digunakan untuk menghitung energi atom!.
Nah, tinggal menghitung energi dari atom Helium sekarang!

$E = 2\sum_{pq}C_pC_qh_{pq} + \sum_{pqrs}C_pC_qC_rC_s$

Jika keenam langkah tersebut dituliskan dalam bahasa pemprograman Fortran 90, salah satunya akan menjadi :

PROGRAM HeliumHF
!-------------------------------------------------------------
! Kalkulasi total energi Helium dengan pendekatan Hartree-Fock
! by Cica, 12 Nov 2011
!-------------------------------------------------------------
IMPLICIT NONE

REAL(8), DIMENSION(:,:),ALLOCATABLE :: H,T,A,Sov,F,v
REAL(8), DIMENSION(:),ALLOCATABLE :: E,eps,c
REAL(8), DIMENSION(:,:,:,:), ALLOCATABLE :: Qpqrs
INTEGER :: n,p,q,r,s,i,j,info,l
REAL(8) :: pi,E2,E1
!koefisien basis
REAL(8) :: alfa(4)= (/ 0.298073d0, 1.242567d0, 5.782948d0, 38.474970d0 /)

n = 4
pi = 4.d0*ATAN(1.d0)   !nilai PI
!------Elemen-elemen matriks------------!
!--- Tensor Q_pqrs, ditulis perkolom
ALLOCATE (Qpqrs(n,n,n,n))
DO s=1,n
  DO r=1,n
    DO q=1,n
      DO p=1,n
         Qpqrs(p,q,r,s) = 2.d0*(pi**2.5d0)/(alfa(p)+alfa(q))/(alfa(r)+alfa(s))/&
                          sqrt(alfa(p)+alfa(q)+alfa(r)+alfa(s))
      ENDDO
    ENDDO
  ENDDO
ENDDO

!-----Matriks hamiltonian, H_pq=T_pq+A_pq

ALLOCATE (H(n,n), T(n,n), A(n,n))
DO q=1,n
  DO p=1,n
      T(p,q)=3.d0*alfa(p)*alfa(q)*(pi**1.5d0)/(alfa(p)+alfa(q))**2.5d0
      A(p,q)=-4.d0*pi/(alfa(p)+alfa(q))
      H(p,q)=T(p,q)+A(p,q)
   ENDDO
ENDDO

!----Overlap Matriks Sov
ALLOCATE (Sov(n,n))
DO q=1,n
 DO p=1,n
   Sov(p,q)=pi/(alfa(p)+alfa(q))
   Sov(p,q)=Sov(p,q)**1.5d0
 ENDDO
ENDDO 

!--- Self-consistency Hartee-Fock
ALLOCATE (F(n,n),c(n), eps(n), v(n,n))    

!---- menebak harga awal keofisien c(n) dahulu
c(:)=0.d0
c(1)=1.d0
E1 = 0.d0
!---- self consistency--------------------------------!

DO l=1,300

! ----- fock marix
 DO p=1,n
  DO q=1,n
    F(p,q) = H(p,q)
   DO s=1,n
    DO r=1,n
      F(p,q) = F(p,q) + c(r)*c(s)*Qpqrs(p,q,r,s)
    ENDDO
   ENDDO
  ENDDO
 ENDDO

   call diag(n,n,F,Sov,eps,v)

 c(:) = v(:,1)
!---- menghitung energi
  E2 = E1
  E1 = 0.d0
  DO p=1,n
   DO q=1,n
     E1 = E1 + 2.d0*c(p)*c(q)*H(p,q)
    DO r=1,n
     DO s=1,n
       E1 = E1 + Qpqrs(p,q,r,s)*c(p)*c(q)*c(r)*c(s)
     ENDDO
    ENDDO
   ENDDO
  ENDDO

PRINT*,E1
   ! do j=1,n
   !    write(*,'(1x,4F18.10)') F(:,i)
   !  enddo

ENDDO  

DEALLOCATE(eps,Qpqrs,c,Sov,T,A,H)

END PROGRAM HeliumHF

!----------------------------------------------------------------------
subroutine diag ( n, ldh, h, s, e, v ) !
!    diambil dari : http://www.fisica.uniud.it/~giannozz/
!
!       Finds eigenvalues and eigenvectors of the generalized problem
!       Hv=eSv, where H=hermitian matrix, S=overlap matrix
!       On return S and H are unchanged.
!       Uses level-2 blas dgemm for matrix-matrix multiplication
!       Uses lapack dsyev for matrix diagonalization, checks for
!       too small eigenvalues of S, removes corresponding eigenvectors
!
!-----------------------------------------------------------------------
!
  implicit none
  integer, parameter :: dp = selected_real_kind(14,200)

  integer, intent(in) :: n, &! dimension of the matrix to be diagonalized
                         ldh ! leading dimension of h,s and v, as declared
                             ! in the calling pgm unit
  real (dp), intent(in) :: h(ldh,n), & ! matrix to be diagonalized
                           s(ldh,n)    ! overlap matrix
  real (dp), intent(out):: e(n), &     ! eigenvalues
                           v(ldh,n)    ! eigenvectors (column-wise)
  !
  ! local variables
  integer lwork, info, i, j, nn
  !
  real (dp), parameter :: small=1.d-10
  real (dp), allocatable :: work(:), aux(:,:), h1(:,:)
  !
  info = 0
  lwork=3*n
  allocate (work(lwork), aux(ldh,n))
  !
  !       Copy S into an auxiliary matrix because dsyev destroys the matrix
  !
  aux = s
  !
  !       Diagonalize S
  !
  call dsyev ( 'V', 'U', n, aux, ldh, e, work, lwork, info )
  if (info /= 0) stop 'S-matrix diagonalization failed '
  !
  !   Keep only linearly independent combinations (within a given threshold)
  !   store into matrix "aux" the eigenvectors of S divided by the squares
  !   of the eigenvectors
  !
  nn = 0
  do i=1,n
     if (e(i) > small) then
        nn = nn + 1
        aux(:,nn) = aux(:,i) / sqrt(e(i))
     end if
  end do
  if ( nn < n) write(*,*) " # of linearly independent vectors =", nn, n
  !
  !   Transform H using the "aux" matrix
  !
  !   V(i,j) = \sum_{k=1}^{n} H(i,k) aux(k,j),  i=1,n, j=1,nn
  !
  call dgemm ( 'N', 'N', n, nn, n, 1.0_dp, h, ldh, aux, ldh, 0.0_dp, v, ldh )
  !
  !    h1(i,j) = \sum_{k=1}^{n} aux(k,i) v(k,j),  i=1,nn, j=1,nn
  !    H' = transpose(aux)*H*aux
  !
  allocate (h1(nn,nn) )
  call dgemm ( 'T', 'N', nn, nn, n, 1.0_dp, aux, ldh, v, ldh, 0.0_dp, h1, nn )
  !
  !    Diagonalize transformed H
  !
  info = 0
  call dsyev ('V', 'U', nn, h1, nn, e, work, lwork, info)
  if (info /= 0) stop 'H-matrix diagonalization failed '
  !
  !    Back-transform eigenvectors
  !
  call dgemm ( 'N', 'N', n, nn, nn, 1.0_dp, aux, ldh, h1, nn, 0.0_dp, v, ldh )
  !
  deallocate (h1, aux, work)
  !
end subroutine diag

Jika program dikompilasi dan dijalankan, output akan menjadi :

Selamat mencoba

Like this:

Be the first to like this post.

4 Komentar leave one →

Kasyful Fuadi permalink

November 16, 2011 2:15 pm

mantap betul… glek glek… baca ‘Schrodinger’ saja sudah merinding …

Balas
- チ茶 permalink*
  
  November 17, 2011 2:10 am
  
  Silakan dicoba kk~
  
  Balas
Indra permalink

November 17, 2011 12:42 pm

Mantap caaa. Calon scientist sejati!

Balas