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.30 No.4 pp.161-172
DOI : https://doi.org/10.9727/jmsk.2017.30.4.161

A Molecular Dynamics Simulation Study of Trioctahedral Clay Minerals

Jiyeon Lee1,2, Jin-Yong Lee1,2, Kideok D. Kwon1,2*
1Critical zone Frontier Research Laboratory, Kangwon National University, Chuncheon 24341, Korea
2Department of Geology, Kangwon National University, Chuncheon 24341, Korea
Corresponding author : +82-33-250-8553, kkwon@kangwon.ac.kr
20171019 20171212 20171218

Abstract

Clay minerals play a major role in the geochemical cycles of metals in the Critical Zone, the Earth surface-layer ranging from the groundwater bottom to the tree tops. Atomistic scale research of the very fine particles can help understand the fundamental mechanisms of the important geochemical processes and possibly apply to development of hybrid nanomaterials. Molecular dynamics (MD) simulations can provide atomistic level insights into the crystal structures of clay minerals and the chemical reactivity. Classical MD simulations use a force field which is a parameter set of interatomic pair potentials. The ClayFF force field has been widely used in the MD simulations of dioctahedral clay minerals as the force field was developed mainly based on dioctahedral phyllosilicates. The ClayFF is often used also for trioctahedral mineral simulations, but disagreement exits in selection of the interatomic potential parameters, particularly for Mg atom-types of the octahedral sheet. In this study, MD simulations were performed for trioctahedral clay minerals such as brucite, lizardite, and talc, to test how the two different Mg atom types (i.e., ‘mgo’ or ‘mgh’) affect the simulation results. The structural parameters such as lattice parameters and interatomic distances were relatively insensitive to the choice of the parameter, but the vibrational power spectra of hydroxyls were more sensitive to the choice of the parameter particularly for lizardite.


삼팔면체 점토광물에 대한 분자동역학 시뮬레이션 연구

이 지연1,2, 이 진용1,2, 권 기덕1,2*
1강원대학교 크리티컬존 선도 연구실
2강원대학교 자연과학대학 지질학과

초록

