The Impact of Descriptor Characteristics on the Accuracy of Neural Network Potentials for Predicting Material Properties

Descriptor 특성이 신경망포텐셜의 소재 물성 예측 정확도에 미치는 영향에 관한 연구

Abstract

In this study, we aim to derive the descriptor vector conditions that can simultaneously achieve the efficiency and accuracy of artificial Neural Network Potentials (NNP). The material system selected is silicon, a highly applicable material in various industries. Atomic structure-dependent energy data for training artificial neural networks were generated through density functional theory calculations. Behler-Parrinello type atomic-centered symmetric functions were employed as descriptors, and various length vector NNPs were generated. These NNPs were applied to reproduce the structure and mechanical properties of silicon materials in molecular dynamics simulations. In our findings, the minimum vector length for achieving both learning and computational efficiency while maintaining property reproducibility is approximately 50. It was also observed that, for the same conditions, incorporating more angle-dependent symmetric functions into the descriptor vector, could enhance the accuracy of NNP. Our results can provide guidelines for optimizing the conditions of descriptor vectors to achieve both efficiency and accuracy of NNP, simultaneously.

본 연구에서는 신경망포텐셜(Neural Network Potential)의 효율성과 정확도를 동시에 달성할 수 있는 기술자 벡터 조건을 도출하고자 한다. 소재 시스템은 단원소 소재인 실리콘으로 선정하였으며, 인공신경망 학습을 위한 원자 구조별 에너지 데이터는 밀도범함수이론 계산을 통하여 생성하였다. Behler-Parrinello 타입의 원자중심대칭함수를 기술자 벡터로 사용하였고, 다양한 벡터 길이에 대한 신경망포텐셜 생성 후 분자동역학 시뮬레이션에 적용하여 실리콘 소재의 구조 및 기계적 물성 재현성을 평가하였다. 실험 결과, 물성 재현 정확도를 유지하면서 학습 및 계산 효율성을 동시에 달성할 수 있는 기술자벡터의 최소 길이는 약 50이고, 소재의 기계적 물성이 이 길이에 더 큰 영향을 받으며, 같은 길이의 조건에서는 방사 대비 각도 방향 대칭함수를 더 반영하면 신경망포텐셜의 정확도가 올라감을 발견하였다. 이를 토대로 신경망포텐셜의 효율성과 정확도 동시 달성을 위한 최적의 기술자벡터 설정 가이드라인 제공이 가능할 것으로 기대된다.

1. 서론

최근 컴퓨터 기술의 급속한 발전과 더불어 소재 연구에 대한 원자단위 시뮬레이션(Atomistic simulation) 기법의 활용 또한 확대되고 있다. 특히 분자동역학(Molecular dynamics) 시뮬레이션은 계를 구성하는 모든 원자들의 시간에 따른 거동을 탐색할 수 있는 강력한 도구로써 다양한 유무기 소재의 물성 예측 및 화학 반응 메커니즘 연구 등에 널리 활용되고 있다[1,2]. 분자동역학 수행 중 각 원자에 작용하는 힘은 해석식으로 구성된 원자간 포텐셜에 의해 계산된다. 분자동역학 시뮬레이션의 정확도는 원자간 포텐셜에 의해 정해지므로, 주어진 물질 군을 정확히 기술할 수 있는 포텐셜 개발이 지난 수 십 년간 이루어져 왔다[3]. 이러한 포텐셜 개발은 개발 과정에 대한 연구자의 오랜 경험뿐만 아니라 해당 물질계에 대한 깊은 이해를 필요로 하므로, 분자동역학 적용 범위 확대에 많은 제약을 가해왔다.

이러한 문제점을 해결할 수 있는 수단으로써 머신러닝 기법을 적용한 포텐셜 개발이 최근 주목을 받고 있다 [4,5]. 특히 인공신경망을 활용한 포텐셜 모델은 양자계산데이터(ab initio data)에 기반을 둔 학습을 통하여 포텐셜 파라미터를 도출하므로, 낮은 개발 난이도 및 높은 정확성으로 인해 그 활용도가 매우 높다. 신경망포텐셜(Neural Network Potential-이후 NNP)의 효율성과 정확도는 원자 주변의 다양한 화학 환경을 반영할 수 있는 입력층의 기술자(Descriptor) 벡터에 크게 의존한다. 시스템 내 각 원자들은 다양한 환경에 놓이게 되므로, 시스템의 모든 정보를 충분히 반영하기 위해서는 시스템의 자유도 대비 더 많은 수의 기술자벡터 요소가 필요하다. 반면, 많은 수의 기술자벡터 요소는 입력층의 노드 수 증가를 의미하므로 인공신경망의 학습 효율을 떨어뜨리며, 분자동역학 계산 시간 증가를 초래한다. 따라서 본 논문에서는 기술자벡터 종류 및 길이 변화에 따른 NNP 성능 변화에 대한 분석을 통하여, 효율성과 정확도를 동시에 달성할 수 있는 기술자벡터의 조건을 도출하기 위한 연구를 진행하였다.

