ProsidingPertemuanI/miah SainsMateri 1997
/SSN/4/0-2897
DINAMlKA MOLEKULER KW ANTUM TIGHT-BINDING SEBAGAI ALA T BANTU STUDI MATERIAL: STUDI KASUS SISTEM SILIKON1 Hermawan Kresno Dipojono2
ABSTRAK DINAMIKA MOLEKULER KWANTUM TIGHT-RINDING SEBAGAI ALAT BANTU STUDI MATERIAL: STUDI KASUS SISTEM SILIKON. Metodologi dasar dinamika molekuler kwantum tight-binding sebagai alat untuk mempelajari sifat sifat materi akan dibahas dalam makalah ini. Contoh contoh yang disajikan di sini diambil daTi eksperimen komputer kami mensimulasikan sistem atom dan elektron. Makalah ini memilih sistem sistem silikon sebagai objek penelitian. Banyaknya hasil hasil penelitian, baik eksperimental maupun teoritis-analitis, dari sistem sistem silikon menyebabkan pemilihan ini merupakan kewajaran dalam rangka menguji keabsahan metoda yang digunakan di sini. Hasil hasil simulasi kwantum sistem silikon yang telah kami lakukan dalam lingkungan pada temperatur tertentu akan disajikan di sini.
ABSTRACT TIGH}'-BINDING MOLECULAR DYNAMICS AS A TOOL FOR STUDYING MATERIAL PROPERTIES: SILICON SYSTEM AS A CASE STUDY. The basic methodology of tight-binding quantum molecular dynamics as a tool for studying material properties is described. Examples are taken from our experience using tight-binding molecular dynamics to simulate system of atoms and electrons. This paper chooses the silicon systems as the object of study. The vast available results, both experimental and theoretical, of silicon systems makes it a natural choice to test the correctness of the method to study material properties. Results from our quantum simulations of silicon system at a finite temperature will be presented.
KEY WORD Molecular dynamics,Tight-binding,Simulations
PENDAHULUAN
Oinamika molekuler kwantum tight-binding memungkinkan tidak saja simulasi sistem dengan Dinamika molekuler pada dasamya jumlah atom yang lebih besar tetapi juga simulasi merupakan solusi numerik dari persoalan dinamika untuk pengamatan fisis yang lebih lama. Karena klasik N-buah atom, meskipun bisa saja gaya gaya alasan inilah kami memilih untuk menggunakan yang bekerja antar partikel itu diturunkan secara Hamiltonian tight-binding parametrik sebagai alat kwantum. Secara umum, dalam menggunakan bantu untuk melakukan studi struktur elektronika metoda dinamika molekuler ini, terdapat dua material. persoalan utama: (a) memilih energi potensial antar atom yang tepat, (b) memilih metoda daD algoritma ENERGI TOTAL TIGHT-BINDING numerik yang effisien. Dengan memilih potensial antar atom berupa pendekatan kwantum tightKami menggunakan energi total sistem binding maka makalah ini memang akan lebih silikon dari Tomanek daD Schluter [5,6] yang menitik beratkan pembahasanpada persoalankedua. didasarkan pada tight-binding parametrik dari Chadi Dalam sepuluh tahun terakhir ini, [7]. Agar dapat digunakan untuk maksud maksud kemungkinan menggunakan Hamiltonian tightdinamika molekuler kami mengikuti Khan dan binding untuk mempelajari sifat dinamika condensed Broughton [I] yaitu memasukkan suatu fungsi matter terns dikembangkan (lihat, misalkan di potong (cutoff) kedalam energi total tersebut. [1,2,3]). Dalam metoda ini, pengaruh elektron Dengan demikian untuk sistem silikon yang terdiri digabungkan kedalam dinamika molekuler melalui dari N buah atom dengan koordinat {R,} , snafu Hamiltonian tight-binding parametrik. Di sini, (f.=1,2,...,N), yang masing masingnya memiliki 4 berbeda dengan Car-Parrinello di mana fungsi buah elektron valensi akan mempunyai status bebas gelombang elektron diuraikan menjadi ribuan {lfIj}' (i=I,2,...,N). Selanjutnya status ini diuraikan gelombang bidang, perhitungan yang melibatkan elektron hanya memerlukan sejumlah kecil orbital kedalam suahl ba~js {f/Jea}yang terdiri atas empat untuk setiap atom. Dengan demikian, kelebihan orbital s, Px, py. dan P. yang terletak di s~tiap atom; pendekatan tight-binding terhadap metoda Cardengan kata lain, f. menyatakan indeks atom dan a Parrinello [4] yang menggunakan formalisma menyatakan indeks orbital sebagaimana yang fungsional kepadatan (dengan pendekatan kepadatan ditunjukkan oleh persamaan(I). lokal -local density approximation) dari aspek waktu komputasi amatjelas. Ener?i total sistem atom dan elektron £'0' \{R f }, ~~a }), yang dinyatakan oleh persamaan (2), dibangun berdasarkan energi total Tomanek dan
~Dipresentasikan pada Pertemuan Ilmiah Sains Materi 1997 2Jurusan Teknik Fisika Institut Teknologi Bandung
445
Schluter, dan dimodifikasi untuk maksud dinamika molekuler berdasarkan Khan dan Broughton [9].
IfI;(r)= f fC~atjJfa(r)
(1)
f=la=1
(2)
£'01 = £bs +£2 -E)
(3)
(4)
[
£3 = N I/JI( N nb
2
( nb
)
dapat dianggap sebagai semata-mata problema elektron dengan posisi atom diasumsikan tetap (fIXed) dalam rentang pengamatan. Keabsahan pendekatan ini sangat didukung oleh kenyataan bahwa massa atom jauh lebih besar dibanding massa elektron. Hal ini berarti bahwa energi total barns diminimisasikan terhadap derajat kebebasan elektron untuk setiap konfigurasi posisi atom. Proses ini dalam literatur disebut sebagai proses quenching elektron pada permukaan Born-Oppenheimer (B08Born-Oppenheimer Surface). 8etelah proses ini diselesaikan, E'°/({ R, }.{c~} )1 {c} di 80S merupakan energi pptensial effektif dari atom atom dan trajektorinya diturunkan berdasarkan Lagrangian yang diberikan oleh persamaan(8).
(5)
+1/12N 1-4'3
(8)
L }is;., tot
E,.{R) = ESi (R)fc{R-Rb)-Esi 2
N I-I
(I
nb = L Lfc\jRI
-Rrl-Rb
)
bs
(R)
(6)
(9)
(7)
1~2r~1
g." adalah energi struktilr pita yang dibentuk oleh Chadi melalui model tight-binding parametrik tetangga terdekat (nearest neighbor). Elemen matrik (If! ia IHIIf! i'a') untilk orbital orbital pada atom atom e dan e' yang berjarak sebesar R mengandung fungsi potong, !c(R). £1,\,;adalah energi atom netral terisolir dan besamya sarna dengan 2(E,,+Ep), di mana
E,\ = (t/Ji,' IHIt/Ji.') dan Ep = (t/JiXYZ IHI t/J,xyz) .
E2 adalah energi tolak (repulsive energy), diperlukan untuk menetralisir pengaruh penghitilngan duakali (double counting) interaksi Coulomb pada suku struktilr pita g.s, dan juga mengandung kontribusi energi energi korelasi serta exchange dan interaksi atom-atom. Fungsi Er yang muncul dalam E] diberikan oleh persamaan (6). E~' daD E~:2 berturut-tilrut adalah energi total dan energi struktilr pita dari silikon diatomik, dan Rb adalah lokasi titik
Berdasarkan Lagrangian ini maka persamaan persamaangerak dari atom atom diberikanoleh persamaan(9). Trajektori atom atom akan memenuhi hukum kekekalan energi sebagaimana diberikanoleh persamaan (10). Dengan demikian sifat sifat makroskopis yang diperoleh dalam dinamika mo1ekulerberdasarkanLagrangian di atasmemiliki (E. V; N) yang konstan,dengankata lain merupakanensemblemikrokanonikal,di manaE adalahenergi total, V adalahvolume sistem,dan N merupakanjumlah atom. INTEGRASI PERSAMAAN GERAK Solusi persamaan persamaan gerak ditentukan diantaranya oleh jenis fungsi potong yang digunakan. Secara garis besar fungsi potong dapat dibedakan menjadi dua: (a) fungsi potong yang landai, daD (b) fungsi potong yang berupa fungsi step. Contoh fungsi potong yang landai misalnya yang digunakan oleh Khan daD Broughton [I]
potong. Untilk E.~~ kami menggunakan energi total molekul diatomik dari Raghavachari [8]. Suku ketiga E3 dimasukkan kedalam energi total agar sistem ini fc (R-Rh ) =21 I-tanh ~R-Rh )] (11) dapat menghasilkan struktilr kubik, intan maupun body-centered-cubic. Koreksi ini perlu dilakukan di mana Rbmenyatakan titik potong (panjang ikatan), untilk mengantisipasi perubahan energi akibat !:!. adalah ukuran kelandaian fungsi potong, daD R perubahan jumlah ikatan kimia nb yang dinyatakan menyatakan jarak antar atom yang berinteraksi. Perlu oleh persamaan (7) dan parameter parameter cIIl, c112, ditambahkan di sini bahwa harga Rbberkisar diantara daD <1>3berturut-turut memiliki harga 0,225 eV, harga jarak tetangga terdekat (nearest neighbor) (Ro) 1,945 eV, dan -1,03 eV. daD tetangga terdekat berikutnya (second-nearestMETODA DINAMlKA MOLEKULER neighbor), yang dalam hat silikon dengan struktur Berdasarkan pendekatan Bom- intan adalah 1,63 Ro. Sedangkan untuk harga !:!. fiipilih sekecil mungkin sehingga harga fungsi Oppenheimer, problema sistem atom clan elektron
[
(
ProsidingPertemuanl/miah SainsMateri 1997
/SSN/4/0"2897
potong itu dapat diabaikan jika harga R mendekati jarak tetangga terdekat kedua. Meskipun demikian perlu dicatat bahwa untuk maksud simulasi dinamika molekuler harga d diinginkan cukup besar agar ~ntuk waktu pengamatan proses yang wajar tidak diperlukan terlalu banyak langkah integrasi. Kompromi harus dicari diantara kedua persyaratan yang saling bertentangan ini. Di samping fungsi potong landai seperti di atas dikenal pula fungsi
yang berpusat di sebuah atom yang tertentu. Secara numerik g( r ) dihitung dengan cara mengarnbil harga rata-rata untuk semua atom di dalam boks simulator sebagai titik pusat. Pada padatan kristal g(r) sarna dengan jumlah fungsi fungsi delta yang berpusat di kisi kisi kristal. Mengembangnya kaki kaki fungsi delta di garnbar 1 disebabkan oleh pergeseranacak dari posisi keseimbangannya.
potong landai lain yang juga banyak digunakan dalam simulasi dinamika molekuler, lihat misalnya Goodwin et al [9]. Kami tidak akan membahas lebih banyak lagi mengenai fungsi potong landai ini. Fungsi potong berupa step untuk maksud simulasi dinamika molekuler pertama kali diperkenalkan oleh Dipojono clan Khan [10]. Solusi ini ditawarkan. untuk mengatasi dua hal yang saling bertentangan bila menggunakan fungsi potong yang landai seperti diuraikan di atas. Integrasi persamaan gerak yang ditawarkan mereka temyata merupakan bentuk umum dari algoritma Verlet. Dalam makalah ini kami tidak akan membahas lebih rinci mengenai penurunan algoritma mereka. Cukup bila disampaikan di sini bahwa persamaan integrasi merekalah yang kami gunakan dalam studi ini.
1500
~ "5
0
Gambar
2
4 IisIanca
&
8
rIAl
Korelasi pasangan kristal 64 atom silikon dengan struktur intan. Kedua puncak pertama menunjukkan tetangga terdekat dan tetangga terdekat kedua.
HASIL SIMULASI DAN PEMBAHASAN SimuJasi dinamika molekuJer dilakukan didalam kubus simulator dengan sisi L=I,48 A. SeJanjutnya sebanyak 64 atom silikon dengan struktur intan ditempatkan di daJamkubus tersebut. Syarat batas periodik diberJakukan dalam simulasi ini sehingga setiap atom yang berada di kubus simulator memiliki 26 citra yang berkoordinat sarna dengan koordinat aslinya ditambah atau dikurangi L. Kepadatan sistem terjaga konstan karena setiap ada atom yang keJuar dari satu sisi kubus simulator maka bersamaan dengan itu akan masuk pula atom barn dari sisi yang berlawanan. SeJang waktu integrasi yang digunakan adalah 10-15 detik kecuali jika disebutkan lain. Source code untuk semua perhitungan di sini kami tulis dengan FORTRAN clan kami proses di IBM RISC 6000/550. Gambar 1 menunjukkan koreJasi pasangan (pair correlation) dari padatan kristaJ 64 atom silikon berstruktur intan. Kedua puncak pertama menunjukkan jarak tetangga terdekat clan tetangga kedua. Oistribusi pasangan ini dihitung dengan cara yang digunakan oJeh Rahman [11] di mana g( r ) didetinisikan sebagai
(12) ~] di mana n( r ) adalah jumlah atom yang berada
:f)[4J1r2Ar
dalam jarak antara r dan r+Ar dari suatu atom, N adalah jumlah atom dalam volume V. Atau dengan kala lain, terdapat sebanyak (V/N) g( r ) 47t r2 ,1r atom di dalam cell bola antara jari-jari r clan r+Ar
447
~ ~
> ~
S]O
~ I ,~ .1~
"---~
--
.-
0
50
100
rurmer of118..-
'50
200
Gambar2. Konvergensienergiversusjumlah iterasiyang diperlukan oleh steepestdescentuntuk kristal 64 atom silikon dengan struktur intan pada suhurendah. Sebagaimana telah kami bahas di bab sebelumnya, untuk temperatur terbatas elektron elektron dapat dipertahankan di status dasarnya
dengan
menggunakan
pendekatan
Bom-
Oppenheimer. Jadi, di setiap saat, untuk suatu konfigurasi posisi atom atom, energi total sistem hams diminimumkan terhadap koordinat koordinat elektron dengan memperhatikan kendala ortonormal. Untuk hal tersebut semula kami menggunakan metoda steepest descent. Namun demikian perlu kami kemukakan di sini bahwa waktu untuk proses quenching ini pada dasamya adalah perhitungan orde N3. Jadi, untuk 64 atom saja, proses ini sebanyak dua kali saja telah memerlukan waktu kurang lebih 700 detik CPO di IBM RISC 6000/550. Gambar 2
menunjukkan konvergensi energi versus jumlah iterasi yang diperlukan dengan steepest descent untuk 64 atom silikon dalam struktur intan terpadatkan yang posisinya digeser sedikit dari posisi keseimbangannya kurang lebih sekitar 10-7 A. Meskipun pergeseran itu amat kecil temyata memerlukan sekitar 200 iterasi untuk mendekati harga energi yang benar. Oleh karena itu untuk temperatur tinggi, di mana atom atom bergeser semakin jauh dari posisi keseimbangannya, metoda steepest descent mengalami kesulitan untuk mencapai harga energi minimum. Gambar 3 menunjukkan hal ini. Di sini simulasi diadakan pada temperatur tinggi daD nampak bahkan setelah iterasi sebanyak 500 kali konvergensi belum juga dicapai. Dalam situasi seperti ini diagonalisasi langsung matrik Hamiltonian, selama jumlah atom masih memungkinkan untuk hal tersebut, nampaknya merupakan satu satunya solusi.
30
I! !
c;;
r(A)
Gambar4. Fungsikorelasipasangang(r) versusjarak antar atom, Garis titik-titik dengan bulatan hitam adalah data dari difraksi sinar-X, pada suhu 1700K. Garisterputusdari simulasipadasuhu 3300 K. Garis penuh diperoleh dari simulasi padasuhu1800K.
E(eV)
Gambar 5. Status kepadatanenergi versus energi yang diperolehdari simulasisepertigambar4 pada suhu1800K. Dengan diketahuinya konfigurasi (bahkan posisi rinci setiap) atom di setiap saat maka diperoleh kekayaan informasi yang cukup lengkap daTi sistem atom dan elektron ini. Korelasi pasangan atom, misalnya, dapat dengan mudah diperoleh; juga kepadatan status (density of states) dapat dihitung dengan cepat. Gambar 4 dan 5 menunjukkan informasi mengenai hal hal tersebut. Korelasi pasangan atom untuk status cair (temperatur di atas suhu lebumya) ditumpang-tindihkan dengan basil diffraksi sinar x dan menunjukkan basil yang cukup memuaskan, meskipun waktu simulasi yang amat terbatas.
448
KESIMPULAN Oalam makalah ini telah dibahas metoda dinamika molekuler kwantum tight-binding untuk mensimulasikan sistem atom daD elektron. Contoh yang disajikan di sini diambil daTi simulasi sistem silikon pacta temperatur lebur yang telah kami lakukan. Silikon pacta temperatur caimya bersifat logam daD basil simulasi kami menunjukkan hal tersebut dengan cukup baik. Pasangankorelasi untuk sistem silikon cair telah diketahui secara umum tidak memiliki long range order. Hasil simulasi kamipun menunjukkan hal itu dengan baik; sebagai pembanding basil simulasi dibandingkan dengan basil diffraksi sinar x. Oi samping itu, karena sifat logamnya, sistem silikon cair memiliki celah pita nolo Hal ini ditunjukkan pula oleh kepadatan status daTi basil simulasi kami. Oengan demikian metoda simulasi kwantum ini dapat digunakan sebagai alat bantu analisis sistem atom daDelektron.
ProsidingPertemuanI/miah SainsMateri /997
/SSN/4/0-2897
UCAPAN TERIMA KASIH Penelitian ini sebagian diantaranya dibiayai oleh Departemen Pendidikan dan Kebudayaan Pemerintah Republik Indonesia melalui program New-S 1 Ilmu Material di Institut Teknologi Bandung. Biaya juga diberikan oleh Kementrian Riset dan Teknologi Republik Indonesia melalui Riset Unggulan Terpadu V 1997/1998. Penulis menyampaikan terima kasih atas bantuan dana yang memungkinkan penelitian dan penulisan ini dapat dilakukan.
DAFTARPUSTAKA [1] KHAN, F. S, and BROUGHTON, J. Q., Simulation of Silicon Clusters and Surfaces via Tight-Binding Molecular Dynamics, Phys. Rev. 8,43 (1991), 11754. [2] WANG, CHAN, C. T., and HO, K. M., Tight Binding Molecular Dynamics Study of Defects in Silicon, Phys. Rev. B, 45 (1992),12227. [3] VIRKKUNEN, LAASONEN, K., and NIEMINEN, R. M., Molecular Dynamics Using the Tight-Binding Approximation: Application to Liquid Silicon, J. Phys: Condens. Matter, 3 (1991),7455.
449
[4] CAR and PARRINELLO, M., Unified Approach for Molecular Dynamics and Density Functional Theory, Phys. Rev. Lett., 55 (1985),
2471. [5] TOMANEK and SCHLUTER, Calculation of Magic Number and the Stability of Small Si Clusters, Phys. Rev. Lett., 56 (1986), 1055. [6] TOMANEK and SCHLUTER, Structure and Bonding of Small Semiconductor Clusters, Phys. Rev. B, 36 (1987),1208. [7] CHAD!, Theoretical Study of the Atomic Structure of Silicon (2] I), (31 I), and (33]) Surfaces. Phys. Rev. B, 29 (1984), 785. [8] RAGHA VACHARI, Theoretical Study of Small Silicon Clusters: Cyclic and Ground State Structure of Si3, J. Chern. Phys., 83 (]985),
3520. [9] GOODWIN, SKINNER, A. J., and PElTITFOR, D. G., Generating Transferable Tight-Binding Parameters: Application to Silicon, Europhys. Lett., 9 (1989), 70]. []O]D!POJONO and KHAN, F. S., A Regional Seminar on Computer Method and Simulation in Engineering, Bandung (1997) (accepted paper for presentation). [] ]]RAHMAN, Correlations in the Motion of Atoms in Liquid Argon, Phys. Rev., 136 (]964), A405.