점토광물은 지하수 바닥부터 산림에 이르는 지구의 얇은 표면에 해당하는 ‘크리티컬존(critical zone)’에 존재하는 금속원소의 지구화학적 순환에 중요한 역할을 한다. 입자 크기가 매우 작은 점토광 물에 대한 원자 수준(atomistic scale)의 연구는 지구화학적 순환 과정에 대한 정확한 기작(mechanism) 을 규명할 수 있을 뿐만 아니라 재료개발과 같은 산업분야에도 응용될 수 있다. 원자 간의 페어 퍼텐 셜(pair potential)을 파라미터화한 힘 장(force field)을 사용하는 분자동역학(molecular dynamics) 컴퓨 터 시뮬레이션은 원자 수준의 정보를 제공할 수 있기 때문에 실험과 함께 점토광물의 결정구조와 반 응도 연구에 사용된다. 점토광물 시뮬레이션을 위한 힘 장으로는 이팔면체(dioctahedral) 광물을 기반 으로 만들어진 ClayFF 힘 장이 보편적으로 사용된다. 삼팔면체(trioctahedral) 광물 시뮬레이션에도 ClayFF를 사용하는 연구가 보고되고 있으나, 같은 광물을 계산하더라도 각 연구마다 다른 파라미터 값을 사용하고 있기 때문에 파리미터 선택이 시뮬레이션의 정확도에 어떤 영향을 미치는지 체계적인 테스트가 필요하다. 이번 연구에서는 삼팔면체 광물인 수활석, 리자다이트, 활석을 대상으로 팔면체 마그네슘(Mg)의 원자간 페어 퍼텐셜을 나타내는 파라미터 ‘mgo’와 ‘mgh’를 각각 사용하여 분자동역 학 시뮬레이션 계산결과를 비교하였다. 격자상수, 원자 간의 거리 등 삼팔면체 점토광물의 결정구조는 주어진 두 가지 파라미터에 관계없이 거의 일정한 결과를 보여주었지만, 진동 파워 스펙트럼(vibrational power spectrum)으로 계산한 수산기의 진동수는 파라미터에 따라 상대적으로 뚜렷한 차이를 보였다.


    Ministry of Education, Science and Technology
    2014-3rd-10C-M001Kangwon National University
    520160469

    서 론

    암석의 풍화에 의해 생성된 점토광물은 토양이 나 퇴적물에 다양한 형태로 존재한다(Grim, 1968; Moore and Reynolds, 1997). 일반적으로 점토광물 은 2 μm 이하의 입자를 갖는 층상규산염광물 (phyllosilicate minerals)에 해당한다(Bailey, 1979). 콜로이드 크기의 입자 때문에 넓은 표면적을 갖으 며, 양이온 치환에 의한 구조전하(layer charge)를 띠는 성질 때문에 토양이나 지하수에 존재하는 물 분자, 중금속, 방사핵종이나 Ca2+, Mg2+, K+ 같은 양이온을 잘 포착한다(Tucker, 1999; Dube et al., 2001). 이러한 특징들 때문에 점토광물은 다양한 반응의 촉매나 오염물을 제거하는 흡착제로 주로 사용되고 고준위 방사성 폐기물 처리에 방습제 (barrier materials)로도 제안되고 있다(Soma and Soma, 1989; Savage, 1995; Costanzo, 2001).

    점토광물의 반응은 원자 수준(atomistic scale)으 로 일어나기 때문에, 점토광물에 대한 원자-분자 단위의 연구가 매우 중요하다(Cygan, 2001). 원자 단위의 구조를 측정하기 위해서는 싱크로트론 방 사광과 같은 고-에너지의 소스(source)를 제공하는 시설을 이용지만, 시설이용의 경제적 공간적 제약 이 존재한다. 그 외에도 작은 입자 크기와 낮은 결 정도를 갖는 점토광물은 원자 수준의 실험연구가 어렵다. 컴퓨터 분자 모델링(molecular modeling) 은 이러한 실험의 제약을 보완하는 방법으로 광물 의 결정구조, 물리적, 열역학적 성질에 대하여 분 자에서 원자 수준까지 연구할 수 있다(Cygan et al., 2004). 이번 연구는 분자 모델링의 한 방법인 전통 분자동역학(classical molecular dynamics, MD) 시뮬레이션을 삼팔면체(trioctahedral) 점토광 물의 구조연구에 적용하였다.

    전통 분자동역학 시뮬레이션은 원자 간의 페어 퍼텐셜(pair potential) 에너지를 묘사하는 힘 장 (force field)을 사용한다(Cygan, 2001). 힘 장은 원 자 간의 상호에너지를 경험적인 가변변수(empirical parameters)를 이용한 함수의 형태로 나타낸 것으 로 실험이나 양자역학 계산을 통해 얻어진다. 따라 서 시뮬레이션의 정확도는 힘 장의 정확도에 크게 좌우된다. 점토광물을 대상으로 한 힘 장은 미국 샌디아 국립연구소에서 개발한 ClayFF가 널리 사 용된다(Cygan et al., 2004). 이 힘 장은 점토광물 의 벌크(bulk) 결정구조, 광물 표면과 수용액 이온 사이의 상호작용의 연구에도 널리 사용될 정도로 점토광물 연구에서 그 신뢰성을 인정받고 있다 (Greathouse and Cygan, 2005; Newton et al., 2016; Kwon and Newton et al., 2016).

    점토광물은 사면체판(tetrahedral sheet)과 팔면체 판(octahedral sheet)이 적층된 구조를 가지며, 팔면 체판의 구조에 따라 이팔면체(dioctahedral) 또는 삼팔면체(trioctahedral) 광물로 구분된다. 이팔면체 점토광물의 팔면체판은 주로 3가 양이온(i.e., Al3+, Fe3+)으로 이루어지고 팔면체 세 자리 중 한 자리 가 비어 있으며, 삼팔편체 점토광물의 팔면체판은 주로 2가 양이온(i.e., Mg2+, Fe2+)으로 구성되며 모 든 팔면체 자리가 채워져 있는 구조를 갖는다 (Moore and Reynolds, 1997). 일반적으로 암석의 풍화과정에서 수화에너지가 작은 3가 양이온이 지 표면에 많이 잔류하기 때문에 이팔면체 점토광물 이 삼팔면체 점토광물보다 흔하게 발견된다(Velde, 1992). 그러나 다양한 삼팔면체 광물들이 종이, 페 인트, 플라스틱이나 고무 제조업 등의 다양한 분야 에 응용되고 전략금속인 리튬을 함유하고 있는 헥 토라이트는 리튬 광석으로 사용되기도 한다(Lien and Kramer, 1985; Komadel et al, 1996; U.S. Geological Survey, 2017).

    삼팔면체 광물 시뮬레이션에도 ClayFF가 사용 되지만(Braterman and Cygan, 2006; Wang et al., 2006; Larentzos et al., 2007; Greathouse et al., 2009; Zeitler et al., 2014), 이팔면체 광물을 주 대 상으로 개발되었기 때문에(Bougeard et al., 2000; Cygan et al., 2004; Sainz-Diaz et al., 2001; Teppen et al., 1997), 삼팔면체 광물 시뮬레이션을 위한 힘 장의 파라미터 사용이 명확하지 않다. 퍼 텐셜 에너지 계산에서 이온의 부분전하(partial charge, qi)가 매우 중요하며, ClayFF의 힘 장에서 는 팔면체 마그네슘 이온에 대한 부분전하(qMg)로 각각 ‘mgo’ 또는 ‘mgh’의 서로 다른 부분전하가 제안되어있다. 그런데, 마그네슘의 부분전하로 각 연구 마다 서로 다른 qMg를 사용한다. 활석 [Mg3Si4O10(OH)2]의 시뮬레이션에서는 주로 mgo 를 사용하여 계산된다(Wang et al. 2006; Du and Miller, 2007; Larentzos et al, 2007; Greathouse et al, 2009). 수활석[Mg(OH)2]의 시뮬레이션에서 는 Wang et al. (2004)는 mgo를 사용하지만, Braterman and Cygan (2006)Zeitler et al. (2014)는 mgh를 사용한다. 삼팔면체 1 : 1 광물인 리자다이트[Mg3Si2O5(OH)4]에 대한 MD 시뮬레이 션 결과는 아직 문헌에서 찾아볼 수 없으며, 리자 다이트 시뮬레이션을 통해 어떤 qMg이 적합한지에 대한 판단이 필요하다.

    이번 연구에서는 동형치환으로 인한 구조적 전 하를 띄지 않는 삼팔면체 점토광물인 수활석 (brucite), 리자다이트(lizardite), 활석(talc)을 대상 으로 ClayFF에서 제안된 마그네슘에 대한 두 종류 의 부분전하를 사용하여 분자동역학 시뮬레이션을 실시하였다. 이 부분전하들에 따른 시뮬레이션 결 과를 실험결과와 비교하였고, 삼팔면체 점토광물에 대한 시뮬레이션을 위해 어떤 파라미터를 사용해 야하는지 판단하였다.

    연구 방법

    결정 구조

    광물들의 결정구조 모델은 Fig. 1과 같이 American Mineralogist Crystal Structure Database에 업로드 된 자료를 사용하였다(Perdikatsis and Burzlaff, 1981; Mellini, 1982; Catti, 1995). 수활석은 단일 팔면체판으로 구성되며 육방정계에 속한다. X선 회절(X-ray diffraction, XRD) 분석법으로 측정된 모델의 격자상수는 a = b = 3.15 Å, c = 4.75 Å, 그리고 α = β = 90°, γ = 120°이다(Redfern and Wood, 1992). 리자다이트는 육방정계나 가상 (pseudo-)의 사방정계에 속하며, 사면체판과 팔면 체판이 1 : 1의 비율을 가진다. 중성자 회절(neutron diffraction) 분석법으로 측정된 리자다이트 모델의 격자상수는 a = b = 5.33 Å, c = 7.27 Å, 그리고 α = β = 90°, γ = 120°이다(Gregorkiewitz et al, 1996). 활석은 삼사정계나 단사정계에 해당하며, 활석은 사면체판과 팔면체판이 2 : 1의 비율로 결 합된 구조다. XRD로 측정된 활석 모델의 격자상수 는 a = 5.29 Å, b = 9.17 Å, c = 9.46 Å, 그리고 α = 90.5°, β = 98.7°, γ = 90.1°이다(Perdikatsis and Burzlaff, 1981).

    각 광물 모델들은 동형치환이 일어나지 않았음 을 가정하며 Mg, Si, O, OH의 원소로만 구성된 광물 단위포의 supercell을 사용하였다. 각 시뮬레 이션에서 사용된 supercell은 수활석 7 × 7 × 5 (1225 atoms), 리자다이트 4 × 4 × 3 (864 atoms), 활석 4 × 2 × 2 (672 atoms)이었다. Supercell의 결정축의 길이는 a, b, c 축 순서대로 수활석 22.05 Å × 22.05 Å × 23.85 Å, 리자다이트 21.33 Å × 21.33 Å × 21.70 Å, 활석 21.16 Å × 18.35 Å × 18.92 Å이었다. 이 광물들의 시스템이 supercell에 대한 주기적인 경계 조건을 넘어 대칭에 의한 제한 이 없도록 P1 대칭성(symmetry)을 가지도록 설정 한 후 계산하였다(Cygan et al., 2004).

    ClayFF 힘 장

    분자동역학 시뮬레이션은 주어진 퍼텐셜 에너지 (potential energy, U)를 이용해서 시간에 따른 구 성 원자들의 위치와 힘을 구하기 위해 뉴턴 방정식 을 수치적으로 계산한다(Allen, 2004). 원자(i = 1,2,3,…N)로 이루어진 시스템의 퍼텐셜 에너지는 식 (1)과 같이 각 원자들의 위치(ri)의 함수로 표현 된다.

    U = U ( r i ) = U ( r 1 , , r N )
    (1)

    원자에 가해지는 힘(f)이 보존된다고 할 때, 이 힘 은 퍼텐셜 에너지와 식 (2)와 같은 관계를 갖는다.

    f = m a = d U d r U = f d r
    (2)

    따라서 가속도(a)와 원자의 질량(m)의 뉴턴 방정 식을 풀어 힘의 크기를 알면 시간에 따른 원자의 궤적(trajectory)이 획득된다.

    원자들의 상호작용으로 생성되는 퍼텐셜 에너지 를 나타내기 위해 필요한 힘 장은 경험적인 가변변 수의 세트로 이루어진다. 이번 연구의 분자동역학 시뮬레이션에 사용되는 힘 장인 ClayFF는 광물 전 체 시스템의 전체의 퍼텐셜 에너지를 식 (3)과 같 이 표현한다(Cygan et al., 2004).

    U t o t a l = U C o u l + U V D W + U b o n d s t r e t c h + U a n g l e b e n d
    (3)

    식 (3)에서 UCoul은 쿨롱 에너지, UVDW는 반데르 발스 에너지를 나타내며, 이는 분자 간(intermolecular) 의 비결합 에너지(nonbonded energy)를 의미한다. Ubond-stretch는 공유 결합 된 분자 내 (intra-molecular)의 신축진동(stretching vibration) 을 표현하며, Uangle-bend은 분자의 굽힘 진동(bending vibration)을 표현한다. 이는 결합 에너지(bonded energy)를 의미한다(Cygan et al., 2004).

    쿨롱 에너지(UCoul)는 장거리(long-rang)까지 발 생하는 전하 간의 정전기적 상호작용 에너지로 서 로 다른 두 원자(i, j)의 부분전하의 곱에 비례하고 원자 간의 거리(rij)에 반비례한다(Cygan et al., 2004).

    U C o u l = e 2 4 π ε 0 i j q i q j r i j
    (4)

    식 (4)에서 e는 전자의 전하량(e = -1.602 × 10-19 C)이고 ε0은 진공상태의 유전율(ε0 = 8.85419 × 10-12Fm-1)이다. qiqj는 양자역학적으로 계산 된 원자 i와 j의 부분전하를 나타낸다.

    반데르발스 에너지(UVDW)는 다른 두 원자가 서 로 접근함으로써 증가하는 단거리(short-range) 에 너지와 관련된 척력(repulsive force)과 인력 (attractive force) 에너지를 나타낸다. 두 원자의 거 리(rij)가 서로 가까워질 때 일어나는 전자의 중첩 (electronic overlap)은 척력(repulsive force)을 발 생시키며 인접한 원자의 전자 밀도의 변동은 인력 (attractive force)을 발생시킨다. 일반적으로 원자 간의 거리와 반데르발스 에너지의 관계는 다음과 같이 Lennard-Jones (12-6) 함수(Jones, 1924; Allen and Tildesley, 2017)로 표현된다.

    U V D M = i j D 0 [ ( R 0 r i j ) 12 2 ( R 0 r i j ) 6 ]
    (5)

    식 (5)에 D0R0는 경험적인 가변변수로 D0는 원자 ij가 이루는 힘의 세기를 나타내는 가변변 수로 퍼텐셜 에너지 우물의 깊이(well depth)를 나 타낸다. R0는 최소 에너지를 갖는 원자 ij사이의 거리를 산술 평균(arithmetic mean)으로 계산한 가 변변수로 퍼텐셜 에너지 우물의 위치를 의미한다 (Lee, 2006). 지수함수 rij-12는 척력을 의미하고, rij-6는 인력을 의미하는 함수다.

    ClayFF에서는 수산기(hydroxyl group, OH)의 결합의 신축 진동(Ubond stretch)을 묘사를 위해 조화 퍼텐셜(harmonic potential)이나 모스 퍼텐셜(morse potential)을 사용한다. 조화 퍼텐셜을 사용한 수산 기의 산소와 수소 사이의 결합 에너지는 다음과 같 이 표현된다(Cygan, 2001).

    U b o n d s t r e c h = k 1 ( r r 0 ) 2
    (6)

    식 (6)에서 r은 수산기의 결합 거리이고 r0는 수 산기 내 작용하는 힘이 평형상태에 도달했을 때의 결합 거리이다. k1은 원자 내의 힘의 상수(force constant) 혹은 훅의 법칙(Hooke’s law)의 용수철 상수(spring constant)를 나타낸다. 이것은 원자의 결합 길이가 진동하면서 발생하는 진동 스펙트럼 을 분석하면 얻어진다. 적외선 분광법(infrared spectroscopy)이나 라만 분광법(raman spectroscopy) 에 의해 수산기의 진동 스펙트럼은 분석된 다. 반면에, 모스 퍼텐셜을 사용한 수산기의 결합 에너지는 다음과 같이 표현된다(Cygan, 2001).

    U b o n d s t r e c h = D o [ 1 exp { 1 α ( r r 0 ) } ] 2
    (7)

    식 (7)에서 rr0는 식 (6)에서와 같은 개념이고 D0은 결합 된 수산기가 이루는 퍼텐셜 에너지 우 물의 깊이(well depth)이고, α는 퍼텐셜 에너지 우 물의 폭(well width)이다. α는 평형상태인 퍼텐셜 에너지 우물의 힘 상수(k0)와 D0를 이용하여 식 (8) 처럼 계산된다.

    α = k 0 2 D 0
    (8)

    일반적으로 사용되는 조화 퍼텐셜로 계산된 원 자 간의 거리는 약 10% 내외의 오차 범위를 갖지 만, 계산이 단순하여 근삿값을 구하기 적합하다 (González, 2011). 모스 퍼텐셜은 보다 사실적인 (realistic) 결합 에너지의 묘사가 가능하지만, 계산 에 더 많은 시간이 소요되는 단점이 있다(Cygan, 2001). 이번 연구에서는 모스 퍼텐셜을 사용해서 MD 시뮬레이션을 계산하였다.

    서로 다른 세 원자 i, jk이 이루는 각도(θ)가 갖는 굽힘 에너지(Uangle-bend)는 일반적으로 조화 퍼 텐셜을 사용할 때 이용된다.

    U a n g l e b e n d = k 2 ( θ θ 0 ) 2
    (9)

    식 (9)에서 k2는 각도 굽힘의 힘의 상수를 나타 내고, θ0 는 평형상태의 원자 i, jk의 결합 각도 이다. 주로 광물의 수화된 표면에 금속이온 흡착묘 사를 위해 사용되며, 벌크 단위의 광물 시뮬레이션 에서는 Uangle-bend를 고려하지 않는다(Cygan et al., 2004).

    부분전하

    ClayFF에서 선택 가능한 마그네슘의 부분전하 (qMg)는 두 가지로 팔면체 마그네슘(octahedral magnesium, mgo, qmgo = -1.36 e)과 수산화물 마 그네슘(hydroxide magnesium, mgh, qmgh = -1.05 e) 이 존재한다(Cygan et al., 2004). 앞서 언급한대로 ClayFF 힘장은 이팔면체판 점토광물을 기반으로 개발되어, 팔면체 Al3+ 자리를 동형치환하는 Mg2+ 를 계산하기 위해서는 mgo가, 금속 수산화 광물 (hydroxide minerals)의 팔면체판의 Mg를 나타내 기 위해서는 mgh가 사용되어 왔다. 하지만, 최근 삼팔면체 광물 MD 시뮬레이션에서 수산화물 광물 인 수활석에 mgo와 mgh 모두 사용되었으며, 부분 전하에 따른 MD 시뮬레이션 차이 및 정확도가 명 확하게 정리되지 않았다(Braterman and Cygan, 2006; Wang et al., 2006; Zeitler et al., 2014). 따 라서 이번 연구에서는 삼팔면체 광물의 팔면체판 에 각각 mgo와 mgh를 설정하여 MD 시뮬레이션 결과의 차이점을 비교하여 계산의 정확도를 판단 하였다.

    마그네슘의 부분전하에 선택한 뒤, supercell의 전체전하가 “0” 되도록 일부 산소들의 부분전하를 수정하였다. 약 3% 이하의 범위의 전하 수정은 시 뮬레이션 결과에 중요한 영향을 미치지 않는 것으 로 예측된다(Wang et al., 2006; Cygan et al., 2004). 각 원자를 쉽게 구분하기 위해서 결정 구조 상의 위치에 따라 원자들의 이름을 세부적으로 다 르게 명명하였다(Fig. 2). 사면체판의 정점 산소 (apical oxygen)는 팔면체판과 연결되는 산소로 ‘Oa’, 서로 이웃한 Si-사면체들을 연결하는 산소 (basal oxygen)를 ‘Ob’, 팔면체에 연결된 상태로 층 의 바깥쪽(outer)을 향하는 수산기의 수소를 ‘Hout’, 이와 연결된 산소를 ‘Oout’, 팔면체에 연결된 상태 로 사면체판 빈 공간의 안쪽(inner)을 향하는 수산 기의 수소를 ‘Hin’, 이와 연결된 산소를 ‘Oin’의 기 호로 표현하였다(Fig. 2). 각 광물에서 사용한 자세 한 원소들의 부분전하는 Table 1에 정리하였다.

    시뮬레이션 방법

    이번 연구에서는 Materials Studio 2016 (BIOVIA Inc.)의 Forcite 모듈을 사용해서 시뮬레이션을 실 시하였다. 각 광물의 전체 구조를 안정화시키고 계 산의 오차를 최소화하기 위해서 구조최적화 (geometry optimization)를 우선 진행하였다. 이는 일반적으로 사용되는 스마트 알고리즘(smart algorithm)을 이용해 결정의 내부에너지가 2 × 10-5 kcalmol-1, 힘이 0.001 kcalmol-1-1, 응력이 0.001 Gpa, 거리가 1 × 10-5 Å 이내로 수렴할 때까지 구 조최적화를 계산하였다. UCoul는 Ewald summation 을 사용하였고, UVDW의 절단거리는 수활석, 리자다 이트에서 10 Å, 활석에서 9 Å으로 설정하여 계산 하였다.

    MD 시뮬레이션은 Wang et al. (2004)와 Greathouse et al. (2009)의 방법을 참고하여 정준 앙상블(canonical ensemble, NVT), 등온-등압 앙상 블(isothermal-isobaric ensemble, NPT)의 순서로 계산하였다. 온도(T), 압력(P), 부피(V), 구성 입자 (N), 에너지(E), 엔탈피(H) 등과 같은 상태 변수들 (state variables)이 주어진 분자 시스템의 반복적인 계산을 통해서 통계역학적인 평균 또는 앙상블 평균 에 도달한다. 온도는 Nosé-Hoover thermostat (Nosé, 1984a, 1984b, 1991)를 이용하여 300 K를 유지하고, 압력은 Parrinello-Rahman barostat (Parrinello and Rahman, 1981)를 이용하여 101,325 Pa (1기압)을 유지하였다. 각 광물 모델의 퍼텐셜 에너지가 평형 화(equilibration)되는 시뮬레이션 스텝이 달랐기 때문에 온도가 300 ± 20 K로 안정화될 때까지 계 산을 진행하였다. 대략적으로 광물마다 NVT, NPT 순서대로 각각 200 × 10-12 sec (ps) 이상의 시간이 소요되었다. 시간 간격(time step)은 1.0 × 10-15 sec (fs)로 설정하였고, Verlet 알고리즘에 적용해 서 원자의 다음 위치와 속도를 예측해 MD 시뮬레 이션을 수행하였다.

    온도와 압력이 일정하도록 NVT, NPT 앙상블을 계 산한 이후에 소정준 앙상블(microcanonical ensemble, NVE)을 1.0 fs의 시간 간격으로 100 ps 동안 계산 하였다. NVE 앙상블까지 계산된 광물 구조는 더 이상 온도와 압력에 영향을 받지 않는 일정한 구조 로 안정화되었다. 마지막 20 ps 동안 계산 된 구조 를 대상으로 원자 간의 거리나 수산기의 진동수와 같은 결정 구조를 분석하였다.

    동경분포함수

    동경분포함수(radial distribution function, RDF) 를 이용해서 시뮬레이션 결과로부터 원자 간의 거리 를 분석하였다. RDF는 쌍분포함수(pair correlation function)로 광물 시스템에서 하나의 원자를 선택 하고 일정한 거리(r)만큼 떨어져 있는 동심원 상의 영역에 존재하는 다른 원자들의 수를 계산한다 (Hansen and McDonald, 1990). 균질한 형태의 원 자 구조를 가진 경우 식 (10)과 같이 델타 함수 (delta function, δ ( r ) )로 RDF 계산이 가능하다.

    x ν x μ ρ g ν μ ( r ) = 1 N i = 1 N ν i = 1 N μ δ ( r + r i r j )
    (10)

    식 (10)에서 일정한 거리인 r 내에 존재하는 v와 μ는 서로 다른 화학적 결합을 갖는 원자들이고 x 는 해당 원자의 몰분율(mole fraction)이다. N은 시 스템의 모든 원자의 개수고 ρ는 시스템의 원자 밀 도이다. rirj는 기준 원자의 위치를 나타내며 vμ가 같을 때, ij가 같기 때문에 이 항(term)은 계산에서 제외된다.

    파워 스펙트럼

    일반적으로 원자들의 속도자기상관 함수(velocity autocorrelation function, VACF)를 주파수 함수로 푸리에 변환하여 파워 스펙트럼(power spectrum) 을 얻는다. 시간의존 상관 함수인 VACF는 시스템 의 분자동역학적 과정의 특성을 표현하는데 사용 되며 식 (11)과 같이 시간에 따른 원자 속도의 평 균을 이용한다.

    C ν ( t ) = ν i ( t 0 ) · ν i ( t 0 + n Δ t )
    (11)

    식 (11)에서 vi는 원자 i의 속도를 나타내는 스칼 라 값이며, t0는 초기 시간이고, △t은 시간 간격을 의미한다. 시간은 무한대로 계산될 수 있지만, 주 로 일정한 값(n)까지를 계산하여 Cv(t)를 구한다. VACF를 푸리에 변환한 결과를 파워 스펙트럼 밀 도(power spectrum density)라고도 하며 해당 주파 수 마다 단위 시간 당 해당하는 원자에 작용하는 에너지의 크기를 나타낸다.

    NVE 앙상블의 MD 시뮬레이션을 1.0 fs의 시간 간격(time step)으로 총 100 ps를 계산하였다. 그 후 마지막 20 ps 동안 시뮬레이션 궤적(trajectory) 을 이용하여 VACF를 구하고 푸리에 변환으로 파 워 스펙트럼을 계산하였다. 파워 스펙트럼의 해상 도는 약 2 cm-1이다.

    연구 결과 및 토의

    격자상수

    수활석, 리자다이트, 활석의 격자상수에 대한 시 뮬레이션 계산결과와 실험결과를 Table 2에 정리 하였다. 시뮬레이션으로 얻어진 격자상수는 qMg의 종류에 따라 미세한 차이가 있었지만, 특정한 전하 가 더 좋은 계산결과를 보인다고 단정할 수는 없었 다. 수활석의 경우, ab축은 qmgo를 사용할 때, c 축은 qmgh를 사용했을 때 실험결과와 더 유사하였 다. 각도는 qMg에 상관없이 실험결과와 일치하였 다. 리자다이트의 경우 a, b, c축 모두 qmgh일 때 실험결과와 유사했지만, qmgh로 시뮬레이션을 수행 했을 때 실험결과와 각도의 차이가 존재했다. 활석 의 경우 a, b축은 어떤 qMg에 상관없이 실험결과와 유사했고 c축은 qmgh의 시뮬레이션 결과가 실험결 과에 더 일치하였지만, 각도의 경우에는 qmgo를 사 용해서 계산할 때 실험결과와 오차가 더 작았다.

    광물의 종류에 상관없이, c축의 경우 ab축에 비해 실험결과와 더 큰 차이를 보였다. 각도의 경 우는 γ를 제외하고 각도 αβ의 시뮬레이션 결과 가 변동 폭이 컸고 실험결과와 더 큰 차이를 보였 다. 이런 양상은 광물마다 다르지만, 시뮬레이션 결과와 실험결과는 약 0-4%의 오차 범위에서 일치 하였다. 결과적으로 qMg의 종류에 상관없이 시뮬레 이션 결과는 실험결과와 비슷했고 어느 한쪽의 부 분전하가 더 실험결과와 일치하는 결과를 계산하 는지는 명확하게 판단되지 않았다.

    원자 간의 길이

    원자 간의 길이를 RDF로 분석한 결과 MD 시뮬 레이션의 부분전하 값에 따라 원자 간의 길이가 조 금씩 다르게 계산되었다(Table 3). Mg-O 결합거리 는 qmgo를 사용했을 때, 그리고 Si-O 결합거리는 qmgh를 사용했을 때 실험결과와 더 유사하게 계산 되었다. 수산기의 결합거리(O-H)는 qMg의 종류에 상관없이 일정하게 계산되었다. Mg-O 및 O-H는 실험결과나 밀도 범함수 이론(density functional theory, DFT) 계산결과보다 크게 계산되었지만, Si-O는 MD 시뮬레이션 결과에서 더 작게 계산되 었다. 이러한 경향성은 ClayFF를 사용한 다른 MD 시뮬레이션 결과에도 나타난다(Larentzos et al., 2007; Zeitler et al., 2014).

    모든 삼팔면체 점토광물에서 계산된 Mg-O 결합 거리는 qmgo일 때보다 qmgh일 때 더 길게 계산되었 다. qmgo의 부분전하가 qmgh보다 크기 때문에, 거리 의 변화는 폴링의 법칙(Pauling’s rule)으로 설명할 수 있다. Si-Ob 결합거리는 qMg의 종류에 상관없이 거의 동일하게 계산되었지만, Si-Oa 결합거리는 부 분전하가 작은 qmgh을 사용했을 때 조금 더 길게 계산되었다. 팔면체판의 Mg에 직접 결합한 Oa까 지는 qMg의 영향을 매우 작게나마 받지만, 사면체 판의 Ob까지는 qMg의 영향이 미치지 않은 것으로 보인다.

    수산기는 모스 퍼텐셜로 계산했을 때 qMg의 종류 에 상관없이 0.97 Å의 일정한 길이로 계산되었다. 이는 실험결과보다는 작은 값이다. 수산기의 거리는 Ubond stretch에 의해 크게 영향을 받기 때문에 qMg의 종류에 따른 영향은 거의 받지 않는 것을 확인할 수 있었다. 단, qmgo로 계산한 리자다이트의 경우 수산 기의 위치에 따라 O-Hin는 0.97 Å이고, O-Hout는 0.99 Å으로 계산되었다. 하지만 모든 수산기에 대 한 평균 O-H 거리는 0.97 Å로 계산되었다.

    수산기 진동 파워 스펙트럼

    모스 퍼텐셜을 사용해서 계산한 수활석, 리자다 이트 그리고 활석의 수산기 파워스펙트럼은 Fig. 3 과 같이 나타났다. 광물마다 차이가 있지만, 계산 한 진동수는 모든 광물에서 qmgo를 사용했을 때가 qmgh를 사용했을 때 보다 평균 90 cm-1 정도 작게 계산되었다. 이를 통해 부분전하의 값이 클수록 진 동수가 작게 계산되는 것을 확인하였다.

    실험결과와 수산기 진동수의 계산결과를 qMg에 관해서 비교하였다. 수활석의 경우, 적외선(IR) 분 광분석 실험에서 수산기 신축 진동모드는 3698 cm-1이지만(Frost and Kloprogge, 1999), 이번 qmgo 계산에서는 3492 cm-1, qmgh 계산에서는 3561cm-1 로 분석되었다(Fig. 3a). 실험값에 비해 100-200 cm-1 정도 작은 4-5%의 차이가 있었다. 활석의 경 우, IR 실험 수산기 진동수가 3670 cm-1인데 (Farmer, 1958), qmgo 계산에서는 실험값 보다 약 40 cm-1 작은 3626 cm-1로, qmgh 계산에서는 3720 cm-1로 50 cm-1 큰 결과가 얻어졌다(Fig. 3c).

    리자다이트에는 두 종류의 수산기(OHin과 OHout) 가 존재한다(Fig. 2). 이번 연구에서 계산한 수소의 진동 파워 스펙트럼에서도 두 종류의 수산기를 잘 구분해 보여준다(Fig. 3b). qmgo를 사용한 경우 수 산기의 진동수가 각각 ~3530 cm-1과 ~3700 cm-1 로, qmgh를 사용한 경우에는 ~3670 cm-1과 ~3760 cm-1로 다소 넓은 두 정점으로 계산되었다. 두 정 점 중 낮은 진동수는 OHout의 진동수에 해당하는 데 이는 Hout은 적층된 옆 사면체층의 산소와 수소 결합으로 OHout 결합세기가 약하기 때문이다. 높은 진동수는 OHin의 진동수에 해당한다(Fig. 3b). Balan et al. (2002)의 IR 실험에서 3584 cm-1과 3684 cm-1을 중심으로 다소 넓은 두 진동 피크가 관찰되었는데, 이는 이번 연구에서 qmgo를 이용하 여 얻은 진동 양상과 일치한다. 하지만, 실험에서 관찰된 3584 cm-1 정점은 팔면체판에 동형치환된 Fe2+에 의한 것이라(Fuchs et al., 1998)는 가정 아 래, Balan et al. (2002)는 순수 리자다이트 IR 실 험결과에 3645 cm-1과 3684 cm-1 (평균 3664 cm-1) 을 OHout 진동모드로, 3703 cm-1을 OHin의 진동모 드로 해석하였다. 이 실험결과를 고려하면, OHout 진동모드는 qmgh를 사용했을 때(실험 평균 3664 cm-1과의 차이 : -134 (mgo) vs. +6 (mgh) cm-1) 실험결과와 더 일치하였으나, OHin의 진동모드는 qmgo를 사용했을 때(실험 3703 cm-1과의 차이 : -3 (mgo) vs. +57 (mgh) cm-1) 실험결과와 더 일치하 였다. 이는 OHout만을 가지는 수활석과 OHin만을 가지는 활석에 대한 계산 경향도 일치한다.

    결과적으로 모스 퍼텐셜로 계산한 수산기의 파 워 스펙트럼은 qMg의 종류에 의한 진동수의 차이 가 있었다. 수활석은 qmgh를 사용한 경우 실험결과 에 근접하게 계산되었고 약 4%의 오차를 나타내었 다. 활석은 qmgo를 사용해서 계산한 경우가 약 1% 의 오차를 보이며 실험결과와 유사하게 계산되었 다. 리자다이트 경우, OHout 진동모드는 qmgh을 사 용했을 때(< 1% 오차), OHin 진동모드는 qmgo을 사용했을 때(< 1% 오차) 실험 결과에 더 근접하게 계산되었다. 따라서, MD 시뮬레이션에서 모스 퍼 텐셜을 이용하여 수산기를 계산하기 위하여 qMg의 종류에 대한 신중한 선택이 필요하다. 이번 연구와 같이 단일 팔면체판의 수산기(OHout)를 갖는 수활 석은 mgh를, 사면체판과 양쪽으로 결합하는 팔면 체판의 수산기(OHin)를 갖는 활석은 mgo를 사용할 수 있다. 리자다이트와 같이 팔면체판의 수산기가 두 가지(OHout & OHin)인 경우 mgh와 mgo 모두 고려될 수 있기에 우선하는 수산기의 종류에 따라 마그네슘 부분전하를 선택하여 사용할 수 있다.

    결 론

    삼팔면체 점토광물인 수활석, 리자다이트, 활석 에 대한 MD 시뮬레이션을 ClayFF 힘 장(force field)을 이용하여 실시하였다. 팔면체판에 위치한 Mg 이온의 전하로써 두 가지 부분전하(qmgoqmgh)가 사용될 수 있다. 다른 부분전하를 사용했을 경우 격자상수, 원자 간의 거리, 수산기의 진동 파 워 스펙트럼 시뮬레이션 결과의 차이점을 체계적 으로 테스트하였다. 부분전하에 따라 광물의 격자 상수나 원자 간의 거리 등의 광물구조 결과에 수치 적인 차이는 있었지만, 그 차이가 시뮬레이션 정확 도에 영향을 줄 정도로 크지 않았고 뚜렷한 규칙성 을 보여주지 않았다. 전체 에너지의 일차 미분(first derivative)을 사용하는 결정구조 계산과는 달리, 진동수 계산은 전체 에너지의 이차 미분(second derivative)에 의존하기 때문에 마그네슘의 부분전 하에 따라 계산결과가 크게 차이가 날 수 있다. 예 상한대로 모스 퍼텐셜을 사용해서 수산기의 진동 수를 계산했을 때는 마그네슘의 부분전하에 따라 결과에 뚜렷한 차이가 있었다. 따라서 삼팔면체 점 토광물의 MD 시뮬레이션에서 결정구조는 Mg 부 분전하의 종류에 따라 크게 관계가 없지만, 수산기 의 진동모드 계산은 부분전하에 대한 신중한 고려 가 필요하다. 이번 연구의 계산 결과는 이팔면체 점토광물을 대상으로 개발된 ClayFF 힘 장이 삼팔 면체 점토광물에 활용 가능성을 보여주며, 헥토라 이트의 리튬 치환/치윤 기작 연구와 같은 삼팔면체 광물에 대한 다양한 시뮬레이션 연구에 가장 기본 이 되는 파라미터 선정에 기준이 될 수 있다.

    사 사

    이 연구는 2017년도 정부(미래창조과학부)의 한국연 구재단(No. NRF-2015R1A4A1041105)과 2016년도 강 원대학교 대학회계 학술연구조성비(520160469)를 지원 받아 수행되었다. 원고 향상에 도움을 주신 익명의 심자 자께 감사드린다.

    Figure

    JMSK-30-161_F1.gif

    Crystal structures of (a) brucite, (b) lizardite, (c) talc.

    JMSK-30-161_F2.gif

    Atom types (Hin vs. Hout; Oin vs. Oout; Oa vs. Ob) used in the current study.

    JMSK-30-161_F3.gif

    Calculated hydrogen vibrational power spectra of (a) brucite, (b) lizardite, and (c) talc using an atom type, mgo or mgh, for Mg in the MD simulations. ▼ represents the experimental IR data for hydroxyls (Frost and Kloprogge, 1999; Balan et al., 2002; Farmer, 1958).

    Table

    Two type of partial charges in |e| (qmgo or qmgh) for trioctahedral clay minerals

    *Wang et al. (2004).
    †Zeitler et al. (2014).
    ‡Wang et al. (2006).

    Lattice parameters of brucite, lizardite and talc calculated with different partial charges

    exp : experimental result; * : standard deviation.
    *Redfern et al. (1992).
    †Gregorkiewitz et al. (1996).
    ‡Perdikatsis and Burzlaff (1981).

    Bond distances (in Å) obtained by the MD simulations using different charge groups (qmgo or qmgh) along with experimental and other MD simulation results

    exp : experimental result; DFT : density functional theory; ± : full width at half maximum.
    aCatti et al. (1995).
    bZeitler et al. (2014) : MD (mgh, morse potential), DFT (GGA).
    cBaranek et al. (2001) : DFT (B3LYP).
    dMellini and Viti (1994).
    ePrencipe et al. (2009) : DFT (B3LYP).
    fPerdikatsis and Burzlaff (1981) : DFT (B3LYP).
    gLarentzos et al. (2007) : MD (mgo, harmonic potencial), DFT (GGA).

    Reference

    1. AllenM.P. , AttigN. (2004) Computational Soft Matter: From Synthetic Polymers to Proteins, John von Neumann Institute for Computing, ; pp.1-28
    2. AllenM.P. TildesleyD.J. (2017) Computer simulation of liquids., Oxford university press,
    3. BaileyS.W. (1980) Summary of recommendations of AIPEA nomenclature committee on clay minerals. , Am. Mineral., Vol.65 ; pp.1-7
    4. BalanE. SaittaA. MauriF. LemaireC. GuyotF. (2015) First-principles calculation of the infrared spectrum of lizardite. , Am. Mineral., Vol.87 ; pp.1286-1290
    5. BougeardD. SmirnovK.S. GeidelE. (2000) Vibrational spectra and structure of kaolinite: A computer simulation study. , J. Phys. Chem. B, Vol.104 ; pp.9210-9217
    6. BratermanP.S. CyganR.T. (2006) Vibrational spectroscopy of brucite: A molecular simulation investigation. , Am. Mineral., Vol.91 ; pp.1188-1196
    7. CattiM. FerrarisG. HullS. PaveseA. (1995) Static compression and H disorder in brucite, Mg (OH) 2, to 11 GPa: a powder neutron diffraction study. , Phys. Chem. Miner., Vol.22 ; pp.200-206
    8. CostanzoP.M. (2001) Baseline studies of the clay minerals society source clays: Introduction. , Clays Clay Miner., Vol.49 ; pp.372-373
    9. CyganR.T. (2001) Molecular modeling in mineralogy and geochemistry. , Rev. Mineral. Geochem., Vol.42 ; pp.1-35
    10. CyganR.T. LiangJ.J. KalinichevA.G. (2004) Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field. , J. Phys. Chem. B, Vol.108 ; pp.1255-1266
    11. DuH. MillerJ.D. (2007) A molecular dynamics simulation study of water structure and adsorption states at talc surfaces. , Int. J. Miner. Process., Vol.84 ; pp.172-184
    12. DubeA. ZbytniewskiR. KowalkowskiT. CukrowskaE. BuszewskiB. (2001) Adsorption and migration of heavy metals in soil. , Pol. J. Environ. Stud., Vol.10 ; pp.1-10
    13. EwaldP.P. (1921) Die Berechnung optischer und elektrostatischer Gitterpotentiale. , Ann. Phys., Vol.369 ; pp.253-287
    14. FarmerV.C. (1958) The infrared spectra of talc, saponite and hectorite. , Mineral. Mag., Vol.31 ; pp.829-845
    15. FrostR.L. KloproggeJ.T. (1999) Infrared emission spectroscopic study of brucite. , Spectrochim. Acta A Mol. Biomol. Spectrosc., Vol.55 ; pp.2195-2205
    16. FuchsY. LinaresJ. MelliniM. (1998) Mössbauer and infrared spectrometry of lizardite-1T from Monte Fico, Elba. , Phys. Chem. Miner., Vol.26 ; pp.111-115
    17. GonzA lezM. A. (2011) Force fields and molecular dynamics simulations. , École thématique de la Société Française de la Neutronique, Vol.12 ; pp.169-200
    18. GreathouseJ.A. CyganR.T. (2005) Molecular dynamics simulation of uranyl (VI) adsorption equilibria onto an external montmorillonite surface. , Phys. Chem. Chem. Phys., Vol.7 ; pp.3580-3586
    19. GreathouseJ.A. DurkinJ.S. LarentzosJ.P. CyganR.T. (2009) Implementation of a Morse potential to model hydroxyl behavior in phyllosilicates. , J. Chem. Phys., Vol.130 ; pp.134713
    20. GregorkiewitzM. LebechB. MelliniM. VitiC. (1996) Hydrogen positions and thermal expansion in lizardite-1T from Elba: A low-temperature study using Rietveld refinement of neutron diffraction data. , Am. Mineral., Vol.81 ; pp.1111-1116
    21. GrimR.E. (1968) Clay Mineralogy., McGraw-Hill,
    22. HansenJ.P. McDonaldI.R. (1990) Theory of simple liquids., Academic Press,
    23. JonesJ.E. (1924) On the determination of molecular fields. II. From the equation of state of a gas. , Proceedings of the Royal Society, Vol.106 ; pp.463-477
    24. KomadelP. BujdA kJ. MadejovA J. uchaV. ElsassF. (1996) Effect of non-swelling layers on the dissolution of reduced-charge montmorillonite in hydrochloric acid. , Clay Miner., Vol.31 ; pp.333-345
    25. KwonK.D. NewtonA.G. (2016) Structure and stability of pyrophyllite edge surfaces: Effect of temperature and water chemical potential. , Geochim. Cosmochim. Acta, Vol.190 ; pp.100-114
    26. LarentzosJ.P. GreathouseJ.A. CyganR.T. (2007) An ab initio and classical molecular dynamics investigation of the structural and vibrational properties of talc and pyrophyllite. , J. Phys. Chem. C, Vol.111 ; pp.12752-12759
    27. LeeJ.G. (2006) Computational Materials Science: Introduction, Young., Uiwang,
    28. LienR.H. KramerD.A. (1985) Recovery of Lithium from a Montmorillonite-type Clay., US Department of the Interior, Bureau of Mines,
    29. MelliniM. (1982) The crystal structure of lizardite 1T: hydrogen bonds and polytypism. , Am. Mineral., Vol.67 ; pp.587-598
    30. MelliniM. VitiC. (1994) Crystal structure of lizardite-1T from Elba, Italy. , Am. Mineral., Vol.79 ; pp.1194-1198
    31. MooreD.M. ReynoldsR.C. (1989) X-ray Diffraction and the Identification and Analysis of Clay Minerals., Oxford university press,
    32. NewtonA.G. KwonK.D. CheongD.K. (2016) Edge structure of montmorillonite from atomistic simulations. , Minerals (Basel), Vol.6 ; pp.25
    33. NosA(c)S. (1984) A molecular dynamics method for simulations in the canonical ensemble. , Mol. Phys., Vol.52 ; pp.255-268a
    34. NosA(c)S. (1984) A unified formulation of the constant temperature molecular dynamics methods. , J. Chem. Phys., Vol.81 ; pp.511-519b
    35. NosA(c)S. (1991) Constant temperature molecular dynamics methods. , Prog. Theor. Phys. Suppl., Vol.103 ; pp.1-46
    36. U.S. Geological SurveyU.S. Geological Survey (2017) https://minerals.usgs.gov/minerals/pubs/mcs/
    37. ParrinelloM. RahmanA. (1981) Polymorphic transitions in single crystals: A new molecular dynamics method. , J. Appl. Phys., Vol.52 ; pp.7182-7190
    38. PerdikatsisB. BurzlaffH. (1981) Strukturverfeinerung am Talk Mg3[(OH)2Si4O10]. , Zeitschrift für Kristallographie- Crystalline Materials, Vol.156 ; pp.177-186
    39. RedfernS.A. WoodB.J. (1992) Thermal expansion of brucite, Mg(OH)2. , Am. Mineral., Vol.77 ; pp.1129-1132
    40. Sainz-DiazC.I. Hernandez-LagunaA. DoveM.T. (2001) Modeling of dioctahedral 2 : 1 phyllosilicates by means of transferable empirical potentials. , Phys. Chem. Miner., Vol.28 ; pp.130-141
    41. SavageD. (1995) The scientific and regulatory basis for the geological disposal of radioactive waste., Wiley,
    42. SomaY. SomaM. (1989) Chemical reactions of organic compounds on clay surfaces. , Environ. Health Perspect., Vol.83 ; pp.205
    43. TeppenB.J. RasmussenK. BertschP.M. MillerD.M. SchäferL. (1997) Molecular dynamics modeling of clay minerals. 1. Gibbsite, kaolinite, pyrophyllite, and beidellite. , J. Phys. Chem. B, Vol.101 ; pp.1579-1587
    44. TuckerM.R. (1999) http://www.ncagr.gov/agronomi/uyrst.htm
    45. VeldeB. (1992) Introduction to clay minerals: chemistry, origins, uses and environmental significance., Chapman and Hall Ltd.,
    46. WangJ. KalinichevA.G. KirkpatrickR.J. (2004) Molecular modeling of water structure in nano-pores between brucite (001) surfaces. , Geochim. Cosmochim. Acta, Vol.68 ; pp.3351-3365
    47. WangJ. KalinichevA.G. KirkpatrickR.J. (2006) Effects of substrate structure and composition on the structure, dynamics, and energetics of water at mineral surfaces: A molecular dynamics modeling study. , Geochim. Cosmochim. Acta, Vol.70 ; pp.562-582
    48. ZeitlerT.R. GreathouseJ.A. GaleJ.D. CyganR.T. (2014) Vibrational analysis of brucite surfaces and the development of an improved force field for molecular simulation of interfaces. , J. Phys. Chem. C, Vol.118 ; pp.7946-7953