Many-body system (in liquid phase) Second approximation: Do not neglect environment (solvent). Remarks: o DFT or HF for the entire system (solute and solvent) is computationally very demanding. o In addition, one computation is not representative by statistical mechanical 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. o Adding a simplified model to account for the solvent, and in particular for the interaction with the solute, is a feasible strategy: solvation models as part of a multi-scale strategy.
Many-body system (in liquid phase) Second approximation: Do not neglect environment (solvent). Remarks: o DFT or HF for the entire system (solute and solvent) is computationally very demanding. o In addition, one computation is not representative by statistical mechanical 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. o Adding a simplified model to account for the solvent, and in particular for the interaction with the solute, is a feasible strategy: solvation models as part of a multi-scale strategy. Example: QM/MM-approach (Nobelprize 2013) o Account for the solute by a quantum mechanical (QM) description. o Account for the solvent by a molecular mechanical (MM) description.
Solvation models o Most of the physical and chemical phenomena of interest in chemistry and biology take place in the liquid phase. o Solvent e ff ects play a crucial role in these processes.
Solvation models o Most of the physical and chemical phenomena of interest in chemistry and biology take place in the liquid phase. o Solvent e ff ects play a crucial role in these processes. solvent molecules solute
Solvation models o Most of the physical and chemical phenomena of interest in chemistry and biology take place in the liquid phase. o Solvent e ff ects play a crucial role in these processes. solvent molecules solute 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): o The electrostatic energy carried by the solute is modified in the presence of the solvent. o An extra term, called the electrostatic contribution to the solvation energy, and denoted here by E s , must be added to the electrostatic energy computed in vacuo .
Solvation models o Most of the physical and chemical phenomena of interest in chemistry and biology take place in the liquid phase. o Solvent e ff ects play a crucial role in these processes. solvent molecules solute Goal of this work: Compute the electro- static contribution (long-range interaction) to the solvation energy in an e ffi cient way. 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): o The electrostatic energy carried by the solute is modified in the presence of the solvent. o An extra term, called the electrostatic contribution to the solvation energy, and denoted here by E s , must be added to the electrostatic energy computed in vacuo .
Polarizable Continuum Model (PCM) Chatracteristics: o The solute is accommodated in a molec- solute ular cavity Ω and � is supported in Ω . o The environment (solvent) is represented by an infinite, continuum polarizable dielectric situated in R 3 \ Ω , described � s by its macroscopic dielectric suscepti- bility � s .
Polarizable Continuum Model (PCM) Chatracteristics: o The solute is accommodated in a molec- solute ular cavity Ω and � is supported in Ω . o The environment (solvent) is represented by an infinite, continuum polarizable dielectric situated in R 3 \ Ω , described � s by its macroscopic dielectric suscepti- bility � s . Energy contribution: E s = 1 � R 3 � ( r ) [ V ( r ) � Φ ( r )] d r , 2 where V solves � ( r � ) � in R 3 \ � Ω , div( � � V ) = 4 �� , | r � r � | d r � , Φ ( r ) = and R 3 [ V ] = 0 , on � Ω , � 1 in Ω � � = � ∂ V � � = 0 , on � Ω , in R 3 \ ¯ � � s Ω ∂ n �
dielectric situated in R 3 \ Ω , described � s by its macroscopic dielectric suscepti- Polarizable Continuum Model (PCM) bility � s . Energy contribution: E s = 1 � R 3 � ( r ) [ V ( r ) � Φ ( r )] d r , 2 where V solves � ( r � ) � in R 3 \ � Ω , div( � � V ) = 4 �� , | r � r � | d r � , Φ ( r ) = and R 3 [ V ] = 0 , on � Ω , � 1 in Ω � � = � ∂ V � � = 0 , on � Ω , in R 3 \ ¯ � � s Ω ∂ n � 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 E s = 1 � R 3 � ( r ) W ( r ) d r , 2 in R 3 \ � Ω , and W solves: div( � � W ) = 0 , [ W ] = 0 , on � Ω , � � �� W = ( � s � 1) � Φ on � Ω , � n , � n
dielectric situated in R 3 \ Ω , described � s by its macroscopic dielectric suscepti- Polarizable Continuum Model (PCM) bility � s . Energy contribution: E s = 1 � R 3 � ( r ) [ V ( r ) � Φ ( r )] d r , 2 where V solves � ( r � ) � in R 3 \ � Ω , div( � � V ) = 4 �� , | r � r � | d r � , Φ ( r ) = and R 3 [ V ] = 0 , on � Ω , � 1 in Ω � � = � ∂ V � � = 0 , on � Ω , in R 3 \ ¯ � � s Ω ∂ n � 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 E s = 1 � R 3 � ( r ) W ( r ) d r , 2 in R 3 \ � Ω , and W solves: div( � � W ) = 0 , [ W ] = 0 , on � Ω , � � �� W = ( � s � 1) � Φ on � Ω , � n , � n Rem: If � s = ∞ , then V = 0 on R 3 \ ¯ Ω .
COSMO Chatracteristics: COSMO (COnductor-like Screening MOdel) approximation: the dielec- tric is replaced by a conductor, i.e. ✏ s = ∞ . Also called C-PCM. � s
COSMO Chatracteristics: COSMO (COnductor-like Screening MOdel) approximation: the dielec- tric is replaced by a conductor, i.e. ✏ s = ∞ . Also called C-PCM. Energy contribution: � s C = 1 � E s 2 f ( � s ) ⇥ ( r ) W ( r ) d r , Ω � s − 1 where f ( � s ) = � s +0 . 5 is an empirical function of � s . W is solution to: Find W ∈ H 1 ( Ω ) such that � − ∆ W = 0 in Ω , W = − Φ on ⇤ Ω . Φ is the potential generated by ρ in vacuo : ρ ( r � ) � | r − r � | d r � . Φ ( r ) = R 3
COSMO Chatracteristics: COSMO (COnductor-like The di ffi culty is not to solve Screening MOdel) approximation: the dielec- Laplace’s equation, it is the tric is replaced by a conductor, i.e. ✏ s = ∞ . extremely complex geometry in- Also called C-PCM. cluding thousands of atoms. Energy contribution: � s Ω C = 1 � E s 2 f ( � s ) ⇥ ( r ) W ( r ) d r , Ω � s − 1 where f ( � s ) = � s +0 . 5 is an empirical function of � s . In addition, this energy needs W is solution to: Find W ∈ H 1 ( Ω ) such that to be computed a large num- ber of times! � − ∆ W = 0 in Ω , W = − Φ on ⇤ Ω . Φ is the potential generated by ρ in vacuo : ρ ( r � ) � | r − r � | d r � . Φ ( r ) = R 3
Notion of Molecular Cavity Note: a cavity of a molecule is an artificial construct. Di ff erent cavities are possible: o Van der Waals cavity o Solvent Accessible Surface (SAS) o Solvent Excluded Surface (SES)
Domain Decomposition Approach for the COSMO
Schwarz’ Domain Decomposition approach Decompose the domain into overlapping subdomains: Find W ∈ H 1 ( Ω ) such that M � Ω = Ω i � − ∆ W = 0 in Ω , i =1 W = − Φ on ∂ Ω .
Schwarz’ Domain Decomposition approach Decompose the domain into overlapping subdomains: Find W ∈ H 1 ( Ω ) such that M � Ω = Ω i � − ∆ W = 0 in Ω , i =1 W = − Φ on ∂ Ω .
Schwarz’ Domain Decomposition approach Decompose the domain into overlapping subdomains: Find W ∈ H 1 ( Ω ) such that M � Ω = Ω i � − ∆ W = 0 in Ω , i =1 W = − Φ on ∂ Ω . Formulate a local problem on each subdomain: Ω 3 For each i = 1 , . . . , M : find W i ∈ H 1 ( Ω i ) such that Ω 1 Ω 2 � − ∆ W i = 0 in Ω i , W i = g i on ∂ Ω i . Ω 4
Schwarz’ Domain Decomposition approach Decompose the domain into overlapping subdomains: Find W ∈ H 1 ( Ω ) such that M � Ω = Ω i � − ∆ W = 0 in Ω , i =1 W = − Φ on ∂ Ω . Formulate a local problem on each subdomain: Ω 3 For each i = 1 , . . . , M : find W i ∈ H 1 ( Ω i ) such that Ω 1 Ω 2 � − ∆ W i = 0 in Ω i , W i = g i on ∂ Ω i . Ω 4 What to impose on the g i ’s in order that W | Ω i = W i ?
Local boundary conditions Local boundary: Γ i = ∂ Ω i Γ e Ω 3 2 Splitting of local boundary into exterior and interior boundary: Γ i = Γ i i ∪ Γ e i Ω 2 Ω 1 Γ e Ω 4 2
Local boundary conditions Local boundary: Γ i = ∂ Ω i Γ e Ω 3 2 Splitting of local boundary into exterior and interior boundary: Γ i = Γ i i ∪ Γ e i Ω 2 Ω 1 Γ e Ω 4 2 We impose that on Γ e � − Φ i , � g i = � on Γ i Average of neighbour solutions. W i, N i , � where 1 � ∀ s ∈ Γ i W i, N ( s ) = W j ( s ) i . | N( i, s ) | j ∈ N( i, s ) Set of neighbouring/intersecting spheres at point s ∈ Γ i i .
Local boundary conditions Local boundary: Γ i = ∂ Ω i Γ e Ω 3 2 Splitting of local boundary into exterior and interior boundary: Γ i = Γ i i ∪ Γ e i Ω 2 Ω 1 Γ e Ω 4 2 We impose that on Γ e � − Φ i , � g i = � on Γ i Average of neighbour solutions. W i, N i , � where 1 � ∀ s ∈ Γ i W i, N ( s ) = W j ( s ) i . | N( i, s ) | j ∈ N( i, s ) Set of neighbouring/intersecting spheres at point s ∈ Γ i i . Then, solve in an iterative manner...
Domain Decomposition approach as integral equations For each i = 1 , . . . , M , we need to find W i ∈ H 1 ( Ω i ) such that � − ∆ W i = 0 in Ω i , W i = g i on ∂ Ω i .
Domain Decomposition approach as integral equations For each i = 1 , . . . , M , we need to find W i ∈ H 1 ( Ω i ) such that � − ∆ W i = 0 in Ω i , W i = g i on ∂ Ω i . The harmonic function W i can be represented by the single-layer potential σ i ( s � ) � W i ( s ) = ˜ | s − s � | d s � , S i σ i ( s ) := ∀ s ∈ Ω i , Γ i where σ i ∈ H � 1 2 ( Γ i ) solves the integral equation S i σ i = g i , on Γ i .
Domain Decomposition approach as integral equations For each i = 1 , . . . , M , we need to find W i ∈ H 1 ( Ω i ) such that � − ∆ W i = 0 in Ω i , W i = g i on ∂ Ω i . The harmonic function W i can be represented by the single-layer potential σ i ( s � ) � W i ( s ) = ˜ | s − s � | d s � , S i σ i ( s ) := ∀ s ∈ Ω i , Γ i where σ i ∈ H � 1 2 ( Γ i ) solves the integral equation S i σ i = g i , on Γ i . Here, S i : H � 1 1 2 ( Γ i ) → H 2 ( Γ i ) denotes the trace of the single-layer potential onto the surface Γ i , defined by σ ( s � ) � | s − s � | d s � . S i σ ( s ) = ∀ s ∈ Γ i , Γ i
Domain Decomposition approach as integral equations Thus, in the spirit of the domain decomposition approach, we solve for each i ∈ H − 1 2 ( Γ i ) such that i = 1 , . . . , M : find σ k Ω i S i σ k i = g k on Γ i , i , where on Γ e � − Φ (s) i , 1 2 ( Γ i ) . g k � g k i (s) = observe that i ∈ H j ∈ N( i, s) ˜ 1 S j σ k − 1 on Γ i � P (s) i , � j | N( i, s) | σ k j Ω j s s σ k i Ω i s σ k � s Ω �
Discretization
Discretization using Spherical Harmonics Ingredients of discretization: o Spherical harmonics: Maximum angular momentum N o Lebedev quadrature on S 2 : Number of integration points N g
Discretization using Spherical Harmonics Ingredients of discretization: o Spherical harmonics: Maximum angular momentum N o Lebedev quadrature on S 2 : Number of integration points N g Spherical harmonics: o Complete set of orthonormal functions for L 2 ( S 2 ). o S S 2 is a diagonal operator for the basis { Y lm } and: � Y lm ( s � ) 4 π | · − s � | d s � = S S 2 Y lm = 2 l + 1 Y lm S 2 o Appropriate scaling can be used to define basis functions on each Γ m . o We approximate σ k i by a truncated series of spherical harmonics � s − x i N l � � � Σ k ( Σ k i ( s ) = i ) lm Y lm | s − x i | l =0 m = � l
Discretization using Spherical Harmonics Ingredients of discretization: o Spherical harmonics: Maximum angular momentum N o Lebedev quadrature on S 2 : Number of integration points N g Spherical harmonics: o Complete set of orthonormal functions for L 2 ( S 2 ). o S S 2 is a diagonal operator for the basis { Y lm } and: � Y lm ( s � ) 4 π | · − s � | d s � = S S 2 Y lm = 2 l + 1 Y lm S 2 o Appropriate scaling can be used to define basis functions on each Γ m . o We approximate σ k i by a truncated series of spherical harmonics � s − x i N l � � � Σ k ( Σ k i ( s ) = i ) lm Y lm | s − x i | l =0 m = � l
Discretization using truncated Spherical Harmonics At each iteration k and for each i = 1 , . . . , M : find Σ k i 2 span { Y lm } such that S i Σ k G k S i Σ k i = Π N G k � � � � Γ i = , i , Y lm i , Y lm Γ i , i where on Γ e � � Φ ( s ) i , 1 G k G k 2 ( Γ i ) . � i ( s ) = observe that i 62 H j ∈ N( i, s ) ˜ S j Σ k − 1 1 on Γ i � P ( s ) i , � j | N( i, s ) |
Discretization using truncated Spherical Harmonics At each iteration k and for each i = 1 , . . . , M : find Σ k i 2 span { Y lm } such that S i Σ k G k S i Σ k i = Π N G k � � � � Γ i = , i , Y lm i , Y lm Γ i , i where on Γ e � � Φ ( s ) i , 1 G k G k 2 ( Γ i ) . � i ( s ) = observe that i 62 H j ∈ N( i, s ) ˜ S j Σ k − 1 1 on Γ i � P ( s ) i , � j | N( i, s ) | Remark 1: o If the expansion of G k i in terms of spherical harmonics is know, i.e. N m X X Π N G k ( G k i = i ) lm Y lm , l =0 m = − l then i ) lm = (2 l + 1) ( Σ k ( G k i ) lm . 4 π R i
Discretization using truncated Spherical Harmonics At each iteration k and for each i = 1 , . . . , M : find Σ k i 2 span { Y lm } such that S i Σ k G k S i Σ k i = Π N G k � � � � Γ i = , i , Y lm i , Y lm Γ i , i where on Γ e � � Φ ( s ) i , 1 G k G k 2 ( Γ i ) . � i ( s ) = observe that i 62 H j ∈ N( i, s ) ˜ S j Σ k − 1 1 on Γ i � P ( s ) i , � j | N( i, s ) | Remark 2: o It remains to compute the coe ffi cients ( G k i ) lm at each iteration. By definition, the coe ffi cients of G k i are defined by Z G k i ( s ) Y lm ( s ) d s . Γ i We replace the integral by the numerical integration to obtain N g X ( G k ω n G k i ) lm = i ( s n ) Y lm ( s n ) d s . n =1
Σ k − 1 One iteration j s n Ω j s n y n Σ k At each iteration k and for each sphere i = 1 , . . . , M , i Ω i we need to access the neighbouring potential at the interior integration point s n ∈ Γ i i : s n Ω � 1 Σ k − 1 ˜ X S j Σ k − 1 G k i (s n ) = (s n ) . � j | N( i, s n ) | j ∈ N( i, s n ) Compute y n = s n − x j | s n − x j | N m X X Σ k − 1 ( Σ k − 1 (y n ) = ) lm Y lm (y n ) , j j l =0 m = − l N m 2 l + 1 X X S j Σ k − 1 ( Σ k − 1 (y n ) = ) lm Y lm (y n ) , j j 4 π R j l =0 m = − l ◆ l 2 l + 1 N m ✓ | s n − x j | ˜ X X S j Σ k − 1 ( Σ k − 1 (s n ) = ) lm Y lm (y n ) . j j 4 π R j R j l =0 m = − l
Σ k − 1 One iteration j s n Ω j s n y n Σ k At each iteration k and for each sphere i = 1 , . . . , M , i Ω i we need to access the neighbouring potential at the interior integration point s n ∈ Γ i i : s n Ω � 1 Σ k − 1 ˜ X S j Σ k − 1 G k i (s n ) = (s n ) . � j | N( i, s n ) | j ∈ N( i, s n ) Compute y n = s n − x j | s n − x j | N m X X Σ k − 1 ( Σ k − 1 (y n ) = ) lm Y lm (y n ) , j j l =0 m = − l N m 2 l + 1 X X S j Σ k − 1 ( Σ k − 1 (y n ) = ) lm Y lm (y n ) , j j 4 π R j l =0 m = − l ◆ l 2 l + 1 N m ✓ | s n − x j | ˜ X X S j Σ k − 1 ( Σ k − 1 (s n ) = ) lm Y lm (y n ) . j j 4 π R j R j l =0 m = − l
Σ k − 1 One iteration j s n Ω j s n y n Σ k At each iteration k and for each sphere i = 1 , . . . , M , i Ω i we need to access the neighbouring potential at the interior integration point s n ∈ Γ i i : s n Ω � 1 Σ k − 1 ˜ X S j Σ k − 1 G k i (s n ) = (s n ) . � j | N( i, s n ) | j ∈ N( i, s n ) Compute y n = s n − x j | s n − x j | N m X X Σ k − 1 ( Σ k − 1 (y n ) = ) lm Y lm (y n ) , j j l =0 m = − l N m 2 l + 1 X X S j Σ k − 1 ( Σ k − 1 (y n ) = ) lm Y lm (y n ) , j j 4 π R j l =0 m = − l ◆ l 2 l + 1 N m ✓ | s n − x j | ˜ X X S j Σ k − 1 ( Σ k − 1 (s n ) = ) lm Y lm (y n ) . j j 4 π R j R j l =0 m = − l
Σ k − 1 One iteration j s n Ω j s n y n Σ k At each iteration k and for each sphere i = 1 , . . . , M , i Ω i we need to access the neighbouring potential at the interior integration point s n ∈ Γ i i : s n Ω � 1 Σ k − 1 ˜ X S j Σ k − 1 G k i (s n ) = (s n ) . � j | N( i, s n ) | j ∈ N( i, s n ) Compute y n = s n − x j | s n − x j | N m X X Σ k − 1 ( Σ k − 1 (y n ) = ) lm Y lm (y n ) , j j l =0 m = − l N m 2 l + 1 X X S j Σ k − 1 ( Σ k − 1 (y n ) = ) lm Y lm (y n ) , j j 4 π R j l =0 m = − l ◆ l 2 l + 1 N m ✓ | s n − x j | ˜ X X S j Σ k − 1 ( Σ k − 1 (s n ) = ) lm Y lm (y n ) . j j 4 π R j R j l =0 m = − l
Example of caffeine
Example of caffeine
Global formulation The proposed iterative scheme is a way to solve the global block-sparse non- symmetric linear system: L 11 L 1 M f 1 Σ 1 · · · . . . . ... . . . . = . . . . L M 1 L MM f M Σ M · · · 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.
Global formulation The proposed iterative scheme is a way to solve the global block-sparse non- symmetric linear system: L 11 L 1 M f 1 Σ 1 · · · . . . . ... . . . . = . . . . L M 1 L MM f M Σ M · · · 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. 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 di ff erent techniques: pre- conditioning, DIIS (Direct Inversion in the Iterative Subspace), .... The DIIS post-processing works really good in practice.
Global formulation The proposed iterative scheme is a way to solve the global block-sparse non- symmetric linear system: L 11 L 1 M f 1 Σ 1 · · · . . . . ... . . . . = . . . . L M 1 L MM f M Σ M · · · 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. In consequence: o Well-suited for parallel implementation o Linear memory requirement with respect to the number of atoms, compared to a BEM-implementation (dense solution matrix). o Linear scaling in number of operations with respect to M the number of atoms if number of iterations is constant (see next slide).
Linear scaling? Counter-example: N spheres It takes at least O ( N ) iterations until the Dirichlet BC’s arrive at the center!
Linear scaling? Counter-example: N spheres It takes at least O ( N ) iterations until the Dirichlet BC’s arrive at the center! This is a slightly di ff erent point of view than in an usual DD-paradigm where the domain Ω is constant and the number of subdomains to cover Ω is increased.
Linear scaling? Counter-example: N spheres It takes at least O ( N ) iterations until the Dirichlet BC’s arrive at the center! This is a slightly di ff erent 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).
Linear scaling? Counter-example: N spheres It takes at least O ( N ) iterations until the Dirichlet BC’s arrive at the center! This is a slightly di ff erent 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).
Regularization 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.
Regularization 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. o Integration points suddenly become exterior points. o Since the function on Γ e � − Φ (s n ) i , G k � i (s n ) = j ∈ N( i, s n ) ˜ S j Σ k − 1 1 on Γ i � � (s n ) i , � | N( i, s n ) | j is not continuous (only its exact counterpart), this implies a non-smooth dependency.
Regularization 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. o Integration points suddenly become exterior χ 0 points. o Since the function on Γ e � − Φ (s n ) i , G k � i (s n ) = j ∈ N( i, s n ) ˜ S j Σ k − 1 1 on Γ i � � (s n ) i , � | N( i, s n ) | j is not continuous (only its exact counterpart), this implies a non-smooth dependency. We now write χ 0 ˜ � S j Σ k − 1 G k i (s n ) = − (1 − χ 0 ) Φ (s n ) + (s n ) , j | N( i, s n ) | j ∈ N( i, s n )
Regularization 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. η o Integration points suddenly become exterior χ η points. o Since the function on Γ e � − Φ (s n ) i , G k � i (s n ) = j ∈ N( i, s n ) ˜ S j Σ k − 1 1 on Γ i � � (s n ) i , � | N( i, s n ) | j is not continuous (only its exact counterpart), this implies a non-smooth dependency. We now write χ 0 ˜ � S j Σ k − 1 G k i (s n ) = − (1 − χ 0 ) Φ (s n ) + (s n ) , j | N( i, s n ) | j ∈ N( i, s n ) and replace it by χ η ˜ � S j Σ k − 1 G k i (s n ) = − (1 − χ η ) Φ (s n )+ (s n ) , j | N( i, s n ) | j ∈ N( i, s n )
Force computations o The solvation energy is not the only desired quantity. o The derivative of the energy with respect to the molecular coordinates (posi- tion of each atom) is an important quantity for geometry optimization or time-dependent simulations.
Force computations o The solvation energy is not the only desired quantity. o The derivative of the energy with respect to the molecular coordinates (posi- tion of each atom) is an important quantity for geometry optimization or time-dependent simulations. Let x q denote the coordinate of the q -th nuclei. Then we would like to compute F q = �r q E s C where r q denotes the derivative with respect to the position x q .
Force computations o The solvation energy is not the only desired quantity. o The derivative of the energy with respect to the molecular coordinates (posi- tion of each atom) is an important quantity for geometry optimization or time-dependent simulations. Let x q denote the coordinate of the q -th nuclei. Then we would like to compute F q = �r q E s C where r q denotes the derivative with respect to the position x q . The energy can be written as M X E s C = h Ψ i , Σ i i i =1 for some given multi-polar expansion Ψ and where N l X X h a, b i = ( a ) lm ( b ) lm . l =0 m = − l
Force computations: analytical derivatives Consider the global system by M X L Σ = f L ji Σ i = f j , j = 1 , . . . , M. , i =1 Take the derivative to obtain: M M M X X X r q f j = r q L ji Σ i = r q ( L ji ) Σ i + L ji r q ( Σ i ) . i =1 i =1 i =1
Force computations: analytical derivatives Consider the global system by M X L Σ = f L ji Σ i = f j , j = 1 , . . . , M. , i =1 Take the derivative to obtain: M M M X X X r q f j = r q L ji Σ i = r q ( L ji ) Σ i + L ji r q ( Σ i ) . i =1 i =1 i =1 Thus M M X X L ji r q ( Σ i ) = r q f j � r q ( L ji ) Σ i , i =1 i =1 | {z } (closed-form expression possible) =( h q ) j , and globally M X r q Σ = L − 1 h q L − 1 L ( r q Σ ) = h q r q Σ j = ji ( h q ) i , j = 1 , . . . , M. , , i =1
M M M X X X r q f j = r q L ji Σ i = r q ( L ji ) Σ i + L ji r q ( Σ i ) . Force computations: analytical derivatives i =1 i =1 i =1 Thus M M X X L ji r q ( Σ i ) = r q f j � r q ( L ji ) Σ i , i =1 i =1 | {z } (closed-form expression possible) =( h q ) j , and globally M X r q Σ = L − 1 h q L − 1 L ( r q Σ ) = h q r q Σ j = ji ( h q ) i , j = 1 , . . . , M. , , i =1 Then, M M M * M + X X X X L − 1 F q = �r q h Ψ j , Σ j i = � h Ψ j , r q Σ j i = � ji ( h q ) i Ψ j , j =1 j =1 j =1 i =1 * M + M M X X X L − T = � ji Ψ j , ( h q ) i = � h s i , ( h q ) i i i =1 j =1 i =1 where s solves the adjoint problem L T s = Ψ .
Numerical results
Two spheres: influence of N η = 0: d=1.25 d=1.25 d=1.25 35 d=1.5 d=1.5 d=1.5 d=1.75 d=1.75 d=1.75 0.01 30 number of iterations 0.1 relative error relative error 25 0.001 20 15 0.01 10 0.0001 1 10 100 1 10 100 0 10 20 30 40 50 N N N H − 1 / 2 σ η = 0 . 1: d=1.25 d=1.25 d=1.25 20 0.01 d=1.5 d=1.5 d=1.5 d=1.75 d=1.75 d=1.75 18 0.1 0.001 number of iterations 16 relative error relative error 0.0001 14 12 0.01 0.00001 10 1 x 10 -6 8 0 10 20 30 40 50 0 10 20 30 40 50 0 10 20 30 40 50 N N N H − 1 / 2 σ
Comparison with the state-of-the-art CSC: Continuous Surface Charge by York and Karplus • smooth energy profiles • not systematically improvable • bad conditioning: unfortunately FMM doesn’t help here! System Atoms CSC CSC/FMM ddCOSMO Vancomycin 377 20” 43” < 1” Hiv-1-GP41 530 1’26” 57” < 1” l-Plectasin 567 1’35” 1’17” < 1” Glutaredoxin 1277 8’54” 3’02” < 1” Glutaredoxin* 28’43” - - UBCH5B 2360 94’ 6’18” 1” Carboxylase 6605 777’ 24’42” 3” Table: Timings for the solution of the C-PCM/COSMO linear equations on 2*Xeon E2560 2GHz.. *The computation was repeated without the N 2 storage for CSC.
Hemoglobin subunits: globular molecules Hemoglobin and its subunits 22 350 direct system adjoint system 20 Forces (matrix part) Memory - cavity and scratch 300 Memory - Geometrical parameters 18 16 Memory (Integer MegaWords) 250 14 Wall Time (s) 200 12 10 150 8 100 6 4 50 2 0 0 2000 3000 4000 5000 6000 7000 8000 9000 Number of atoms
Alanine chains: “linear molecules” Alanine chains 90 1600 direct system adjoint system Forces (matrix part) 80 1400 Memory - cavity and scratch Memory - Geometrical parameters 70 1200 Memory (Integer MegaWords) 60 1000 Wall Time (s) 50 800 40 600 30 400 20 200 10 0 0 0 5000 10000 15000 20000 25000 30000 35000 40000 Number of atoms
Acceleration techniques 160 140 120 Number of iterations 100 Jacobi iterations 80 DIIS iterations 60 40 20 0 0 10000 20000 30000 40000 50000 60000 70000 Number of spheres
Parallel implementation using MPI on MeSU-alpha (ICS) 350 Ideal Speedup 300 250 performance gain 200 150 100 50 0 0 50 100 150 200 250 300 350 number of cores
DNA-Origami DNA Origami (PDB ref. 2YMG) 74708 Atoms N g = 302, N = 10 Timings and memory requirement: 5.5GB RAM needed for the computation 147s for the linear system, 121 for the adjoint one 39 s for the (PCM related) contributions to the forces
This opens the door to a new universe: Applications
Coupling with polarizable force-field models 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 F. Lipparini, L. Lagardère, Ch. Raynaud, B. Stamm, M. Schnieders, B. Mennucci, P. Ren,Y. Maday, J.-P. Piquemal, J. Chem. Theory Comput., 2015 Molecular dynamics simulation of ubiquitin in water (100ps), with the solvent described by a polarizable continuum, using a time step of 1fs
Coupling with polarizable force-field models 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 F. Lipparini, L. Lagardère, Ch. Raynaud, B. Stamm, M. Schnieders, B. Mennucci, P. Ren,Y. Maday, J.-P. Piquemal, J. Chem. Theory Comput., 2015 Molecular dynamics simulation of ubiquitin in water (100ps), with the solvent described by a polarizable continuum, using a time step of 1fs
Coupling with polarizable force-field models 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 F. Lipparini, L. Lagardère, Ch. Raynaud, B. Stamm, M. Schnieders, B. Mennucci, P. Ren,Y. Maday, J.-P. Piquemal, J. Chem. Theory Comput., 2015 Molecular dynamics simulation of ubiquitin in water (100ps), with the solvent described by a polarizable continuum, using a time step of 1fs Table 1. Stability of NVE Simulations, as a Function of the Convergence Threshold for the Coupled Polarization Equations and of the Time Step (AMBER/Wang Variational Force Field) a 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 a For each combination, the short-time average fl uctuation (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 of the Time Step (AMOEBA Force Field) a 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 Dual Xeon Model E5-2650 cluster 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 node (16 cores, 2 GHz) equipped with a For each combination, the short-time average fl uctuation (STF), the long-time drift (LTD), and the average computational time per time step 64GB of DDR3 memory (on TACC) (ACT) are reported. All the energies are given in terms of kcal/mol, all timings are given in seconds.
Coupling with polarizable force-field models 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 F. Lipparini, L. Lagardère, Ch. Raynaud, B. Stamm, M. Schnieders, B. Mennucci, P. Ren,Y. Maday, J.-P. Piquemal, J. Chem. Theory Comput., 2015 Molecular dynamics simulation of ubiquitin in water (100ps), with the solvent described by a polarizable continuum, using a time step of 1fs Table 1. Stability of NVE Simulations, as a Function of the Convergence Threshold for the Coupled Polarization Equations and of the Time Step (AMBER/Wang Variational Force Field) a 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 a For each combination, the short-time average fl uctuation (STF), the long-time drift (LTD), and the average computational time per time step Figure 1. Scaling of the MPI code for the coupled induced dipoles/ (ACT) are reported. All the energies are given in terms of kcal/mol, all timings are given in seconds. ddCOSMO problem, with respect to a single-node (16 cores) computation. Table 2. Stability of NVE Simulations, as a Function of the Convergence Threshold for the Coupled Polarization Equations and of the Time Step (AMOEBA Force Field) a 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 Dual Xeon Model E5-2650 cluster 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 node (16 cores, 2 GHz) equipped with a For each combination, the short-time average fl uctuation (STF), the long-time drift (LTD), and the average computational time per time step 64GB of DDR3 memory (on TACC) (ACT) are reported. All the energies are given in terms of kcal/mol, all timings are given in seconds.
Multi-scale model: 3 layer QM/MM/Continuum 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.
Multi-scale model: 3 layer QM/MM/Continuum 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. Achieving linear scaling in computational cost for a fully polarizable MM/Continuum embedding, S. Caprasecca, S. Jurinovich, L. Lagardère, B. Stamm, F. Lipparini, J. Chem. Theory Comput., 2015
Perspectives
Solvent Excluded Surface (Thesis of Chaoyu Quan)
Solvent Excluded Surface (Thesis of Chaoyu Quan)
Solvent Excluded Surface (Thesis of Chaoyu Quan)
PCM (and not COSMO) The PCM-energy contribution writes E s = 1 Z R 3 ⇢ ( r ) W ( r ) d r , 2 where W solves: in R 3 \ @ Ω , div( ✏ r W ) = 0 , [ W ] = 0 , on @ Ω , � ✏@ W = ( ✏ s � 1) @ Φ on @ Ω . @ n , @ n The outside medium couples with all the spheres touching the interface. We will need the Fast Multipole Method (FMM) for achieving linear scaling. ε = ε s Ω ε = 1
Conclusions Characteristics: o The ddCOSMO is robust, smooth (energy profile) and incredibly fast. o The approximation setting is fully under control: systematically improvable. o Linear scaling for all molecules (including globular ones where FMM is performing poorly). Extensions available: o Extension to more general charge distributions (beyond point-charges: Multipoles from po- larizable force-field models (MD), charge densities from Quantum mechanics: semi-empirical methods, DFT, Hartree-Fock) o Implementation in large codes: Gaussian and Tinker. Special remarks: o This project shows how to build bridges between di ff erent disciplines. o It needs some e ff ort and time until chemists and mathematicians can talk to each other: common research seminars, Rosco ff summer school, etc.. Special thanks to the Institute of Scientific Computing (ICS) at Paris 6 who is pushing such interdisciplinary work. o This algorithm is about to be implemented in standard chemistry codes: a wonderful example how a mathematical idea has a large impact in an other (applied) community.
References Domain decomposition for implicit solvation models E Cancès, Y Maday, B Stamm ddCOSMO J. Chem. Phys., 2013 A Fast Domain Decomposition Algorithm for Continuum Solvation Models: Energy and First Derivatives F Lipparini, B Stamm, E Cancès, Y Maday, B Mennucci J. Chem. Theory Comput., 2013 Quantum, Classical and Hybrid QM/MM Calculations in Solution: General Implementation of the ddCOSMO Linear Scaling Strategy F. Lipparini, G. Scalmani, L. Lagardère, B. Stamm, E. Cancès, Y. Maday, J.-P. Piquemal, M. Frisch and B. Mennucci, J. Chem. Phys., 2014 Scalable evaluation of the polarization energy and associated forces in polarizable molecular dynamics: I. towards massively parallel direct space computations. F. Lipparini, L. Lagardère, B. Stamm, E. Cancès, M. Schnieders, P. Ren, Y. Maday, J.-P. Piquemal, J. Chem. Theory Comput., 2014 PFF Scalable Evaluation of Polarization Energy and Associated Forces in Polarizable Molecular Dynamics: II. Towards Massively Parallel Computations using Smooth Particle Mesh Ewald L. Lagardère, F. Lipparini, E. Polack, B. Stamm, E. Cancès, M. Schnieders, P. Ren, Y. Maday, J.-P. Piquemal, J. Chem. Theory Comput., 2015 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, J. Phys. Chem. Lett., 2014 Coupling Polarizable Molecular Dynamics in a Polarizable Continuum Solvent, F. Lipparini, L. Lagardère, Ch. Raynaud, B. Stamm, M. Schnieders, B. Mennucci, P. Ren,Y. Maday, J.-P. Piquemal, J. Chem. Theory Comput., 2015 Achieving linear scaling in computational cost for a fully polarizable MM/Continuum embedding, S. Caprasecca, S. Jurinovich, L. Lagardère, B. Stamm, F. Lipparini, J. Chem. Theory Comput., 2015
References Domain decomposition for implicit solvation models E Cancès, Y Maday, B Stamm ddCOSMO J. Chem. Phys., 2013 A Fast Domain Decomposition Algorithm for Continuum Solvation Models: Energy and First Derivatives F Lipparini, B Stamm, E Cancès, Y Maday, B Mennucci J. Chem. Theory Comput., 2013 Quantum, Classical and Hybrid QM/MM Calculations in Solution: General Implementation of the ddCOSMO Linear Scaling Strategy F. Lipparini, G. Scalmani, L. Lagardère, B. Stamm, E. Cancès, Y. Maday, J.-P. Piquemal, M. Frisch and B. Mennucci, J. Chem. Phys., 2014 Scalable evaluation of the polarization energy and associated forces in polarizable molecular dynamics: I. towards massively parallel direct space computations. F. Lipparini, L. Lagardère, B. Stamm, E. Cancès, M. Schnieders, P. Ren, Y. Maday, J.-P. Piquemal, J. Chem. Theory Comput., 2014 PFF Scalable Evaluation of Polarization Energy and Associated Forces in Polarizable Molecular Dynamics: II. Towards Massively Parallel Computations using Smooth Particle Mesh Ewald L. Lagardère, F. Lipparini, E. Polack, B. Stamm, E. Cancès, M. Schnieders, P. Ren, Y. Maday, J.-P. Piquemal, J. Chem. Theory Comput., 2015 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, J. Phys. Chem. Lett., 2014 Coupling Polarizable Molecular Dynamics in a Polarizable Continuum Solvent, F. Lipparini, L. Lagardère, Ch. Raynaud, B. Stamm, M. Schnieders, B. Mennucci, P. Ren,Y. Maday, J.-P. Piquemal, J. Chem. Theory Comput., 2015 Achieving linear scaling in computational cost for a fully polarizable MM/Continuum embedding, S. Caprasecca, S. Jurinovich, L. Lagardère, B. Stamm, F. Lipparini, J. Chem. Theory Comput., 2015 Thank you for your attention!
Comparison Elapsed time to solve COSMO equations 80 Alanine chains Water clusters Hemoglobin subunits 70 60 50 Wall Time (s) 40 30 20 10 0 0 5000 10000 15000 20000 25000 30000 Number of atoms
Coupling with QM, SE, hybrid QM/MM 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
Scalable evaluation of the polarization energy System of N atoms Permanent atomic multipoles : q i , µ s,i , Θ i Polarisability tensors (3*3) on atomic sites : α i Induced dipoles on atomic sites : µ i , associated global induced dipole vector (3N) : µ E i : electric field created by the permanent multipoles on site i , associated global vector (3N) E T i,j operators, i 6 = j : � T i,j µ j : eletric field created by the dipole µ j on site i X 8 i, µ i = α i ( E i � T i,j µ j ) j 6 = j Matrix representation : α − 1 T 12 T 13 . . . T 1 N 1 α − 1 T 21 T 23 . . . T 2 N 2 ... T = T 31 T 32 . . . . . . . . . α − 1 T N 1 T N 2 . . . N Linear system to solve : T µ = E
Scalable evaluation of the polarization T is symetrical positive definite iterative methods : JSOR, Jacobi+DIIS, Conjugate gradient (with preconditionner) 0.1 JSOR CG JI/DIIS PCG Threshold 0.01 Norm of the increment/gradient 0.001 0.0001 1e-05 1e-06 1e-07 0 5 10 15 20 25
Coupling with Polarizable Force-Fields (PFF) Coupled Equations : L σ = g (0) + g ( µ ) L ∗ s = ψ (0) + ψ ( µ ) T µ = E (0) + E ( σ , s )
Recommend
More recommend
