Dynamika molekularna

(1)

Elementy dynamiki molekularnej

AB

(2)

(3)

Co to jest DM?

Dynamika molekularna (MD lub DM) zajmuje się czastkami masywnymi,

poruszającymi się pod działaniem sił newtonowskich. Pozwala symulować

układy cząstek w ciałach stałych, cieczach i gazach, opisywać ruch gwiazd i

galaktyk.

1

(4)

Co to jest DM?

W przypadku układów chaotycznych, jak np. oscylator Duffinga

m¨

x (t) = −γ ˙

x (t) + 2ax − 4bx

3

+ F

0

cos(ωt) .

rozwiązania numeryczne są praktycznie niewykonalne. Małe błedy numeryczne

w danych oraz błędy numerycznych zaokrągleń prowadzą do rozbieżności

(badanie wykładników Liapunowa).

(5)

Przestrzeń fazowa. Mapy Poincare’go.

Analizę wyników prowadzi się najczęściej w przestrzeni fazowej położeń (x(t)) i

prędkości (v (t)).

W przypadku zachowań chaotycznych przestrzeń fazową zastępuje się mapą

Poincare’go, która jest zbiorem punktów przestrzeni fazowej, odpowiadających

chwilom czasu t = 0, T , 2T , ..., gdzie T jest okresem związanym z

charakterystyczną częstością w układzie cząstek, np. częstością w członie

wymuszającym w przypadku oscylatora Duffinga, a więc ω. W przypadku

układów okresowych otrzymamy jeden punkt w przestrzeni fazowej. Dla

układów chaotycznych mogą to być dowolne zbiory punktów. Mogą pojawiać

się tzw. bifurkacje (rozgałęzienia). Układy z tzw. chaosem deterministycznym

mogą prowadzić do struktur fraktalnych, dziwnych atraktorów, których wymiar

(6)

Metoda Runge-Kutta I

dx

dt

=

v

(1)

dv

dx

=

−γv + 2ax − 4bx

3

+ F

0

cos(ωt)

(2)

Program całkuje układ równań ruchu metodą Runge’go-Kutty 4go rzędu:

~

k

1

=

h~

f (t, ~

x (t))

~

k

2

=

h~

f (t +

h

2 , ~

x (t) +

1

2 ~

k

1

)

~

k

3

=

h~

f (t +

h

2 , ~

x (t) +

1

2 ~

k

2

)

~

k

2

=

h~

f (t + h, ~

x (t) + ~

k

3

)

~

x (t + h)

=

~

x (t) +

1

6 (~

k

1

+ 2~

k

2

+ 2~

k

3

+ ~

k

4

) .

(3)

(7)

Program I

! Duffing Oscillator program duffing

implicit none

! real kind with decimal precision() >= 14 and exponent range() >= 300 integer, parameter :: dp = selected_real_kind(14, 300)

! constants

real(dp), parameter :: zero = 0._dp, one = 1._dp, two = 2._dp real(dp), parameter :: half = one / two, six = (one + two) * two ! physics parameters

real(dp) :: omega = 2.4_dp, gamma = 0.1_dp real(dp) :: a = 0.5_dp, b = 0.25_dp real(dp) :: F0 = 2.0_dp

! variables

real(dp), dimension(2) :: xv, dxvdt real(dp) :: t, t_max, dt

(8)

Program II

print *, ’Enter initial x and v: ’ read *, xv(1), xv(2)

print *, ’Enter time step dt’ read *, dt

print *, ’Enter integration time: ’ read *, t_max

print *, ’Enter output file name: ’ read *, file_name

open(unit=2, file=file_name) t = zero

write(2, *) t, xv(1), xv(2) ! main integration loop do

if (t >= t_max) exit ! 4th order Runge-Kutta step call find_dxvdt(t, xv, dxvdt) k1 = dt * dxvdt

xv_step = xv + k1 / two

call find_dxvdt(t + half * dt, xv_step, dxvdt) k2 = dt * dxvdt

xv_step = xv + k2 / two

call find_dxvdt(t + half * dt, xv_step, dxvdt) k3 = dt * dxvdt

xv_step = xv + k3

(9)

Program III

k4 = dt * dxvdt

