Enabling Ab-Initio Molecular Dynamics under Bias: The CP2K+SMEAGOL Interface for Integrating Density Functional Theory and Non-Equilibrium Green Functions

기술 요약: CP2K+SMEAGOL 인터페이스를 통한 바이어스 하의 Ab-Initio 분자 동역학 구현

문제 제기
밀도 범함수 이론(DFT)과 비평형 그린 함수(NEGF)의 결합은 외부 바이어스 전위 하에서의 양자 수송을 모델링하는 표준적인 접근 방식이지만, 기존 방법론은 정적 시스템에 국한되어 있습니다. 기존의 구현 방식들은 대개 고정된 원자 위치에서 수송 특성을 계산하는 정준(canonical) 프레임워크 내에서 작동합니다. 이는 이론적 이해와 거시적 실험 결과 사이의 상당한 간극을 만듭니다. 특히 전기화학(EC) 시스템과 같이 변환이 본질적으로 그랜드 캐노니컬(grand canonical)인 경우, 이러한 간극은 더욱 두드러집니다. 전기화학적 변환의 자연스러운 환경을 정확하게 묘사하기 위해서는 바이어스 하에서의 분자 동역학(MD)과 계산 셀을 드나들 수 있는 전자의 명시적인 개방 경계 기술이 필요합니다. 또한, 기존의 DFT-NEGF 구현은 구조 최적화 및 동역학에 필요한 전류 유도 힘(current-induced forces)을 계산하는 기능이 부족하거나, 응축상 시뮬레이션에 최적화된 최첨단 패키지의 기능들을 갖추지 못한 특정 DFT 코드(예: SIESTA)에 국한되어 있는 경우가 많습니다.

방법론
저자들은 대중적인 혼합 가우시안/평면파 전자 구조 패키지인 CP2K와 NEGF 코드인 SMEAGOL 사이의 새로운 인터페이스를 제시합니다. CP2K+SMEAGOL이라 명명된 이 인터페이스는 원래의 SMEAGOL 코드를 최대한 재사용하면서, 어떤 DFT 소프트웨어와도 연결 가능한 독립형 라이브러리로 설계되었습니다.

주요 기술적 구현 사항은 다음과 같습니다:

• 전류 유도 힘(Current-Induced Forces): 이 인터페이스는 바이어스 하에서 이온에 작용하는 비보존적 힘(non-conservative forces)의 계산을 구현합니다. 이 힘은 원자의 운동량에 대한 전체 해밀토니안의 기울기를 활용하여 원자의 운동량 시간 미분으로부터 유도됩니다. 결정적으로, 이 구현은 밀도 행렬($\Omega_{ij}$)의 미분에서 발생하는 힘 항을 포함하며, 이는 "전자 바람(electron wind)" 성분(전자가 이온에 전달하는 운동량)을 설명합니다. 이 항은 다른 구현체에서 종종 무시되지만, 질적인 정확성을 확보하는 데 필수적임을 보여줍니다.
• 밀도 행렬 구축: 바이어스 윈도우 내의 결합 상태(bound states)를 올바르게 점유하기 위해, 인터페이스는 밀도 행렬 계산을 업데이트합니다. 단일 컨투어 평가나 좌/우 기여도의 단순 평균을 사용하는 대신, Brandbyge 등의 방식을 따라 비평형 기여도에 의해 결정되는 가중치를 가진 가중 이중 컨투어(weighted double-contour) 스킴을 채택하여 양쪽 리드(lead)로부터의 적분을 합산합니다.
• 시뮬레이션 설정: 이 방법은 바이어스 하에서의 Born-Oppenheimer 분자 동역학을 허용합니다. 시스템은 반무한 전극에 부착된 중앙의 "확장된 분자(extended molecule)" 영역으로 모델링됩니다. 리드의 자기 에너지(self-energies)는 중앙 영역의 전하 밀도와 독립적이며 정적인 것으로 취급되는데, 이는 중앙 영역이 전하 변화를 차폐할 만큼 충분히 큰 경우 유효한 가정입니다.

주요 결과 및 검증
저자들은 여러 시스템에 대해 SIESTA+SREAGOL을 이용한 참조 계산을 통해 CP2K+SMEAGOL 인터페이스를 검증하였습니다:

