Abstract:
A computational method to determine electrostatic interaction by performing direct numerical integration. The method recasts the Poisson equation and approximates the integral by using numerical integration schemes. Multi-dimensional integrals are separated into a coupled product of one-dimensional integrals. Linear transformations are performed and the total electrostatic potential is obtained as a sum of potential contributions for each integration point. The method is computationally efficient and well suited for parallel computers.

Description:
CROSS REFERENCE TO RELATED APPLICATIONS  
       [0001]     This application claims priority to U.S. Provisional application Ser. No. 60/580,205, filed Jun. 16, 2004, the disclosure of which is hereby incorporated by reference herein. 
     
    
     TECHNICAL FIELD  
       [0002]     This patent relates generally to molecular modeling and modeling of semiconductor devices. Specifically, the patent relates to the determination of electrostatic potential by using a direct computational approach.  
       BACKGROUND OF THE INVENTION  
       [0003]     Electrostatic potential in molecular systems is created by the nuclei and the electrons. Chemical reactivity and molecular interactions depend on the electrostatic potential. Electrostatic potentials are of fundamental importance in simulations of charging processes of semiconductor structures and devices. Electrostatic potential is a measurable physical quantity, but it is more commonly obtained in computer simulations. Nuclear contributions to the electrostatic potential can be obtained analytically, whereas calculation of the electronic contribution is more involved as it is created by the static electron distribution rather than by point charges.  
         [0004]     Because electrostatic potentials can be obtained analytically only for some simple systems, determination largely relies on numerical methods. Calculation of electrostatic potentials is straightforward, since electrostatic potentials-are generated as the integral of the reciprocal interparticle distance, |r 1 -r 2 |, times a general charge distribution. Newtonian gravitational potentials may be obtained using similar expressions as used in the calculation of electrostatic potentials.  
         [0005]     Such a direct approach has seldom been used as it is considered to be very complicated and time consuming due to the potential involving a six-dimensional space with singularities at r 1 =r 2 . The preferred method of solving for the electrostatic potentials is to numerically solve the Poisson equation in three dimensions. The Poisson equation is a second-order elliptical partial differential equation describing the electrostatic potential caused by a fixed charge distribution. In three dimensions, the Poisson equation is shown as: 
 
