1. 서 론

분자동역학은 전산모사 기술의 하나로써 원자를 전하, 부피 및 질량을 가지고 있는 점으로 간주하고 이들 사 이에 화학적인 결합이 스프링처럼 연결이 되어있다고 가 정한 후에 이러한 점들이 역장(力場), 즉 force-field 안 에서 움직이는 것을 계산하여 원자의 움직임을 모사하는 방법이다. 이를 통하여, 원자 및 분자의 운동과 관련된 diffusion, configuration, morphology 및 crystallization 등의 다양한 특성들을 계산할 수 있다[1-3]. 이러한 특 징으로 인하여, 고분자 분리막의 주요 성능 지표인 diffusivity, ion conductivity, solubility 등을 계산하는데 활 용될 수 있고[4-6], 실험만으로는 정확히 파악하기 힘든 주쇄의 conformation이나 configuration 같은 구조 특성 [7,8], 그리고 고분자 분리막 내부에 형성되는 free volume, water channel과 같은 형태적 특성[9-12] 등을 분 석할 수 있기 때문에 점점 더 많은 연구결과 들이 발표 되고 있다. 따라서 이러한 유용한 기술인 분자동역학의 연구결과들을 실험에 적용하고자 하는 경우나, 혹은 직 접 분자동역학 연구를 진행해 보고자 하는 연구자들에 게 있어서 분자동역학의 근간이 되는 force-field에 대한 이해는 반드시 선행되어야 한다.

앞서 언급하였듯이, 분자동역학에서는 force-field를 바 탕으로 원자의 움직임을 결정하기 때문에, 이러한 움직 임에 영향을 주는 여러 가지 요소들이 미리 정의가 되어 있어야 한다[13]. 일반적으로 크게 (1) 화학결합을 매개 로 한 원자들 간의 상호작용으로 표현되는 bond term 과 (2) 서로 결합하지 않은 원자들이 2차 간력을 통하 여 상호작용하는 non-bond term으로 구분되며, 다음 식 (1)에 제시된 고분자에서 널리 사용되고 있는 대표적인 force-field 중의 하나인 COMPASS (condensed-phase optimized molecular potentials for atomistic simulation studies)를 통하여 이를 살펴볼 수 있다[13-15]. 이러한 force-field의 존재는 분자동역학의 매우 중요한 특징으 로써, 일일이 전자의 상태를 풀어서 계산해 나가는 양 자역학과 달리 미리 입력된 변수를 이용해서 운동방정 식만 계산하면 되므로, 최소계산 단위가 양자역학 보다 더 큰 원자라는 것을 감안하더라도 분자동역학의 계산 속도가 월등히 빨라지는 장점으로 작용한다. 반면에 force-field가 정의되어 있지 않은 조건에서는 분자동역 학 계산이 불가능하다는 심각한 단점이 되기도 하는데, 전체 시스템에서 일부만 이런 경우가 발생해도 전체 계 산결과의 신뢰도가 급격하게 하락하거나 아예 계산이 되지 않는 결과를 초래하기도 한다.

(1)

본 연구에서는 이러한 force-field의 종류가 고분자 분리막 내부에 존재하는 기체 분자의 확산 현상에 미치 는 영향을 분석하고자 하였다. 특히, 고분자 분리막 소 재로서 상용화되어 널리 이용되고 있는 폴리이미드인 Kapton^®을 이용하여, force-field에 정의되어 있지 않은 원자가 존재하여 구조가 깨지거나 계산이 실패하는 상 황을 배제하고 장기간 분자동역학 계산을 수행하여, 각 각의 force-field가 분자동역학 계산에 미치는 영향을 좀 더 자세히 살펴보고자 하였다.

2. 전산모사

2.1. 기체 투과 전산모사를 위한 Kapton^® 모델 제작