1. 제로 바이어스 힘 (무한 Au 와이어): 인터페이스는 제로 바이어스에서의 원자력을 성공적으로 재현하였으며, 이는 CP2K 단독 계산과 일치합니다. 또한 에너지 밀도 행렬($\Omega$) 항을 무시할 경우 질적으로 잘못된 힘이 발생함을 입증하였습니다.
2. 평행판 커패시터: 인터피스는 커패시터 전반에 걸친 기대 전위 강하를 재현하며, 리드에서는 일정한 전위를, 진공 영역에서는 선형적인 전위 강하를 보여주었습니다. 또한 가중 이중 컨투어 스킴이 대칭적인 전하 밀도와 올바른 하트리(Hartree) 포텐셜을 생성하는 데 필수적임을 확인하였으며, 단일 컨투어 평가는 비대칭적이고 잘못된 결과를 낳음을 보여주었습니다.
3. Au-H$_2$-Au 접합 구조 최적화: 인가된 바이어스(최대 1.5 V) 하에서, 인터페이스는 H-H 결합의 신장과 원자 변위를 정확하게 예측하였으며, 이는 이전의 SIESTA+SMEAGOL 결과와 일치합니다. 이러한 변화는 전기장이 아닌 전류에 기인한 것으로 분석되었습니다.
4. 용매화된 Au 와이어의 분자 동역학: 저자들은 대규모 응축상 시스템(166개의 물 분자로 용매화된 단원자 Au 와이어)에 대해 최초의 DFT-NEGF 분자 동역학 시뮬레이션을 수행했습니다.
   • 안정성: 개방형 시스템임에도 불구하고, 역학은 최소한의 장기 에너지 드리프트와 함께 안정적으로 유지되었습니다.
   • 전자 이동: 1 V 바이어스에서의 200 fs 궤적 동안, 시뮬레이션은 전자 이동 효과를 관찰했습니다. 즉, 국부 전하 밀도의 최소점이 중앙 와이어 원자에서 오른쪽 끝 원자로 이동하였는데, 이는 오른쪽 전극의 양전하를 차폐하는 현상과 일치합니다.
   • 투과(Transmission): 용매화된 와이어의 투과 스펙트럼은 이상적인 무한 와이어와 질적으로 유사한 모습을 보였으며, 차이점은 물 분자에 의해 추가로 열린 터널링 채널 때문인 것으로 나타났습니다.

성능
주된 계산 병목 지점은 그린 함수의 적분으로서 전자 밀도를 평가하는 것이며, 이는 단일 해밀토니안 대각화를 대신하여 SCF 단계당 96개의 에너지 포인트 계산(실수축 64개, 복소축 32개)을 요구합니다.

• 병렬화: 코드는 실공간 계산 및 k-포인트를 위해 MPI를 사용하고, DO 루프 및 LAPACK 루틴을 위해 OpenMP를 사용합니다. 성능은 32개의 OpenMP 스레드까지 잘 확장됩니다.
• 비용: 용매화된 Au 와이어 시스템의 경우, CP2K+SMEAGOL 계산은 8배 더 많은 코어(2048 vs 256)를 사용했음에도 불구하고 표준 CP2K 계산보다 SCF 단계당 약 3.2배 느렸습니다. 저자들은 이것이 한 차원 높은 비용 증가를 의미하며, 추가적인 가속화(예: ScaLAPACK 또는 머신 러닝 힘 필드) 없이는 장기 MD 궤적 수행이 어려울 수 있다고 언급했습니다.

의의 및 주장
본 논문은 CP2K+SMEAGOL이 현재 CP2K에서 사용할 수 있는 가장 기능이 풍부하고 완전한 DFT-NEGF 구현체라고 주장하며, 이는 $\Gamma$-지점 수송(transport)에 국한되고 힘(force) 계산이 불가능한 기존의 네이티브 CP2K 구현을 능가하는 것입니다.

이 연구의 주요 의의는 DFT-NEGF를 사용하여 실제 작동 조건 하의 대규모 응축상 시스템에 대한 분자 동역학 시뮬레이션을 수행할 수 있음을 입증했다는 데 있습니다. 저자들이 아는 바로는, 이는 이러한 계산이 수행된 첫 번째 사례입니다. 이 인터페이스는 인가된 전위, 전류 흐름, 그리고 원자 동역학 사이의 상호작용이 결정적인 전기 촉매 반응 및 전기화학 시스템 연구를 가능하게 하여, 미시적 이론 모델과 거시적 실험 조건 사이의 간극을 메워줍니다.