xv = xv + (k1 + two * k2 + two * k3 + k4) / six t = t + dt

write(2, *) t, xv(1), xv(2) end do

print *, ’Output in file: ’, file_name close(2) contains ! Duffing equation subroutine find_dxvdt(t, xv, dxvdt) implicit none real(dp), intent(in) :: t

real(dp), dimension(2), intent(in) :: xv real(dp), dimension(2), intent(out) :: dxvdt dxvdt(1) = xv(2)

dxvdt(2) = - gamma * xv(2) + 2 * a * xv(1) - 4 * b * xv(1)**3 dxvdt(2) = dxvdt(2) + F0 * cos(omega * t)

(10)

Przykładowe wyniki

��

(11)

Przykładowe wyniki

��

(12)

Przykładowe wyniki

��

(13)

Przykładowa mapa Poincare’go

��

(14)

Podobna mapa dla wahadła tłumionego

3

2

1

0

1

2

3 Angle

1.0

0.5

0.0

0.5

1.0

1.5

2.0 AngularVelocity

d

2

θ/dt

2

+ sin(θ(t)) + β ˙

θ(t) + A ∗ sin(Ωt) = 0

(15)

Zadanie

Proszę dopisać w programie fragment obliczeń mapy Poincarego dl oscylatora

Duffinga. Narysować mapę dla kilku wybranych przypadków danych

początkowych. a) tylko punkty. b) punkty i linie. Porównać i wyciągnąć

wnioski. Wykonać to samo dla prostego oscylatora harmonicznego.

(16)

Algorytm Verleta I

Algorytm Verleta jest tzw. algorytmem symplektycznym i zachowuje energię

całkowitą cząstek spełniających równania ruchu Hamiltona dokładniej niż

algorytm Runge’go-Kutty.

Załóżmy, że jednowymiarowy ruch cząstki dany jest równaniwm F = ma.

Rozwińmy położenie x cząstki w szereg Taylora

x (t + h) = x (t) + v (t) +

1

2 a(t)h

2

+

1

6 a(t)h

˙

3

+ O(h

4

)

x (t − h) = x (t) − v (t) +

1

2 a(t)h

2

−

1

6 a(t)h

˙

3

+ O(h

4

)

Dodając te równania otrzymamy prosty algorytm Verleta:

x (t + h) = 2x (t) − x (t − h) + a(t)h

2

+ O(h

4

)

v (t) =

x (t + h) − x (t − h)

2h

