Comparison of results of the Ewald Particle Grid method with the standard Ewald method.
Project proposal PAPIIT, DGAPA, or CONACYT 2003.
MS. Roberto Pérez
INTRODUCTION
In a typical MD simulation, about 95% of the CPU time, it is spent examining the complete set of N (N-1) / 2 pairs, identifying these pairs separated by at least the cutting distance rc, and computing the internal forces for this subset. The interactions of remaining pairs do not influence the dynamics of the system. For each molecule the set of neighbors within the cutting distance changes over time; The problem of identifying and rejecting the others is CPU time consumed. Any other increase in speed can reduce the time spent currently evaluating peer forces within the cutting radius.
Multi-step temporary method
(MULTIPLE TIME STEP METHOD)
The multi-step method (multi-step method), was designed for this, the neighborhood of an atom i is split into 2 groups, which are the primary neighbors, which rest at a distance rp from the atom i, and the secondary neighbors , which are at a distance between rp and rc of i. In this way, the total fixed force at i is separated into a primary fip component, due to the closest neighbors and a second fis component due to the furthest neighbors. Typically for a Lennard Jones fluid, rp will be chosen to be between and 1.5.
Where is the length parameter in a potential pair.
Comparación de resultados de
el método de Rejilla de Partículas de Ewald con el método de Ewald estándar.
Propuesta
de proyecto PAPIIT, DGAPA, ó CONACYT 2003.
MS.
Roberto Pérez
INTRODUCCIÓN
En una simulación típica de
MD, cerca del 95% del tiempo de CPU, se gasta en examinar el conjunto completo
de N(N-1)/2 pares, identificando estos pares separados por lo menos la
distancia de corte rc, y
computando las fuerzas internas para este subconjunto. Las interacciones de
pares remanentes no influyen en la dinámica del sistema. Para cada molécula el conjunto
de vecinos dentro de la distancia de corte cambia con el tiempo; el problema de
identificar y rechazar las otras es tiempo de CPU consumido. Cualquier otro
incremento en la velocidad puede alcanzarse reduciendo el tiempo gastado
actualmente evaluando fuerzas entre pares dentro del radio de corte.
Método
de múltiples pasos temporales
(MULTIPLE TIME STEP METHOD)
El método de múltiples pasos (multiple time step method), fue diseñado
para esto, la vecindad de un átomo i es partida en 2 grupos, los cuales son
los vecinos primarios, que descansan a una distancia rp del átomo i,
y los vecinos secundarios, los cuales se encuentran a una distancia entre rp y rc de i. De
esta manera, la fuerza total fi
en i es separada en una componente
primaria fip, debida a los
vecinos más cercanos y una segunda
componente fis debido a
los vecinos más lejanos. Típicamente para un fluido de Lennard Jones, rp se escogerá para
encontrarse entre s y 1.5s .
Donde s es el parámetro
de longitud en un potencial par.
El movimiento del átomo i es dominado por un cambio rápido de la
fuerza primaria resultado de las colisiones con el vecino más cercano. La
fuerza secundaria es más pequeña y cambia lentamente con el tiempo. El método de múltiples pasos obtiene
ventajas de esta separación, y tiene el fin de calcular las interacciones pares
secundarias de manera menos frecuente que las primarias que son más importantes.
El método es
mejor si se usa conjuntamente con el predictor de Grear el cual es un algoritmo
corrector de distancias. Al tiempo t
las fuerzas primarias y secundarias en cada átomo i son calculadas de manera usual. Al mismo tiempo se calculan las
derivadas temporales de las fuerza secundarias,
, …
hasta el orden m, son calculadas, y una lista de
vecinos primarios es compilada. A cada
uno del siguiente tmts-1 pasos, la fuerza primaria
es calculada explícitamente. La fuerza secundaria es estimada usando la serie
de Taylor de orden m.
Siguiendo estos tmts pasos (uno que
envuelve el cálculo completo de fuerzas y el segundo las derivadas temporales,
y tmts-1 usando la ecuación anterior para las fuerzas
secundarias), el proceso es repetido íntegramente, empezando con un nuevo
cálculo de la lista primaria. Entonces el método usa dos pasos temporales, dt para las
interacciones primarias y tmts dt para las secundarias.
Expresiones convenientes para
las derivadas de fuerza pueden obtenerse en términos de las derivadas del
tiempo de las posiciones de las partículas y las derivadas especiales del
potencial. Por ejemplo,
Y las sumas
son sobre todos los vecinos secundarios j.
Estos cálculos toman la forma en el loop
de fuerza; derivadas temporales como
Son calculadas de las
derivadas temporales de posiciones atómicas disponibles para el predictor de
escenario de algoritmos de Gear.
En el estudio de Streett et al, una serie de tercer orden de
Taylor, con tmts =10, rc
= 2.5 y rp = 1.1 fueron usados.
En la simulación de un fluido de Lennard Jones a una densidad reducida r
=
0.8, los números promedio de vecinos primarios y secundarios son del orden de
1.4 y 24, respectivamente. El método de
múltiples tiempos en este caso resultó en una mejora de velocidad de un factor
entre 7 y 10 sobre un programa convencional de MD que no tenían listas de vecinos y un factor
de 3 a 5 sobre un programa usando la lista de vecinos de Verlet. Estos tiempos
relativos dependerán de detalles precisos de operación llevados a cabo en el
loop de fuerza, y será dependiente del tipo de máquina.
El método ha sido felizmente
aplicado a fluidos moleculares [Swindoll y Haile 1984]. Se ha extendido al caso
cuando hay 3 regiones dentro de la esfera de radio rc [Nicolas et al].
Una aplicación interesante de este método es en la simulación de fluidos
con fuerzas de tres cuerpos [Haile 1978]. Para el potencial de Axilrod-Teller
la contribución de tres cuerpos a la fuerza varía lentamente con el tiempo, y
es posible tratar la componente entera de tres cuerpos como una interacción
secundaria. De esta manera la sumatoria prohibitiva y cara sobre el triplete
puede ser evitada, para la mayoría de los pasos en la simulación. Todavía
existen dificultades para este método.
Primero.- Hay un esfuerzo de programación considerable
requerido para incorporar el código relevante en el programa convencional.
Segundo.- El método ha sido usado
solamente con un tiempo dt más pequeño que el usado en
simulaciones convencionales (5 X 10 -15
s) en oposición al más usual 10 -14.
Se encuentra que usar un tiempo más grande requiere un incremento en rp y un decremento en el
intervalo tmts, entonces se compensaran algunas de las ventajas
ganadas usando el método de múltiples pasos. La eficiencia de este método es
claramente sensitiva a la elección de rp,
sin embargo, la mejor manera de probar este método es fijar el tiempo deseado,
decidir sobre un nivel de energía satisfactorio, y ajustar rp y el intervalo tmts, para maximizar la
eficiencia de la simulación sujetas a estas variables. Para fluidos
moleculares, donde un tiempo menor se requiere, estas dificultades se hacen
menos evidentes.
COMO
LIDIAR CON FUERZAS DE LARGO ALCANCE
En las secciones previas, se
ha discutido el núcleo del programa cuando las fuerzas involucradas son de
corto alcance. En esta sección cambiamos nuestra atención para lidiar con
fuerzas de largo alcance. Una fuerza de largo alcance es frecuentemente
definida como una fuerza en la cual la interacción espacial decae no más rápido
que r-d donde d es la dimensionalidad del sistema.
En esta categoría se encuentran las
interacciones carga-carga entre los
iones (Vzz(r)~ r -1) y la interacción
dipolo-dipolo, la cual es básicamente una interacción molecular (Vmm(r)~ r -3). Estas fuerzas son un serio
problema para las simulaciones computacionales, debido a que su rango es mayor
que la mitad de la longitud de la caja para las simulaciones típicas de 500
moléculas. La solución por medio de fuerza bruta para este problema
podría incrementar el tamaño de la caja central L en cientos de nanómetros, lo
cual generaría que el apantallamiento de los vecinos disminuya el rango
efectivo de los potenciales.
Incluso con computadoras
modernas esta solución es impractica, debido a que el tiempo requerido para
ejecutar dicha simulación es aproximadamente proporcional a N2, i.e.
L6.
¿Cómo
solucionar este tipo de problemas?
El problema es particularmente
agudo para Vzz(r).De esta
manera el corte esférico del potencial puede desecharse. La esfera resultante
alrededor de un ión dado puede desecharse. La esfera resultante alrededor de un
ión dado puede cargarse, debido a que el número de aniones y cationes no
necesitan balancearse a cualquier instante. La tendencia de los iones a migrar hacia atrás
y adelante sobre la superficie esférica creará efectos artificiales en r = rc. Esto puede contabilizarse
distribuyendo una nueva carga sobre la superficie de la esfera, igual en
magnitud y opuesta en signo a la carga neta de la esfera, y esto se efectúa
para garantizar la electroneutralidad local. Esto es como cambiar el potencial,
Adams [1983b] mostró que los resultados de esta aproximación son sistemas
dependientes del tamaño, pero para sistemas de 512 iones los resultados están
en buen acuerdo con los resultados obtenidos de las sumas de Ewald.
En contraste, el método básico
de mínima imagen correspondiente a un corte de potencial en la superficie del
cubo rodeado del ión en cuestión. Este cubo será eléctricamente neutro. De
cualquier manera, la desventaja es que iones cargados similares tenderán a
ocupar posiciones en las esquinas opuestas del cubo: la estructura de imágenes
periódicas será impuesta directamente en lo que debería de ser un líquido
isotrópico, y esto tiene como resultado una distorsión en la estructura del
líquido. Este efecto podría reducirse en las cajas no cúbicas. De manera
similar, y con efectos menos importantes se obtiene al efectuar el corte
esférico o de mínima imagen a un sistema polar (Vmm(r)).
En esta trabajo nos
concentraremos en 2 métodos que pueden usarse para resolver el problema de las
fuerzas de largo alcance. El método de las sumas de Ewald, incluye la interacción
de un ión o molécula con todas sus imágenes periódicas. Este método tenderá a
enfatizar la naturaleza periódica del fluido modelo.
MÉTODO
DE PARTICLE MESH EWALD(PPME) vs. Método Estándar de Ewald[8]
Este método tenderá a
enfatizar la naturaleza continua de un
fluido polar y requiere una estimación a priori
de la permitividad relativa. Ambos métodos usan las ideas bien conocida de la
teoría electrostática. En particular, una distribución de carga con una cavidad
esférica polariza el medio que lo rodea. Esta polarización que depende de la
permitividad relativa del medio, tiene un efecto en la distribución de carga en
la cavidad.
En la siguiente parte del
desarrollo de este proyecto se explica como lidiar con las fuerzas
electrostáticas de largo alcance y cómo combinar el método de SPME con el
multinivel de integración de las ecuaciones de movimiento a fin de obtener algoritmos
de simulación eficientes.
La evaluación de interacciones
de largo alcance entre partículas es extremadamente cara si usamos métodos convencionales (es
decir las sumas de Ewald), su costo computacional escala típicamente como N2
(siendo N el número de partículas), rápidamente excede cualquier límite
razonable. Como se comentó anteriormente, las interacciones de largo alcance es uno de los mayores problemas con los que
hay que lidiar en las simulaciones atomísticas de Dinámica Molecular. El problema es particularmente
agudo en el caso de biopolímeros donde la presencia de carga neta distribuida
hace que el potencial local oscile. La naturaleza condicionalmente convergente de
las series de energía electrostática para un sistema periódico tal como la caja de Dinámica Molecular en
condiciones de contorno periódicas (Periodic Boundary Conditions) hace
de cualquier método posterior de corte es esencialmente incosteable.
El campo de reacción es en
principio un método físico que corrige las contabilidades de los efectos de
largo alcance y requiere solamente un limitado esfuerzo computacional. La
técnica asume interacciones electrostáticas explícitas con una esfera limitada
y rodeada por un medio dieléctrico el cual ejerce sobre la esfera una
polarización o un “campo de reacción”.
El medio dieléctrico tiene una constante dieléctrica, que es el mismo que el de
la esfera interna. El método se ha probado y ha dado resultados idénticos a los
obtenidos con el método exacto de
Ewald en una simulación de Monte Carlo
de esferocilindros bipolares donde la constante dieléctrica que introdujeron en
el campo de reacción es cambiada periódicamente de acuerdo al valor encontrado
en la esfera. El método del campo de reacción de cualquier manera sufre de dos
grandes desventajas que limita fuertemente su uso.
Hemos descrito como obtener
interacciones de largo alcance, ahora en la separación de la escala de tiempo
de los potenciales modelo de sistemas moleculares complejos. Adicionalmente,
proveemos un subdivisión del potencial general aplicándolo a sistemas
biológicos, al igual que en otros sistemas químicos interesantes incluyendo
cristales líquidos. Este tipo de sistema son típicamente caracterizados por una
gran flexibilidad y por fuertes interacciones Coulombianas intermoleculares.
Esquemáticamente, podemos escribir el potencial V debido a dos contribuciones:
V = Vbnd
+ Vnbn
Aquí, la parte ligada o
intramolecular Vbnd es
rápida y es responsable de la flexibilidad del sistema. La parte no ligada o
intermolecular (intergrupal) Vnbn
es dominada por las interacciones coulombianas. El animo de la sección
siguiente es describir un protocolo general para las subdivisiones de tales
formas de potencial de interacción y para mostrar como obtener razonablemente
eficientes y transferibles integradores de múltiples pasos temporales válidos
para cualquier sistema molecular complejo.
Por el momento nos vamos a
concentrar en el potencial Vnbn
, que se encarga de determinar y modelar las interacciones coulombianas, este
potencial se determina de manera exacta con las sumas de Ewald.
El
método de Ewald de partículas en redes suaves
Antes de discutir la separación de múltiples pasos temporales no ligadas, es útil describir con detalle una
de las técnicas más avanzadas para tratar las fuerzas de largo alcance. En
lugar de esto, este tipo de fuerzas no ligadas son las más engorrosas para
lidiar y merecen un escrutinio detallado.
En la literatura reciente, pocas
técnicas están disponibles para tratar este problema de las interacciones de largo
alcance en las simulaciones
computacionales de partículas cargadas a diferente nivel de aproximación [1, 2,
3]. El
método de Ewald provee el resultado exacto para la energía electrostática de un
sistema periódico consistente de cajas infinitas replicadas neutras de partículas cargadas. El
método es la elección natural en simulaciones de Dinámica Molecular de sistemas
moleculares complejos.
El potencial de Ewald [2] está
dado por :
(1)
(2)
con
(3)
(4)
donde, ri es el vector de posición de la carga atómica qi, rij = ri-rj, rn es el vector de la caja,
es la función de error
complementario, erf(x) = 1 – erfc(x), V el volumen de la celda unitaria, m el vector de la caja reciproca y
a es el parámetro de
convergencia de Ewald. En la parte de la distancia directa a la caja, la ecuación (1), se omiten los contactos intramoleculares. En adición, en la Ec. (2)
el término Vintra sustrae,
en el espacio directo, la energía intra-molecular entre pares ligados, la cual
es incluida automáticamente en la parte derecha
de esta ecuación. En
consecuencia, la sumatoria sobre i y j en la Ec. (4) es sobre todos los
contactos intramoleculares excluidos. Debemos remarcar que en el potencial de
Ewald [7] hemos asumido las condiciones de frontera también llamadas
condiciones de contorno laminar:
1)La esfera de Ewald [7] es
inmersa en un medio perfectamente conductor y entonces el termino bipolar en la
superficie de la esfera de Ewald es cero. Para sistemas grandes el costo computacional
de la sumatoria estándar de Ewald, la cual se incrementa como N2,
se vuelve muy larga para aplicaciones practicas. Algoritmos alternativos los
cuales se incrementan en una potencia menor que N que el estándar de Ewald se
han propuesto en el pasado. Entre los algoritmos más rápidos diseñados para
sistemas periódicos es el algoritmo de mallas de partículas de Ewald (PME,
Particle Mesh Ewald) [4,5], inspirados por el método de malla de partículas
desarrollado por Hockney y Eastwood [6]. Aquí una aproximación multidimensional
de piezas discretas interpoladas son
usadas para computar la energía de la caja reciproca, Vqr, de la
ecuación (2), donde la parte directa, Vqd, es computada
directamente. El bajo costo computacional de el método de PME permite la
elección de grandes valores del parámetro de convergencia de Ewald a, comparado con el usado en el
método convencional de Ewald. Correspondientemente, menores distancias de
corte en el espacio directo de Ewald son adoptados. Si uj es la coordenada de escala fraccional de la i-ésima
partícula, el factor de estructura S(m)
en la Ec. (3), puede ser reescrito como :
(5)
Donde, N es el número de partículas, K1,
K2, K3 y m1,
m2, m3 son enteros. La componente a de la coordenada de
escala fraccional del i-ésimo átomo
puede escribirse como:
(6)
donde , = 1,2,3 son los vectores base
de las cajas recíprocas.
S(m) en la Ec. (5)
puede mirarse como una transformada discreta de Fourier (FT) de un conjunto de
cargas puestas irregularmente dentro de
la celda unitaria. Se han desarrollado técnicas en el pasado para aproximar S(m) con expresiones que involucren la
transformada de Fourier de una malla regular de puntos. Tales aproximaciones de
factor de estructura son ventajosas desde el punto de vista computacional
debido a que estas pueden evaluarse de las transformadas de Fourier rápidas
(FFT). Todas estas FFT, son en cierto sentido aproximaciones para juntar las
cargas sobre una malla cercana de puntos para producir una distribución regular de carga. El método de PME soluciona este problema por
interpolación. Entonces, las exponenciales complejas
Computadas en la posición de
la i-ésima carga en la Ec. (5), son rescritas como una suma de coeficientes
interpolados multiplicados por sus valores en los puntos cercanos a la reja. En
la versión suave de PME (SMPE) [5], la suma es multiplicada por un factor
apropiado, llamado:
(7)
donde n es el orden de la interpolación , define los
coeficientes de la coordenada de interpolación en la coordenada de escala . En la Ec. (7) la suma sobre k, representan los puntos de la reja, es
solo sobre un rango finito de enteros, debido a que las funciones son cero fuera del intervalo 0 £ u £ n. Debe hacerse
hincapié en que los coeficientes complejos b(mi)
son independientes de las coordenadas de las cargas ui y se necesitan computar sólo al inicio de la
simulación. Una derivación detallada de las funciones y de los coeficientes ba son dados en la referencia [33].
Insertando la ecuación (7) en (5), S(m)
puede rescribirse como:
(8)
donde es la transformada de Fourier
discreta (FT) en el punto de la reja del arreglo con 1£ ki £ Ki, i = 1,2,3. El arreglo de carga enrejada , es definida como:
(9)
Insertando el
factor de estructura de la Ec. (2) y usando el hecho de que , la energía de la caja
reciproca SPME puede escribirse como:
= (10)
con:
(11)
(12)
(13)
Usando el teorema de convolución para FFT
la energía (10) puede rescribirse como:
(14)
Ahora usamos
la identidad
para llegar a
(15)
Primero
notamos que no depende de las posiciones de las cargas y
que es diferenciable para n >2
(la cual es siempre el caso en aplicaciones practicas). Entonces la fuerza
sobre cada carga puede obtenerse tomando la derivada de la Ec. (15), llamado
=
(16)
En la práctica,
el cálculo es llevado a cabo de acuerdo
al siguiente esquema:
i)
A cada paso de la simulación se calcula
las coordenadas fraccionales de la reja y llena un arreglo con Q de acuerdo a la Ec. (9). En este
momento, las derivadas de las funciones Mn
son también calculadas y puestas en la memoria.
ii)
El arreglo que contiene a Q es entonces sobrescrito con , i.e. la transformada tridimensional de
Fourier de las Q’s.
iii)
Subsecuentemente, la energía electrostática es calculada vía la ecuación
(10). Al mismo tiempo, el arreglo conteniendo es sobrescrito por el producto
de el mismo con el arreglo conteniendo BC (calculado al inicio de la corrida).
iv)
El arreglo resultante es entonces la
transformada de Fourier para obtener la convolución .
v)
Finalmente, las fuerzas son calculadas vía la ecuación (16) usando las
derivadas calculadas con anterioridad de las funciones Mn, para formar .
La memoria requerida para el
método de SPME es limitada. Las variables de doble precisión son necesarios para el arreglo de la
reja de carga Q, cuando los cálculos de las funciones y sus derivadas requieren solamente 6 X n X N números de doble precisión. Los enteros determinan la finura de la
reja a lo largo del vector a-ésimo de la celda
unitaria. La precisión de salida de la
energía y fuerzas dependen de los parámetros del SPME:
-El parámetro de convergencia a, el espaciado de
la reja y el orden
n de la interpolación. Para un
valor típico a »
0.4A-1
una precisión
relativa entre 10-4 y 10-5 para
la energía electrostática son
obtenidos cuando el espaciado de la reja es alrededor de 1A a lo
largo de cada eje, y el orden n de la interpolación B es 4 o 5.
La potencia
del algoritmo SPME, comparado con la implementación estándar del método de
Ewald, es muy grande. En la Fig 2. se reporta el tiempo de CPU obtenido en una
computadora Pentium 3 para la evaluación de la energía de la caja y las fuerzas
obtenidas vía SPME como una función del número de átomos en el sistema. Se
usaron rutinas de dominio público de la rutina de 3-D FFT fueron usadas. El
algoritmo es prácticamente lineal y para 12 000 partículas el SPME toma solamente 2
segundos del CPU para llevar a cabo los cálculos. Una simulación estándar de
Ewald para una caja de 64 X 64 X 64 A3 (i.e. con una espaciado de
reja en el espacio k de » 0.01 A-1)
para el mismo ejemplo y al mismo nivel de precisión podría tomar varios
minutos.
REFERENCIAS
[1] P. Ewald. Ann.
Phys., 64:253, 1921.
[2]
S.W. deLeeuw, J.W. Perram, and E. R Smith. Proc.
R. Soc. London A, 373:27, 1980
[3]
J.P. Hansen. Molecular-dynamics simulation of coulomb systems in two and three
dimensions. In Molecular Dynamics
Simulation of Statistical-Mechanics Systems, Proceedings of the
International School of Physics “Enrico Fermi”. North Holland Physics, 1986.
[4]
T. Darden, D. York, and L. Pedersen. J.
Chem. Phys., 98:10089, 1993.
[5]
U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L.G. Pedersen. J. Chem. Phys., 101:8577, 1995.
[6]
R. W. Hockney. Computer Simulation Using
Particles. McGraw-Hill, New York, 1989.
[7]
M.P. Allan, D.J. Tildesley, COMPUTER SIMULATIONS OF liquids, Oxford Science
Publications.
[8]
Massimo Marchiand Piero Procacci , Orac
Manual and Guide,
A Molecular Dynamics Program to Simulate
Complex Molecular Systems with Realistic Interactions.
FORTRAN
PRINT i1 , PARA ENTEROS INDICANDO EL NUMERO DE DIGITOS
No hay comentarios:
Publicar un comentario