분자동역학 계산을 위하여 사용한 전산모사 프로그램 은 Materials Studio (Dassault Systemes, BIOVIA Corp., USA)로써 아래에 언급되는 모든 모듈들은 동 프로그램 내에 이식되어 있는 것을 사용하였다. Kapton^®의 화학 구조를 바탕으로 Visualizer 모듈을 이용하여 repeat unit 모델을 생성하였고, 이렇게 생성된 repeat unit은 분자동 역학 모듈인 Forcite를 이용하여 Geometry optimization 을 수행하였다. 이때 force-field가 고분자 구조에 미치 는 영향을 파악하기 위하여, Forcite 모듈에서 선택 가 능한 COMPASS, pcff, Dreiding, Universal, cvff forcefield[ 13]들을 repeat unit에 적용하여 Geometry optimization을 수행하고 각각의 구조를 비교하였다(Fig. 1). 기 체 투과 전산모사를 위한 3D 분리막 모델은 기존 선행 연구 결과와 비교하기 위하여, repeat unit 120개를 갖는 linear chain을 생성하여 이를 Amorphous cell 모듈의 입 력 구조로 활용하여 298K 조건에서 생성한 후 compress relaxation protocol[11,16,17]을 이용하여 최종 완성하였 다. 이렇게 생성된 3D 모델의 기체 투과 전산모사를 수 행하기 위하여, 기체 투과 실험에 널리 사용되고 있는 H₂, O₂, N₂, 및 CO₂의 네 종류 기체를 Sorption 모듈의 Fixed loading 기능을 이용하여 각각 모델에 도입하였 다. Sorption 모듈은 GCMC (grand canonical monte carlo) 기법에 기반한 Metropolis 알고리즘[18]을 통하여 흡착 분자의 삽입/삭제 과정을 계산해 나가는 방법이다[19]. 삽입/삭제 과정에서의 확률은 다음 식 (2)와 같이 계산 한다.

이때, U는 쿨롬 상호작용과 반데르발스 상호작용과 같은 non-bond에너지의 합이며, N_s는 삽입되는 기체 분 자의 수이다. 이때 에너지 변화량 ΔU가 음의 값을 갖 거나 Boltzmann factor, e ( − Δ U k b T )가 0과 1 사이에서 생성 된 특정값보다 크면 기체 분자가 삽입이 되게 된다.

고분자 linear chain의 초기 분포 및 기체 분자의 초 기 위치가 조금만 달라도 기체 확산에 있어서 매우 큰 영향을 끼치기 때문에, 이를 해결하기 위하여 가능한 초기 상태를 동일한 조건으로 유지시킬 수 있도록, 앞 선 과정은 단일 모델을 이용하여 진행하였으며, 기체 분 자를 삽입한 이후에 switching step[13]을 적용하여 각 force-field별 구조를 최종 제작하였다.

2.2. 전산모사 결과 분석

분자동역학을 통해 얻어진 각각의 3D 모델의 기체 확 산도와 고분자 linear chain의 움직임은 다음과 같은 식 (3)을 이용하여 계산되었다[4,20].

여기서, N_α는 확산하고 있는 기체 α의 개수를 나타내 고, r_i(0)과 r_i(t)는 기체 분자가 t시간 동안 이동하였을 때 각각 초기위치와 시간 t에서의 위치를 나타낸다. 이 때 〈 [ r i ( t ) − r i ( 0 ) ] 2 〉를 mean-square displacement (MSD) 라고 한다.

다음으로 고분자 linear chain이 서로 다른 force-field 하에서 어떻게 공간상에 분포하고 있는지를 계산하기 위 해, radius of gyration을 계산하였다. 식 (4)는 radius of gyration이 어떻게 정의되는지를 보여주고 있다.

여기서, N은 고분자 linear chain 내의 원자의 개수를 나타내고 있으며, r_i와 r_c,m은 각각 i번째 원자와 고분자 linear chain의 무게중심의 위치 벡터를 나타내고 있다.

3. 결과 및 고찰

3.1. Force-field에 따른 고분자 화학구조 변화