+ O(h

2

(17)

Algorytm Verleta II

Nie jest to algorytm samostartujący (jak np. algorytm RK). Samostartujący

algorytm Verleta o podobnej dokładności jest

x (t + h)

=

x (t) + v (t)h +

1

2 a(t)h

2

+ O(h

4

)

v (t + h)

=

v (t) +

1

2 [a(t) + a(t + h)] h + O(h

2

)

(5)

Jest to algorytm dwukrokowy. Położenie x(t + h) wylicza się na podstawie

x (t), v (t) i a(t). Stąd a(t + h) i v (t + h).

(18)

Schemat algorytmu Verleta I

! samostartujący algorytm Verleta

!

subroutine sVerlet(dt, x, dxdt, d2xdt2, n, alpha)

dimension x(n), dxdt(n), d2dxt(n)

do i=1,n

x(i) = x(i) + dxdt(i) * dt + 0.5 * d2xdt2(i) * dt**2

dxdt(i) = dxdt(i) + 0.5 * d2xdt2(i) * dt

end do

compute_accelerations(x, d2xdt2, alpha, n)

do i=1,n

dxdt(i) += 0.5 * d2xdt2(i) * dt

end do

(19)

Schemat algorytmu Verleta II

subroutine compute_accelerations(x, d2xdt2, alpha, n)

dimension x(n), d2xdt2,

real alpha

do i = 2, n-1

dx_plus = x[i+1] - x[i]

dx_minus = x[i] - x[i-1]

d2xdt2[i] = dx_plus - dx_minus + &

alpha * (dx_plus2 - dx_minus2.0);

end do

(20)

(21)

Modelowanie argonu

Rozpatrujemy N atomów szlachetnego argonu (m ∼ 6, 7 × 10

−26

kg). Atomy w

przybliżeniu zachowują się jak sztywne kule przyciągające się siłami van der

Waalsa, które będziemy modelować potencjałem Lennarda-Jonesa’:

V (r ) = 4

_σ

r

12

−

σ

r

6

.

Tutaj r jest odległością centrów atomów, = 1, 65 × 10

−21

J,

a σ = 3, 4 × 10

−10

m odpowiada r dla którego energia jest zerem. Człon r

−12

opisuje odpychający twardy rdzeń, a r

−6

reprezentuje przyciąganie typu

dipol-dipol.

(22)

Potencjał L-J

-6 -4 -2 0 2 4 6 3 4 5 6 7 8 9 10 V(r) [meV] r [σ]

(23)

Potencjał L-J

Zadania

1. Obliczyć V w jego minimum r

min

.

2. Obliczyć siłę F działającą na cząstki umieszczone w potencjale L-J.

3. Narysować F (r ) dla 0 < r ≤ 10.

(24)

Przygotowanie programu MD (java) I

• Ustalenie jednostek. m = σ = = 1. Jednostką czasu jest w tym układzie

τ =

q

mσ2

= 2, 17 × 10

−12

s. Czas charakterystyczny jest więc rzędu

picosekund.

• Ustalenie liczby cząstek N i rozmiarów sieci (pojemnika) L

• Zadanie położeń początkowych cząstek

• Zadanie prędkości początkowych (metoda initialize())

• Obliczenie sił i przyspieszeń. Siła działająca na atom i -ty jest

~

F

i

=

P

j =1,j 6=i

F

ij

, gdzie F

ij

dostajemy z potencjału L-J. Stąd,

przyspieszenie ~

a

i

(t) = d ~

v

i

(t)/dt = d

2

r

i

(t)/dt

2

= ~

F

i

/m (metoda

computeAccelerations()).

• Całkowanie 3N równań różniczkowych drugiego rzędu z użyciem algorytmu

Verleta, z zadanym krokiem czasowym (metoda velocityVerlet()).

(25)

Przygotowanie programu MD (java) II

• Ponieważ liczba cząstek N w układzie jest stała i stała jest objętość L

3

_{, a}

siły L-J są konserwatywne, więc całkowita energia jest również zachowana.

Dla układu w równowadze termicznej spełniona jest zasada ekwipartycji

energii:

3(N − 1) ×

1

2 k

B

T =

*

m

2

N

X

i

v

i2

+

.

Nawias h· · · i oznacza średnią termiczną po zespole. Wielkość 3(N − 1)

jest liczbą temperaturowych stopni swobody (wszystkie stopnie swobody

minus trzy translacyjne). Stąd obliczamy temperaturę T .

• Program główny main() wywołuje metodę inicjalizacji (initialize()),

następnie oblicza krok za krokiem kolejne stany układu (velocityVerlet())

i na koniec oblicza temperaturę T (metoda instantenousTemperature).

(26)

Schemat algorytmu MD I

// Basis:

// http://www.physics.buffalo.edu/gonsalves/ public class MD {

static int N = 64; // number of particles

double[][] r = new double[N][3]; // positions

double[][] v = new double[N][3]; // velocities

double[][] a = new double[N][3]; // accelerations

double L = 10; // linear size of cubical volume

double vMax = 0.1; // maximum initial velocity component void initialize() {

// initialize positions

int n = (int)(Math.ceil(Math.pow(N, 1.0/3)));

// ^ number of atoms in each direction

double a = L / n; // lattice spacing

int p = 0; // particles placed so far

for (int x = 0; x < n; x++) for (int y = 0; y < n; y++)

for (int z = 0; z < n; z++) { if (p < N) {

r[p][0] = (x + 0.5) * a; r[p][1] = (y + 0.5) * a;

(27)

Schemat algorytmu MD II

r[p][2] = (z + 0.5) * a; } ++p; } // initialize velocities

for (int part = 0; part < N; part++) for (int i = 0; i < 3; i++)

v[part][i] = vMax * (2 * Math.random()); }

