Journal Search Engine
Search Advanced Search Adode Reader(link)
Download PDF Export Citaion korean bibliography PMC previewer
ISSN : 1225-309X(Print)
ISSN : 2288-7172(Online)
Journal of the mineralogical society of korea Vol.29 No.4 pp.209-220
DOI : https://doi.org/10.9727/jmsk.2016.29.4.209

A Molecular Dynamics Simulation Study of Hydroxyls in Dioctahedral Phyllosilicates

Sangbo Son, Kideok D. Kwon*
(Department of Geology, Kangwon National University, Chuncheon 24341, Korea)

Associate Editor: Sung Keun Lee

Corresponding author: +82-33-250-8553, kkwon@kangwon.ac.kr
September 20, 2016 December 16, 2016 December 20, 2016

Abstract

Clay minerals are a major player to determine geochemical cycles of trace metals and carbon in the critical zone which covers the atmosphere down to groundwater aquifers. Molecular dynamics (MD) simulations can examine the Earth materials at an atomic level and, therefore, provide detailed fundamental-level insights related to physicochemical properties of clay minerals. In the current study, we have applied classical MD simulations with clayFF force field to dioctahedral clay minerals (i.e., gibbsite, kaolinite, and pyrophyllite) to analyze and compare structural parameters (lattice parameter, atomic pair distance) with experiments. We further calculated vibrational power spectra for the hydroxyls of the minerals by using the MD simulations results. The MD simulations predicted lattice parameters and interatomic distances respectively deviated less than 0.1~3.7% and 5% from the experimental results. The stretching vibrational wavenumber of the hydroxyl groups were calculated 200-300 cm-1 higher than experiment. However, the trends in the frequencies among different surface hydroxyl groups of each mineral was consistent with experimental results. The angle formed by the surface hydroxyl group with the (001) plane and hydrogen bond distances of the surface hydroxyls were consistent with experimental result trends. The inner hydroxyls, however, showed results somewhat deviated from reported data in the literature. These results indicate that molecular dynamics simulations with clayFF can be a useful method in elucidating the roles of surface hydroxyl groups in the adsorption of metal ions to clay minerals.


분자동역학 시뮬레이션을 이용한 이팔면체 점토광물 수산기 연구

손 상보, 권 기덕*
강원대학교 자연과학대학 지질학과

초록

