電子擴散系統的快速計算法(2/3)

行政院國家科學委員會專題研究計畫 期中進度報告

電子擴散系統的快速計算法(2/3)

中 華 民 國 93 年 5 月 31 日

Fast Algorithms for Electro-Diffusional Systems

(Progress Report)

I-Liang Chern2

April 28, 2004

1

Background

Most biochemical processes involve macromolecules in solution. Usually, the macro-molecule is surrounded by a hydration shell and is immersed in mobile ionic solvent. The corresponding electrostatic potential (in or around) is of central importance for understanding their functions [16]. Two typical approaches are the molecular dy-namics simulation and the continuum modeling. The molecular dydy-namics simulation provides more adequate description at microscopic scale, but it requires a large num-ber of adjustable parameters and computing resources. Indeed, it is impractical for the spatial and time scale we are interested in. On the other hand, the continuum approach is simple, adjustable parameter free, relatively inexpensive, and has been widely used [16].

The continuum model for molecules in solution was pioneered by Debye-H¨uckel

in 1924 [2, 10]. The macromolecule such as protein or nucleotide represented by

a structured and polarized clusters of charges. It sits inside a region Ωin with low

dielectric constant 1, and is surrounded by a hydration shell Γ, which is immersed in

a mobile ionic solvent in region Ωout with high dielectric constant 2. The hydration

shell prevents the ionic solvent move into the inside region. Sometimes, the hydration shell could occupy a layer region (denoted by Ωlayer) with dielectric constant 2. In Ωout,

there could be several kinds of ions. For easy presentation of our numerical model, we assume there are only two kinds of ions with opposite sign and same charge unit, and with total neutralized solution. The corresponding electrostatic potential φ satisfies the Poisson-Boltzmann equation [2, 18, 19]:

−∇ · [(x)∇φ(x)] −X

i

K(x)ziexp(−ziφ(x)/kT ) = Q(x). (1.1)

Here, i is the index for the ith species of ions, (x) is the dielectric function taking value 1 in Ωin and 2 in Ωout∪ Ωlayer, K(x) = 0 for x ∈ Ωin∪ Ωlayer and K(x) = K for

1_{Project number NSC 93-2115-M-002-004}

x ∈ Ωout, a modified Debye-H¨ukel parameter, represents the ionic strength for solvent, the function Q(x) = 4π m X k=1 qkδ(x − xk) (1.2)

represents the charge distribution of the macromolecule in Ωin.

In the community of computational chemistry, there were decades of efforts of improvement on modeling by various groups of peoples such as Hornig (DelPhi pack-age), McCammon, etc. Previous numerical models can be classified into grid-based methods such as finite difference, finite element and finite volume [7, 20, 13], and boundary integral methods such as boundary element for integrating the discretized linear systems such as SOR [7], multigrid [6] have also attracted attentions. The nonlinear counterpart was solved by direct iteration [20], or Newton’s iteration [6]. However, most of above methods are low order, or not fast enough. Some of them indeed create severe numerical errors for large macromolecular such as DNA. For instance, one major error source comes from the treatment of singular point charge source. This point charge source term in a finite difference method is usually han-dled by direct discretization, or by distributing them into neighboring grid points.

This creates O(1/h1_{) error in L}∞ _{near macromolecules. In addition, the treatment}

of discontinuous dielectric constant is also naive. The coefficient smoothing method commonly used in this community has O(1) error in potential. The error may be larger for the force. These numerical problems leave room for mathematicians to do improvement.

2

What we have done so far

In these two years, we have done the following subprojects.

• “Accurate Evaluation of Electrostatics for Macromolecules in Solution,” (with Jian-Guo Liu and Wei-Cheng Wang), Methods and Applications of Analysis, Vol. 10, No. 2, pp. 309 (2003). The main point is to use multipole method to handle the singular source and to use a body-fitting grid method to handle interface discontinuities.

• “New Formulation and Fast Poisson Solvers for Interface Problems in Polar coordinates,” (with Zhilin Li, Wei-Cheng Wang and Ming-Chih Lai), The main thing addressed is to use a distance function to remove jump singularity of the source terms on interfaces. SIAM J. Sci. Comp., Vol. 25, No. 1, pp. 224-245 (2003).(SCI)

• “A spectral method for two-dimensional Poisson-Boltzmann equation,” (with Chien-Chang Yen, Jann-Guo Liu and Ming-Chen Shiu), in preparation. The work here is to develop a spectral method when the interface is a circle or a ring.

• “A New Immersed Interface Method for Elliptic Equations with Discontinuous Ciefficients,” (with Yu-Chen Shu), in preparation. In this research, we have develop a robust and simple second order finite method for elliptic equations with discontinuous coefficients. The underlying grid is regular, say Cartesian, and fixed. The interfaces may move and are represented by level sets. Interface problem is a well-konwn problem since the work of Peskin. The issue is to design a simple and high-order method under Cartesian grid and to allow interfaces

dynamically move. We have compared our method with previous works of