2. 이론적 배경

2.1 원자중심대칭함수

원자가 포텐셜에 의해 도출되는 총 에너지는 전체계의 병진 및 회전, 그리고 같은 화학종 원자간 치환에 대하여 변하지 않고 일정하게 유지되어야 한다. 따라서 기술자벡터 또한 이러한 연산에 대한 불변성(invariant)을 가지고 있어야 하기 때문에, 시스템 내 각 원자의 좌표를 기술자벡터에 직접 적용하는 것은 부적절하다. 따라서 본 연구에서는 Behler-Parrinello 타입의 원자중심대칭함수[6]을 기술자벡터로 사용하였다. 적용한 대칭함수의 방사(Radial) 성분 및 각도(Angular) 성분은 식(1), 식(2)로 각각 표현된다.

\(\begin{align}G_{i}^{2}=\sum_{j} e^{-\eta\left(R_{i j}-R_{s}\right)^{2}} \cdot f_{c}\left(R_{i j}\right)\end{align}\)      (1)

\(\begin{align}\begin{aligned} G_{i}^{4}= & 2^{1-\zeta} \sum_{j, k \neq i}\left(1+\lambda \cos \theta_{i j k}\right)^{\zeta} \cdot e^{-\eta\left(R_{i j}^{2}+R_{i k}^{2}+R_{j k}^{2}\right)} \\ & \cdot f_{c}\left(R_{i j}\right) \cdot f_{c}\left(R_{i k}\right) \cdot f_{c}\left(R_{j k}\right)\end{aligned}\end{align}\)       (2)

식(1), 식(2)에서 i는 중심 원자의 인덱스, j와 k는 이웃 원자들의 인덱스, 그리고 R_ij, R_ik, R_jk는 이들 원자간 거리이다. 또한 θ_ijk는 원자 i에서 원자 j로 이어지는 벡터와 원자 i에서 원자 k로 이어지는 벡터 사이의 각도를 의미한다. 가우시안 함수의 폭과 중심은 η와 R_s에 의하여 각각 결정되며, 각도 함수의 개형은 ζ와 λ에 따라 바뀐다. 함수 f_c는 원자간 거리가 R_c에 접근함에 따라 대칭함수가 식(3)에 의하여 점차 0으로 접근할 수 있도록 해준다.

