Molecular Dynamics를 기반으로 한 분자시뮬레이션과

Polyelectrolyte Complex Fluid의 미세구조특성 연구

Theories and Applications of Chem. Eng., 2002, Vol. 8, No. 1

이 현 수, 전 명 석

한국과학기술연구원 복잡유체연구팀

Molecular Simulations Based on Molecular Dynamics

and Study on the Microstructural Properties of Polyelectrolyte Complex Fluids

Hyun-Su Lee, Myung-Suk Chun

Complex Fluids Team, Korea Institute of Sci. and Tech., Seoul 130-650

서론

Fig. 1과 같이 직경이 수십 마이크론 이하인 미세한 흐름채널에서 complex fluid가 나타내는 다양하고 복잡한 현상은, physicochemical micro-hydrodynamics나 micro- rheology 분야에서 학술적 가치가 높은 연구주제라 여겨진다.  실제응용적 측면에서도, 모세관 전기영동 및 모세관 크로마토그래피에 의한 초미세분리, 멤브레인 기공, 젤, 혈관을 통한 물질이동, lab-on-a-chip과 같은 microfluidic devices의 설계분야에서 중요한 비중을 차지한다.

구형 콜로이드 입자와 함께 전형적인 complex fluid로서, 하전된 monomer들이 사슬로 연결된 구조인 polyelectrolyte 서스펜션의 static 및 dynamic 특성은 전부터 관심의 대상이 되어 왔다[1].  화학공학적 현상을 미세한 수준에서 이해하고자 하는 최근의 추세에 따라 polyelectrolyte 시스템에 대한 분자시뮬레이션 기법의 적용이 확장되고 있다[2].  본 연구에서는, polyelectrolyte의  conformation 특성 예측을 목적으로 기존에 연구된 molecular dynamics와 Monte Carlo의 hybrid 기법에서 한층 효율적으로 개선된 hybrid 기법을 개발하였다.  미세구조특성과 함께 polyelectrolyte chain의 dynamic 거동을 파악하고자 확산계수를 다루었고, 산출된 이들 결과에 대해 토론해 보고자 한다.

Fig. 1. Schematic of complex fluids in microchannel.

화학공학의 이론과 응용 제8권 제1호 2002년

Molecular Dynamics(MD)와 Monte Carlo(MC) 기법의 Hybrid Scheme

1. MD (Langevin) equation of motion

주어진 polyelectrolyte의 i-th monomer에 대해 표현되는 (1)식은, Newton의 운동방정식에 solvent와의 상호작용에 의한 frictional force와 random force Wi(t)를 추가한 Langevin 식이다.

[eqn: m {d^2 r_i} over { d t^2} ~=~ - del V_i ~-~GAMMA {d r_i} over dt ~+~W_i (t)]

(1)