Peskin (immersed boundary method, see the review article [21]), LeVeque and Zhilin Li [25] (immersed interface method), Liu, Fedkiw and Kang [22] (ghost fluid method), our new method is the best. It is simple and second-order for both potentials and its gredients.

3

What we plan to do in the last year

3.1

The research subprojects

• We plan to write the paper for the new immersed interface method. • We can do one of the following researches in the third year:

– We plan to extend the new immersed interface method to 3-d.

– We plan to implement the new immersed interface method to Poisson-Boltzmann equations. The new element will be some nonlinear precondi-tioner.

– We plan to import some charge distribution data for some prototype pro-tein from the propro-tein data bank and to test our method for real drug design problems.

3.2

The difficulties and the approaches

• In 3-d, the difficulty needed to be resolved is the tangential derivative correction on the interfaces. Our 2-d method can be extended directly to 3-d without difficulty. However, the algebraic multigrid method for solving the resulting linear system is needed to tune up.

• For solving Poisson-Boltzmann equation, we plan to find some nonlinear pre-condotioner to speed up the nonlinear iterations.

• The real problem is an implementation problem.

3.3

Expected results

After this work, we should have a method which is better than any current popular methods (Hornig, McCommon, Holst, etc. see references in the bibliography).

References

[1] C. F. Anderson and M. T. Thomas, Jr., Ion distributions around DNA and other cylindrical polyions: theoretical descriptions and physical implications, Annual Review Biophysics and Biophysical Chemistry, 1990.

[2] Debye, P. and E. H¨uckel, Physik. Z., 24, 185 (1923).

[3] Gilson, M., K. Sharp and B. Honig, Calculating the electrostatic potential of molecules in solution: method and error assessment, J. Comp. Chem., 9, 327 (1987).

[4] Golub, G.H. and C.F. Van Loan, Matrix Computations, 3rd ed., Johns Hopkins University Press, Baltimore, 1996.

[5] Greengard, L. and V. Rokhlin, A fast algorithm for particle simulation, J. Com-put. Phys., 73 325 (1987).

[6] Holst, M. and F. Saied, Multigrid solution of the Poisson-Boltzmann equation, J. of Comp. Chem., 14, 105 (1993).

[7] Nicholls, A. and B. Honig, A rapid finite difference algorithm, utilizing successive over-relaxation to solve the Poisson-Boltzmann equation, J. Comp. Chem., 12, 435 (1991).

[8] Jackson, J.D., Classical Electrodynamics, 2nd Ed., Wiley, New York, 1975. [9] Kellogg, O.D., Foundations of Potential Theory, Dover, New York, 1953. [10] Kirkwood, J.G., J. Chem., Phys., 2, 351 (1934).

[11] LeVeque, R.J. and Z. Li, The immersed interface method for elliptic equations with discontinuous coefficients and singular sources, SIAM J. Numer. Anal., 31 1019 (1994).

[12] LeVeque, R.J. and Z. Li, SIAM J. Numer. Anal. (1997)

[13] Luty B., M. Davis and J. McCammon, Electrostatic energy calculations by a finite-difference method: rapid calculation of charge-solvent interaction energies, J. Comp. Chem., 13, 768 (1992).

[14] Peskin, C.S., Numerical analysis of blood flow in the heart, J. Comput. Phys., 25 220 (1977)

[15] Press, W.H., B.P. Flannery, S.A. Teukolsky and W.T. Vetterling, Numerical Recipes in C : The Art of Scientific Computing, 2nd ed., Cambridge University Press, New York, 1996.

[16] Rashin, A., Hydration phenomena, classical electrostatics, and the boundary element method, J. Phys. Chem., 94, 1725 (1990).

[17] Schaefer, M. and C. Froemmel, A precise analytical method for calculating the electrostatic energy of macromolecules in aqueous solution, J. Mol. Biol., 216, 1045 (1990).

[18] Setlow, R.B. and E.C. Pollard, Molecular Biophysics, Addison-Wesley, Boston, 1962.

[19] Sharp K. A. and B. Honig, Annu. Rev. Biophys. Biophys. Chem., 94, 7684 (1990). [20] Yoon, B. and A. Lenhoff, A boundary element method for molecular

electrostat-ics with electrolyte effects, J. Comp. Chem., 11, 1080 (1990).

[21] C. S. Peskin, The immersed boundary method, Acta Numerica, 1-39(2002). [22] X. D. Liu, R. Fedkiw, and M. Kang, A boundary condition capturing method for

Poisson’s equation on irregular domains, J. Comput. Phys, 160, 151-178(2000). [23] F. Gibou, R. Fedkiw, L. T. Cheng, and M. Kang, A second-order-accurate

sym-metric discretization of the poisson equation on irregular domains, J. Comput. Phys. 176, 205-227(2002)

[24] Z.Li, W.C. Wang, I-Liang, Chern, and M.C. Lai, New formulation and fast poisson solvers for interface problem in polar coordinates. NCSU CRSC-TR01-24(2001).

[25] R. LeVeque and Z. Li, The immersed interface method for elliptic equations with discontinuous coefficients and singular sources, SIAM J. Numer. Anal. 31, 1019-1044(1994).