서론
유리는 오랜 기간 동안 다양한 산업 분야에서 핵심적인 소재로 활용되어 왔다. 특히, 우수한 열적 안정성, 화학적 내구성, 광학적 특성 덕분에 디스플레이, 광섬유, 바이오소재, 에너지 저장 장치 등 첨단 기술 분야에서도 널리 사용되고 있다.[1,2] 최근에는 유리가 기존의 반도체 기판을 대체할 수 있는 소재로 주목받으며, 차세대 전자소자 및 패키징 기술에서 그 중요성이 점차 부각되고 있다.[3] 유리는 고체전해질(Solid electrolyte)로서의 가능성 측면에서도 많은 관심을 받고 있다. 유리계 고체 전해질은 유연한 조성 설계가 가능하고, 높은 기계적⋅화학적 안정성과 더불어 입계 저항(Grain boundary resistance)이 없다는 점에서 뛰어난 계면 안정성을 갖는다. 이러한 장점으로 인해, 차세대 리튬 이차전지, 나트륨 전지, 전고체 배터리(Solid-state battery) 등의 응용 분야에서 활발히 연구되고 있으며, 유리 상태에서 높은 이온 전도도를 확보할 수 있는 다성분 산화물계, 황화물계, 할로겐계 유리에 대한 연구도 꾸준히 확대되고 있다.[4-6] 이러한 배경 속에서 유리 재료에 대한 이해는 과거보다 더욱 중요한 연구 주제로 자리 잡고 있다. 고기능성 유리 소재의 설계 및 응용을 위해서는 유리의 물성과 기능을 결정짓는 유리 구조에 대한 완벽한 이해가 선행되어야 한다. 유리는 장범위 규칙도(Long-range order)가 존재하지 않는 비정질(Amorphous) 상태이며,[7] 원자 배열이 무질서하게 분포되어 있기 때문에 결정질 재료처럼 X-ray diffraction 등의 실험 기법만으로는 그 구조를 명확히 규명하는 데 한계가 존재한다. 유리의 구조는 완전히 무작위적으로 분포되어 있지는 않으며, 단범위 규칙도(Short-range) 및 중간범위 규칙도(Medium-range)를 포함하고 있으며, 이는 유리의 물리적 특성에 중요한 영향을 준다.[8]
이에 따라, 유리 구조를 예측하기 위한 다양한 이론적 예측 모델이 제안되었고, 구조 분석을 위한 측정 기법도 지속적으로 발전해왔다. M. B. Volf 등은 1988년, 유리의 화학 조성만을 이용하여 물리적 특성을 예측할 수 있는 1차 및 2차 혼합 모델을 제안하였다.[9] 1차 모델은 적용 가능한 조성 범위가 제한적이지만, 해당 범위 내에서는 높은 정확도를 보이며, 2차 모델은 조성 성분 간의 상호작용을 고려함으로써 보다 넓은 조성 영역에서의 예측이 가능하다. 이러한 모델들은 유리 내 국소 원자 구조(Local structure units)를 기반으로 하며, 구조 단위(Structure unit)의 개념을 활용하여 인접 원자 간의 배치를 정의한다.[10] 구조 단위는 열역학 계산, 구조적 평형 가정, 핵자기공명(NMR, Nuclear magnetic resonance) 분석 등을 통해 도출된다. K. Budhwani와 S. Feller는 이러한 구조 단위를 활용하여 유리의 밀도를 예측할 수 있는 모델을 개발했으며,[10] 이후 H. Inoue et al.은 M. B. Volf의 모델을 구조 단위 기반으로 확장하여 Na2O-B2O3-SiO2 유리 시스템의 다양한 물리적 특성을 예측할 수 있는 모델을 제시하였다.[11] 해당 모델은 W. J. Dell과 P. J. Bray가 1983년에 최초로 제안한 구조 단위 정의에 기반하고 있으며,[12,13] K. Budhwani와 S. Feller의 계산법을 사용해 BO3/2, BO4/2Na, BO2/2ONa, BO1/2O2Na2, SiO4/2 (Q4), SiO3/2ONa (Q3), SiO2/2O2Na2 (Q2), SiO1/2O3Na3 (Q1)의 구조 단위를 정의하였다. [SiO2]/[B2O3] = K, [Na2O]/[B2O3] = R이라는 조성비를 바탕으로 0<R<Rmax = 0.5 + 0.0625 K(영역 1), Rmax<R<R1=0.5+0.25 K(영역 2), R1<R<R2=1.5+0.75 K (영역 3), R2<R<R3=2+ K(영역 4)를 정의하였다.[11] 각 영역별로 비가교 산소(NBO, Nonbridging oxygen)의 수를 계산하고 NMR 분석 결과와 조합하여, Qn unit(Q4, Q3, Q2, Q1)의 분포를 추정할 수 있도록 하였다.[11] 해당 계산법을 통해 각 구조영역에 대한 구조 단위의 몰 분율과 SiO4 사면체당 비가교 산소의 수를 계산할 수 있으며, 각 계산에 대한 방정식은 Table 1에 요약하였다. L. S. Du와 J. F. Stebbins는 2005년에 Na2O-B2O3-Al2O3-SiO24성분계에 대한 예측 모델을 제안하였고, W. J. Dell and P. J. Bray 모델의 핵심 지표인 BO3/BO4 비율(N4)을 유지하면서 Al2O3가 BO4 단위 형성에 기여함을 고려하여 Fig. 1에서 나타난 바와 같이 N4=[B2O3]× N4 +[Al2O3] / [B2O3]+[Al2O3]로 확장하였다.[14] L. S. Du와 J. F. Stebbins 모델은 ZnO, ZrO2 등 Na2O 이외의 modifier 가 포함될 경우 예측 정확도가 낮아지는 한계를 가진다. 이러한 문제를 보완하기 위해 X. Lu 등은 2021년에 기존 L. S. Du 와 J. F. Stebbins 모델을 확장하여 Na2O 이외에도 Li2O, CaO, ZnO 등의 modifier와 Al2O3, ZrO2 등의 network former를 고려할 수 있는 일반화된 모델을 제안하였다.[15] 이 모델은 조성 정보만으로도 Borate glass에서의 BO4 unit 비율(N4)을 예측할 수 있으며, 관련 구조 인자를 머신 러닝 모델의 입력 변수로 사용할 수 있도록 정량화하였다. 이와 같이 유리 구조의 예측을 위한 다양한 이론적 모델이 제시되었으나, 조성이 복잡해질수록 예측의 범용성과 신뢰도가 저하된다는 한계가 존재한다. 대부분의 모델들은 실험적으로 도출된 인자에 기반하고 있으며, 구조 단위를 간접적으로 예측할 수 있어 실제 유리의 정밀한 구조를 완전히 규명하기에는 한계가 따른다.
Fig. 1.
Comparison of calculated N’4 values ([B2 O3] × N4 + [Al2O3] / [B2O3] + [Al2O3]) of R’. Reproduced from L.S. Du et al. J. Non-Cryst. Solids. 2005;351:43-45 with permission of ELSEVIER.[14]

