...

aree witoelar - Witoelar.com

by user

on
Category: Documents
0

views

Report

Comments

Transcript

aree witoelar - Witoelar.com
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−mm
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 <Classes.hpp>
#include <Controls.hpp>
#include <StdCtrls.hpp>
#include <Forms.hpp>
#include <ExtCtrls.hpp>
#include <Chart.hpp>
#include <Series.hpp>
#include <TeEngine.hpp>
#include <TeeProcs.hpp>
#include <ComCtrls.hpp>
//-------------------------------------------------------------------------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 <vcl.h>
#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 <vcl.h>
#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 <vcl.h>
#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 <vcl.h>
#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 <vcl.h>
#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 <vcl.h>
#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 <vcl.h>
#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 <vcl.h>
#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
Fly UP