Numerical analysis of the DMC method in a simple case. Benjamin - PowerPoint PPT Presentation

Numerical analysis of the DMC method in a simple case . Numerical analysis of the DMC method in a simple case. Benjamin Jourdain Joint work with Eric Canc` es, Mohamed El Makrini, Tony Leli` evre HIM, Bonn, April 7-11th 2008 Benjamin Jourdain (ENPC CERMICS) HIM, Bonn, April 7-11th 2008 1 / 34

Numerical analysis of the DMC method in a simple case . Structure of the talk Introduction 1 The Diffusion Monte Carlo method : a variance reduction technique 2 Analysis of the bias : the fixed-node approximation 3 Numerical implementation 4 Benjamin Jourdain (ENPC CERMICS) HIM, Bonn, April 7-11th 2008 2 / 34

Numerical analysis of the DMC method in a simple case . Introduction The Born-Oppenheimer approximation For most applications, systems of limited size (e.g. molecules) are described by M nuclei with electric charges z k ∈ N ∗ and positions ¯ x k ∈ R 3 , k ∈ { 1 , . . . , M } . slow variables → modelled by classical mechanics. N electrons with positions x i ∈ R 3 , i ∈ { 1 , . . . , N } and charge − 1. Very light particles. fast variables → modelled by (nonrelativistic) quantum mechanics. In the Born-Oppenheimer approximation, for electronic computations, the nuclei positions are supposed to be constant. Benjamin Jourdain (ENPC CERMICS) HIM, Bonn, April 7-11th 2008 3 / 34

Numerical analysis of the DMC method in a simple case . Introduction Electronic state To simplify, we leave the spin variables aside. Electronic state → modelled by the Hamiltonian N N M 1 z k 1 � � � � H = − 2 ∆ x i − x k | + | x i − ¯ | x i − x j | i = 1 i = 1 k = 1 1 ≤ i < j ≤ N where the positions ¯ x k of the nuclei are supposed to be fixed. Electronic ground state energy : � � � R 3 N | ψ | 2 = 1 E 0 = inf � ψ, H ψ � , ψ ∈ D ( H ) , (1) with � � N N � � L 2 ( R 3 ) L 2 ( R 3 ) D ( H ) = ψ ∈ , H ψ ∈ i = 1 i = 1 Benjamin Jourdain (ENPC CERMICS) HIM, Bonn, April 7-11th 2008 4 / 34

Numerical analysis of the DMC method in a simple case . Introduction Electronic state Antisymmetrized tensor product space N � L 2 ( R 3 ) = { ψ ∈ L 2 (( R 3 ) N ) : ∀ σ ∈ S n , ∀ x ∈ ( R 3 ) N , ψ ( x σ ) = ( − 1 ) ε ( σ ) ψ ( x ) } i = 1 with x σ = ( x σ ( 1 ) , . . . , x σ ( N ) ) for x = ( x 1 , . . . , x N ) and ε ( σ ) signature of σ . Justified by exchangeability of the electrons (in terms of the density | ψ | 2 : | ψ ( x σ ) | 2 = | ψ ( x ) | 2 ) Pauli’s exclusion principle for fermions ( | ψ | 2 vanishes when two positions x i are equal) Benjamin Jourdain (ENPC CERMICS) HIM, Bonn, April 7-11th 2008 5 / 34

Numerical analysis of the DMC method in a simple case . Introduction Ground state properties Ground state = element ψ 0 of D ( H ) which minimizes the energy (1) < ψ 0 , H ψ 0 > = E 0 with � ψ 0 � 2 = 1 . Theorem 1 When N ≤ Z = � M k = 1 z k (neutral molecule or positive ion), then H is self-adjoint in D ( H ) and there exists a ground state ψ 0 ( Zhislin 1960 ). Any ground state ψ 0 belongs to C θ ( R 3 N ) for θ ∈ ( 0 , 1 ) and to C ∞ ( R 3 N \ Γ) where Γ = { ( x 1 , . . . , x N ) ∈ ( R 3 ) N : ∃ i � = j s . t . x i = x j or ∃ i , k s . t . x i = ¯ x k } . Any ground state ψ 0 solves the Euler-Lagrange equation H ψ 0 = E 0 ψ 0 . (2) Benjamin Jourdain (ENPC CERMICS) HIM, Bonn, April 7-11th 2008 6 / 34