Table 1.
Mole fractions of the Structural unit and the number of Nonbridging Oxygen per SiO4 Tetrahedron. Reproduced from Inuoe et al. J. Am. Ceram. Soc: 2012;95:211-16, with permission of The American Ceramic Society.[11]
최근 Synchrotron XRD[16] 및 Neutron diffraction,[17] NMR,[18] Extended X-ray Absorption Fine Structure [19] 등 고분해능 3차원 분석 기법들의 빠른 발전으로 유리 구조 분석의 정밀도는 크게 향상되었다. 이러한 기법들은 유리 내 원자 배치와 국소 구조에 대한 정보를 실험적으로 획득하는데 매우 유용하지만, 유리의 무질서한 네트워크 구조나 구성 원소 간의 상호작용 등 미시적 수준의 복잡한 구조적 특성을 완전히 규명하기에는 여전히 한계가 존재한다.[20] 실험적 분석기법과 이론적 모델은 각각의 강점을 지니고 있으나, 독립적으로는 유리 구조를 완전하게 설명하기에는 불충분하다.[21]
이에 따라 유리 구조 해석의 대안으로 원자 수준의 시뮬레이션(Atomistic simulation)이 주목받게 되었다. 이 시뮬레이션 기법은 유리 재료의 통계역학적 상태 공간(Statistical mechanical phase space)을 수치적으로 탐색할 수 있는 수단을 제공하며, 이 과정을 통해 생성된 구조와 물성은 상태 공간상의 여러 지점들을 평균한 결과로 간주할 수 있다.[22] 분자동역학 시뮬레이션(Molecular dynamics simulation)과 몬테카를로 시뮬레이션(MC, Monte carlo simulation)은 유리 연구에 실질적으로 많이 활용되고 있으며, 비록 MC 시뮬레이션 기법이 MD 시뮬레이션 기법보다 먼저 개발되었지만, 동역학 정보를 함께 제공할 수 있는 MD 시뮬레이션 기법이 유리 연구에 더 널리 사용되고 있다.[2] 유리에 대한 최초의 MD 시뮬레이션은 1976년 L. V. Woodcock 등에 의해 수행되었으며, Born-Mayer-Huggins 포텐셜을 사용하여 실리카 유리를 대상으로 연구되었다.[23] 이후 실리카 및 알칼리 실리케이트 유리에 대한 다수의 MD 연구들이 발표되었으며, SiO2-Na2O,[24] Na2O-CaO-SiO2,[25] Na2O-Al2O3-SiO2,[26] Al2O3-Y2O3-SiO2,[27] SiO2-Na2O-CaO-P2O5,[28] SiO2-B2O3-Na2O-CaO,[29] Al2O3-P2O5-SiO2-CeO2,[30] Li2O-P2O5,[31,33] Na2O-B2O3[32] glass system 등 보다 복잡하고 현실적인 유리 시스템에 대한 M D 시뮬레이션이 가능해졌다. A. Masuno 등은 Born-Mayer 포텐셜을 사용하여 Na2O-B2O3-SiO2 유리를 대상으로 MD 시뮬레이션을 수행하였다. 이 연구에서는 Na2O 함량에 따른 BO4 단위의 비율(N4) 변화를 성공적으로 재현하였으며, W. J. Dell의 모델과 비교하였다. 또한, 비가교 산소 형성 경향 또한 실험 결과와 비교하여 합리적으로 설명하였다.[34] Fig. 2에 나타난 바와 같이, 이 연구에서는 SiO4 단위 내 Qn 구조(Q1∼Q4)의 분포 변화, B–O 및 Si–O 네트워크 내 산소 연결성 변화를 분석함으로써, 유리의 중간 범위 구조 특성까지 체계적으로 규명하였다. 이러한 진전은 시뮬레이션 기법의 발전뿐만 아니라 다양한 포텐셜 모델(부분 전하 모델, 분극화 모델, 반응형 포텐셜 등)의 개발과 고성능 컴퓨터 자원의 향상된 접근성 덕분이다.[22] MD 시뮬레이션은 일반적으로 주기적 경계 조건이 적용된 원자 혹은 분자 집합체를 대상으로 수행되며, 원자 간 상호작용은 경험적 포텐셜을 통해 계산된다. 원자의 초기 위치는 무작위로 배치하거나 결정 구조로부터 시작하며, 속도는 볼츠만 분포에 따라 초기화된다. 그 이후 원자들은 뉴턴의 운동 방정식을 매 femtosecond 단위로 계산해 가며 이동하고, 이 과정을 반복하며 전체 시뮬레이션 시간을 구성하게 된다. 시뮬레이션은 일반적으로 NVT, NPT, NVE와 같은 다양한 앙상블 조건하에서 수행되며, 이로부터 압력, 열용량 등 다양한 열역학적 물성을 도출할 수 있다. 유리 구조 생성을 위한 가장 일반적인 방식은 가열-용융-냉각 과정을 거치는 것으로, 시뮬레이션에서는 표면에서부터 녹는 자연스러운 용융이 어려우므로 일반적인 녹는점보다 훨씬 높은 5000 K 이상의 고온이 필요하다.[22,35] 이후 선형 또는 비선형 냉각을 통해 유리 구조가 형성되며, 이때 저장된 원자 궤적은 구조 분석 및 확산 분석 등에 활용된다. 현재 MD 시뮬레이션은 학계와 산업계 모두에서 실질적인 문제 해결을 위해 활발히 활용되고 있으며, 그 사례로 디스플레이용 화학 강화 유리의 이온 교환 과정을 Mean square displacement 분석을 통해 이해하려는 연구가 있다.[36-38] 대부분의 산업 및 기술적으로 중요한 유리는 다성분 산화물 유리(Multicomponent oxide glasses)이며, 유리는 넓은 조성 범위에서도 유리 상태를 유지할 수 있는 특징이 있어 조성 조절을 통한 물성 최적화가 가능하다. 이러한 재료에 대한 시뮬레이션 모델 구축의 중요성에 따라 다양한 산화물계 유리에 적용하려는 연구가 활발히 진행되고 있다.[39]
Fig. 2.
(a) Comparison of N4 values from MD simulations and Dell's model for various K = [SiO2]/[B2O3] ratio (0.5, 1.0, 1.5). Reproduced from Inuoe et al. J. Phys. Chem B: 2012; 116:12325-31, with permission of American Chemical Society.[11]

