泊松─能斯特─普朗克類型方程的高效能計算及KcsA離子運輸模擬

Abstract

本論文針對泊松-能斯特-普朗克（Poisson-Nernst-Plank）類型的時變偏微分方程及其穩定態，提出高精度的數值格式，並實作到平行化、高效能的三維程式——NekCEM上，應用到生物上離子通道的大型模擬。 研究主要分成三部分，第一部分為解數值泊松-能斯特-普朗克方程。我們考慮一至三維度空間，選用高精度譜方法，其有限元素的特性將空間分成多域，使得資料結構適合平行化。時間方面，首先採用半顯隱式一階有限差分法;再來利用理察森外推法(Richardson extrapolation)，提高對時間差分的精度;最後我們也研究去雅可比-牛頓-克里夫法（Jacobian-free Newton-Krylov），加速得到穩定態解的過程。我們對這些方法進行數值實驗，觀察在不同邊界條件下，空間上近似解的誤差呈指數收斂，時間上達到二階收斂，符合理論預期。 第二部份我們將前述的數值格式應用到離子通道問題，我們拿KcsA鉀離子通道的蛋白質結構，在形狀複雜的三維通道上解泊松-能斯特-普朗克方程，從動態解到接近穩態，以模擬離子通過通道時的傳輸行為。由於離子通道的形狀曲折，對幾何形狀的不規則性做了預處理，而為了有足夠多的解析度，我們使用的網格點數超過一億，整個數值模擬的計算量和數據量都非常龐大。這部份使用了美國Argonne國家實驗室的超級電腦Mira跟Cetus，用了超過上千個節點，以及用其旗下的另一臺電腦Cooley上做資料的視覺化與後處理。 第三部份我們討論如何基於泊松-能斯特-普朗克方程的格式，加上位阻項及高階位阻項。這些高階非線性項需對顯隱式做額外的差值處理，以降低非線性項影響，使格式穩定進而收斂。而不同離子之間的耦合則必須解特別大的線性系統，這使得程式的編寫必須特別小心。我們會討論可能的邊界條件及呈現初步的數值實驗結果。

This thesis presents a high-order accurate scheme to solve time-dependent and time-independent Poisson-Nernst-Plank(PNP) type equations, including the implementation into a highly efficient parallel code--NekCEM, and large scale simulations on ion channel. For the first part, we focus on the numerical PDE methods. We use high-order spectral element method to discretize the space. As for the temporal discretization, first, we use first order finite difference with implicit-explicit methods. Then we use Richardson extrapolation to achieve second order accurate in time. Finally, we study Jacobian-free Newton-Krylov method to build a steady state solver. We provide numerical experiments to those methods under different kinds of boundary conditions, and successfully observing the exponential convergent in space and second order in time. For the second part, we apply our PNP schemes and codes to simulate ions transport in a potassium channel--KcsA. The channel structure comes from real biological data, which has to be taken care for its the complex three dimensional protein structure in pre-process. We solve time-dependent Poisson-Nernst-Plank equations on that complex geometry toward near steady state. To have enough the resolution, the total numbers of grid points are over a hundred million. The simulations are performed on supercomputers in Argonne Leadership Computing Facility (ALCF). We use over thousands nodes of Mira and Cetus, and using Cooley to do the visualizations. For the third part, we will discuss how to add the steric terms and the higher order steric terms into our PNP schemes. These leading nonlinear terms need an extra-extrapolation to increase the stability, and the coupling cross terms need a larger linear solver. We will discuss the implementation of some possible boundary conditions and some early state results will be shown.