\(\begin{align}f_{c}\left(R_{i j}\right)=\left\{\begin{array}{ll}0.5\left[\cos \left(\frac{\pi R_{i j}}{R_{c}}\right)+1\right] & \left(R_{i j} \leq R_{c}\right) \\ 0 & \left(R_{i j}>R_{c}\right)\end{array}\right.\end{align}\)       (3)

2.2 고차원신경망

NNP는 입력층에 기술자벡터의 형태로 주어지는 시스템 정보를 반영하여 에너지를 도출한다. 이때 시스템의 자유도는 입력층의 노드 수, 즉 기술자벡터의 길이에 의해 결정되는데, 시스템 내 원자의 수가 증가할수록 자유도는 증가하므로 기술자벡터의 길이는 시스템의 크기에 비례한다. 따라서 일반적인 feed-forward 인경신경망으로 포텐셜을 구성하면 한번 결정된 기술자벡터의 길이를 변경하는 것이 불가능하므로, 특정 시스템 크기에 대하여 학습된 NNP는 시스템 크기가 달라질 경우에 적용이 불가능해지는 치명적인 한계가 존재한다.

이러한 한계를 극복하기 위해, 원자 신경망으로 구성된 고차원 신경망을 적용하기도 한다 [7,8]. 고차원신경망에서는 각 원자 주변 상황을 반영한 원자에너지(E_i)가 원자신경망에 의하여 도출되며, 시스템의 총 에너지(E)는 식(4)에 의해 원자에너지의 총 합으로 나타내어진다.

\(\begin{align}E=\sum_{i} E_{i}\end{align}\)       (4)

총 N개의 원자가 포함된 시스템에 대한 고차원신경망의 구성을 그림 1에 도식으로 나타내었다. 이러한 고차원신경망으로 구성된 포텐셜에서 원자에 작용하는 힘은 식(5)에 의해 해당 원자의 위치에 대하여 총 에너지를 편미분 함으로써 계산된다.

그림 1. (a) N 개의 원자가 포함된 시스템에 대한 고차원신경망 포텐셜의 구성 및 (b) 고차원 신경망포텐셜 내 원자신경망의 구조

Fig. 1. (a) Structure of a high-dimensional neural network potential for a system with N atoms and (b) schematic diagram of atomic neural network in the high-dimensional neural network potential

\(\begin{align}F_{k, \alpha}=-\frac{\partial E}{\partial R_{k, \alpha}}=-\sum_{i=1}^{N} \sum_{j=1}^{M} \frac{\partial E_{i}}{\partial G_{i, j}} \frac{\partial G_{i, j}}{\partial R_{k, \alpha}}\end{align}\)       (5)

식(5)에서 F_k,α와 R_k,α는 원자 k에 대한 α-성분(α = x,y,z)의 힘과 좌표를 의미하며, M은 기술자벡터의 길이를 나타낸다

3. 개발 및 분석 방법

3.1 기술자벡터 설정

기술자벡터 길이에 따른 분자동역학 정확도 분석을 위하여 본 연구에서는 다양한 길이의 기술자벡터를 설정하였다. 식(1) 및 식(2)로 표현되는 대칭함수 내 파라미터 중, R_s는 중심 원자의 spherical shell에 대한 분석을 제외한 상황에서는 일반적으로 0으로 고정되며, λ 또한 +1 혹은 –1로 고정된다. 따라서 기술자벡터의 길이는 방사 방향의 정보를 반영하는 η와 각도 방향의 정보를 반영하는 ζ의 설정수치 개수에 따라 변하게 된다.

본 연구에서는 기술자벡터의 길이를 16, 28, 40, 52, 64로 변화시키면서 NNP의 성능을 분석하였다. 방사 방향 및 각도 방향 정보에 대한 상대적 중요성의 분석하기 위하여, η의 설정수치 개수 조절을 통하여 기술자벡터의 길이를 조정한 경우 (표 1) 와 ζ의 설정 수치 개수 조절을 통하여 기술자벡터의 길이를 조정한 경우(표2)로 나누어 비교 분석하였다.

표 1. 기술자벡터 길이에 따른 대칭함수 내 η 설정수치 개수 변화

Table 1. The number of η parameter values in symmetry functions for various descriptor lengths

표 2. 기술자벡터 길이에 따른 대칭함수 내 ζ 설정수치 개수 변화

Table 2. The number of ζ parameter values in symmetry functions for various descriptor lengths

3.2 데이터 생성 및 학습

단원소 소재 중 정보전자 및 화학산업에서 응용도가 높은 실리콘을 본 연구의 소재 시스템으로 선정하였다. 인공신경망 학습을 위한 데이터는 Quantum Espresso package [9]를 이용한 밀도범함수이론(Density functional theory-이후 DFT) 계산으로 생성하였다. 교환-상관성 범함수는 GGA-PBE [10], 역격자 적분을 위하여 1 × 1 × 1 k-point mesh를 각각 적용하였다.

결정질과 비정질 고체 실리콘 (그림 2), 그리고 액체 상태의 실리콘 대상으로 학습 데이터를 생성하였다. 결정질 고체 실리콘과 액체 실리콘은 각각 500 K과 2000 K에서, 그리고 비정질 고체 실리콘은 2000 K에서 300 K 까지 70 K/ps 의 속도로 냉각시키면서 ab-initio 분자동역학을 수행하였으며 20 fs 간격으로 시스템의 원자 구조 및 에너지를 추출하여 90%는 인공신경망의 학습 데이터로 사용하였고 나머지 10%는 학습 완료된 NNP의 검증 과정에 활용하였다. 원자신경망은 두 개의 은닉층를 배치하였으며, 각 은닉층은 30개의 노드를 설정하였다. SIMPLE-NN package [11] 를 통하여 신경망 학습, 검증 및 포텐셜 생성을 진행하였으며, 최적화 과정에서 에너지와 힘에 대한 에러를 함께 고려하였다.

그림 2. (a) 결정질 및 (b) 비정질 고체 실리콘 원자 구조

Fig. 2. Atomic structure of (a) crystalline and (b) amorphous solid silicon

3.3 분자동역학 시뮬레이션

학습이 완료된 NNP 적용하여 분자동역학을 수행 후 그 정확도를 분석하였다. 본 연구에서 수행된 모든 분자동역학 시뮬레이션은 LAMMPS package [12]를 이용하였으며, 시스템 온도와 압력은 NPT 앙상블에서 300K 밎 1 atm으로 유지되었다. 각 시뮬레이션은 1 ft의 시간 간격으로 1 ns동안 적분 되었으며, 시뮬레이션 결과로부터 시스템 물성을 계산하였다. 구조 물성 예측 성능 분석을 위하여 방사분포함수(Radial distribution function)을 측정하였으며, 기계적 물성 예측 성능 분석을 위하여 부피탄성계수(Bulk modulus)를 측정하였다.

4. 결과 및 고찰

임의의 원자 구조에 대한 에너지 및 힘을 학습 완료된 NNP로 계산하고, 이를 DFT로 계산한 결과와 비교를 통하여 성능을 평가하였다. 기술자벡터의 길이가 64인 경우와 η의 설정 개수 변화를 통해 기술자벡터의 길이를 16으로 줄인 경우, η 설정 개수 변화를 통해 기술자벡터의 길이를 16으로 줄인 경우에 NNP와 DFT 계산 결과 사이의 선형 관계도로 나타내었다(그림 3). 기술자벡터의 길이가 64일 경우 에너지와 힘에 대한 계산 모두에서 NNP와 DFT 결과 사이에 뛰어난 선형 관계를 보여주고 있다 (그림 3a). 반면, 기술자벡터의 길이가 16으로 줄어들게 되면 64길이 대비 선형 관계에서 벗어나는 경우가 증가하는 것을 확인할 수 있다 (그림 3b,c).

그림 3. NNP 및 DFT 통한 에너지(E) 및 힘(F) 계산 비교. NNP는 기술자벡터의 길이를 (a) 64로 설정한 경우, (b) η개수 변화를 통하여 16으로 설정한 경우, (c) ζ 개수 변화를 통하여 16으로 설정한 경우에 대하여 계산

Fig. 3. Comparison of NNP and DFT calculations for energy(E) and force(F). For NNP calculations, descriptor vector lengths are (a) 64, (b) 16 by reducing the number of η values, and (c) 16 by reducing the number of ζ values

기술자벡터 길이에 따른 NNP 와 DFT 계산 결과 차이를 정량적으로 측정하기 위하여, 에너지와 힘에 대한 NNP 와 DFT 계산 결과의 평균적인 차이를 식(6)과 식(7)로 각각 정의하였다.

\(\begin{align}\Delta E_{\text {arg }}=\frac{1}{N_{s}} \sum_{i=1}^{N_{s}}\left(E_{i}^{\mathrm{NNP}}-E_{i}^{\mathrm{DFT}}\right)^{2}\end{align}\)       (6)

\(\begin{align}\Delta F_{\mathrm{avg}}=\frac{1}{N_{s}} \sum_{i=1}^{N_{s}}\left(F_{i}^{\mathrm{NNP}}-F_{i}^{\mathrm{DFT}}\right)^{2}\end{align}\)       (7)

위 식에서 E^α_i와 F^α_i(α = NNP, DFT)는 α를 통하여 계산된 원자 구조 i에 대한 에너지와 힘, N_s는 계산에 사용된 원자 구조의 개수를 의미한다. 그림 4는 기술자벡터 길이에 따른 ∆E_avg와 ∆F_avg의 계산 결과이다.

그림 4. 기술자벡터 길이에 따른 (a) 에너지와 (b) 힘에 대한 NNP와 DFT의 계산 결과 차이

Fig. 4. Difference of NNP and DFT calculation results of (a) energy and (b) force based on descriptor length

기술자벡터 길이가 줄어들수록 ∆E_avg와 ∆F_avg는 증가하는 것을 확인할 수 있다. 구체적으로, 기술자벡터의 길이가 64에서 52로 줄었을 경우 큰 변화를 관찰할 수 없었지만 40으로 줄었을 경우 NNP와 DFT 계산 결과 차이가 눈에 띄기 시작하며, 28에서 16으로 줄었을 경우 ∆E_avg와 ∆F_avg 모두 급격히 증가하는 것을 관찰하였다. 또한, 기술자벡터 길이가 같은 경우에 ζ 개수가 적은 경우가 η 개수가 적은 경우 대비 ∆E_avg와 ∆F_avg 모두 더 큰 폭으로 증가시키는 것을 확인하였다. 이것은 기술자벡터의 길이가 같다면, 방사 방향의 대칭함수보다 각도 방향의 대칭함수를 기술자벡터에 더 많이 반영하는 것이 NNP 정확도 증가에 유리하다는 의미이다.

에너지와 힘 계산에 대한 NNP의 정확도가 분자동역학 시뮬레이션의 소재 물성 예측 성능에 미치는 영향을 확인하기 위하여 방사분포함수 및 부피탄성계수를 분자동역학 시뮬레이션 결과로부터 도출하였다.

그림 5는 결정질 및 비정질 고체 실리콘에 대한 방사분포함수 산출 결과이다. ∆E_avg및 ∆F_avg와는 달리, NNP의 기술자벡터 길이에 따른 방사분포함수 산출 결과는 거의 변화가 없는 것을 확인할 수 있다. 비정질 실리콘의 방사분포함수의 경우 거의 동일하게 유지되었으며, 비록 기술자벡터 길이가 가장 짧은 경우(즉, 길이가 16인 경우)는 조금의 차이를 보였지만 그 변화폭은 크지 않았다 (그림 5a). 특히 결정질 실리콘의 방사분포함수는 모든 기술자벡터 길이에서 동일한 결과를 보였다(그림 5b). 즉, 방사분포함수로 대표되는 소재의 구조 물성은 NNP의 에너지 및 힘 계산 정확도에 크게 영향을 받지 않는다는 것을 알 수 있다.

그림 5. 기술자벡터 길이에 따른 (a) 비정질 및 (b) 결정 질 실리콘에 대한 방사분포함수

Fig. 5. Radial distribution functions of (a) amorphous and (b) crystalline silicon for various descriptor lengths

반면, 분자동역학을 통한 소재의 기계적 물성 산출은 NNP의 정확도에 크게 영향을 받는 것으로 확인 된다. 그림 6은 결정질 실리콘의 부피탄성계수를 다양한 길이의 기술자벡터를 적용하여 산출한 결과이다. 산출된 부피탄성계수는 i) 기술자벡터 길이가 64에서 52로 줄었을 경우 변화가 거의 없고, ii) 길이가 40 이하에서는 길이가 짧아짐에 따라 기술자벡터의 변화 폭이 급격히 증가하며, iii) 같은 길이에서 ζ 개수가 적은 경우가 η 개수가 적은 경우 대비 변화 폭이 크다는 측면에서 ∆E_avg및 ∆F_avg의 변화 패턴과 정확히 일치한다. 따라서 NNP의 에너지 및 힘 계산 정확도는 분자 동역학을 통한 소재의 기계적물성 산출에 직접적인 영향을 주는 것을 확인할 수 있었다.