본 리뷰에서는 Molecular Dynamics 시뮬레이션에 사용되는 포텐셜 함수, 앙상블 조건, 시뮬레이션 알고리즘 및 계산 기법, 시뮬레이션 프로그램과 분석 도구에 대해 체계적으로 설명하고, 다성분계 oxyfluoride glass system의 구조 분석 사례를 소개하고자 한다.[40] 시뮬레이션을 통해 얻어진 구조적 특성을 실험 결과와 비교함으로써 구조 모델의 타당성과 신뢰성을 평가하였으며, 원자 수준에서의 구조 분석 및 시뮬레이션 기반 구조 해석의 유용성을 제시하고자 한다.
본론
2.1 Fundamentals of MD simulation analysis
2.1.1 Interatomic potentials and force field
Molecular dynamics 시뮬레이션에서 경험적 포텐셜 또는 힘 장(force field)은 원자 간 상호작용을 기술하는 핵심요소로, 시뮬레이션의 정확성과 신뢰도를 좌우하는 중요한 입력 변수이다. 이 포텐셜 함수는 원자 또는 이온의 위치에 따라 시스템의 포텐셜 에너지를 계산하며, 이 값을 바탕으로 뉴턴의 운동 방정식을 수치적으로 적분하여 시간에 따른 원자의 움직임을 추적할 수 있다. 따라서 포텐셜 함수의 수학적 형태, 포함된 항, 피팅 방식에 따라 시뮬레이션 결과의 정밀도와 물리적 타당성은 큰 영향을 받게 된다. 특히, 특정 유리 조성 시스템이 시뮬레이션 가능 여부 또한 해당 조성에 적합한 포텐셜의 존재에 달려 있다.[22,35] 일반적으로 시스템의 총 포텐셜 에너지는 원자의 위치 벡터의 함수로 나타내며, 이는 1체(one-body), 2체(two-body), 3체(three-body), 또는 필요 시 그 이상의 고차 항의 합으로 구성된다.[35] 이 중 1체 항은 외 부 전기장 등에 의해 발 생 하는 장거리 쿨롱 상호작용(Long-range Coulombic interaction)을 포함한다. 쿨롱 상호작용은 전하를 띤 입자들 사이의 전기적인 힘을 의미하며, 두 입자 사이의 거리 r이 멀어질수록 그 세기는 1/r에 비례하여 천천히 감소한다. 이처럼 천천히 감소하는 거리 의존성 때문에, 쿨롱 상호작용을 정확하게 계산하기 위해서는 모든 입자 쌍 간의 상호작용을 고려해야 하므로 계산량이 너무 많아지고, 그 결과 계산이 느리게 수렴한다. 이러한 문제를 해결하기 위해 Ewald summation이나 Particle Mesh Ewald 와 같은 효율적인 계산 기법이 활용된다.[22] 2체 항은 원자 쌍 사이의 결합 상호작용을 설명하며, 짧은 거리에서는 반발력, 특정 거리 이상에서는 반데르발스(van der Waals) 인력 등의 효과를 포함한다. 이 항은 원자 간 거리 r에 따라 급격히 변하는 것이 특징이다. 3체 이상의 항은 결합각이나 다면체의 기하학적 제약 조건을 유지하기 위해 포함된다. 이온성이 높은 비정질 고체인 유리의 경우, 이러한 상호작용 중에서도 이온 간 결합과 전하 간 인력이 지배적이기 때문에, 장거리 쿨롱 상호작용과 단거리 상호작용을 모두 포함한 포텐셜이 주로 사용된다.[35] 대표적으로 Lennard-Jones,[41,42] Morse,[43] Buckingham,[44] Born-Mayer-Huggins[45] 포텐셜 등이 있으며, 각각은 다양한 유리 시스템에 적용되어 왔다. Lennard-Jones 포텐셜은 초기 2체 상호작용 모델 중 하나로, 단순한 비극성 분자나 기체 상태에 적합하도록 설계되었으며, 식 1과 같이 표현된다.
ε은 potential well의 깊이이며, σ는 포텐셜 에너지 E가 0이 될 때의 거리를 나타내고, rij는 원자 간 거리를 정의하는 매개변수이다. r12 항은 반발력, r6 항은 반데르발스 인력을 나타낸다. Morse 포텐셜은 분자의 진동 특성을 보다 현실적으로 반영하기 위해 개발된 2체 포텐셜로, 결합 거리 근처의 에너지를 정밀하게 묘사할 수 있으며, 식 2에 나타내었다.[22,35]
반발력은 지수 함수 형태로 표현되며, Dij는 potential well 깊이, αij는 포텐셜 폭 변수이며, r0는 평형 결합 거리 변수, rij는 원자 간 거리를 나타낸다. Buckingham 포텐셜은 식 3에 나타난 바와 같이 Lennard-Jones 포텐셜과 Morse 포텐셜의 특성을 결합한 형태로, 지수형 반발력과 r-6 인력항을 결합하여 나타낸다.[22,35]
A는 반발상수, ρ는 반발 범위 변수, C는 분산상수, r은 원자 간 거리로 나타낸다. 이 포텐셜은 산화물 유리나 세라믹 같은 이온성 고체에 널리 사용되지만, 짧은 거리에서 인력 항이 과도하게 증가하여 원자가 비물리적으로 융합되어 시뮬레이션이 불안정해지는 문제가 발생한다. 이를 해결하기 위해 Ziegler-biersack-littmark (ZBL) 포텐셜,[45] Splice function[46] 등이 도입되었다. 특히 ZBL 포텐셜은 고에너지 이온 충돌을 다루는 시뮬레이션에서 효과적인 짧은 거리 반발력을 제공한다. J. Du 등은 Buckingham 포텐셜의 안정화를 위해 splice function을 도입하여, 특정 거리 r0에서 포텐셜과 그 도함수가 연속되도록 설계하였다. 이는 시뮬레이션의 물리적 안정성을 크게 향상시킨다. Born-Mayer-Huggins 포텐셜은 식 4에 나타난 바와 같이 Buckingham potential에 r-8 항을 추가하여 인력을 강화한 형태로 복잡한 이온성 구조에 대해 더욱 정확한 묘사가 가능하다.
보다 복잡한 유리 시스템을 정확히 시뮬레이션하기 위해서는 단순한 2체 포텐셜을 넘어서 조성비, 배위수, 결합각 등의 물질의 화학적 특성을 반영한 정교한 모델이 필요하다. D. M. Teter가 개발하고 A. N. Cormak과 J. Du 등이 수정한 포텐셜은 이러한 요구를 반영하여 실리카, 알루미노실리카, 인산염 유리 등 보다 복잡한 유리 시스템에 성공적으로 적용된 바가 있으며, 많은 구조 모델들에 활용이 되고 있다.[39,47] 이 포텐셜 모델은 부분 이온 전하를 기반으로 하며, 산소 이온에는 −1.2의 전하를, 양이온에는 조성비에 따라 비례적으로 조정된 전하가 할당한다. 이 모델의 단거리 상호작용은 식 5에 나타난 바와 같이 Buckingham 포텐셜의 형태를 따르며, 장거리 상호작용은 쿨롱 항으로 구현된다.
A는 반발상수, ρ는 반발 범위 변수, C는 분산상수, r는 원자 간 거리, Zij는 입자 i, j의 전하, ε0는 자유공간 내 유전율로 정의되며, 본 포텐셜의 매개변수들은 Table 2[39]에 요약하였다.
Table 2.
Atomic charge and Buckingham potential parameters for oxide glasses. Reproduced from Du et al. Springer, Cham: 2015;157-80, with permission of Springer Nature.[39]
A. Pedone 등은 기존 모델에 바탕으로 단거리 반발력을 보완하기 위해 Morse 포텐셜에 r-12항을 추가한 포텐셜을 개발하였다.[48] 이 포텐셜은 실리케이트 광물의 결정 구조 및 물성 데이터를 기반으로 피팅되어 다양한 다성분 유리 조성에 적합하게 적용된다. 복잡한 다성분 유리의 경우, 조성에 따라 이온의 배위수가 변할 수 있다. 이러한 조성과 배위수의 상관관계를 반영하기 위해, A. N. Cormack과 B. Park는 배위 의존 포텐셜을, L. Huang과 J. Kieffer는 전하 평형 기반의 포텐셜을 제안하였다.[50] L. H. Kieu 등의 연구에서는 고정 전하 방식으로도 조성 의존성을 반영할 수 있음을 제시하였으며,[51] H. Inoue 등은 단순한 2체 포텐셜과 고정 전하 모델만으로도 신뢰성 있는 결과를 얻을 수 있음을 보였다.[34] 그러나 여전히 유리형성제가 여러 개 존재하는 시스템에 대해 보편적으로 적용 가능한 고정밀 포텐셜 개발은 진행 중인 연구 과제로 남아 있다.
한편, 화학 반응, 결합 형성/파괴, 계면 반응 등을 정확히 모사하기 위해서는 전하의 시간적 변화를 반영하는 가변 전하 및 결합 차수 기반 포텐셜이 필요하다. 대표적으로 Reactive force-field (ReaxFF)[52,53]와 Charge optimized many body potential [54] 포텐셜이 있으며, 이들은 시뮬레이션 중 원자의 전하를 실시간으로 업데이트하여 동적 화학 반응을 재현할 수 있도록 설계되었다. Rea x FF는 전기음성도 평형화 기법(Elec-tronegativity equalization method)을 사용하여 원자의 부분 전하를 주변 환경에 따라 조정함으로써 복잡한 화학 반응 메커니즘을 정확하게 재현할 수 있다. 이온 분극은 유리의 전기적, 구조적 특성에 큰 영향을 미치는 요소이며, 이를 설명하기 위해 코어-셸 모델이 자주 사용된다.[55] 이 모델에서는 중심 이온(코어)과 질량이 없는 셸(shell)이 스프링으로 연결되어 있으며, 셸의 전하 및 스프링 상수는 이온의 분극율을 반영하도록 조정된다. A. Tilocca 등은 이 모델을 사용하여 소다라임 실리케이트 및 생체활성 유리의 Qn 분포 및 실리콘 배위 결함을 정밀하게 재현하였다.[56] 마지막으로, Ab initio 기반의 분자동역학 시뮬레이션(Ab Initio Molecular Dynamics, AIMD)은 경험적 포텐셜 없이 양자역학적 계산을 기반으로 원자 간 상호 작용을 직접 계산하기 때문에 가장 보편적이고 정확한 시뮬레이션 방법으로 간주된다.[57] 그러나 AIMD는 계산 비용이 매우 크기 때문에 수백 개의 원자 이하, 수십 피코 초(ps) 수준 이내의 시간 범위에만 적용 가능하며, 다성분 유리 시스템에는 적용이 제한적이다. 그럼에도 불구하고 AIMD는 고전적 MD 시뮬레이션 결과를 검증하는데 중요한 기준 역할을 하며, 특히 복잡한 결합 특성을 지닌 시스템에서는 거의 유일한 시뮬레이션 도구로 활용된다.[39]
2.1.2 Computational algorithms
분자동역학 시뮬레이션은 주어진 원자들의 초기 위치와 속도를 바탕으로 시간에 따른 시스템의 상태를 추적해 나가는 방식이다. 이는 뉴턴의 제2운동 방정식에 기반하며, 특정한 시간 간격 Δt(타임 스텝) 동안, 원자에 작용하는 힘이 주어지면, 그 시간 동안의 운동을 수치적으로 계산할 수 있다. 이 과정을 반복함으로써 원자의 시간에 따른 궤적을 생성하게 되며, 이러한 정보를 기반으로 다양한 물리적 및 열역학적 성질을 추출할 수 있다.[35] 이러한 궤적 생성을 위해서는 운동 방정식의 수치적 적분이 필수적으로 요구되며, 이 역할을 수행하는 것이 바로 적분 알고리즘이다. 일반적으로 원자 X개로 이루어진 시스템의 경우 총 3 X개의 2차 미분방정식을 해결해야 하며, 정확도 뿐만 아니라 계산 효율성, 에너지 보존, 안정성을 고려한 다양한 알고리즘이 개발되어왔다.[2] 대표적인 알고리즘으로는 Verlet,[22,58] Leapfrog,[22,59] Velocity Verlet[22,60] 알고리즘이 있으며, 이 외에도 목적에 따라 다양한 변형 알고리즘이 존재한다.
Verlet 알고리즘은 가장 오래되고 널리 사용되어온 MD 적분 방식 중 하나로, 위치에 대한 테일러 급수 전개 기반으로 한다. 이 방식에서는 t시점과 t-Δt 시점의 위치 정보를 이용하여 t+Δt 시점의 위치를 계산하므로, 속도를 직접 계산하지 않아도 된다. 결과적으로 구현이 간단하여 메모리 사용량이 적고 계산적으로도 매우 안정적이라는 장점을 갖는다. 그러나, 속도가 직접적으로 계산되지 않기 때문에, 운동 에너지나 점도와 같은 속도 기반 물성의 경우, 후속 단계에서 변형된 식을 통해 추정해야 하며, 이로 인해 속도의 정확도가 떨어지고 오차가 누적될 수 있는 한계가 존재한다. 또한 두 시점의 위치를 저장해야 하므로, 과거 컴퓨팅 자원이 제한된 환경에서 메모리 부담이 문제되기도 했다.[22] Leapfrog 알고리즘은 Verlet 방식의 변형으로, 속도와 위치를 서로 반 타임스텝 차이로 계산하는 구조를 가진다. 속도는 t+1/2Δt 시점에서 먼저 계산되고, 이를 바탕으로 위치는 t+Δt 시점에서 계산된다. 그 후, 속도는 t+3/2Δt 지점에서 업데이트되는 형식으로 반복된다. 이 방식은 속도를 명시적으로 계산하므로 운동 에너지 보존 특성이 우수하며, 상대적으로 큰 타임스텝에서도 안정적으로 동작하는 특징이 있다. 속도와 위치가 동기화되어 있지 않기 때문에, 특정 시점에서의 총 에너지(운동 에너지+포텐셜 에너지)를 정확하게 계산하기 어렵다는 제한점이 있으며, 총 에너지를 정밀하게 추적해야 하는 경우 사용에 제약이 있다.[59] Velocity Verlet 알고리즘은 기존 Verlet 방식의 단점을 보완하기 위해 고안된 알고리즘으로, 위치와 속도, 힘을 동일한 타임 스텝에서 명시적으로 계산할 수 있다는 점에서 현재 가장 널리 사용되는 방식이다. 이 알고리즘의 핵심은 하나의 타임 스텝 내에서 위치와 속도를 함께 업데이트하면서, 힘 계산은 한 번만 수행한다는 점이다. 즉, 시간 t에서의 위치와 속도를 바탕으로 새로운 위치 r(t+Δt) 계산을 하고, 이 위치에서의 새로운 힘 F(t+Δt) 또는 가속도 a(t+Δt)를 계산한다. 이후 t와 t+Δt 시점의 가속도를 평균화하여 속도 v(t+Δt)를 업데이트하는 순서로 계산된다.[60] 이 방식은 운동 에너지와 포텐셜 에너지를 동일 시간에서 계산할 수 있으며, 시간 가역성과 에너지 보존 측면에서도 뛰어난 특성을 가진다. 따라서 이 방식은 LAMMPS (Large-scale atomic molecular massively parallel simulator),[61] DL_POLY,[62] GROMACS[63] 등 대표적인 MD 패키지에서 기본 옵션으로 채택되어 있으며, 속도 기반 물성을 다루는 시뮬레이션에서 특히 효과적이다.[22]
이 외에도 여러 알고리즘이 개발되어 있으며, 그중 하나인 Beeman 알고리즘[64]은 Verlet 알고리즘의 변형 중 하나로, 특히 속도 계산의 정밀도 향상과 에너지 보존에 초점을 맞춘 방식이다. 위치 계산 시 현재와 이전 시점의 가속도를 함께 활용하고, 속도 계산에는 t, t-Δt, t+Δt 시점의 가속도를 모두 반영한다. 그 결과, 속도 계산의 정밀도가 높고 장시간 시뮬레이션에도 에너지 누적 오차가 작은 장점이 있다. 그러나 추가적인 가속도 저장 및 계산이 필요하며, 알고리즘 구현이 상대적으로 복잡하기에, 계산 비용이 높은 유리 시스템이나 다성분계 시뮬레이션에서는 거의 사용되지 않는다. Predictor-Corrector 알고리즘은[65] 다음 시간 스텝의 위치, 속도, 가속도를 테일러 전개를 통해 먼저 예측한 뒤, 해당 위치에서의 실제 힘을 계산하여 결과값을 보정하는 방식이다. 초기에는 높은 정확도를 위해 사용되었으나, 예측-보정 과정이 반복적이며 계산 시간이 길고, 시스템 조건에 따라 수치적 불안정성이 나타날 수 있는 단점이 있다.[22] 현재는 대부분의 MD 시뮬레이션에서는 사용되지 않으며, 특정한 고정밀 계산 환경에서만 제한적으로 활용된다. 최근 주목받고 있는 Stochastic integration[66]은 운동 방정식 적분과 동시에 온도 및 압력 제어 기능을 통합한 방식으로, 열역학 앙상블 조건을 직접 샘플링 할 수 있도록 설계되었다. 이러한 방식은 소규모 시스템, 유체 및 생체 시스템, 혹은 정확한 앙상블 분포 유지가 필요한 경우에 효과적이며, 최근에는 GULP[67]와 같은 프로그램에서 기본 적분기로 채택되기도 한다. 다만, 확률적 요소로 인해 궤적의 부드러움이 떨어지고, 동역학적 특성 분석에는 다소 제약이 따를 수 있다.
2.1.3 Ensemble
앞서 언급하 바와 같이 분자동역학 시뮬레이션은 통계역학에 기반한 계산 기법으로, 원자 또는 분자의 궤적을 생성하여 시스템의 미시 상태를 샘플링하고, 이를 통해 다양한 거시적 물성을 도출한다. 이때 시뮬레이션이 수행되는 열역학적 조건을 정의하는 것이 바로 앙상블이며, 이는 동일한 입자 수, 에너지, 압력, 온도 등의 거시 변수를 공유하는 미시 상태들의 통계적 집합으로 이해할 수 있다. MD 시뮬레이션은 단일 궤적을 계산하지만, Ergodic hypothesis[68]에 따라 충분한 시간 동안 생성된 궤적의 시간 평균값이 앙상블 평균값과 일치한다고 간주한다. 이러한 전제 하에서, 시뮬레이션 시스템은 평형에 도달한 이후 가능한 상태 공간을 충분히 탐색해야 하며, 이를 통해 도출된 시간 평균값은 통계적으로 의미 있는 물리량으로 해석된다.[35] 따라서 연구 목적 및 시스템 특성에 따라 적절한 앙상블의 선택은 매우 중요한 고려 요소이며, 일반적으로 MD 시뮬레이션에서는 NVE, NVT, NPT 앙상블이 활용된다.
NVE 앙상블은 입자 수(N), 체적(V), 에너지(E)가 일정하게 유지되는 앙상블로, 외부 환경과의 에너지 교환이 없는 고립계를 가정한다. 뉴턴의 운동 방정식에 기반한 에너지 및 운동량 보존이 자연스럽게 이루어지기 때문에, 대부분의 MD 시뮬레이션 세팅에서 기본적으로 적용되며, 보통 시스템이 열역학적 평형에 도달한 후 단계에서 구조 및 확산 특성 등의 계산에 사용된다.[22] 현실적인 물리 조건을 고려할 때는, 에너지를 기준으로 제어하지 않고 대체로 온도를 지표로 사용하기 때문에, 입자 수(N), 체적(V), 온도(T)가 일정한 NVT 앙상블이 가장 널리 사용된다.[22] 이 앙상블에서는 시스템이 가상 열 저장소와 연결되어 평균 운동 에너지를 조절함으로써 일정한 온도가 유지된다. 가장 간단한 온도 제어 방식은 입자 속도를 일정한 비율로 재조정하는 속도 재스케일링 방식이며, 대표적으로 Berendsen thermostat이 사용된다.[69,70] 이 방식은 시스템의 온도를 목표 온도 T0에 수렴시키기 위해 식 6과 같이 온도 변화율을 조절한다.
T(t)는 현재 시스템 온도, τ는 시스템의 온도 완화를 결정하는 완화 시간이다. 이 방식은 빠른 수렴 속도라는 장점을 가지나, 생성되는 속도 분포가 엄밀한 NVT 앙상블 분포에서 벗어날 수 있다는 한계가 있다. 보다 이론적으로 정밀한 온도 제어를 위해서는 Nosé-Hoover thermostat가 사용된다.[71,72] 이 방법은 시스템에 확장된 자유도를 추가하여, 식 7과 같은 형태의 운동 방정식으로 온도 제어를 수행한다.
pξ는 온도조절계 변수의 운동량이며, Q는 열 저장소와 시스템 간 에너지 교환의 강도를 결정하는 관성 매개변수이다. 이 방식은 이론적으로 정확한 NVT 앙상블을 구현할 수 있으며, 특히 장시간 안정적인 유리 시스템에 자주 사용된다. 확률론적 방식의 Andersen thermostat 는 일정 시간 간격으로 무작위로 선택된 일부 입자의 속도를 맥스웰-볼츠만 분포에 따라 재할당하는 방식으로, 이론적으로 NVT 앙상블을 정확히 재현할 수 있다.[70,73] 하지만 궤적의 연속성이 저해되므로, 동역학 특성 분석에는 다소 부적절할 수 있다. NPT 앙상블은 입자 수(N), 압력(P), 온도(T)가 일정하며, 시스템의 부피가 자유롭게 변화할 수 있는 조건을 가진 앙상블이다. 이는 구조 변화, 상전이, 열팽창 등 압력 의존적 현상을 모사하기 위한 시뮬레이션에서 필수적인 앙상블이며, 유리 시스템의 밀도 조정이나 응력 상태 분석에서 중요한 역할을 한다.[22] 압력은 일반적으로 비리얼 상태 방정식(Virial equation)을 기반으로 계산되며, 포텐셜 에너지의 거리 및 힘에 대한 미분 항에 민감하게 작용하기 때문에, 온도나 에너지보다 더 큰 시간적 변동성을 갖는다. 압력 제어 방식으로는 Berendsen barostat이 대표적이며,[69] 식 8과 같이 압력 변화율이 표현된다.
τp는 압력 완화 시간, P0는 목표 압력을 의미한다. 이때 시스템의 좌표 μij 는 식 9와 같은 식으로 동적으로 선형 스케일링 된다.
β는 시스템의 등온 압축률을 나타내며, 시스템이 목표 압력에 점진적으로 평형 상태에 수렴하도록 조정하는 역할을 한다. 시뮬레이션에서 사용되는 포텐셜 모델에 따라 시스템이 수렴하는 최종 밀도가 다를 수 있으므로, 실험적으로 측정된 밀도 조건에 맞춰 시뮬레이션을 수행하는 경우가 많다. 특히 비정질 유리의 경우, 결합각이나 원자 간 평균 거리와 같은 밀도 의존적 구조 특성이 크기 때문에, 보다 정확한 비교와 예측을 위해 실험 밀도 조건에서 시뮬레이션을 고정시키는 전략이 중요하게 고려된다.[22]
2.1.4 Simulation programs & tools
분자동역학 시뮬레이션의 발전은 계산 이론 및 물리 모델링의 정립뿐만 아니라, 이를 실질적으로 구현하고 가속화할 수 있는 소프트웨어 및 하드웨어의 기술적 진보에 크게 의존해왔다. 유리 시스템과 같은 비정질 재료는 높은 구성 복잡도와 다성분성을 가지기 때문에, 이를 효과적으로 다룰 수 있는 고성능 계산 도구 및 환경이 필수적이다. 현재 유리 시스템의 원자 구조와 동역학적 특성을 모사하기 위해 가장 보편적으로 사용되는 MD 시뮬레이션 패키지로는 LAMMPS, DL_POLY, GROMACS, GULP 등이 있다. LAMMPS는 미국 샌디아 국립연구소에서 개발된 오픈소스 기반 병렬 MD 코드로, 다양한 유형의 포텐셜 함수와 앙상블 조건을 유연하게 지원한다. 대규모 원자 시스템에 대한 병렬 연산이 가능하여,유리의 구조 형성, 이온 확산, 계면 거동 등 광범위한주제에 활용되고 있다.[61] DL_POLY는 영국 Daresbury 연구소에서 개발된 범용 MD 패키지로, 고전적인 포텐셜모델을 기반으로 강체 및 연체 시스템 모두를 지원한다.다양한 출력 옵션과 효율적인 재시작 기능을 제공하며,이온성 유리의 열역학적 특성 분석에 유리한 플랫폼으로 평가된다.[66] GROMACS는 원래 생체분자 시뮬레이션을 목적으로 개발되었으나, 최근에는 무기 비정질재료로의 적용 범위가 확대되고 있다. 높은 계산 효율성과 직관적인 사용자 인터페이스를 바탕으로, 비교적 단일 목적의 시뮬레이션이나 구조적 초기화에 많이 활용된다.[63] GULP는 결정성 고체 및 격자 기반 시스템에특화된 시뮬레이터로, 격자 진동 분석 및 탄성 계수 계산등에서 우수한 성질을 보이며, 결정-비정질 전이 모델링에 활용 가능하다.[67] 이러한 프로그램들은 대부분오픈소스 또는 학술 연구용 무료 라이선스로 배포되며, MPI (Message passing interface)를 활용한 병렬 연산을 통해 수천∼수만 개 원자 규모의 시뮬레이션을 효율적으로 수행할 수 있다.
MD 시뮬레이션으로부터 얻어진 궤적 데이터를 시각적으로 분석하고 구조 정보를 정량화하기 위해 다양한 후처리 및 시각화 도구가 사용된다. 대표적인 예시로는 VMD (Visual molecular dynamics), VESTA, OVITO (Open visualization tool) 그리고 Materials studio가 있다. VMD는 시간에 따른 원자 구조의 변화를 3 D로 시각화 할 수 있으며, 원자 간 거리 계산, 결합각 분석, 표면 구조 렌더링 등의 기능을 통해 구조적 특성 파악에 효과적이다.[74] VESTA는 결정 및 비정질 구조 모두를 고해상도로 시각화 할 수 있는 툴로, 결합 길이, 다면체 시각화 등 각종 구조 분석에 유용하며, 실험 데이터와 비교에도 자주 활용된다.[75] OVITO는 오픈소스 기반의 시각화 도구로, 시뮬레이션 결과에서 특정 구조 단위를 자동 식별하고 시간에 따라 시각화 할 수 있는 기능을 제공한다. 특히 비정질 구조 분석에 최적화되어 있어 유리 시뮬레이션 후처리에서 자주 활용된다.[76] Materials studio는 Biovia에서 제공하는 상용 통합 시뮬레이션 플랫폼으로, 모델링, 포텐셜 설정, 시뮬레이션 수행, 결과 해석까지 하나의 Graphical user interface (GUI) 환경에서 일괄 처리할 수 있는 상용 통합 플랫폼으로 상당히 활용성이 높다. 특히, GUI 기반 워크플로우 설정의 용이성과 내장된 고정밀 포텐셜 라이브러리를 통해 산업 및 응용 연구에서 자주 활용된다.[77] 이러한 도구들은 계산 결과의 직관적 해석과 실험 데이터와의 비교를 위한 시각적 근거를 제공하며, 시뮬레이션 결과에 대한 해석의 정밀도 확보에 핵심적인 역할을 수행한다. 계산 효율성, 데이터 처리 속도, 사용자 편의성 등도 시뮬레이션 수행의 전반적인 품질과 생산성을 결정짓는 중요한 요소로 작용하므로, 연구 목적에 따라 적절한 시뮬레이션 프로그램과 분석 도구를 선택하는 것이 필수적이다.
2.2 Structure analysis of Multicomponent Oxyfluoride glass system
본 리뷰에서는 다성분계 유리의 MD 시뮬레이션을 통한 구조분석 사례로, BaF2-BaO-La2O3-B2O3계 oxyfluoride 다성분계 glass system에 대한 구조 분석 결과를 소개하고자 한다.
2.2.1 Oxyfluoride glass
Oxyfluoride glass는 산소(O)와 불소(F)를 동시에 음이온으로 포함하는 복합 비정질 소재로, 전통적인 산화물 유리(Oxide glass)와 불화물 유리(Fluoride glass)의 장점을 모두 가지는 고기능성 유리로 주목받고 있다.[78-84] 일반적으로 F 원자는 O 원자보다 높은 전기음성도와 낮은 분극성을 가지므로, 유리의 굴절률 분산 특성(Abbe 수)을 향상시킬 수 있으며, 이는 파장 분산 제어나 고굴절 광학 설계가 요구되는 분야에서 유리한 특성으로 작용한다.[78-81] 1993년 Y. Wang과 H. Ohwaki는 PbxCd1-xF2 결정상을 포함하고 Er3+및 Yb3+ 이온이 도핑된 oxyfluoride glass ceramic을 처음으로 보고하였다.[76] 이 시스템은 산화물 기반 매트릭스의 우수한 열적, 기계적 안정성과 불화물 결정상의 우수한 발광 특성을 동시에 갖춤으로써 차세대 광기능성 소재로서 큰 관심을 받아왔다.[81] 이후 LaF3, PbF2, CaF2 등 다양한 불화물 결정상을 유리 내 나노 크기로 안정화하고 희토류 이온을 도핑함으로써, 고효율 광 증폭기, 레이저, 업컨버전 소자 등으로의 응용 가능성이 확장되고 있다.[79]
특히 결정 크기가 가시광 파장 이하로 유지될 경우 산란 손실 없이 높은 투명성이 확보되며, 낮은 독성과 높은 화학적 안정성, 수분 흡착에 의한 발광 소멸 억제 특성 등은 실용적인 광학 소재로서의 강점을 뒷받침한다.[79-85] 그러나 제조 과정 중 SiF4 및 GeF4와 같은 휘발성 화합물의 형성으로 인해 조성 제약이 존재하며, 이를 극복하기 위한 방안으로 P2O5, B2O3, Nb2O5 등의 산화물이 유리 안정화 및 유리 형성을 위해 효과적으로 활용되고 있다. 최근 연구들은 불소와 산소를 동시에 포함하는 유리계의 가능성에 주목하고 있으며, 일부 조성에서는 산소가 오히려 유리 형성 안정성을 향상시키는 역할을 수행함이 밝혀졌다. S. A Polishchuk 등은 LnF3– MFₓ– SiO2 (또는 GeO2) 기반의 시스템에서는 희토류 이온의 농도가 60 mol% 이상에 이를 수 있으며, 이는 유리의 굴절률 및 발광 특성 향상에 기여할 수 있다고 정의하였다. 이러한 조성에서는 M–F 결합보다 M–O 결합의 이온성이 낮기 때문에, 희토류 이온은 기존 산화물 또는 불화물 기반 시스템과는 다른 고유의 화학 환경을 형성하게 된다.[74] 이러한 배경 속에서 BaO-La2O3-B2O3 기반 산화물 유리 시스템에 BaF2를 첨가한 fluoroborate 유리 시스템의 분자동역학 시뮬레이션을 통한 구조 모델을 리뷰 하였다. X-ray 회절로부터 얻어진 총 상관함수 (Total correlation fuction)와 19F MAS-NMR 분석을 기반으로, BaF2의 첨가가 유리 구조에 미치는 영향을 평가하였다. 특히 B-O 네트워크의 연결성, BO4 단위 비율(N4 값), 결합각 분포, 비가교 산소의 수, Qn 구조 단위의 분포 및 중간범위 질서에서 나타나는 간극 구조 등에 미치는 영향을 평가한 사례를 소개하고자 한다.[40]
2.2.2 Experiment and Molecular dynamics simulation details
MD 시뮬레이션을 통해 35 Ba O–15 La2O3–50B2O3 (BLB0), 20BaF2–15 BaO–15 La2O3–50B2O3 (BLB20) 그리고 35 BaF2–15 La2O3–50B2O3 (BLB35)(mol%) 3종의 유리에 대한 구조 모델을 확보하였으며, 각 유리조성별 구조적 특성을 분석하였다. 19F MAS-NMR 측정을 위해 20BaF2–30BaO–50B2O3 유리(BB20)가 추가로 제조되었다. 시료의 원자간 평균 구조 분석을 위해, X-ray 회절 측정을 기반의 총 상관함수 분석이 수행되었으며, 실험적 총 상관함수 Texp(r)과 분자동역학 시뮬레이션으로부터 얻은 Tcal(r)를 비교하여 모델의 신뢰도를 검증하였다.
분자동역학 시뮬레이션은 Born-Mayer 형태의 pairwise 포텐셜을 기반으로 수행되었으며, 원자 간 상호작용 에너지 Φij는 식 10과 같이 정의되며, 해당 식에 사용된 포텐셜 파라미터는 Table 3에 요약되어 있다.
Table 3.
Parameter of the Born-Mayer Potentials. Reproduced from Chung et al. J. Phys. Chem B: 2023;127:3091-9, with permission of American Chemical Society.[40]
Zi, Zj는 유효 전하량, Bij, ρij는 각각의 반발 상수와 연질성 매개변수이며, 장거리의 쿨롱 상호작용은 Ewald summation 방법을 이용하여 계산되었다. 각 원소에 대한 부분 전하는 실험적 데이터를 바탕으로 Ba, La, B, F 원자에 대해 +1.7, +.255, +1.8, −0.8 그리고 O 원자는 전하 중성을 유지되도록 설정되었다.
시뮬레이션에 사용된 구조 모델은 약 3000개의 원자로 구성된 정육면체 셀을 기반으로 하였으며, 주기적 경계 조건(Periodic boundary conditions)이 적용되었다. 셀의 크기는 실험적 밀도를 반영하여 설정되었으며, 초기 원자 배치는 무작위로 구성되었다. 시뮬레이션은 4,000 K에서 40,000 step 동안 평형화를 수행한 후, 총 50,000 step에 걸쳐 293 K까지 천천히 냉각하였으며, 293 K에서 추가적인 40,000 step의 평형화를 진행하였다. 모든 과정은 NVT 앙상블 조건에서 수행되었다. 각 조성에 대해 상이한 초기 조건에서 5개의 독립 시뮬레이션을 수행하고, 이를 평균하여 대표 구조 모델을 도출하였다. 시뮬레이션으로 얻어진 구조 모델을 바탕으로 부분 상관 함수 Tij(r)를 계산하였으며, 이를 합산하여 총 상관함수 Tcal(r)를 구했다. 시뮬레이션 구조 모델의 정확도는 실험으로 도출된Texp(r)과의 일치도를 평가하는 RX factor를 통해 검증되었다. RX factor는 식 11과 같이 계산되며, 1∼10 Å 범위에서 산출되었다.[40]
2.2.3 Structure analysis
Fig. 3(a)는 BLB0, BLB20, BLB35 유리의 19F MAS-NMR 스펙트럼을 나타낸 것이다. 세 유리 모두 −2.2, −3.9, −4.3 ppm 부근에서 하나의 넓은 피크를 보였으며, 반칙폭 (Full width at half maximum, FWHM)은 각각 47.9, 54.3, 54.0 ppm으로 측정되었다. BaF2 결정에서는 19F의 피크가 −15.5 ppm에서 나타나는 반면, B-F 결합이 존재하는 BO1.3F0.4 화합물에서는 −108∼-139 ppm 사이에서 피크가 형성된다고 보고된 바 있다. 또한 F 원자와 B 원자가 결합하는 B2O3-PbF2 및 B2O3-PbO-P bF 유리 시스템에서 −110 ppm 부근에서 19F 피크가 관찰된다고 알려져 있다. 본 연구에서 BB20 유리의 스펙트럼에서는 −110 ppm 부근에 명확한 피크가 관찰되지 않았으며, 이는 F 원자가 B 원자와 결합하기보다는 Ba 원자와 주로 결합하고 있음을 시사한다. La F3 결정에서는 세 가지 결정학적 위치에 대응하여 −23.5, 17.1, 24.9 ppm에서 피크가 나타난다. BLB20 및 BLB35 시료에서는 이러한 피크들이 상방 이동되었으며, FWHM은 BLB0 대비 약 6 ppm 넓어졌다. 이는 F 원자가 Ba 뿐만 아니라 La 원자와도 결합하고 있음을 나타내며, La2O3가 포함된 유리에서도 F-B 결합은 극히 드물다는 점이 추가적으로 확인되었다. X-ray 회절 측정을 통한 총 상관함수 결과에서는 BLB0 조성에서 1.33, 2.66, 4.19, 6.65 Å 부근에 네 개의 주요 피크가 관찰되었다. 첫 번째 피크는 B-O 결합에 대응되나, B와O 원자의 산란 인자가 작아 피크 위치 및 면적의 정확성은 제한적이다. 두 번째 및 세 번째 피크는 각각 Ba-O, La-O 및 Ba-Ba, Ba-La, La-La 결합쌍을 반영하는 것으로 해석된다. BaO를 BaF2로 치환해도 총 상관함수 곡선의 형태는 유지되었지만, 6.7 Å 부근의 네 번째 피크는 BaF2 농도가 증가할수록 점차 넓어지는 경향을 보였다. 이러한 차이는 Ba-Ba, Ba-La, La-La 결합의 세 번째 이웃 간 변화와 관련된 것으로 해석되나, 변화가 3차 이웃에만 국한되는 구체적인 원인은 아직 명확히 규명되지 않았다. Fig. 3(b)에서는 시뮬레이션을 통해 계산된 Tcal(r) 결과와 실험적 Texp(r) 결과가 일치함을 보였으며, BLB0, BLB20, BLB35 조성 각각에서 1∼10 Å 범위의 계산된 RX 값은 0.039, 0.031, 0.035로 나타나 구조 모델의 신뢰성을 입증하였다.[40]
Fig. 3.
(a) 19F MAS NMR spectra of BB20, BLB20, and BLB35. Reproduced from Chung et al. J. Phys. Chem B: 2023; 127:3091-9, with permission of American Chemical Society,[40] (b) Calculated and experimental total correlation functions of BLB0, BLB20, and BLB35. Reproduced from Chung et al. J. Phys. Chem B: 2023; 127:3091-9, with permission of American Chemical Society,[40] (c) Pair distribution of BLB0, BLB20, and BLB35. (d) Comparison of BO4 unit fractions (N4) obtained from MD simulation and11 B MAS NMR. Reproduced from Chung et al. J. Phys. Chem B: 2023; 127:3091-9, with permission of American Chemical Society.[40]