void computeAccelerations() {

for (int i = 0; i < N; i++) // set all accelerations to zero for (int k = 0; k < 3; k++)

a[i][k] = 0;

for (int i = 0; i < N-1; i++) // loop over all distinct pairs i,j for (int j = i+1; j < N; j++) {

double[] rij = new double[3]; // position of i relative to j double rSqd = 0;

for (int k = 0; k < 3; k++) { rij[k] = r[i][k] - r[j][k];

(28)

Schemat algorytmu MD III

a[j][k] -= rij[k] * f; } } } void velocityVerlet(double dt) { computeAccelerations(); for (int i = 0; i < N; i++)

for (int k = 0; k < 3; k++) {

r[i][k] += v[i][k] * dt + 0.5 * a[i][k] * dt * dt; v[i][k] += 0.5 * a[i][k] * dt;

}

computeAccelerations(); for (int i = 0; i < N; i++)

for (int k = 0; k < 3; k++) v[i][k] += 0.5 * a[i][k] * dt; }

double instantaneousTemperature() { double sum = 0;

for (int i = 0; i < N; i++) for (int k = 0; k < 3; k++)

(29)

Schemat algorytmu MD IV

return sum / (3 * (N - 1)); }

public static void main(String[] argv) { MD md = new MD();

md.initialize(); double dt = 0.01;

for (int i = 0; i < 2000; i++) { md.velocityVerlet(dt);

System.out.println("" + md.instantaneousTemperature()); }

}

(30)

Potencjał L-J

��

(31)

Dalsze ulepszenia algorytmu MD

• Periodyczne warunki brzegowe

-możliwość ucieczki cząstek z układu.

• Warunki początkowe dla położeń i

prędkości powinny być bardziej

fizyczne (np. sieć fcc - face centered

cubic lattice; prędkości zgodne z

rozkładem Maxwella-Boltzmanna).

• Układ powinien dążyć do stanu

równowagi przy zadanej

temperaturze.

• Pomiar (obliczenia) różnych

Niebieskie atomy (4) tworzą bazę komórki elementarnej. Ich położenia w jednostkach a są: (0, 0, 0), (0.5, 0.5, 0), (0.5, 0, 0.5), (0, 0.5, 0.5). W programie przesuwamy bazę o (0.5, 0.5, 0.5) tak, że atomy znajdują się

(32)

Vdistribution.java I

/**

The method gasdev() from Numerical Recipes returns random numbers with a Gaussian probability distribution

P(x) = (1/(2\pi\sigma^2)) * e^{-(x-x0)^2/(2*sig^2)},

with center x0 = 0 and unit variance \sigma^2 = 1. This function uses the Box-M\"uller algorithm.

*/

import java.util.Random; class Vdistribution {

static boolean available = false; static double gset;

Random ran = new Random(); double gasdev () { double fac, rsq, v1, v2; if (!available) { do { v1 = 2.0 * ran.nextDouble(); v2 = 2.0 * ran.nextDouble(); rsq = v1 * v1 + v2 * v2;

(33)

Vdistribution.java II

} while (rsq >= 1.0 || rsq == 0.0); fac = Math.sqrt(-2.0 * Math.log(rsq) / rsq); gset = v1 * fac; available = true; return v2*fac; } else { available = false; return gset; } }

void initVelocities(double[][] v, int N) { // Gaussian with unit variance

for (int n = 0; n < N; n++) // number of particles

for (int i = 0; i < 3; i++)

v[n][i] = gasdev(); // x, y, z components of velocity

}

// Test method; plot histogram public static void main(String[] argv){

(34)

Vdistribution.java III

for(int i=0;i<N;i++) { double g = vd.gasdev(); int j = (int)(g / hbin); if (j < nbins) ++bins[j];

//System.out.println(" "+available); }

int mmax = bins[0]; for(int i=1;i<nbins;i++)

mmax = bins[i]>mmax ? bins[i] : mmax; // plot hist

for(int i=0;i<nbins;i++){

double hist = (50*bins[i])/mmax; if(hist>0){ System.out.print("|"); for(int j=1;j<hist;j++) System.out.print("-"); } if(50*bins[i]/mmax>0) System.out.println("*"); } } }

(35)

Liczby gaussowskie

!" " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " # !" " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " # !" " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " # !" " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " # !" " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " # !" " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " # !" " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " # !" " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " # !" " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " # !" " " " " " " " " " " " " " " " " " " " " " " " " " " " " " # !" " " " " " " " " " " " " " " " " " " " " " " " " " " # !" " " " " " " " " " " " " " " " " " " " " " " " # !" " " " " " " " " " " " " " " " " " " " " # !" " " " " " " " " " " " " " " " " " " # !" " " " " " " " " " " " " " " " # !" " " " " " " " " " " " " # !" " " " " " " " " " " # !" " " " " " " " " " # !" " " " " " " # !" " " " " " # !" " " " # !" " " # !" " " # !" " # !" # !# !# !#

(36)

Schemat algorytmu MD I

import java.util.Random; class MD2 {

// simulation parameters

int N = 64; // number of particles

double rho = 1.0; // density (number per unit volume)

double T = 1.0; // temperature

double L; // linear size of cubical volume

/* Functions

void initialize(); // allocates memory, calls following 2 functions

void initPositions(); // places particles on an fcc lattice

void initVelocities(); // initial Maxwell-Boltzmann velocity distribution void rescaleVelocities(); // adjust the instanteous temperature to T

double gasdev(); // Gaussian distributed random numbers

*/

Random gen = new Random(); // random number generator

double[][] r = new double[N][3]; // positions

double[][] v = new double[N][3]; // velocities

(37)

Schemat algorytmu MD II

void initialize() { initPositions(); initVelocities(); } /**

Calculations of lattice constant from number of particles and number density, and positions; fcc latice

*/

void initPositions() {

// compute side of cube from number of particles and number density L = Math.pow(N / rho, 1.0/3);

// find M large enough to fit N atoms on an fcc lattice int M = 1;

while (4 * M * M * M < N) ++M;

(38)

Schemat algorytmu MD III

int n = 0; // atoms placed so far

for (int x = 0; x < M; x++) for (int y = 0; y < M; y++)

for (int z = 0; z < M; z++) for (int k = 0; k < 4; k++) if (n < N) { r[n][0] = (x + xCell[k]) * a; r[n][1] = (y + yCell[k]) * a; r[n][2] = (z + zCell[k]) * a; ++n; } }

/**

The method gasdev() from Numerical Recipes returns random number with a Gaussian probability distribution

P(x) = (1/(2\pi\sigma^2)) * e^{-(x-x0)^2/(2*sig^2)}, with center x0 = 0 and unit variance \sigma^2 = 1. This function uses the Box-M\"uller algorithm.

(39)

Schemat algorytmu MD IV

double gasdev () { double fac, rsq, v1, v2; if (!available) { do { v1 = 2.0 * gen.nextDouble(); v2 = 2.0 * gen.nextDouble(); rsq = v1 * v1 + v2 * v2; } while (rsq >= 1.0 || rsq == 0.0); fac = Math.sqrt(-2.0 * Math.log(rsq) / rsq); gset = v1 * fac; available = true; return v2*fac; } else { available = false; return gset; } } /**

Initial velocities of all particles (Maxwell-Boltzmann dist) */

(40)

Schemat algorytmu MD V

v[n][i] = gasdev();

// Adjust velocities so center-of-mass velocity is zero double vCM[] = {0, 0, 0};

for (int n = 0; n < N; n++) for (int i = 0; i < 3; i++)

vCM[i] += v[n][i]; for (int i = 0; i < 3; i++)

vCM[i] /= N;

v[n][i] -= vCM[i];

// Rescale velocities to get the desired instantaneous temperature rescaleVelocities();

} /**

Rescales velocities of particles */

void rescaleVelocities() { double vSqdSum = 0; for (int n = 0; n < N; n++)

for (int i = 0; i < 3; i++) vSqdSum += v[n][i] * v[n][i];

double lambda = Math.sqrt( 3 * (N-1) * T / vSqdSum ); for (int n = 0; n < N; n++)

(41)

Schemat algorytmu MD VI

for (int i = 0; i < 3; i++)

v[n][i] *= lambda; }

/**

Calculates accelerations of particles */