∇ 1   2 φ( x   1   , y   1   , z   1 )=−4πρ( x   1   , y   1   , z   1 ) 
 
         [0006]     Solving the Poisson equation numerically is a complicated task since it involves large linear matrix equations with crucial system specific boundary conditions. Because of its importance, much effort in the scientific community has been dedicated towards developing efficient methods of solving the Poisson equation.  
         [0007]     Other approaches have been useful for obtaining the electrostatic potential; however, these methods often provide qualitative rather than quantitative accuracy. Several numerical methods of determining the electrostatic potential in chemical and biological systems have been developed in computer programs such as, APBS, DelPhi, ITPACT and Manifold Code.  
         [0008]     Equally important applications of Poisson solvers are real space computational methods for electronic structure, the electrostatic potentials caused by the electrons are determined with high accuracy using the Poisson equation. However, this is one of the most time consuming steps in real space computations.  
         [0009]     Electrostatic potential is also of great importance in studies of semiconductor devices. Solutions to the Poisson equation for semiconductor structures and quantum dots provide information about their properties and physical insights into the single electron charging processes.  
         [0010]     What is needed is a technique that accurately determines electrostatic potential for complex molecular and semiconductor systems.  
       SUMMARY OF THE DISCLOSURE  
       [0011]     A method for solving differential equations for electrostatic potential is presented, in which the electrostatic potential may be estimated in a system without knowing boundary conditions. The method includes separating a multi-dimensional integral into a coupled product of multiple one dimensional (1D) integrals by applying an integral transformation and using numerical tensorial basis functions, and constructing matrices containing one-dimensional auxiliary integrals for each dimension. Further, the method approximates the auxiliary integral of the integral transformation by using numerical quadrature. Matrix multiplications are performed for each dimension and for each integration point in an auxiliary dimension. The differential contributions are numerically integrated to obtain the electrostatic potential for the system.  
         [0012]     A computer system is described that executes a software program in estimating the electrostatic potential for the system. The computer system includes a processor, a memory and a program stored in the memory and executable by the processor. The program incorporates the logic described above in estimating the electrostatic potential.  
         [0013]     The described embodiments solve the Poisson equation to computationally estimate the electrostatic potential. 
     
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
       [0014]      FIG. 1  depicts a logic diagram consistent with the teachings of the disclosure.  
         [0015]      FIG. 2  is a schematic diagram of a system capable of performing computations in accordance with the logic diagram of  FIG. 1 . 
     
    
     DETAILED DESCRIPTION  
       [0016]     The described embodiments recast an equation for calculating electrostatic potential into a more usable format and incorporate this equation into a computer software program. The electrostatic potential is defined as an integrated average of a charge distribution multiplied by a reciprocal distance between a position of a charge causing the potential and a potential coordinate. Mathematically expressed as:  
               ϕ   ⁡     (       x   1     ,     y   1     ,     z   1       )       =       ∫     -   ∞     ∞     ⁢       ∫     -   ∞     ∞     ⁢       ∫     -   ∞     ∞     ⁢       1     r   12       ⁢           ⁢     ρ   ⁡     (       x   2     ,     y   2     ,     z   2       )       ⁢     ⅆ     x   2       ⁢     ⅆ     y   2       ⁢     ⅆ     z   2                       (   1   )             
 
 where ρ(x 2  , y 2 , z 2  ) is a charge density,  
         r   12     =              r   1     -     r   2            =           (       x   1     -     x   2       )     2     +       (       y   1     -     y   2       )     2     +       (       z   1     -     z   2       )     2               
 
 is a distance, and φ(x 1 , y 1 , z 1 ) is the electrostatic potential. Thus, determination of φ(x 1 , y 1 , z 1 ) using Eq. (1) involves six spatial dimensions e.g. (x 1 , y 1 , z 1 ) and (x 2  , y 2 , z 2 ), and a singular function, because r 12  appears in the denominator. Singularities may be removed; by applying an integral transformation to recast the mathematical expression.  
                 ∫     -   ∞     ∞     ⁢       1     r   12       ⁢           ⁢     ρ   ⁡     (       x   2     ,     y   2     ,     z   2       )       ⁢     ⅆ     x   2       ⁢     ⅆ     y   2       ⁢     ⅆ     z   2           =       2     π       ⁢       ∫   0   ∞     ⁢       ∫     -   ∞     ∞     ⁢       ⅇ     -         t   2     ⁡     (       r   1     -     r   2       )       2         ⁢   ρ   ⁢     (       x   2     ,     y   2     ,     z   2       )     ⁢     ⅆ     x   2       ⁢     ⅆ     y   2       ⁢     ⅆ     z   2       ⁢     ⅆ   t                     (   2   )             
 
         [0017]     The integral transformation in Eq. (2) has been used in deriving an efficient recursion relation for the calculation of two electron integrals over Gaussian functions. Next the charge density ρ(x 2 , y 2 , z 2 ) may be expanded in a numerical tensorial basis. The tensorial basis consists of basis functions constructed as an outer product of one dimensional basis functions.  
               ρ   ⁡     (       x   2     ,     y   2     ,     z   2       )       =       ∑   αβγ     ⁢           ⁢         ⅆ     αβ   ⁢           ⁢   γ       ⁢       χ   α     ⁡     (     x   2     )         ⁢       χ   β     ⁡     (     y   2     )       ⁢       χ   γ     ⁡     (     z   2     )                   (   3   )             
 
         [0018]     In the three dimensional case, substitution of the density in Eq. (3) into Eq. (2) yields a separation of a three dimensional integral into a coupled product of three one dimensional integrals. This substitution also derives an expression for calculation of a potential φ(x 1 , y 1 , z 1 ) at selected points in space. Coordinates of the chosen potential points are shown in the exponent of the Gaussian function; the one dimensional integrals involving the Gaussian function times the basis function have to be calculated analytically or numerically for each potential point and basis function yielding a computational scaling that is proportional to N x   2 +N y   2 +N z   2 , where N x , N y , and N z  are the number of grid points in each dimension. Integration in the t direction may be performed numerically using Gaussian quadrature. Thus, an expression for the calculation of the potential in points (x 1 , y 1 , z 1 ) may be written as:  
               ϕ   ⁡     (       x   1     ,     y   1     ,     z   1       )       =       2     π       ⁢       ∑     α   t       ⁢           ⁢       w     α   t       ⁢       ∑   αβγ     ⁢           ⁢       ⅆ   αβγ     ⁢       ∫     -   ∞     ∞     ⁢       ⅇ     -         t     α   t     2     ⁡     (       x   1     -     x   2       )       2         ⁢           ⁢       χ   α     ⁡     (     x   2     )       ⁢       ∫     -   ∞     ∞     ⁢       ⅇ     -         t     α   t     2     ⁡     (       y   1     -     y   2       )       2         ⁢       χ   β     ⁡     (     y   2     )       ⁢           ⁢       ∫     -   ∞     ∞     ⁢       ⅇ     -         t     α   t     2     ⁡     (       z   1     -     z   2       )       2         ⁢       χ   γ     ⁡     (     z   2     )       ⁢     ⅆ     x   2       ⁢     ⅆ     y   2       ⁢     ⅆ     z   2                                     (   4   )             
 
 where integration points t α     t    and corresponding weights w α     t    have been introduced. The weights w α     t    are integration weights of the Gauss integration. Other numerical schemes may yield other weight factors. By denoting the integrals of the Gaussian function times the basis function χ γ     x    (x 2 ) for the calculation of the potential in point x α     x    by  
               F       γ   x     ⁢     α   x         x   ,     α   t         =       ∫     -   ∞     ∞     ⁢       ⅇ     -         t     α   t     2     ⁡     (       x     α   x       -     x   2       )       2         ⁢       χ     γ   x       ⁡     (     x   2     )       ⁢           ⁢     ⅆ     x   2                   (   5   )             
 
 and similar expressions for the y and z terms, the final expression can be written as  
               v       α   x     ⁢     α   y     ⁢     α   z         =       2     π       ⁢       ∑     α   t       ⁢           ⁢       w     α   t       ⁢       ∑     γ   z       ⁢           ⁢       F       γ   z     ⁢     α   z         z   ,     α   t         ⁢       ∑     γ   y       ⁢           ⁢       F       γ   y     ⁢     α   y         y   ,     α   t         ⁢       ∑     γ   x       ⁢           ⁢       F       γ   x     ⁢     α   x         x   ,     α   t         ⁢     ⅆ       γ   x     ⁢     γ   y     ⁢     γ   z                                   (   6   )             
 
 where ν α     x     α     y     α     z    denotes the electrostatic potential values for selected grid points. The evaluation of Eq. (6) includes three coupled matrix multiplications, the matrix size of which are N x ×N y , N y ×N z , and N z ×N y , respectively. Auxiliary integrals in F x,α     t   , F y,α     t   , and F y,α     t   , at can be calculated analytically using error function, but for small t values, the analytical expression suffers from numerical instabilities. However, for small t values the auxiliary integrals can be accurately obtained numerically by using, e.g., Gaussian quadrature. The matrix multiplications in Eq. (6) are performed for each grid point in the remaining direction (i.e. z, y, x) and for each t value. This leads to a computational scaling of  
         (         N   x   2     ⁢     N   y     ⁢     N   z       +       N   x     ⁢     N   y   2     ⁢     N   z       +       N   x     ⁢     N   y     ⁢     N   z   2         )     ⁢     N   t     ⁢           ⁢   or       
         3   ⁢     N   x   4     ⁢     N   t       =       3   ⁢     N     4   3       ⁢     N   t     ⁢           ⁢   when   ⁢           ⁢     N   x       =       N   y     =     N   z             
 
 N x N y N z =N are assumed. Thus, the method scales almost linearly with the total number of grid points as N t  is independent of grid size used. Two outer loop indices can be used for the distribution of the computational efforts to the processors of one or more parallel computers. Thus, the increase in speed should be substantially linear because the distributed tasks consist of matrix multiplications with no requirement for communication between the processors. 
 
         [0019]     Referring now to  FIG. 1 . The software implementation of Equation 6 is a series of calculation loops, several smaller loops at  24 ,  28 , and  32  within one larger loop at  38 . The number of smaller loops is dependent on the number of dimensions in the system. The larger loop at  38  represents computations for each integration point in the t dimension. In this operation, the Einstein summation convention is used. The Einstein summation convention implies that when an index occurs more than once in the same expression, the expression is implicitly summed over all possible values for the index. The auxiliary integrals are constructed at  20  for each dimension and at each grid point in space. For large t α     t   , values, the F x,α     t   , the F y,α     t    and the F z,α     t    matrices are band dominant, a property that may be used for acceleration of the computational speed. Linear transformation in the x dimension is performed at  24 . The matrix multiplications for the external indices α 1 , and y z  may be performed on one or more parallel processors, thus accelerating computational speed.  
         [0020]     An optional reorder function is performed at  22 ; this function reorganizes the data into a more rapidly accessible form. This reorder function is introduced to make the matrix multiplications as fast as possible, any other manner of reordering the terms or no reordering at all would be acceptable. Additional optional reorder functions at  26 ,  30  and  34  similarly reorganize the data.  
         [0021]     Linear transformation in the y dimension is performed at  28 . In this operation, the Einstein summation convention is used again. The α t , and α x  indices are external indices and the corresponding matrix multiplications may be performed on one or more parallel processors. Linear transformation in the z dimension is performed at  32 . In this operation, the Einstein summation convention is also used. The α t  and α x  are external indices and the corresponding matrix multiplications may be performed on one or more parallel processors.  
         [0022]     Contributions to the electrostatic potential for each integration point in the t dimension are multiplied by the integration weight factor and added to obtain the total electrostatic potential at  36 .  
         [0023]      FIG. 2  is a schematic diagram of one embodiment of a computer system used to calculate electrostatic potential. The computer  60  is operatively connected to one or more processors  68 , a memory  62 , an output device  66  and an input device  64 . An executable program  70  is stored in the memory  62  and accessible by the computer  60 . The executable program  70  may use the logic described in  FIG. 1 . The executable program  70  calculates, in the spatial grid points, the contribution to the electrostatic potential of each t α     t    value and these contributions are summed to obtain the total electrostatic potential for the system. A notation is used wherein matrix elements with increasing first lower indices lie subsequently in the computer memory  62 , thus allowing the computer  60  to keep the values as long as necessary and possible in the cache memory, yielding an accelerated computational speed. The computer  60  may divide the matrix multiplications for the outer indices α t  and γ z  at  24  and for the outer indices at and ax at  28  and  32  between the processors  68 , thereby accelerating the computational process.  
         [0024]     In one embodiment, a tensorial product of Lagrange interpolation functions may be used as a basis function. In this numerical representation, expansion coefficients of the functions are amplitudes of the functions in the grid points. Element functions of arbitrary order may be employed; however, in the described embodiment, second, fourth and sixth order Lagrange interpolation functions have been used. Other kinds of local basis functions, such as, for example, wavelets, splines or any other basis set which can be expressed as a tensor product of the one-dimensional (1D) basis functions may be used. One benefit of this embodiment is that in solving more complicated systems, where higher-order element functions result in more accurate potentials, the higher-order element functions may be performed with almost no additional computational costs (i.e., without slowing down the program).  
         [0025]     The described embodiments are directed at solving the Poisson equation in three dimensions. The disclosure could be adapted in other embodiments to solve the Poisson equation in any number of dimensions as well as solving other types of nonlinear Poisson-Boltzmann equations and other related differential equations such as, for example, Schrödinger equations. In the case of non-linear differential equations, the method may be used iteratively to obtain a solution. For example, in the case of the Schrödinger equation, the method begins with an arbitrary initial guess for wavefunction and energy. Integration. may be performed and new energy values used in a subsequent integrations. The process may be repeated until the energy value converges. One skilled in the art will realize the wide range of equations that may be solved by various implementations of the disclosure.