Fig. 3(c)는 BLB0, BLB20, BLB35 구조 모델에서의 BO3, BO4 단위 내 B-O 결합 길이 분포를 보여준다. 구조 모델을 통해 BO3 단위에서는 1.34∼1.35 Å, BO4 단위에서는 1.44 Å으로 나타났으며, 이는 MB5O8 (M= K, Rb, Cs, Tl)와 Na3ZnB5O10 결정에서 BO3 1.34∼1.39 Å, BO4 144∼1.48 Å와 비교했을 때 매우 유사하며, 시뮬레이션을 통해 확보된 구조 모델의 높은 신뢰성을 의미한다. MD 시뮬레이션의 결과는 배위수에 따른 결합 거리 변화를 잘 재현하였으며, B-O 결합의 cut-off 값은 2.0 Å으로 설정되었다. Fig. 3(d)는 구조 모델에서의 BO4 단위 비율(N4)과 11 B MAS-NMR 스펙트럼으로부터 얻은 BO4 분율을 비교한 것이다. N4값은 BaF2 함량 증가에 따라 상승하는 경향을 보였다. O 원자는 결합 환경에 따라 고립 산소(Oiso), 비가교산소(NBO), 가교산소(BO), 삼중결합 산소(Otri)로 분류하였다. BLB0 조성에서는 전체 O 원자의 약 2.5%에 해당하는 46.6개가 B 원자와 결합하지 않은 고립 산소로 존재하였다. F 원자 도입에 따라 O 원자가 치환되면서 NBO 수가 감소하고, BO 수가 증가하는 경향을 나타냈다. 전통적인 이성분계 알칼리 붕산염 유리에서는 유리 형성 산화물 농도 증가에 따라 O/B 비율이 감소하며, O/B가 약 1.75일 때 N4값이 최대에 도달한 이후 감소하고 NBO 수가 증가하는 것으로 알려져 있다. BLB0 조성의 O/B 비율은 2.3으로, N4값이 0.23에 불과하고 NBO 수가 BO 수보다 많았다. BaF2 첨가로 O/B 비율이 점차 감소함에 따라, BLB35 조성에서는 1.95까지 감소하였다. 이는 F 원자가 B원자와 직접 결합하지 않고, O 원자만 치환함을 나타내며, 결과적으로 N4값의 증가로 이어졌다. 여기서 주목할 점은 F 원자와 NBO의 총합이 BLB0, BLB20, BLB35 조성 모두에서 각각 1003.0, 1007.6, 1013.8로 거의 일정하게 유지되었다는 사실이다. 이는 F 원자와 NBO가 Ba 및 La 양이온에 대한 전하 보상 역할을 함께 수행함을 의미한다.[40] 시뮬레이션 모델에서는 Ba 및 La 원자 수가 고정되어 있으며, O와 F 원자의 전하는 각각 −1.35 및 −0.8로 설정되어 있다. 따라서 F 원자와 NBO는 거의 동일한 전하 보상 능력을 갖는 것으로 간주할 수 있으며, 이는 비록 고정 전하 모델을 사용했더라도 실제 유리 네트워크 내 전하 분포가 국지적으로 민감하게 조정될 수 있음을 시사한다. Ba-O, Ba-F, La-O, La-F의 부분 상관함수 결과를 통해 Ba-(O, F) 결합의 cutoff 거리는 3.6 Å, La-(O, F) 결합의 cutoff 거리는 3.1 Å으로 설정되었으며, Ba 및 La 주변의 O 원자 및 F 원자의 평균수를 Table 4에 요약하였다. F 원자 도입이 증가할수록 Ba 및 La 주변의 비가교 산소 수는 감소하고, 가교 산소 수는 소폭 증가하는 경향을 보였다. 특히 La 주변의 F 및 NBO 합계가 Ba 주변보다 약 10% 많았으며, 이는 La 양이온이 Ba 양이온에 비해 더 높은 전하를 가지기 때문으로 해석된다. 또한 La-O, La-F 결합 거리가 Ba-O, Ba-F 결합보다 짧아 더 효율적인 전하 보상이 가능함이 구조적으로 확인되었다.[40]
Table 4.
Distribution of O and F Atoms coordinated to La and Ba. Reproduced from Chung et al. J. Phys. Chem B:2023;127:3091-9, with permission of American Chemical Society.[40]
B 원자 주변의 O-B-O Bond angle distribution (BAD)은 유리 네트워크의 기본적인 구조 특성을 반영한 다. Fig. 4(a)는 BLB0 유리에서 BO3 및 BO4 단위의 O-B-O BAD를 각각 분리하여 나타낸 것으로, BO3 단위는 약 120º를 중심으로, BO4 단위는 약 109º 중심으로 분포하여 이상적인 삼각형 및 사면체 구조를 잘 반영한다. 전체 BAD는 이 두 단위의 조합으로 설명할 수 있으며, 이를 통해 유리 매트릭스 내 B-O 네트워크가 BO3와 BO4 단위로 구성되어 있음을 확인할 수 있다. BLB0, BLB20, BLB35 조성 간의 비교를 위해 Fig. 4(b)에 O-B-O Bond angle distribution을 나타냈다. Ba F2 함량 증가에 따라 전체 BAD는 더 낮은 각도(109º 부근)으로 치우치는 경향을 보였으며, 이는 BO4 단위의 증가를 시사한다. 한편, Ba 및 La 이온은 직접적으로 B-O-B 네트워크에 속해 있지는 않지만, 주변 O 원자와 결합하여 네트워크 간 간극을 조절하는 역할을 한다. Fig. 4(c)와 (d)는 각각 O-Ba-O, O-La-O Bond angle distribution을 나타내며, Ba 주변의 O 원자 간 결합각은 주로 40º∼60º 범위에 하는 것을 확인했으며, La 주변의 O 원자 간 결합각은 60º∼80º 범위에 분포하는 것을 확인할 수 있었다.
Fig. 4.
(a) Decomposed plots of bond angle distribution of O-B-O angle in BLB0. (b) Bond angle distribution of O-B-O angle of BB20, BLB20, and BLB35. (c) Bond angle distribution of O-Ba-O angle of BB20, BLB20, and BLB35. (d) Bond angle distribution of O-La-O angle of BB20, BLB20, and BLB35.

