Revista Colombiana de F ´ ´ısica, Vol. 43, No. 2 de 2011.
Solucion o´ n De La Ecuacion o´ n De Poisson-Boltzmann Usando Diferencias Finitas Solution Of The Poisson-Boltzmann Equation Using Finite Differences F. Fonseca * a , R. Mart´ Mart´ınez ınez a , P. Teher´ Teheran a´ n a a
Departamento de F´ısica, ısica, Univ ersidad Nacional de Colombia, Sede Bogot´ a.
Recibido 15.03.11; Aceptado 14.04.11; Publicado en l´ l´ınea ınea 04.09.11.
Resumen Se presenta la aplicaci´ aplicacion o´ n del m´ metodo e´ todo de diferencias finitas finitas a la soluci´ solucion o´ n de la ecuaci´ ecuacion o´ n de Poisson-Boltzmann en dos dimensiones. Se estudian tres posibles m´etodos etodos de soluci´on, on, Jacobi, Gauss-Sediel Gauss-Sediel y sobrerelajaci´ sobrerelajaci´ on sucesiva, encontr´andose andose que la diferencia entre los tres algoritmos yace en la velocidad de convergencia.
Palabras Clave: Poisson-Boltzmann; Diferencias finitas.
Abstract We present the application of finite difference method to solve the Poisson-Boltzmann equation in two dimensions. We studied three possible methods of solution, Jacobi, Gauss-Sediel and successive overrelaxation. We found that the difference between the three algorithms lies in the speed of convergence.
Keywords: Poisson-Boltzmann; Finite differences. PACS: 31.15.xf; 31.15.E-. c 2011. Revista Colombiana de F´ F´ısica. ısica. Todos los derechos reservados.
Intro Introduc ducci ci´on o´ n
1.
Los modelo modeloss contin continuos uos de mol´ moleculas, e´ culas, en solucione solucioness ionicas, o´ nicas, son herramie herramientas ntas muy importan importantes tes en el estudio estudio de las interacciones electrost aticas. a´ ticas. El primer modelo propuesto, fue elaborado por primera vez en 1923 por Debye y Huckel, u¨ ckel, [1]. [1]. Las interacciones electrost´aticas aticas juegan un papel central en la determinaci´on on de la estructura y comportamiento de prote´ prote´ınas ınas y sistemas complejos como los virus y las enzimas. Por otra parte, se ha realizado investigaci on o´ n en forma amplia en poliiones, incluyendo esferas de poliestireno, part´ part´ıculas ıculas elongadas como el virus del mosaico del tabaco, en cadenas r´ıgidas ıgidas de polielectrolitos, polielectrolitos, como el ADN , en membranas flexibles, en esferoides de silicato de arcilla, etc. *
[email protected]
Cuando las biomol eculas e´ culas residen en un electrolito acuo´ y funcionamiento debido so, este afecta su conformaci on a los efectos de apantallamiento diel ectrico. e´ ctrico. Adem´ Ademas, a´ s, las mol´eculas eculas polarizables u orientables de un medio hu´esped, sirven para disminuir las fuerzas el ectricas e´ ctricas que determinan la estructura y funci on ´ de macromol eculas. e´ culas. Por lo anterior, un c alculo a´ lculo exacto del potencial electrost atico a´ tico puede mejorar la comprensi on o´ n del comportamiento y la estructura de las macromol´eculas. eculas. Por otra parte, desde el punto de vista matem atico, a´ tico, la ´ diferencial ecuaci´ ecuacion o´ n de Poisson-Boltzmann es una ecuaci on parcial, altamente no-lineal, para un potencial electrost atico a´ tico local, la cual puede ser resuelta anal´ anal ´ıticamente ıticamente para unos pocos casos. Sin embargo, la teor´ıa ıa de Poisson-Boltzmann Poisson-Boltzmann se vuelve inmanejable para geometr´ geometr ´ıas ıas no triviales, y por lo ge-
F. Fonseca, R. Mart ´ınez, P. Teher ´ an: Soluci´ on De La Ecuaci´ on De Poisson-Boltzmann Usando Diferencias Finitas
neral debe ser linealizada para tratarse m a´ s f a´ cilmente desde un punto de vista anal´ıtico o en una aproximaci o´ n num´ericamente m´as estable.
Cambiando el potencial termodin a´ mico y la funci o´ n independiente, la condici o´ n de minimizaci o´ n (2) se puede expresar como
δ Ω[ρ] = u(r). δρ(r)
En este trabajo se presenta la soluci´on de la ecuaci´on de Poisson-Boltzmann usando discretizaci o´ n en diferencias finitas, en dos dimensiones. Se exploran los m e´ todos de Jacobi, Gauss-Seidel y sobrerelajaci o´ n sucesiva. En la secci o´ n (2) se presenta una revisi o´ n de la formulaci o´ n funcional de la teor´ıa de Poisson-Boltzmann, el cual representa una aproximaci´on te´orica muy interesante al problema. En la secci o´ n (3) se muestra el m e´ todo de discretizaci o´ n en diferencias finitas y se aplica a la ecuaci o´ n de Poisson-Boltzmann. En la seccio´ n (4), ofrecemos los resultados y los comparamos seg´un los diferentes m´etodos. Por u´ ltimo la secci´on (5) se discuten las conclusiones y en la ( 6) se dan los agradecimientos.
(4)
Por lo tanto, un sistema con un potencial externo dado es ´ asociado con un unico densidad de equilibrio y viceversa. Se considera una aproximaci o´ n de campo medio significa para un sistema que contiene una densidad de iones ρ± y una densidad de campo externo q (r), de tal manera que la densidad de carga es eq (r). Se considera como sistema de referencia el gas ideal. La interacci o´ n de larga distancia entre dos iones de la misma clase es v1,α,α ( r, r ) y entre los iones y la distribuci´on de carga externa es denotado v1,α,q (r, r ) = v1,q,α (r, r ) y se asume que ambas distribuciones obedecen la ley de Coulomb,
2.
La Ecuaci´on de Poisson-Boltzmann
v1,α,α
e2 αq v1,α,q (r, r ) = |r − r | La funcional de la energ´ıa es
(7)
α
1,α,α
α
1,α,q
(r, r )ρα (r )
α,α =±
+
1 2
(r, q )q (r );
α,α =±
usando la funcional de las propiedades intr´ınsecas de la energ´ıa libre de Helmholtz de un gas ideal como [4]
F id [ρ] =
drρ(r)(ln( ρ(r) ) − 1), ρs
α,α =±
β drρ(r)(ln( ρ(r) ) − 1) ρ 1
(8)
entonces la ecuaci´on (7) se convierte en
(1)
βF [ρ] =
(9)
s
α,α =±
+
2
drC (r)φ(r),
donde ρs es la densidad de iones del bloque, y la densidad de r ) = −ρ− (r) + ρ+ (r) + q (r) y el potencial carga total es C ( electrost´atico φ es definido como:
(2)
donde ρ es la densidad de equilibrio. La energ´ıa libre intr´ınseca de Helmholtz se puede definir con la siguiente transformaci´on de Legendre.
drρ0 (r)u(r).
dr ρ (r)v dr dr ρ (r)v
βF [ρ] = F id [ρ] 1 + dr 2
donde µ es el potencial qu´ımico. Usando la derivada funcional, la variaci o´ n del gran potencial con respecto al potencial intr´ınseco es:
F [ρ0 ] = Ω[u] +
(6)
(5)
Partiendo de la densidad funcional Ω para un sistema de N part´ıculas id´enticas, interactuando con un por1 , r2 , r3 , ..., rN ) y un potencial externo V = tencial φ( N ri ), se define el potencial qu´ımico local intr´ınsei=1 V ext ( co como:
δ Ω[u] = −ρ0 (r), δu(r)
y
La teor´ıa de la densidad Funcional (DFT) se puede usar para determinar la estructura de un sistema de muchos cuerpo utilizando funcionales, tales como la densidad electr o´ nica dependiente del espacio. Esta teor´ıa de forma general se basa en suponer que las propiedades de equilibrio del sistema est a´ n determinadas u´ nicamente por la densidad de part´ıculas, y por lo tanto el sistema de muchos cuerpos, de 3N coordenadas espaciales, se reduce a 3 coordenadas espaciales a trav´es del uso de funcionales. De la misma forma, se define y minimiza un funcional de energ´ıa para el sistema, y por lo tanto emerge la densidad de part ´ıculas, que es la densidad de equilibrio del sistema. Un tratamiento m a´ s profundo a la teor´ıa del DFT, puede ser hallado en la siguiente bibliograf´ıa [2] y [3].
u(r) = µ − V ext (r),
e2 αα (r, r ) = , |r − r |
φ(r) = λB
C (r ) dr − , |r r |
(10)
con λB la longitud Bjerrum que se define como la distancia a la que la interacci o´ n electrost´atica entre dos cargas es
(3) 237
Rev.Col.F ´ ıs., Vol. 43, No. 2 de 2011.
igual a la unidad K B T . En la funcional (9), hacemos el proceso de minimizaci on ´ con respecto a las densidades i onicas ´ ρ± , lo cual conduce a la ecuaci o´ n diferencial parcial el´ıptica no-lineal de Poisson-Boltzmann:
conectar estos valores iniciales se produce un valor al lado derecho V 1i, j , el cual se propaga a la red. Este proceso se repite un n u´ mero finito de veces de tal forma que la soluci o´ n converge. El esquema es:
− ∇ · ((r)∇φ(r)) + κ2 sinh(φ(r)) = 4πλB q (r) en R2 ,
+1 φm i,j
(11) r) es donde κ es la longitud de apantallamiento de Debye, ( el diel´ectrico que en general puede depender del espacio, en r ) es la densidad de carga nuestro caso ser a´ uniforme, y q ( externa. Para los iones m o´ viles en el solvente estar a´ dada por una distribuci´on de Boltzmann, y en este caso esta representada por dos clases de iones de carga inversa en equilibrio estad´ıstico, el cual est a´ dado por el sinh.
3.
En el m´etodo iterativo de Gauss-Seidel se hace una modificaci´on del m´etodo de Jacobi, demostr´andose que converge m´as r´apido. Supongamos que barremos la red con el fin de incrementar i y j , luego los t´erminos de la izquierda y m´as bajo de cada sitio en la red se han actualizado. Esto da lugar a una forma del algoritmo de Gauss-Seidel:
´ en diferencias finitas para dos dimenDiscretizacion siones
+1 φm i,j
El proceso de discretizaci o´ n se lleva a cabo sobre las r). Por lo tancoordenadas espaciales x y y , y la funci o´ n φ( to, se establece una escalas en el sistema, h, asociada a la discretizaci´on espacial, la cual se asume igual, en ambas direcciones x y y. Entonces, las derivadas espaciales quedan definidas como: m φm ∂φ i+1,j − φi,j = , ∂x h
(12)
m φm ∂φ i,j +1 − φi,j = , ∂y h
(13)
ω m+1 (φ + φm i,j −1 (19) 4 i,j+1 +1 2 2 2 m m m +φm i+1,j + φi−1,j + h 4πλ B q (i,j ) − κ h sin(φi,j )). +1 φm = (1 − ω)φm i,j + i,j
Como se necesita definir cu a´ ndo el proceso iterativo ha sido lo suficientemente exitoso, o cu a´ ndo la convergencia se ha alcanzado, un criterio es solicitar que la soluci´on aproximada deje de cambiar significativamente de una iteraci´on a la siguiente. Una forma de estimar este error relativo es calcular el promedio sobre todos los sitios de la red que experimentan un cambio distinto de cero. Por lo tanto se define:
(14)
m m φm ∂ 2 φ i,j +1 − 2φi,j + φi,j −1 (15) = . ∂y 2 h2 Los sub´ındices i y j representan un ro´ tulo para la discretizaci´on espacial y el super´ındice m representa la discretizaci o´ n
(|1 − error =
temporal, de la misma forma se tomar a´ como la cantidad de iteraciones necesarias para validar el m e´ todo de solucio´ n. Aplicando este esquema de discretizaci´on a la ecuaci´on (1), tenemos:
φm i,j
1 m+1 m+1 m (φ + φm i,j −1 + φi+1,j + φi−1,j (18) 4 i,j+1 2 2 m + h2 4πλB q (m i,j ) − κ h sin(φi,j )). =
Los m´etodos de Jacobi y Gauss-Seidel no utilizan el valor de V ij en el mismo punto de la red para actualizar el valor de V i,j . Resulta que el convergencia de la iteraci´on se puede mejorar considerablemente mediante el uso de una combinaci´on lineal de las soluciones nuevas y viejas de la siguiente manera:
y la segunda derivada m m φm ∂ 2 φ i+1,j − 2φi,j + φi−1,j = , ∂x 2 h2
1 m m m (φi,j+1 + φm i,j −1 + φi+1,j + φi−1,j (17) 4 2 2 m + h2 4πλB q (m i,j ) − κ h sin(φi,j )). =
L
i,j
φi,j |) φ(Nuevo)i,j
(20)
=0 para φ(Nuevo)i,j
1 m m m (φ + φm i,j −1 + φi+1,j + φi−1,j (16) 4 i,j +1 2 2 m + h2 4πλB q (m i,j ) − κ h sin(φi,j )). =
Por u´ ltimo, para lograr mayor estabilidad en el sistema, la no linealidad, es decir el t´ermino sin(φm i,j )), en la ecuaci´on de Poisson-Boltzmann, se convierte en:
Para resolver el sistema computacionalmente se explicar´an, brevemente, las que pueden ser consideradas tres t´ecnicas b´asicas para resolver el sistema.
1 (sin(φm i+1,j )) + 4 m m sin(φm i−1,j )) + sin(φi,j +1 )) + sin(φi,j −1 ))), sin(φm i,j )) →
La primera es el m´etodo iterativo de Jacobi para resolver ecuaciones simult a´ neas. Para ello se empieza con un valor inicial V 0i,j para la soluci o´ n de los nodos internos. Al
(21)
o un promedio espacial sobre sus vecinos m a´ s cercanos. 238
F. Fonseca, R. Mart ´ınez, P. Teher ´ an: Soluci´ on De La Ecuaci´ on De Poisson-Boltzmann Usando Diferencias Finitas
4.
Resultados
Fig. 4: Campo gradiente y estructura de las curvas de nivel para la ecuaci´on (19). Fig. 1: Soluci´on computacional usando el m´etodo de sobrerelajaci´on, ec. (19).
En las figuras (1), (2) y (3) se presentan los resultados para una simulaci´on tridimensional, con un tama˜no del sistema de L = 104, escalado a 1, generando una malla de 10816 puntos internos. Las condiciones de frontera de Dirichlet que se imponen en el sistema son:
φ(0, y) = 0, φ(L, y) = 8, φ(x, 0) = −5, φ(x, L) = −10.
(22)
De la misma forma el sistema se inicializa con 4 puntos iniciales de potencial que sirven como fuentes:
φ(L/4, L/4) = −12, φ(L/4, 3L/4) = 12, φ(3L/4, L/4) = −12, φ(3L/4, 3L/4) = 12.
Fig. 2: Estructura del campo gradiente superpuesto a las curvas de nivel para la ecuaci´on (19).
(23)
Para los tres m´etodos de soluci´on planteados, Jacobi, Gauss-Seidel y sobrerelajaci o´ n sucesiva, los resultados gr´aficos cualitativamente no son muy diferentes. Por lo anterior, so´ lo se presentan los resultados obtenidos con el m e´ todo de sobrerelajaci´on sucesiva, figuras (1), (2) y (3). Las diferencias aparecen en el n u´ mero de iteraciones para cada uno de los m e´ todos. En los tres casos se define un error de 1 % para alcanzar la convergencia en el m e´ todo num´erico. Para Jacobi la convergencia se alcanza para un n´umero de 249 iteraciones, para Gauus-Seidel el n´umero de iteraciones es 203 y finalmente para la sobrerelajaci on ´ sucesiva el numero ´ de iteraciones fue 216, con un valor del par´ametro de sobre-relajaci o´ n de ω = 0,56.
Fig. 3: Estructura tridimensional de las curvas de nivel para la ecuacio´ n (19).
239
Rev.Col.F ´ ıs., Vol. 43, No. 2 de 2011.
La figura (2) presenta la estructura del campo gradiente, superpuesta a las curvas de nivel, bidimensionalmente. Los m´etodos resuelven con suficiente eficiencia la distribuci o´ n de las curvas de nivel sobre las fronteras y sobre las cargas fuente.
les, tales como c´ırculos, elipses y estadios.
6.
Agradecimientos
F. Fonseca agradece la financiaci o´ n de este trabajo por parte de la Universidad Nacional de Colombia, en su Divisi´on de Investigaci o´ n sede Bogot´a, con n´umero de proyecto (DIB-8003355).
La figura (3), presenta la estructura tridimensional de las curvas de nivel, y la figura (4) muestra la distribuci o´ n de valores de φ donde la escala de colores representa los valores de campo escalar sobre la malla.
Referencias 5.
Conclusiones [1] P. Debye and E. H¨uckel, Physik. Z., 24, 185 (1923).
Se ha implementado la soluci´on computacional de la ecuaci´on de Poisson-Boltzmann usando diferencias finitas, para los m´etodos de Jacobi, Gauss-Seidel y sobrerelajaci o´ n sucesiva en dos dimensiones. Los tres m e´ todos resuelven el sistema de forma efectiva, diferenci a´ ndose en su velocidad de convergencia.
[2] R. van Roij, Soft Condensed Matter Theory, course book (2009). [3] J.-L. Barrat, and J.-P Hansen, Basic Concepts for Simple and Complex Liquids, Cambridge University Press (2003).
Como trabajo futuro se puede desarrollar la soluci o´ n de la ecuacio´ n de Poisson-Boltzmann sobre dominios no trivia-
[4] McQuarrie, D. Statistical Mechanics. 2nd ed. University Science Books, 2000. ISBN: 9781891389153.
240
