Suma de Ewald

La suma de Ewald , llamada así por Paul Peter Ewald , es un método para calcular interacciones de largo alcance (por ejemplo, interacciones electrostáticas ) en sistemas periódicos. Se desarrolló primero como método para calcular las energías electrostáticas de cristales iónicos y ahora se usa comúnmente para calcular interacciones de largo alcance en química computacional . La suma de Ewald es un caso especial de la fórmula de suma de Poisson , que reemplaza la suma de energías de interacción en el espacio real con una suma equivalente en el espacio de Fourier . En este método, la interacción de largo alcance se divide en dos partes: una contribución de corto alcance y una contribución de largo alcance que no tiene una singularidad . La contribución de corto alcance se calcula en el espacio real, mientras que la contribución de largo alcance se calcula utilizando una transformada de Fourier . La ventaja de este método es la rápida convergencia de la energía en comparación con la de una suma directa. Esto significa que el método tiene una alta precisión y una velocidad razonable al calcular interacciones de largo alcance, y por lo tanto es el método estándar de facto para calcular interacciones de largo alcance en sistemas periódicos. El método requiere neutralidad de carga del sistema molecular para calcular con precisión la interacción coulombiana total. Kolafa y Perram proporcionan un estudio de los errores de truncamiento introducidos en los cálculos de energía y fuerza de sistemas de carga puntual desordenados. ^[1]

Derivación

