Skip to content

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:

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

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

  5. 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 :

  6. \sum_{pq}C_pS_{pq}C_q = 1

  7. 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!.
  8. 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😀

15 Komentar leave one →
  1. Kasyful Fuadi permalink
    November 16, 2011 2:15 pm

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

  2. Indra permalink
    November 17, 2011 12:42 pm

    Mantap caaa. Calon scientist sejati!

  3. Maret 18, 2012 3:43 am

    wow! :-O

  4. Desember 9, 2014 12:41 pm

    mbak, nanti itu keluarannya dalam bentuk file o bukan? cara liat hasilnya gimana yah?

    • Desember 13, 2014 10:09 pm

      File .o itu keluaran binary yg belum di-linker.
      Tergantung dari OS, kalau ga salah keluaran di Windows biasanya .exe.
      Kalau misalnya pakai linux, karena ini code cuman kecil dan satu file bisa begini :

      $ gcc file.c -o hartreefock
      $ ./hartreefock

      ekstensi ga masalah dan linker ga perlu kalau cuman satu file

  5. Fika permalink
    Maret 19, 2015 7:48 am

    Maaf saya mau tanya, saya kan sudah coba source code diatas tapi kenapa pas di compile di fortran 90 (for windows) gak bisa ada, dan muncul sbb :

    C:\Users\Rafika\AppData\Local\Temp/ccEX2oEJ.o:HeliumHF.f90:(.text+0x1adf): undef
    ined reference to `_dsyev_’
    C:\Users\Rafika\AppData\Local\Temp/ccEX2oEJ.o:HeliumHF.f90:(.text+0x1f34): undef
    ined reference to `_dgemm_’
    C:\Users\Rafika\AppData\Local\Temp/ccEX2oEJ.o:HeliumHF.f90:(.text+0x200f): undef
    ined reference to `_dgemm_’
    C:\Users\Rafika\AppData\Local\Temp/ccEX2oEJ.o:HeliumHF.f90:(.text+0x2063): undef
    ined reference to `_dsyev_’
    C:\Users\Rafika\AppData\Local\Temp/ccEX2oEJ.o:HeliumHF.f90:(.text+0x2109): undef
    ined reference to `_dgemm_’

    itu maksudnya bagaimana ya? mohon pencerahannya🙂

    • Maret 19, 2015 12:11 pm

      Hai Fika,
      Jadi kesalahannya fungsi _dsyev, _dgemm tidak ditemukan.
      _dsyev itu pake dari LAPACK dan _dgemm itu paket dari BLAS,
      Fika harus install BLAS dan LAPACK dulu.

      • Fika permalink
        Maret 19, 2015 10:40 pm

        oh oke. Makasih.
        Oya utk perhitungan atom elektron banyak (misal lithium, berilium, dll) sudah pernah coba juga kah?

      • Maret 20, 2015 1:36 am

        Jelas bisa dicoba untuk elektron banyak, tapi akurasi tidak menjamin, berhubung metoda HF melakukan simplifikasi dengan satu Slater Determinant sebagai basis

  6. Fika permalink
    Maret 20, 2015 1:13 am

    ohya maaf saya mau tanya lagi, saya belum berhasil utk install BLAS & LAPACK, kalau boleh tau bagaimana ya prosedur instalasinya? Terima kasih

Trackbacks

  1. Calculate Helium Ground State Energy using Hartree-Fock « Celestial Aviary

Tinggalkan Balasan

Isikan data di bawah atau klik salah satu ikon untuk log in:

Logo WordPress.com

You are commenting using your WordPress.com account. Logout / Ubah )

Gambar Twitter

You are commenting using your Twitter account. Logout / Ubah )

Foto Facebook

You are commenting using your Facebook account. Logout / Ubah )

Foto Google+

You are commenting using your Google+ account. Logout / Ubah )

Connecting to %s

%d blogger menyukai ini: