PDF de programación - Dinámicas Híbridas QM/MM aceleradas en GPUs

Dinámicas Híbridas QM/MM aceleradas en GPUs

Mariano C. González Lebrero, Matías A. Nitsche y Darío A.

Estrin

IQUIFIB-CONICET, INQUIMAE-CONICET, FCEN-UBA.

1 de agosto de 2012

Denición: Simular es experimentar con un modelo.

Simulación

TeoríamodeloExperimentoSimulación Modelos.

Simulación

MM

Simulación

Ya lo vieron ayer...

Sin detalle elctrónico

Potenciales parametrizados.

Aplicable a cientos de miles de átomos y en escalas de los
nanosegundos o más.

QM

Simulación

ˆHψ = Eψ

Parece simple, pero...

No es razonable, intuitiva ni fácilmente comprensible.

No se puede resolver exactamente para más de 1 electrón.

Soluciones autoconsistentes mediante iteraciones.

Soluciones aproximadas, aunque pueden ser muy buenas.

Las soluciones mejores sólo se pueden aplicar a unos pocos
átomos en el vacío, para todo lo demas existe...

Métodos QM

Simulación

Hartree Fock

Semiempíricos.

HF

Post Hartree Fock (MP2,
CC, CI, etc.)

Híbridos (mezcla de DFT y HF)

B3LYP, PBE0, etc.

DFT

Semiempíricos

Local Density Aproximation
(LDA)

Generalized Gradient
Aproximation (GGAs: PBE,
BP86, etc.)

Simulación

QM-MM
Métodos Híbridos y Dinámica Molecular

Para sistemas complejos:
Métodos híbridos o QM-MM: porción clásica (soluto) +
porción cuántica (solvente)

QM/MM y dinámica molecular

Simulación

E = EQM + EM M + EQM/M M



N cargas

i=1

EQM/M M =

qi

Clas

cuan

i=1

j=1

ρ(r)

| Ri − r | dr+

[VLJ (R)+

qiZj

| Ri − rj | ]

Para hacer dinámica necesitos las fuerzas, o sea las derivadas de la
energía.

F = −∇E

Haciendo un DFT eciente...

Simulación

Partimos de un programa existente, escrito en el año 1992 por
Darío Estrin, en fortran 77.

1 Nuevos algoritmos, basados en el trabajo de Stratmann et al.

Menor complejidad
Variaciones novedosas, teniendo en cuenta: Procesador grácos
y los tipos de sistemas moleculares que se busca tratar

2 Nueva implementación

Extiende el programa existente (Fortran77)

Se modica la porción más signicativa del cálculo (XC) para
hacerlo en la GPU.

Código para procesador gráco

Particularidades del hardware y evaluación de alternativas de
implementación

Optimizaciónes a mansalva (mayormente reordenamientos de
loops) y uso de memoria dinámica.
Empleo de bibliotecas MKL de intel (parelelas).
otros.

La cuenta

Método de Kohn-Sham
La energía

E[ρ] = Ts[ρ]+Vne[ρ]+

ρ(r1)ρ(r2)

r12

1
2

d r1dr2+Exc[ρ]

Exc : energía de intercambio y correlación

Noc
M

i=1

µ=1

ρ(r) =

χi(r) =

|χi(r)|2

Ciµφk(r)

... se llega a que

→ ρ(r) =

M

k=1

Noc

i=1

2

Cikφk(r)

Las φk son funciones llamadas
base, en nuestro caso son
Gaussianas.

Se trata de encontrar los
coecientes Cik que minimicen la
energía.

Esquema de cálculo

La cuenta

InicializaciónFunción aproximadaDensidad electrónicaMatriz de FockHasta convergenciaEnergía y fuerzasmuevo los núcleosInicializaciónFunción aproximadaDensidad electrónicaMatriz de FockHasta convergenciaEnergía y fuerzas Energía de Intercambio y Correlación

La cuenta

Múltiples métodos de aproximación, el más simple es LDA:



Exc[ρ] =

ρ(r)εxc(ρ(r))dr

εxc : funcional de intercambio y correlación
Si sólo depende de la densidad.
Otros método más preciso: GGA (Global Gradient Approach).



Exc[ρ] =

ρ(r)εxc(ρ(r),∇ρ(r))dr

No existe solución analítica de esta integral, se integra
numéricamente sobre una grilla.


I

I =

i

F (r)dr

ωiF (ri)

La Grilla

La cuenta

Se toma una grilla esférica cómo base.
Superposición de grillas base, c/escalamiento: capas
concéntricas, sobre cada átomo
En sistemas poliatómicos: unión de las grillas de c/átomo

Las funciones.

La cuenta

La función de onda se describe como combinación lineal de un
grupo de funciones, llamado Base.
La base puede ser de funciones Gaussianas, ondas planas, etc.
En nuestro caso son Gaussianas, es decir tienen esta pinta:

φ = (x − x0)n(y − y0)m(z − z0)leα|−→r −−→r0|2

Donde n, m, y l son número enteros que dan el momento angular a
la función.
Usar funciones Gaussianas es venajoso porque reproducen
adecuadamente la distribución electrónica y son matemáticamente
amigables.

La densidad.

La cuenta

Los electrones estan en orbitales moleculares:

χi(r) =

Ciµφk(r)

µ=1

M
M
Noc
M
M

i=1

k=1

Entonces la densidad electrónica total es la suma de las
contribuciones de cada orbital:
→ ρ(r) =

Cikφk(r)

O dando vuelta la sumatoria se puede llegar a que:

2

→ ρ(r) =

Pjkφk(r)φj(r)

Donde Pijes la llamada matriz densidad.

j=1

k=1

Implementación

Empezando a pensar algoritmos ecientes.