본 연구 사례에서는 Ba O를 BaF2로 치환함에 따라 BO 비율이 증가하고 O/B 비율이 감소하는 경향이 관찰되었다. 또한, F 원자가 Ba 및 La 양이온의 전하를 보상함에 따라, 이에 대응하는 비가교 산소 수도 감소하는 것으로 나타났다. 이를 보다 상세히 분석하기 위해 가교산소 수에 따른 Qn 구조 단위 비율을 Fig. 5에 나타냈다. BLB0 조성에서는 BO3 단위 내 Q2 단위, BO4 단위 내 Q3 단위의 비율이 가장 높게 나타났으며, 대부분의 BO3 및 BO4 단위는 하나의 NBO를 가지고 있었다. BLB20 조성에서는 BO3 내 Q3 및 BO4 단위 내 Q4 단위의 비율이 증가하였으며, 특히 NBO가 없는 Q3 및 Q4 비율이 가장 높게 나타났다. 이는 BaF2 첨가가 BO3 및 BO4 단위 내 NBO 수를 감소시키는 효과를 가져온다는 것을 의미한다. 실제 BO3 단위에서의 평균 NBO 수는 BLB0, BLB20, BLB35 조성에서 각각 1.31, 0.92, 0.62개였으며, BO4 단위에서는 각각 1.08, 0.72, 0.48개로 나타났다. 시뮬레이션 구조 모델에서는 BO4 단위가 BO3 단위보다 NBO 수가 항상 적었는데, 이는 BO4 단위가 음전하를 띄고 있어 추가적인 NBO 형성이 불리하기 때문으로 해석된다. Fig. 6는 BLB0, BLB20, BLB35 유리의 Ring size distribution을 나타낸다. BaF2 함량이 증가할수록 ring size가 줄어드는 결과를 얻었는데, 이는 O가 F로 치환됨에 따라 유리 네트워크의 연결성(network connectivity)이 감소했음을 의미한다.
Fig. 5.
(a) Qn Distribution of BLBO glass associated with BO3 and BO4 units. (b) Qn Distribution of BLB2O glass associated with BO3 and BO4 units. (c) Qn Distribution of BLB35 glass associated with BO3 and BO4 units. (d) Compositional trend of Qn Distribution contributions from BO3 and BO4 units in BLB glasses.