그림 6. 기술자벡터 길이에 따른 결정질 실리콘에 대한 부피탄성계수

Fig. 6. Bulk modulus of crystalline silicon as a function of descriptor length

5. 결론

적절한 조건의 기술자벡터를 설정하는 것은 NNP의 학습 효율성 및 시뮬레이션 정확도를 동시에 달성하기 위한 필수 조건이다. 본 논문에서는 다양한 조건의 기술자벡터를 적용하여 NNP를 개발하였고, 이를 적용한 분자동역학 시뮬레이션을 실리콘 소재에 적용하여 다음과 같은 결과를 얻었다.

1. 단원소 소재시스템에 NNP를 적용할 경우, 기술자벡터 길이가 약 50 이상에서는 NNP의 정확도 차이를 거의 관찰할 수 없었다. 따라서 정확도를 유지하면서, 학습 및 계산 효율성을 동시에 달성할 수 있는 벡터의 최소 길이는 약 50인 것으로 나타났다.

2. 기술자벡터의 길이가 50보다 짧아지면 NNP의 정확도가 감소하기 시작하며, 특히 28보다 짧아질 경우 정확도가 급격히 감소함을 확인하였다. 따라서 효율성 향상을 위해 길이를 50 이하로 설정하더라도 최소한의 정확도 확보를 위해서는 약 30 이상으로는 유지하는 것이 필요할 것으로 판단된다.

3. 같은 길이의 기술자벡터 조건에서는 방사 방향의 대칭함수보다 각도 방향의 대칭함수를 기술자벡터에 더 많이 반영하는 것이 NNP의 정확도를 높일 수 있는 것으로 확인되었다.

4. 분자동역학 시뮬레이션을 통하여 소재의 물성을 도출할 경우, NNP 정확도는 소재의 기계적 물성에 직접적인 영향을 미치지만 구조 물성에는 큰 영향을 미치지 않는다. 따라서 소재의 기계적 물성 예측을 위한 분자동역학 수행의 경우, 충분한 길이의 기술자벡터를 적용하여 NNP의 정확도를 높일 필요가 있다. 본 연구에서 분석한 결과를 적용하여 상황별 맞춤형 기술자벡터 조건을 설정한다면, 효율성과 정확도를 동시에 갖춘 NNP를 개발할 수 있을 것으로 기대된다.