깁스 앙상블 몬테-카를로 모사를 통한 상평형 연구

안 지 철, 김 순 호, 김 화 용

서울대학교 응용화학부

A study of Gibbs Ensemble Monte-Carlo Simulation for Phase Equilibria

Jichul An, Sunho Kim, Hwayong Kim

School of Chemical Engineering, Seoul National University

서론

Gibbs Ensemble Monte Carlo(GEMC)는 유체의 상거동을 파악하는 하나의 방법으로서 최근에 주목을 받고 있다. GEMC는 분자 모델에 기인한 모사를 통하여 상평형도를 직접적으로 바로 구할 수 있는 장점이 있다. 최근에 이를 통하여 여러 유체 모델에 대한 상평형도가 그려지고 있다. 궁극적으로는 실제 유체를 모사하기 위한 방법으로 Gibbs Ensemble Monte Carlo 모사를 통하여 순수한 Lennard-Jones 유체 monomer와 8-mer에 대해 각각 상 거동을 살펴보았고 이성분 monomer 혼합물에 대해서도 모사를 통해 상 거동을 살펴보았다. 8-mer를 모사할 때에는 steric overlap을 방지하기 위한 방법으로 GEMC에 Configurational-bias 방법을 결합시켜 사용하였다. 이것은 여러 방향으로 사슬을 자라게 하여 가장 최적의 위치에 사슬을 옮기는 방법으로 이 방법을 사용하면 어느 정도 긴 사슬 유체의 모사도 가능하게 만든다. 또한 Square-well 유체에 대한 모사를 실행하여 상 거동을 확인하고 well-width가 다른 Square-well 유체 혼합물의 상 거동을 모사하여 보았다.

이론

GEMC는 분자 모델에 기인한 상평형을 직접적으로 계산하기 위한 방법이다. 이는 열역학적으로 연결된 두 개의 모사 상자로 구성되어 있고 각 상자는 periodic boundary condition을 가진다. 기액 평형의 경우라면 각 상자는 액상과 기상을 나타내는 상이라고 생각할 수 있다. GEMC는 NVT, NPT, □VT모사 방법 세가지를 합친 모사 방법으로 한 번의 모사를 통하여 두 상의 평형치을 예측할 수 있다는 것이 장점이다. GEMC는 두 상자를 합친 전체로서는 NVT를 유지하고 있지만 각 상자는 전체 부피를 유지하는 범위에서 부피를 변화시키고 두 상자간에 입자를 교환할 수 있다. 물론 전체를 NPT로 유지할 수도 있는데 이는 주로 혼합물의 모사에 주로 이용된다.

GEMC의 기본적인 이동 알고리즘은 particle displacement, volume change, particle transfer로 구성되어 있는데 각각 상자내의 내부평형을, 두 상자간의 동일 압력을, 두 상자간의 동일 chemical potential을 유지하기 위함이다.

①  Particle displacement.

임의로 한입자를 선정하여 그것이 속한 상자내에서 U(o)를 계산하고 다음과 같이 임의의 방향으로 움직이게 한다.

[eqn: X -> X + □]	(2.1)

그리고 그 위치에서의 U(n)를 구하고 움직임에 의해 야기된 시스템의 에너지변화(□U)를 계산할 수 있다. 이때 다음과 같은 확률로 이동을 받아들일 수 있다.

[eqn: acc(o->n) = min(1,exp{-□[U(n) - U(o)]})]	(2.2)

□U<0일 때, 즉 시스템이 좀 더 낮은 에너지 상태로 가게 된다면 입자의 움직임은 허용되고, 그 입자는 새로운 위치로 옮겨가게 된다. 만일 □U>0이면, [eqn: exp(-□□E)]의 확률을 가지고 움직이게 된다. 즉 0과 1 사이의 난수 P를 취하여 , [eqn: P<exp(-□□U)]라면 입자는 새로운 위치로 가게 되고, 그렇지 않다면 이전 위치로 돌아가게 된다.

사슬분자의 경우 단순 이동뿐만이 아니라 회전 이동도 필요하다. 이 때는 임의의 분자를 선택하고 그 분자의 분절도 임의로 선택하여 그 분절을 중심으로 회전시킨다. 회전 전후의 에너지를 계산하여 위의 식(2.2)와 동일한 확률로 받아들인다.

② Volume change

[eqn: V_total = V_1 + V_2]인 조건을 유지하면서 다음과 같이 부피를 변화시킨다.

[eqn: ln ({V_1 over V_2}) ->ln ({V_1 over V_2}) + □]	(2.3)

부피를 변화시키기 전의 전체 에너지(U(o))를 구하고 부피를 변화시킨 후의 전체 에너지(U(n))를 구한다. 이때 다음과 같은 확률로 부피 변화를 받아들일 수 있다.

[eqn: acc(o->n) = min(1,({{V_1 (n)} over {V_1 (o)} })^{n_1 +1} ({{ V_2 (n) } over {V_2 (o)}})^{n_2 +1} exp(-□□U))]	(2.4)

③ Particle transfer