Fig. 1은 repeat unit의 화학구조가 force-field에 따라 서 어떻게 달라지는지를 나타내고 있다. 앞서 언급했던 것처럼, 본 연구의 주목적은 force-field가 맞지 않아서 잘못된 구조가 나오는 경우[13]를 최대한 배제하고 구 조적인 차이가 크게 없는 상황에서도 force-field에 따른 특성 차이가 나타날 수 있는지를 살펴보는 것이다. 그 런 의미에서 Fig. 1에 나타난 repeat unit의 geometry optimization결과는 Kapton^® 폴리이미드의 구조가 본 연구의 목적에 적절한 선택이었음을 확인시켜 준다. 벤 젠 구조의 공명구조도 잘 반영이 되어 있으며, 이미드 구조의 경우에도 전체 force-field에서 깨짐 없이 적절히 표시되고 있음을 확인할 수 있다. 일반적으로 이미드와 같은 상용 고분자의 경우 많은 전산모사 연구가 이미 수행되었으며, 이에 따라 각 force-field에도 이러한 이 미드 결합의 특성이 이미 잘 반영되어 파라미터화 되어 있기 때문으로 판단된다.

최종적으로 완성된 기체 투과 전산모사용 고분자 분 리막 3D 구조는 Fig. 2에서 확인할 수 있다. 다만, forcefield별, 기체분자별 전체 모델의 수는 총 20개에 이르 기 때문에, Fig. 2에서는 수소 기체가 포함된 모델만 force-field별로 도시하였다. 본 연구에서 살펴본 총 4종 류의 기체가 포함된 3D 구조는 초기 고분자 분포 및 기체 분자 분포가 최대한 동일하도록 설계가 되었기 때 문에, 기체별로 큰 틀에서 구조차이가 없는 것을 알 수 있다. 그럼에도 불구하고, 초기 구조와 비교하였을 때, 그 리고 각 force-field별로 고분자 linear chain의 torsional angle과 같은 conformation에 차이가 있는 것을 관찰할 수 있었다. 특히, 같은 functional form[13]을 사용하는 COMPASS와 pcff에서도 다른 force-field만큼이나 서로 다른 conformation 차이를 나타내었다.

3.2. Force-field에 따른 고분자 linear chain의 공간 분포 및 움직임 분석

이러한 구조적 차이는 이미지만으로는 명확한 차이를 파악할 수 없기 때문에, 좀 더 자세한 차이를 알아보기 위하여, radius of gyration 분석을 통하여 3D 공간 속에 분포되어 있는 고분자 linear chain의 반경을 분석하였다 (Fig. 3). 각 force-field별로 4종류의 기체 분자가 들어있 는 모델에서의 공간 분포를 살펴보면, 우선 COMPASS와 Universal의 경우 매우 작은 편차를 보여주고 있다. 특히 최대값과 최소값과의 차이가 0.2 Å 이내에 불과하였으 며, 그래프의 개형도 전체적으로 매우 비슷한 분포를 보여주고 있었다. 그러나 cvff와 pcff의 경우는 값의 편 차가 크게 나타났으며, 그래프의 개형도 다른 모양을 보 여주고 있었다. Dreiding의 경우는 편차는 크게 나타났 지만, 개형은 비슷한 또 다른 특징을 보여주고 있다. 다 만, 편차가 가장 큰 경우에도 그 차이가 0.5 Å 이내에 불 과하기 때문에, 앞에서 설명한 바와 같이 큰 틀에서 고 분자의 linear chain의 분포는 큰 차이가 없으며, forcefield별로 설정된 conformation 관련 파라미터들의 차이 로 인하여 조금씩 3D 공간 내에서의 구조차이가 나타 나게 된다는 것을 다시 한 번 확인할 수 있었다.