void computeAccelerations() {

for (int i = 0; i < N; i++) // set all accelerations to zero

for (int k = 0; k < 3; k++) a[i][k] = 0;

for (int i = 0; i < N-1; i++) // loop over all distinct pairs i,j

for (int j = i+1; j < N; j++) {

double[] rij = new double[3]; // position of i relative to j double rSqd = 0;

for (int k = 0; k < 3; k++) { rij[k] = r[i][k] - r[j][k];

(42)

Schemat algorytmu MD VII

rij[k] += L; }

rSqd += rij[k] * rij[k]; }

double f = 24 * (2 * Math.pow(rSqd, -7) - Math.pow(rSqd, -4)); for (int k = 0; k < 3; k++) { a[i][k] += rij[k] * f; a[j][k] -= rij[k] * f; } } } /**

Calculates positions/velocities from the Verlet algorithm */

void velocityVerlet(double dt) { computeAccelerations(); for (int i = 0; i < N; i++)

for (int k = 0; k < 3; k++) {

r[i][k] += v[i][k] * dt + 0.5 * a[i][k] * dt * dt; // use periodic boundary conditions

if (r[i][k] < 0) r[i][k] += L; if (r[i][k] >= L)

(43)

Schemat algorytmu MD VIII

r[i][k] -= L;

v[i][k] += 0.5 * a[i][k] * dt; }

computeAccelerations(); for (int i = 0; i < N; i++)

for (int k = 0; k < 3; k++) v[i][k] += 0.5 * a[i][k] * dt; } /** Calculate temperature */ double instantaneousTemperature() { double sum = 0;

sum += v[i][k] * v[i][k]; return sum / (3 * (N - 1)); }

(44)

Schemat algorytmu MD IX

public static void main(String[] argv) { MD2 md = new MD2();

md.initialize(); double dt = 0.01;

for (int i = 0; i < 2000; i++) { md.velocityVerlet(dt); System.out.println("" + md.instantaneousTemperature()); if (i % 100 == 0) md.rescaleVelocities(); } } }

(45)

Wyniki MD2

��

(46)

(47)

MD wydajniej

• Obcięcie potencjału – dla potencjałów krótkozasięgowych i szybko

malejących z odległością (np. r /σ > 1) można dokonać obcięcia

• .

(48)

Schemat algorytmu MD I

import java.util.Random; public class MD3 {

// simulation parameters

int N = 864; // number of particles

double rho = 1.0; // density (number per unit volume)

double T = 1.0; // temperature

double L; // will be computed from N and rho

double[][] r, v, a; // positions, velocities, accelerations

// variables to implement Verlet’s neighbor list

double rCutOff = 2.5; // cut-off on Lennard-Jones potential and force

double rMax = 3.3; // maximum separation to include in pair list

int nPairs; // number of pairs currently in pair list

int[][] pairList; // the list of pair indices (i,j)

double[][] drPair; // vector separations of each pair (i,j)

double[] rSqdPair; // squared separation of each pair (i,j)

int updateInterval = 10; // number of time steps between

(49)

Schemat algorytmu MD II

Random gen = new Random(); void initialize() { r = new double[N][3]; v = new double[N][3]; a = new double[N][3]; initPositions(); initVelocities();

// allocate memory for neighbor list variables nPairs = N * (N - 1) / 2;

pairList = new int[nPairs][2]; // to store indices i and j

drPair = new double[nPairs][3]; // to store components x,y,z

rSqdPair = new double [nPairs]; }

double computeSeparation (int i, int j, double dr[]) { // find separation using closest image convention double rSqd = 0;

(50)

Schemat algorytmu MD III

dr[d] += L; rSqd += dr[d]*dr[d]; } return rSqd; } void updatePairList() { nPairs = 0;

double[] dr = new double[3];

for (int i = 0; i < N-1; i++) // all distinct pairs

for (int j = i+1; j < N; j++) { // of particles i,j

double rSqd=computeSeparation(i, j, dr); if (rSqd < rMax*rMax) { pairList[nPairs][0] = i; pairList[nPairs][1] = j; ++nPairs; } } } void updatePairSeparations() { double[] dr = new double[3]; for (int p = 0; p < nPairs; p++) {

(51)

Schemat algorytmu MD IV

int j = pairList[p][1]; double rSqd = 0; rSqd = computeSeparation(i, j, dr); for (int d = 0; d < 3; d++) drPair[p][d] = dr[d]; rSqdPair[p] = rSqd; } } void computeAccelerations() {

for (int i = 0; i < N; i++) // set all accelerations to zero

for (int k = 0; k < 3; k++) a[i][k] = 0;

for (int p = 0; p < nPairs; p++) { int i = pairList[p][0]; int j = pairList[p][1];

if (rSqdPair[p] < rCutOff*rCutOff) { double r2Inv = 1 / rSqdPair[p]; double r6Inv = r2Inv*r2Inv*r2Inv;

(52)

Schemat algorytmu MD V

}

} }