Ba와 La은 첫 번째 배위권에서 주로 O 및 F 원자와 결합하며, 많은 O 원자들은 B 원자와 결합하여 B-O 네트워크를 형성한다. 따라서 Ba와 La는 B-O 네트워크에 인접해 존재하며, 이들 양이온의 배위 구조에 따라 인접 네트워크 간 간격이 조정된다. Fig. 7은 BLB0, BLB20, BLB35 각각의 시뮬레이션 모델에서 BO3, BO4 단위 및 Ba-O, Ba-F, La-O, La-F 결합 구조를 나타낸다. BLB0에서는 B-O 네트워크가 모델 전체에 일정한 간격을 유지하며 고르게 퍼져 있으며, Ba 및 La 양이온은 가교 산소 또는 비가교 산소를 통해 B-O 네트워크와 연결되어 있다. 전체 산소의 약 2.5%를 차지하는 고립 산소는 B 원자와 결합하지 않고 Ba 또는 La에만 결합하여 네트워크 간 간극을 확장시키고 국지적인 클러스터를 형성하였다. Oiso 주변의 원자 배열을 보면, 4개의 Ba와 1개의 La가 하나의 Oiso를 중심으로 위치해 있으며, 이들은 주변 B− O 네트워크의 BO 또는 NBO와도 결합되어 약간의 네트워크 확장을 유발하고 있음을 알 수 있다. BLB20 모델에서는 160개의 산소가 320개의 F 원자로 치환되었으며, 전체 음이온의 약 16%를 차지하는 F 원자들은 주로 Ba 및 La 배위 다면체에 포함되어 있었다. F 원자들이 주로 B−O 네트워크 간의 간극에 분포하는 모습을 보여주었으며, 평균적으로 Ba 및 La 각각에 약 2개의 F 원자가 결합되어 있었다. 이는 BLB0의 Oiso와 유사한 방식으로 작용하며, 네트워크 간의 간격을 더욱 증가시키는 결과를 가져왔다. BLB35 모델에서는 F 원자 밀도가 더욱 증가하여, 하나의 배위 다면체에 최대 3개의 F 원자가 결합된 경우도 관찰되었다. 특히 F 원자가 다수 존재하는 Ba 및 La 배위 다면체가 형성되면서, B− O 네트워크 간의 간극이 더욱 넓어졌으며, F 원자 주변에 Ba− F− Ba 또는 La− F− La 구조가 형성되어 B− O 네트워크 사이의 간극을 채우는 모습을 확인할 수 있었다.
Fig. 7.
(a) Structure model of BLB0 showing BO3 and BO4 units along with Ba-O, Ba-F, La-O and La-F bonding environments. (b) Structure model of BLB20 showing BO3 and BO4 units along with Ba-O, Ba-F, La-O and La-F bonding environments. (c) Structure model of BLB35 showing BO3 and BO4 units along with Ba-O, Ba-F, La-O and La-F bonding environments.

