분자 전산 모사를 이용한 물성치 예측

장  욱, 이태용

한국과학기술원 화학공학과

Physical Properties Prediction Using Molecular Dynamics Simulation

Wook Chang and Tai-yong Lee

Department of Chemical Engineering, KAIST

서론

분자의 거동과 물성치를 파악하는 연구는 주로 실험을 통하여 진행되어 왔다. 이를 통하여 수많은 물성치와 상관식들을 예측하여 왔으며 이를 설명할 수 있는  모델을 개발하여 왔다. 컴퓨터 계산 능력이 향상됨에 따라서 분자 단위로 계산을 할 수 있게 되었고 많은 시스템에 관한 모사를 통하여 물성치를 실험을 거치지 않고 예측할 수 있다. 본 연구에서는 실험치와 비교 분석을 하기 쉬운 물을 예로 들어 분자 전산 모사를 하였고, 정적 열용량과 자체 확산 계수, 전단 점도를 얻었다.

분자 전산 모사

일반적으로 분자의 정확한 모사를 하기 위해서는 양자 역학에서 소개되는 Schr□dinger 식을 풀어서 각 분자들의 에너지 준위를 얻어야 하지만 수소와 헬륨같은 가벼운 기체가 아니고 어느 정도 온도 (50K) 이상의 경우에는 에너지 수준을 연속적으로 볼 수 있다. [1] 따라서 분자 전산 문제들은 각 분자들 사이에 존재하는 상호 인력, 전기적 인력, 수소 결합등으로 인한 상호 인력 등으로 표현되는 고전 역학 문제로 구성될 수 있으며 수많은 분자간의 조합으로 구성된 Newton의 운동 법칙을 푸는 문제로 바뀌어 진다.

다원자들의 화학 결합으로 구성된 물질의 상호 작용은 각 분자 사이의 상호 작용 (intermolecular interaction)과 화학 결합 내에 존재하는 상호 작용(intramolecular interaction)으로 분류가 가능하며 화학 결합은 각 분자들 사이에 움직임이나 위치에 제한을 줌으로서 표현할 수 있다. [2] 분자 전산 모사는 이러한 상호 작용을 바탕으로 분자가 받는 전체 힘을 계산하고 문제를 풀고 분자의 위치와 속도를 얻는 과정을 거치면서 물성치 예측의 기본 자료를 얻게 된다.

지배 방정식

한 분자 i가 다른 분자 j로부터 받는 전체 힘을 Fi라 하면 Fi는 다음과 같이 표현할 수 있다.

[eqn: bold{F}`-_i = SUM from {j} SUM from {beta} SUM from {alpha} bold{f}`_{i alpha j beta}]

()

분자 j에 존재하는 원자 □로 인하여 i 분자에 있는 원자 □가 받는 힘을 fi□j□이라 하면 전체 힘은 그의 합으로 표현되고 이는 Newton의 운동 법칙으로 설명된다.

[eqn: M_i bold{R}`_i prime prime = bold{F}`_i]

()

이를 이용하여 분자들의 위치와 속도 등을 파악하기 위해서는 fi□j□의 크기를 알아야 한다. 이 값들은 각 상호 작용을 고려한 포텐셜 함수의 구배(gradient)를 구함으로서 가능하다. j 분자에 있는 원자 □와 i 분자에 있는 원자 □사이의 포텐셜 에너지를 □i□j□(ri□j□)라 하면 계의 전체 포텐셜 에너지와 분자내 원자들 간의 힘은 식 (3)으로 표현할 수 있다.

[eqn: &U &=& SUM from {i} SUM from {j>i} SUM from {alpha} SUM from {beta}` phi_{i alpha j beta} (bold{r}`_{-i alpha j beta}`)# &bold{f}`_{-i alpha j beta} &=& - bold{TRIANGLED} phi_{i alpha j beta} (bold{r}`_{-i alpha j beta}`)]

()

분자내 존재하는 회전 효과는 분자 i의 토크를 계산하여 고려할 수 있다. 원자 □로 인해 분자 i가 받는 토크는 분자 i의 질량 중심(Ri)과의 차와 힘을 이용하여 식 (4)와 같이 나타낼 수 있다.