이러한 force-field에 따른 conformation 차이의 원인을 파악하기 위하여, Fig. 4와 같이 고분자 linear chain의 움직임을 mean square displacement (MSD) 계산을 통 하여 살펴보았다. 상대적으로 radius of gyration의 편차 값이 적었던, COMPASS와 Universal의 경우 3D 모델 내 에 존재하는 기체 분자의 종류에 관계없이 모든 모델에 서 linear chain의 움직임이 작은 것을 관찰할 수 있다. 반면에 cvff와 pcff에서는 기체의 종류에 따라서 MSD 결과, 즉 linear chain의 움직임이 매우 다르게 나타났으 며, 특히 CO₂와 O₂ 분자 존재 하에서 각 force-field의 MSD 값이 다른 force-field에 비하여 매우 큰 값을 보이 고 있었다. 또한, Dreiding의 경우 H₂와 O₂ 분자 존재 하 에서 MSD 값이 높게 나타났지만, cvff와 pcff의 값보다 는 작게 나타난 것을 알 수 있다. 이는 앞의 radius of gyration 결과와 정확히 일치하는 경향으로서, 고분자 전체 linear chain의 공간 분포는 크게 차이가 나지 않지만, 그 안에 존재하는 개별 원자들의 움직임은 force-field에 따라서 서로 다른 움직임을 보이고 있으며, 그것이 결과 적으로 linear chain의 공간상에서의 미세한 움직임에 차 이를 나타나게 하는 것을 알 수 있었다.

3.3. Force-field에 따른 고분자 분리막 3D 모델의 기체 확산도 분석

본 연구에 사용된 모델들은 동일 구조에서 출발하였 기 때문에 각각의 기체 종류 및 force-field별로 그 영향 이 확연하게 구분되며, 이는 기체 분자의 MSD 결과를 분 석함으로써 알 수 있다. Fig. 5에서 볼 수 있듯이, 같은 기체 분자가 도입된 같은 종류의 구조를 사용하여 동일 한 조건에서 분자동역학 전산모사를 수행하였음에도 force-field에 따른 MSD 개형이 완전히 다르게 나타나는 것을 확인할 수 있다. H₂ 분자의 경우 COMPASS에서 가장 높은 기울기를 나타내고 있으며 Dreiding이 그 다 음이고, 나머지 force-field의 경우 비교할 수 없을 정도 로 기울기가 낮은 것을 알 수 있었다. 특히, COMPASS 의 H₂ MSD 개형은 기존 연구 결과[17]와 비슷한 범위 값을 나타내고 있으나, 나머지 force-field의 경우 매우 낮은 값을 가지면서 비교 자체가 어려운 낮은 신뢰도의 값을 나타내었다. 특히, 같은 functional form을 사용하는 pcff에서도 이러한 경향이 나타났으며, 이는 특별히 구 조가 깨지는 등의 이상이 없고 구조상으로 유사하더라도 force-field의 차이만으로 완전히 다른 결과가 나올 수 있다는 것을 의미한다. N₂의 경우에도 매우 다른 MSD 개형을 보였으며, H₂에서와 마찬가지로 COMPASS는 기존 결과[17]와 비슷한 개형을 얻을 수 있었지만 다른 force-field는 이와 매우 동떨어진 결과를 보여주었고, 특 히 pcff와 cvff는 통계상으로 의미 없는 값을 나타내었다. 반면에 O₂의 경우에는 전체적으로 값의 편차가 다소 적 었으나, 역시 마찬가지로 force-field별 경향은 차이가 있 었다. COMPASS와 pcff가 비슷한 MSD 개형을 보였지 만 cvff는 상대적으로 낮은 값을 보였으며, Universal과 Dreiding은 통계적으로 무의미한 경향을 보여주었다. CO₂ 의 경우는 앞의 경우와 달리 pcff에서 월등히 높은 기 울기를 갖는 것을 관찰할 수 있었다. 반면에 Dreiding은 통계적으로 무의미한 결과를 보여주었다.