Numerical analysis of the DMC method in a simple case . Introduction The Tiling property For any continuous function ψ on R 3 N , we define U = R 3 N \ ψ − 1 ( 0 ) . For σ ∈ S N and x = ( x 1 , . . . , x N ) ∈ ( R 3 ) N , we set x σ = ( x σ ( 1 ) , . . . , x σ ( N ) ) . When ψ is antisymmetric, for any connected component C of U and ∀ σ ∈ S N , σ ( C ) = { x σ : x ∈ C} is also a connected component of U . Theorem 2 Any ground state ψ 0 satisfies the tiling property : for any connected component C of U 0 = R 3 N \ ψ − 1 0 ( 0 ) , U 0 = � σ ∈S N σ ( C ) ( Ceperley 91 ). Moreover, for any connected component C of U 0 , � 1 � � � � |∇ ψ | 2 + ψ 2 = 1 def V ψ 2 , ψ ∈ H 1 E 0 = E C = inf 0 ( C ) , . 2 C C C Benjamin Jourdain (ENPC CERMICS) HIM, Bonn, April 7-11th 2008 7 / 34

Numerical analysis of the DMC method in a simple case . Introduction Numerical computation of the ground state energy Difficult problem because dimension 3 N with N large, antisymmetry condition due to the fermionic nature of electrons Some numerical methods: Hartree-Fock methods (variational approximation : restriction of D ( H ) to Slater determinants det ( φ i ( x j )) ) , Density Functional Theory (Thomas-Fermi, Kohn-Sham) , Quantum Monte Carlo methods (Variational Monte Carlo, Diffusion Monte Carlo) . see E. Canc` es, M. Defranceschi, W. Kutzelnigg, C. Le Bris and Y. Maday, Computational Quantum Chemistry: a Primer , Handbook of Numerical Analysis, volume X (2003). Benjamin Jourdain (ENPC CERMICS) HIM, Bonn, April 7-11th 2008 8 / 34

Numerical analysis of the DMC method in a simple case . Introduction Trick : Schr¨ odinger in complex time Add a ficticious time variable : φ ( t , x ) = e − tH ψ I ( x ) solves � ∂ t φ = − H φ = 1 2 ∆ φ − V φ φ ( 0 , . ) = ψ I ( . ) ∈ D ( H ) with V ( x ) = − � N � M x k | + � z k 1 | x i − x j | . i = 1 k = 1 | x i − ¯ 1 ≤ i < j ≤ N If E 0 isolated single eigenvalue for H associated with eigenstate ψ 0 ( � ψ 0 � 2 = 1) and < ψ I , ψ 0 > � = 0, then φ ( t , x ) = e − tH ψ I ( x ) ∼ e − E 0 t < ψ I , ψ 0 > ψ 0 ( x ) ( t large) 1 ⇒ E 0 = − lim t log | φ ( t , x ) | . (3) t → + ∞ Benjamin Jourdain (ENPC CERMICS) HIM, Bonn, April 7-11th 2008 9 / 34

Numerical analysis of the DMC method in a simple case . Introduction Feynman-Kac interpretation If ( W t ) Brownian motion in R 3 N , for r ∈ [ 0 , t ] , � R r R r d r e − 0 V ( x + W s ) ds φ ( t − r , x + W r ) = e − 0 V ( x + W s ) ds ∇ x φ ( t − r , x + W r ) . dW r � � � − ∂ t φ + 1 + 2 ∆ φ − V φ ( t − r , x + W r ) dr . � �� � 0 Therefore � t R t R r e − 0 V ( x + W s ) ds ψ I ( x + W t ) = φ ( t , x )+ e − 0 V ( x + W s ) ds ∇ x φ ( t − r , x + W r ) . dW r 0 � 0 V ( x + W s ) ds � R t ψ I ( x + W t ) e − and φ ( t , x ) = E . By (3), � 0 V ( x + W s ) ds �� 1 � R t � ψ I ( x + W t ) e − � ⇒ E 0 = − lim t log � E � t → + ∞ Benjamin Jourdain (ENPC CERMICS) HIM, Bonn, April 7-11th 2008 10 / 34

Numerical analysis of the DMC method in a simple case . Introduction Variance reduction � 0 V ( x + W s ) ds �� 1 � R t � ψ I ( x + W t ) e − � E 0 = − lim t log � E � t → + ∞ Problem : large fluctuations of V in the exponential factor ⇒ variance too large. Need of variance reduction. Principle of the Diffusion Monte Carlo method (importance sampling) choose ψ I ∈ D ( H ) as close as possible to the ground state ψ 0 , modify the dynamics of the Brownian motion by adding the drift term ∇ ψ I ψ I , replace V by E L = H ψ I ψ I in the exponential factor (when ψ I = ψ 0 , E L = E 0 constant). Benjamin Jourdain (ENPC CERMICS) HIM, Bonn, April 7-11th 2008 11 / 34

Numerical analysis of the DMC method in a simple case . Introduction The Diffusion Monte Carlo method Yields very good results and is widely used in the chemistry community. J.C. Grossman, J. Chem. Phys., 117 (2002). Benjamin Jourdain (ENPC CERMICS) HIM, Bonn, April 7-11th 2008 12 / 34

Numerical analysis of the DMC method in a simple case . The Diffusion Monte Carlo method : a variance reduction technique Introduction 1 The Diffusion Monte Carlo method : a variance reduction technique 2 Analysis of the bias : the fixed-node approximation 3 Numerical implementation 4 Benjamin Jourdain (ENPC CERMICS) HIM, Bonn, April 7-11th 2008 13 / 34

Numerical analysis of the DMC method in a simple case . The Diffusion Monte Carlo method : a variance reduction technique The DMC method Let = < e − tH ψ I , H ψ I > E ( t ) def < e − tH ψ I , ψ I > By spectral decomposition, E ( t ) ∼ E 0 < ψ I , ψ 0 > 2 e − E 0 t → E 0 as t → + ∞ . < ψ I , ψ 0 > 2 e − E 0 t For f ( t , x ) = ψ I ( x ) e − tH ψ I ( x ) = ψ I ( x ) φ ( t , x ) and E L = H ψ I ψ I , E ( t ) = < f ( t , . ) , E L ( . ) > . < f ( t , . ) , 1 > Benjamin Jourdain (ENPC CERMICS) HIM, Bonn, April 7-11th 2008 14 / 34

Numerical analysis of the DMC method in a simple case . The Diffusion Monte Carlo method : a variance reduction technique Probabilistic interpretation of E ( t ) The function f ( t , x ) = ψ I ( x ) e − tH ψ I ( x ) = ψ I ( x ) φ ( t , x ) solves ∂ t f = 1 2 ∆ f − ∇ . ( bf ) − E L f , f ( 0 , . ) = ψ 2 I ( . ) (4) where b = ∇ ψ I ∆ ψ I ψ I and E L = H ψ I ψ I = − 1 ψ I + V . Assume � ψ I � 2 = 1. 2 Without the term − E L f , (4) Fokker-Planck equation for the density of X t solving dX t = dW t + b ( X t ) dt , X 0 ∼ ψ 2 I ( x ) dx . (5) With this term, h ( t , x ) defined by � � 0 E L ( X s ) ds � R t ∀ g : R 3 N → R , g ( X t ) e − R 3 N g ( x ) h ( t , x ) dx = E solves (4). Benjamin Jourdain (ENPC CERMICS) HIM, Bonn, April 7-11th 2008 15 / 34