점토광물은 대기부터 지하수에 이르는 크리티컬존(critical zone) 영역에서 금속과 탄소 순환을 결정짓는 역할을 한다. 계산광물학 연구방법 중에 하나인 분자동역학(molecular dynamics) 시뮬레이션 은 지구물질을 원자단위로 계산하기 때문에, 점토광물의 물리화학적 현상들에 대해 원자수준의 자세 한 정보를 제공할 수 있다. 이번 연구에서는 clayFF 힘 장(force field)을 사용한 분자동역학 시뮬레이 션을 이팔면체 점토광물인 깁사이트(gibbsite, Al(OH)3), 카올리나이트(kaolinite, Al2Si2O5(OH)4), 파이 로필라이트(pyrophyllite, Al2Si4O10(OH)2)에 적용하여 300 K, 1기압조건에서 각 광물이 가지는 격자상 수와 원자간 거리를 계산하고 실험결과와 비교하였다. 더불어 수산기의 방향성 및 수소결합의 양상 그리고 파워스펙트럼(power spectrum)을 추가적으로 계산하였다. 계산 결과, 격자상수는 기존의 실험 결과와 약 0.1~3.7% 미만의 오차율을 보였으며, 원자간 거리는 실험결과와 약 5% 미만의 차이를 가 졌다. 깁사이트나 카올리나이트의 팔면체층 표면에 존재하는 수산기가 가지는 신축진동파수(stretching vibrational wavenumber)는 실험값 보다 약 200-300 cm-1 높게 계산되었지만, 팔면체층 표면에 존재하 는 수산기들의 상대적 크기의 경향은 기존 실험 결과와 일치하였다. 팔면체층 표면의 수산기가 (001) 면과 이루는 각도도 기존 실험결과와 상당부분 일치한 반면에 내부 수산기의 경우는 다소 차이를 보 였다. ClayFF를 사용한 분자동역학 시뮬레이션 연구 방법은 이팔면체 점토광물 표면 내(층간) 금속이 온 흡착에 대한 수산기의 역할을 규명하는데 유용한 연구방법이 될 수 있음을 시사한다.


    Ministry of Science, ICT and Future Planning
    National Research Foundation
    NRF-2015R1A4A1041105Kangwon National University
    520150427

    서 론

    이팔면체 점토광물(dioctahedral clay minerals) 은 알루미늄 팔면체층과 규소 사면체층으로 이루 어진 층상구조를 가지는 규산염광물로써 대기부터 지하수에 이르는 크리티컬존(critical zone) 영역 (Brantley et al., 2006; National Research Council, 2001)에서 금속과 탄소의 원소순환에 중요한 역할 을 한다. 또한, 점토광물은 미세한 입자크기로 인 한 넓은 표면적 및 구조 내에서의 이온치환으로 인 해 팽윤성과 뛰어난 이온흡착/교환 성질을 가져 의 약품, 화장품, 석유탐사와 오염물질 처리 등 다양 한 분야에서 응용되고 있다(Sposito et al., 1999; Karaborni et al., 1996). 특히, 이온흡착과 교환은 점토광물의 대표적인 성질로 그 기작에 대한 연구 가 많이 진행되어왔다. 이팔면체 점토광물에는 구 조적 수산기(OH)가 존재하는데, 이들은 층간 (interlayer) 물의 함량 또는 금속이온의 위치에 따 라 방향성을 달리하는 등으로 이온흡착 과정에 관 여한다고 오랫동안 제안되어 왔다(Cygan et al., 1998; Hawkins and Egelstaff, 1980). 따라서, 수산 기와 층간 이온흡착 사이의 관계, 나아가 점토광 물 물성에 대한 수산기의 영향을 이해하기 위해서 는 점토광물 내 수산기의 위치와 방향성을 파악하 는 것이 중요하다. 그러나 점토광물은 주로 입자 크기가 작고 결정도도 매우 낮아서 실험으로 분석 하기 힘들 뿐만 아니라, 특히, 수산기의 위치는 일 반 실험실에서 쓰이는 X-선 실험으로 관찰할 수 없다.

    계산광물학(computation mineralogy)은 컴퓨터 를 이용하여 지구물질을 원자-분자 스케일로 연구 하는 방법으로, 물질에 대한 물리화학적 특성을 규 명하고 예측하는데 널리 사용되고 있다(Cygan et al., 2004; Yi and Lee, 2014). 급속도로 발전하는 컴퓨터의 하드웨어와 계산을 수행하는 소프트웨어 의 발달에 따라 계산광물학의 정확도 및 계산속도 의 능력이 빠르게 진보되고 있다. 계산광물학은 실 험만으로 분석하기 힘든 점토광물의 구조와 점토 광물 표면에서의 오염물질 흡착 현상을 미시적 관 점에서 이해하는데 매우 중요한 연구방법이다.

    계산광물학은 크게 양자역학 계산방법(quantum mechanics)과 분자역학 계산방법(molecular mechanics) 으로 나뉜다. 양자역학 계산방법은 다전자 (many electrons)체계로 이루어진 시스템의 전자 (electron)들에 대한 슈뢰딩거(Schrödinger)파동방 정식의 해를 구하여 물질의 성질을 파악하는 방법 으로, 분자역학 계산방법에 비해 정확도 및 신뢰도 가 더 높지만 계산시간이 오래 걸리기 때문에 상대 적으로 분자역학 계산방법에서 다루는 시스템의 크기보다 매우 작은 크기 밖에 다룰 수가 없다. 이 에 반해 분자역학 계산방법은 양자역학 계산방법 보다 신뢰도 및 정확도는 낮지만, 전자가 아닌 원 자 또는 이온을 다루기 때문에 양자역학 계산방법 에 비해 계산속도가 월등히 빠르며, 시스템의 크기 도 양자역학 계산방법보다 수십~수백 배가 큰 시 스템을 다룰 수 있어 실제 실험에서 다루는 크기에 좀 더 근접하다. 하지만, 분자역학 계산방법은 원 자들의 상호작용에 따라 발생하는 에너지(쿨롱에 너지, 반데르발스 에너지 등)를 나타내기 위해 힘 장(force field)이라는 파라미터 세트를 사용하므로, 계산 결과의 정확도는 힘 장의 질(quality)에 크게 좌우된다. 점토광물에 대한 분자역학 계산에서는 clayFF (Cygan et al., 2004)라는 힘 장이 널리 사 용된다.

    이번 연구에서 사용한 분자동역학(classical molecular dynamics)시뮬레이션은 분자역학 계산 방법을 기반으로 구해진 시스템의 에너지를 이용 하여 원자들에게 작용되는 힘을 계산한 후에 뉴턴 의 운동방정식에 기초해 원자들이 시간에 따라 변 화하는 속도와 위치를 계산하는 통계역학 방법이 다. 이번 연구에서는 clayFF를 사용한 분자동역학 계산결과(격자상수, 원자간 거리, 수산기의 상태)를 기존 실험결과와 비교하여 이팔면체 점토광물에 대한 실험결과를 재현할 수 있음을 확인하고, 이팔 면체 점토광물 내 수산기 연구에 대하여 분자동역 학 시뮬레이션이 적합한지 테스트를 실시하였다.

    연구 방법

    이팔면체 점토광물

    이번 연구에서는 구조적전하가 없는 이팔면체 점 토광물인 단일 팔면체층구조의 깁사이트(gibbsite, Al(OH)3), 사면체층과 팔면체층이 1 : 1 구조를 가 지는 카올리나이트(kaolinite, Al2Si2O5(OH)4), 사면 체층과 팔면체층이 2 : 1 구조를 가지는 파이로필라 이트(pyrophyllite, Al2Si4O10(OH)2)를 계산하였다.

    힘 장(Force field, FF)

    ClayFF (Cygan et al., 2004)를 사용하여 시스템 이 갖는 전체에너지는 다음과 같이 표현된다.

    E total = E coul + E VDW + E bondstretch + E ∠⃒ bend
    (1)

    식 (1)에서 EcoulEVDW는 원자 간의 비결합 (non-bonded) 에너지를 표현하는 항(intermolecular terms)이며, EbondstretchE∠ bend는 물 분자나 수산 기 (OH) 결합의 신축진동(stretching vibration)과 굽힘 진동(bending vibration)에 의한 결합 내 에너 지를 나타내는 항(intramolecular terms)이다.

    Ecoul은 원자 간의 쿨롱에너지를 의미하며, 다음 과 같이 표현된다(식 (2)).

    E coul = e 2 4 π ϵ 0 i j q i q j r ij
    (2)

    식 (2)에서 qiqj는 각각 원자 i와 원자 j의 전 하를 의미하며, rij는 원자 ij의 거리를 의미한 다. e 는 전자의 전하량(-1.602 × 10-19C)을 의미하 고, ε0은 진공상태의 유전율(8.85419 × 10-12F/m) 을 의미한다. Ecoul과 같은 경우 원자들 간의 거리 가 멀리 떨어져 있어도 수렴하지 않는 특징을 보여 이를 계산하기 위해 Ewald summation (Ewald, 1921)을 사용하였으며 계산의 정확도로는 소수점 다섯 번째 자리까지 계산하도록 설정하였다.

    EVDW는 원자들 간에 작용하는 반데르발스 에너 지를 의미하며 다음으로 표현된다.

    E VDW = i j D o , ij R o , ij r ij 12 2 R o , ij r ij 6
    (3)

    식 (3)에 표현된 함수는 Lennard-Jones (12-6) 함 수(Jones, 1924)이며, rij는 두 원자 ij의 거리를 의미한다. (Ro,ij/rij)12항은 두 원자 간의 척력을 의 미하며, (Ro,ij/rij)6항은 두 원자 간의 인력을 나타 낸다. Ro,ij는 두 원자 ij가 평형을 이루었을 때 의 거리를 의미하고, Do,ij는 원자 ij가 이루는 Lennard-Jones (12-6) 함수에서의 퍼텐셜의 깊이를 나타내는 항이다. Ro,ijDo,ij는 Lorentz-Berthelot combination rule (Halgren, 1992)을 이용하여 얻 어진다. EVDW 같은 경우에는 Ecoul보다 상대적으로 빠르게 수렴하는 특징을 가지고 있어 불필요한 계 산을 줄이고 계산 효율을 높이기 위해 절단 거리 (cut off distance)를 설정하여 계산하게 되며, 통상 적으로 시스템이 가지는 격자상수에서 가장 짧은 길이의 반(half)을 절단 거리로 설정하게 된다.

    Ebondstretch는 원자간 공유결합 에너지를 표현하 는 항으로 특히 수산기(O-H)의 신축(stretching) 운 동에 따른 에너지를 표현하는데 사용된다. 종류로 는 조화 퍼텐셜(harmonic potential) (식 (4)) 또는 모스 퍼텐셜(Morse potential)(식 (5))이 있다.

    E bondstretch ij harmonic = k 1 r ij r o 2
    (4)
    E bondstretch ij morse = D 0 1 e α r r 0 2
    (5)

    식 (4)의 k1은 힘의 상수(force constant)로써 훅 의 법칙의 용수철상수(spring constant)와 대응되 며, 단위는 kcalmol-1 Å-2이다. 이는 적외선 분광 법(infrared spectroscopy) 또는 라만 분광법 (Raman spectroscopy)에 의하여 구할 수 있다. rij 는 원자 ij의 거리를 의미한다. r0는 두 원자가 평형을 이룰 때의 거리를 의미하며, X-ray 회절법 과 같은 실험에 의하여 얻어진다. 식 (5)의 r은 원 자 ij사이의 거리를 나타내고, D0는 원자 ij 가 해리되기 위하여 필요한 에너지(dissociation energy)에서부터 원점까지의 퍼텐셜 깊이를 의미 한다. α k e / 2 D 0 로써 표현되며 ke는 퍼텐셜 우 물 바닥에서의 힘의 상수(force constant)를 나타낸 다. 두 원자 간의 신축 진동은 조화 진동자 (harmonic oscillator)의 성질을 갖지 않고 비조화 성 진동자(anharmonic oscillator)의 성질을 가지고 있다. 따라서 모스 퍼텐셜의 경우는 두 원자들 간 의 신축진동양상을 조화 퍼텐셜보다 실제에 가깝 게 표현하기 위하여 제시된 퍼텐셜이다.

    E∠ bend는 세 개의 원자가 이루는 각도에 따른 굽 힘 진동(bending vibration)에 의한 에너지를 의미 하며 다음 식 (6)과 같다.

    E ∠⃒ bend ijk = k 2 θ ijk θ 0 2
    (6)

    식 (6)에서 k2도 힘의 상수(force constant)로써 식 (4)에서 말한 용수철상수와 같은 의미를 지니지 만 단위는 kcalmol-1rad-2이다. θijk는 원자 ij 그리고 k가 이루는 각도를 의미하고, θ0는 세 원자 가 평형일 때 이루는 각도를 의미한다. 이 항은 주 로 점토광물의 측면 표면(lateral edge)이나 외부 기저표면(basal surface)에서 일어나는 수산기와 양 이온 흡착과 같은 연구에 주로 사용되며, 표면의 성질과 관련이 없는 벌크(bulk) 점토광물 시뮬레이 션에서는 통상 사용하지 않아 이번 연구에서는 이 항에 대한 계산을 배제하였다(Cygan et al., 2004). 이번 연구에서는 Cygan et al. (2004)의 수산기에 대한 조화 퍼텐셜을 포함한 clayFF 힘 장을 수정 없이 사용하였고, 모스 퍼텐셜은 Greathouse et al. (2009)의 파라미터를 사용하였다.

    분자동역학 계산 (Molecular dynamics simulation)

    분자동역학 계산방법은 시스템을 구성하는 원자 들에 대하여 상호작용에 따른 에너지 및 힘이 주어 졌을 때 이를 고전역학의 뉴턴 운동방정식을 이용 하여 수치적분법으로 해석함으로써 원자들의 시간 에 따른 거동의 변화를 연구하는 방법이다. 시스템 을 이루고 있는 각각의 원자들은 다음과 같은 고전 역학 운동방정식을 따르게 된다.

    dE dr = m d 2 r dt 2
    (7)

    식 (7)는 시스템을 이루고 있는 각각의 원자에 작용하는 힘을 나타내며, 우변의 m은 각 원자의 질량을 나타내고 d2r/dt2은 각 원자의 가속도를 나 타낸다. 이번 연구에서 원자들의 운동방정식을 해 석하기 위한 알고리즘으로는 Verlet (1967)이 개발 한 Velocity Verlet algorithm을 사용하였다. 분자 동역학계산을 수행할 때에 시스템에 주기적경계조 건(periodic boundary conditions)을 적용하여 계산 속도의 향상 및 시스템 표면에 존재하는 원자들의 상호작용에 의한 왜곡효과를 방지하였다(González, 2011).

    통계 앙상블(Ensemble)

    분자동역학계산을 통하여 얻은 시간에 따라 변 화하는 시스템의 미시적인 정보들을 기존 실험결 과인 거시적인 정보와 비교하기 위해서는 앙상블 이라는 통계열역학에서 사용하는 개념이 필요하다. 앙상블(ensemble)이란, 거시적으로는 동일한 열역 학적 상태를 보일 수 있는 시스템의 미시적인 상태 들에 대한 모음을 의미한다. 앙상블을 결정하는 요 소로는 원자의 수 N, 온도 T, 압력 P, 에너지 E, 부피 V 등이 존재하며, 주어진 시스템을 등온-등 압 앙상블(isothermal-isobaric ensemble, NPT 일 정), 정준 앙상블(canonical ensemble, NVT 일정) 등과 같이 나눌 수 있다(Frenkel and Smit, 2002). 위의 요소들은 실제로 실험을 행할 때에 조절하는 온도와 압력 등과 대응되며 분자동역학 계산도 이 를 위해 온도와 압력을 통제하는 알고리즘을 사용 하게 된다. 이번 연구에서는 온도를 유지하기 위해 Nose-Hoover thermostat (Nosé, 1991)을 이용하였 으며, 압력을 유지하기 위해서는 Parrinello- Rahman barostat (Parrinello and Rahman, 1981)을 사용하였다.

    이번 분자동역학 계산에서는, Saalfeld and Wedde (1974)의 깁사이트(Fig. 1a), Bish (1993)의 카올리나이트(Fig. 1b), Lee and Guggenheim (1981) 의 파이로필라이트(Fig. 1c) 결정구조를 기반으로 2 X 2 X 2의 supercell을 이용하였다. 먼저 시스템 의 원자들에게 초기 속도를 부여하기 위해 NVT (T = 300 K) 하에서 200 ps (깁사이트의 경우 250 ps) 동안 계산한 후에 NPT [300 K, 1기압(1.013 × 10-4 Gpa)] 하에서 100 ps 동안 추가적으로 시뮬레 이션을 실시하여 300 K, 1기압 조건에서의 이팔면 체 점토광물들이 갖는 평형 구조를 얻었다. 이후 수산기의 진동모드 및 방향성을 분석하기 위해 추 가적으로 NVE ensemble 하에서 계산하였다. 원자 들의 위치 및 속도를 계산하는 시간간격(time step) 으로는 1 fs (1 × 10-15 sec)를 사용하였으며, 시뮬 레이션 결과는 2 fs마다 저장하여 계산결과를 분석 하였다. 모든 시뮬레이션 계산을 위해 Materials Studio의 Forcite 모듈을 사용하였다(Accelrys, 2016).

    동경분포함수(Radial distribution function)

    동경분포함수는 하나의 원자를 기준으로 했을 때 일정거리의 반경(r) 내에 다른 원자의 종류(α, β)가 존재하는 확률[g(r)]을 나타낸다(Hansen and McDonald, 1990).

    x α x β ρ g r = 1 N < i = 1 N α j = 1 N β δ r r i + r j >
    (8)

    식 (8) 좌변에서의 x는 원자의 종류를 의미하는 αβ의 몰 분율을 의미하며, ρ는 시스템의 밀도 이다. 우변에서의 N 은 시스템 내에 존재하는 총 원자들의 수, rirj는 기준이 되는 원자 ij의 위치(coordinate)를 나타내며, <>(angular brackets) 시간 평균값을 의미한다. 이번 연구에서는 동경분 포함수를 이용하여 특정 원자들 간의 결합거리 및 수산기가 이루는 수소결합의 길이를 계산하였다.

    파워스펙트럼(Power spectrum)

    속도자기상관 함수(velocity autocorrelation function, VACF)는 시스템을 구성하고 있는 원자들에 관하여 시간에 따른 속도의 변화에 대한 상관관계 를 나타낸다.

    ψ t = < υ i t υ i t 0 + t >
    (9)

    식 (9)에서 υ는 원자 i의 속도를 의미하며, t는 시간간격을 의미한다. 속도자기상관 함수를 푸리에 변환(fourier transform)했을 때 원자가 가지는 진 동모드(vibrational mode)에 따른 파수(wavenumber, cm-1)의 함수를 얻을 수 있으며 이를 파워스펙트럼 (power spectrum)이라 한다. 분자동역학 계산결과 를 저장하는 간격(step)에 따라 파수의 범위가 결 정되는데(예 : 저장 간격이 짧아질수록 높은 파수 의 범위까지 분석 가능), 이번 연구에서는 NVE ensemble 하에서 2 fs (1 fs = 1 × 10-15 sec)마다 분자동역학 시뮬레이션 결과를 저장하여 8,000 cm-1 범위까지 파워스펙트럼을 계산하였다. 파워스 펙트럼의 해상도는 속도자기상관 함수를 분석하는 시간이 증가할수록 높아지는데, 테스트 결과 10 ps (1 ps = 1 × 10-12 sec) 이상이면 수산기(OH)의 파워 스펙트럼을 분석하기 충분하였고, 이번 연구에서는 20 ps 동안 수산기들의 속도자기상관 함수를 분석 하여 파워스펙트럼 결과에 대하여 약 2 cm-1의 해 상도를 얻을 수 있었다.

    연구결과 및 토의

    격자상수 및 원자간 거리

    분자동역학 시뮬레이션으로 얻은 깁사이트, 카올 리나이트, 파이로필라이트의 격자상수와 원자들 간 의 거리를 기존 실험결과와 비교하였다(Table 1). 최대 오차율은 단일 팔면체층 구조가 3.7%, 사면 체층을 포함한 구조(1 : 1 구조, 2 : 1 구조)가 1.7% 의 오차율을 보여 규소 사면체층을 포함한 구조(규 산염광물)가 실험결과와 더 가까웠다. 깁사이트의 경우, 모스 퍼텐셜을 사용하였을 때 300 K, 1기압 조건에서 평형화를 얻을 수 없어 Table 1에서는 조화 퍼텐셜의 결과만 나타내었다.

    조화 퍼텐셜과 모스 퍼텐셜은 수산기의 결합을 다루기 때문에, 계산된 격자상수는 퍼텐셜의 종류 에 크게 다르지 않았다. 다만, 표면 수산기가 있는 팔면체층이 적층되는 광물의 경우, 수산기의 결합 및 수소 결합의 영향으로 c축 상수의 값이 퍼텐셜 에 따라 상이한 값을 보여주었다. 예를 들어, 파이 로필라이트의 c축의 경우에는 퍼텐셜에 따라 계산 된 결과가 크게 다르지 않았다. 그러나, 카올리나 이트 c축 상수의 경우, 모스 퍼텐셜의 결과가 조화 퍼텐셜의 결과 보다 크게 계산되었다(Table 1). 이 는 조화 퍼텐셜을 사용했을 때 수산기의 길이가 모 스 퍼텐셜보다 길게 계산되어 수소결합이 짧아지 는 양상을 보이므로 상대적으로 팔면체층 표면의 수산기와 그 다음 사면체층과의 인력이 강하게 작 용된 결과로 판단된다.

    원자간 거리(Si-O와 Al-O)에 대해서는 실험결과 와의 오차율은 Si-O와 Al-O의 길이가 각각 2.4- 3.7%와 0.5-2.6%를 보였지만, Si-O의 길이는 실험 값 보다 작게 계산된 반면, Al-O의 길이는 실험 보 다 큰 값을 보였다. 특히, 계산된 광물들 중에서는 카올리나이트의 Al-O 길이가 실험값과 거의 근접 한 값으로 계산되었는데, 이는 다른 두 가지 광물 의 구조를 고려했을 때, 팔면체층 표면에 존재하는 수산기와 사면체층의 인력이 Al-O 거리에 영향을 준 것으로 판단된다.

    수산기의 결합길이는 조화 퍼텐셜을 사용했을 때 1.02-1.03 Å으로 계산되었으며, 모스 퍼텐셜을 사용했을 때는 0.97 Å로 계산되어 기존 실험결과 는 모스 퍼텐셜이 더 잘 반영하였다. 깁사이트의 경우, 조화 퍼텐셜을 사용한 결과가 실험결과와 약 25%의 차이가 나는데, 이는 실험에서 발표한 수산 기 길이(Saalfeld and Wedde, 1974)가 수소의 위 치를 직접 측정한 결과가 아니고 X-ray 회절분석 을 통해 얻은 산소와 알루미늄 위치로부터 간단한 에너지 프로그램을 이용해 수소원자의 위치를 유 추했기 때문으로 판단된다. 또한 상대적으로 분자 역학 계산보다 신뢰성이 더 높은 양자역학 계산에 서도 깁사이트의 수산기 길이를 0.98 Å (Gale et al., 2001)로 계산한 것으로 보아 깁사이트 수산기 에 대한 분자동역학 시뮬레이션 결과는 조화 퍼텐 셜이 가지는 수산기 길이에 대한 오차율을 잘 반영 하여 얻어진 것이라고 판단된다. 마지막으로 이번 연구에서 얻은 분자동역학 결과들은 기존의 분자 동역학 계산결과(Cygan et al., 2004; Wang et al., 2006)와도 일치한다.

    깁사이트 수산기 구조

    Saalfeld and Wedde (1974)는 깁사이트의 수산 기를 두 그룹으로 분류하였다. OH(1), OH(2), OH(4)는 팔면체층과 거의 평행하게 위치하고, 나 머지 OH(3), OH(5), OH(6)는 팔면체층과 수직적 으로 배열되어 있으며, 이들은 불규칙한 삼각기둥 모양의 형태를 가진다. 분자동역학 계산결과 역시 수산기들을 수직과 수평으로 방향을 가지는 두 그 룹으로 나눌 수 있었다(Fig. 1a). Frost et al. (1999)은 적외선(infrared) 분광분석에서 얻은 깁사 이트의 3400-3700 cm-1의 파수(wavenumber) 범위 를 수산기의 신축진동파수(stretching vibrational wavenumber)로 해석하였다. 위의 범위에서 3460 cm-1은 수직적으로 배열된 수산기, 3529 cm-1와 3622 cm-1는 수평으로 배열된 수산기라고 구분하 였다. 이번 연구의 속도자기상관 함수으로부터 계 산된 파워스펙트럼을 분석한 결과, 3750-3850 cm-1 의 범위 중, 수평적으로 배열된 수산기가 수직적으 로 배열된 수산기보다 더 높은 파수를 보였다(Fig. 2). 계산된 범위는 실험결과와 정량적으로는 큰 차 이는 있지만, Frost et al. (1999)이 발표한 경향과 는 일치한다. Ruan et al. (2001)은 라만(Raman) 분광분석으로부터 얻은 깁사이트의 3360-3617 cm-1 의 범위 중, 가장 낮은 3364 cm-1는 층간에서 수소 결합을 가지는 수직으로 배열된 수산기에 의한 것 으로 해석하였다. 이 또한 분자동역학 계산으로 얻 은 수산기의 파워스펙트럼 경향과 일치한다. 각 수 산기 그룹들이 이루는 수소결합의 길이를 동경분 포함수를 통하여 계산한 결과, 층과 수직적으로 배 열되어 있는 수산기들[OH(3), (5), (6)](약 1.5~1.9 Å)이 수평적으로 배열되어 있는 수산기들[OH(1), (2), (4)](> 1.6 Å) 보다 수소결합의 길이가 더 짧 은 양상을 보였다(Fig. 3).

    카올리나이트 수산기 구조

    실온에서 측정한 실험데이터는 문헌에서 찾을 수 없지만, Bish (1993)는 1.5 K에서 중성자회절분 석을 통해 카올리나이트의 수산기를 두 그룹으로 나누었다. 팔면체층 표면에 존재하는 수산기[OH(2), (3), (4)] (surface OH)는 (001)면과 수직에 가깝게 위치하여 그 다음 사면체층과 수소결합을 이루고, 팔면체층 내부에 있는 수산기[OH (1)] (inner OH) 는 층과 거의 평행하게 위치한다(Fig. 1b). Bish (1993) 실험에 의하면, 수산기가 (001)면과 이루는 각도는 OH(1) < OH(4) < OH(3) < OH(2) 순으로 증가하였다. 열에너지(thermal energy)를 고려하지 않은 0 K에서 계산한 밀도범함수(density functional theory, DFT) 구조최적화 결과에서는 수산기의 각 도가 OH(1) < OH(4) < OH(2) < OH(3) 순으로 증가한다고 보고하였다(White et al., 2009)(Table 2).

    300 K에서 실시한 분자동역학 계산결과, 수산기 가 (001)면과 이루는 평균 각도는 두 포텐셜 모두, OH(1) < OH(2) < OH(3) < OH(4) 순서로 증가하 는 경향을 보여주었다(Table 2). 낮은 온도(e.g., 1.5 K와 0 K)에서 실시한 기존의 실험 및 계산결 과와는 비록 상이한 부분이 있지만, inner OH인 OH(1)이 가장 낮은 각도를 가지고, surface OH인 OH(2), (3), (4)가 상대적으로 높은 각도를 가지는 것은 기존 결과의 경향과 일치한다. Fig. 4는 실온 (300 K)에서 카올리나이트의 팔면체층 표면 수산 기들이 넓은 분포의 각도를 가질 수 있다는 것을 보여준다. OH(2)는 저각도에서 고각도까지 넓은 분포를 띄고, OH(3)와 OH(4)는 비록 고각도의 범 위에만 분포되어 있지만 약 15~20°의 범위를 가진 다. 하지만, 팔면체층 내부 수산기는 두 개의 퍼텐 셜을 사용한 결과 중에서 조화 퍼텐셜만 0~60°까 지 넓은 범위를 보여 정확한 판단을 내리기에는 무 리가 있다고 판단된다. 이 분자동역학 결과에서 팔 면체층 내부의 수산기에 대해서는 두 개의 퍼텐셜 에 의한 일관성을 보이지 않았지만, 팔면체층 표면 의 수산기들과 내부의 수산기를 구분할 수 있는 분 자동역학 시뮬레이션의 정확도는 보여주었다. 더불 어 두 퍼텐셜 모두 OH(2), OH(3), OH(4) 순으로 긴 수소결합의 양상을 보이고(Fig. 6), 이는 Bish (1993)의 수소결합 길이에 대한 상대적 순서와도 일치한다.

    IR 분광분석 실험에서 카올리나이트 수산기의 신축진동파수는 약 3600-3700 cm-1 범위에서 뚜렷 한 네 개의 밴드로 관찰된다고 보고되었으며 (Johnston et al., 1990), 네 개의 밴드는 각각 OH(1) < OH(4) < OH(3) < OH(2) 순서로 진동파 수가 증가한다(Bish, 1993). 분자동역학 시뮬레이 션에서 수산기들의 파워스펙트럼을 계산한 결과, 조화 퍼텐셜과 모스 퍼텐셜은 각각 약 3794-3823 cm-1, 약 3650-3750 cm-1의 범위를 보였다(Fig. 5). 조화 퍼텐셜 결과에서는 OH(4) < OH(3) < OH(1) < OH(2) 순으로 진동파수의 크기가 증가하였다 (Fig. 5a). 모스 퍼텐셜을 사용한 결과 역시 수산기 들의 진동파수의 상대적 크기가 조화 퍼텐셜 결과 와 일치하였다. 다만, 조화 퍼텐셜과 달리 모스 퍼 텐셜 결과는 상대적으로 넓은 진동파수를 보여주 어 정확도는 다소 떨어졌다(Fig. 5b). Bish (1993) 는 OH(1)의 진동파수가 가장 낮았지만, 분자동역 학 시뮬레이션에서는 OH(1)의 진동파수가 가장 낮 은 진동모드가 아니었다. 하지만, OH(2), (3), (4) 에 대해서는 진동파수의 크기 경향이 실험결과의 경향과 일치한다.

    파이로필라이트 수산기 구조

    파이로필라이트에 존재하는 수산기(inner OH)의 위치에 대한 실험결과는 문헌에서 찾을 수 없으며, 몇몇 분자모델링 연구결과만 존재한다(Refson et al., 2003). Giese (1973)는 원자들 간의 정전기적 에너지(electrostatic energy)를 이용한 계산을 통해 수산기가 (001)면과 26°의 각도를 이룬다고 제시 하였으며, 이 각도는 사면체층의 규소와 수산기의 수소 사이에서 발생하는 척력의 결과로 해석하였 다(Giese, 1973; Lee and Guggenheim, 1981). DFT 구조최적화 계산에 의하면(0 K), 파이로필라 이트의 수산기는 (001)면과 약 24-27°를 이룬다 (Refson et al., 2003). 이번 clayFF를 사용한 분자 동역학 시뮬레이션에서는 수산기가 평균 약 10° (조화 퍼텐셜)와 11° (모스 퍼텐셜)의 각도를 보여 주었다. 카올리나이트의 분자동역학 시뮬레이션에 서는 팔면체층 내부 수산기의 각도가 기존 연구결 과 보다 높은 값을 보여준 것과 대비된다. 아직 파 이로필라이트에 대한 수산기 방향을 실험한 결과 가 없기 때문에, 이번 분자동역학 시뮬레이션에서 계산한 수산기의 방향성에 대한 정확성을 명확히 판단하기 어렵다. 하지만, 상대적으로 신뢰성이 높 은 기존 DFT 결과와 비교했을 때, clayFF가 계산 한 수산기의 방향성은 약 10-20° 정도의 차이를 가 진다고 볼 수 있다.

    Fig. 7은 분자동역학 시뮬레이션결과로부터 계산 한 파이로필라이트의 수산기 신축진동파수에 해당 하는 파워스펙트럼이다. 조화 퍼텐셜로 계산했을 때는 약 3800-3830 cm-1의 범위를 보이며, 모스 퍼 텐셜로 계산했을 때는 약 3700-3760 cm-1의 범위 를 보였다. 이는 기존에 Farmer and Russell (1964)이 적외선 분광분석으로 얻은 수산기의 신 축진동파수에 해당하는 3675 cm-1와 비교하면 상 당히 큰 값이지만, 모스 퍼텐셜을 사용한 계산결과 가 실험결과와 더 근접하였다.

    결 론

    이번 연구에서는 점토광물 분자동역학 시뮬레이 션에 널리 사용되는 clayFF 힘 장(force field)을 이용하여 이팔면체 점토광물인 깁사이트, 카올리나 이트, 파이로필라이트가 300 K, 1기압(1.013 × 10-4 Gpa) 조건에서 가지는 격자상수와 원자들 간 의 결합거리를 계산하였다. 특히, 각 광물들에 존재 하는 수산기들의 경우에는 조화 퍼텐셜과 모스 퍼 텐셜을 각각 적용하여 계산하였다. 분자동역학 결 과를 기존 실험결과와 비교함으로써 점토광물 구조 계산에 대한 clayFF의 특성을 파악하고 그 정확성 을 판단하였다.

    이팔면체 점토광물의 격자상수에 대해서는 실험 결과와 약 0.1%-3.7%의 오차율을 보였다. 계산된 각각의 점토광물 구조에 대한 최대 오차율을 고려 했을 때는 깁사이트보다 사면체층을 포함한 카올 리나이트와 파이로필라이트에 대한 오차가 더 낮 았다. 격자상수 중에서 c축의 길이는, 팔면체층 표 면에 존재하는 수산기의 유무와 각기 다른 퍼텐셜 을 사용하여 얻은 수산기의 길이에 크게 영향을 받 는 것으로 나타났다. Si-O와 Al-O의 길이는 실험 결과와 약 0.5-3.7%의 오차율을 보였다. Si-O의 길 이는 실험값 보다 짧게 계산되었고, Al-O의 거리 는 실험값보다 길게 계산되는 경향을 보였다. 수산 기의 길이는, 조화 퍼텐셜 보다 모스 퍼텐셜을 사 용한 결과가 실험 및 DFT 결과와 더 근접한 결과 를 보여주었다.

    파워스펙트럼으로 계산한 수산기의 신축진동파 수는 두 퍼텐셜을 사용한 결과가 모두 실험값 보다 큰 값을 보였다. 비록 모스 퍼텐셜을 사용한 결과 가 더 실험값에 근접하였다 할지라도 비교적 넓은 범위를 가져 결과에 대한 정확도를 판단하기 힘들 며, 상대적으로 신뢰도가 더 높은 기존의 양자역학 계산(Bougeard et al., 2000; Gale et al., 2001; Greathouse et al., 2009)에서는 정확하게 실험값을 재현한 것으로 보아, 추후의 수산기의 신축진동파 수에 대한 더욱 정확한 분자역학 계산결과를 얻기 위해서는 수산기에 대한 파라미터의 개선이 필요 하다고 판단된다. 하지만, 두 퍼텐셜을 사용한 결 과에서 팔면체층 내부의 수산기(inner OH)를 제외 한 팔면체층 표면 수산기들(surface OH)의 신축진 동파수에 대한 상대적인 크기 순서에 해당하는 경 향은 실험결과와 일치하였으며, 팔면체층 표면의 수산기가 가지는 수소결합 길이를 계산한 결과도 실험에서 제시한 경향과 일치하였다.

    카올리나이트의 수산기들이 (001)면과 이루는 각도에 대한 시뮬레이션 결과는 실험에서 제시한 상대적으로 낮은 각도를 이루는 팔면체층 내의 수 산기(inner OH)와 큰 각도를 이루는 팔면체층 표 면의 수산기(surface OH)들을 명확히 구분할 수 있다는 것을 보여주었다. 300 K에서 계산한 수산 기의 평균 각도는 낮은 온도에서의 실험(1.5 K) 값 과 정량적인 차이를 가지지만, 이는 다시 말해 카 올리나이트의 수산기가 실온(300 K)에서 상당히 넓은 범위를 가질 수 있다는 가능성을 보여준다. 카올리나이트의 팔면체층 표면에 존재하는 수산기 들은 (001)면과 각도를 이루는 양상이 실험 및 양 자역학 계산결과의 경향과 비슷한 반면에, 파이로 필라이트의 팔면체층 내부에 존재하는 수산기들은 (001)면과 이루는 각도가 양자역학 계산과는 상이 한 경향을 보여주었다. 따라서, 위의 파워스펙트럼 및 수소결합의 길이, (001)면과의 각도의 결과는 clayFF를 사용한 분자동역학 계산이 팔면체층 내 부의 수산기보다 팔면체층 표면에 위치한 수산기 를 실험결과에 더 가깝게 재현함을 알 수 있으며, 이는 추후에 분자동역학 시뮬레이션을 통한 점토 광물의 수산기 연구를 할 때 신중히 고려해야 할 시뮬레이션의 정확도 및 변수로 판단된다.

    모스 퍼텐셜을 사용한 결과가 조화 퍼텐셜을 사 용한 결과보다 수산기의 구조적인 특징을 정량적 으로 더 잘 반영하였지만, 수산기들이 가지는 상대 적인 경향에 대해서는 두 퍼텐셜을 사용한 결과가 큰 차이를 보이지 않았다. 따라서, 깁사이트의 경우 와 같이 평형상태의 구조를 얻을 수가 없는 모스 퍼텐셜의 호환성 문제를 고려한다면, 점토광물 내 의 수산기와 관련된 현상을 관찰하기에는 상대적으 로 평형상태를 얻기 쉽고 clayFF와 호환성이 높은 조화 퍼텐셜을 사용하는 것이 계산 효율적으로 더 낫다고 판단된다. 이번 연구결과는 clayFF를 사용 하는 분자동역학 계산방법이 실험으로 확인하기 힘 든 점토광물 구조 및 표면 수산기의 구조를 연구할 수 있다는 것을 보여준다. 분자동역학 계산은 구조 전하를 가지는 점토광물의 수산기와 층간 금속이온 사이의 관계 규명에도 유용한 방법임을 시사한다.

    사 사

    논문은 2016년도 정부(미래창조과학부)의 재원으로 한 국연구재단의 지원(NRF-2015R1A4A1041105)과 2015년도 강원대학교 대학회계 학술연구조성비(관리번호-520150427) 지원을 받아 작성하였다. 분자동역학 시뮬레이션에 도움을 준 Aric Newton 박사께 감사드린다. 원고 향상에 도움을 주신 이용재 교수님과 익명의 두 분의 심사자께도 감사 드린다.

    Figure

    JMSK-29-4-209_F1.gif

    Snapshots of MD simulations (at 300 K, 1 atm) for (a) gibbsite, (b) kaolinite, (c) pyrophyllite. Polyhedra represent Al-octahedra (pink) or Sitetrahedra (yellow). White balls represent hydrogen.

    JMSK-29-4-209_F2.gif

    Vibrational power spectra of hydroxyls in gibbsite obtained from MD simulations performed with the harmonic potential.

    JMSK-29-4-209_F3.gif

    Distribution of H-bond distances in gibbsite obtained from MD simulations with the harmonic potential.

    JMSK-29-4-209_F4.gif

    Hydroxyl orientation angles relative to the (001) plane of kaolinite obtained from MD simulations performed with the (a) harmonic potential and (b) Morse potential.

    JMSK-29-4-209_F6.gif

    Hydrogen-bond distances of hydroxyls in kaolinite obtained from MD simulations performed with the (a) harmonic potential (b) Morse potential.

    JMSK-29-4-209_F5.gif

    Vibrational power spectra of hydroxyls in kaolinite derived from MD simulations performed with the (a) harmonic potential and (b) Morse potential.

    JMSK-29-4-209_F7.gif

    Vibrational power spectra of hydroxyls in pyrophyllite derived from MD simulations.

    Table

    Structural parameters of dioctahedral phyllosilicates obtained by molecular dynamics (MD) simulations

    *EXP = Experimetal results; N.A. = Not available; values in ( ) = Full width at half maximum of a corresponding RDF peak
    aSaalfeld and Wedde (1974)
    bBish (1993)
    cLee and Gugenheim (1981)
    dDensity functional theory calculation result from Gale (2001)

    Averaged angles (°) of hydroxyls of kaolinite with respect to its basal plane

    aExperiment at 1.5 K(Bish, 1993)
    bCalculation at 0 K(White, 2009)
    cCurrent simulation at 300 K

    Reference

    1. Accelrys, Inc (2016) Forcite module , Materials Studio. San Diego,
    2. Bish DL (1993) Rietveld refinement of the kaolinite structure at 1.5 K , Clays and Clay Minerals, Vol.41 ; pp.738-744
    3. Bougeard D , Smirnov KS , Geidel E (2000) Vibrational spectra and structure of kaolinite : A computer simulation study , The Journal of Physical Chemistry B, Vol.104 ; pp.9210-9217
    4. Brantley SL , White TS , White AF , Sparks D , Richter D , Pregitzer K , Derry L , Chorover J , Chadwick O , April R , Anderson S , Amundson R (2006) Frontiers in Exploration of the Critical Zone , Report of a workshop sponsored by the National Science Foundation (NSF). An NSF-Sponsored Workshop. Newark. DE, Vol.30 ; pp.24-26
    5. Cygan RT , Nagy KL , Brady PV (1998) Molecular models of cesium sorption on kaolinite , Adsorption of Metals by Geomedia. Academic, ; pp.383-399
    6. Cygan RT , Liang J , Andrey GK (2004) Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field , Journal of Physcal Chemistry B, Vol.108 ; pp.1255-1266
    7. Ewald PP (1921) The computation of optical and electrostatic lattice potentials , Annals of Physics, Vol.64 ; pp.253(in German without English abstract)
    8. Farmer VT , Russell JD (1964) The infra-red spectra of layer silicates , Spectrochimica Acta, Vol.20 ; pp.1149-1173
    9. Frenkel D , Smit B (2002) Understanding Molecular Simulation , From Algorithms to Applications, Academic,
    10. Frost RL , Kloprogge JT , Russell SC , Szetu JL (1999) Vibrational spectroscopy and dehydroxylation of aluminum (oxo) hydroxides gibbsite , Applied Spectroscopy, Vol.53 ; pp.423-434
    11. Gale JD , Rohl AL , Milman V , Warren MC (2001) An ab initio study of the structure and properties of aluminum hydroxide: gibbsite and bayerite , The Journal of Physical Chemistry B, Vol.105 ; pp.10236-10242
    12. Giese RF (1973) Hydroxyl orientation in pyrophyllite , Nature, Vol.241 ; pp.151
    13. González MA (2011) Force fields and molecular dynamics simulations , Collection SFN, Vol.12 ; pp.169-200
    14. Greathouse JA , Durkin JS , Larentzos JP , Cygan RT (2009) Implementation of a Morse potential to model hydroxyl behavior in phyllosilicates , Journal of Chemical Physics, Vol.130 ; pp.134713
    15. Halgren TA (1992) The representation of van der Waals (vdW) interactions In molecular mechanics force fields: potential form. combination rules. And vdW parameters , Journal of the American Chemical Society, Vol.114 ; pp.7827-7843
    16. Hansen JP , McDonald IR (1990) Theory of Simple Liquids: With Applications to Soft Matter , Academic,
    17. Hawkins RK , Egelstaff PA (1980) Interfacial water structure in montmorillonite from neutron diffraction experiments , Clays and Clay Minerals, Vol.28 ; pp.19-28
    18. Johnston CT , Agnew SF , Bish DL (1990) Polarised single-crystal fourier-transform infra-red microscopy of Ouray dickite and Keokuk kaolinite , Clays and Clay Minerals, Vol.38 ; pp.573-583
    19. Jones JE (1924) On the determination of molecular fields. II. From the equatio n of state of a gas , The Royal Society of London A. Mathematical. Physical and Engineering Sciences, Vol.106 ; pp.463-477
    20. Karaborni S , Smit B , Heidug W , Urai J , VanOort E (1996) The swelling of clays : molecular simulations of the hydration of montmorillonite , Science, Vol.271 ; pp.1102-1104
    21. Lee JH , Guggenheim S (1981) Single crystal X-ray refinement of pyrophyllite-1tc , American Mineralogist, Vol.66 ; pp.350-357
    22. National Research Council (2001) Basic Research Opportunities in Earth Science, The National Academies Press, ; pp.154
    23. Nosé S (1991) Constant temperature molecular dynamics methods , Progress of Theoretical Physics Supplement, Vol.103 ; pp.1-46
    24. Parrinello M , Rahman A (1981) Polymorphic transitions in single crystals crystals: A new molecular dynamics method , Journal of Applied Physics, Vol.52 ; pp.7182-7190
    25. Refson K , Park SH , Sposito G (2003) Ab initio computational crystallography of 2:1 clay minerals. 1. Pyrophyllite-1Tc , The Journal of Physical Chemistry B, Vol.107 ; pp.13376-13383
    26. Ruan HD , Frost RL , Kloprogge JT (2001) Comparison of Raman spectra in characterizing gibbsite Bayerite. diaspore and boehmite , Journal of Raman Spectroscopy, Vol.32 ; pp.745-750
    27. Saalfeld H , Wedde M (1974) Refinement of the crystal structure of gibbsite , Al(OH)3 Zeitschrift for Kristallographie-Crystalline Materials, Vol.139 ; pp.129-135
    28. Sposito G , Skipper NT , Sutton R , Park SH , Soper AK , Greathouse JA (1999) Surface geochemistry of the clay minerals , Proceedings of the National Academy of Sciences, Vol.96 ; pp.3358-3364
    29. Verlet L (1967) Computer “experiments” on classical fluids I. Thermodynamical properties of Lennard. Jones molecules , Physical review, Vol.159 ; pp.98
    30. Wang J , Kalinichev AG , Kirkpatrick RJ (2006) Effects of substrate structure and composition on the structure, Dynamics and energetics of water at mineral surfaces: A molecular dynamics modeling study , Geochimica et cosmochimica acta, Vol.70 ; pp.562-582
    31. White CE , Provis JL , Riley DP , Kearley GJ , van Deventer JS (2009) What is the structure of kaolinite? Reconciling theory and experiment , Journal of Physical Chemistry B, Vol.113 ; pp.6756-6765
    32. Yi YS , Lee SK (2014) Quantum Chemical Calculations of the Effect of Si-O Bond Length on X-ray Raman Scattering Features for MgSiO3 Perovskite , Journal of the Mineralogical Society of Korea, Vol.27 ; pp.1-15(in Korean with English abstract)