Costo de calcular la integral de intercambio y corralacion=
αM 2 × N donde M es el número de funciones base y N el número
de puntos de la grilla.

¾Que hacemos?

Respuesta: funciones signicativas.

Implementación

funciones signicativas.

1 Agrupar los puntos de la grilla.

2 Encontrar las funciones signicativas

φ es signicativa en un grupo si φ(ri)
para algún punto (o
lugar) del grupo.

Ahora hay que denir el tamaño de los grupos y el valor de
.
Grupos muy chicos -> menos funciones x grupo pero más
grupos.

grande -> menos funciones signicativas pero mayor error.

Si hacemos esto para sistemas grandes el número de funciones
signicativas no crece con el tamaño del mismo.

½Escalamiento lineal!

Implementación

Algoritmos ecientes.
Ej: Particionado del sistema y funciones signicativas.

Algoritmo usual basado en cubos ⇒Muchos cubos con pocos
puntos (malo para GPU).
Algoritmo basado en esferas y cubos ⇒Mucho más eciente.

Implementación

El nuevo elefante: Coulomb

ρ(r1)ρ(r2)
CijCikΦiΦj =

Donde

d r1dr2

r12



PijΦiΦj

i(base)

j(base)

ρ(r) =

i(ocup).



j(Base)

ECoul = 1
2


ρ ∼= ˜ρ =

k(base)

ECoul ≈ ρ(r1)˜ρ(r2)

r12

Auxiliar

d r1dr2 − 1

2

CkGk

˜ρ(r1)˜ρ(r2)

d r1dr2

r12

Implementación

Coulomb

Se pueden tirar términos basado en el teorema de prouctos de
Gaussianas.

e−a(−→r −−→ra)2 × e−b(−→r −−→rb )2 = e

a+b (−→ra−−→rb )2
− ab

−(a+b)(−→r − a−→ra+b−→rb
e

a+b

)

Escala cuadráticamente.

Se puede guardar en memoria y no recalcular en cada iteración
(los productos de las Gaussianas no cambian, sólo los
coecientes Pij)
El cálculo de las derivadas (para la fuerza) es lo más
demandante para sistemas medianos (HEMO), junto con la
integral QM/MM.

Implementación

Oootro tema: Precisión.

Energía de formación de un agregado.

24 H2O −−→ (H2O)24

∇E = E((H2O)24) − 24E(H2O)

∇E = −1832,896H − 24 × (−76,352H) = 0,437H

La presición química es del orden de 1 kcal/mol o sea 10−3Hartrees!
Es indispensable usar presición alta siempre, no??

Implementación

Respuesta: no

Se puede usar precisión simple para muchas partes
manteniendo la calidad del resultado.

Para otras se puede usar precisión mixta.

Mejora los tiempos de cálculo (½sobre todo en GPUs!).

Disminuye el uso de memoria (muy útil en Coulomb).

Implementación

uso de presición simple.

En la integral de intercambio y correlación para todo salvo para la
acumulación.
Para las integrales de Coulomb se aplica el mismo criterio basado
enel producto de Gaussianas. Si son muy chicas, se tiran, si tienen
un valos intermedio se hacen en simple, si tienen un valor mayor en
doble.

Implementación

Calidad.

510152025Exponente-0,6-0,5-0,4-0,3-0,2-0,100,1KCal/molGPU-simpleGPU-dobleCPU-simpleCPU-dobleDif. energia de formaciondel agua 24 Implementación

Calidad

Usando precision simple se comete un error en la energía de
0,2 Kcal/mol para los sitemas mas grandes (incluso para la
valinomicina, de 168 átomos).

Con exponentes mayores a 10 ya se obtiene un buen resultado.

GPU en doble es equivalente a CPU en doble.

Implementación

Parametros nuevos.

Tienen impacto en la performance y la calidad.

max_function_exponent => Dene el valor de las funciones
que se tiran.

Tamaño de los cubos.

Tamaxo de las esferas (entre 0. y 1. )

rmax => Dene los pares de Gaussianas que se tiran cuando
hace Coulomb.

rmaxs => Dene los pares de Gausianas que se hacen en
doble.

Escalamiento y aceleración.

Resultados

Taxol-631G

Valinomicina-631G

Hemo-DZVP

Generacion de la grilla

Iteración SCF

Intercambio y correlación

GPU
0.99
3.70
2.41

CPU/GPU

19.53

8.19

11.96

GPU
1.63
6.40
3.5

CPU/GPU

18.88

6.03

10.09

GPU
0.42
2.92
2.3

CPU/GPU

18.79

9.13

11.26

Scalability anlysisexchange-correlation time [s]00.20.40.60.811.21.4# water molecules5101520Low-densityMedium-densityHigh-densityGPU/CPU execution time relationx0246810121416# water molecules5101520Low-densityMedium-densityHigh-density Garchamber

Resultados

Ahora lo bueno, a este código mejorado lo pegué al amber.

Usé buena parte del código de QM/MM de Amber (para
DFTB y PM3).

Neceista un achivo de input extra con las bases y los
parámetros del DFT. (y se tiene que llamar input)

Pase todo a memoria dinámica (digamos un 90 %) y las
variables en un módulo (no mas COMM).

Compilado como biblioteca.

Cut o, no Ewald.

¾Cómo anda?

Resultados

Sistema

segundos por paso de dinámica

agua1 en agua
agua2 en agua
agua3 en agua

Hemo O2

trip en agua

0.56
1.3
2.3
180
20

Resultados para el agua

Resultados

40000 pasos de dinámica en 6 horas y media

24681234g(r)010002000300040000,0050,010,0150,020,025 Conclusiones.

¾Qué se logró?

1 Una implementación de DFT con variaciones novedosas:

Esferas → mejora desempeño en GPU
En sistemas r