여기서, m은 monomer의 질량, □는 frictional coefficient,[eqn: W_i (t) = delta_ij delta(t-t') 6 k_B T GAMMA] 는 white noise로서 구성된다.  Heat bath MC와 결합된 (1)식의 시간에 대한 적분은 velocity verlet 알고리즘에 의한다[3].

사슬에 대한 Hookean 모델과 하전된 monomer들 사이의 Debye-H□ckel(즉, Yukawa) potential에 따른 Hamiltonian은 다음과 같다.

[eqn: H ~=~ sum from {i=1} to {N-1} {3 k_B T} over {2 b^2 } ( {vec r}_i ``- {vec r}_i+1 ) ^2 ~+~ sum from {i=2} to {N} sum from {j=1} to {i-1} k_B `T `lambda_B {exp(- kappa r_ij )} over r_ij]

(2)

여기서, Bjerrum length □B(= e2/4□□kBT)는 electrostatic energy와 Boltzmann thermal energy가 같은 경우에 정의되는 특성길이이고, kB는 Boltzmann constant, □는 inverse Debye length, b는 monomer간의 bond length이다.

2. Hybrid scheme

하전된 monomer들 사이에 Debye-H□ckel potential이 작용하는 경우는 short range와 long range 길이척도 모두에 대해서 평형점을 찾아주어야 한다.  이때의 MC 단계에서 가장 효율적인 off-lattice pivot 알고리즘은 사슬에서 pivot point와 회전각도를 임의로 설정한 후 새로운 위치를 결정하는 것을 근간으로 하고 있다.  새로운 위치의 결정은 configuration 변화 전후의 Coulomb 에너지를 비교하여 결정한다.  Polyelectrolyte는 Fig. 2와 같이 totally stretched state와 임의의 회전각도를 준 random walks state의 두 경우를 초기조건으로 하고, 최종 conformation 값이 일치하는 시점을 평형에 도달한 것으로 간주하였다[4,5].

화학공학의 이론과 응용 제8권 제1호 2002년

계산 알고리즘은 (1)식에 대해서 103 order 횟수의 시간 적분(MD 단계)을 수행한 후에 104 order의 pivot 단계(MC Metropolis)을 통하여 large scale의 사슬구조를 평형상태에 도달하게 한다.  이후 105 order의 MD와 rare MC pivot 단계를 거쳐 계는 평형상태에 도달함으로써 계산수행을 마치게 된다.  계산과정 결과를 나타낸 Fig. 3a에서, 초기의 pivot 단계를 통해 large scale의 변이가 일어나고 이후의 MD 단계를 통해 서서히 평형에 도달함을 볼 수 있다.  초기 조건에 따라서 평형에 도달하는 시간에 차이가 있는데, totally stretched state에서 시작된 경우가 random walks에 비해 약 2배 정도 일찍이 도달함을 알 수 있다.  Fig. 3b에는 MC pivot 단계를 선택하는 tolerance에 따른 평형상태 도달 거동을 나타냈다.

Polyelectrolyte 서스펜션에 대한 시뮬레이션 결과

1. Conformation properties

MD의 velocity verlet 알고리즘의 시간 적분 단계에서 off-lattice pivot 알고리즘을 사용하게 된다.  Rare pivot 방식에 따라서, random number를 발생시켜 일정한 tolerance인 10-1 order 범위 이내로 들어오는 경우에만 pivot 알고리즘을 사용하였다[3-5].  반면에 본 연구에서는, 이같은 기존의 hybrid 방식과는 달리 MC Metropolis에 의한 off-lattice pivot 단계만 약 5×104회 계산하고 이를 초기조건으로 하여 MD 계산 단계에서는 pivot의 영향을 최소화하였다.  이러한 새로운 기법의 정확성을 검증하기 위해 문헌 결과와 비교하였는데, relative bond length 변화인 Fig. 4와 같이 정확히 일치하였다.  Polyelectrolyte의 conformation 특성에 관련하여 다음의 인자들을 정의할 수 있다.

End-to-End distance:        [eqn: R_{rm End} ~=~ sqrt { ({vec r}_N - {vec r}_1 )^2 }]	(3)

Radius of Gyration:     [eqn: R_G ~=~ sqrt{ 1 over N sum from i=1 to N bigg< ({vec r}_i - {vec R}_CM )^2 bigg>}]	(4)

여기서, [eqn: {vec R} _CM]은 center-of-mass vector이다.  한편, coil-to-rodlike microstructure transition 인자인 relative stretching ratio □(=[eqn: R_End^2 / R_G ^2])이 정의되는데, ideal chain (random walk)에서는 6, rodlike limit인 totally stretched state에서는 12, good solvent에서의 유연한 사슬에서는 6.3을 갖는다.

Structure Factor: Structure factor는 모든 length scale에 대한 사슬의 미세구조적인 매우 의미있는 정보로서, 실험 결과치와의 비교가 가능하게 하고 이론 연구자들에게 중요한 기본 데이터를 제공한다[4-6].  Wave number q에 대한 spherically-averaged structure factor Ssp(q)가 정의되는데, 이는 pair correlation function의 Fourier 변환을 의미한다.

[eqn: S_sp (q) &``=~ left < 1 over N left | Sum from i<j to N exp(i ``vec q cdot (vec r_i - vec r_j )) right |^2 right >~=~ 1 over N int d ``vec r ``exp(i ``vec q cdot vec r `) g(vec r `)#&``=~ 1 ``+`` 1 over N {sin(q b) } over {q b}]

(5)

본 연구에서 계산한 결과와 문헌 결과를 Fig. 5에 나타낸 바, monomer 갯수 N은 q값이 2□/REnd보다 작을 때는 structure factor에 크게 영향을 주고 있음을 알 수 있다. Structure factor에 대해 여러 model이 제안되어 있는데, q값이 2□/REnd보다 작은 경우에는 Debye 식이, q값이 증가함에 따라 generalized Debye 식이 적합한 것으로 보고되어 있다[4].

화학공학의 이론과 응용 제8권 제1호 2002년

2. Diffusion coefficient

Stokes-Einstein 관계식에서 polyelectrolyte의 수행된 trajectory에 대한 transport 특성인 diffusion coefficient는 다음으로 산출 가능하다[7].

[eqn: D ~=~ lim from { t -> inf } {left< (R_CM (t)-R_CM (0))^2 right>} over {6 t}]

(6)

여기서, 주어진 polyelectrolyte configuration에 대한 mean square displacement 값은 시간과 선형관계에 있다.

감사문

본 연구는 특정기초연구사업(KOSEF, R01-2001-00411)의 지원으로 수행되었음을 밝힌다.  아울러, 독일연방연구협회(DFG)의 지원에 의한 Max-Planck Institute for Polymer Research, Theory Group의 방문기간 중, Dr. Christian Holm과 Professor Kurt Kremer와의 기본 알고리즘에 관한 토의에 감사드린다.

참고문헌

화학공학의 이론과 응용 제8권 제1호 2002년

[1] Chun, Myung-Suk,  Park, O Ok and Yang, Seung-Man: J. Colloid Interface Sci., 161, 247(1993).

[2] Sadus, R.J.: "Molecular Simulation of Fluids: Theory, Algorithms and Object -Orientation", Elsevier, New York(1999).

[3] Binder, K. and Ciccotti, G. (Eds.): "Monte Carlo and Molecular Dynamics of Condensed Matter Systems", Societa Italiana Di Fisica, Bologna(1996).

[4] Micka, U. and Kremer, K.: Phys. Rev. E, 54, 2653(1996).

[5] Micka, U. and Kremer, K.: Europhys. Lett., 38, 279(1997).

[6] Micka, U., Holm, C. and Kremer, K.: Langmuir, 15, 4033(1999).

[7] Valdez, M.A. and Manero, O.: J. Colloid Interface Sci., 190, 81(1997).

화학공학의 이론과 응용 제8권 제1호 2002년

화학공학의 이론과 응용 제8권 제1호 2002년

Fig. 2. Two limiting initial conditions and their equilibrium configuration.

Fig. 3. Temporal variations of characteristic ratio of end-to-end distance and radius of gyration (N = 128, b = 2.0, □ = 0.01).

Fig. 4. Relative bond length stretched as a function of [eqn: 1/kappa] for N = 128, b = 2.

Fig. 5. Structure factor comparison between Micka and Kremer[4] and our present study (N = 128, b = 2.0, □ = 0.01).