이러한 결과는 실제 기체 투과 전산모사에서 forcefield를 적용할 때 매우 신중을 기해야 한다는 것을 나 타내고 있다. 일반적으로 전산모사 모델의 1차적 검증 은 구조적 유효성과 밀도, 에너지의 안정성 등으로 판 단하게 된다. 그러나 이러한 조건들이 만족되었다고 해 서 기체 투과 결과가 항상 타당성 있게 나오지는 않는 다는 것을 본 연구에서는 보여주고 있다. 동일한 구조 라고 하더라도 기체 분자의 운동은 시간에 따라 지속적 으로 force-field의 영향 하에서 움직이기 때문에, 고분 자 linear chain과 같은 거대 분자에 비하여 그 영향을 훨씬 크게 받으며 따라서 어떤 기체 분자에서는 성공적 인 분자동역학 전산모사가 수행되더라도 반대로 특정 기체 분자에서는 제대로 된 결과를 얻지 못할 수도 있 다는 것을 본 연구에서 보여주고 있다.

4. 결 론

본 연구에서는 고분자 분리막 3D 모델의 기체 투과 를 관찰하기 위한 분자동역학 전산모사에서 force-field 가 고분자 구조와 기체 분자의 움직임에 미치는 영향을 관찰하고자 하였다. 이를 위하여 동일한 구조에서 출발 한 고분자 3D 모델을 이용하여 구조적 변수를 최소화 한 상태에서 분자동역학을 수행한 후에 force-field 종류 에 따른 고분자 구조, 3D 공간 내의 linear chain 분포, 기체 분자의 MSD 결과를 서로 비교하였다. 본 연구에 사용된 Kapton^® polyimide 구조는 여기서 사용된 모든 force-field에서 구조의 깨짐 없이 모두 정상적으로 모사 되었다. 그러나 force-field마다 다른 conformation 변수 들의 영향으로 부분적인 구조에서의 차이는 관찰이 되 었으며, 이는 radius of gyration과 고분자 linear chain의 MSD 결과를 통하여 분석되었다. 특히, 이러한 부분적인 구조의 차이는 force-field의 변화에 따른 원자 및 linear chain의 움직임의 차이에서 오는 것을 알 수 있었다. 이 러한 부분적인 구조차이와 달리 기체 분자의 움직임은 force-field별로 연관성을 찾기 어려울 정도로 매우 다른 경향을 나타냈으며, 이는 같은 functional form을 사용하 는 COMPASS와 pcff에서도 관찰되었다. 이는 고분자 와 같은 거대 구조의 경우 force-field에 따른 변수 값이 차이가 난다고 하더라도 전체 구조에 미치는 영향은 상 대적으로 작고 분자동역학 계산 중에 움직임이 제한적 이지만, 크기가 작은 기체 분자의 경우 긴 시간 동안 지속되는 분자동력학 계산 중에 force-field의 영향을 훨 씬 크게 받기 때문인 것을 알 수 있다.

특히, 분자동역학의 기술적인 관점에 보면, 현재 연 구자들이 사용하는 분자동역학 프로그램의 종류나 선 호도에 따라서 서로 다른 force-field를 사용하게 되는 데, 이로 인하여 동일한 구조를 갖고 동일한 조건에서 계산을 수행해도 서로 다른 결과가 나올 수 있다는 것 을 염두에 두고 충분한 검증과정을 통하여 전산모사 결 과에 대한 신뢰도를 높여야 할 것이다. 또한, 고분자와 같이 복잡한 구조를 갖는 모델의 경우 한 가지 프로그 램만으로는 구조 생성이 어렵기 때문에 다양한 프로그 램을 한 모델에 적용하는 경우가 많다. 예를 들면, 고분 자 구조는 그에 최적화된 특정 프로그램을 이용하여 생 성한 후, 다른 분자동역학 전용 프로그램에 이식하여 계 산을 수행하는 경우가 많은데, 이 경우 각 프로그램별로 지원되는 force-field가 다른 경우가 있으므로 충분한 주 의를 기울여 전산모사를 진행해야 할 것이다.