martes, 20 de agosto de 2019

Comparison of results of the Ewald Particle Grid method with the standard Ewald method.


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) = 1erfc(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 A4, PARA CARACTERES 4 ES EL NUMERO DE CARACTERES
PRINT i1  , PARA ENTEROS INDICANDO EL NUMERO DE DIGITOS





















No hay comentarios:

Publicar un comentario

zen consultora

Blogger Widgets

Entrada destacada

Platzy y el payaso Freddy Vega, PLATZI APESTA, PLATZI NO SIRVE, PLATZI ES UNA ESTAFA

  Platzy y los payasos fredy vega y cvander parte 1,  PLATZI ES UNA ESTAFA Hola amigos, este post va a ir creciendo conforme vaya escribiend...