Mata Kuliah Simulasi dan Pemodelan
ANALISIS PROGRAM MICROCANONICAL MOLECULAR – DYNAMICS PROGRAM
Dosen
: Dr. A. Benny Mutiara
Oleh
: Ernastuti & Nur Yuliani
Program Doktor Teknologi Informasi/Ilmu Komputer Universitas Gunadarma Maret 2006
Mata Kuliah Simulasi dan Pemodelan
Analisis Program Microcanonical Molecular – Dynamics Program Oleh : Ernastuti dan Nur Yuliani
Molecular Dynamics Metode molecular dinamis mengkomputasi phase spase trajectory dari suatu koleksi molekul yang secara individu mengikuti hukum klasik dari motion. Simulasi awal dilakukan system dengan energi yang konstan. Maka property dikalkulasikan dengan microcanonical esembel, dimana N adalah jumlah partikel, V adalah volume dan E adalah energi konstan. Kadang-kadang T adalah temperature konstan. Microcanonical Molecular-Dynamics adalah suatu algoritma dasar dari molecular dinamis dengan energi yang tetap. Point awal adalah Hamiltonian yang mendeskripsikan interaksi antara N partikel. Persamaan motion dengan mengunakan form Hamiltonian, adalah sebagai berikut:
1 1 2 pi + ∑ u (rij ) ∑ 2 i m i< j Dimana rij adalah jarak antara partikel i dan j, N adalah jumlah partikel konstan dan P H=
adalah constraint dari zero linear momentum. Persamaan motion dengan menngunakan form Newtonian, adalah sebagai berikut:
d 2 ri (t ) 1 = ∑ Fi (rij ) m i< j dt 2 Secara analisis, solusi system dari persamaan differensial order kedua didapat dengan melakukan integrasi dua kali dari waktu 0 ke waktu t, untuk mendapatkan velocity dan posisi. Simulasi komputer dari system molecular dapat dibagi menjadi 3 bagian, yaitu: 1. Initialitation Pada bagian ini, kondisi initial ditetapkan. Untuk memulai Algoritma dibutuhkan posisi dan velocity. Pada algoritma ini anggap saja initial posisi adalah lattice dan velocity diambil dari distibusi Boltzman. 2. Equilibration Equilibrium dapat dibangun jika system sudah ditetapkan untuk nilai energi kinetic dan energi potential tertentu. 3. Production Semua kuantitas dikomputasi bersama dengan trajectory dari system pada phase space. Program diawali dengan pendeklarasian variable dan parameter yang digunakan dalam program. Lalu dilanjutkan dengan menyiapkan inisial konfigurasi .
Halaman 1
Mata Kuliah Basis Data dan Sistem berbasis Web
Algoritma pada phase initialitation 1. Spesifikasi posisi r io dan ri 1 2. Komputasi forces pada time step n : Fin. 3. Komputasi posisi pada time step n+1 seperti pada ri n+1 4. Komputasi velocity pada time step n seperti pada vi n Algoritma pada phase equilibrium 1. Integrasikan persamaan dari motion untuk beberapa time step 2. komputasi energi kinetik dan energi potensial. 3. jika energi tidak sama dengan yang diharapkan, maka lakukan skala pada velocity 4. Ulangi dari langkah 1 sampai system mencapai equilibrium,
Kesimpulan Program Microcanonical Molecular-Dynamics yang tertera dibawah mensimulasikan system Lennard-Jones dengan N partikel pada suatu sel MD dengan volume V = σL’ dan energi tetap. Semua parameter diasumsikan berada pada form skala. Untuk posisi selanjutnya dan untuk komputasi velocity, program menggunakan form velocity dari algoritma Verlet. Pada tahap initial, partikel diletakkan pada suatu face-centered-cubic lattice. Jumlah partikel yang dibutuhkan adalah multiple 4. Velocity dipilih dari Distribusi Gaussian. Agar energi yang diinginkan dapat tercapai atau ekuivalen mendekati mean temperature, velocity diskala secara periodik sampai energi masuk kedalam interval yang diterima. Sistem tidak dapat diubah setelah phase equilibrium. Program ini juga dapat digunakan untuk simulasi isokinetik molecular dinamik dengan mensetting parameter IREP sama dengan unity dan ISTP lebih besar dari MD step.
Listing Program 1. ! Microcanonical Molecular-Dynamics Program ! !=================================================================== ! MOLECULAR DYNAMICS ! ! MICROCANONICAL ENSEMBLE ! APPLICATION TO ARGON. THE LENNARD JONES POTENTIAL IS TRUNCATED ! AT RCOFF AND NOT SMOOTHLY CONTINUED TO ZERO. INITIALLY THE ! NPART PARTICLES ARE PLACED ON AN FCC LATICE. THE VELOCITIES ! ARE DRAWN FROM A BOLTZMANN DISTRIBUTION WITH TEMPERATURE TREF. ! ! INPUT PARAMETER ARE AS FOLLOWS ! ! NPART NUMBER OF PARTICLES (MUST BE A MULTIPLE OF 4) ! SIDE SIDE LENGTH OF THE CUBICAL BOX IN SIGMA UNITS ! TREFF REDUCED TEMPERATURE ! RCOFF CUTOFF OF THE POTENTIAL IN SIGMA UNITS ! H BASIC TIME STEP ! IREP VELOCITIES SCALING EVERY IREP'TH TIME STEP ! ISTOP STOP OF SCALING EVERY IREP'TH TIME STEP ! TIMEMX NUMBER OF INTEGRATION STEPS ! ISEED SEED FOR TH
Halaman 2
Mata Kuliah Simulasi dan Pemodelan
! !=================================================================== ! REAL X(1:786), VH (1:786),F(786) ! REAL FY(1:256),FZ(1:256),Y(1:256), Z(1:256) ! REAL H,HSQ,HSQ2,TREF REAL KX,KY,KZ ! INTEGER CLOCK,TIMEMX ! EQUIVALENCE (FY,F(257)),(FZ,F(513)),(Y,X(257)),(Z,X(513)) ! !--------------------------------------------------------------------! ! DEFINITION OF THE SIMULATION PARAMETERS ! !--------------------------------------------------------------------! NPART = 256 SIDE = 6.75284 TREF = 0.722 THEN = 0.83134 RCOFF = 2.5 H = 0.064 IREF = 10 TIMEMX = 3 ISEED = 4711 ! !--------------------------------------------------------------------! ! SET THE OTHER PARAMETERS ! !--------------------------------------------------------------------! WRITE (6,*) 'MOLECULAR DYNAMICS SIMULATON PROGRAM' WRITE (6,*) '------------------------------------' ! WRITE (6,*) WRITE (6,*) 'NUMBER OF PARTICLES IS ',NPART WRITE (6,*) 'SIDE LENGTH OF THE BOX IS ',SIDE WRITE (6,*) 'CUT OFF IS ',RCOFF WRITE (6,*) 'REDUCED TEMPERATURE IS ',TREF WRITE (6,*) 'BASIC TIME STEP IS ',H WRITE (6,*) ! A = SIDE /4.0 SIDEH = SIDE * 0.5 HSQ = H * H HSQ2 = HSQ * 0.5 NPARTM = NPART - 1 RCOFFS = RCOFF * RCOFF TSCALE = 16.0 / (1.0 * NPART - 1.0) VAVER = 1.13 * SQRT (TREF / 24.0) ! N3 = 3 * NPART IOF1 = NPART IOF2 = 2 * NPART !
Halaman 3
Mata Kuliah Basis Data dan Sistem berbasis Web
CALL RANSET (ISEED) ! ! !--------------------------------------------------------------------! ! THIS PART OF THE PROGRAM PREPARES THE INITIAL CONFIGURATION ! !--------------------------------------------------------------------! ! SET UPFCC LATTICE FOR TE ATOMS INSIDE THE BOX ! IJK = 0 DO 10 LG = 0,1 DO 10 I = 0,3 DO 10 J = 0,3 DO 10 K = 0,3 IJK = IJK + 1 X(IJK ) = I * A + LG * A * 0.5 Y(IJK ) = J * A + LG * A * 0.5 Z(IJK ) = K * A ! 10 CONTINUE DO 15 LG = 1,2 DO 15 I=0,3 DO 15 J=0,3 DO 15 K=0,3 IJK = IJK + 1 X(IJK ) = I * A + (2-LG) * A * 0.5 Y(IJK ) = J * A + (LG-1) * A * 0.5 Z(IJK ) = K * A + A * 0.5 ! 15 CONTINUE ! ! ASSIGN VELOCITIES DISTRIBUTED NORMALLY ! ====================================== ! CALL MXWELL(VH,N3,H,TREFF) ! DO 50 I=1,N3 F(I) = 0.0 ! 50 CONTINUE ! !--------------------------------------------------------------------! ! START OF THE ACTUAL MOLECULAR DYNAMICS PROGRAM ! THE EQUATION OF MOTION ARE INTEGRATED USING THE 'SUMMED FORM' ! EVERY IREP'TH STEP VELOCITIES ARE RESCALED SO AS TO GIVE ! THE SPECIFIED TEMPERATURE TREF ! !--------------------------------------------------------------------! DO 200 CLOCK=1,TIMEMX ! 210 !
==ADVANCE POSITIONS ONE BASIC TIME STEP== DO 210 I=1,N3 X(I) = X(I) + VH(I) + F (I) CONTINUE
Halaman 4
Mata Kuliah Simulasi dan Pemodelan
!
215 ! 220 ! ! !
225 ! !
! ! !
!
270 ! ! ! ! ! ! !
==APPLY PERIODIC BOUNDARY CONDITION== DO 215 I=1,N3 IF (X(I) .LT. 0) X(I) = X(I) + SIDE IF (X(I) .GT. SIDE) X(I) = X(I) - SIDE CONTINUE ==COMPUTE THE PARTIAL VELOCITIES== DO 220 I=1,N3 VH(I) = VH(I) + F(I) CONTINUE == THIS PART COMPUTES THE FORCES ON THE PARTICLES == VIR = 0.0 EPOT = 0.0 DO 225 I=1,N3 F(I) = 0.0 CONTINUE DO 270 I=1,NPART XI = X(I) YI = Y(I) ZI = Z(I) DO 270 J=I+1,NPART XX = XI - X(J) YY = YI - Y(J) ZZ = ZI - Z(J) IF (XX IF (XX
.LT. .GT.
-SIDEH) XX SIDEH) XX
= XX = XX
+ SIDE - SIDE
IF (YY IF (YY
.LT. .GT.
-SIDEH) YY SIDEH) YY
= YY = YY
+ SIDE - SIDE
(ZZ .LT. -SIDEH) ZZ = ZZ (ZZ .GT. SIDEH) ZZ = ZZ = XX * XX + YY * YY + ZZ * ZZ ( RD .GT. RCOFF ) GOTO 270
+ SIDE - SIDE
IF IF RD IF
EPOT = EPOT + RD ** (-6.0) - RD ** (-3.0) R148 = RD ** (-7.0) - 0.5 * RD ** (-4.0) VIR = VIR - RD * R148 KX = XX * R148 F(I) = F(I) + KX F(J) = F(J) - KX KY = YY * R148 FY(I) = FY(I) + KY FY(J) = FY(J) - KY KZ = ZZ * R148 FZ(I) = FZ(I) + KZ FZ(J) = FZ(J) - KZ CONTINUE =================================================== = = = END OF THE FORCE CALCULATION = = = =================================================== ==COMPUTE THE VELOCITIES==
Halaman 5
Mata Kuliah Basis Data dan Sistem berbasis Web
300 ! !
305 ! !
306 ! ! !
310
! ! ! ! ! !
!
DO 300 I=1,N3 VH(I) = VH(I) + F(I) CONTINUE ==COMPUTE THE KINETIC ENERGY== EKIN = 0.0 DO 305 I=1,N3 EKIN = EKIN + VH(I) + VH(I) CONTINUE EKIN = EKIN / HSQ ==COMPUTE THE AVERAGE VELOCITY== VEL = 0.0 COUNT = 0.0 DO 306 I=1,NPART VX = VH(I) * VH(I) VY = VH(I+IOF1)* VH(I+IOF1) VZ = VH(I+IOF2)* VH(I+IOF2) SQ = SQRT ( VX + VY + VZ ) SQT = SQ / H IF (SQT .GT. VAVER ) COUNT = COUNT + 1 VEL = VEL + SQ CONTINUE VEL = VEL / H IF ( CLOCK .LT. ISTOP ) THEN ==NORMALIZE THE VELOCITIES TO OBTAIN THE== ==SPECIFIED REFERENCE TEMPERATURE == IF ( (CLOCK/IREP)*IREP .EQ. CLOCK) THEN WRITE(6,*) 'VELOCITY ADJUSTMENT' TS = TSCALE * EKIN WRITE(6,*) 'TEMPERATURE BEFORE SCALING IS ', TS SC = TREF / TS SC = SQRT( SC ) WRITE(6,*) 'SCALE FACTOR IS ', SC DO 310 I=1,N3 VH(I) = VH(I) * SC CONTINUE EKIN = TREF / TSCALE END IF END IF =================================================== = = = COMPUTE VARIOUS QUANTITIES = = = =================================================== IF ( (CLOCK/2)*2 .EQ. CLOCK ) THEN EK = 24.0 * EKIN EPOT = 4.0 * EPOT ETOT = EK + EPOT TEMP = TSCALE * EKIN PRES = DEN * 16.0 * ( EKIN - VIR ) / NPART VEL = VEL / NPART RP = (COUNT / 256.0) * 100.0 WRITE(6,6000) CLOCK,EK,EPOT,ETOT,TEMP,PRES,VEL,RP END IF
Halaman 6
Mata Kuliah Simulasi dan Pemodelan
200 6000 !
CONTINUE FORMAT(1I6,7F15.6) STOP END
! !===================================================================== = ! ! M X W E L L RETURN UPON CALLING A MAXWELL DISTRIBUTED DEVIATES ! FOR THE SPECIFIED TEMPERATURE TREF. ALL DEVIATED ARE SCALED BY ! A FACTOR. ! ! CALLING PARAMETERS ARE AS FOLLOW ! ! VH VECTOR OF LENGTH NPART ( MUST BE A MULTIPLE OF 2) ! NH3 VECTOR LENGTH ! H SCALING FACTOR WITH WHICH ALL ELEMENTS OF VH ARE ! MULTIFIED. ! TREF TEMPERATURE ! !===================================================================== = ! SUBROUTINE MXWELL (VH,N3,H,TREF) ! REAL VH(1:N3) ! REAL U1,U2,V1,V2,S,R ! NPART = N3 / 3 IOF1 = NPART IOF2 = 2 * NPART TSCALE = 16.0 / (1.0 * NPART - 1.0) ! DO 10 I=1,N3,2 ! ! 1 U1 = RAND() U2 = RAND() ! V1 = 2.0 * U1 - 1.0 V2 = 2.0 * U2 - 1.0 S = V1 * V1 + V2* V2 ! IF ( S .GE. 1.0 ) GOTO 1 R = -2.0 * ALOG(S) / S VH(I) = V1 * SQRT( R ) VH(I+1) = V2 * SQRT( R ) 10 CONTINUE ! EKIN = 0.0 SP = 0.0 DO 20 I=1,NPART SP = SP + VH(I) 20 CONTINUE SP = SP / NPART DO 21 I=1,NPART VH(I) = VH(I) - SP
Halaman 7
Mata Kuliah Basis Data dan Sistem berbasis Web
21
22
23
24
25
30
EKIN = EKIN + VH(I) * VH(I) CONTINUE WRITE(6,*) 'TOTAL LINEAR MOMENTUM IN X SP = 0.0 DO 22 I = IOF1 + 1,IOF2 SP = SP + VH(I) CONTINUE SP = SP / NPART DO 23 I=IOF1 + 1,IOF2 VH(I) = VH(I) - SP EKIN = EKIN + VH(I) * VH(I) CONTINUE WRITE(6,*) 'TOTAL LINEAR MOMENTUM IN Y SP = 0.0 DO 24 I = IOF2 + 1,N3 SP = SP + VH(I) CONTINUE SP = SP / NPART DO 25 I = IOF2 + 1,N3 VH(I) = VH(I) - SP EKIN = EKIN + VH(I) * VH(I) CONTINUE WRITE(6,*) 'TOTAL LINEAR MOMENTUM IN Z WRITE(6,*) 'VELOCITY ADJUSTMENT' TS = TSCALE * EKIN WRITE(6,*) 'TEMPERATURE BEFORE SCALING SC = TREF / TS SC = SQRT ( SC ) WRITE(6,*) 'SCALE FACTOR IS ' , SC SC = SC * H DO 30 I=1,N3 VH(I) = VH(I) * SC CONTINUE END
DIRECTION IS ' , SP
DIRECTION IS ' , SP
DIRECTION IS ' , SP IS ' , TS
Output Program 1. Molecular Dynamics Simulation Program
Halaman 8
