Mathematics meets chemistry: a new paradigm for implicit solvation - - PowerPoint PPT Presentation
Mathematics meets chemistry: a new paradigm for implicit solvation models Journes inaugurales de la machine de calcul, 11/6/2015 P6 Benjamin Stamm Laboratoire Jacques-Louis Lions Universit Pierre et Marie Curie and CNRS Mathematical
Search phrase: Google Scholar MathSciNet ratio Navier Stokes 537 000 8 220 65.3 Maxwell equations 194 000 1 892 102.5 Molecular dynamics 1 550 000 325 4 769.2 Density functional theory 1 200 000 124 9 677.4 Particle Mesh Ewald 24 200 1 24 200 Continuum solvation model 6 320 Polarizable force field 3 670
Search phrase: Google Scholar MathSciNet ratio Navier Stokes 537 000 8 220 65.3 Maxwell equations 194 000 1 892 102.5 Molecular dynamics 1 550 000 325 4 769.2 Density functional theory 1 200 000 124 9 677.4 Particle Mesh Ewald 24 200 1 24 200 Continuum solvation model 6 320 Polarizable force field 3 670
Mechanics & Engineering 13.3% Physics 14.0% Earth & Environmental Science 6.5% Life Science 6.1% Others 0.9% Chemistry & Materials 59.2%
Search phrase: Google Scholar MathSciNet ratio Navier Stokes 537 000 8 220 65.3 Maxwell equations 194 000 1 892 102.5 Molecular dynamics 1 550 000 325 4 769.2 Density functional theory 1 200 000 124 9 677.4 Particle Mesh Ewald 24 200 1 24 200 Continuum solvation model 6 320 Polarizable force field 3 670
Mechanics & Engineering 13.3% Physics 14.0% Earth & Environmental Science 6.5% Life Science 6.1% Others 0.9% Chemistry & Materials 59.2%
Search phrase: Google Scholar MathSciNet ratio Navier Stokes 537 000 8 220 65.3 Maxwell equations 194 000 1 892 102.5 Molecular dynamics 1 550 000 325 4 769.2 Density functional theory 1 200 000 124 9 677.4 Particle Mesh Ewald 24 200 1 24 200 Continuum solvation model 6 320 Polarizable force field 3 670
Mechanics & Engineering 13.3% Physics 14.0% Earth & Environmental Science 6.5% Life Science 6.1% Others 0.9% Chemistry & Materials 59.2%
ematics can intervene and have an impact.
Before we start serious things… a little bit of fun: First article of three “naive” mathematicians: Domain decomposition for implicit solvation models E Cancès, Y Maday, B Stamm The Journal of Chemical Physics, 2013
Before we start serious things… a little bit of fun: First article of three “naive” mathematicians: Domain decomposition for implicit solvation models E Cancès, Y Maday, B Stamm The Journal of Chemical Physics, 2013 Review: “There's a useful numerical method here, at the heart of it all, but the presentation is a convoluted mess …, if it is to be useful to (and understood by) the chemical physics community.” “..., however the manuscript suffers greatly from the reliance on the authors trying to write and act like mathematicians” “It's just a 3-d differential equation that they are mapping onto a boundary-element method and solving via domain decomposition. It's not rocket science, it's just linear equations -- so why present the material in a complicated (and tediously long) way?”
Goal: model liquid!
Goal: model liquid! Schr¨
HΨ = EΨ.
state, exicted states, etc.
number of electrons (Born-Oppenheimer approxima- tion).
compute an approximate solution due to the curse
Goal: model liquid! Schr¨
HΨ = EΨ.
state, exicted states, etc.
number of electrons (Born-Oppenheimer approxima- tion).
compute an approximate solution due to the curse
H = −
N
X
j=1
1 2∆j +
M
X
α=1
Zα |Rα − xj| ! +
N
X
i,j=1 i6=j
1 |xi − xj|,
Goal: model liquid! Schr¨
HΨ = EΨ.
state, exicted states, etc.
number of electrons (Born-Oppenheimer approxima- tion).
compute an approximate solution due to the curse
Small molecule to be studied Solvent Focussed model: Iden- tify chemically interest- ing portion of the sys- tem.
H = E is still too complex.
approximation or the Hartree-Fock (HF) approximation. Based on mini- mizing a simplified energy functional: min EDF T
min EHF Thus, approximates only the ground state.
Second approximation: Do not neglect environment (solvent). Remarks:
demanding.
reasons: jumping from one local minimum to another is very easy and there- fore sampling is required. Thus, several DFT or HF computation are needed to describe the system.
the interaction with the solute, is a feasible strategy: solvation models as part of a multi-scale strategy.
Second approximation: Do not neglect environment (solvent). Remarks:
demanding.
reasons: jumping from one local minimum to another is very easy and there- fore sampling is required. Thus, several DFT or HF computation are needed to describe the system.
the interaction with the solute, is a feasible strategy: solvation models as part of a multi-scale strategy. Example: QM/MM-approach (Nobelprize 2013)
biology take place in the liquid phase.
biology take place in the liquid phase.
biology take place in the liquid phase.
Given the charge distribution ρ (classical point charges, electric dipoles and multi- poles in force-field models, classical nuclear charges and quantum electronic charge density in first-principle or semi-empirical models):
the solvent.
and denoted here by Es, must be added to the electrostatic energy computed in vacuo.
biology take place in the liquid phase.
Given the charge distribution ρ (classical point charges, electric dipoles and multi- poles in force-field models, classical nuclear charges and quantum electronic charge density in first-principle or semi-empirical models):
the solvent.
and denoted here by Es, must be added to the electrostatic energy computed in vacuo. Goal of this work: Compute the electro- static contribution (long-range interaction) to the solvation energy in an efficient way.
s Chatracteristics:
ular cavity Ω and is supported in Ω.
by an infinite, continuum polarizable dielectric situated in R3\Ω, described by its macroscopic dielectric suscepti- bility s.
s Chatracteristics:
ular cavity Ω and is supported in Ω.
by an infinite, continuum polarizable dielectric situated in R3\Ω, described by its macroscopic dielectric suscepti- bility s.
Energy contribution: Es = 1 2
where V solves div(V ) = 4, in R3\Ω, [V ] = 0,
∂n
and Φ(r) =
(r) |r r| dr, =
in Ω s in R3\¯ Ω
s dielectric situated in R3\Ω, described by its macroscopic dielectric suscepti- bility s.
Energy contribution: Es = 1 2
where V solves div(V ) = 4, in R3\Ω, [V ] = 0,
∂n
and Φ(r) =
(r) |r r| dr, =
in Ω s in R3\¯ Ω Let W = V Φ denote the reaction-field potential generated by the charge distri- bution in presence of the dielectric continuum. Then, the energy contribution writes Es = 1 2
and W solves: div(W) = 0, in R3\Ω, [W] = 0,
n
n ,
s dielectric situated in R3\Ω, described by its macroscopic dielectric suscepti- bility s.
Energy contribution: Es = 1 2
where V solves div(V ) = 4, in R3\Ω, [V ] = 0,
∂n
and Φ(r) =
(r) |r r| dr, =
in Ω s in R3\¯ Ω Let W = V Φ denote the reaction-field potential generated by the charge distri- bution in presence of the dielectric continuum. Then, the energy contribution writes Es = 1 2
and W solves: div(W) = 0, in R3\Ω, [W] = 0,
n
n ,
Rem: If s = ∞, then V = 0 on R3\¯ Ω.
s Chatracteristics: COSMO (COnductor-like Screening MOdel) approximation: the dielec- tric is replaced by a conductor, i.e. ✏s = ∞. Also called C-PCM.
Φ is the potential generated by ρ in vacuo: Φ(r) =
ρ(r) |r − r|dr. Energy contribution: Es
C = 1
2f(s)
⇥(r)W(r) dr, where f(s) =
s−1 s+0.5 is an empirical function of s.
W is solution to: Find W ∈ H1(Ω) such that
in Ω, W= −Φ
s Chatracteristics: COSMO (COnductor-like Screening MOdel) approximation: the dielec- tric is replaced by a conductor, i.e. ✏s = ∞. Also called C-PCM.
Φ is the potential generated by ρ in vacuo: Φ(r) =
ρ(r) |r − r|dr. Energy contribution: Es
C = 1
2f(s)
⇥(r)W(r) dr, where f(s) =
s−1 s+0.5 is an empirical function of s.
W is solution to: Find W ∈ H1(Ω) such that
in Ω, W= −Φ
s Chatracteristics: COSMO (COnductor-like Screening MOdel) approximation: the dielec- tric is replaced by a conductor, i.e. ✏s = ∞. Also called C-PCM. Ω The difficulty is not to solve Laplace’s equation, it is the extremely complex geometry in- cluding thousands of atoms. In addition, this energy needs to be computed a large num- ber of times!
Note: a cavity of a molecule is an artificial construct. Different cavities are possible:
Decompose the domain into overlapping subdomains: Find W ∈ H1(Ω) such that
in Ω, W= −Φ
Ω =
M
Ωi
Decompose the domain into overlapping subdomains: Find W ∈ H1(Ω) such that
in Ω, W= −Φ
Ω =
M
Ωi
Decompose the domain into overlapping subdomains: Formulate a local problem on each subdomain: Ω1 Ω2 Ω3 Ω4 Find W ∈ H1(Ω) such that
in Ω, W= −Φ
For each i = 1, . . . , M: find Wi ∈ H1(Ωi) such that
in Ωi, Wi= gi
Ω =
M
Ωi
Decompose the domain into overlapping subdomains: Formulate a local problem on each subdomain: Ω1 Ω2 Ω3 Ω4 Find W ∈ H1(Ω) such that
in Ω, W= −Φ
What to impose on the gi’s in order that W|Ωi = Wi? For each i = 1, . . . , M: find Wi ∈ H1(Ωi) such that
in Ωi, Wi= gi
Ω =
M
Ωi
Ω1 Ω2 Ω3 Ω4 Γe
2
Splitting of local boundary into exterior and interior boundary: Γi = Γi
i ∪ Γe i
Local boundary: Γi = ∂Ωi Γe
2
Ω1 Ω2 Ω3 Ω4 Γe
2
Splitting of local boundary into exterior and interior boundary: Γi = Γi
i ∪ Γe i
Local boundary: Γi = ∂Ωi Average of neighbour solutions. We impose that gi =
i,
Wi,N
i,
where Wi,N(s) = 1 |N(i, s)|
Wj(s) ∀s ∈ Γi
i.
Set of neighbouring/intersecting spheres at point s ∈ Γi
i.
Γe
2
Ω1 Ω2 Ω3 Ω4 Γe
2
Then, solve in an iterative manner... Splitting of local boundary into exterior and interior boundary: Γi = Γi
i ∪ Γe i
Local boundary: Γi = ∂Ωi Average of neighbour solutions. We impose that gi =
i,
Wi,N
i,
where Wi,N(s) = 1 |N(i, s)|
Wj(s) ∀s ∈ Γi
i.
Set of neighbouring/intersecting spheres at point s ∈ Γi
i.
Γe
2
For each i = 1, . . . , M, we need to find Wi ∈ H1(Ωi) such that
in Ωi, Wi= gi
For each i = 1, . . . , M, we need to find Wi ∈ H1(Ωi) such that
in Ωi, Wi= gi
The harmonic function Wi can be represented by the single-layer potential ∀s ∈ Ωi, Wi(s) = ˜ Siσi(s) :=
σi(s) |s − s| ds, where σi ∈ H 1
2 (Γi) solves the integral equation
Siσi = gi,
For each i = 1, . . . , M, we need to find Wi ∈ H1(Ωi) such that
in Ωi, Wi= gi
Here, Si : H 1
2 (Γi) → H 1 2 (Γi) denotes the trace of the single-layer potential onto
the surface Γi, defined by ∀s ∈ Γi, Siσ(s) =
σ(s) |s − s| ds. The harmonic function Wi can be represented by the single-layer potential ∀s ∈ Ωi, Wi(s) = ˜ Siσi(s) :=
σi(s) |s − s| ds, where σi ∈ H 1
2 (Γi) solves the integral equation
Siσi = gi,
Ωj s s s Ω s Ωi σk
i
σk
j
σk
Thus, in the spirit of the domain decomposition approach, we solve for each i = 1, . . . , M: find σk
i ∈ H− 1
2 (Γi) such that
Siσk
i = gk i ,
where gk
i (s) =
i, 1 |N(i,s)|
P
j∈N(i,s) ˜
Sjσk−1
j
(s)
i,
gk
i ∈ H
1 2 (Γi).
Ingredients of discretization:
Ingredients of discretization:
Spherical harmonics:
SS2Ylm =
Ylm(s) | · −s| ds = 4π 2l + 1Ylm
i by a truncated series of spherical harmonics
Σk
i (s) = N
l
(Σk
i )lmYlm
s − xi |s − xi|
Ingredients of discretization:
Spherical harmonics:
SS2Ylm =
Ylm(s) | · −s| ds = 4π 2l + 1Ylm
i by a truncated series of spherical harmonics
Σk
i (s) = N
l
(Σk
i )lmYlm
s − xi |s − xi|
At each iteration k and for each i = 1, . . . , M: find Σk
i 2 span{Ylm} such that
i , Ylm
i , Ylm
, SiΣk
i = ΠNGk i
where Gk
i (s) =
i, 1 |N(i,s)|
P
j∈N(i,s) ˜
SjΣk−1
j
(s)
i,
Gk
i 62 H
1 2 (Γi).
At each iteration k and for each i = 1, . . . , M: find Σk
i 2 span{Ylm} such that
i , Ylm
i , Ylm
, SiΣk
i = ΠNGk i
where Gk
i (s) =
i, 1 |N(i,s)|
P
j∈N(i,s) ˜
SjΣk−1
j
(s)
i,
Gk
i 62 H
1 2 (Γi).
Remark 1:
i in terms of spherical harmonics is know, i.e.
ΠNGk
i = N
X
l=0 m
X
m=−l
(Gk
i )lmYlm,
then (Σk
i )lm = (2l + 1)
4πRi (Gk
i )lm.
Remark 2:
i )lm at each iteration. By definition,
the coefficients of Gk
i are defined by
Z
Γi
Gk
i (s) Ylm(s) ds.
We replace the integral by the numerical integration to obtain (Gk
i )lm = Ng
X
n=1
ωn Gk
i (sn) Ylm(sn) ds.
At each iteration k and for each i = 1, . . . , M: find Σk
i 2 span{Ylm} such that
i , Ylm
i , Ylm
, SiΣk
i = ΠNGk i
where Gk
i (s) =
i, 1 |N(i,s)|
P
j∈N(i,s) ˜
SjΣk−1
j
(s)
i,
Gk
i 62 H
1 2 (Γi).
Ωi Ωj Ω sn sn sn
Σk
i
Σk−1
j
yn
At each iteration k and for each sphere i = 1, . . . , M, we need to access the neighbouring potential at the interior integration point sn ∈ Γi
i:
Gk
i (sn) =
1 |N(i, sn)| X
j∈N(i,sn)
˜ SjΣk−1
j
(sn). Compute yn = sn − xj |sn − xj| Σk−1
j
(yn) =
N
X
l=0 m
X
m=−l
(Σk−1
j
)lm Ylm (yn) , SjΣk−1
j
(yn) =
N
X
l=0 m
X
m=−l
2l + 1 4πRj (Σk−1
j
)lm Ylm (yn) , ˜ SjΣk−1
j
(sn) =
N
X
l=0 m
X
m=−l
✓|sn − xj| Rj ◆l 2l + 1 4πRj (Σk−1
j
)lm Ylm (yn) .
Ωi Ωj Ω sn sn sn
Σk
i
Σk−1
j
yn
At each iteration k and for each sphere i = 1, . . . , M, we need to access the neighbouring potential at the interior integration point sn ∈ Γi
i:
Gk
i (sn) =
1 |N(i, sn)| X
j∈N(i,sn)
˜ SjΣk−1
j
(sn). Compute yn = sn − xj |sn − xj| Σk−1
j
(yn) =
N
X
l=0 m
X
m=−l
(Σk−1
j
)lm Ylm (yn) , SjΣk−1
j
(yn) =
N
X
l=0 m
X
m=−l
2l + 1 4πRj (Σk−1
j
)lm Ylm (yn) , ˜ SjΣk−1
j
(sn) =
N
X
l=0 m
X
m=−l
✓|sn − xj| Rj ◆l 2l + 1 4πRj (Σk−1
j
)lm Ylm (yn) .
Ωi Ωj Ω sn sn sn
Σk
i
Σk−1
j
yn
At each iteration k and for each sphere i = 1, . . . , M, we need to access the neighbouring potential at the interior integration point sn ∈ Γi
i:
Gk
i (sn) =
1 |N(i, sn)| X
j∈N(i,sn)
˜ SjΣk−1
j
(sn). Compute yn = sn − xj |sn − xj| Σk−1
j
(yn) =
N
X
l=0 m
X
m=−l
(Σk−1
j
)lm Ylm (yn) , SjΣk−1
j
(yn) =
N
X
l=0 m
X
m=−l
2l + 1 4πRj (Σk−1
j
)lm Ylm (yn) , ˜ SjΣk−1
j
(sn) =
N
X
l=0 m
X
m=−l
✓|sn − xj| Rj ◆l 2l + 1 4πRj (Σk−1
j
)lm Ylm (yn) .
Ωi Ωj Ω sn sn sn
Σk
i
Σk−1
j
yn
At each iteration k and for each sphere i = 1, . . . , M, we need to access the neighbouring potential at the interior integration point sn ∈ Γi
i:
Gk
i (sn) =
1 |N(i, sn)| X
j∈N(i,sn)
˜ SjΣk−1
j
(sn). Compute yn = sn − xj |sn − xj| Σk−1
j
(yn) =
N
X
l=0 m
X
m=−l
(Σk−1
j
)lm Ylm (yn) , SjΣk−1
j
(yn) =
N
X
l=0 m
X
m=−l
2l + 1 4πRj (Σk−1
j
)lm Ylm (yn) , ˜ SjΣk−1
j
(sn) =
N
X
l=0 m
X
m=−l
✓|sn − xj| Rj ◆l 2l + 1 4πRj (Σk−1
j
)lm Ylm (yn) .
The presented scheme can be considered as a block-type Jacobi iteration scheme. Other iterative solvers (GMRes, Gauss-Seidel, . . . ) can also be considered to accelerate the convergence rate. The proposed iterative scheme is a way to solve the global block-sparse non- symmetric linear system: L11 · · · L1M . . . ... . . . LM1 · · · LMM Σ1 . . . ΣM = f1 . . . fM
The presented scheme can be considered as a block-type Jacobi iteration scheme. Other iterative solvers (GMRes, Gauss-Seidel, . . . ) can also be considered to accelerate the convergence rate. The proposed iterative scheme is a way to solve the global block-sparse non- symmetric linear system: L11 · · · L1M . . . ... . . . LM1 · · · LMM Σ1 . . . ΣM = f1 . . . fM Jacobi acceleration technique: Jacobi-like iterations are nevertheless well-suited for a parallel implementation: More, but more parallel work will result in a faster computation! In addition, convergence of iterative scheme can be accelerated by different techniques: pre- conditioning, DIIS (Direct Inversion in the Iterative Subspace), .... The DIIS post-processing works really good in practice.
The presented scheme can be considered as a block-type Jacobi iteration scheme. Other iterative solvers (GMRes, Gauss-Seidel, . . . ) can also be considered to accelerate the convergence rate. The proposed iterative scheme is a way to solve the global block-sparse non- symmetric linear system: L11 · · · L1M . . . ... . . . LM1 · · · LMM Σ1 . . . ΣM = f1 . . . fM In consequence:
to a BEM-implementation (dense solution matrix).
atoms if number of iterations is constant (see next slide).
Counter-example: N spheres It takes at least O(N) iterations until the Dirichlet BC’s arrive at the center!
Counter-example: N spheres It takes at least O(N) iterations until the Dirichlet BC’s arrive at the center! This is a slightly different point of view than in an usual DD-paradigm where the domain Ω is constant and the number of subdomains to cover Ω is increased.
Counter-example: N spheres It takes at least O(N) iterations until the Dirichlet BC’s arrive at the center! This is a slightly different point of view than in an usual DD-paradigm where the domain Ω is constant and the number of subdomains to cover Ω is increased. This luckily does not happen for chemically/biologically interesting molecules (proteins), in contrast to crystals (physics/material science).
Counter-example: N spheres It takes at least O(N) iterations until the Dirichlet BC’s arrive at the center! This is a slightly different point of view than in an usual DD-paradigm where the domain Ω is constant and the number of subdomains to cover Ω is increased. This luckily does not happen for chemically/biologically interesting molecules (proteins), in contrast to crystals (physics/material science).
The energy depends in a discontinuous way on the distance between two spheres! This is due to the numerical integration and may lead to non-physical local minima for the energy.
points.
Gk
i (sn) =
i, 1 |N(i,sn)|
SjΣk−1
j
(sn)
i,
is not continuous (only its exact counterpart), this implies a non-smooth dependency. The energy depends in a discontinuous way on the distance between two spheres! This is due to the numerical integration and may lead to non-physical local minima for the energy.
χ0
points.
Gk
i (sn) =
i, 1 |N(i,sn)|
SjΣk−1
j
(sn)
i,
is not continuous (only its exact counterpart), this implies a non-smooth dependency. We now write Gk
i (sn) = −(1 − χ0)Φ(sn) +
χ0 |N(i, sn)|
˜ SjΣk−1
j
(sn), The energy depends in a discontinuous way on the distance between two spheres! This is due to the numerical integration and may lead to non-physical local minima for the energy.
χη η
points.
Gk
i (sn) =
i, 1 |N(i,sn)|
SjΣk−1
j
(sn)
i,
is not continuous (only its exact counterpart), this implies a non-smooth dependency. We now write Gk
i (sn) = −(1 − χ0)Φ(sn) +
χ0 |N(i, sn)|
˜ SjΣk−1
j
(sn), and replace it by Gk
i (sn) = −(1−χη)Φ(sn)+
χη |N(i, sn)|
˜ SjΣk−1
j
(sn), The energy depends in a discontinuous way on the distance between two spheres! This is due to the numerical integration and may lead to non-physical local minima for the energy.
tion of each atom) is an important quantity for geometry optimization or time-dependent simulations.
tion of each atom) is an important quantity for geometry optimization or time-dependent simulations. Let xq denote the coordinate of the q-th nuclei. Then we would like to compute Fq = rqEs
C
where rq denotes the derivative with respect to the position xq.
tion of each atom) is an important quantity for geometry optimization or time-dependent simulations. The energy can be written as Es
C = M
X
i=1
hΨi, Σii for some given multi-polar expansion Ψ and where ha, bi =
N
X
l=0 l
X
m=−l
(a)lm(b)lm. Let xq denote the coordinate of the q-th nuclei. Then we would like to compute Fq = rqEs
C
where rq denotes the derivative with respect to the position xq.
Consider the global system by LΣ = f ,
M
X
i=1
LjiΣi = fj, j = 1, . . . , M. Take the derivative to obtain: rqfj = rq
M
X
i=1
LjiΣi =
M
X
i=1
rq(Lji)Σi +
M
X
i=1
Ljirq(Σi).
Consider the global system by LΣ = f ,
M
X
i=1
LjiΣi = fj, j = 1, . . . , M. Take the derivative to obtain: rqfj = rq
M
X
i=1
LjiΣi =
M
X
i=1
rq(Lji)Σi +
M
X
i=1
Ljirq(Σi). Thus
M
X
i=1
Ljirq(Σi) = rqfj
M
X
i=1
rq(Lji)Σi | {z }
=(hq)j,
(closed-form expression possible)
, and globally L(rqΣ) = hq , rqΣ = L−1hq , rqΣj =
M
X
i=1
L−1
ji (hq)i,
j = 1, . . . , M.
rqfj = rq
M
X
i=1
LjiΣi =
M
X
i=1
rq(Lji)Σi +
M
X
i=1
Ljirq(Σi). Thus
M
X
i=1
Ljirq(Σi) = rqfj
M
X
i=1
rq(Lji)Σi | {z }
=(hq)j,
(closed-form expression possible)
, and globally L(rqΣ) = hq , rqΣ = L−1hq , rqΣj =
M
X
i=1
L−1
ji (hq)i,
j = 1, . . . , M. Then, Fq = rq
M
X
j=1
hΨj, Σji =
M
X
j=1
hΨj, rqΣji =
M
X
j=1
* Ψj,
M
X
i=1
L−1
ji (hq)i
+ =
M
X
i=1
* M X
j=1
L−T
ji Ψj, (hq)i
+ =
M
X
i=1
hsi, (hq)ii where s solves the adjoint problem LT s = Ψ.
10 20 30 40 50
N
1x10-6 0.00001 0.0001 0.001 0.01
relative error
d=1.25 d=1.5 d=1.75 10 20 30 40 50
N
0.01 0.1
relative error
d=1.25 d=1.5 d=1.75
σ H−1/2
10 20 30 40 50
N
8 10 12 14 16 18 20
number of iterations
d=1.25 d=1.5 d=1.75 1 10 100
N
0.0001 0.001 0.01
relative error
d=1.25 d=1.5 d=1.75 1 10 100
N
0.01 0.1
relative error
d=1.25 d=1.5 d=1.75
σ H−1/2
10 20 30 40 50
N
10 15 20 25 30 35
number of iterations
d=1.25 d=1.5 d=1.75
η = 0: η = 0.1:
CSC: Continuous Surface Charge by York and Karplus
2 4 6 8 10 12 14 16 18 20 22 2000 3000 4000 5000 6000 7000 8000 9000 50 100 150 200 250 300 350 Wall Time (s) Memory (Integer MegaWords) Number of atoms Hemoglobin and its subunits direct system adjoint system Forces (matrix part) Memory - cavity and scratch Memory - Geometrical parameters
10 20 30 40 50 60 70 80 90 5000 10000 15000 20000 25000 30000 35000 40000 200 400 600 800 1000 1200 1400 1600 Wall Time (s) Memory (Integer MegaWords) Number of atoms Alanine chains direct system adjoint system Forces (matrix part) Memory - cavity and scratch Memory - Geometrical parameters
20 40 60 80 100 120 140 160 10000 20000 30000 40000 50000 60000 70000 Number of iterations Number of spheres Jacobi iterations DIIS iterations
50 100 150 200 250 300 350 50 100 150 200 250 300 350 performance gain number of cores Ideal Speedup
Molecular dynamics simulation of ubiquitin in water (100ps), with the solvent described by a polarizable continuum, using a time step of 1fs
The (ddCOSMO) linear system needs to be inverted ~500’000 times. Whether you have 3 seconds or 25 minutes decides whether one makes such a simulation or not!
Polarizable Molecular Dynamics in a Polarizable Continuum Solvent
Molecular dynamics simulation of ubiquitin in water (100ps), with the solvent described by a polarizable continuum, using a time step of 1fs
The (ddCOSMO) linear system needs to be inverted ~500’000 times. Whether you have 3 seconds or 25 minutes decides whether one makes such a simulation or not!
Polarizable Molecular Dynamics in a Polarizable Continuum Solvent
Molecular dynamics simulation of ubiquitin in water (100ps), with the solvent described by a polarizable continuum, using a time step of 1fs
The (ddCOSMO) linear system needs to be inverted ~500’000 times. Whether you have 3 seconds or 25 minutes decides whether one makes such a simulation or not!
Table 1. Stability of NVE Simulations, as a Function of the Convergence Threshold for the Coupled Polarization Equations and
Convergence = 10−4 Convergence = 10−5 Convergence = 10−6 time step (fs) STF LTD ACT STF LTD ACT STF LTD ACT 1.00 0.42 20.7 0.35 0.46 18.5 0.46 0.43 20.7 0.75 0.50 0.11 8.7 0.35 0.10 2.1 0.46 0.10 4.4 0.75 0.25 0.03 3.6 0.35 0.02 0.0 0.46 0.02 0.5 0.75
aFor each combination, the short-time average fluctuation (STF), the long-time drift (LTD), and the average computational time per time step
(ACT) are reported. All the energies are given in terms of kcal/mol, all timings are given in seconds.
Table 2. Stability of NVE Simulations, as a Function of the Convergence Threshold for the Coupled Polarization Equations and
Convergence = 10−4 Convergence = 10−5 Convergence = 10−6 time step (fs) STF LTD ACT STF LTD ACT STF LTD ACT 1.00 0.46 −1.6 1.18 0.46 3.1 1.34 0.45 9.0 1.72 0.50 0.11 10.4 1.18 0.11 1.5 1.34 0.11 1.2 1.72 0.25 0.03 2.1 1.18 0.03 0.2 1.34 0.03 −0.1 1.72
aFor each combination, the short-time average fluctuation (STF), the long-time drift (LTD), and the average computational time per time step
(ACT) are reported. All the energies are given in terms of kcal/mol, all timings are given in seconds.
Polarizable Molecular Dynamics in a Polarizable Continuum Solvent
Dual Xeon Model E5-2650 cluster node (16 cores, 2 GHz) equipped with 64GB of DDR3 memory (on TACC)
Molecular dynamics simulation of ubiquitin in water (100ps), with the solvent described by a polarizable continuum, using a time step of 1fs
The (ddCOSMO) linear system needs to be inverted ~500’000 times. Whether you have 3 seconds or 25 minutes decides whether one makes such a simulation or not!
Table 1. Stability of NVE Simulations, as a Function of the Convergence Threshold for the Coupled Polarization Equations and
Convergence = 10−4 Convergence = 10−5 Convergence = 10−6 time step (fs) STF LTD ACT STF LTD ACT STF LTD ACT 1.00 0.42 20.7 0.35 0.46 18.5 0.46 0.43 20.7 0.75 0.50 0.11 8.7 0.35 0.10 2.1 0.46 0.10 4.4 0.75 0.25 0.03 3.6 0.35 0.02 0.0 0.46 0.02 0.5 0.75
aFor each combination, the short-time average fluctuation (STF), the long-time drift (LTD), and the average computational time per time step
(ACT) are reported. All the energies are given in terms of kcal/mol, all timings are given in seconds.
Table 2. Stability of NVE Simulations, as a Function of the Convergence Threshold for the Coupled Polarization Equations and
Convergence = 10−4 Convergence = 10−5 Convergence = 10−6 time step (fs) STF LTD ACT STF LTD ACT STF LTD ACT 1.00 0.46 −1.6 1.18 0.46 3.1 1.34 0.45 9.0 1.72 0.50 0.11 10.4 1.18 0.11 1.5 1.34 0.11 1.2 1.72 0.25 0.03 2.1 1.18 0.03 0.2 1.34 0.03 −0.1 1.72
aFor each combination, the short-time average fluctuation (STF), the long-time drift (LTD), and the average computational time per time step
(ACT) are reported. All the energies are given in terms of kcal/mol, all timings are given in seconds.
Polarizable Molecular Dynamics in a Polarizable Continuum Solvent
Dual Xeon Model E5-2650 cluster node (16 cores, 2 GHz) equipped with 64GB of DDR3 memory (on TACC)
Figure 1. Scaling of the MPI code for the coupled induced dipoles/ ddCOSMO problem, with respect to a single-node (16 cores) computation.
Figure 7. Total elapsed time (in seconds) for the solution of the coupled MMPol/ddCOSMO equations as a function of the system size.
At each Self-Consistent Field (SCF) iteration devoted to the QM-part, one needs to solve the coupled polarization equation.
Figure 7. Total elapsed time (in seconds) for the solution of the coupled MMPol/ddCOSMO equations as a function of the system size.
Achieving linear scaling in computational cost for a fully polarizable MM/Continuum embedding,
At each Self-Consistent Field (SCF) iteration devoted to the QM-part, one needs to solve the coupled polarization equation.
The PCM-energy contribution writes Es = 1 2 Z
R3 ⇢(r)W(r) dr,
where W solves: div(✏rW) = 0, in R3\@Ω, [W] = 0,
✏@W @n
@n ,
ε = εs ε = 1
The outside medium couples with all the spheres touching the interface. We will need the Fast Multipole Method (FMM) for achieving linear scaling.
Extensions available:
larizable force-field models (MD), charge densities from Quantum mechanics: semi-empirical methods, DFT, Hartree-Fock)
Special remarks:
research seminars, Roscoff summer school, etc.. Special thanks to the Institute of Scientific Computing (ICS) at Paris 6 who is pushing such interdisciplinary work.
how a mathematical idea has a large impact in an other (applied) community. Characteristics:
Domain decomposition for implicit solvation models E Cancès, Y Maday, B Stamm
A Fast Domain Decomposition Algorithm for Continuum Solvation Models: Energy and First Derivatives F Lipparini, B Stamm, E Cancès, Y Maday, B Mennucci
Quantum, Classical and Hybrid QM/MM Calculations in Solution: General Implementation of the ddCOSMO Linear Scaling Strategy
Scalable evaluation of the polarization energy and associated forces in polarizable molecular dynamics: I. towards massively parallel direct space computations.
Scalable Evaluation of Polarization Energy and Associated Forces in Polarizable Molecular Dynamics: II. Towards Massively Parallel Computations using Smooth Particle Mesh Ewald
Quantum calculations in solution for large to very large molecules: a new linear scaling QM/continuum approach F Lipparini, L Lagardère, G Scalmani, B Stamm, E Cancès, Y Maday, J-P Piquemal, M J Frisch, and B Mennucci,
Polarizable Molecular Dynamics in a Polarizable Continuum Solvent,
Achieving linear scaling in computational cost for a fully polarizable MM/Continuum embedding,
Domain decomposition for implicit solvation models E Cancès, Y Maday, B Stamm
A Fast Domain Decomposition Algorithm for Continuum Solvation Models: Energy and First Derivatives F Lipparini, B Stamm, E Cancès, Y Maday, B Mennucci
Quantum, Classical and Hybrid QM/MM Calculations in Solution: General Implementation of the ddCOSMO Linear Scaling Strategy
Scalable evaluation of the polarization energy and associated forces in polarizable molecular dynamics: I. towards massively parallel direct space computations.
Scalable Evaluation of Polarization Energy and Associated Forces in Polarizable Molecular Dynamics: II. Towards Massively Parallel Computations using Smooth Particle Mesh Ewald
Quantum calculations in solution for large to very large molecules: a new linear scaling QM/continuum approach F Lipparini, L Lagardère, G Scalmani, B Stamm, E Cancès, Y Maday, J-P Piquemal, M J Frisch, and B Mennucci,
Polarizable Molecular Dynamics in a Polarizable Continuum Solvent,
Achieving linear scaling in computational cost for a fully polarizable MM/Continuum embedding,
10 20 30 40 50 60 70 80 5000 10000 15000 20000 25000 30000 Wall Time (s) Number of atoms Elapsed time to solve COSMO equations Alanine chains Water clusters Hemoglobin subunits
TABLE I. Absolute timings (in seconds) for the various ddCOSMO and CSC tasks for a pure QM, pure Semiempirical (SE) and hybrid QM/MM computation. Keys as in figure 2. Task SE QM QM/MM dd CSC dd CSC dd CSC Phi 0.0530 0.0703 0.8134 0.6969 1.6740 1.3666 Psi 0.0002 1.9835 1.9852 X/Solve 0.0494 4.9328 0.0522 5.4797 0.6652 170.2033 S 0.0519 0.0557 0.7274 F1/Fock 0.0406 0.0405 0.9942 0.9400 1.8419 1.6333 F2 0.0001 3.1696 3.1350 Total 0.1951 5.0436 7.0687 7.1165 10.0288 173.2032
System of N atoms Permanent atomic multipoles : qi, µs,i, Θi Polarisability tensors (3*3) on atomic sites : αi Induced dipoles on atomic sites : µi, associated global induced dipole vector (3N) : µ Ei : electric field created by the permanent multipoles on site i , associated global vector (3N) E Ti,j operators, i 6= j : Ti,jµj : eletric field created by the dipole µj on site i 8i, µi = αi(Ei X
j6=j
Ti,jµj)
Matrix representation : T = α−1
1
T12 T13 . . . T1N T21 α−1
2
T23 . . . T2N T31 T32 ... . . . . . . . . . TN1 TN2 . . . α−1
N
Linear system to solve : Tµ = E
1e-07 1e-06 1e-05 0.0001 0.001 0.01 0.1 5 10 15 20 25 Norm of the increment/gradient JSOR CG JI/DIIS PCG Threshold
Figure 2: Computed IR spectra (black) and projected vibrational density of states (blue) on C and O atoms of all amide groups for a α-helix polypeptide model in vacuo 1000 2000 3000 4000
ω (cm
0.2 0.4 0.6 0.8 1
Intensity (arb. u.)
Figure 3: Computed IR spectra (black) and projected vibrational density of states (blue) on C and O atoms of all amide groups for a α-helix polypeptide model in water 1000 2000 3000 4000
ω (cm
0.2 0.4 0.6 0.8 1
Intensity (arb. u.) Figure 2: UV/Vis spectrum of the PCP light-harvesting complex as computed with the ZINDO method in vacuo and with ddCOSMO. The spectra have been normalized with respect to the maxima. The inset shows the experimental spectrum extracted from ref.[40]