La suma de Ewald reescribe el potencial de interacción como la suma de dos términos, donde representa el término de corto alcance cuya suma converge rápidamente en el espacio real y representa el término de largo alcance cuya suma converge rápidamente en el espacio de Fourier (recíproco). La parte de largo alcance debe ser finita para todos los argumentos (en particular r  = 0), pero puede tener cualquier forma matemática conveniente, la más típica es una distribución gaussiana . El método supone que la parte de corto alcance se puede sumar fácilmente; por lo tanto, el problema se convierte en la suma del término de largo alcance. Debido al uso de la suma de Fourier, el método supone implícitamente que el sistema en estudio es infinitamente periódico (una suposición sensata para los interiores de los cristales). Una unidad repetitiva de este sistema periódico hipotético se llama celda unitaria . Una de esas celdas se elige como la "celda central" para referencia y las celdas restantes se denominan imágenes . $${\displaystyle \varphi (\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \varphi _{sr}(\mathbf {r} )+\varphi _{\ell r}(\mathbf {r} ),}$$ ${\displaystyle \varphi _{sr}(\mathbf {r} )}$ ${\displaystyle \varphi _{\ell r}(\mathbf {r} )}$

La energía de interacción de largo alcance es la suma de las energías de interacción entre las cargas de una celda unitaria central y todas las cargas de la red cristalina. Por lo tanto, se puede representar como una doble integral sobre dos campos de densidad de carga que representan los campos de la celda unitaria y la red cristalina, donde el campo de densidad de carga de la celda unitaria es una suma sobre las posiciones de las cargas en la celda unitaria central y el campo de densidad de carga total es la misma suma sobre las cargas de la celda unitaria y sus imágenes periódicas. $${\displaystyle E_{\ell r}=\iint d\mathbf {r} \,d\mathbf {r} ^{\prime }\,\rho _{\text{TOT}}(\mathbf {r} )\rho _{uc}(\mathbf {r} ^{\prime })\ \varphi _{\ell r}(\mathbf {r} -\mathbf {r} ^{\prime })}$$ ${\displaystyle \rho _{uc}(\mathbf {r} )}$ ${\displaystyle \mathbf {r} _{k}}$ ${\displaystyle q_{k}}$ $${\displaystyle \rho _{uc}(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{\mathrm {charges} \ k}q_{k}\delta (\mathbf {r} -\mathbf {r} _{k})}$$ ${\displaystyle \rho _{\text{TOT}}(\mathbf {r} )}$ ${\displaystyle q_{k}}$ $${\displaystyle \rho _{\text{TOT}}(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{n_{1},n_{2},n_{3}}\sum _{\mathrm {charges} \ k}q_{k}\delta (\mathbf {r} -\mathbf {r} _{k}-n_{1}\mathbf {a} _{1}-n_{2}\mathbf {a} _{2}-n_{3}\mathbf {a} _{3})}$$

Aquí, es la función delta de Dirac , , y son los vectores de red y , y abarcan todos los números enteros. El campo total se puede representar como una convolución de con una función de red ${\displaystyle \delta (\mathbf {x} )}$ ${\displaystyle \mathbf {a} _{1}}$ ${\displaystyle \mathbf {a} _{2}}$ ${\displaystyle \mathbf {a} _{3}}$ ${\displaystyle n_{1}}$ ${\displaystyle n_{2}}$ ${\displaystyle n_{3}}$ ${\displaystyle \rho _{\text{TOT}}(\mathbf {r} )}$ ${\displaystyle \rho _{uc}(\mathbf {r} )}$ ${\displaystyle L(\mathbf {r} )}$ $${\displaystyle L(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{n_{1},n_{2},n_{3}}\delta (\mathbf {r} -n_{1}\mathbf {a} _{1}-n_{2}\mathbf {a} _{2}-n_{3}\mathbf {a} _{3})}$$

Como se trata de una convolución , la transformación de Fourier de es un producto donde la transformada de Fourier de la función reticular es otra suma sobre funciones delta donde se definen los vectores espaciales recíprocos (y permutaciones cíclicas) donde es el volumen de la celda unitaria central (si es geométricamente un paralelepípedo , lo que suele ser el caso, pero no necesariamente). Nótese que tanto y son funciones reales, pares. ${\displaystyle \rho _{\text{TOT}}(\mathbf {r} )}$ $${\displaystyle {\tilde {\rho }}_{\text{TOT}}(\mathbf {k} )={\tilde {L}}(\mathbf {k} ){\tilde {\rho }}_{uc}(\mathbf {k} )}$$ $${\displaystyle {\tilde {L}}(\mathbf {k} )={\frac {\left(2\pi \right)^{3}}{\Omega }}\sum _{m_{1},m_{2},m_{3}}\delta (\mathbf {k} -m_{1}\mathbf {b} _{1}-m_{2}\mathbf {b} _{2}-m_{3}\mathbf {b} _{3})}$$ ${\displaystyle \mathbf {b} _{1}\ {\stackrel {\mathrm {def} }{=}}\ 2\pi {\frac {\mathbf {a} _{2}\times \mathbf {a} _{3}}{\Omega }}}$ ${\displaystyle \Omega \ {\stackrel {\mathrm {def} }{=}}\ \mathbf {a} _{1}\cdot \left(\mathbf {a} _{2}\times \mathbf {a} _{3}\right)}$ ${\displaystyle L(\mathbf {r} )}$ ${\displaystyle {\tilde {L}}(\mathbf {k} )}$

Para abreviar, defina un potencial efectivo de una sola partícula. $${\displaystyle v(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \int d\mathbf {r} ^{\prime }\,\rho _{uc}(\mathbf {r} ^{\prime })\ \varphi _{\ell r}(\mathbf {r} -\mathbf {r} ^{\prime })}$$

Dado que esto también es una convolución, la transformada de Fourier de la misma ecuación es un producto donde la transformada de Fourier está definida $${\displaystyle {\tilde {V}}(\mathbf {k} )\ {\stackrel {\mathrm {def} }{=}}\ {\tilde {\rho }}_{uc}(\mathbf {k} ){\tilde {\Phi }}(\mathbf {k} )}$$ $${\displaystyle {\tilde {V}}(\mathbf {k} )=\int d\mathbf {r} \ v(\mathbf {r} )\ e^{-i\mathbf {k} \cdot \mathbf {r} }}$$

La energía ahora se puede escribir como una integral de campo único . $${\displaystyle E_{\ell r}=\int d\mathbf {r} \ \rho _{\text{TOT}}(\mathbf {r} )\ v(\mathbf {r} )}$$

Utilizando el teorema de Plancherel , la energía también se puede sumar en el espacio de Fourier. $${\displaystyle E_{\ell r}=\int {\frac {d\mathbf {k} }{\left(2\pi \right)^{3}}}\ {\tilde {\rho }}_{\text{TOT}}^{*}(\mathbf {k} ){\tilde {V}}(\mathbf {k} )=\int {\frac {d\mathbf {k} }{\left(2\pi \right)^{3}}}{\tilde {L}}^{*}(\mathbf {k} )\left|{\tilde {\rho }}_{uc}(\mathbf {k} )\right|^{2}{\tilde {\Phi }}(\mathbf {k} )={\frac {1}{\Omega }}\sum _{m_{1},m_{2},m_{3}}\left|{\tilde {\rho }}_{uc}(\mathbf {k} )\right|^{2}{\tilde {\Phi }}(\mathbf {k} )}$$

donde en la suma final. ${\displaystyle \mathbf {k} =m_{1}\mathbf {b} _{1}+m_{2}\mathbf {b} _{2}+m_{3}\mathbf {b} _{3}}$

Este es el resultado esencial. Una vez calculado, la suma/integración es sencilla y debería converger rápidamente. La razón más común de la falta de convergencia es una celda unitaria mal definida, que debe ser neutral en cuanto a carga para evitar sumas infinitas. ${\displaystyle {\tilde {\rho }}_{uc}(\mathbf {k} )}$ ${\displaystyle \mathbf {k} }$

Método de Ewald de malla de partículas (PME)

La suma de Ewald se desarrolló como un método en física teórica , mucho antes de la llegada de las computadoras . Sin embargo, el método de Ewald ha disfrutado de un uso generalizado desde la década de 1970 en simulaciones por computadora de sistemas de partículas, especialmente aquellos cuyas partículas interactúan a través de una ley de fuerza del cuadrado inverso, como la gravedad o la electrostática . Recientemente, PME también se ha utilizado para calcular la parte del potencial de Lennard-Jones con el fin de eliminar artefactos debido al truncamiento. ^[2] Las aplicaciones incluyen simulaciones de plasmas , galaxias y moléculas . ${\displaystyle r^{-6}}$

En el método de malla de partículas, al igual que en la suma de Ewald estándar, el potencial de interacción genérico se separa en dos términos . La idea básica de la suma de Ewald de malla de partículas es reemplazar la suma directa de las energías de interacción entre partículas puntuales con dos sumas, una suma directa del potencial de corto alcance en el espacio real (esta es la parte de partículas de la suma de Ewald de malla de partículas ) y una suma en el espacio de Fourier de la parte de largo alcance. ${\displaystyle \varphi (\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \varphi _{sr}(\mathbf {r} )+\varphi _{\ell r}(\mathbf {r} )}$ $${\displaystyle E_{\text{TOT}}=\sum _{i,j}\varphi (\mathbf {r} _{j}-\mathbf {r} _{i})=E_{sr}+E_{\ell r}}$$ ${\displaystyle E_{sr}}$ $${\displaystyle E_{sr}=\sum _{i,j}\varphi _{sr}(\mathbf {r} _{j}-\mathbf {r} _{i})}$$ $${\displaystyle E_{\ell r}=\sum _{\mathbf {k} }{\tilde {\Phi }}_{\ell r}(\mathbf {k} )\left|{\tilde {\rho }}(\mathbf {k} )\right|^{2}}$$

donde y representan las transformadas de Fourier del potencial y la densidad de carga (esta es la parte de Ewald ). Dado que ambas sumas convergen rápidamente en sus respectivos espacios (real y de Fourier), se pueden truncar con poca pérdida de precisión y una gran mejora en el tiempo computacional requerido. Para evaluar la transformada de Fourier del campo de densidad de carga de manera eficiente, se utiliza la transformada rápida de Fourier , que requiere que el campo de densidad se evalúe en una red discreta en el espacio (esta es la parte de malla ). ${\displaystyle {\tilde {\Phi }}_{\ell r}}$ ${\displaystyle {\tilde {\rho }}(\mathbf {k} )}$ ${\displaystyle {\tilde {\rho }}(\mathbf {k} )}$

Debido a la suposición de periodicidad implícita en la suma de Ewald, las aplicaciones del método PME a sistemas físicos requieren la imposición de simetría periódica. Por lo tanto, el método es más adecuado para sistemas que pueden simularse como infinitos en extensión espacial. En simulaciones de dinámica molecular , esto normalmente se logra construyendo deliberadamente una celda unitaria de carga neutra que puede "mosaicarse" infinitamente para formar imágenes; sin embargo, para tener en cuenta adecuadamente los efectos de esta aproximación, estas imágenes se reincorporan nuevamente a la celda de simulación original. El efecto general se denomina condición de contorno periódica . Para visualizar esto con mayor claridad, piense en un cubo unitario; la cara superior está efectivamente en contacto con la cara inferior, la derecha con la cara izquierda y la parte delantera con la cara trasera. Como resultado, el tamaño de la celda unitaria debe elegirse cuidadosamente para que sea lo suficientemente grande como para evitar correlaciones de movimiento inapropiadas entre dos caras "en contacto", pero lo suficientemente pequeño como para que sea computacionalmente factible. La definición del límite entre interacciones de corto y largo alcance también puede introducir artefactos.

La restricción del campo de densidad a una malla hace que el método PME sea más eficiente para sistemas con variaciones "suaves" de densidad o funciones de potencial continuas. Los sistemas localizados o aquellos con grandes fluctuaciones de densidad pueden tratarse de manera más eficiente con el método multipolar rápido de Greengard y Rokhlin.

Término dipolar

La energía electrostática de un cristal polar (es decir, un cristal con un dipolo neto en la celda unitaria) es condicionalmente convergente , es decir, depende del orden de la suma. Por ejemplo, si las interacciones dipolo-dipolo de una celda unitaria central con celdas unitarias ubicadas en un cubo cada vez mayor, la energía converge a un valor diferente que si las energías de interacción se hubieran sumado esféricamente. En términos generales, esta convergencia condicional surge porque (1) el número de dipolos interactuantes en una capa de radio crece como ; (2) la fuerza de una única interacción dipolo-dipolo cae como ; y (3) la suma matemática diverge. ${\displaystyle \mathbf {p} _{uc}}$ ${\displaystyle R}$ ${\textstyle R^{2}}$ ${\textstyle 1/{R^{3}}}$ ${\textstyle \sum _{n=1}^{\infty }{\frac {1}{n}}}$

Este resultado algo sorprendente se puede conciliar con la energía finita de los cristales reales porque dichos cristales no son infinitos, es decir, tienen un límite particular. Más específicamente, el límite de un cristal polar tiene una densidad de carga superficial efectiva en su superficie donde es el vector normal a la superficie y representa el momento dipolar neto por volumen. La energía de interacción del dipolo en una celda unitaria central con esa densidad de carga superficial se puede escribir ^[3] donde y son el momento dipolar neto y el volumen de la celda unitaria, es un área infinitesimal en la superficie del cristal y es el vector desde la celda unitaria central hasta el área infinitesimal. Esta fórmula resulta de integrar la energía donde representa el campo eléctrico infinitesimal generado por una carga superficial infinitesimal ( ley de Coulomb ) El signo negativo deriva de la definición de , que apunta hacia la carga, no lejos de ella. ${\displaystyle \sigma =\mathbf {P} \cdot \mathbf {n} }$ ${\displaystyle \mathbf {n} }$ ${\displaystyle \mathbf {P} }$ ${\displaystyle U}$ $${\displaystyle U={\frac {1}{2V_{uc}}}\int {\frac {\left(\mathbf {p} _{uc}\cdot \mathbf {r} \right)\left(\mathbf {p} _{uc}\cdot \mathbf {n} \right)}{r^{3}}}\,dS}$$ ${\displaystyle \mathbf {p} _{uc}}$ ${\displaystyle V_{uc}}$ ${\displaystyle dS}$ ${\displaystyle \mathbf {r} }$ ${\displaystyle dU=-\mathbf {p} _{uc}\cdot d\mathbf {E} }$ ${\displaystyle d\mathbf {E} }$ ${\displaystyle dq\ {\stackrel {\mathrm {def} }{=}}\ \sigma dS}$ $${\displaystyle d\mathbf {E} \ {\stackrel {\mathrm {def} }{=}}\ \left({\frac {-1}{4\pi \epsilon }}\right){\frac {dq\ \mathbf {r} }{r^{3}}}=\left({\frac {-1}{4\pi \epsilon }}\right){\frac {\sigma \,dS\ \mathbf {r} }{r^{3}}}}$$ ${\displaystyle \mathbf {r} }$

Historia

La suma de Ewald fue desarrollada por Paul Peter Ewald en 1921 (ver referencias a continuación) para determinar la energía electrostática (y, por lo tanto, la constante de Madelung ) de los cristales iónicos.

Escalada

En general, los distintos métodos de suma de Ewald dan diferentes complejidades temporales . El cálculo directo da , donde es el número de átomos en el sistema. El método PME da . ^[4] ${\displaystyle O(N^{2})}$ ${\displaystyle N}$ ${\displaystyle O(N\,\log N)}$

Véase también

• Pablo Peter Ewald
• Constante de Madelung
• Fórmula de suma de Poisson
• Modelado molecular
• Suma de lobos

Referencias

1. Kolafa, Jiri; Perram, John W. (septiembre de 1992). "Errores de corte en las fórmulas de suma de Ewald para sistemas de carga puntual". Molecular Simulation . 9 (5): 351–368. doi :10.1080/08927029208049126.
2. Di Pierro, M.; Elber, R.; Leimkuhler, B. (2015), "Un algoritmo estocástico para el conjunto isobárico-isotérmico con sumas de Ewald para todas las fuerzas de largo alcance.", Journal of Chemical Theory and Computation , 11 (12): 5624–5637, doi :10.1021/acs.jctc.5b00648, PMC 4890727 , PMID  26616351 
3. Herce, HD; Garcia, AE; Darden, T (28 de marzo de 2007). "El término de superficie electrostática: (I) sistemas periódicos". The Journal of Chemical Physics . 126 (12): 124106. Bibcode :2007JChPh.126l4106H. doi :10.1063/1.2714527. PMID  17411107.
4. Darden, Tom; York, Darrin; Pedersen, Lee (15 de junio de 1993). "Ewald de malla de partículas: un método N ⋅log( N ) para sumas de Ewald en sistemas grandes". The Journal of Chemical Physics . 98 (12): 10089–10092. Bibcode :1993JChPh..9810089D. doi :10.1063/1.464397. ISSN  0021-9606.

• Ewald, P (1921). "Die Berechnung optischer und elektrostatischer Gitterpotentiale". Ana. Física . 369 (3): 253–287. Código bibliográfico : 1921AnP...369..253E. doi : 10.1002/andp.19213690304.
• Darden, T; Perera, L; Li, L; Pedersen, L (1999). "Nuevos trucos para modeladores del conjunto de herramientas de cristalografía: el algoritmo de malla de partículas de Ewald y su uso en simulaciones de ácidos nucleicos". Estructura . 7 (3): R55–R60. doi : 10.1016/S0969-2126(99)80033-1 . PMID  10368306. S2CID  40964921.
• Frenkel, D., y Smit, B. (2001). Comprensión de la simulación molecular: de los algoritmos a las aplicaciones , Academic press.