[eqn: bold{N}`_i &= SUM from {alpha} ` left(bold{r}`_{i alpha} - bold{R}`_i right) times bold{f}`_{i alpha}# bold{R}`_i &= 1/M_i SUM from {alpha} ` m_{i alpha} bold{r}`_{i alpha}#]

()

Newton-Euler 식을 이용하여 intertia tensor (Ii)를 계산하면 다음과 같이 쓸 수 있다.[eqn: &bold{I}`_i ~ CDOT bold{omega}`_i prime - bold{omega}`_i ~ times bold{I}`_i ~ CDOT bold{omega}`_i = bold{N}`_i# &bold{I}`_i = SUM from {alpha} ` m_{i alpha} left( (r_{i alpha} - R_i `)^2 bold {1} - (bold{r}`_{i alpha} - bold{R}`_i `) (bold{r}`_{i alpha} - bold{R}`_i `) right)]

()

여기서 □i는 분자 i의 각속도를 의미한다.

포텐셜 함수

액체 상태의 물에 대하여 여러 가지 포텐셜 함수들이 제안되었는데 본연구에서는 기본적으로 산소원자 사이의 Lennard-Jones 포텐셜과 수소 원자에 양전하가 존재하고 HOH의 결합각의 이등분선 위에 음전하가 존재하는 TIP4P 모델을 이용하여 수행하였다.

대상 시스템과 모사 조건

본 연구에서는 256개의 물분자에 대하여 밀도는 0.997g/□, 온도는 298K과 318K, canonical ensemble하에서 분자 스케일의 전산 모사를 수행하였다. 온도를 일정하게 유지하기 위해서 우리가 다루는 cell 주위에 가상의 heat bath를 둠으로서 시스템내의 운동 에너지를 일정하게 유지시키는 Nos□-Hoover 방법 [3]을 이용하였다. 컴퓨터 계산 효율을 높이기 위한 방법으로 이 분야에서 널리 사용되는 cut-off 방법과 주기 경계 조건을 사용하였다. 본 연구에서는 열역학적 물성치인 정적 열용량과 물질의 이동 현상에 관련된 자체 확산 계수와 전단 점도(shear viscosity)를 계산하였다.

결과 및 토론

1. 정적 열 용량

Canonical ensemble에서 내부 에너지 변동의 평균은 다음과 같이 정의되며 이는 다음과 같이 정적 열용량의 식으로 유도할 수 있다. [4]

[eqn: <(E-<E>^2 )> = <E^2 > - <E>^2 = k_B T^2 {partial E} over {partial T}= k_B T^2 C_v]

()

여기서 <> 기호는 시간에 대한 평균을 의미하며 모사 시간에 따른 정적 열용량의 수렴 과정을 그림 1에 나타내었다. 시간에 지남에 따라 일정한 값으로 점차 수렴하는 경향을 보이며 수렴된 값 두 경우 거의 차이가 크게 나지 않았으며 문헌치와 대체로 잘 맞음을 표 1에서 확인할 수 있다.

2. 자체 확산 계수

이동 현상을 표현하는 물성치들은 Einstein 식 [5]으로 표현이 가능하며 자체 확산 계수는 각각 식 (7)과 같이 위치의 제곱 평균과 직접 연결됨을 알 수 있다.

[eqn: D = lim from { t -> inf } {LEFT < (bold{r}(t) - bold{r}(0))^2 RIGHT >} over {6t}]

()

식 (7)은 Green-Kubo 관계식 [6]을 이용하여 적분식으로 바꿀 수 있으며 본 연구에서는 적분식을 이용하여 계산을 하였다. 그림 2는 모사 시간에 대한 자체 확산 계수를 나타내었는데 100 ps를 넘어 가면서 거의 일정한 값을 유지하며 수렴된 값은 표1에 비교하였다. 계산된 확산 계수들은 보편적인 액체 물질들의 확산 계수와 비슷한 크기를 가지며, 온도가 높아짐에 따라 확산 계수가 높아진다는 사실을 확인할 수 있었다.

그림  정적 열용량의 변화

그림  자체 확산 계수의 변화

3. 전단 점도

전단 점도도 확산 계수와 마찬 가지로 Einstein 식과 Green-Kubo 관계식을 이용하여 표현을 할 수 있다. 그러나, 전단 응력은 항상 전체 시스템과 연관이 되기 때문에 분자 전산 모사에서 사용하는 적은 수의 분자에 대한 평균을 통계적으로 처리하는 방법으로는 정확히 구할 수 없다. 동일 상이고 균일한 유체에 대해서는 다음과 같이 각 방향에 대해서 계산을 함으로서 구할 수 있다.

[eqn: eta_xy &= 1 over {Vk_B T} lim from { t -> inf } {LEFT < (bold{A}(t) - bold{A}(0))^2 RIGHT >} over {2t} # bold{A}(t) &= m SUM from { {i}=0} to N dot{bold{x}}_i (t) bold{y}_i (t) # eta &= 1 over 6 LEFT ( eta_xy + eta_xz + eta_yx + eta_yz + eta_zx + eta_zy RIGHT)]

()

시간 변화에 따른 계산된 점도의 값들을 그림 3에 표현하였는데 앞의 경우와 유사하게 일정한 값으로 수렴을 하고 있으며 수렴된 값은 문헌치와 비슷한 값을 얻었다는 사실을 확인할 수 있었고 온도가 높아짐에 따라 작아지는 사실을 알 수 있었다.

결론

본 연구에서는 물을 대상으로 TIP4P 모델을 이용하여 포텐셜 함수를 구성하고 지배 방정식을 세워 256개의 물분자와 canonical ensemble의 조건하에서 분자 전산 모사를 하였고, 이를 이용하여 정적 열용량과 자체 확산 계수, 전단 점도를 Einstein 식과 Green-Kubo 관계식을 이용하여 값을 예측하였다. 수렴된 값과 실험을 통한 문헌치와도 대체로 잘 맞음을 확인하였고, 앞으로 분자 전산 모사를 통하여 다성분 물질의 물성치 예측을 하는 과정은 물성치 파악과 실험, 계산에 필요한 노력들을 대체할 수 있는 좋은 방법이 되리라 생각된다.

표  계산된 물성치들의 수렴값과 비교

물성치	298 K 모사 결과	298 K 실험값 [7]	318 K 모사 결과	318 K 실험값 [7]
Cv (J Kg-1 K-1)	4151.52	4180.25	4131.56	4176.10
D (m2 s-1)	1.72×10-9	1.71×10-9	2.36×10-9	2.40×10-9
□ (Kg m-1 s-1)	8.44×10-4	9.09×10-4	6.87×10-4	6.12×10-4

그림  시간에 따른 전단 점도의 값

감사의 글

이 논문은 2000년도 두뇌한국 21 사업에 의하여 지원되었기에 감사드립니다.

참고 문헌

[1]	Rowley R. L. : "Statistical Mechanics for Thermophysical Property Calculations", PTR Prentice-Hall, Engleweood Cliffs, NJ (1994).

[2]	Allen M. P. and D. J. Tildesley : "Computer Simulation of Liquids", Clarendon·Oxford, New York, NY (1987).

[3]	Nos□ S. : Mol. Phys., 52, 255 (1984).

[4] Plischke M. and B. Bergersen : "Equilibrium Statistical Mechamics", 2nd ed. World Scientific, River Edge, NJ (1994).

[5] Einsten A. : "On the Movement of Small Particles Suspended in a Stationary Liquid Demanded by the Molecular-Kinetic Theory of Heat", Investigations on the Theory of the Brownian Movement, Dover, New York (1956).

[6] Haile J. M. : "Molecular Dynamics Simulation - Elementary Methods", John Wiley & Sons, Inc., New York (1997).

[7] Reid R. C., J. M. Prausnitz and B. E. Poling : "The Properties of Gases and Liquids", 4th ed. McGraw-Hill, New York (1986).