void velocityVerlet(double dt) {

// assume accelerations have been computed for (int i = 0; i < N; i++)

for (int k = 0; k < 3; k++) {

r[i][k] += v[i][k] * dt + 0.5 * a[i][k] * dt * dt; // use periodic boundary conditions

if (r[i][k] < 0) r[i][k] += L; if (r[i][k] >= L) r[i][k] -= L; v[i][k] += 0.5 * a[i][k] * dt; } updatePairSeparations(); computeAccelerations(); for (int i = 0; i < N; i++)

for (int k = 0; k < 3; k++) v[i][k] += 0.5 * a[i][k] * dt; }

(53)

Schemat algorytmu MD VI

void initPositions() {

// compute side of cube from number of particles and number density L = Math.pow(N / rho, 1.0/3);

// find M large enough to fit N atoms on an fcc lattice int M = 1;

while (4 * M * M * M < N) ++M;

double a = L / M; // lattice constant of conventional cell

// 4 atomic positions in fcc unit cell double[] xCell = {0.25, 0.75, 0.75, 0.25}; double[] yCell = {0.25, 0.75, 0.25, 0.75}; double[] zCell = {0.25, 0.25, 0.75, 0.75};

int n = 0; // atoms placed so far

for (int x = 0; x < M; x++) for (int y = 0; y < M; y++)

for (int z = 0; z < M; z++) for (int k = 0; k < 4; k++)

(54)

Schemat algorytmu MD VII

} }

double gasdev () { double fac, rsq, v1, v2; if (!available) { do { v1 = 2.0 * gen.nextDouble(); v2 = 2.0 * gen.nextDouble(); rsq = v1 * v1 + v2 * v2; } while (rsq >= 1.0 || rsq == 0.0); fac = Math.sqrt(-2.0 * Math.log(rsq) / rsq); gset = v1 * fac; available = true; return v2*fac; } else { available = false; return gset; } }

(55)

Schemat algorytmu MD VIII

void initVelocities() {

// Gaussian with unit variance for (int n = 0; n < N; n++)

for (int i = 0; i < 3; i++) v[n][i] = gasdev();

// Adjust velocities so center-of-mass velocity is zero double[] vCM = {0, 0, 0};

vCM[i] += v[n][i]; for (int i = 0; i < 3; i++)

vCM[i] /= N;

v[n][i] -= vCM[i];

// Rescale velocities to get the desired instantaneous temperature rescaleVelocities();

(56)

Schemat algorytmu MD IX

vSqdSum += v[n][i] * v[n][i];

double lambda = Math.sqrt( 3 * (N-1) * T / vSqdSum ); for (int n = 0; n < N; n++)

for (int i = 0; i < 3; i++) v[n][i] *= lambda; }

double instantaneousTemperature() { double sum = 0;

sum += v[i][k] * v[i][k]; return sum / (3 * (N - 1)); }

public static void main(String[] argv) { MD3 md = new MD3(); md.initialize(); md.updatePairList(); md.updatePairSeparations(); md.computeAccelerations(); double dt = 0.01;

for (int i = 0; i < 1000; i++) { md.velocityVerlet(dt);

(57)

Schemat algorytmu MD X

System.out.println("" + md.instantaneousTemperature()); if (i % 200 == 0) md.rescaleVelocities(); if (i % md.updateInterval == 0) { md.updatePairList(); md.updatePairSeparations(); } } } }

(58)

Wyniki MD2

��

(59)