Resumen de Ewald

La suma de Ewald , que recibe su nombre de Paul Peter Ewald , es un método para calcular interacciones de largo alcance (por ejemplo, interacciones electrostáticas ) en sistemas periódicos. Se desarrolló inicialmente para calcular las energías electrostáticas de cristales iónicos y actualmente se utiliza 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 por 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 sin singularidad . La contribución de corto alcance se calcula en el espacio real, mientras que la de largo alcance se calcula mediante una transformada de Fourier . La ventaja de este método radica en la rápida convergencia de la energía en comparación con la suma directa. Esto significa que el método posee alta precisión y una velocidad razonable al calcular interacciones de largo alcance, convirtiéndose así en el método estándar de facto para el cálculo de 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 presentan un estudio sobre los errores de truncamiento introducidos en los cálculos de energía y fuerza de sistemas de carga puntual desordenados. ^{[ 1 ]}

La suma de Ewald reescribe el potencial de interacción como la suma de dos términos, $${\displaystyle \varphi (\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \varphi _{sr}(\mathbf {r} )+\varphi _{\ell r}(\mathbf {r} ),}$$ dónde${\displaystyle \varphi _{sr}(\mathbf {r} )}$representa el término de corto alcance cuya suma converge rápidamente en el espacio real y${\displaystyle \varphi _{\ell r}(\mathbf {r} )}$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, más típicamente 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 razonable para los interiores de los cristales). Una unidad repetitiva de este sistema periódico hipotético se llama celda unitaria . Una de estas celdas se elige como la "celda central" para referencia y las celdas restantes se llaman imágenes .

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, puede representarse como una integral doble sobre dos campos de densidad de carga que representan los campos de la celda unitaria y la red cristalina. $${\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 })}$$ donde el campo de densidad de carga de la celda unitaria${\displaystyle \rho _ {uc}(\mathbf {r} )}$es una suma sobre las posiciones${\displaystyle \mathbf {r} _{k}}$de los cargos${\displaystyle q_{k}}$en la celda unitaria central $${\displaystyle \rho _{uc}(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{\mathrm {cargas} \ k}q_{k}\delta (\mathbf {r} -\mathbf {r} _{k})}$$ y el campo de densidad de carga total${\displaystyle \rho _{\text{TOT}}(\mathbf {r} )}$es la misma suma sobre las cargas de la celda unitaria${\displaystyle q_{k}}$y sus imágenes periódicas $${\displaystyle \rho _{\text{TOT}}(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{n_{1},n_{2},n_{3}}\sum _{\mathrm {cargas} \ 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í,${\displaystyle \delta (\mathbf {x} )}$es la función delta de Dirac ,${\displaystyle \mathbf {a} _{1}}$,${\displaystyle \mathbf {a} _{2}}$y${\displaystyle \mathbf {a} _{3}}$son los vectores de la red y${\displaystyle n_{1}}$,${\displaystyle n_{2}}$y${\displaystyle n_{3}}$abarca todos los números enteros. El campo total${\displaystyle \rho _{\text{TOT}}(\mathbf {r} )}$puede representarse como una convolución de${\displaystyle \rho _ {uc}(\mathbf {r} )}$: $${\displaystyle \rho _{\text{TOT}}(\mathbf {r} )=\int d\mathbf {r'} \rho _{uc}(\mathbf {r'} )L(\mathbf {r} -\mathbf {r'} ),}$$ con una función de red${\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})}$$

Dado que se trata de una convolución , la transformada de Fourier de${\displaystyle \rho _{\text{TOT}}(\mathbf {r} )}$es un producto $${\displaystyle {\tilde {\rho }}_{\text{TOT}}(\mathbf {k} )={\tilde {L}}(\mathbf {k} ){\tilde {\rho }}_{uc}(\mathbf {k} )}$$ donde la transformada de Fourier de la función reticular es otra suma sobre funciones delta $${\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})}$$donde se definen los vectores del espacio recíproco${\displaystyle \mathbf {b} _{1}\ {\stackrel {\mathrm {def} }{=}}\ 2\pi {\frac {\mathbf {a} _{2}\times \mathbf {a} _{3}}{\Omega }}}$(y permutaciones cíclicas) donde${\displaystyle \Omega \ {\stackrel {\mathrm {def} }{=}}\ \mathbf {a} _{1}\cdot \left(\mathbf {a} _{2}\times \mathbf {a} _{3}\right)}$es el volumen de la celda unitaria central (si geométricamente es un paralelepípedo , lo cual suele ser así, pero no necesariamente el caso). Nótese que ambos${\displaystyle L(\mathbf {r} )}$y${\displaystyle {\tilde {L}}(\mathbf {k} )}$son reales, incluso funciones.

Para abreviar, definamos un potencial efectivo de partícula única. $${\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 $${\displaystyle {\tilde {V}}(\mathbf {k} )\ {\stackrel {\mathrm {def} }{=}}\ {\tilde {\rho }}_{uc}(\mathbf {k} ){\tilde {\Phi }}(\mathbf {k} )}$$ donde se define la transformada de Fourier $${\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} )}$$

dónde${\displaystyle \mathbf {k} =m_{1}\mathbf {b} _{1}+m_{2}\mathbf {b} _{2}+m_{3}\mathbf {b} _{3}}$en el resumen final.

Este es el resultado esencial. Una vez${\displaystyle {\tilde {\rho }}_{uc}(\mathbf {k} )}$se calcula la suma/integración sobre${\displaystyle \mathbf {k} }$Es sencillo 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 eléctricamente neutra para evitar sumas infinitas.

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 tenido 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 inversa al cuadrado, como la gravedad o la electrostática . Recientemente, el PME también se ha utilizado para calcular la${\displaystyle r^{-6}}$parte del potencial de Lennard-Jones para eliminar artefactos debidos a la truncación. ^{[ 2 ]} Las aplicaciones incluyen simulaciones de plasmas , galaxias y moléculas .

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.${\displaystyle \varphi (\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \varphi _{sr}(\mathbf {r} )+\varphi _{\ell r}(\mathbf {r} )}$La idea básica de la suma de Ewald de malla de partículas es reemplazar la suma directa de energías de interacción entre partículas puntuales. $${\displaystyle E_{\text{TOT}}=\sum _{i,j}\varphi (\mathbf {r} _{j}-\mathbf {r} _{i})=E_{sr}+E_{\ell r}}$$ con dos sumas, una suma directa${\displaystyle E_{sr}}$del potencial de corto alcance en el espacio real $${\displaystyle E_{sr}=\sum _{i,j}\varphi _{sr}(\mathbf {r} _{j}-\mathbf {r} _{i})}$$ (esta es la parte de partículas de la malla de partículas de Ewald ) y una suma en el espacio de Fourier de la parte de largo alcance $${\displaystyle E_{\ell r}=\sum _{\mathbf {k} }{\tilde {\Phi }}_{\ell r}(\mathbf {k} )\left|{\tilde {\rho }}(\mathbf {k} )\right|^{2}}$$

dónde${\displaystyle {\tilde {\Phi }}_{\ell r}}$y${\displaystyle {\tilde {\rho }}(\mathbf {k} )}$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), pueden truncarse con poca pérdida de precisión y una gran mejora en el tiempo de cálculo requerido. Para evaluar la transformada de Fourier${\displaystyle {\tilde {\rho }}(\mathbf {k} )}$Para medir eficientemente el campo de densidad de carga, 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 la malla ).

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 de extensión espacial infinita. En las simulaciones de dinámica molecular , esto se logra normalmente mediante la construcción deliberada de una celda unitaria eléctricamente neutra que puede ser "teselada" infinitamente para formar imágenes; sin embargo, para tener en cuenta adecuadamente los efectos de esta aproximación, estas imágenes se reincorporan 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 izquierda y la frontal con la posterior. Como resultado, el tamaño de la celda unitaria debe elegirse cuidadosamente para que sea lo suficientemente grande como para evitar correlaciones de movimiento impropias entre dos caras "en contacto", pero aún lo suficientemente pequeño como para que sea computacionalmente factible. La definición del límite entre las 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 potenciales 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.

La energía electrostática de un cristal polar (es decir, un cristal con un dipolo neto)${\displaystyle \mathbf {p} _{uc}}$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${\displaystyle R}$crece como${\textstyle R^{2}}$; (2) la fuerza de una sola interacción dipolo-dipolo disminuye como${\textstyle 1/{R^{3}}}$; y (3) la suma matemática${\textstyle \sum _{n=1}^{\infty }{\frac {1}{n}}}$diverge.

Este resultado, algo sorprendente, puede conciliarse con la energía finita de los cristales reales, ya que 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.${\displaystyle \sigma =\mathbf {P} \cdot \mathbf {n} }$dónde${\displaystyle \mathbf {n} }$es el vector normal de la superficie y${\displaystyle \mathbf {P} }$representa el momento dipolar neto por volumen. La energía de interacción${\displaystyle U}$del dipolo en una celda unitaria central con esa densidad de carga superficial se puede escribir ^{[ 3 ]}$${\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}$$ dónde${\displaystyle \mathbf {p} _{uc}}$y${\displaystyle V_{uc}}$son el momento dipolar neto y el volumen de la celda unitaria,${\displaystyle dS}$es un área infinitesimal en la superficie del cristal y${\displaystyle \mathbf {r} }$es el vector desde la celda unitaria central hasta el área infinitesimal. Esta fórmula resulta de la integración de la energía${\displaystyle dU=-\mathbf {p} _{uc}\cdot d\mathbf {E} }$dónde${\displaystyle d\mathbf {E} }$representa el campo eléctrico infinitesimal generado por una carga superficial infinitesimal${\displaystyle dq\ {\stackrel {\mathrm {def} }{=}}\ \sigma dS}$( Ley de Coulomb ) $${\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}}}}$$ El signo negativo deriva de la definición de${\displaystyle \mathbf {r} }$, que apunta hacia la carga, no en dirección opuesta a ella.

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

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

• Pablo Peter Ewald
• Thomas A. Darden
• Constante de Madelung
• fórmula de suma de Poisson
• Modelado molecular
• Resumen de Wolf