p~
~
N~
H,t..,
Nt..t.- M.. ~
X KL4,ISSN 1410,6g6
SIMULASI DINAMIKA MOLEKUL (SEBUAH PENGANTAR) HermawanK. Dipojono LaboratoriumKomputasiMaterial, DepartemenTeknikFisiko ITB .n. Ganesha10 Bandung40132 e-mail:
[email protected]
ABSTRAK SIMULASI DINAMIKA MOLEKUL (SEBUAHPENGANTAR). Makalah iniakanmembahas simulasidinamika molekulsebagaisebuah metodaalternatif yangdapatdigunakan untukmempelajari sifatsifatataukarakteristik suatumateri.Beberapa hasilsimulasiuntuksistematom danelektronakandisajikanpulasebagaisuatucontohkegunaan dankesederhanaan metoda ini.
ABSTRACT MOLECULER DYNAMICS SIMULATION (AN INTRODUCTORY). Moleculer dynamics simulation as analtemative methodethatcanbe for studyingmaterialscharacteristics is presented in thispaper.A numberof simulation resultsfor systemos atomandelectronis discussed and presented asexamplesofits spectrum andsimplicity.
PENDAHULUAN Simulasi dinamika molekul merupakansebuah metodayang dapatdigunakanuntukmelakukanprediksi terhadap sifat sifat statik maupun dinamik yang diturunkansecaralangsungdari interaksiditingkat atom atau molekul. Mengingat belum actaaltematif lainnya yang dapatdigunakanuntuk memecahkanpersoalanitu sampai ketingkatan yang cukup rinci maka metoda simulasidinamika molekul ini merupakansesuatuyang tidak terhindarkan dalam penelitian baik untuk ilmu mumi maupunrekayasa.Metoda ini, sejak pertamakali dipopulerkan kembali oleh Alder dan Wainwright [1], telah digunakandenganamatluas oleh berbagaidisiplin ilmu, termasukdalam ilmu dan rekayasamateri [2-6]. Makalah ini akan mengupas secara rinci mengenai
metodaini dan menyajikansejumlahbasil simulasiyang telahdilakukanuntuk sistematomdanelektron. Pertanyaanyang sering munculdalam berbagai konteks penelitian ilmu dan rekayasa materi adalah adakah relasi antara sifat sifat makroskopis materi denganinteraksiyang actadi tingkat atom atau molekul penyusunmateri tersebut. Sudah barang tentu melalui eksperimen,yang pasti tidak murah dan tidak pula sederhana,dapat dideduksi sifat sifat mikroskopisyang berpengaruhterhadap sifat sifat makroskopismateri. Dinamika molekul memberikansebuahaltematif untuk mereproduksi tingkah laku mikroskopis itu dengan
~,
menggunakan sistem sistem model. Kemajuan teknologi komputer telah memungkinkan untuk mencari jawab persoalan yang lebih kompleks disertai meningkatnya harapan untuk bisa mendapatkan jawaban yang lebih
realistis. Terdapat dua buah persoalan utama yang hams diselesaikan dalam menggunakan metoda simulasi dinamika molekul ini. Pertama adalah menentukan model energi potensial yang mengatur hubungan antara atom atom atau molekul molekul yang saling berinteraksi di dalam sistem. Dari model energi potensial inilah gayagaya yang mempengaruhi dinamika sistem dapat ditentukan. Kedua, menentukan algoritma daD metoda numerik yang tepat daD efisien untuk digunakan dalam memecahkan persoalan dinamika yang amat kompleks ini. Makalah ini akan memulai pembahasannya dengan menyoroti persoalan pertama terlebih dahulu. Pada bagian ini akan dibahas terlebih dahulu model model energi potensial klasik-empirik daD diakhiri dengan model yang menggunakan sepenuhnya pendekatan mekanika kwantum. Kemudian makalah ini akan membahas pula sejumlah algoritma daD metoda numerik yang banyak digunakan dalam studi materi dengan menggunakan metoda ini. Termasuk dalam bagian ini juga akan dikemukakan algoritma untuk melakukan simulasi dengan temperatur daDtekanan yang terkontrol.
6J~ 2001
~
/)~"~"--~H..llIt..l(S1t1.4 P~)
H~~ MODELENERGIPOTENSIAL Keabsahanmodel materi di tataranmikroskopis sangat ditentukan oleh pemahaman komprehensip mengenai komponen komponen penyusun materi itu. Meskipundeskripsimengenaihal ini secaraprinsipharus didasarkanpada mekanika kwantum namunpendekatan mekanika klasik, yaitu memperlakukanatom atom atau molekul molekul sebagaisuatutitik massa,dalam batas batas tertentu masih dapat digunakan. Model paling sederhanayang dapat digunakan untuk menjelaskan suatu materi didasarkanpada konsep partikel partikel berbentuk bola yang saling berinteraksi yang untuk kemudahandalam pembahasanselanjutnyapartikel itu disebut atom. lnteraksi itu, dalam tataran yang paling sederhana,terjadi antar pasanganatom yang memenuhi dua buah kriteria dasar hubunganantar atom. Pertama, interaksi itu harus mampu menahantekananpasangan atomyang saling berinteraksidaDkonsekuensinya adalah pasanganyang saling mendekatakan menghadapigaya tolak menolak. Kedua, interaksi itu mempunyai sifat mengikat pasanganatom sehinggapasanganitu akan tarik menarikjika saling menjauh.Model sederhanaini kemudian terns berkembang hingga mencapai kemantapannya dengan dilibatkannya mekanika kwantum untuk menanganihadirnya partikel elementer sepertielektron. Dalam bagian ini kita akan membahas model energi potensialpaling sederhanaterlebih dahulu sebelum membahas model ab initio yang lebih komprehensip.
K.~
sehingga waktu komputasi dapat dihemat dengan cara mengabaikan interaksi pasangan atom yang mempunyai jarak lebih besar daTiRc. Dalam banyak kasus biasanya dipilih Rc=2,500; di mana U(Rc)=-0,0163E dan besarnya gaya adalah -0,039E /0". Parameter E, mempunyai satuan energi sedangkan parameter 0", mempunyai satuan panjang. Untuk tujuan dinamika molekul berbasis potensial Lennard-Jones ini sering juga digunakan satuan satuan tanpa dimensi yaitu dengan menggunakan 0; untuk satuan panjang m, untuk satuan massa, dan E, untuk satuan energi. Dengan kata lain mengganti r ~ rO" untuk panjang, e ~ eE untuk energi, dan t~ t w"d/E untuk waktu. Jika model energi ini digunakan untuk keperluan memodelkan sistem argon cair [7] maka relasi antara satuan satuan tanpa dimensi dan satuan fisis dinyatakan dalam besaran berikut. Panjang dinyatakan dalam E = 3,4A. Satuan energi diberikan oleh E/kB = 120K, yang berarti bahwa E = 120 x S,314Joule/mole. Dengan menggunakanmassa atom argon adalah m=39,95 x 1,6747 x 10-24gram, maka satuan waktu dinamika molekul setara dengan 2,161 x 1012,detik; jadi sebuah selang waktu yang umum dipakai sebesar At = 0,005, adalah setara dengan kurang lebih 10-24, detik. Selanjutnya jika Na, adalah jumlah atom yang berada di dalam boks simulator dengan panjang sisi sisinya adalah L, maka untuk rapat jenis cairan argon sebesar 0,942 grarn/cm3, akan berarti bahwa L = 4,142 Nal/3 A, yang dalam satuan Yal1g telah direduksi akan menjadi L = I,2ISNaI/3
PotensialLennard-Jones Parameterisasi Ikatan Kuat Model energi potensial yang dianggappaling Parameterisasi ikatan kuat (tight binding sederhana,yang pada mulanya dimaksudkan sebagai parameterized energy model) merupakan salah satu model cairan argon, adalah potensial Lennard-Jones. model yang cukup populer karena di samping dapat Dalam model ini untuk pasanganatom denganindeks I digunakan untuk memodelkan sistem yang melibatkan daDl' yang terletakberturut-turutdi koordinat~ daDRi, elektron elektron juga komputasiny relatif lebih hemat energipotensialnyadinyatakanoleh :
U(Rul) = 4£ (~)12 Ru'
-(~)6
(1)
Rll'
di mana Rti' = I Ri _Ri' I yaitu jarak antara atom / clanr. Parameter e , menyatakan kekuatan interaksi sedangkan parameter cr, menyatakan suatu skala panjang; interaksi itu akan tolak menolak pada jarak yang dekat clan akan tarik menarik pada jarak yang relatif jauh, serta akhirnya akan tidak berpengaruh sama sekali pada suatu jarak yang lebih besar dari jarak potong (cut oj]). Potensial interaksi Lennard-Jones di atas dapat disederhanakan lebih lanjut dengan mengabaikan gaya tarik menarik yang terjadi di bagian yang mendekati Dol yaitu dengancara merubah persamaan (I) menjadi :
U(Rtll) =
4£ 0
2
[(~j12 Ru'
~ {~)8 Rll'
] +a
Rei' ~ Rc Ru' > Rc
~/6J.,.,
waktu (dibandingkan dengan model ab initio) yang akan dibahas di bab selanjutnya. Model energi ini banyak digunakan untuk studi bahan bahan semikonduktor baik untuk kristal tunggal maupun amorf. Parameter parameter di dalam model ini ditentukan, misalnya, dengan menggunakan metoda fitting terhadap data daTi eksperimen maupun data daTi model yang lebih akurat (seperti data model energi ab initio ). Dalam kesempatan ini akan dibahas secara rinci sebuah contoh dari model energi ikatan kuat untuk sistem silikon. Sengaja dipilih silikon karena sistem ini mempunyai sifat sifat yang menarik, yaitu pada temperatur OaK, bersifat isolator, sedangkan pacta rasa caimya bersifat konduktor, dan pada temperatur diantara keduanya sistem silikon bersifat semikonduktor. Di samping itu data data sistem silikon amat lengkap sehingga merupakan sistem yang tepat untuk dipilih sebagaicontoh. Misalkan sistem silikon tersebut terdiri atas N, buah atom dengan koordinatnya dinyatakan oleh ~ dimana I = 1,2, N menyatakan identitas atom atom tersebut.
2001
~
Mengingat masing masing atom tersebut mempunyai empat buah elektron valensi maka keadaanbebasnya dapatdinyatakanoleh fungsi gelombangIIIi di mana i = 1,2, ,2N. Kemudiankeadaanbebasini diuraikandalam suatukombinasilinier basis {J2I;a}dimanat menyatakan indeks atom dan a, menyatakanindeks orbital.Dengan demikiandapatdituliskanbahwa: N
4
'i';(r) = E E C~a(j)la(r)
(2)
l=la=1
di mana, C famerupakanbobot dari fungsi gelombang denganindeksi untuk atomdenganindekst, pada orbital a. Adapun energi total sistemelektrondan atom untuk sistem silikon U({~}{ Cj}), dapat dibangun dengan menggunakanmodel Tomanek dan Schluter [8], yang untuk kepentingandinamika molekul perlu dilakukan sedikit modifikasi yaitu dengan memasukkanfaktor berupasebuahfungsipotong (cuto.fJ), sebagaiberikut: u = co + E] -E3
'3)
di mana Ebs,adalah energi struktur pita (band structure energy) yang diturunkan berdasarkan parameterisasi tetangga terdekat atom terkait sebagaimana diusulkan oleh Chadi [9], Adapun elemen matriks < VJi',aIHI 'fila] untuk orbital orbital pada atom atom I', clan I, yang terpisah oleh jarak sebesar R. akan mengandung fungsi potong yang disebutkan di atas, Secara lebih terinci fungsi energi tersebut di atas dapat dituliskan sebagai berikut
N t-l
~ = EE
c}'a' c}a(~l'a/IHI~la)
-N Esi
(4) C}'a'c}a (fa'
l'a/
(5)
olehpersamaan persamaansebagaiberikut Er(R} = E~?;(R}f(lRi
-Rill
'- Rc} -E~,S(R)
(7)
N i-I
f(lRi-
Rl'l-
Rc}
(8)
l=2 l'
i l/ a' ta 2N i
Er aRt-Rill)
Fungsi Er daD nbyang muncul dalam E2 dan E3 diberikan
-N Esi
i 2N
= 2 LLL
H~- K.~
t=2t'=1 Sedangkan E3, dimasukkan sebagai salah satu faktor dari energi total sistem dengan alasan agar sistem tidak terpaku pacta satu kemungkinan struktur atom saja.Dengan adanya faktor ini maka sistem dapat mengantisipasi perubahan energi akibatadanya perubahanjumlah ikatan kimia nb. E3, ini diberikan oleh persarnaan
nb = ~~
= 2 LLL
HM..l (S4.4P~)
potong yang biasanya dapat diarnbil sarna dengan panjang ikatan atom (bond length). Sedangkan EoSi, adalah energi sebuah atom netral (terisolasi) yang besarnya sarna dengan 2(Es+Ep), dimana Es=[f2JlsIHIf2Jls] dan Es=[QJixyz IHIQJ/xyz]. Faktor energi total lainnya, yaitu E2, adalah energi tolak menolak yang diperlukan untuk menetralisir pengaruh penghitungan dua kali (double counting) dari interaksi Coulomb yang terdapat dalarn suku energi ikatan kuat. Energi E2, ini juga mengandung kontribusi energi korelasi (exchange correlation) dan interaksi antar atom. Secara terinci energi ini diberikan oleh
2N
Ebs =2 L(l/JiIHll/Ji)
~
IHlla) -N Es;
Adapun parameter parameter e/. f2fzdon ej berturut-turut mempunyai harga 0,225, eV, 1,945eV, clan -1,03eV. Sedangkan harga harga parameter lainnya secara lebih terinci dapat dilihat dalam Khan[ 11].
la
Adapun elemen matriks [1'a1Hlla] dapat ditentukan dengan memperhatikan kedua kasus berikut ini. Pertarna, untuk kasus di mana atom I sarna dengan atom r, maka dapat digunakan persarnaanberikut ini
(ta' IHlla)= { ~
if a = 1, if a = 2,3, 4, if a:f:a'
Kedua, untuk kasus di mana I * r elemen matriks [I'a1Hlla] bagi orbital orbital pacta atom yang berbeda diasumsikan mengandung fungsi potong .f(IRr ~~ -Rc)
sehinggadiperoleh
dengan Ro=2,35A, adalah jarak antar tetangga terdekat di dalam struktur intan silikon. Adapun Rc, adalah jarak
~,
Energi Total (Ab initio) Perbedaan utama antara metoda (ab initio) dan metoda semi empiris (seperti metoda parameterisasi ikatan kuat yang telah kita bahas sebelumnya) adalah bahwa metoda yang pertama tidak memerlukan data data empiris sedangkan metoda memerlukan data data tersebut. Dengan demikian energi total di sini lebih bersifat generik. Meskipun demikian perlu ditekankan di sini bahwa sifat generik ini tidak berarti bahwa metoda ini sepenuhnya analitis dan eksak karena pada dasamya selalu ada aproksimasi yang harus dilakukan juga dalam usaha mencari solusi sistem atom dan elektron yang amat kompleks ini. Aproksimasi yang dilakukan itu khususnya adalah bahwa sistem atom dan elektron dapat dilihat sebagai problema elektron semata-mata (BornOppenheimer approximations) mengingat bahwa massa atom jauh lebih besar dibanding massa elektron. Di samping itu juga sering digunakan aproksimasi kerapatan lokal (local density approximations) untuk menentukan energi korelasi dan exchange.
~J~ 2001
~
Oi dalammetoda(ab inito) ini energipotensial pada umumnya dihitung secara konsisten (self consistently) dengan menggunakan rapat elektron mengikuti suatu skema teori fungsional rapat elektron (density functional theory) atau dft. Namun ada juga sejumlah metoda yang secara langsungmenyelesaikan persamaanSch"odinger dalam suatu medan potensial yang secaracerdasdiusulkan sehinggadiperoleh basil basil yang dapatdipertanggung-jawabkan. Metoda yang kedua ini biasanya sangat berguna untuk sistem nonperiodik. Oalam makalah ini hanya akan dibahas metodapertamayaitu metodayang berbasispada teori fungsionalrapatelektron. Persamaan-persamaan yang biasa digunakan dalamformalisma dftadalahpersamaanKohn-Sham[12] untuk satu himpunankeadaanefektif elektron elektron Pada dasarnyaini menyatakanorbital orbital elektron yang memenuhikaidahorthonormalitas
f1jJ*cr)1jJj(r)dV={ ~
:;~
(9)
daDenergipotensialtotal dapatdinyatakanoleh U[{Iji(r)}] = ~k;Jljij(r)(-~V2Ijii(r))dV+JUnuc(r)n(r)dV
(10)
I
+ ~j;(r)n(r)dV+ jlxc(n(r))dV + Unuc-nuc di mana fnuc-nuc, adalah suatu fungsi yang diketahui bentuknya, Unuc adalah medan potensial yang ditimbulkan oleh suatu inti atom yang diberikan oleh 2~ Zl Unuc(r):; -kc qe L., R l
(11)
l
di mana kc, adalah konstanta Coulomb, Rt, adalah jarak dari titik r ke inti atom I clan qc, adalah muatan elektron.Unuc-nuc merupakan energi interaksi statik antar inti atom yang dinyatakan oleh
1 -2 ~ ZlZl' Unuc-nuc(r)= 2kC!leL., ~ l~l'
(12)
U
dengan Z/, menyatakan jumlah atom dari spesies I. Sedangkan kj menyatakan jumlah elektron yang mengisi orbital i (biasanya sarnadengan dua) clan
sedangkan0(r ), menyatakanpotensialyang ditimbulkan olehmuatansebuahelektrondi r, yang dapatditentukan dengan menyelesaikan persamaan Poisson sebagai berikut
V2tf>(r)= -41\"~ n(r) Oari persamaan (10) terlihat dengan jelas bahwa energi total merupakan suatu fungsional dari orbital orbital elektron. Pilihan terhadap orbital orbital yang benar dilakukan dengan cara melakukan minimisasi energi total ini dengan memperhatikan syarat orthonormalitas yang
4
~,
~
H..ilII..l (~
P~)
H~- K.~
diberikan oleh persamaan(9). Secara garis besar ada dua cara untuk memperoleh energi' total minimum itu. Pertama, menghitungnya langsung dengan menganggapnya sebagai problema minimisasi numerik daD menyelesaikannya dengan menggunakan teknik teknik minimisasi numerik berkendala. Kedua, dengan mencari persamaan persamaan pengali Lagrange bagi suatu minimisasi berkendala dan selanjutnya menggunakanmetoda numerik untuk mencari solusi dari persamaan persamaan tersebut. Persamaan pengali Lagrange untuk teori fungsional rapat elektron inilah yang dikenal luas sebagai persamaanKohn-Sham. Dapat ditunjukkan bahwa persamaan ini mempunyai bentuk sebagaimanahalnya persamaanSchrodinger yaitu
Ii
-2m V2'I/Ji(r) + U(r)'l/Ji(r) = Eit/Ji(r) di mana bentuk potensial U(r) berikut
(13)
, didetinisikan sebagai
U(r) = Unuc(r) + q,(r) + f:.:cc(n(r)) (14) daD dalam hal ini &t= At/ ki. Mengingat bahwa pengali Lagrange di awal perhitungan merupakan konstanta yang tidak diketahui besarnya maka di sini digunakan &j, yang dapat dilihat sebagai energi untuk setiap orbital. Strategi umum yang digunakan untuk mencari energi total biasanya adalah sebagaiberikut. Untuk dapat menghitung semua suku daripersamaan (10) maka diperlukan orbital orbital {VIi (r)}yang tepat daD ini dapat diperoleh dengan menyelesaikan persamaan (13) Untuk menentukan U(r), berdasarkan (14) maka semua komponen persamaan ini hams dapat ditentukan terlebih dahulu. Harga minimum U[{(\j/,{r)}], dapat diperolehjika semua persamaan persamaan itu dipenuhi secara simultan. Khususnya rapat elektron hams konsisten (selfconsistent): harga rapat elektron masukan hams menghasilkan U(r), sedemikian rupa sehingga diperoleh sebuah himpunan orbital orbital \j/,{r), yang memberikansebuah rapat elektron sarna dengan rapat elektron masukannya. Makalah ini tidak akan membahas lebih rinci lagi mengenai kedua metoda tersebut karena kedua metoda itu sudah merupakan sesuatuyang bersifat baku daD telah cukup banyak buku acuan yang membahasnyasecara rinci (Iihat misalnya di [13]). Sebelum menutup bab ini perlu disampaikan di sini mengenai masalah saWall saWall atomis. Hal ini menjadi amat penting mengingat kesalahan kecil dalam komputasi bisa jadi memerlukan waktu berjam-jam untuk mencari penyebabnya. Salah satu cara untuk menghindari hal seperti ini adalah dengan menggunakan analisis dimensional (Iihat misalnya di [14]). Dalam hal teori fungsional rapat elektron dapat dilihat ada empatbuah konstanta fisis: konstanta Planck II massa elektron m, konstanta Coulomb Kc, daD muatan elektron q.. Oleh karena itu dari pada menggunakan saWall satuan standar meter, kilogram, daD detik dapat digunakan saWall satuan L. M, daD T, untuk saWall fundamental
6J~ 2001
~
~
M..it/G..l (5tt t
di manadenganmenggunakanpersamaan persamaan(4), (5) daan(6) dapatditunjukkanbahwa
panjang, massa, clan waktu. Ketiga alternatif satuan ini sangat membantu mengingat bahwa konstanta konstanta fisis itu selalu muncul dalarn dua bentuk kombinasi saja yaitu h2/m, atau Kclle2,yang selanjutnya dapat direduksi menjadi sarna dengan satujika digunakan satuan tertentu.
8{t'a'IHjla) 8B.t
-t'
E
0'
,ta
.!!lL~+J(Ru-).!!lL~[Et'o',la] R~ 8B.t
Jadi J(Rti') Et'a',to~
')
Ji-
=
m
lEL2
kcq;
di mana E, menyatakanenergisatuanyang dalam hal ini E=IML2/r. Dari persamaan(15) dapatdenganmudah diperoleh solusi untuk L dan E, sebagaisatuanatomis standar untuk panjang dan energi. Jadi untuk tujuan simulasidinamika molekul denganmenggunakanmodel energi (ab initio ), dan juga untuk energi ikatan kuat, biasanyadigunakansatuanpanjangdanenergi. L=-1;1/m 1 baht= -,
=O,52917725 A
E=1 hattt..=-:,-- / ,
(k qlj1 1;- m
~
8Rt
(17)
aU(Ru/}
(18)
aRJ;
di mana ~r= I~ -Rrl menyatakanjarak antara atom / dan /'. Untuk sistem Lennard-Jones seperti yang diberikan oleh persamaan (1) dengan mudah dapat kita turunkan persamaangerak atom I, adalah -N
28'3
,,"
MI~=24.}::: 1'..1["Dir-m] ..", /I' [a.(Rt.-R".J+i,(R"-R,.,J+i.(R,.-R,,.J] di mana posisi atom I, dinyatakan oleh vektor Rl = a,. R1x +Iiy R1y + ~ R1z
Adapun untuk sistem yang menggunakan energi total dengan parameterisasi ikatan kuat seperti yang diberikan oleh persamaan (3) maka persamaan geraknya akan diberikan oleh
Rtl'R"
aEbs --a-R::;"
aEl -a"R;
[
] p(Ru./R"jc,-l)_q(Rtl,/R"j(l-l)
"I
(B-Ruo/R")"
{(Ru,/R-jP-(Rtl'/R")'}{,-
(
"IV
j ~l)}](Rt'-Rtj
n 'n.\(~l) B-Ru'/R"
ikatan kuat terparameterisasi. Persamaan gerak untuk sistem yang menggunakan energi total ab initio dicari dengan menggunak!\n metoda Car-Parrinello [15]. Metoda ini menawarkan perhitungan yang paling murah dari segi waktu komputasi. Oleh karena itu metoda ini banyak digunakan untuk simulasi dinamika molekul ab initio. Pacta dasarnya metoda ini memperlakukan derajat kebebasan (degree of freedom) elektron sebagai suatu massa semu (fictitous mass) daD dengan demikian mempunyai dinamika tersendiri. Jadi orbital elektron daD posisi atom divariasikan secara bersama-sama untuk menentukan energi minimum. Dalam formalisma ini sistem akan mempunyai dua jenis massa dengan dinamikanya masing masing. Lagrangian dari sistem tidak saja dapat digunakan untuk menentukan energi minimum sistem namun juga dapat digunakan untuk keperluan simulasi dinamika molekul. Jika massa semu elektron adalah J.l, daD parameter pengali Lagrange dinyatakan oleh A, maka persamaan gerak dari massa semu elektron akan diberikan oleh
aU[{Ri},{\II.}]+ L'" J"'l1k = -I3'IIk
aE) a~
Ablr\ll",
(20)
Sedangkan persamaan gerak dari atom atom diberikan oleh persamaansebagai berikut
MtRt = _aU[{~},{\IIJ}1 Mt~ ..+- =
8 ["R;i: ~ ]
persamaangerak seandainyadigunakan model energi
=27,211396.V
Gaya yang bekerja pacta sebuah atom akibat suatu distribusi energi U(R), pacta dasarnya diturunkan dari negatif gradien distribusi energi ini. Dengan demikian untuk suatu sistem yang terisolasi, yaitu sistem denganjumlah atom N, dan volume sistem V, yang tetap, dengan distribusi energi U(R), akan mempunyai persamaangerak untuk atom /, sebagai berikut
l~t
8Rt
Sedangkan turunan dari energi struktur pita untuk molekul diatomik silikon dapat diperoleh dengan cara sarna. Dengan demikian lengkaplah sudah seluruh perhitungan yang diperlukan untuk menurunkan
PERSAMAAN GERAK
MtRt = -f
[ = --.!!:.-exp +
(16)
k,q;
R:',
Jika digunakan energi total molekul diatomik silikon daTi Raghavachari [10] maka dapat diperoleh
(15)
= lEt
')
P~)
H~- K.. ~
aRt
(19)
+ L I\,rt,.r a {\IIkl\ll"r}} kNr
(21)
aRt
Sebagaimanahalnya denganenergi ikatan kuat, dalam komputasinya, orbital elektron biasanya diuraikan
~I
6J~ 2001
5
~
sebagai kombinasi linier suatu basis. Dengan demikian yang tidak diketahui dari persarnaan persarnaan diatas adalah koordinat atom (Rt) koordinat massa elektron semu yang dalarn hal ini dinyatakan oleh koefisien koefisien basis dalarn kombinasi linier itu, dan konstanta konstanta pengali Lagrange. Perhatikan bahwa matriks Akm merupakan matriks simetri yang besarnya sarna
dengan setengah keadaan yang diduduki (occupied states) oleh elektron elektron dalarn sistem. Penghitungan gaya merupakan proses yang paling menyita waktu dalam simulasi dinarnika molekul mengingat sumasi harns dilakukan untuk setiap atom. Meskipun demikian ada sejumlah cara untuk mengurangi konsumsi waktu tersebut. Salah satu di antaranya adalah memanfaatkan hukum ke tiga Newton, yaitu aksi sarna dengan negatifreaksi atau dengan kata lain FI/'= -FI/' .Di sarnping itu penggunaan fungsi potong, daftar tetangga, dan untuk kasus kasus yang memungkinkan
digunakannya pembagian ruang simulator menjadi beberapa sel simulator juga dapat mereduksi konsumsi waktu dalarn komputasi persarnaangerak ini.
SYARA T BATAS PERIODIK Oinarnika molekul biasanya dipakai untuk sistem sistemyang terdiri atas beberapapuluh sarnpai beberaparibu atomdenganelektronelektronnya.Sistem yang sedemikian kecil akan didominasi oleh efek permukaan,yaitu interaksi antara atom dan dinding simulator. Untuk simulasi di mana efek permukaanini tidak menjadi perhatian utama maka dapat digunakan syaratbatasperiodiksbp. Oalarnmenggunakansbp untuk mensimulasikan sistemyang terdiri atas N, buah atomyang terkurungdi dalarn simulatordenganvolume V, diasumsikan bahwa volume V, ini hanya merupakan bagian kecil dari padatanyang menjadiobjek studioVolume V, ini disebut sel utama; diasumsikanbahwa sel utama ini mewakili objek materi yang dibentuk oleh sel utama dan replika replikanya.Replika ini merupakansetcitra ( imagecell ), yang mempunyaiukuran dan bentuk yang tepat sarna dengansel utama. Jadi sel utama itu diasumsikansecara periodik diulang ke seluruh arab sehingga terbentuk sarnpelobjekmakroskopisnya. Ada dua konsekuensidari syaratbatasperiodik ini. Pertama,sebuah atom yang meninggalkanruang simulator melalui satu dinding akan digantikan oleh sebuahatomcitranya dari dinding dihadapandinding di manaatomitu keluar. Kedua,jarak antar atomditentukan oleh jarak terdekat, baik antar sesarnaatom dalarn sel utamaatau denganatomcitra di sel citra. Efek ini harus diperhatikan dalarn menghitung gaya dan integrasi persarnaanpersarnaangeraknya. Setiap tahap integrasi numerik harus diikuti dengan mengecek kembali koordinatkoordinatyang diperolehnya.Jika sebuahatom temyatadiketahuitelah keluar dari ruang simulatormaka koordinatnyaharusdisesuaikan kembali.
6
~,
/)~:-~.,.~HM..i (.s«t.4P~)
H~- K.~
Misal koordinat x dari atom atom itu didetinisikan terletak dalam selang -L/2, dan Lx/2 di mana Lx merupakan panjang ruang simulator dalam arab x, maka tes yang harns dilakukan adalah: jika Rtx ~ L/2, maka ganti R/xdengan Rtx -Lx, dan jika Rtx < -L/2 , maka ganti R/x dengan R/x + Lx. Efek syarat batas dalam menentukan jarak pasangan atom dapat ditentukan dengan tes: jika RI/' > L/2 maka ganti R/xdengan RI/'- Lx, seandainya Rtx > 0 atau ganti dengan Rtx + Lx, seandainya R/x < O. Syarat batas periodik paling mudah ditangani untuk ruang simulator yang berbentuk kubus. Bentuk ruang simulator sembarangpunsebenarnyatidak menjadi hala!1gan bagi diberlakukannya syarat batas periodik namun tentunya efek sbp menjadi tidak sesederhanajika bentuknya kubus. Motivasi utama untuk menggunakan ruang simulator yang non-kubus adalah agar diperoleh rasio volume terhadap permukaan sebesar mungkin sehingga memperbesar jarak maksimum yang bisa dicapai sebelum pengaruh periodisitas kemudian mensyaratkanuntuk menggunakanjarak terkecil. Alasan lain dalam menggunakan bentuk yang lebih kompleks adalah agar dapat digunakan untuk mensimulasikan struktur kristal yang mempunyai sumbu sumbu nonorthogonal.
DAFTARTETANGGA Bagian yang paling menyita waktu dalarn sebuah simulasi dinarnika molekul adalah menghitung gaya gaya yang bekerja pacta sebuah atom. Sebagai contoh jika di dalarn sistem terdapat N, buah atom maka gaya yang bekerja pacta atom atom pacta dasarnya akan melibatkan sebanyak y~ (N-I), buah pasangan atom. Konsekuensinya jika seluruh pasangan atom dilibatkan dalarn menghitung gaya maka jumlah perhitungan akan sebanding dengan kuadrat jumlah atom. Narnun jika digunakan fungsi potong (cut off) sedemikian rupa sehingga gaya itu besarnya akan sarna dengan nol jika jarak pasangan atom itu melebihi suatu harga maka waktu yang dibutuhkan untuk menghitung gaya menjadi
berkurang. Salah satu cara yang banyak untuk mereduksi Waktu komputasi dalarn simulasi dinarnika molekul, khususnya dalarn mengevaluasi jarak antar pasangan atom yang saling berinteraksi, adalah dengan menggunakan daftar tetangga. Untuk setiap atom l, metoda ini menawarkan sebuah daftar atom atom yang berada dalam jarak interaksi RL, yang umumnya sedikit lebih besar dari jarak potong Rc. Untuk sistem LennardJones biasanya digunakan RL=Rc+O,30cr. Dengan demikian daftar ini mengidentifikasi atom atom yang memberi kontribusi dalarn penentuan gaya yang bekerja pacta atom l, tersebut. Daftar ini diperbaharui di setiap atau setelah beberapa langkah integrasi numerik. Sebagai ilustrasi misalnya sistem yang disimulasikan adalah sistem Lennard-Jones dengan rapat jenis pcr3=O,8dan jarang potong sekitar Rc=2,5cr. Untuk
~J~ 2001
~
~
/1..itJ...l /5tt t
P~)
H~- K.~ precision akan membantumemperkecilkesalahanjenis kedua.Dalam bab di baw~ ini akan disajikan secara tiga jenis metoda integrasi numerik yang banyak digunakandalamsimulasidinamikamolekul.
sistem seperti ini setiap atom akan mempunyai tetangga sekitar 75, buah atom dalam radius sekitar R L=2,8cr. Tetapi di dalam daftar tetangga untuk atom ell, hanya disimpan atom atom dengan indek f, di mana f > I, karena untuk I' < I, atom I akan muncul dalam daftar tetangga dari atom f. Jadi secara rata rata, daftar tetangga itu memerlukan sekitar 75/2 7 40, lokasi penyimpanan per atom. Dengan demikian untuk simulasi sistem 256, atom dapat digunakan sebuah l(array), sebut saja misalnya DFTR, dengan sekitar 11000, elemen. Kemudian dengan menggunakan sebuah array lainnya, sebut saja misalnya ISI, untuk menentukan lokasi di dalam DFTR tetangga tetangga daTi suatu atom. Dengall menggunakan dua buah array satu dimensi, bukan sebuah array dua dimensi, dapat dihindari biaya komputasi yang harus dibayar karena menggunakan variabel dengan dua indek. Pembahasanyang lebih rinci dapat dilihat misalnya di [16].
Algoritma Verlet Algoritma inilah yang paling banyak digunakan untuk keperluan dinamika molekul. Ide dasarnya adalah menguraikan posisi atom, misal atom dengan indek I, R[ dalam deret Taylor sampai orde ketiga, baik secara maju ( forward) maupun mundur (backward) dalam waktu. Jadi dapat dituliskan R,(t + .1t) =
R,(t) + ~(t).1t + ~~(t).1t' + i ~(t).1"
R,(t -.1 t) = R,(t) -~(t).1
t + ~~(t).1 t' -i ItI(t).1t3 + 0(.1 t4)
Jika kedua persarnaandi alas dijumlahkan maka akan diperolehbentukdasardari algoritmaVerlet yaitu
ALGORITMA INTEGRASI WAKTU
Rt{t +M) =2 Rt{t)-Rt{t-M)
Algoritrna ini diperlukan untuk mengintegrasi secaranumerik persarnaan persarnaangerak atom atom yang saling berinteraksi.Algoritrna ini didasarkanpada metoda beda hingga (finite difference),dimana waktu didiskritisasipadasebuahgrid terhingga,denganjangkah waktu (time step) At, merupakanjarak antaradua titik berturutanpada grid tersebut. Jika diketahui posisi daD turunanwaktunyapadawaktu t, makaalgoritrnaini akan memberikanbesaranyang sarnapada waktu berikutnya, yaitu t+At. Selanjutnya dengan mengulang (iterasi) prosedur itu maka akan diperoleh evolusi waktu dari sistemyang sedangdisimulasikan. Sudah barang tentu solusi yang diperoleh merupakansuatu solusi pendekatandaD berarti selalu mengandungkesalahanbawaan.Kesalahanitu sekurangkurangnyadapatdiklasifikasikansebagaiitemizeitem. Truncation errors, hal ini mengingatbahwametoda beda hingga biasanyadidasarkanpada uraian deret Taylor yang dipotong sampaipada suku tertentu. Kesalahanjenis ini tidak tergantung pada faktor implementasitetapi merupakankesalahanintrinsik dari algoritmayangdigunakan.
+ ~{t)Atl + 0 (At4) (22)
Mengingat bahwa persamaan di atas diperoleh dari integrasi numerik persamaan gerak Newton maka RAt) = -(I/MU VU ({RAt)}). Seperti terlihat dalam persamaan (22) maka kesalahanjenis pertama dari algoritma ini hila sistem ber-evolusi sebesar Llt adalah dalam order Llt4, meskipun turunan ketiga tidak muncul secara eksplisit. Pada saat yang sarna algoritma ini amat sederhana dan mudah implementasinya sehingga amat populer penggunaannya dalam simulasi dinamika molekul. Persoalan yang timbul dalam menggunakan algoritma Verlet versi ini adalah bahwa kecepatan atom tidak langsung tersedia. Meskipun kecepatan itu tidak diperlukan untuk mengetahui evolusi trayektori namun pengetahuan mengenai kecepatan ini kadang kadang diperlukan, misalnya untuk menghitung energi kinetik yang amat diperlukan untuk menguji aspek konservasi energi total sistem. Kecepatan itu dapat saja dihitung dengan menggunakanpersamaan Rt(t) = Rt(t + ~ t) -Rt(t -~ t)
2~t
Round-off errors, hal ini berkaitan dengan implementasi dari algoritma yang digunakan. Misamya, sampai seberapa banyak angka yang digunakandalamkomputasiaritmatikanya. Keduajenis kesalahandi atas dapat diperkecil dengan memperkecil ~t. Untuk harga M, yang besar maka kesalahanakan didominasi oleh jenis pertama namun kesalahan ini akan mengecil dengan cepat bila ~t diperkecil. Kesalahanjenis keduamengecillebih lambat dibanding dengan kesalahan jenis pertama dengan mengecilnya ~t. Penggunaan komputasi double
~,
+ 0(.1 t4)
Namun kesalahan dalam persamaan untuk kecepatan ini adalah dalam order Llr, bukan lagi dalam order Llt4, sebagaimana halnya persamaan (22). Untuk mengatasi persoalan ini telah dikembangkan beberapa variasi dari algoritma Verlet ini yang salah satu di antaranya adalah Rt(t + ~ t)
=
Rt(t) + Rt(t) ~ t + (1/2) Rt(t) ~ t2
Rt(t + ~ t/2)
;
Rt(t) + (1/2) Rt(t) ~ t
Rt(t+~t)
;
-(l/MvVU(Rt(t+~t))
Rt(t + ~ t)
=
Rt(t + ~ t/2) + (1/2) Rt(t + ~ t) ~ t
Predictor-Corrector Algoritma ini pertama kali digunakan untuk keperluan simulasi dinamika molekul oleh Rahman[17].
6J~ 2001
~
Sedangkan versi yang kemudian sering digunakan biasanya diambil dari metoda yang koleksi oleh Gear [18]. Sesuai dengan namanya metoda ini secara garis besar terdiri atas dua tahap yaitu prediksi dan kemudian disusul dengan koreksi terhadap prediksi itu. Prediksi posisi atom R" pada waktu t+Llt dilakukan dengan menggunakan deret Taylor sampai orde ke lima yang di dasarkan atas pengetahuan mengenai posisi dan turunan turunannya pada waktu t. Dengan demikian turunan turunan R'" R"" R(iii),dan R(vi), diperlukan di setiap tahap dan nilalinya pada waktu t+Llt, dicari dengan menggunakan deret Taylor sebagaiberikut ~('+4t)
(t)~ 41 + R(.i l'( J ~ 51 + R(,o> t t B.,(t+4t)
t 3! ) ~+~ +Rii")(t
R,(t + At)
n~;;')(t+4t) ~ ~iji)(t)+R~;.)(~)4t+R~II)(~)~ R~w)(t+6t) = R}W)(t)+~O)(t)4t R}°)(t+.1t) ~ Ri~)(t) Gaya yang bekerja pada setiap atom pada waktu t, diperoleh dengan menggunakan prediksi koordinat koordinat atom tersebut. Koreksi terhadap prediksi koordinat clan turunan turunannya dilakukan dengan menggunakan adanya perbedaan antara percepatan basil prediksi clan yang diperoleh dari perhitungan gaya. Jika didefinisikan bahwa =
Rt{t
+ &}
-R.~{t
1/6
5/6
~
q=5'
19/120 3/4
1/3
1/2
3/6 251/360 1 11/18
1/12
1/6 1/60
PENGENDALIAN TEMPERA TUR
~ R.,(t)+R}4H)(t) 4£+ R}ill)(t)~~R}V)(t)~
~ == ~A Rt{At}1 2
~
<Xi
<Xo-
as
4! .
P~)
H~- ~.~
Tabel 1. Harga parameter <Xi.dalam algoritma Gear untuk persamaan diferensial orde dua dengan menggunakan prediktor orde q.
~
= Rot(t)+~(t)4l+R}H')(t)~
HJ4d (~
pada orde dari deret Taylor yang digunakan dalam proses prediksi. Gear mendapatkan harga harga parameter ini dengan menggunakan setiap algoritmanya pada persamaan persamaan diferensial linier daD analisis matrik stabilitasnya.
~ ~ ~
= R,r;(t)+~(t)4t...Iit(t)~+R}tU)(t)¥
~
+At}
(23)
Terdapatcukup banyak mekanismetermostatik yang tersediadalamliteratur. Banyak peneliti melakukan pengendalian temperatur dengan cara menyesuaikan harga kecepatansecarateratur clanterns menerusselama simulasi berlangsung. Formulasi metoda dinamika molekul pada temperaturkonstanyang paling lengkap pertamakali diperkenalkanoleh Nose. Ia menunjukkan bahwa metoda metoda terdahulu dapat diturunkan daTi metoda yang diperkenalkannya.Karena alasan inilah maka hanya metoda ini yang akan disampaikandalam makalah ini. Meskipun demikian di sini akan disampaikan metoda Nose [19] yang telah disederhanakan oleh Hoover [20] sehinggakomputasinya menjadilebihmurah. Berdasarkanalgoritma Nose-Hoover maka persamaan gerakperlu dimodiflkasimenjadi
di mana R~{t+t.t} , diperoleh melalui prediksi maka denganmenggunakan algoritma Gear koreksi dilakukan sebagaiberikut
Rt = Rf +~ao ~&
=
Sukutambahanzeta,dot(BR)ell, bertindaksebagaimana halnya gaya gesek. Koefisien zeta, adalah variabel dinamis baru denganevolusinya diatur oleh persamaan sebagaiberikut:
~t\t +~!Xl
Q ~ = L Mt (Rr -gkeT
(25)
t di mana Q, adalah masa dinamis dari tennostat, k B, adalah konstanta Botlzmann, daD T, adalah temperatur atom yang diinginkan, serta g=3N, dengan N, adalah jurnlah atom yang acta di d~lam simulator. Penggerak Adapun harga parameter parameter C1j,dalam algoritma Gear di atas dapat dilihat dalam tabel di bawah ini. Parameter ini rnernelihara stabilitas nurnerik dari algoritrna ini. Harga parameter (Xi,tergantung pada orde dari persamaan diferensial yang akan dipecahkan clan
8
~,
perubahan dari
koefisien ~
adalah
adanya
ketidakseimbangan antara energi kinetik daD harga rata rata gkBT. Umpan balik negatif akan menjamin energi kinetik berosilasi di sekitar harga konstan yang diberikan
6J~ 2001
~
~
M..itIt..l(Stt t P~)
H,-- 1(. ~ oleh gkBT. Jika harga energi kinetik ini lebih besar dari gk BT, turunan terhadap waktu dari ~, berharga positif, atau ~' > 0, dan ~ bertarnbah besar harganya sehingga
tekanan dapat dikendalikan. Sebagairnana dalam pengendalian temperatur di sini pun digunakan sebuah massa semua Mv, yang merepresentasikan massa piston yang bergerak sehingga memungkinkan volume bervariasi besarnya. Namun untuk mengurangi kompleksitas jika harus menganalisis dinamika piston maka diasumsikan bahwa perubahan volume bersifat uniform. Kedalam sistem ditambahkan dua buah derajat kebebasan sehingga persamaan persamaan gerak yang terlibat adalah
rnenjadi positif. Di daerah ~ positif sistern akan rnengalarni pendinginan. Jika harga energi kinetik lebih rendah dari gkBT rnekanisrna umpan balik negatif bekerja sebaliknya; dalarn hal ini ~ akan rnengecil sarnpai rnenjadi negatif, dan di daerah ~ negatif sistern akan rnengalarni pernanasan.Dengan tara seperti inilah energi kinetik akan berosilasi di sekitar harga rata rata gkBT. Jadi energi kinetik rata rata akan rnernpunyai harga yang sarnadengan hasil teori ekuipartisi gkBT. Jika didefmisikan suatu variabel s sedernikian rupa sehingga
.,
~
.?
,=~ s
S
V
sehingga dapat dituliskan persamaandinamis untuk s sebagaiberikut
maka Hamilton sistem dengan kebebasan itu menjadi
1
= -M;Vl73VU{{~})=
=
S-
-+S
sV -+. S
(s:;+w?V ) R.t
G IS Ms
(31)
G')Sl 3M\,V
di mana
Gl
=
Gl
= mV1/3L Ri+Vl/3 L Rtt .Fu'- 3PV
tambahan derajat
Mt.yl/3I:
.
Rl-gkBT
.l
di mana RD' = R, -R', dan FI/' adalah gaya yang bekerja pacta atom I akibat interaksinya dengan atom r. Sedangkan untuk harga tekanan P, dan temperatur T, merupakan harga yang diinginkan dan oleh karena itu merupakan harga awal yang diberikan. Adanya tambahan dua derajat kebebasan maka perataan dinamika sistem fisis ini akan menghasilkan perataan yang berada dalam ensambel (NPT). Hamilton dari sistem dengan tarnbahan dua derajat kebebasanmenjadi
Walaupun Hamilton ini lestari namun tidak mempunyai arti fisis. Mengingat sifat kelestariannya Hamilton ini dapat digunakan untuk mengecekkeabsahandari proses simulasi untuk ensambel (NVT) ini. Algoritma untuk menyelesaikan persamaan (24) secara numerik dapat diturunkan sebagaimana halnya algoritma Verlet. Jika didefinisikan /3= (1/2),l;(t)~t, maka dapat diperoleh evolusi posisi atom sebagai berikut
11= ~MtV1/J 1;: ai+U{{~}}+
~(t+&)= ~
~M,{lls}1 +~M"{V Is}1+tk6nns+PV
(32)
['2RtCt)-(I-I3)~(t-&)+~(,)&?] +0(&' (27) Walaupun Hamilton ini lestari namun tidak mempunyaiarti fisis. Mengingat sifat kelestariannya Hamilton ini dapat digunakan untuk mengecek keabsahandaTiproses simulasiuntukensambel(NPT) ini.
Sedangkan evolusi untuk variabel dinarnis l;; dapat dicari dengan menggunakan persarnaanberikut
+ O(At]) (28)
SIMULASI SISTEM SILIKON PENGENDALIAN TEKANAN DAN TEMPERA TUR
Walaupuncukup banyak metodapengendalian tekanan namun gabungan pengendaliantekanan dan temperaturcukup menarikuntuk disoroti di sini. Metoda ini pada dasarnyaditurunkan dari ensambelisotermalisobarik(NPT). Dalam makalah ini secararingkas akan disampaikanalgoritma (NPT) yang secara lebih rinci dapat dilihat misalnya di [21]. Diasumsikan bahwa volume sistem dapat divariasikanuntuk menjagaagar
~,
Sebagai contoh dari metodologi simulasi dinarnika molekul di sini disajikan simulasi sistem 64atom silikon dengan struktur intan yang ditempatkan dalarn suatu simulator kubus dengan sisi L=lO,48A Dengan demikian kepadatannya adalah 2,59gr/cm3 yaitu tepat sarna dengan kepadatan silikon cair pada temperatur leburnya. Selanjutnya diberlakukan syarat batas periodik sehingga di setiap saat sebuah atom dalarn kubus akan memiliki 26 citra. Kepadatan sistem akan
6J~ 2001
9
~
konstan karena diberlakukannya syarat batas periodik ini. Selang waktu integrasi yang digunakan di sini adalah 1,0 X10-15detik. Energi potensial sistem yang kami gunakan adalah ikatan kuat terparameterisasi sebagaimana diberikan oleh persamaan (3). Untuk mendapatkan energi minimum sekaligus menurunkan persamaan gerak kami menggunakan formalisme Carr-Parrinello. Dengan demikian pada dasarnya digunakan persamaan persamaan (20) clan (21) hanya energi total kami turunkan dari ikatan kuat terparameterisasi. Juga kami sajikan di sini penggunaan mekanisme termostatika Nose-Hoover baik pada temperatur rendah maupun pada temperatur tinggi, yaitu temperatur lebur silikon.
~
H..lt/...l(Slt.4 P~)
H~- K.~
dimampatkan. Mengembangnya kaki kaki fungsi delta
disebabkanoleh karena adanyapergeseranacak posisi atom dari posisi keseimbangannya. Posisi posisi puncak pertama dan kedua menunjukkan jarak tetangga tetangga terdekat pertama dan kedua, yaitu si 2,27Adan 3,72,A. Sebagai pembanding kami sampaikan disini bahwa jarak tetangga terdekat pertama dan kedua dari kristal silikon dengan struktur intan (normal dan tidak dimampatkan) adalah 2,35Adan 3,84,A,. Dalam studi ini kami memilih untuk menggunakan fungsi pasangan korelasi g(r), dibanding dengan faktor struktur statis S(k), mengingat terbatasnya ukuran sel yang kami gunakan. Tabel 2 menunjukkan sejumlah parameter penting basil studi mengenai silikon pada rasa cair dari para peneliti
terdahulu. 1500 I
'
1000 ~ '.:' '5\
500
nL-0
-2
4
a
-6
jarak r(A)
Gambar 1. Korelasi pasangan korelasi (pair correlation) untuk sistem 64 atom silikon yang dimampatkan (compressed). Posisi puncak-puncak pertama dan kedua berturut-turut menunjukkan jarak tetangga terdekat pertama dan kedua O'Oe+OO
Konfigurasi atom atom pada rasa cair sangatlah berbeda dengan konfigurasi atom atom padatan. Oleh karena itu perlu digunakan suatu kuantitas yang mampu memberikan deskripsi kuantitatif mengenai distribusi jarak antar atom. Besaran ini dapat digunakan sebagai ukuran konfigurasi atom atom dalam rasa cairo Untuk tujuan ini di sini digunakan fungsi distribusi pasangan g(r), sebagaimanahalnya [] yang didefinisikan sebagai
~ -6.00-06
I. .
1.2e.05
-2.4e-Q5
V n(r)Ar ] g(r) = N [47fr2
.-05
'1.Oe...oo
dimana n(r), adalah jumlah atom yang berada pada jarak antara r daD r+Ar dari suatu atom yang dimaksud, N, adalah jumlah atom di dalam volume V. Persamaan di atas pada dasarnya juga menyatakan bahwa n(r) sarna dengan jumlah atom yang berada di dalam elemen bola konsentris dengan jari jari r, daD r+L\r dengan atom yang dimaksud sebagai titik pusatnya. Dalam prakteknya perhitungan numerik dilakukan dengan cara perataan dari semua atom di dalam kubus simulator sebagai titik pusatnya. Gambar
korelasi untuk
10
1.4e-12
2.1e-12
2.~lJ -" ~-
3.5e-lJ
waktu (detik)
1. menunjukkan fungsi pasangan sistem 64 atom silikon yang
~I
7.0e-1J
Gambar 3. Evolusi ketelitian; garis tebal merupakan evolusi ketelitian jika sistem mengalami quenching satu kali di awal simulasi dengan melakukan diagonalisasi matrik Hamilton, sedangkan garis terputus adalah evolusi jika sistem mengalami quenching sebanyak due kali dengan menggunakan metoda steepest descent di awal simulasi . Evolusi temperatur elektron elektron semu. Garis tebal merupakan evolusi jika sistem mengalami quenching satu kali dengan melakukan diagonalisasi matrik Hamilton. Garis terputus merupakan evolusi temperatur jika sister:n mengalami quenching sebanyak
6J~ 2001
~
~
/1..id..l(~
P~)
H~- K.~
parameterQ=2850, daD didapatkan temperatur sistem berosilasidi sekitar temperaturyang diinginkan dengan amplitudo yang wajar. Potensial sistem dimana temperatur dikendalikan menggunakan mekanisme termostat melalui penyesuaian kembali kecepatan ditumpangtindihkansebagaiperbandingan.
dua kali dengan menggunakan metoda steepest descent di awal simulasi Gambar 2 daD 3 menunjukkan temperatur dan ketelitian elektron elektron semujika sistem mengalami quenching satu kali di awal simulasi dengan menggunakan diagonalisasi langsung matrik Hamilton. Hasil basil proses sejenisjika dilakukan quenching sampai dua kali ditampilkan secara bersama sarna untuk membandingkan kinerja kedua metoda tersebut. Di sini kami tidak menyajikan evolusi potensial daD temperatur ion ion untuk kedua kasus itu karena harga harganya nyaris sarna.
2200
2000
g 5 ~
190
~ E
1!.
184 .\
.,
.c
.." .'
g ~
\
I .. ..
".I
/,' ,':"
1400
..
l~~+OO
8-
E
~
..' I' I 166 160I O.Oe+OO
., 4.De-13 .waktu
.; " ..
", ,
.S.oOe:13
7.sOe:13
1.00e:12
1.2ie-12
Gambar 5. Temperatur sistem versus waktu dengan termostat diatur pads T=17000K dan Q=2400. Posisi awal atom atom diambil dari keadaan cair yang amat panas
., :'
I
.2.sOe:13
waktu (detik)
" " .-, « " .' ..
.. B.De-13
1.2e-12 (detik)
1.6e-12
2.De-12
Gambar 4. Temperatur atom sebagai fungsi waktu. Garis tebal berasal dari Lagrangean semu sedangkan garis terputus-putus dari Lagrangean fisis. Dalam kedua kasus pengamatan ini temperatur sistem dikendalikan dengan menggunakan mekanisme Nose-Hoover
Temperaturatom atomadalahsekitar 175°Kyang berarti jauh di atas temperaturelektron semutersebut.Dengan demikiandapatdisimpulkandi sini bahwapadadasarnya melakukan proses quenching satu kali dengan diagonalisasi matrik Hamilton adalah sarna baiknya denganmelakukanproses quenchingdua kali dengan metoda minimisasi steepestdescent.Oleh karena itu dalam penelitian yang kami sajikan di sini telah digunakanmetodadiagonalisasimatrik Hamilton di awal simulasisistem. Dalam simulasipada temperaturkonstan175°K, tersebutkami menggunakanparameterQ=205, dengan satuan energi x waktu2,dan titik potong Rc=2,80A. Gambar4 menunjukkantemperatursistemsebagaifungsi waktu untuk suatu simulasi selama2 pikodetik. Dalam gambar tersebutditunjukkan keadaantemperatursistem untuk kedua metoda simulasi, yaitu baik untuk Lagrangeansemumaupun Lagrangeanfisis (minimisasi energi sistem dilakukan di setiap langkah integrasi. Sedangkandi gambar5 kami tunjukkan potensialatom atom sebagaifungsiwaktu. Dari keduagambarini dapat kita lihat bahwameskipuntemperatursistemLagrangean semu lebih dekat ke temperatur yang diinginkan dibandingkan Lagrangean fisis namun potensial keduanyanyarissarna. Gambar5 menunjukkansistempadatemperatur leburnya. Untuk simulasi ini kami menggunakan
~,
Gambar 6 menunjukkan mean square displacement ~ msd yang dihitung dengan menggunakan posisi atom R (t) dalam selang waktu t, sebagai berikut
(R2(t») =
di mana 't menyatakanhimpunan keadaan awal daD tanda kurung miring menyatakan perataan. Selanjutnya koefisien difusi D, dihitung dengan menggunakan bagian linier dari msd ini. Simulasi pada gambar 6 berlangsung selama 0,15, pikodetik. Sedangkan gambar 7 menunjukkan korelasi pasangan g(r),terhadap jarak r. Ketelitian darisimulasi ini ditunjukkan oleh gambar 8.
6J~ 2001
l "
~
~
H..ltiL..l 15tt...tP~)
H~~ Harga Q, yang terlalu kecil akan menyebabkan tiadanya interaksi antara massa sumberpanas dan massaatom, masing masing mencapaikeadaanseimbangnyasendiri sendiri. Untuk harga Q=2400, kami mencapaiketelitian dalamorde 10.3. 3.0
2.0 -;:-
'51 .0
0.0 ...J! 0.0
1.2
.I 2.4
..I 3.6
4.8
6.0
jarak r(A)
Gambar 7. Korelasi pasangan g(r), versus jarak r, pada temperatur lebur sistem.Simulasi berlangsung dalam waktu 0,15 pikodetik. Bulatan hitam berasal dari difraksi sinar X [22] sedangkan garis vertikal terputus-putus
merupakan lokasi titik potong
KESIMPULAN Dalam makalah ini telah dikemukakansecara cukup rinci mengenai metoda simulasi dinamika molekul. Kerincian dari makalah ini dimaksudkanagar dapat digunakan sebagai sumber bagi para pemula dalam menggunakanmetoda ini untuk keperluannya masingmasing. Sebagaicontoh aplikasisederhanatelah disajikan basil simulasi sistem silikon dengan menggunakanenergi ikatan kuat dan formalisma CarrParrinello.Sengajadipilih sistemini karenatelah banyak basil basil mengenaisistem ini dari peniliti terdahulu yang menggunakanberbagai metoda lainnya. Dengan demikian sistem silikon ini dapat digunakan untuk mengujikeabsahansimulasiyang telah kami lakukanitu. Walaupunsimulasiyang kami lakukan bersifat sangat dini dan sederhananamunhasilnyatemyatatidak banyak berbedadenganbasilpara penelititerdahulu. PENGHARGAAN Kami ingin menyampaikan ucapan terima kasih kepada Menteri Riset dan Teknologi Republik Indonesia yang telah memberikan dukungan keuangan bagi penelitian ini melalui RUT V 1997-1999. Tanpa
12
~I
1(.~
dukungan tersebut penelitian dan hasilnya tidak dapat kami sajikan dan laporkan di sini. '
DAFTARPUSTAKA [1]. B.J. ALDER AND T.E. WAINWRIGHT, J. Chern.Phys. 27,1208 (1957) [2]. T.E. KAP.AKASIDIS AND M. MEYER, Modeling Simul. in Mat. Sci. and Eng. (8 ),117 (2000) [3]. O.M. CABARCOS ADN J.M. LISY, Inti. Journal of Mass Spectroscopy(185 186187),883 (1999) [4]. S. PALUCHA AND Z. GBURSKI, J. Mol. Structure(480-481),343 (1999) [5]. YOUNG-KYUN KWON, D. TOMANEK, ADN S. IIJIMA, Phys.Rev.Lett. (82), 1470(1999) [6]. M. J. MEHL AND D. A. PAPACON STANTOPOULOS,Phys. Rev. B (54), 4519 (1996) [7]. L. VERLET, Phys.Rev. (159),98 (1967) [8]. TOMANEK AND SCHL"UTER,Phys.Rev. Lett. (56), 1055(1986) [9]. D.J. CHADI, Phys.Rev.B (29),785 (1984) 9 [10]. K. RAGHAVACHARI, J. Chern. Phys. (83), 3520(1985 ) [11]. F.S. KHAN AND J.Q. BROUGHTON,Phys.Rev. (39 ), 3688(1989 ) [12]. W. KOHN AND L.J. SHAM, Phys.Rev. (140), Al133 (1965) [13]. J.M. THIJSSEN, Computational Physics, CambridgeUniversityPress,Cambridge1999 [14]. TOM'AS A. ARIAS, Notes on the (ab initio) theory of moleculesandsolids: Densityfunctional theory, Physics 680 class notes, Departmentof Physics,CornellUniversity,Cornell2000 [15]. R. CAR AND M. PARRINELLO, Phys. Rev. Lett. (55),2471 (1985 ) [16]. J.M. HAILE, Molecular DynamicsSimulation: Elementary Methods, John Wiley, New York, 1992 16
[17]. RAHMAN, Phys.Rev. (136) (2A), 405 (1964) [18]. C.W. GEAR, Numerical Initial Value Problems in Ordinary Differential Equations,Prentice-Hall, EnglewoodCliffs, NJ 1971 [19]. S. NOS'E, Molec. Phys. (52),255 (1984) (20]. W.G. HOOVER,Phys.Rev. A (31 ),1695 (1985) [21]. D.C. RAPAPORT, The Art of Molecular Dynamics Simulation, Cambridge University Press,Cambridge1995 [22]. STICH, R. CAR, M. PARRINELLO, Phys.Rev. Lett. 63,2240 (1989 )
6 J~ 2001
Ke Daftar Isi
