Laporan Tugas Akhir
PERANCANGAN DAN ANALISA SIMULASI DINAMIKA MOLEKUL ENSEMBLE MIKROKANONIKAL DAN KANONIKAL DENGAN POTENSIAL LENNARD JONES Diajukan untuk memenuhi syarat kelulusan tahap sarjana di departemen Teknik Fisika Institut Teknologi Bandung
Oleh:
AREE WITOELAR 13396018
Pembimbing: Hermawan K. Dipojono, Dr., Ir. Nugraha, Dr., Ir.
DEPARTEMEN TEKNIK FISIKA FAKULTAS TEKNOLOGI INDUSTRI INSTITUT TEKNOLOGI BANDUNG 2002
i
PERANCANGAN DAN ANALISA SIMULASI DINAMIKA MOLEKUL ENSEMBLE MIKROKANONIKAL DAN KANONIKAL DENGAN POTENSIAL LENNARD JONES
Oleh:
AREE WITOELAR 13396018
Disahkan oleh:
Pembimbing I
Pembimbing II
Hermawan K. Dipojono, Dr., Ir.
Nugraha, Dr., Ir.
DEPARTEMEN TEKNIK FISIKA FAKULTAS TEKNOLOGI INDUSTRI INSTITUT TEKNOLOGI BANDUNG 2002
ii
1.1
ABSTRAK
Perancangan material baru dengan sifat-sifat tertentu sangat berguna untuk berbagai bidang. Namun perancangan material dengan eksperimen membutuhkan biaya mahal. Salah satu solusi adalah dengan melakukan eksperimen komputer dengan teknik dinamika molekul. Dinamika Molekul adalah teknik simulasi komputer yang mengamati pergerakan molekul-molekul yang saling berinteraksi. Dari molekul-molekul ini, dapat diperoleh sifat material yang diamati. Dinamika molekul banyak digunakan dalam berbagai bidang ilmu untuk mempelajari material. Berdasarkan itu, maka dirancang suatu simulasi dinamika molekul untuk sistem terisolasi (ensemble mikrokanonikal) dan sistem dengan pengendalian temperatur (ensemble kanonikal). Untuk menguji kebenaran simulasi, maka digunakan model potensial Lennard-Jones dan dilakukan analisa hasil simulasi dengan berbagai parameter. Perangkat lunak simulasi dikembangkan dengan bahasa pemrograman C++ dengan berorientasi obyek. Data yang diperoleh menunjukkan bahwa simulasi dinamika molekul untuk ensemble mikrokanonikal dan kanonikal memenuhi spesifikasi ketelitian yang dibutuhkan. Simulasi dinamika molekul dapat digunakan atau dikembangkan dalam penelitian selanjutnya.
iii
DAFTAR ISI DAFTAR ISI ......................................................................................................... i DAFTAR GAMBAR........................................................................................... iv DAFTAR NOTASI ............................................................................................. vi BAB 1 PENDAHULUAN.....................................................................................1 1.1
Latar Belakang Masalah ........................................................................1
1.2
Perumusan Masalah...............................................................................2
1.3
Tujuan Penelitian...................................................................................2
1.4
Pembatasan Masalah .............................................................................2
1.5
Sistematika Pembahasan .......................................................................3
BAB 2 DASAR TEORI........................................................................................5 2.1
Dinamika Molekul.................................................................................5
2.1.1
Ensemble ...........................................................................................6
2.1.1.1
Ensemble Mikrokanonikal (E, V, N)..........................................6
2.1.1.2 Ensemble Kanonikal (T, V, N)...................................................7 2.1.1.3
Ensemble Isobarik-Isotermal (P, T, N).......................................7
2.1.2
Langkah-Langkah Simulasi Dinamika Molekul................................7
2.1.3
Syarat Batas Periodik ........................................................................8
2.2
Mekanika Klasik..................................................................................10
2.3
Mekanika Statistik ...............................................................................11
2.3.1
Energi Total .....................................................................................11
2.3.2
Temperatur ......................................................................................12
2.3.3
Distribusi Kecepatan........................................................................12
2.3.4
Tekanan ...........................................................................................13
2.3.5
Entalpi..............................................................................................14
2.4
Model Interaksi Antar Molekul ...........................................................14
2.4.1
Potensial Lennard-Jones ..................................................................15
2.4.2
Gaya.................................................................................................16
2.4.3
Pemotongan Potensial .....................................................................17
2.5
Pengendalian Temperatur ....................................................................17
i
2.5.1
Velocity Scaling ..............................................................................18
2.5.2
Termostat Nose-Hoover ..................................................................18
2.5.3
Termostat Berendsen .......................................................................20
2.6
Pengendalian Tekanan.........................................................................21
2.7
Pengembangan Persamaan Gerak Diskrit............................................21
2.7.1
Algoritma Verlet..............................................................................22
2.7.2
Algoritma Leap-Frog.......................................................................23
2.7.3
Algoritma Verlet dengan Pengendalian Temperatur .......................24
2.7.4 Pendiskritan ζ untuk Termostat Nose-Hoover dan Berendsen ......25 BAB 3 PERANCANGAN SIMULASI DINAMIKA MOLEKUL.....................26 3.1
Garis Besar ..........................................................................................26
3.1.1
Prosedur Simulasi ............................................................................26
3.1.2
Parameter Simulasi ..........................................................................29
3.1.2.1
Kotak Simulasi dan Syarat Batas Periodik ...............................29
3.1.2.2
Pemotongan Potensial ..............................................................30
3.1.2.3
Posisi Awal ...............................................................................30
3.1.3
Paralelisasi .......................................................................................31
3.1.3.1
Dinamika molekul paralel dengan dekomposisi atom..............32
3.1.3.2
Dinamika molekul paralel dengan dekomposisi geometrik .....32
3.2
Perangkat Lunak ..................................................................................33
3.2.1
Kandidat Kelas ................................................................................33
3.2.2
Atribut Kelas ...................................................................................33
3.2.2.1
Simulation.................................................................................34
3.2.2.2
Molekul.....................................................................................34
3.2.2.3
Parameter ..................................................................................35
3.2.2.4
Property ....................................................................................36
3.2.2.5
Initialize ....................................................................................37
3.2.2.6
Verlet ........................................................................................37
3.2.2.7
LennardJones ............................................................................37
3.2.3
Operasi Kelas...................................................................................37
3.2.3.1
Simulation.................................................................................37
ii
3.2.3.2
Molecule ...................................................................................38
3.2.3.3
Parameter ..................................................................................38
3.2.3.4
Property ....................................................................................38
3.2.3.5
Initialize ....................................................................................38
3.2.3.6
Verlet ........................................................................................40
3.2.3.7
LennardJones ............................................................................40
3.2.4
Hubungan Antar Kelas ....................................................................41
3.2.5
Interface (Antar Muka)....................................................................42
BAB 4 ANALISA ...............................................................................................43 4.1
Simulasi Sistem Mikrokanonikal ........................................................43
4.1.1
Uji Kesalahan Energi Simulasi........................................................43
4.1.2
Momentum ......................................................................................46
4.1.3
Uji Simulasi Terhadap Temperatur Awal........................................48
4.2
Simulasi Ensemble Kanonikal.............................................................50
4.2.1
Termostat Nose Hoover...................................................................50
4.2.1.1
Pengujian Terhadap Parameter Nose Hoover Q.......................50
4.2.1.2
Uji Hamiltonian ........................................................................52
4.2.1.3
Error Temperatur pada Perubahan Tset .....................................53
4.2.1.4
Penurunan Matematis ...............................................................54
4.2.2
Termostat Berendsen .......................................................................57
4.2.2.1
Pengujian Terhadap Parameter Berendsen Q ...........................57
4.2.2.2
Error Temperatur pada Perubahan Tset .....................................58
4.2.2.3
Penurunan Matematis ...............................................................59
4.3
Tekanan terhadap Temperatur .............................................................61
4.4
Entalpi..................................................................................................63
4.5
Posisi....................................................................................................64
BAB 5 KESIMPULAN DAN SARAN ...............................................................66 5.1
Kesimpulan..........................................................................................66
5.2
Saran ....................................................................................................66
DAFTAR PUSTAKA..........................................................................................68 LAMPIRAN A LISTING PROGRAM ...............................................................70
iii
DAFTAR GAMBAR Gambar 2.1 Sistem yang diamati dalam simulasi dinamika molekul................................5 Gambar 2.2 Langkah-langkah simulasi dinamika molekul. [2] ........................................8 Gambar 2.3 Pembagian materi menjadi sel primer dan sel citra. ......................................9 Gambar 2.4 Distribusi kecepatan molekul pada berbagai temperatur [2] .......................13 Gambar 2.5 Potensial Lennard-Jones ..............................................................................16 Gambar 2.6 Pemasangan reservoir kalor pada sistem .....................................................19 Gambar 2.7 Pemasangan reservoir kalor dan reservoir tekanan pada sistem..................21 Gambar 3.1 Penghasilan trajektori molekul. [1]..............................................................26 Gambar 3.2 Prosedur simulasi dinamika molekul...........................................................28 Gambar 3.3 Syarat batas periodik untuk molekul yang keluar dari kotak simulasi. .......29 Gambar 3.4 Syarat batas periodik untuk jarak antar molekul .........................................30 Gambar 3.5 Penggunaan linked-list untuk molekul ........................................................35 Gambar 3.6 Penambahan molekul dalam linked-list [11] ...............................................39 Gambar 3.7 Hubungan antar kelas ..................................................................................42 Gambar 3.8 Contoh tampilan Simulasi Dinamika Molekul ............................................42 Gambar 4.1 Grafik perbandingan energi kinetik dan energi potensial............................44 Gambar 4.2 Grafik error energi total pada ∆t = 0,005 ....................................................44 Gambar 4.3 Grafik error energi total pada ∆t = 0,01 ......................................................45 Gambar 4.4 Grafik error energi total pada ∆t = 0,02 ......................................................45 Gambar 4.5 Grafik error energi total pada ∆t = 0,05 ......................................................46 Gambar 4.6 Grafik momentum total pada sumbu x, y dan z...........................................47 Gambar 4.7 Grafik temperatur terhadap waktu ...............................................................48 Gambar 4.8 Error energi total besar pada temperatur tinggi ...........................................49 Gambar 4.9 Lompatan energi potensial...........................................................................49 Gambar 4.10 Grafik temperatur dengan termostat Nose-Hoover Q = 1 .........................50 Gambar 4.11 Grafik temperatur dengan termostat Nose-Hoover Q = 5 .........................51 Gambar 4.12 Grafik temperatur dengan termostat Nose-Hoover Q = 15 .......................51 Gambar 4.13 Grafik Error Hamiltonian terhadap waktu.................................................52
iv
Gambar 4.14 Temperatur pada perubahan Tset dengan termostat Nose-Hoover..............53 Gambar 4.15 Perubahan ζ dan ζ& terhadap waktu.........................................................54 Gambar 4.16 Grafik temperatur dengan termostat Berendsen Q = 0,1 ...........................57 Gambar 4.17 Grafik temperatur dengan termostat Berendsen Q = 1 ..............................58 Gambar 4.18 Grafik temperatur dengan termostat Berendsen Q = 10 ............................58 Gambar 4.19 Temperatur pada perubahan Tset dengan termostat Berendsen ..................59 Gambar 4.20 Perbandingan tekanan dan temperatur sistem............................................61 Gambar 4.21 Perbandingan komponen Pm dan Pf dari tekanan sistem ...........................62 Gambar 4.22 Perbandingan entalpi dan tekanan sistem ..................................................63
v
DAFTAR NOTASI A
Luas
ai
Percepatan molekul i
E
Energi total
F
Gaya
H
Hamiltonian
H
Entalpi
K
Energi kinetik
kB
Konstanta Boltzmann
P
Tekanan
p
Momentum
L
Rusuk kotak simulasi
m
Massa
N
Jumlah molekul
Q
Parameter pengendalian temperatur
Rc
Radius cutoff
Ri
Posisi molekul i
Rij
Jarak antar molekul
T
Temperatur
Tset
Temperatur set
U
Energi potensial
V
Volume
vi
Kecepatan molekul i
∆t
Time step
ε
Parameter kekuatan interaksi
σ
Parameter panjang interaksi
ζ
Variabel pengendali temperatur
vi
1 2
β
β ≡ .ζ .∆t
s
Variabel Nose
vii
BAB 2 PENDAHULUAN
2.1
LATAR BELAKANG MASALAH
Teori-teori baru mengenai material pada skala atomik mempermudah peneliti untuk memprediksi perilaku material pada skala makroskopik dan memberikan kemampuan untuk merancang material-material baru dengan sifatsifat tertentu yang diinginkan. Namun analisa dan rancangan material dahulu hanya dapat dilakukan dengan eksperimen berkali-kali yang memerlukan biaya yang sangat mahal dan waktu yang lama. Selain itu, ada berbagai kondisi yang sulit atau tidak bisa diimplementasikan, antara lain eksperimen pada suhu yang sangat tinggi atau tekanan yang sangat besar. Dengan adanya kemajuan dalam ilmu komputer dan kemampuan komputasi yang jauh lebih kuat daripada dahulu, terbuka kemungkinan baru yaitu eksperimen komputer. Eksperimen komputer adalah jembatan antara teori dan eksperimen yang telah diterima sebagai salah satu metoda penelitian dan pengembangan material. Suatu eksperimen material secara fisik dapat didahului oleh eksperimen komputer untuk menentukan kondisi yang dibutuhkan atau memprediksi
hasilnya.
Eksperimen
komputer
dapat
berfungsi
untuk
mengkomplemen eksperimen [1]. Keuntungan eksperimen komputer adalah harga yang relatif lebih murah dan mudah meskipun untuk sistem yang kompleks. Selain itu, kemampuan eksperimen komputer meningkat sejalan dengan kemajuan komputer. Dinamika Molekul adalah teknik simulasi komputer yang mengamati pergerakan molekul-molekul yang saling berinteraksi. Teknik ini mensimulasikan molekul-molekul yang saling menarik dan mendorong dan menabrak satu sama lain. Simulasi dinamika molekul memberikan informasi statik dan dinamik pada skala atomik, seperti posisi dan kecepatan. Informasi ini lalu dapat diolah menjadi informasi pada skala makroskopis seperti tekanan, temperatur dan lain-lain.
1
Berbeda dengan beberapa teknik simulasi lainnya (lihat di [2]), dinamika molekul bersifat deterministik. Jika keadaan suatu materi diketahui pada waktu tertentu, maka keadaan materi tersebut pada waktu berbeda dapat ditentukan dengan sempurna.
2.2
PERUMUSAN MASALAH
Terdapat dua masalah utama saat merancang simulasi dinamika molekul. Pertama adalah pemodelan sistem yang terdiri dari model interaksi antar molekul dan interaksi antara sistem dengan lingkungan. Kedua adalah pemilihan algoritma dan metoda numerik. Pemodelan sistem akan menentukan kebenaran hasil simulasi dari segi fisis. Namun hasil yang lebih akurat memerlukan model yang lebih kompleks dan sulit dipecahkan dan membutuhkan waktu komputasi lebih lama. Pemodelan dan metoda numerik yang semakin akurat memerlukan komputasi yang lebih mahal. Karena salah satu hambatan terbesar adalah keterbatasan kemampuan komputasi, maka suatu kompromi dibutuhkan antara ketelitian simulasi dengan kecepatan simulasi.
2.3
TUJUAN PENELITIAN
Penelitian ini memiliki tujuan: 1. Memahami prinsip dinamika molekul 2. Mengembangkan simulasi dinamika molekul dengan algoritma yang tepat dan efisien 3. Melakukan uji coba simulasi dinamika molekul dengan berbagai parameter dan kondisi yang berbeda
2.4
PEMBATASAN MASALAH
Penelitian dengan simulasi dinamika molekul di laporan ini dibatasi: 1. Materi yang diamati berada dalam fasa padat sepanjang simulasi dengan formasi kristal kubus.
2
2. Simulasi dilakukan dibatasi pada energi total sistem konstan dan temperatur sistem konstan. 3. Model potensial intermolekuler yang digunakan adalah potensial LennardJones [2]. 4. Algoritma persamaan gerak yang digunakan adalah algoritma Verlet [2]. 5. Model interaksi sistem dengan lingkungan dibatasi pada metoda Velocity Scaling [1], thermostat Nose-Hoover [8] dan termostat Berendsen [9].
2.5
SISTEMATIKA PEMBAHASAN
Penelitian dilakukan dengan langkah-langkah sebagai berikut: 1. Studi literatur mengenai termodinamika statistik dan simulasi dinamika molekul. 2. Studi literatur mengenai teknik dan bahasa pemrograman C++. 3. Perancangan perangkat lunak simulasi dinamika molekul. 4. Menjalankan simulasi dengan berbagai parameter. 5. Analisa hasil simulasi. 6. Penulisan laporan. Laporan tugas akhir ini disusun dalam lima bab. Bab I
Pendahuluan Bab ini menjelaskan latar belakang, perumusan masalah, tujuan penelitian, pembatasan masalah penelitian, dan sistematika pembahasan.
Bab II
Dasar teori Bab ini menjelaskan landasan teori yang digunakan dalam penelitian ini, yaitu dasar dinamika molekul, model interaksi antar molekul, model interaksi sistem dengan lingkungan dan dasar mekanika statistik yang digunakan untuk mengolah informasi.
Bab III
Perancangan Bab ini membahas implementasi teori yang dijelaskan pada bab II ke dalam perangkat lunak. Perancangan perangkat lunak menggunakan bahasa pemrograman C++ dengan pemrograman berorientasi obyek.
Bab IV
Analisa
3
Bab ini memberikan hasil uji coba simulasi dinamika molekul untuk melihat kesesuaian dengan spesifikasi yang diinginkan dan analisa hasil yang diperoleh. Bab V
Kesimpulan dan Saran Bab ini memberikan kesimpulan penelitian dan saran-saran untuk penelitian selanjutnya.
4
BAB 3 DASAR TEORI
3.1
DINAMIKA MOLEKUL
Dinamika molekul mengamati molekul-molekul dalam suatu sistem tertutup, di mana jumlah materi (molekul) dalam sistem tidak berubah. Energi dapat keluar atau masuk sistem, tergantung dari jenis simulasi yang dilakukan. Suatu sistem adalah suatu kuantitas materi atau volume yang dipilih untuk diamati, sedangkan materi dan volume di luar sistem disebut lingkungan [3]. Pemisah antara sistem dan lingkungan disebut batas, yang secara teoritis tidak memiliki massa ataupun volume tersendiri. Sistem terbagi atas sistem tertutup, jika materi tidak dapat menembus batas, dan sistem terbuka, jika materi dapat menembus batas.
LINGKUNGAN
SISTEM
BATAS MOLEKUL
Gambar 3.1 Sistem yang diamati dalam simulasi dinamika molekul
Materi pada skala makroskopis terdiri dari molekul-molekul berjumlah sangat besar (bilangan Avogadro berorde 1023). Karena keterbatasan komputasi, maka simulasi dinamika molekul hanya dapat melakukan perhitungan untuk ratusan atau ribuan molekul. Semakin banyak molekul dalam simulasi, semakin realistis hasil yang diperoleh, tetapi biaya komputasi semakin mahal.
5
Tujuan pertama simulasi dinamika molekul adalah menghasilkan trajektori molekul-molekul sepanjang suatu jangka waktu terhingga. Pada setiap waktu, molekul-molekul dalam simulasi memiliki suatu posisi dan momentum tertentu untuk masing-masing sumbu. Untuk N molekul dalam ruang 3 dimensi, terdapat ruang posisi berdimensi 3N dan ruang momentum berdimensi 3N, sehingga terbentuk ruang fasa berdimensi 6N. Suatu konfigurasi posisi dan momentum molekul-molekul dapat diartikan sebagai koordinat dalam ruang fasa tersebut. Menurut Boltzmann [2], dengan jangka waktu yang cukup, suatu sistem akan pernah memiliki semua konfigurasi yang mungkin terjadi. Dengan kata lain, sistem tersebut akan pernah berada pada setiap koordinat dalam ruang fasa. Ruang fasa ini merupakan informasi keadaan mikroskopis sistem. Informasi ini lalu digunakan untuk memperoleh keadaan makroskopis sistem dan memprediksi ruang fasa pada waktu selanjutnya. Ruang fasa yang berbeda dapat menghasilkan keadaan makroskopis yang sama. Dengan kata lain, suatu keadaan makroskopis belum pasti berasal dari satu keadaan mikroskopis yang unik.
3.1.1
Ensemble
Suatu ensemble adalah koleksi dari keadaan sistem yang mungkin yang memiliki keadaan mikroskopis berbeda tetapi memiliki keadaan makroskopis sama [4]. Contohnya adalah sistem dengan konfigurasi posisi atau momentum yang berbeda namun memiliki temperatur yang sama. Beberapa ensemble yang sering digunakan dalam dinamika molekul adalah ensemble mikrokanonikal, ensemble kanonikal dan ensemble isobarikisotermal.
3.1.1.1 Ensemble Mikrokanonikal (E, V, N)
Ensemble mikrokanonikal adalah ensemble yang memiliki karakteristik jumlah molekul N dan volume yang tidak berubah serta energi total yang konstan. Ensemble ini diperoleh dari sistem yang terisolasi sehingga tidak ada interaksi sistem dengan lingkungan. Dengan demikian energi tidak dapat keluar dan masuk ke sistem dan energi total memiliki harga konstan. Ensemble ini biasa dinamakan
6
ensemble (E, V, N). Ensemble mikrokanonikal adalah ensemble yang paling sederhana untuk simulasi dinamika molekul, namun kurang praktis mensimulasi keadaan eksperimen dalam laboratorium. Ini disebabkan energi total sistem sulit dipertahankan konstan dalam eksperimen.
3.1.1.2 Ensemble Kanonikal (T, V, N)
Ensemble kanonikal adalah ensemble dengan keadaan makroskopis suhu yang tetap. Selain itu jumlah molekul N dan volume tidak berubah, maka dinamakan ensemble (T, V, N). Dalam laboratorium, temperatur sistem lebih mudah dikendalikan daripada energi total sistem, maka eksperimen sering dilakukan pada temperatur konstan. Ensemble kanonikal mendekati keadaan eksperimen pada temperatur konstan.
3.1.1.3 Ensemble Isobarik-Isotermal (P, T, N)
Dinamika molekul juga dapat dilakukan dengan mempertahankan tekanan dan temperatur sistem pada harga yang konstan. Tekanan dan temperatur adalah sifat makroskopis yang mudah dikendalikan dalam eksperimen. Dalam ensemble isobarik-isotermal, volume sistem dapat berubah atau menjadi suatu variabel. Jumlah molekul tidak berubah, maka ensemble ini juga dinamakan ensemble (P, T, N).
3.1.2
Langkah-Langkah Simulasi Dinamika Molekul
7
Pemodelan Dinamika Molekul
Pengembangan Model
Pemodelan Interaksi Antar Molekul
Simulasi Dinamika Molekul
Pemodelan Interaksi SistemLingkungan
Pengembangan Persamaan Gerak
Generasi Trajektori
Analisa Trajektori
Inisialisasi
Ekuilibrasi
Produksi
Gambar 3.2 Langkah-langkah simulasi dinamika molekul. [2] Dinamika molekul dilakukan dengan langkah-langkah berikut: 1. Pengembangan Model Pengembangan model harus dilakukan sebagai persiapan simulasi dinamika molekul. Model dapat diperoleh dari teori atau eksperimen. Model interaksi antar molekul dibutuhkan jika molekul berinteraksi satu sama lain. Model interaksi antara sistem dan lingkungan dibutuhkan jika sistem
tidak
terisolasi.
Interaksi
antara
sistem
dan
lingkungan
dipergunakan untuk pengendalian temperatur sistem dan pengendalian tekananc sistem. 2. Simulasi Dinamika Molekul Tahap simulasi dilakukan dengan penghasilan trajektori molekul-molekul dalam simulasi, lalu analisa yang dapat dilakukan bersamaan atau setelah simulasi. Generasi trajektori dibagi menjadi 3 tahap, yang akan dibahas lebih mendalam pada bab 4.1.
3.1.3
Syarat Batas Periodik
Molekul-molekul dekat batas sistem atau permukaan akan menimbulkan masalah karena memiliki molekul tetangga yang lebih sedikit daripada molekul yang berada di tengah sistem. Pada materi berskala makroskopis, molekul dekat
8
permukaan jauh lebih sedikit daripada molekul total sistem sehingga efek permukaan sangat kecil. Simulasi dinamika molekul sangat terpengaruh oleh efek permukaan. Ini menyebabkan informasi yang didapatkan adalah sifat materi dekat permukaan, padahal yang lebih penting diamati adalah sifat materi itu sendiri. Untuk menghindari ini, maka interaksi molekul dengan batas dihilangkan dengan menggunakan syarat batas periodik. [5] Material yang diamati akan dibagi menjadi sel-sel yang identik satu sama lain. Sel yang diamati disebut sebagai sel primer atau sel komputasi, sedangkan sel lain yang tidak diamati disebut sebagai sel citra. Di dalam sel citra terdapat citra-citra dari molekul-molekul dalam sel utama. Sel citra memiliki semua informasi (misal jumlah, posisi relatif dan kecepatan molekul) sama dengan sel primer. Sel citra adalah hasil replika dari sel primer secara periodik ke segala arah.
Sel Citra
Sel Citra
Sel Citra
Sel Citra
Sel Primer
Sel Citra
Sel Citra
Sel Citra
Sel Citra
Gambar 3.3 Pembagian materi menjadi sel primer dan sel citra. Ada dua implikasi dari syarat batas periodik. Pertama, jika sebuah molekul meninggalkan sel primer, molekul tersebut akan digantikan oleh citranya yang masuk ke dalam sel primer secara bersamaan. Posisi molekul yang keluar dari sel primer tersebut diganti dengan posisi baru yaitu posisi citranya yang masuk ke dalam sel primer. Kondisi ini akan menyebabkan jumlah atom N yang berada dalam sel primer akan konstan.
9
Implikasi kedua syarat batas periodik adalah konvensi citra minimum (minimum-image
convention).
Molekul-molekul
dalam
sel
utama
dapat
berinteraksi baik dengan molekul-molekul lain dalam sel utama maupun dengan citra-citranya dalam sel-sel citra. Interaksi antara dua molekul hanya diperhitungkan satu kali, dengan memilih interaksi molekul pertama dengan molekul kedua yang terdekat, baik dalam sel utama atau citranya dalam sel citra [6].
3.2
MEKANIKA KLASIK
Dalam dinamika molekul digunakan ketiga Hukum Newton: 1. Suatu partikel akan tetap diam atau bergerak dengan kecepatan tetap kecuali jika menerima gaya-gaya eksternal dengan resultan tidak sama dengan nol.
r 2. Jika partikel dengan massa m menerima gaya F , maka partikel itu akan mengalami percepatan sebesar
r r F a= m
(3.1)
r 3. Jika partikel i memberikan gaya pada partikel j sebesar Fij , maka partikel j r memberikan gaya pada partikel i sebesar − Fij . r r Fij = − F ji (3.2)
Hukum
Newton
ini
memberikan
konsekuensi
hukum
kekekalan
momentum. Dalam suatu sistem terisolasi, momentum masing-masing molekul dapat berubah-ubah akibat interaksi satu sama lain, namun momentum total tidak akan berubah. Momentum total sistem dapat diamati untuk memeriksa kebenaran simulasi ensemble mikrokanonikal.
d d ( ∑ p i ) = ( ∑ mi v i ) = 0 dt i dt i
(3.3)
di mana m adalah massa molekul dan p adalah momentum molekul.
10
3.3
MEKANIKA STATISTIK
Mekanika statistik atau termodinamika statistik dibutuhkan untuk mengkonversikan informasi pada skala atomik menjadi informasi pada skala makroskopik
[5].
Konfigurasi
posisi
dan
momentum
molekul-molekul
menentukan sifat-sifat yang dimiliki materi tersebut. Sifat-sifat itu antara lain adalah energi, temperatur, tekanan dan entalpi. Menurut mekanika statistik, kuantitas fisis diperoleh sebagai rata-rata konfigurasi tersebut terhadap waktu.
3.3.1
Energi Total
Energi total suatu sistem tersusun dari energi potensial sistem dan energi kinetik sistem. Energi potensial adalah jumlah dari semua energi potensial molekul-molekul dalam sistem.
U = ∑U i (R N )
(3.4)
i
RN adalah set posisi titik pusat massa atom atau molekul, RN = {R1, R2, R3, … , RN}. Energi kinetik sistem adalah jumlah dari energi kinetik setiap molekul.
K = ∑ K i (vi )
(3.5)
i
2
dengan K i =
1 p mi (vi ) 2 = i . 2 2mi
Untuk sistem terisolasi di mana tidak ada energi yang menembus batas, sistem bersifat konservatif atau energi sistem konstan. Konservasi energi ini adalah salah satu cara untuk memperiksa kebenaran simulasi ensemble mikrokanonikal.
11
3.3.2
Temperatur
Menurut termodinamika statistik, temperatur tidak lain adalah suatu skala dari energi kinetik molekul-molekul penyusunnya. Untuk tiga dimensi, hubungan antara energi kinetik dengan temperatur dinyatakan oleh K=
3 Nk BT 2
(3.6)
T=
2 K 3 Nk B
(3.7)
Atau,
di mana K adalah energi kinetik total sistem, N adalah jumlah molekul sistem, kB adalah konstanta Boltzmann dan T adalah temperatur.
3.3.3
Distribusi Kecepatan
Molekul-molekul dalam materi dapat memiliki kecepatan yang berbedabeda, sehingga terbentuk suatu distribusi kecepatan. Secara statistik dapat diperoleh bahwa molekul-molekul akan paling banyak berada pada suatu kecepatan tertentu, dan semakin berkurang jumlah molekulnya dengan semakin jauh kecepatan dari kecepatan tersebut. Salah satu penyebabnya adalah karena molekul-molekul dalam materi akan saling bertabrakan dan berinteraksi. Interaksi ini menyebabkan adanya pemerataan energi kinetik, karena molekul yang bergerak lebih cepat memberikan tambahan momentum pada molekul yang bergerak lebih lambat dan sebaliknya. Distribusi kecepatan yang terjadi berbentuk distribusi normal, dan dinamakan Distribusi Maxwell-Boltzmann [7]. Distribusi Maxwell-Boltzmann bergantung dari temperatur, dirumuskan: f (v)dv =
mv 2 1 exp − C 2kT
dv
(3.8)
12
Gambar 3.4 Distribusi kecepatan molekul pada berbagai temperatur [2]
3.3.4
Tekanan
Tekanan didefinisikan sebagai gaya yang bekerja tegak lurus pada suatu satuan luas. Px =
Fx A
(3.9)
Dengan menggunakan hukum Newton kedua, Px =
1 d (mvx ) A dt
(3.10)
13
Maka tekanan adalah suatu fluks momentum atau momentum yang menembus suatu satuan luas dalam suatu satuan waktu [2]. Menurut termodinamika statistik, ini terdiri dari dua bagian yaitu: 1. Pm, adalah fluks momentum akibat molekul yang menembus suatu permukaan luas selama dt. Pm =
2N K 3V
(3.11)
2. Pf, adalah fluks momentum akibat gaya yang bekerja antara dua molekul yang berada pada sisi yang berbeda dari permukaan luas. Pf =
1 3V
r r F ∑∑ ij ⋅ Rij i
(3.12)
j
Maka tekanan total menurut termodinamika statistik adalah P=
3.3.5
2N 1 K + 3V 3V
r
∑∑ F
ij
i
r ⋅ Rij
(3.13)
j
Entalpi
Entalpi didefinisikan sebagai gabungan dari energi total sistem dengan energi aliran [3]: H = E + PV
3.4
(3.14)
MODEL INTERAKSI ANTAR MOLEKUL
Model interaksi antar molekul yang diperlukan adalah hukum gaya antar molekul, yang ekivalen dengan fungsi energi potensial antar molekul. Pemilihan fungsi energi potensial harus dilakukan sebelum simulasi apa pun dapat dikerjakan. Pemilihan model interaksi antar molekul sangat menentukan kebenaran simulasi dari sudut pandang fisika. Karena berada dalam skala atomik, interaksi secara prinsip harus diturunkan secara kuantum, di mana berlaku prinsip ketidakpastian Heisenberg. Namun kita dapat melakukan pendekatan mekanika klasik di mana atom atau molekul dianggap sebagai suatu titik massa.
14
Model interaksi itu harus memenuhi dua buah kriteria. Pertama, molekulmolekul harus mampu menahan tekanan pasangan molekul yang saling berinteraksi. Ini dapat diartikan bahwa ada gaya tolak-menolak antar molekul. Kedua, molekul-molekul itu harus saling mengikat, atau ada gaya tarik-menarik antara molekul. Jika molekul-molekul terlalu dekat, gaya resultan adalah gaya tolak-menolak. Sebaliknya, jika terlalu jauh, maka gaya resultan adalah gaya tarik-menarik. Pada suatu jarak tertentu, kedua gaya tersebut saling meniadakan sehingga gaya resultannya sama dengan nol. Untuk N jumlah atom dalam suatu simulasi maka fungsi energi potensial adalah U(RN) di mana RN adalah set posisi titik pusat massa atom atau molekul, RN = {R1, R2, R3, … , RN}. Model energi potensial yang sederhana dan umum digunakan adalah pair-wise, di mana potensial adalah jumlah dari interaksi antara dua molekul yang diisolasikan.
U = ∑∑ U ( Rij ) ; i < j i
3.4.1
(3.15)
j
Potensial Lennard-Jones
Salah satu model energi potensial antara dua molekul yang dikembangkan adalah Potensial Lennard-Jones. Model ini dianggap paling sederhana, namun memiliki ketelitian yang baik untuk simulasi. Model Potensial ini dirumuskan: σ U ( Rij ) = kε Rij
n
σ − R ij
m
(3.16)
n dan m adalah bilangan bulat positif yang dipilih di mana n > m. n n k= n−mm
m /( n − m )
dan n > m
Dan i dan j adalah indeks dari molekul, Rij ≡ Ri − R j atau jarak antara molekul i dan j. Sedangkan σ adalah parameter jarak, dan ε adalah parameter yang menyatakan kekuatan interaksi. Pilihan yang umum untuk m dan n adalah m=6 dan n=12, sehingga persamaan (2.16) menjadi:
15
σ U ( Rij ) = 4ε Rij
3.4.2
12 6 σ − R ij
(3.17)
Gaya
Gaya adalah negatif dari gradien potensial. Untuk potensial LennardJones, besar gaya adalah:
ε d F ( Rij ) = − U ( Rij ) = 24 dr σ
σ 2 Rij
13 7 σ − R ij
(3.18)
Menurut konvensi, gaya positif adalah gaya tolak-menolak, dan gaya negatif adalah gaya tarik-menarik. Model ini menggambarkan adanya gaya tolak-
σ menolak dengan suku r
13
yang mendominasi pada jarak dekat dan gaya tarik-
7
σ menarik dengan suku yang mendominasi pada jarak jauh. r 4 3
Potensial
2
Gaya
1 0
ε
-1 -2
σ
-3 0.6
0.8
1
1.2
1.4
1.6
1.8
2
2.2
2.4
2.6
2.8
Gambar 3.5 Potensial Lennard-Jones
Dapat diturunkan untuk masing-masing sumbu:
F = −∇ ⋅ U ( Rij )
16
12 6 σ rx σ ∂ 2 − Fx = − U ( Rij ) = 24ε R ∂x Rij Rij ij 12 6 σ ry σ ∂ 2 − Fy = − U ( Rij ) = 24ε R ∂y Rij Rij ij
(3.19)
12 6 σ rz σ ∂ 2 − Fz = − U ( Rij ) = 24ε R ∂z Rij Rij ij
3.4.3
Pemotongan Potensial
Penghitungan gaya-gaya yang terjadi antar molekul adalah proses yang paling lama dalam simulasi dinamika molekul. Dalam prakteknya, sering kali potensial diberikan jarak cutoff Rc dan interaksi antar atom yang berjarak lebih besar dari Rc diabaikan. Cutoff ini dipasang pada jarak di mana interaksi seperti potensial dan gaya sudah kecil dan dapat diabaikan. 4ε σ U ( Rij ) = Rij 0
3.5
12 6 σ − R ij
Rij ≤ Rc
(3.20)
Rij > Rc
PENGENDALIAN TEMPERATUR
Metoda dinamika molekul dengan ensemble mikrokanonikal tidak selalu adalah cara terbaik untuk mendapatkan rata-rata statistik tertentu. Eksperimen laboratorium lebih sering dilakukan pada temperatur konstan (ensemble kanonikal T,V,N) daripada energi konstan (ensemble kanonikal E,V,N), karena temperatur lebih mudah dikendalikan pada skala makroskopis. Namun, untuk ensemble kanonikal diperlukan adanya fungsi tambahan untuk menjaga temperatur. Fungsi tambahan ini merupakan interaksi sistem dengan lingkungan di sekitarnya. Beberapa
model
interaksi
sistem-lingkungan
yang
dikembangkan
melibatkan adanya heat bath atau reservoir kalor yang berinteraksi dengan sistem. Reservoir adalah sesuatu di luar sistem yang memiliki temperatur tertentu yang
17
tidak akan berubah meskipun menerima dari atau memberikan kalor pada sistem. Kalor berpindah dari temperatur tinggi ke temperatur rendah, maka jika reservoir lebih panas daripada sistem, kalor berpindah dari reservoir menuju sistem, dan sebaliknya.
3.5.1
Velocity Scaling
Pengendalian temperatur dapat dilakukan dengan menyesuaikan harga kecepatan masing-masing molekul terus menerus atau pada periode tertentu. Tset ′ vi vi = T
(3.21)
Dengan metoda ini dapat dipastikan bahwa temperatur sistem sesaat akan tepat berada pada temperatur set point. Namun teknik ini tidak menjelaskan bagaimana interaksi sistem-lingkungan bekerja. Kelemahan lain teknik ini adalah sifat-sifat materi tidak dapat diturunkan dengan benar, karena pengambilan informasi harus dilakukan pada saat sistem bebas dan tanpa paksaan. Velocity scaling berguna dalam kecepatan awal atau jika simulasi suatu saat perlu menaikkan atau menurunkan temperatur. Cara lain adalah dengan secara periodik, namun tidak setiap saat, melakukan velocity scaling sehingga sistem lambat laun akan berkisar pada temperatur yang diinginkan. Teknik ini disebut quench method [1].
3.5.2
Termostat Nose-Hoover
Shuichi
Nose
[8]
melakukan
pengendalian
temperatur
dengan
memasangkan reservoir kalor pada sistem. Reservoir kalor dapat memberi atau menerima energi kalor dari sistem. Perpindahan energi dari lingkungan ke sistem bersifat non-lokal, atau terjadi di semua lokasi dalam sistem secara bersamaan. Molekul-molekul dalam sistem lalu akan dipercepat jika menerima kalor atau diperlambat jika memberi kalor sehingga distribusi kecepatannya mendekati bentuk distribusi kecepatan Maxwell-Boltzmann pada temperatur
yang
diinginkan.
18
ISOLASI
RESERVOIR KALOR
ENERGI
SISTEM
SISTEM
Gambar 3.6 Pemasangan reservoir kalor pada sistem
Dengan termostat Nose, molekul-molekul dalam sistem menerima gaya tambahan yang merupakan fungsi dari kecepatannya. Gaya ini dapat diumpamakan sebagai gaya gesek. Nose membuat suatu variabel baru s yang mengatur pertukaran kalor sistem dengan lingkungan. Hoover menyederhanakan ini dengan mendefinisikan variabel dinamis ζ :
ζ ≡
s& s
(3.22)
Pengendalian temperatur ini lalu dinamakan termostat Nose-Hoover dengan persamaan gerak baru dan penentuan harga ζ sebagai berikut: &r&i (t ) = −
1 ∂U ({R}) − ζ (t )r&i (t ) Mi ∂R
(
)
Qζ& (t ) = ∑ mi r&i (t ) 2 − 3Nk B Tset
(3.23) (3.24)
i
Dengan definisi temperatur dalam persamaan (2.7), persamaan dapat dituliskan Qζ& = 3(k B T − k B Ts )
(3.25)
di mana Q adalah parameter termostat yang menunjukkan kekuatan pengendalian temperatur dan Ts adalah temperatur set-point.
19
Termostat Nose-Hoover mengendalikan temperatur sistem dengan mengatur ζ& (perubahan ζ terhadap waktu). Jika temperatur sistem melebihi Tset,
ζ& akan menjadi positif dan ζ bertambah. Karena ζ berkisar di sekitar nol, maka setelah beberapa waktu ζ akan menjadi positif. Jika ζ positif maka umpan balik negatif bekerja dengan memperlambat molekul-molekul, sehingga temperatur sistem turun. Jika temperatur sistem lebih rendah daripada Tset, ζ& akan menjadi negatif dan ζ berkurang. Setelah beberapa waktu, ζ menjadi bernilai negatif dan molekul-molekul dipercepat, sehingga temperatur sistem naik. Termostat Nose-Hoover mengumpamakan perpindahan kalor memiliki ‘inersia termal’, dan seluruh sistem dengan reservoir kalor mengkonservasi nilai Hamiltonian H. Hamiltonian tidak memiliki nilai fisis namun dapat dipergunakan untuk memerika kebenaran simulasi. Hamiltonian ini dirumuskan H =
3.5.3
1 1 M i ( sR& i ) 2 + U ({R}) + Qs& 2 + gk B T ln s ∑ 2 2 i
(3.26)
Termostat Berendsen
Metoda pengendalian temperatur lain yang memiliki banyak kesamaan dengan termostat Nose-Hoover dirumuskan oleh Berendsen et al. [9]. Termostat Berendsen menggunakan persamaan gerak yang sama dengan termostat NoseHoover, &r&i (t ) = −
1 ∂U ({R}) − ζ (t )r&i (t ) Mi ∂R
Perbedaannya adalah termostat Berendsen bukan mengatur harga ζ& melainkan harga ζ dengan persamaan:
ζ =
1 (k B T − k B Ts ) 2Qk B T
(3.27)
di mana Q di sini adalah parameter termostat yang berupa konstanta waktu. Metoda ini mengarahkan temperatur sistem lebih langsung daripada metoda Nose-Hoover, sehingga kurva temperatur akan lebih tajam.
20
3.6
PENGENDALIAN TEKANAN
Tekanan merupakan fungsi dari temperatur dan volume sistem. Pengendalian tekanan dilakukan untuk simulasi ensemble isobarik-isotermal (P, V, N) dengan mengatur perubahan volume sistem sehingga tekanan tidak berubah. Sistem dijadikan piston yang dapat bergerak sehingga volume membesar atau mengecil. Sebuah reservoir tekanan lalu ditambahkan di luar sistem yang menggerakkan piston jika terdapat perbedaan tekanan.
ISOLASI
RESERVOIR KALOR
ENERGI
SISTEM
RESERVOIR TEKANAN
SISTEM
VOLUME
Gambar 3.7 Pemasangan reservoir kalor dan reservoir tekanan pada sistem
3.7
PENGEMBANGAN PERSAMAAN GERAK DISKRIT
Dinamika molekul mengamati posisi dan derivasinya, antara lain kecepatan dan percepatan. Maka persamaan gerak dapat dituliskan: dri (t ) = vi (t ) dt dvi (t ) = ai (t ) dt
(3.28)
Persamaan gerak ini harus dijadikan diskrit agar dapat dipecahkan secara numerik. Pendiskritan ini dilakukan dengan metoda beda hingga (finitedifference). Waktu yang diskrit memiliki langkah waktu (time step) ∆t yang merupakan selisih antara dua waktu berturutan. Metoda beda hingga dilakukan dengan melakukan ekspansi Taylor. Suatu metoda beda hingga yang mengaproksimasi solusi sebuah persamaan diferensial,
21
akan selalu menimbulkan truncation error (error pemotongan). Truncation error dihitung dari suku bukan nol pertama yang dihilangkan dari suatu ekspansi. Metoda yang dipilih harus memberikan error kecil namun tidak terlalu kompleks sehingga memerlukan waktu komputasi yang lama. Beberapa metoda integrasi persamaan gerak antara lain adalah algoritma Verlet, algoritma LeapFrog dan algoritma Velocity Form [9].
3.7.1
Algoritma Verlet
Algoritma Verlet sering digunakan karena algoritmanya yang sederhana namun memiliki ketelitian yang baik. Caranya adalah dengan menggunakan ekspansi Taylor untuk t + ∆t dan t − ∆t sebagai berikut: dri (t ) 1 d 2 ri (t ) 1 d 3 ri (t ) 2 (∆t ) + ∆t + (∆t ) 3 + O((∆t ) 4 ) ri (t + ∆t ) = ri (t ) + 3 2 2 dt 6 dt dt
(3.29)
dri (t ) 1 d 2 ri (t ) 1 d 3 ri (t ) 2 (∆t ) − ∆t + (∆t ) 3 + O((∆t ) 4 ) ri (t − ∆t ) = ri (t ) − 3 2 2 dt 6 dt dt
(3.30)
Penjumlahan persamaan (2.28) dan (2.29) menghasilkan: ri (t + ∆t ) + ri (t − ∆t ) = 2ri (t ) +
d 2 ri (t ) (∆t ) 2 + O((∆t ) 4 ) dt 2
ri (t + ∆t ) = 2ri (t ) + ri (t − ∆t ) +
d 2 ri (t ) (∆t ) 2 + O((∆t ) 4 ) 2 dt
atau ri (t + ∆t ) = 2ri (t ) + ri (t − ∆t ) + ai (t ).(∆t ) 2 + O((∆t ) 4 ) ri (t + ∆t ) = 2ri (t ) + ri (t − ∆t ) +
Fi (t ) .(∆t ) 2 + O((∆t ) 4 ) mi
(3.31) (3.32)
Maka diperoleh posisi molekul pada t+dt dengan truncation error berorde (∆t ) 4 . Sedangkan kecepatan pada t diperoleh dengan mengurangkan persamaan (2.29) dengan persamaan (2.30), menjadi ri (t + ∆t ) − ri (t − ∆t ) = 2
dri (t ) ∆t + O((∆t ) 3 ) dt
(3.33)
yang dapat dituliskan vi (t ) =
dri (t ) ri (t + ∆t ) − ri (t − ∆t ) + O((∆t ) 3 ) = dt 2∆t
(3.34)
22
Kelemahan dari algoritma Verlet adalah penanganan kecepatan yang kurang praktis, karena harus memprediksi posisi berikut sebelum dapat menghitung kecepatan sesaat. Selain itu, posisi sama sekali tidak ditentukan oleh kecepatan pada saat t, maka algoritma ini tidak mudah mempergunakan velocity scaling untuk simulasi pada T konstan.
3.7.2
Algoritma Leap-Frog
Algoritma Leap-Frog adalah metoda di mana kecepatan melompati percepatan, dan posisi melompati kecepatan. Menggunakan ekspansi Taylor sampai orde kedua untuk t + 12 ∆t dan t − 12 ∆t : dri (t ) 1 1 d 2 ri (t ) 1 ( 2 ∆t ) 2 + O((∆t ) 3 ) ri (t + ∆t ) = ri (t ) + 2 ∆t + 2 2 dt dt
(3.35)
dri (t ) 1 1 d 2 ri (t ) 1 ∆ t + ( 2 ∆t ) 2 + O((∆t ) 3 ) 2 2 2 dt dt
(3.36)
1 2
ri (t − 12 ∆t ) = ri (t ) −
Dengan cara yang sama dengan penghitungan kecepatan Verlet yaitu pengurangan (2.35) dengan (2.36) diperoleh ri (t + 12 ∆t ) − ri (t − 12 ∆t ) = 2 ri (t + 12 ∆t ) = ri (t − 12 ∆t ) +
dri (t ) 1 ( 2 ∆t ) + O((∆t ) 3 ) dt
dri (t ) ∆t + O((∆t ) 3 ) dt
(3.37)
Dengan mentranslasi seluruh persamaan sebesar 12 ∆t , diperoleh: ri (t + ∆t ) = ri (t ) +
dri (t + 12 ∆t ) ∆t + O((∆t ) 3 ) dt
atau ri (t + ∆t ) = ri (t ) + vi (t + 12 ∆t )∆t + O((∆t ) 3 )
(3.38)
Dengan cara yang sama diperoleh persamaan untuk kecepatan: vi (t + 12 ∆t ) − vi (t − 12 ∆t ) = 2 vi (t + 12 ∆t ) = vi (t − 12 ∆t ) +
dvi (t ) 1 ( ∆t ) + O((∆t ) 3 ) dt 2
dvi (t ) ∆t + O((∆t ) 3 ) dt
(3.39)
23
atau vi (t + 12 ∆t ) = vi (t − 12 ∆t ) +
3.7.3
Fi (t ) .∆t + O((∆t ) 3 ) mi
(3.40)
Algoritma Verlet dengan Pengendalian Temperatur
Algoritma Verlet dapat digunakan dengan pengendalian temperatur dengan modifikasi dengan menggabungkan persamaan (2.23), persamaan (2.32) dan persamaan (2.34):
&r&i (t ) = −
1 ∂U ({R}) − ζ (t )r&i (t ) Mi ∂R
ri (t + ∆t ) = 2ri (t ) + ri (t − ∆t ) + ai (t ).(∆t ) 2 + O((∆t ) 4 ) vi (t ) =
dri (t ) ri (t + ∆t ) − ri (t − ∆t ) + O((∆t ) 3 ) = dt 2∆t
maka persamaan gerak menjadi 1 ∂U ({R}) − ζ (t ) r&i (t ) .(∆t ) 2 + O((∆t ) 4 ) ri (t + ∆t ) = 2ri (t ) + ri (t − ∆t ) − ∂R Mi
(3.41)
1 ∂U ({R}) ri (t + ∆t ) − ri (t − ∆t ) + O((∆t ) 3 ) . − ζ (t ) ri (t + ∆t ) = 2ri (t ) + ri (t − ∆t ) − 2 ∂ M R ∆ t i 2 4 (∆t ) + O((∆t ) ) r (t + ∆t ) − ri (t − ∆t ) 1 ∂U ({R}) ri (t + ∆t ) = 2ri (t ) + ri (t − ∆t ) − ζ (t ) i .(∆t ) 2 − .(∆t ) 2 2∆t Mi ∂R
+ O((∆t ) 3 )
Didefinisikan 1 2
β (t ) ≡ ζ (t )∆t
(3.42)
sehingga persamaan (2.41) menjadi ri (t + ∆t ) = 2ri (t ) + ri (t − ∆t ) − β (t )( ri (t + ∆t ) − ri (t − ∆t )) −
1 ∂U ({R}) .( ∆t ) 2 + O ((∆t ) 3 ) Mi ∂R
(3.43)
24
ri (t + ∆t )(1 + β ) = 2ri (t ) + ri (t − ∆t )(1 − β ) − ri (t + ∆t ) =
3.7.4
1 ∂U ({R}) .(∆t ) 2 + O((∆t ) 3 ) Mi ∂R
1 1 ∂U ({R}) .(∆t ) 2 + O ((∆t ) 3 ) (3.44) 2ri (t ) + ri (t − ∆t )(1 − β ) − Mi 1+ β ∂R
Pendiskritan ζ untuk Termostat Nose-Hoover dan Berendsen
Dengan melakukan ekspansi Taylor untuk ζ pada t + ∆t dan t − ∆t , maka diperoleh persamaan-persamaan: 1 2
(3.45)
1 2
(3.46)
ζ (t + ∆t ) = ζ (t ) + ζ& (t ).∆t + ζ&&(t ).( ∆t ) 2 + O (( ∆t ) 3 ) ζ (t − ∆t ) = ζ (t ) − ζ& (t ).∆t − ζ&&(t ).(∆t ) 2 + O ((∆t ) 3 ) Persamaan (2.41) dikurangi persamaan (2.42) menjadi
ζ (t + ∆t ) − ζ (t − ∆t ) = 2ζ& (t ).∆t + O((∆t ) 3 )
(3.47)
Dengan definisi ς& :
(
)
Qζ& (t ) = ∑ mi r&i (t ) 2 − 3Nk B Ts i
maka diperoleh
ζ (t + ∆t ) = ζ (t − ∆t ) +
2∆t mi r&i (t ) 2 − 3Nk B Ts + O((∆t ) 3 ) ∑ Q i
(
)
(3.48)
25
BAB 4 PERANCANGAN SIMULASI DINAMIKA MOLEKUL
4.1
GARIS BESAR
4.1.1
Prosedur Simulasi
Simulasi dinamika molekul terdiri atas tiga tahap: inisialisasi, ekuilibrium, dan produksi. Dalam penelitian ini bagian ekuilibrium tidak dilakukan secara spesifik.
Gambar 4.1 Penghasilan trajektori molekul. [1]
26
1. Tahap inisialisasi terdiri dari penentuan sistem unit, algoritma dan parameter simulasi dan inisialisasi molekul-molekul. Inisialisasi molekul melibatkan penentuan posisi awal dan kecepatan awal molekul-molekul. 2. Tahap ekuilibrium diperlukan agar keadaan awal simulasi tidak dominan mempengaruhi analisa dari simulasi. Dalam penelitian ini bagian ekuilibrium sudah termasuk ke dalam tahap produksi. 3. Tahap produksi adalah tahap utama dalam simulasi dinamika molekul, saat memperoleh informasi dalam simulasi. Algoritma finite-difference digunakan berulang kali sambil mengumpulkan trajektori ruang fasa. Tahap produksi harus dilakukan dengan jangka waktu yang cukup untuk mengurangi ketidakpastian statistik sampai memenuhi spesifikasi.
27
Inisialisasi Parameter dan Molekul
Maju Satu Iterasi
Hitung Interaksi Antar Molekul
Kalkulasi Posisi dan Kecepatan
Tidak
Ya
Peroleh Sifat Materi
Pengendalian Temperatur
Iterasi Selesai?
Ya
Ulangi Simulasi
Tidak
Gambar 4.2 Prosedur simulasi dinamika molekul.
28
4.1.2
Parameter Simulasi
4.1.2.1 Kotak Simulasi dan Syarat Batas Periodik
Simulasi dinamika dilakukan dalam kubus dengan dimensi L x L x L, karena bentuk ini paling memudahkan syarat batas periodik. Maka syarat batas periodik diimplementasikan sebagai berikut: 1. Jika sebuah molekul melebihi batas kotak simulasi, maka xi = xi + αL dengan memilih α = …,-1,0,1,… sehingga molekul masih berada dalam kotak simulasi.
Gambar 4.3 Syarat batas periodik untuk molekul yang keluar dari kotak simulasi.
2. Jarak suatu molekul ke molekul kedua ditentukan sebagai jarak terdekat di dalam sel utama atau dengan citranya. Maka jarak dua molekul pada suatu sumbu adalah Rijx = Rijx + αL dengan memilih α = …,-1,0,1,… untuk mendapatkan Rijx terkecil.
29
Sel Primer
Sel Citra
Sel Citra
Sel Citra
Gambar 4.4 Syarat batas periodik untuk jarak antar molekul
4.1.2.2 Pemotongan Potensial
Penghitungan interaksi antar molekul memakan waktu paling lama dalam simulasi dinamika molekul. Agar komputasi tidak mahal, maka interaksi molekul hanya dihitung dengan tetangga terdekat saja. Teknik ini banyak digunakan dalam dinamika molekul materi dalam fasa padat. Pemotongan potensial dilakukan dengan memasang jarak cutoff di antara tetangga pertama dan tetangga kedua molekul. Jarak cutoff dipasang di antara jarak molekul ke tetangga pertama dan tetangga keduanya. Jarak cutoff dipasang di antara tetangga pertama yang berjarak L dan tetangga kedua yang berjarak L 2 , yaitu sebesar 1,2L.
4.1.2.3 Posisi Awal
Penempatan molekul-molekul secara acak dapat menyebabkan adanya molekul-molekul yang saling menumpuk atau berjarak terlalu dekat sehingga terjadi gaya tolak-menolak yang terlalu besar. Untuk mencegah hal ini, maka
30
molekul-molekul disusun dalam bentuk kubus sebagai posisi awal. Struktur ini sesuai dengan struktur materi dalam fasa padat. Panjang rusuk kubus adalah jarak di mana molekul-molekul berada pada potensial terendah. Untuk model potensial Lennard-Jones, maka jarak tersebut kalau hanya memperhitungkan tetangga terdekatnya adalah
Fx = − 0 = 2( 2(
(
σ Rij
Rij
σ
r σ ∂ σ U ( Rij ) = 24ε x 2( )12 − ( ) 6 = 0 ∂x Rij Rij Rij
σ Rij
)12 − (
)12 = (
σ Rij
σ Rij
)6
)6
)6 = 2
Rij = σ 6 2
4.1.3
(4.1)
(4.2)
Paralelisasi
Simulasi dinamika molekul memerlukan kemampuan komputasi yang besar untuk jumlah molekul banyak dan algoritma yang kompleks. Salah satu cara meningkatkan kemampuan komputasi adalah dengan menggunakan paralelisasi komputer [10]. Paralelisasi simulasi dinamika molekul dengan P prosesor terdiri dari mengambil sistem N partikel yang berinteraksi dan mendekomposisi menjadi P sistem yang lebih kecil yang dapat diproses secara paralel. Dengan P proses bekerja simultan maka kecepatan komputasi secara ideal akan berlipat P kali dari kecepatan satu prosesor. Dalam praktek, komunikasi antar komputer dan beban yang berbeda pada masing-masing prosesor akan menambah waktu simulasi. Dua teknik utama adalah dekomposisi molekular dan dekomposisi geometrik.
31
4.1.3.1 Dinamika molekul paralel dengan dekomposisi atom
Teknik ini adalah cara paling mudah untuk membagi-bagi suatu sistem. Untuk P prosesor, masing-masing prosesor menangani N/P molekul. Alokasi molekul pada satu prosesor tidak berdasar pada kondisi tertentu, seperti lokasi yang dekat. Molekul itu akan ditangani oleh prosesor tersebut sepanjang simulasi. Untuk mengetahui suatu molekul berinteraksi dengan molekul lain yang mana, setiap prosesor memiliki daftar tetangga dari molekul yang ditanganinya. Daftar tetangga adalah daftar molekul-molekul lain yang berada dalam bola dengan radius Rc dan berpusat pada molekul. Daftar tetangga perlu disusun kembali setiap periode tertentu. Teknik ini mempunyai kelemahan yaitu komunikasi antar prosesor untuk mempertukarkan informasi molekulnya dapat memakan waktu banyak.
4.1.3.2 Dinamika molekul paralel dengan dekomposisi geometrik
Pada teknik ini, ruang dipartisi menjadi banyak kotak yang disusun sehingga seluruh ruang kotak simulasi terpenuhi. Kotak partisi lalu didistribusikan pada banyak prosesor sehingga kotak partisi yang berdekatan secara fisik akan berada pada prosesor yang berdekatan. Pertukaran data jarak jauh akan diminimisasikan. Kelemahan teknik ini adalah jika molekul-molekul berpindah kotak partisi sehingga beban yang ditangani prosesor menjadi tidak seimbang. Prosesor dengan beban lebih kecil akan selesai lebih dahulu dan harus menunggu prosesor dengan beban terbesar selesai. Alternatif lain adalah dengan menyesuaikan distribusi ruang agar distribusi beban prosesor merata, namun ini menambah kompleksitas simulasi. Dekomposisi geometrik baik untuk digunakan dalam simulasi molekul dengan perpindahan kecil sehingga molekul-molekul kurang lebih berada tetap dalam kotak partisinya.
32
4.2
PERANGKAT LUNAK
Simulasi dinamika molekul dirancang dengan bahasa pemrograman C++ berorientasi objek dengan Borland C++ Builder Version 4.0. Simulasi dibentuk dengan berbagai kelas yang masing-masing memiliki tugas tertentu dalam simulasi. Penggunaan kelas dalam program ini mempermudah modifikasi program andaikata perlu mempergunakan model atau algoritma lain.
4.2.1
Kandidat Kelas
Nama Kelas
Type
Tanggung Jawab
Simulation
Process
Mengatur
langkah-
langkah
simulasi
dinamika molekul Molekul
Actor
-
Parameter
Process
Menginisialisasi parameter simulasi
Property
Process
Memperoleh sifat materi dari simulasi
Initialize
Process
Mengatur jalan inisialisasi simulasi
dinamika
molekul Verlet
Process
Menjalankan produksi
tahap dengan
algoritma Verlet LennardJones
Process
Menghitung antar
molekul
interaksi dengan
potensial LennardJones
4.2.2
Atribut Kelas
33
4.2.2.1 Simulation
Atribut
Data Type
Keterangan
t
integer
Nomor time step simulasi
Atribut
Data Type
Keterangan
index
integer
Nomor indeks molekul
4.2.2.2 Molekul
dalam kotak simulasi m
double
Massa molekul
R
array[3] of double
Posisi
molekul
pada
sumbu x, y dan z pada t Rnext
array[3] of double
Posisi
molekul
pada
sumbu x, y dan z pada t + ∆t
Rprev
array[3] of double
Posisi
molekul
pada
sumbu x, y dan z pada t − ∆t
v
array[3] of double
Kecepatan molekul pada sumbu x, y dan z pada t
F
array[3] of double
Gaya total yang diterima molekul pada sumbu x, y dan z
next
Molecule*
Pointer
ke
berikut
dalam
molekul kotak
simulasi Atribut kelas molekul bersifat public dan akan dimanipulasi oleh kelaskelas lain. Algoritma Verlet membutuhkan posisi pada waktu t − ∆t , t, dan t + ∆t , kecepatan dan percepatan, namun percepatan diggantikan dengan gaya yang dialami molekul agar lebih fleksibel untuk simulasi ensemble mikrokanonikal dan ensemble kanonikal.
34
Molekul-molekul dalam simulasi disusun dalam bentuk linked list (daftar terkait) di mana suatu molekul memiliki pointer (penunjuk) untuk molekul berikutnya. Pointer ini diberi nama next.
Gambar 4.5 Penggunaan linked-list untuk molekul
4.2.2.3 Parameter
Atribut
Data Type
Keterangan
N
integer
Jumlah
molekul
dalam
simulasi Naxis
integer
Jumlah
molekul
sepanjang rusuk kubus SIG
double
Parameter panjang dalam interaksi
EPS
double
Parameter
kekuatan
interaksi L0
double
Panjang rusuk satu kubus molekul
L
double
Panjang kubus simulasi
RC
double
Jarak cutoff
DT
double
Besar satu time step
T0
double
Temperatur awal sistem
Kb
double
Konstanta Boltzmann
TMAX
integer
Nomor iterasi maksimum
35
Q
double
Parameter kekuatan untuk pengendalian temperatur
Tset
double
Temperatur sistem yang diinginkan
4.2.2.4 Property
Atribut
Data Type
Keterangan
U
double
Energi potensial sistem
K
double
Energi kinetik sistem
E
double
Energi total sistem
E0
double
Energi total awal sistem
error
double
momentum
array[3] of double
error =
E (t ) − E (0) E (0)
Momentum pada
total
sistem
masing-masing
sumbu zeta
double
ζ
zetadot
double
ζ&
beta
double
β ≡ ζ∆t
s
double
s
sdot
double
s&
T
double
Temperatur sistem
P
double
Tekanan
Pf
double
Komponen tekanan yang
1 2
bergantung pada gaya Pm
double
Komponen tekanan yang bergantung
pada
momentum H
double
Entalpi
Ham
double
Hamiltonian
36
Ham0
double
HamError
double
Hamiltonian awal sistem HamError =
H (t ) − H (0) H (0)
4.2.2.5 Initialize
Kelas Initialize tidak memiliki atribut.
4.2.2.6 Verlet
Atribut
Data Type
Keterangan
zetaprev
double
ζ (t − ∆t )
zetanext
double
ζ (t + ∆t )
sprev
double
s (t − ∆t )
snext
double
s (t + ∆t )
Atribut
Data Type
Keterangan
Raxis
array[3] of double
Jarak dua molekul pada
4.2.2.7 LennardJones
sumbu x, y dan z Rij
double
Jarak dua molekul
force
array[3] of double
Gaya yang terjadi antar molekul i dan j atau Fij
4.2.3
Operasi Kelas
4.2.3.1 Simulation
Operasi
Return Type
Keterangan
Simulation()
-
Konstruktor kelas
Run()
void
Menjalankan
simulasi
37
dinamika molekul Reset()
void
Menghapus hasil simulasi
NextStep()
void
Memindahkan
variabel-
variabel t + ∆t menjadi t
4.2.3.2 Molecule
Operasi
Return Type
Keterangan
Molecule()
-
Konstruktor kelas
Operasi
Return Type
Keterangan
Parameter()
-
Konstruktor kelas
Set()
void
Memasang
4.2.3.3 Parameter
harga
parameter sesuai dengan harga input
4.2.3.4 Property
Operasi
Return Type
Keterangan
Property()
-
Konstruktor kelas
Run()
void
Menghitung
sifat-sifat
materi dari simulasi Kinetic()
void
Menghitung
energi
kinetik sistem Energy()
void
Menghitung energi total sistem dan besar error
4.2.3.5 Initialize
Operasi
Return Type
Keterangan
Initialize()
-
Konstruktor kelas
38
Run()
void
Menjalankan inisialisasi simulasi
dinamika
molekul CreateMolecule()
void
Menambahkan molekulmolekul baru ke dalam simulasi
InitialPosition()
void
Memasang
molekul-
molekul pada posisi awal InitialVelocity()
void
Memberikan awal
pada
kecepatan molekul-
molekul
head
NULL
Molekul N
head
NULL
Molekul N-1
head
Molekul N
NULL
...
Molekul N
Molekul 1
head
NULL
Gambar 4.6 Penambahan molekul dalam linked-list [11]
39
4.2.3.6 Verlet
Operasi
Return Type
Keterangan
Verlet()
-
Konstruktor kelas
Run()
void
Mengatur jalan algoritma Verlet pada setiap iterasi
RunOne()
void
Mengatur jalan algoritma Verlet
untuk
iterasi
pertama PredictPosition()
void
Memprediksi
posisi
molekul-molekul
pada
iterasi berikut PredictPositionOne()
void
Memprediksi
posisi
molekul-molekul setelah iterasi pertama Velocity()
void
Menghitung
kecepatan
molekul-molekul PredictZeta()
void
Memprediksi ζ (t + ∆t )
PredictS()
void
Memprediksi s (t + ∆t )
VelocityScale()
void
Melakukan
velocity
scaling jika diperlukan
4.2.3.7 LennardJones
Operasi
Return Type
Keterangan
LennardJones()
-
Konstruktor kelas
Run()
void
Mengatur interaksi antar molekul
ClearValues()
void
Mengosongkan
harga
potenial dan gaya yang dialami molekul Distance(p1:
Molecule*, void
Menghitung
jarak
dua
40
p2: Molecule*)
molekul pada
Force(p1: Molecule*, p2: void
Menghitung gaya yang
Molecule*)
terjadi antar dua molekul
Potential(p1:
Molecule*, void
Menghitung
p2: Molecule*)
potensial
yang terjadi antar dua molekul void
PressureByForce(p1: Molecule*,
Menghitung tekanan
p2:
komponen
yang
berupa
fungsi gaya antar dua
Molecule*)
molekul
4.2.4
Hubungan Antar Kelas Verlet -zetaprev : double -zetanext : double -snext : double -sprev : double +Verlet() +Run() : void +RunOne() : void -NextStep() : void -PredictPosition() : void -PredictPositionOne() : void -Velocity() : void -PredictZeta() : void -PredictS() : void -VelocityScale() : void
Simulation +t : int *
+Simulation() +Run() : void +Reset() : void -NextStep() : void *
*
* LennardJones
*
-RAxis[3] : double -Rij : double -force : double +LennardJones() +Run() : void -ClearValues() : void -Distance(in p1 : Molecule*, in p2 : Molecule*) : void -Force(in p1 : Molecule*, in p2 : Molecule*) : void -Potential(in p1 : Molecule*, in p2 : Molecule*) : void -PressureByForce(in p1 : Molecule*, in p2 : Molecule*) : void
* *
*
*
*
* *
*
* *
*
*
Property
*
Parameter
*
* *
* Initialize
+Initialize() +Run() : void -CreateMolecules() : void -InitialPosition() : void -InitialVelocity() : void
+N : int +NAxis : int +SIG : double +EPS : double +L0 : double +L : double +RC : double +DT : double +T0 : double +Kb : double +TMAX : int +Q : double +Tset : double +Parameter() +Set() : void
Molecule +index : int +m : double +R[3] : double +Rnext[3] : double +Rprev[3] : double +v[3] : double +F[3] : double +next : Molecule* +Molecule() *
*
*
* *
*
* * *
*
+U : double +K : double +E : double +E0 : double +error : double +momentum[3] : double +zeta : double +zetadot : double +beta : double +s : double +sdot : double +T : double +P : double +Pf : double +Pm : double +H : double +Ham : double +Ham0 : double +HamError : double +Property() +Run() : void -Kinetic() : void -Energy() : void
41
Gambar 4.7 Hubungan antar kelas
4.2.5
Interface (Antar Muka)
Tampilan perangkat lunak simulasi dinamika molekul dalam penelitian ini:
Gambar 4.8 Contoh tampilan Simulasi Dinamika Molekul
42
BAB 5 ANALISA
5.1
SIMULASI SISTEM MIKROKANONIKAL
5.1.1
Uji Kesalahan Energi Simulasi
Energi diamati untuk melihat kebenaran simulasi, karena dalam sistem terisolasi atau ensemble mikrokanonikal, energi tidak berubah sepanjang simulasi. Namun karena simulasi dinamika molekul menggunakan persamaan diskrit atau metode beda hingga, maka tidak terhindarkan timbulnya error. Perubahan energi potensial akan diimbangi oleh perubahan energi kinetik sehingga energi total tidak berubah.
43
Gambar 5.1 Grafik perbandingan energi kinetik dan energi potensial
Kesalahan energi total dihitung dengan error =
E (t ) − E (0) E (0)
(5.1)
Hasil perhitungan error dalam beberapa simulasi 512 molekul dengan berbagai besar time step ∆t : 1. ∆t = 0,005, 100000 time step
Gambar 5.2 Grafik error energi total pada ∆t = 0,005
2. ∆t = 0,01, 50000 time step
44
Gambar 5.3 Grafik error energi total pada ∆t = 0,01 3. ∆t = 0,02, 25000 time step
Gambar 5.4 Grafik error energi total pada ∆t = 0,02 4. ∆t = 0,05, 10000 time step
45
Gambar 5.5 Grafik error energi total pada ∆t = 0,05 Time step ( ∆t )
Orde Error
0,005
10-9
0,01
10-8
0,02
10-8
0,05
10-7 Tabel 5.1 Orde error terhadap besar time step ∆t
Dapat terlihat bahwa simulasi memiliki error lebih besar dengan ∆t yang lebih besar. Namun simulasi dengan ∆t lebih kecil memerlukan lebih banyak time step untuk menjalankan simulasi dengan waktu yang sama. Dari analisa ini kita pilih
∆t =0,02 untuk simulasi lain dalam penelitian agar dapat
mensimulasikan waktu lebih lama namun error masih berorde 10-8.
5.1.2
Momentum
Dalam simulasi ensemble mikrokanonikal, perubahan kecepatan dan momentum molekul hanya terjadi akibat gaya antar molekul. Karena Fij = -Fji, maka perubahan momentum molekul i akan diimbangi molekul j, sehingga momentum total tidak berubah.
46
Gambar 5.6 Grafik momentum total pada sumbu x, y dan z Terlihat momentum total pada ketiga sumbu hanya mengalami perubahan kecil (berorde 10-10) dan dapat diabaikan. Simulasi berjalan dengan baik dari sudut pandang momentum.
47
5.1.3
Uji Simulasi Terhadap Temperatur Awal
Uji simulasi ini mengamati perkembangan temperatur sejalan waktu dan kemudian menghitung temperatur rata-rata. Hasil simulasi 512 molekul dengan temperatur awal 0,05 selama 2500 langkah:
Gambar 5.7 Grafik temperatur terhadap waktu T = 0,02563
Terlihat bahwa temperatur rata-rata berbeda dengan temperatur awal. Dalam simulasi ensemble mikrokanonikal, temperatur rata-rata tidak dapat diprediksi atau diatur sebelum simulasi berakhir. Perolehan data untuk suatu temperatur yang ditentukan sebelum simulasi menjadi kurang efisien, karena sistem bisa jadi jarang berada pada temperatur yang diinginkan. Ini adalah kerugian simulasi ensemble mikrokanonikal. Dari hasil simulasi di atas dapat dilihat bahwa temperatur rata-rata ( T = 0,02563) lebih rendah daripada temperatur awal (T(0) = 0,05) dan temperatur sepanjang simulasi tidak pernah melebihi temperatur awal. Ini disebabkan karena pada awal simulasi, molekul-molekul sudah berada pada posisi berenergi potensial terendah. Perubahan posisi apa pun, baik saling mendekati atau saling menjauhi, selama semua molekul masih berada dalam kotak simulasi, akan
48
menyebabkan potensial yang lebih tinggi daripada potensial awal. Karena sistem bersifat konservatif, semua perubahan energi potensial akan selalu diimbangi dengan negtatif perubahan energi kinetik. Maka potensial yang lebih tinggi menjadikan energi kinetik lebih rendah, dan temperatur lebih rendah. Simulasi dengan temperatur awal terlalu tinggi dapat menyebabkan terjadi error energi total yang besar, seperti pada simulasi berikut:
Gambar 5.8 Error energi total besar pada temperatur tinggi Lompatan error energi total pada t ≈ 25 dan t ≈ 30 terjadi akibat lompatan energi potensial yang tidak diimbangi energi kinetik:
Gambar 5.9 Lompatan energi potensial
49
Kasus ini terjadi karena pada temperatur terlalu tinggi, terdapat molekul yang bergerak terlalu cepat, sehingga dapat menembus dari luar ke dalam radius cut-off. Dua molekul yang sebelumnya tidak berinteraksi lalu memiliki potensial baru berharga negatif dan gaya baru yang berupa gaya tarik.
5.2
SIMULASI ENSEMBLE KANONIKAL
5.2.1
Termostat Nose Hoover
5.2.1.1 Pengujian Terhadap Parameter Nose Hoover Q
Simulasi dengan pengendalian temperatur Nose Hoover dilakukan dengan parameter Tset sebagai temperatur yang diinginkan dan T sebagai temperatur sistem dan selisihnya sebagai error temperatur. Simulasi dilakukan untuk 512 molekul dan Tset = 0,02 selama 25000 langkah waktu, dengan Q divariasikan. 1. Q = 1
Gambar 5.10 Grafik temperatur dengan termostat Nose-Hoover Q = 1 2. Q = 5
50
Gambar 5.11 Grafik temperatur dengan termostat Nose-Hoover Q = 5 3. Q = 15
Gambar 5.12 Grafik temperatur dengan termostat Nose-Hoover Q = 15 Q
T
1
0,02004
5
0,02001
15
0,01981 Tabel 5.2 Temperatur rata-rata terhadap parameter termostat Nose Hoover Q Termostat Nose-Hoover menjadikan sistem berosilasi di sekitar harga Tset
= 0,02. Temperatur rata-rata T hampir sama dengan Tset, maka pengendalian temperatur berhasil dilakukan dengan baik.
51
Koreksi error temperatur bergantung pada parameter termostat NoseHoover Q. Hasil simulasi menunjukkan bahwa semakin besar Q, semakin besar perioda temperatur berosilasi, dengan kata lain koreksi lebih lambat. Ini sesuai dengan rumus:
(
)
Qζ& (t ) = ∑ mi r&i (t ) 2 − 3Nk B Ts i
Jika Q besar, ζ& menjadi kecil, sehingga perubahan ζ terhadap waktu menjadi lebih kecil. Maka perubahan temperatur oleh termostat Nose Hoover juga lebih lambat.
5.2.1.2 Uji Hamiltonian
Termostat Nose-Hoover tidak mempertahankan energi total sistem pada harga yang konstan melainkan harga total Hamiltonian sistem dan reservoir kalor. Maka simulasi ensemble kanonikal dengan termostat Nose-Hoover dapat diuji kebenarannya dengan melihat Hamiltonian sepanjang simulasi. Error Hamiltonian dihitung dengan error =
H (t ) − H (0) H (0)
(5.2)
Simulasi dijalankan dengan 512 molekul, temperatur awal 0,02 dan termostat Nose Hoover dengan Tset 0,02.
Gambar 5.13 Grafik Error Hamiltonian terhadap waktu Error Hamiltonian sepanjang simulasi memenuhi spesifikasi (berorde 10-4).
52
5.2.1.3 Error Temperatur pada Perubahan Tset
Pada bagian ini simulasi mulai tanpa pengendalian temperatur, lalu pada t=100 diberikan pengendalian temperatur dengan termostat Nose-Hoover dengan Tset = 0,02.
Gambar 5.14 Temperatur pada perubahan Tset dengan termostat Nose-Hoover Temperatur sistem sebelum pengendalian temperatur berkisar di sekitar 0,01. Setelah sistem diberikan pengendalian temperatur, temperatur sistem berosilasi di sekitar Tset = 0,02 dengan amplitudo kurang lebih (0,02-0,01)=0,01. Meskipun sistem dijalankan beberapa waktu, amplitudo kurang lebih tidak berubah, atau tidak teredam. Ini dapat terjadi karena termostat Nose-Hoover dianggap memiliki ‘massa’ dan ‘momentum’, sehingga perpindahan kalor ke dalam atau ke luar sistem memiliki ‘inersia thermal’ [9]. Momentum termal menyimpan energi kinetik yang ditambah atau dikurangi dari sistem ke lingkungan, sehingga bertindak seperti kapasitor. Termostat Nose-Hoover mengontrol ζ& dan bukan ζ itu sendiri. Maka dapat terjadi kasus di mana error temperatur dan ζ& adalah nol namun ζ bukanlah nol, menyebabkan adanya koreksi lebih meskipun pada saat itu temperatur sistem tidak perlu dikoreksi.
53
0.2
ζ&
ζ
0.1
0
-0.1
-0.2
Gambar 5.15 Perubahan ζ dan ζ& terhadap waktu
Jika Q besar, ζ& kecil dan ζ lambat berubah. Maka perubahan temperatur oleh termostat Nose Hoover juga lebih lambat. Perubahan temperatur yang lambat menyebabkan perubahan ζ& juga lambat. Maka perubahan ζ , ζ& dan T tersebut menjadi lebih lambat, namun hubungan satu sama lain tidak berubah.
5.2.1.4 Penurunan Matematis
Sistem dianggap berosilasi di sekitar temperatur rata-rata T = Tset dengan ~ error temperatur T . Kecepatan rata-rata molekul dalam sistem dinamakan v dan error kecepatan v~ . Andai tidak ada interaksi antar molekul, atau potensial dan gaya antar molekul adalah nol, maka percepatan molekul hanya diatur oleh termostat.
a = −ζv
(5.3) ~
ζ& berbanding lurus dengan T , atau: ~
ζ& = k1T dengan k1 ∝
(5.4)
1 . Q
ζ dapat dihitung dari integrasi ζ& :
54
ζ = ∫ ζ&dt
(5.5) ~
ζ = ∫ (k1T )dt ~
ζ = k1 ∫ T dt sehingga percepatan menjadi ~ a = −(k1 ∫ T dt ).v
dv ~ = −(k1 ∫ T dt ).v dt
(5.6)
~ Hubungan antara T dan v~ diperoleh sebagai berikut:
T = k 2 .v 2
(5.7)
dT = 2k 2 v dv ~ dT ~ .v T= dv ~ T = 2k 2 vv~
(5.8)
Persamaan (4.6) disubstitusi ke dalam persamaan (4.4)
dv = − k1 ∫ (2k 2 vv~ )dt.v dt dv dv~ + = −k1 ∫ (2k 2 (v + v~ )v~ )dt.(v + v~ ) dt dt Karena temperatur rata-rata tidak berubah terhadap waktu, maka
dv = 0. dt
dv~ + k1 ∫ (2k 2 (v + v~ )v~ )dt.(v + v~ ) = 0 dt dv~ + 2k1 k 2 ∫ ((v + v~ )v~ )dt.(v + v~ ) = 0 dt Jika dianggap osilasi temperatur sangat kecil, v~ << v maka v + v~ ≈ v dan persamaan dapat dilinierkan: dv~ + 2k1 k 2 ∫ (v v~ )dt.(v ) = 0 dt dv~ + 2k1 k 2 v 2 ∫ v~dt = 0 dt
55
Setelah diturunkan terhadap t, diperoleh persamaan diferensial linier orde 2:
d 2 v~ + 2k1 k 2 v 2 v~ = 0 2 dt Persamaan ini memiliki solusi umum v~ = Ae − λt
(5.9)
(5.10)
dengan λ = a + ib sehingga persamaan menjadi
λ2 Ae − λt + 2k1 k 2 v 2 Ae − λt = 0 maka
λ 2 + 2 k1 k 2 v 2 = 0
(5.11)
karena k1 , k 2 dan v 2 selalu positif, λ2 adalah negatif dan λ adalah suatu bilangan imajiner.
λ = i 2k1 k 2 v 2 , a = 0 dan b = v 2k1 k 2
(5.12)
Persamaan ini adalah persamaan harmonis sederhana
ω = v 2 k1 k 2
(5.13)
v~ = Ae iωt = A. cos(ωt + φ )
(5.14)
v~ berosilasi tanpa redaman. Untuk temperatur, ~ T = 2k 2 vv~
~ T = 2 Ak 2 v cos(ωt )
(5.15)
Error temperatur akan berosilasi tanpa redaman dengan perioda yang bergantung dari harga v , k1 dan k 2 . Karena k1 ∝
1 , maka dapat disimpulkan, Q
untuk error temperatur kecil, semakin besar parameter Nose-Hoover Q, semakin lambat osilasi temperatur. Dalam penelitian ini, error temperatur tidak kecil sehingga tidak berlaku v~ << v dan persamaan menjadi tidak linier, namun masih menunjukkan sifat-sifat
yang hampir sama.
56
Selain itu dapat disimpulkan bahwa semakin tinggi temperatur rata-rata juga osilasi semakin lambat.
5.2.2
Termostat Berendsen
5.2.2.1 Pengujian Terhadap Parameter Berendsen Q
Simulasi dengan pengendalian temperatur Berendsen dilakukan dengan 512 molekul dan Tset = 0,02 selama 25000 langkah waktu, dengan parameter Q divariasikan. 1. Q = 0,1
Gambar 5.16 Grafik temperatur dengan termostat Berendsen Q = 0,1 2. Q = 1
57
Gambar 5.17 Grafik temperatur dengan termostat Berendsen Q = 1 3. Q = 10
Gambar 5.18 Grafik temperatur dengan termostat Berendsen Q = 10
Q
T
0,1
0,02000
1
0,01996
10
0,01967 Tabel 5.3 Temperatur rata-rata terhadap parameter Berendsen Q Temperatur sistem berosilasi di sekitar Tset, berarti pengendalian
temperatur berhasil dijalankan dengan baik. Dari beberapa hasil simulasi dengan
Q bervariasi, terlihat bahwa error temperatur lebih besar pada Q besar. Ini terjadi karena Q lebih besar berarti pengendalian temperatur lebih lemah, dan menyebabkan sistem dapat berada pada temperatur yang lebih jauh dari Tset.
5.2.2.2 Error Temperatur pada Perubahan Tset
Simulasi dimulai tanpa pengendalian temperatur, lalu pada t=100 diberikan pengendalian Berendsen dengan Tset = 0,02.
58
Gambar 5.19 Temperatur pada perubahan Tset dengan termostat Berendsen Temperatur sistem sebelum pengendalian temperatur berkisar di sekitar 0,01. Setelah sistem diberikan pengendalian temperatur, temperatur sistem berosilasi di sekitar 0,02 dengan amplitudo error temperatur yang semakin kecil terhadap waktu. Setelah beberapa saat, temperatur sistem berkisar pada Tset dengan error kecil. Maka termostat Berendsen lebih baik beradaptasi pada perubahan Tset daripada termostat Nose-Hoover.
5.2.2.3 Penurunan Matematis
Seperti pada termostat Nose-Hoover, sistem dianggap berosilasi di sekitar ~ temperatur rata-rata T = Tset dengan error temperatur T . Kecepatan rata-rata molekul dalam sistem dinamakan v dan error kecepatan v~ . Andai tidak ada interaksi antar molekul, atau potensial dan gaya antar molekul adalah nol, maka percepatan molekul hanya diatur oleh termostat.
a = −ζv namun dalam termostat Berendsen digunakan: ~ T ζ = k1 T di mana k1 ∝
(5.16)
1 . Q
Sehingga percepatan menjadi
59
~ T a = − k1 .v T ~ dv T = −k1 .v dt T ~ Hubungan antara T dan v~ sudah diperoleh:
(5.17) (5.18)
T = k 2 .v 2 ~ T = 2k 2 vv~
(5.19)
Persamaan (4.19) disubstitusi ke persamaan (4.18) 2k 2 vv~ dv v = − k1 dt k 2 .v 2
(5.20)
dv dv~ + = −2k1v~ dt dt Temperatur rata-rata tidak berubah terhadap waktu, maka
(5.21)
dv = 0 . Diperoleh dt
persamaan diferensial linier berorde satu: dv~ + 2k1v~ = 0 dt Persamaan ini memiliki solusi umum v~ = Ae − λt
(5.22)
(5.23)
dengan λ = a + ib sehingga persamaan menjadi
− λAe − λt + 2k1 Ae − λt = 0 maka
λ = 2k1
(5.24)
karena k1 adalah bilangan real positif, λ adalah suatu bilangan real positif.
v~ = Ae −2 k1t Error temperatur dihitung ~ T = 2k 2 vv~ ~ T = 2 Ak 2 ve −2 k1t
(5.25)
60
Error temperatur berubah terhadap waktu menuju nol secara eksponensial. Semakin besar k1 , semakin cepat error temperatur menuju nol. Karena k1 ∝
1 , Q
dapat disimpulkan bahwa semakin besar Q, koreksi error temperatur berlangsung lebih lambat.
5.3
TEKANAN
Simulasi ini dilakukan dengan 512 molekul, σ =10,0 ε =0,625 ∆t =0,02 sebanyak 2500 iterasi.
Gambar 5.20 Perbandingan tekanan dan temperatur sistem Tekanan yang dihitung dari simulasi berfluktuasi hampir sama dengan temperatur sistem. Menurut teorema virial,
61
P=
2N 1 K + 3V 3V
r
∑∑ F
ij
i
r ⋅ Rij
j
Komponen pertama (Pm) dari tekanan berkaitan dengan temperatur dan komponen kedua (Pf) berkaitan dengan gaya antar molekul. Perbandingan antara kedua komponen tersebut dalam simulasi:
Gambar 5.21 Perbandingan komponen Pm dan Pf dari tekanan sistem Komponen Pm yang berorde 10-4 jauh lebih besar daripada Pf yang berorde 10-10. Ini disebabkan karena molekul-molekul berada pada posisi potensial terendah sehingga gaya antar molekul sangat kecil. Maka komponen kedua memberikan kontribusi yang jauh lebih kecil daripada komponen pertama.
62
Kondisi seperti ini juga dialami oleh gas ideal, di mana interaksi antar molekul sangat kecil hingga diasumsikan nol, sehingga persamaan tekanan dari teorema virial P=
2N 1 K + 3V 3V
r
∑∑ F
ij
i
r ⋅ Rij
j
menjadi persamaan gas ideal
PV = NkT
5.4
(5.26)
ENTALPI
Entalpi sistem diperoleh dari simulasi ensemble mikrokanonikal (sistem terisolasi) 512 molekul:
Gambar 5.22 Perbandingan entalpi dan tekanan sistem
63
Perubahan entalpi sistem mengikuti perubahan tekanan sistem, karena seperti dirumuskan
H = E+PV Pada simulasi ensemble mikrokanonikal, volume sistem pada simulasi ini adalah konstan dan energi total sistem juga dipertahankan konstan, maka faktor yang mempengaruhi perubahan entalpi sistem hanya adalah tekanan sistem.
5.5
POSISI
Pengamatan posisi molekul-molekul dalam simulasi 64 molekul pada beberapa temperatur tertentu menunjukkan: 1. T = 0,1
Pada temperatur ini, molekul-molekul bergetar di sekitar posisi awal yang berbentuk kubik. Melihat formasi molekul yang teratur ini, materi dapat dianggap berada dalam fasa padat. Tidak ada molekul yang memasuki atau keluar dari radius cutoff molekul lain, sehingga kesalahan simulasi kecil. 2. T = 1
64
Temperatur yang lebih tinggi menyebabkan molekul-molekul mulai lepas dari posisi awal. Terdapat molekul yang memasuki atau keluar dari radius cutoff molekul lain. 3. T = 10
Pada temperatur ini, molekul-molekul sudah lepas sama sekali dari posisi awal yang teratur. Ini dapat diartikan bahwa molekul sudah tidak berada pada fasa padat lagi melainkan fasa cair/gas. Keadaan ini harus dihindarkan, karena beberapa algoritma seperti algoritma Verlet tidak dapat dipergunakan untuk materi cair atau gas.
65
BAB 6 KESIMPULAN DAN SARAN
6.1
KESIMPULAN
1. Simulasi ensemble mikrokanonikal mendapat hasil yang baik, dengan error semakin kecil jika time step kecil. Time step terbesar dengan error yang masih memenuhi spesifikasi adalah 0,02 dan dapat digunakan untuk simulasi selanjutnya. 2. Dalam simulasi ensemble mikrokanonikal, temperatur sistem sesaat dan rata-rata tidak pernah melebihi temperatur awal jika simulasi dimulai dengan posisi kubik berpotensi terendah. Temperatur akhir dan rata-rata sistem tidak dapat diprediksi dari temperatur awal. 3. Temperatur sistem yang terlalu tinggi akan menyebabkan adanya molekul yang masuk ke dalam jangkauan molekul lain, dan akan menimbulkan error yang besar. 4. Termostat Nose Hoover dan termostat Berendsen berhasil mengendalikan temperatur sistem sehingga berkisar pada temperatur yang diinginkan. 5. Parameter Q pada termostat Nose Hoover mempengaruhi frekuensi osilasi error temperatur. Pada termostat Berendsen, Q mempengaruhi besar error temperatur. 6. Termostat Nose Hoover tidak melakukan peredaman error temperatur sistem dari temperatur set. Termostat Berendsen lebih baik dalam mengendalikan temperatur. 7. Tekanan hampir seluruhnya bergantung pada temperatur, seperti gas ideal. Entalpi hampir seluruhnya bergantung pada tekanan.
6.2
SARAN
1. Pengujian simulasi dari penelitian ini dengan parameter-parameter spesifik untuk material tertentu, misal Argon.
66
2. Pengujian simulasi dengan jumlah molekul lebih banyak dengan bantuan pembuatan daftar tetangga (neighbour list) dan komputasi paralel. 3. Modifikasi simulasi agar dapat melakukan pengendalian tekanan. 4. Modifikasi simulasi agar dapat melakukan simulasi dinamika molekul pada fasa cair dan gas. 5. Visualisasi tiga dimensi.
67
DAFTAR PUSTAKA 1. Ceperly,
D.,
Techniques,
Simulation
http://archive.ncsa.uiuc.edu/Apps/CMP/lectures/cms01/lectures.ppt 2. Haile, J.M., Molecular Dynamics Simulation: Elementary Methods, John Wiley & Sons, Inc., 1992 3. Cengel, Yunus A., dan Boles, Michael A., Thermodynamics: An Engineering Approach, McGraw-Hill Inc, 1994
4. Stote, R., Dejaegere, A., Kuznetsov, D., dan Falquet, L., CHARMM Molecular
Simulations,
Dynamics
http://www.ch.embnet.org/MD_tutorial/ 5. Ercolessi,
Furio,
A
Molecular
Dynamics
Primer,
http://www.sissa.it/furio/md 6. Refson,
Keith,
Moldy
User’s
Manual,
http://www.earth.ox.ac.uk/~keith/moldy-manual/moldy.html 7. Huang, Kerson, Statistical Mechanics, John Wiley & Sons, Inc., 1987 8. Nose, Shuichi, A Molecular Dynamics Method for Simulations in the Canonical Ensemble, Molecular Physics, 1984, Vol 52, No. 2, 255-268
9. Allen, M.P. and Tildesley, D.J., Computer Simulation of Liquids, Oxford University Press, 1987 10. Cornwell, C.F dan Wille, L.T., Parallel Molecular Dynamics Simulations for Short-Ranged Many-Body Potentials, Computer
Physics Communications 128 (2000), 477-491 11. Kruse, R. L., Leung, B. P. dan Tondo, C. L., Data Structures and Program Design in C, Prentice-Hall Inc., 1998
12. Fahmi, Mahuddin, Perancangan Perangkat Lunak Simulasi Dinamika Molekular dengan Model Potensial Lennard Jones, Tugas Akhir, 1999
13. Aitken, P. dan Jones, B., Teach Yourself C in 21 Days, SAMS Publishing, 1994 14. Lee, R.C. dan Tepfenhart, W. M., UML and C++, Prentice-Hall Inc., 1997
68
15. Nose, Shuichi, A Unified Formulation of the Constant Temperature Molecular Dynamics Methods, Journal of Chemical Physics 81, 1984
16. Reif, F., Fundamentals of Statistical and Thermal Physics, McGrawHill Publishing Company, 1965 17. Reisdorph, Kent, Teach Yourself Borland C++ Builder in 21 Days, SAMS Publishing, 1997 18. Steinbach, Peter J., Introduction to Macromolecular Simulation, http://cmm.cit.nih.gov/intro_simulation 19. Walpole, Ronald E. dan Myers, Raymod H., Probability and Statistics for Engineers and Scientists, Macmillan Publishing Company, 1993
20. Woodcock, L.V., Isothermal Molecular Dynamics Calculations for Liquid Salts, Chemical Physics Letters, Agustus 1971, Vol 10, No. 3
21. Stadler, J., Mikulla, R. dan Trebin, H.-R., IMD: A Software Package for Molecular Dynamics Studies on Parallel Computers, International
Journal of Modern Physics C, Vol. 8, No. 5 (1997)
69
LAMPIRAN A LISTING PROGRAM MainForm.h //-------------------------------------------------------------------------#ifndef MainFormH #define MainFormH //-------------------------------------------------------------------------#include
#include #include <StdCtrls.hpp> #include #include <ExtCtrls.hpp> #include #include <Series.hpp> #include #include #include //-------------------------------------------------------------------------class TForm1 : public TForm { __published: // IDE-managed Components TPanel *Panel1; TPanel *Panel2; TLabel *Label5; TEdit *Time; TButton *ButtonRun; TEdit *U; TEdit *K; TEdit *E; TEdit *error; TLabel *Label12; TLabel *Label13; TLabel *Label14; TLabel *Label15; TEdit *beta; TLabel *Label16; TEdit *Tavg; TLabel *Label18; TEdit *T; TLabel *Label20;
70
TEdit *Ham; TLabel *Label21; TEdit *P; TLabel *Label22; TLabel *Label23; TEdit *Pavg; TGroupBox *GroupBox1; TRadioButton *RadioButton64; TRadioButton *RadioButton125; TRadioButton *RadioButton216; TRadioButton *RadioButton343; TRadioButton *RadioButton512; TGroupBox *GroupBox2; TLabel *Label2; TEdit *SIG; TLabel *Label3; TEdit *EPS; TLabel *Label4; TEdit *L0; TLabel *Label6; TEdit *DT; TLabel *Label7; TEdit *TMAX; TLabel *Label8; TEdit *RC; TLabel *Label19; TEdit *Kb; TLabel *Label9; TEdit *T0; TGroupBox *GroupBox3; TLabel *Label10; TEdit *Q; TLabel *Label11; TEdit *TSet; TLabel *Label17; TEdit *TSetStart; TRadioButton *RadioButtonNone; TRadioButton *RadioButtonNoseHoover; TRadioButton *RadioButtonBerendsen; TButton *ButtonReset; TButton *ButtonStop; TCheckBox *DisableLennardJones; TPanel *PanelEnergy; TChart *ChartPotential; TFastLineSeries *Series1; TChart *ChartKinetic; TFastLineSeries *Series2; TButton *Button1; TButton *Button2;
71
TButton *Button3; TButton *Button4; TButton *Button5; TPanel *PanelTemperature; TChart *ChartTemperature; TChart *ChartError; TFastLineSeries *Series3; TFastLineSeries *Series9; TChart *ChartZeta; TChart *ChartZetaDot; TFastLineSeries *Series10; TFastLineSeries *Series11; TChart *Chart1; TFastLineSeries *Series4; TButton *Button6; TPanel *PanelPressure; TChart *ChartPressure; TFastLineSeries *Series17; TChart *ChartEnthalpy; TFastLineSeries *Series18; TPanel *PanelPosition; TChart *ChartPosition; TPointSeries *Series21; TRadioButton *RadioButtonXY; TRadioButton *RadioButtonYZ; TRadioButton *RadioButtonZX; TPanel *PanelHamiltonian; TChart *ChartHamiltonian; TChart *ChartS; TFastLineSeries *Series13; TChart *ChartSDot; TPanel *PanelMomentum; TChart *ChartMomentumX; TChart *ChartMomentumY; TChart *ChartMomentumZ; TFastLineSeries *Series5; TFastLineSeries *Series6; TFastLineSeries *Series7; TChart *ChartPf; TChart *ChartPm; TFastLineSeries *Series19; TFastLineSeries *Series20; TChart *ChartHamiltonianError; TLabel *Label1; TEdit *HamError; TFastLineSeries *Series14; TFastLineSeries *Series15; TFastLineSeries *Series16;
72
void __fastcall ButtonRunClick(TObject *Sender); void __fastcall ButtonStopClick(TObject *Sender); void __fastcall ButtonResetClick(TObject *Sender); void __fastcall TSetChange(TObject *Sender); void __fastcall Button1Click(TObject *Sender); void __fastcall Button3Click(TObject *Sender); void __fastcall Button6Click(TObject *Sender); void __fastcall Button5Click(TObject *Sender); void __fastcall Button4Click(TObject *Sender); void __fastcall Button2Click(TObject *Sender); private: // User declarations public: // User declarations __fastcall TForm1(TComponent* Owner); void Output(); bool stop; }; //-------------------------------------------------------------------------extern PACKAGE TForm1 *Form1; //-------------------------------------------------------------------------#endif
73
MainForm.cpp //-------------------------------------------------------------------------#include #pragma hdrstop #include #include #include #include
"MainForm.h" "Simulation.h" "Property.h" "Parameter.h"
//-------------------------------------------------------------------------#pragma package(smart_init) #pragma resource "*.dfm" TForm1 *Form1; Simulation *Sim = new Simulation; //-------------------------------------------------------------------------__fastcall TForm1::TForm1(TComponent* Owner) : TForm(Owner) { stop = false; } //-------------------------------------------------------------------------void __fastcall TForm1::ButtonRunClick(TObject *Sender) { ButtonRun->Enabled = false; ButtonStop->Enabled = true; Sim->Run(); ButtonStop->Enabled = false; ButtonReset->Enabled = true;
} //-------------------------------------------------------------------------void __fastcall TForm1::ButtonStopClick(TObject *Sender) { stop = 1; } //--------------------------------------------------------------------------
74
void __fastcall TForm1::ButtonResetClick(TObject *Sender) { ButtonReset->Enabled = false; Sim->Reset(); stop = 0; Time->Text = U->Text = K->Text = E->Text = error->Text = beta->Text = T->Text = Ham->Text = HamError->Text P->Text = Tavg->Text = Refresh();
0; 0; 0; 0; 0; 0; 0; 0; = 0; 0; 0;
while (Series1->Count() != 0) { Series1->Delete(0); Series2->Delete(0); Series3->Delete(0); Series4->Delete(0); Series5->Delete(0); Series6->Delete(0); Series7->Delete(0); //Series8->Delete(0); Series9->Delete(0); Series10->Delete(0); Series11->Delete(0); //Series12->Delete(0); Series13->Delete(0); Series14->Delete(0); Series15->Delete(0); Series16->Delete(0); Series17->Delete(0); Series18->Delete(0); Series19->Delete(0); Series20->Delete(0); } while (Series21->Count() != 0) { Series21->Delete(0);
75
} ButtonRun->Enabled = true; } //-------------------------------------------------------------------------void TForm1::Output() { Time->Text = Sim->t; U->Text = Prop->U; K->Text = Prop->K; E->Text = Prop->E; error->Text = Prop->error; beta->Text = Prop->beta; T->Text = Prop->T; Ham->Text = Prop->Ham; HamError->Text = Prop->HamError; P->Text = Prop->P; double a = StrToFloat(Tavg->Text); Tavg->Text = (Sim->t*a + Prop->T)/(Sim->t+1); double b = StrToFloat(Pavg->Text); Pavg->Text = (Sim->t*b + Prop->P)/(Sim->t+1); if (Sim->t%(Par->TMAX/500) == 0) { Series1->AddXY(Sim->t*Par->DT, >U/Par->N, ' ' , clBlue ); Series2->AddXY(Sim->t*Par->DT, >K/Par->N * pow(10,3), ' ' , clBlue ); Series3->AddXY(Sim->t*Par->DT, >E/Par->N, ' ' , clBlue ); Series4->AddXY(Sim->t*Par->DT, >error * pow(10,9), ' ' , clBlue );
PropPropPropProp-
Series5->AddXY(Sim->t*Par->DT, Prop>momentum[0] * pow(10,9), ' ' , clBlue ); Series6->AddXY(Sim->t*Par->DT, Prop>momentum[1] * pow(10,9), ' ' , clBlue ); Series7->AddXY(Sim->t*Par->DT, Prop>momentum[2] * pow(10,9), ' ' , clBlue ); Series9->AddXY(Sim->t*Par->DT, Prop->T, ' ' , clBlue ); Series10->AddXY(Sim->t*Par->DT, Prop>zeta, ' ' , clBlue ); Series11->AddXY(Sim->t*Par->DT, Prop>zetadot, ' ' , clBlue );
76
Series13->AddXY(Sim->t*Par->DT, >Ham, ' ' , clBlue ); Series14->AddXY(Sim->t*Par->DT, >HamError * pow(10,3), ' ' , clBlue ); Series15->AddXY(Sim->t*Par->DT, >s, ' ' , clBlue ); Series16->AddXY(Sim->t*Par->DT, >sdot, ' ' , clBlue );
Prop-
Series17->AddXY(Sim->t*Par->DT, * pow(10,6), ' ' , clBlue ); Series18->AddXY(Sim->t*Par->DT, >H, ' ' , clBlue ); Series19->AddXY(Sim->t*Par->DT, >Pf * pow(10,6), ' ' , clBlue ); Series20->AddXY(Sim->t*Par->DT, >Pm * pow(10,6), ' ' , clBlue ); }
Prop->P
PropPropProp-
PropPropProp-
if (ChartPosition->Visible == true) { while (Series21->Count() != 0) Series21->Delete(0); for (Molecule *p = head; p != NULL; p = p>next)
{
, clBlue);
if (RadioButtonXY->Checked == true) Series21->AddXY(p->R[0], p->R[1], ' ' else if (RadioButtonYZ->Checked == true) Series21->AddXY(p->R[1], p->R[2], ' '
, clBlue);
else if (RadioButtonZX->Checked == true) Series21->AddXY(p->R[2], p->R[0], ' '
, clBlue);
} } Refresh();
} //-------------------------------------------------------------------------//--------------------------------------------------------------------------
77
void __fastcall TForm1::TSetChange(TObject *Sender) { if (TSet->Text != "") Par->Tset = StrToFloat(TSet->Text); } //-------------------------------------------------------------------------void __fastcall TForm1::Button1Click(TObject *Sender) { PanelEnergy->Visible = true; PanelMomentum->Visible = false; PanelTemperature->Visible = false; PanelHamiltonian->Visible = false; PanelPressure->Visible = false; PanelPosition->Visible = false; } //-------------------------------------------------------------------------void __fastcall TForm1::Button2Click(TObject *Sender) { PanelEnergy->Visible = false; PanelMomentum->Visible = true; PanelTemperature->Visible = false; PanelHamiltonian->Visible = false; PanelPressure->Visible = false; PanelPosition->Visible = false; } //-------------------------------------------------------------------------void __fastcall TForm1::Button3Click(TObject *Sender) { PanelEnergy->Visible = false; PanelMomentum->Visible = false; PanelTemperature->Visible = true; PanelHamiltonian->Visible = false; PanelPressure->Visible = false; PanelPosition->Visible = false; } //-------------------------------------------------------------------------void __fastcall TForm1::Button4Click(TObject *Sender) { PanelEnergy->Visible = false; PanelMomentum->Visible = false;
78
PanelTemperature->Visible PanelHamiltonian->Visible PanelPressure->Visible PanelPosition->Visible
= = = =
false; true; false; false;
} //-------------------------------------------------------------------------void __fastcall TForm1::Button5Click(TObject *Sender) { PanelEnergy->Visible = false; PanelMomentum->Visible = false; PanelTemperature->Visible = false; PanelHamiltonian->Visible = false; PanelPressure->Visible = true; PanelPosition->Visible = false; } //-------------------------------------------------------------------------void __fastcall TForm1::Button6Click(TObject *Sender) { PanelEnergy->Visible = false; PanelMomentum->Visible = false; PanelTemperature->Visible = false; PanelHamiltonian->Visible = false; PanelPressure->Visible = false; PanelPosition->Visible = true; } //--------------------------------------------------------------------------
79
Molecule.h //-------------------------------------------------------------------------#ifndef MoleculeH #define MoleculeH class Molecule { public: Molecule(); int index; //molecule index number double m; //mass of molecule double R[3]; //position at t double Rnext[3]; //position at t+dt double Rprev[3]; //position at t-dt double v[3]; //velocity at t double F[3]; //force received at t Molecule* next; //pointer to next molecule }; //-------------------------------------------------------------------------#endif
80
Molecule.cpp //-------------------------------------------------------------------------#include #pragma hdrstop #include "Molecule.h" //-------------------------------------------------------------------------#pragma package(smart_init) Molecule::Molecule() { m = 1; next = NULL; }
81
Simulation.h //-------------------------------------------------------------------------#ifndef SimulationH #define SimulationH #include #include #include #include #include #include
<math.h> "MainForm.h" "Molecule.h" "LennardJones.h" "Property.h" "Parameter.h"
//-------------------------------------------------------------------------class Simulation { public: Simulation(); int t; void Run(); void Reset(); private: void NextStep(); }; extern extern extern extern
Molecule *head; LennardJones *LJ; Property *Prop; Parameter *Par;
#endif
82
Simulation.cpp //-------------------------------------------------------------------------#include #pragma hdrstop #include "Simulation.h" #include "Verlet.h" #include "Initialize.h" //-------------------------------------------------------------------------#pragma package(smart_init) Initialize *Init = new Initialize; Verlet *Ver = new Verlet; Molecule *head = NULL; LennardJones *LJ = new LennardJones; Property *Prop = new Property; Parameter *Par = new Parameter; Simulation::Simulation() { t = 0; } void Simulation::Run() { for (t = 0; t < 1; t++) { Init->Run(); LJ->Run(); Ver->RunOne(); } Form1->ChartPosition->LeftAxis->Minimum = >L/2; Form1->ChartPosition->LeftAxis->Maximum = Form1->ChartPosition->BottomAxis->Minimum >L/2; Form1->ChartPosition->BottomAxis->Maximum >L/2;
-ParPar->L/2; = -Par= Par-
for (t = 1; t < (Par->TMAX+1); t ++) { NextStep();
83
LJ->Run(); Ver->Run(); Form1->Output(); Application->ProcessMessages(); if (Form1->stop) break; }
}
void Simulation::Reset() { t = 0; while (head != NULL) { Molecule* p = head; head = p->next; delete (p); } } void Simulation::NextStep() { for (Molecule *p = head; p != NULL; p = p->next) { for (int i = 0; i < 3; i++) { p->Rprev[i] = p->R[i]; p->R[i] = p->Rnext[i]; } } }
84
Parameter.h //-------------------------------------------------------------------------#ifndef ParameterH #define ParameterH //-------------------------------------------------------------------------class Parameter { public: Parameter(); void Set(); int N; //number of molecules int NAxis; //number of molecules along length of cube double SIG; //length parameter of interaction double EPS; //strength of interaction double L0; //length of molecular lattice double L; //length of simulation box double RC; //cutoff distance double DT; //iteration time step double T0; //average initial velocity of each molecule double Kb; //Boltzmann constant int TMAX; //maximum iteration number double Q; //Nose-Hoover thermostat parameter double Tset; //temperature set point }; //-------------------------------------------------------------------------#endif
85
Parameter.cpp //-------------------------------------------------------------------------#include #pragma hdrstop #include "Parameter.h" #include "Simulation.h" //-------------------------------------------------------------------------#pragma package(smart_init) Parameter::Parameter() { } void Parameter::Set() { if (Form1->RadioButton64->Checked) { N = 64; NAxis = 4; } if (Form1->RadioButton125->Checked) { N = 125; NAxis = 5; } if (Form1->RadioButton216->Checked) { N = 216; NAxis = 6; } if (Form1->RadioButton343->Checked) { N = 343; NAxis = 7; } if (Form1->RadioButton512->Checked) { N = 512; NAxis = 8; } SIG = StrToFloat(Form1->SIG->Text); EPS = StrToFloat(Form1->EPS->Text); L0 = StrToFloat(Form1->L0->Text)*SIG; L = NAxis * L0;
86
DT = StrToFloat(Form1->DT->Text); TMAX = StrToFloat(Form1->TMAX->Text); RC = StrToFloat(Form1->RC->Text)*L0; T0 = StrToFloat(Form1->T0->Text); Kb = StrToFloat(Form1->Kb->Text); Q = StrToFloat(Form1->Q->Text); Tset = StrToFloat (Form1->TSet->Text); }
87
Property.h //-------------------------------------------------------------------------#ifndef PropertyH #define PropertyH //-------------------------------------------------------------------------class Property { public: Property(); void Run(); double U; double K; double E; double E0; double momentum[3]; double error; double zeta; double zetadot; double beta; double s; double sdot; double T; double P; double Pf; double Pm; double H; double Ham; double Ham0; double HamError; private: void Kinetic(); void Energy(); }; //-------------------------------------------------------------------------#endif
88
Property.cpp //-------------------------------------------------------------------------#include #pragma hdrstop #include "Property.h" #include "Simulation.h" //-------------------------------------------------------------------------#pragma package(smart_init) Property::Property() { U = 0; K = 0; E = 0; E0 = 0; for (int i = 0; i<3; i++) momentum[i] = 0; error = 0; zeta = 0; zetadot = 0; beta = 0; s = 1; sdot = 0; T = 0; P = 0; H = 0; Ham = 0; } void Property::Run() { Kinetic(); Energy(); T = 2 * Prop->K /(3*Par->N*Par->Kb); Pm = 2*Par->N/(3*pow(Par->L,3))*Prop->K; P = Pf+Pm; H = E + P * pow(Par->L,3); if (s>0) Ham = U + K*pow(s,2) + 0.5*Par->Q*pow(sdot,2) + 3*Par->Kb * T *log(s);
89
} void Property::Kinetic() { K = 0; for (int i = 0; i<3; i++) momentum[i] = 0; for (Molecule *p = head; p != NULL; p = p->next) for (int i = 0; i < 3; i++) { K += 0.5 * p->m * pow(p->v[i],2); momentum[i] += p->m * p->v[i]; } } void Property::Energy() { E = U + K; if (E0 != 0) error = (E-E0)/E0; if (Ham0 != 0) HamError = (Ham-Ham0)/Ham0; }
90
Initialize.h //-------------------------------------------------------------------------#ifndef InitializeH #define InitializeH //-------------------------------------------------------------------------class Initialize { public: Initialize(); void Run(); private: void CreateMolecules(); void InitialPosition(); void InitialVelocity(); }; //-------------------------------------------------------------------------#endif
91
Initialize.cpp //-------------------------------------------------------------------------#include #pragma hdrstop #include <stdlib.h> #include <math.h> #include #include #include #include #include
"Initialize.h" "Parameter.h" "MainForm.h" "LennardJones.h" "Property.h"
#include "Simulation.h" //-------------------------------------------------------------------------#pragma package(smart_init) //-------------------------------------------------------------------------Initialize::Initialize() { } void Initialize::Run() { Par->Set(); CreateMolecules(); InitialPosition(); InitialVelocity(); } void Initialize::CreateMolecules() { Molecule* p; for (int i = 0; i < Par->N; i++) { if ((p = new Molecule) != 0) { p->index = Par->N - i; p->next = head; head = p; } } }
92
void Initialize::InitialPosition() { //CUBE FORMATION double x0 = 0.5*(1-Par->NAxis)*Par->L0; double y0 = x0; double z0 = x0; Molecule* p = head; double z = z0; for (int i = 0; i < Par->NAxis; i++) { double y = y0; for (int j = 0; j < Par->NAxis; j++) { double x = x0; for (int k = 0; k < Par->NAxis; k++) { p->R[0] = x; p->R[1] = y; p->R[2] = z; p = p->next; x += Par->L0; } y += Par->L0; } z += Par->L0; } } void Initialize::InitialVelocity() { double K = 0; double T = 0; Randomize(); for (int i = 0; i < 3; i++) { double V = 0;
>next)
//ASSIGN RANDOM VALUES for (Molecule* p = head; p != NULL; p = p{ p->v[i] = 1; for (int j = 0; j < 200; j++) p->v[i] += 0.01*(random(2)-0.5); //WEIGHTED RANDOM NUMBER
93
if (random(2) == 0) p->v[i] = -p->v[i]; //DIRECTION V += p->v[i]; K += 0.5*p->m*pow(p->v[i],2); } //READJUSTING MOMENTUM ON AXIS for (Molecule* p = head; p != NULL; p = p>next) }
p->v[i] -= V/Par->N;
T = 2*K/(3*Par->N*Par->Kb); for (int i = 0; i < 3; i++) { for (Molecule* p = head; p != NULL; p = p>next) { p->v[i] *= sqrt(Par->T0/T); } } }
94
Verlet.h //-------------------------------------------------------------------------#ifndef VerletH #define VerletH class Verlet { public: Verlet(); void Run(); Simulation void RunOne(); private: //double zeta; double zetaprev; double zetanext; //double beta; double sprev; double snext; void NextStep(); void PredictPosition(); void PredictPositionOne(); void Velocity(); void PredictZeta(); void PredictS(); void VelocityScale(); };
//Run the
//-------------------------------------------------------------------------#endif
95
Verlet.cpp //-------------------------------------------------------------------------#include #pragma hdrstop #include <math.h> #include #include #include #include #include
"Verlet.h" "Parameter.h" "MainForm.h" "LennardJones.h" "Property.h"
#include "Simulation.h" //-------------------------------------------------------------------------#pragma package(smart_init) Verlet::Verlet() { zetanext = 0; zetaprev = 0; snext = 1; sprev = 1; } void Verlet::Run() { NextStep(); PredictPosition(); Velocity(); Prop->Run();
}
PredictZeta(); PredictS(); Prop->zetadot = (zetanext-zetaprev)/(2*Par->DT);
void Verlet::RunOne() { PredictPositionOne(); Prop->Run(); Prop->E0 = Prop->E; Prop->Ham0 = Prop->Ham; } void Verlet::NextStep()
96
{
zetaprev = sprev = Prop->zeta Prop->s Prop->beta Prop->sdot
Prop->zeta; Prop->s; = zetanext; = snext; = 0.5 * Prop->zeta * Par->DT; = Prop->zeta * Prop->s;
} void Verlet::PredictPosition() { for (Molecule *p = head; p != NULL; p = p->next) { for (int i = 0; i < 3; i++) { if (Form1->RadioButtonNone->Checked ==
true)
p->Rnext[i] = 2 * p->R[i] - p->Rprev[i] + p->F[i]/p->m * pow(Par->DT,2); else p->Rnext[i] = (2 * p->R[i] - (1-Prop>beta) * p->Rprev[i] + p->F[i]/p->m * pow(Par>DT,2))/(1+Prop->beta); if (p->Rnext[i] > 0.5*Par->L) { p->Rnext[i] -= Par->L; p->R[i] -= Par->L; p->Rprev[i] -= Par->L; } if (p->Rnext[i] < -0.5*Par->L) { p->Rnext[i] += Par->L; p->R[i] += Par->L; p->Rprev[i] += Par->L; } }
}
} void Verlet::PredictPositionOne() { for (Molecule *p = head; p != NULL; p = p->next) for (int i = 0; i < 3; i++) p->Rnext[i] = p->R[i] + p->v[i] * Par->DT + p->F[i]/p->m * pow(Par->DT,2); }
97
void Verlet::Velocity() { for (Molecule *p = head; p != NULL; p = p->next) { for (int i = 0; i < 3; i++) { p->v[i] = (p->Rnext[i] - p->Rprev[i]) / (2 * Par->DT); } } } void Verlet::PredictZeta() { if (Form1->RadioButtonBerendsen->Checked == true) zetanext = (Prop->T-Par->Tset)/(2*Par->Q*Prop>T); else if (Form1->RadioButtonNoseHoover->Checked == true) zetanext = zetaprev + 2 * Par->DT * ((2*Prop>K)-(3*Par->N*Par->Kb*Par->Tset))/ Par->Q; } void Verlet::PredictS() { snext = sprev + 2 * Prop->sdot*Par->DT; } void Verlet::VelocityScale() { for (int i = 0; i < 3; i++) { for (Molecule* p = head; p != NULL; p = p>next) { p->v[i] *= sqrt(Par->Tset/Prop->T); p->Rnext[i] = p->Rprev[i] + p->v[i] * 2 * Par->DT; } } }
98
LennardJones.h #ifndef LennardJonesH #define LennardJonesH #include "Molecule.h" class LennardJones { public: LennardJones(); void Run(); private: double RAxis[3], Rij; double force[3]; void ClearValues(); void Distance (Molecule *p1, Molecule *p2); void Force(Molecule *p1, Molecule *p2); void Potential(Molecule *p1, Molecule *p2); void PressureByForce(Molecule *p1, Molecule *p2); };
//-------------------------------------------------------------------------#endif
99
LennardJones.cpp //-------------------------------------------------------------------------#include #pragma hdrstop #include <math.h> #include "LennardJones.h" #include "Simulation.h" //-------------------------------------------------------------------------#pragma package(smart_init) LennardJones::LennardJones() { Rij = 0; for (int i = 0; i < 3; i++) { RAxis[i] = 0; force[i] = 0; } } void LennardJones::Run() { ClearValues(); for (Molecule *p1 = head; p1 != NULL; p1 = p1>next) { for (Molecule *p2 = p1->next; p2 != NULL; p2 = p2->next) { Distance(p1,p2); if (Rij < Par->RC) { if (Form1->DisableLennardJones->Checked == false) { Force(p1,p2); Potential(p1,p2); PressureByForce(p1,p2); } } } } } void LennardJones::ClearValues()
100
{
Prop->U = 0; Prop->Pf = 0; for (Molecule *p = head; p != NULL; p = p->next) for (int i = 0; i < 3; i++) p->F[i] = 0;
} void LennardJones::Distance (Molecule *p1, Molecule *p2) { for (int i=0; i<3; i++) { RAxis[i] = p1->R[i] - p2->R[i]; if (RAxis[i] > 0.5*Par->L) RAxis[i] -= Par->L; if (RAxis[i] < -0.5*Par->L) RAxis[i] += Par->L; } Rij = sqrt(pow(RAxis[0],2)+pow(RAxis[1],2)+pow(RAxis[2],2)); } void LennardJones::Force (Molecule *p1, Molecule *p2) { for (int i=0; i<3; i++) { force[i] = 24 * Par->EPS * RAxis[i] * pow(Par>SIG,6) / pow(Rij,8) * (2 * pow(Par->SIG/Rij,6) - 1); p1->F[i] += force[i]; p2->F[i] -= force[i]; //Newton's Third Law } } void LennardJones::Potential (Molecule *p1, Molecule *p2) { Prop->U += 4 * Par->EPS *(pow(Par->SIG/Rij,12)pow(Par->SIG/Rij,6)); } void LennardJones::PressureByForce(Molecule *p1, Molecule *p2) { for (int i=0; i<3; i++) Prop->Pf += force[i] * RAxis[i]; Prop->Pf = Prop->P/(3*pow(Par->L,3));
101
}
102