임의의 입자를 선정하여 그 입자가 속한 상자내에서 입자가 받는 에너지(U(o))를 계산하고 다른 상자의 임의의 위치에 그 입자를 옮겨 그 위치에서 그 입자가 받는 에너지(U(n))를 계산한다. 이때 다음과 같은 확률로 입자의 이동을 받아 들일수 있다.

[eqn: acc(o->n) = min({1, {n_1 V_2} over {(n_2 +1)V_1} exp(-□□U)})]	(2.5)

ㅇ Configurational-bias Monte Carlo

일반적 GEMC에서 입자들은 임의의 위치로 다른 상자로 삽입된다. 일반적으로 Lennard-jones유체에서 한 입자를 액상에 삽입하기 위한 시도에 있어 다른 입자와 겹쳐지지 않을 확률이 0.005정도이다. 유사한 조건에서 n개의 분절을 가진 분자를 삽입할 때에는 [eqn: 0.005^n]이 될 것이다. 그러므로 분자 전이가 상당히 적은 확률로 이루어지기가 힘들어 질 것이다. 이를 해결하기 위해 GEMC에 Configurational-bias Monte Carlo를 결합시켜 long chain에 대한 전이를 수월하게 할 수 있다.

n개의 분절 입자를 상자 II에서 상자 I로 전이한다고 할 때 첫 번째 분절을 상자 I의 임의의 위치에 놓고 다음의 무게인자를 구한다.

[eqn: w_I (1) = exp(-□U_I,1 )]	(2.6)

여기서 [eqn: U_I,1`]은 첨가된 분절과 상자 I에 있는 나머지 입자와의 에너지이다. 다음 2번째 분절의 삽입에 있어 m번의 다른 방향으로의 성장을 시도하고 그 중 다음의 확률을 가지는 k번째 성장시도를 선택한다.

[eqn: P_k ={ exp(-□U_I,2^k )} over {w_I (2)}]	(2.7)

여기서 [eqn: P_k]는 k번째 시도가 선택될 확률이고, [eqn: U_I,2^k]는 k번째 시도후 자라난 두 번째 분절과 나머지와의 에너지이다. 그리고 [eqn: w_I (2)]는 다음과 같이 주어진다.

[eqn: w_I (2) = SUM from { { i}=1} to m exp(-□U_I,2^i )]	(2.8)

n분절까지 같은 방법을 적용하여 상자 I에서의 전체 Rosenbluth 무게 인자를 구한다.

[eqn: W_I =PI_j=1 ^n w_I (j)]	(2.9)

상자 II에서의 n분절의 제거에도 같은 방법이 쓰이나 한가지 다른점은 m번 시도 중 실제 분자의 위치가 포함되어야 한다는 것이다. 그래서 상자 II에서의 전체 Rosenbluth 무게인자는 위식과 같이 나타나며 [eqn: W_II]로 나타낸다.

이때에 n분절 입자의 다른 상자로의 전이는 다음의 확률로 받아들여진다.

[eqn: acc(o->n) = min LEFT { 1, {N_II V_I W_I} over {(N_I + 1)V_II W_II}RIGHT }]	(2.10)

결과 및 토론

1. Gibbs Ensemble Monte Carlo Simulation를 통하여 순수 Lennard Jones monomer의 상거동을 확인하였고 임계온도와 임계 밀도가 각각 1.35, 0.35에 근접함을 확인할 수 있었다.

2. GEMC에서 사슬분자의 모사시 steric overlap를 줄이기 위한 방법으로 Cofigurational-bias방법을 채용하여 freely-jointed Lennard-jones 유체 8-mer 에 대한 상거동을 확인하여 보았다. 이로서 실질적인 포텐셜 모델이 주어진다면 수많은 사슬형 실제 유체의 모사 가능성을 확인할 수 있었다.

3. 혼합물에 있어 순물질보다 수렴이 어렵고 큰 fluctuation이 생겼지만 혼합물 상평형 모사에 있어 가장 효과적인 방법인 것을 알수 있었다.

참고문헌

[1] F. A. Escobedo and J. J. De Pablo, Simulation and prediction of vapour-liquid equilibria for chain molecules, Molecular Physics, 87, 347, (1996)

[2] Panagiotopoulos, A.Z. (1987) Direct determination of phase coexistence properties of fluids by Monte Carlo simulation in a new ensemble, Molec. Phys. 61, 813-826.

[3] Panagiotopoulos, A.Z. (1992) Direct determination of fluid phase equilibria by simulation in the Gibbs ensemble: A review, Molec. Simulation 9, 1-23.

[4] Allen, M. P., and Tildesley, D.J., 1987, Computer Simulation of Liquids (Oxford: Clarendon Press).

Fig.3.1 단일성분 Lennard-jones monomer와 8mer

Fig.3.2. 이성분 Lennard-jones monomer

([eqn: sigma_1 /sigma_2 = 1], [eqn: epsilon_11 = 1], [eqn: epsilon_12 = 0.75], [eqn: epsilon_22 = 1]인 Mixture)

Fig.3.3 Square-well 단일성분 monomer 와 mixture

([eqn: sigma_1 /sigma_2 = 1], [eqn: epsilon_11 = epsilon_12 = epsilon_22 = 1], [eqn: lambda_11 = 1.4], [eqn: lambda_22 = 1.5]인 Mixture)