결론
유리 구조에 대한 정밀한 이해는 고기능성 유리 소재의 설계 및 응용에 필수적인 요소이다. 그러나 유리는 결정질 재료와 달리 장거리 질서를 가지지 않는 비정질 상태로 존재하며, 단순한 실험적 분석이나 이론적 예측 모델만으로는 전체 구조를 정확히 규명하는 데 한계가 따른다. 특히 다성분계 유리 시스템의 경우, 조성의 복잡성과 원자 간 상호작용의 다양성으로 인해, 실험적 기법만으로는 short-range order와 medium-range order에 대한 정밀한 해석이 어렵다. 이러한 배경 속에서 분자동역학 (Molecular dynamics, MD) 시뮬레이션은 유리 구조 해석의 중요한 도구로 부각되고 있다. MD 시뮬레이션은 실험적 분석의 한계를 보완하고, 구조와 물성 간의 관계를 정량적으로 제시할 수 있는 강력한 수단을 제공한다. 본 리뷰에서는 MD 시뮬레이션의 핵심 구성 요소인 포텐셜 함수, 수치 적분 알고리즘, 열역학 앙상블 조건에 대해 체계적으로 정리하고, 유리 시스템에 적합한 시뮬레이션 기법의 이론적 기반을 고찰하였다. 포텐셜 함수는 원자 간의 결합 특성과 화학적 환경을 수치적으로 표현하는 핵심 인자로서, 단순한 2체 모델부터 전하 평형 기반의 고차원 포텐셜까지 그 특성과 적용 가능 범위를 상세히 논의하였다. 또한 Verlet 계열 적분기법을 중심으로 알고리즘 선택의 기준과 에너지 및 운동량 보존 특성, 계산 안정성 측면에서의 장단점을 비교하고, NVT 및 NPT 등 다양한 앙상블 조건의 실제적 적용 방식과 그 물리적 의미를 명확히 제시하였다. 이러한 이론적 토대를 바탕으로 다성분계 glass system의 구조 분석 사례로서, Ba F2가 첨가된, Ba O-La2O3-B2O3계 다성분 fluoroborate 유리계에 대한 분자동역학 시뮬레이션 기반 구조 모델을 리뷰하였다. F 원자의 도입에 따라 BO4 단위 비율(N4)의 증가, 비가교 산소(NBO)의 감소, Qn 구조 단위의 재배열, 네트워크 링 사이즈 감소 등 다양한 구조적 변화 양상이 실험 결과와 일치하도록 재현되었으며, F 원자가 O 원자를 치환하면서도 유사한 전하 보상 효과를 갖는다는 점을 구조적으로 설명하였다. 또한, La3⁺ 및 Ba2⁺ 이온의 배위 환경 분석을 통해, 이들 이온이 유리 네트워크 간 간극을 조절하고 전하 보상에 기여한다는 구조적 해석을 제시하였다. 결론적으로, 본 연구는 MD 시뮬레이션이 유리 구조의 복잡한 미시구조를 예측하고 해석하는 데 있어 매우 유효한 방법론임을 입증하였다. 향후 포텐셜 모델의 고도화, 시뮬레이션 스케일의 확대, 실험–계산 융합 분석 기법의 발전을 통해, MD 시뮬레이션 기반 유리 구조 분석은 더욱 정밀하고 신뢰성 높은 예측 도구가 될 것으로 기대된다.