1 / 45

A quantum Monte Carlo study of noncovalent interactions

A quantum Monte Carlo study of noncovalent interactions. Kenta Hongo. Assistant Professor, School of Information Science, JAIST (Prof. Ryo Maezono’s group). in collaboration with: A. Aspuru-Guzik, M.A. Watson, R.S. Sánchez-Carrera , T. Iitaka (for molecular crystals)

gizela
Download Presentation

A quantum Monte Carlo study of noncovalent interactions

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. A quantum Monte Carlo study of noncovalent interactions Kenta Hongo Assistant Professor, School of Information Science, JAIST (Prof. Ryo Maezono’s group) in collaboration with: A. Aspuru-Guzik, M.A. Watson, R.S. Sánchez-Carrera, T. Iitaka (for molecular crystals) R. Maezono, Nguyen Thanh Cuong (for DNA stacking) Quantum Monte Carlo in the Apuan Alps VII @ TTI, Vallico Sotto, Tuscany, Italy, 08/03/2012

  2. Outline • Motivation of our studies • Noncovalent interactions • DMC and DFT studies of • Polymorphism in para-diiodobenzene (DIB) • Base-pair-step stacking in B-DNA (AT-AT)

  3. Noncovalent interactions:Challenges for theoretical methods accuracy: • The energy scale of the interaction is quite small (order of 1 mHa) • Noncovalent bond: hydrogen / dispersion-dominated complexes/mixture of them computational cost: • Week interactions between complicated molecules • containing hundreds number of electrons… • Accurate quantum chemistry methods (eg. CCSD(T)) are not routinely applicable • On the other hand, commonly-used DFTs mostly fail to describe noncovalent interactions. • A few QMC studies have been performed so far • eg. For the S22 benchmark: Martin Korth, et. al., J. Phys. Chem. A, 112, 2104, (2008). Noncovalent systems are good challenges to theoretical methods.

  4. Molecular crystals Many applications to… pharmaceutical materials, semiconductors, superconductors, etc. Picene DNTT Aspirin R. Mitsuhashi, et al, Nature 464, 76 (2010). R. S. Sánchez-Carrera, et al, J. Phys. Chem. C 114, 2334 (2010). P. Veshweshwar, et al, J.Am. Chem. Soc. 127, 16802 (2005)

  5. Polymorphism in molecular crystals • The existence of more than one form of a compound • Each polymorph has different physical and chemical properties e.g. effects of polymorphism on charge transport in organic semiconductors Fluorinated 5,11-bis(triethylsilylethynyl)antradihiophene (diF TES ADT) O. D. Jurchescu, et al, Phys. Rev. B 80, 085201 (2009).

  6. para-diiodobenzene (p-DIB) π-conjugated organic semiconductors as (opt)electronic devices -Typical room-temperature charge-carrier mobility ≈ 1 cm2/(V s) para-diiodobenzene: Room-temperature hole mobility = 12 cm2/(V s) L. M. Schwartz and J. F. Horning, Mol. Cryst. 2, 379 (1967). R. S. Sánchez-Carrera, et al. Chem. Mater. 20, 5832 (2008).

  7. Polymorphism in p-DIB known to exhibit polymorphism (from experiment): α-phase (Pbca) β-phase (Pccn) Transition from α- to β-phase occurs at 326 K (The α-phase is more stable than the β-phase at zero temperature)

  8. Previous study of DIB • A. Brillante, et al., JACS 127, 3038 (2005) DFT calculation: Exchange-correlation functional = BLYP functional Pseudopotential = Troullier-Martins The α-phase is less stable by ≈0.002 eV/atom than the β-phase at 0 K. Contrary to experimental results! Good challenge for DMC!

  9. Collaborators & HARIKEN project Dr. Toshiaki Iitaka (RIKEN in Japan) Prof. Alán Aspuru-Guzik (AG group) Dr. Mark A. Watson (Princeton University) Dr. Roel S. Sánchez-Carrera (Former member) = + (hurricane) Harvard RIKEN The first application of QMC to molecular crystal systems: K. Hongo, M. A. Watson, R. S. Sánchez-Carrera, T. Iitaka, A. Aspruru-Guzik, J. Phys. Chem. Lett. 1, 1789 (2010).

  10. Computational details: QMC • DMC calc. using QMCPACK • Pseudo potential = Trail-Needs type (CASINO library) • The total number of electrons = 168 • Trial wavefunction = Slater-Jastrow type • Slater determinant: LDA (cut off energy = 70 Ry) • Jastrow factor: one- & two-body terms (6 & 4 parameters) • 1x1x1 simulation cell with a finite size correction • KZK scheme (Kwee, Zhang, Krakauer, PRL, 100, 126404, 2008) • computational conditions of DMC • # of walkers = 16384, # of MC steps = 3.2×107

  11. Relative stability between the two phases 68 29 25 B3LYP+D DMC LDA ΔE = E(α) – E(β) (<0: Experiment) PBE B3LYP PW91 -48+/-24 -123 -155 • The DMC result is consistent with experiment, i.e., a negative ΔE. • GGAs and B3LYP predict a positive ΔE. • LDA and B3LYP+D correctly predict the relative stability, but strongly overestimate the magnitude of ΔE, compared to DMC. • It is well-known that LDA frequently overbinds. • The Grimme semiempirical dispersion correction is important.

  12. Biological systems Adenine (A) Guanine (G) Cytosine (C) Thymine (T) B-DNA • Noncovalent bonding: Hydrogen bond, dispersion, and their mixture…

  13. Present study: Base-pair-step stacking in B-DNA Adenine Adenine Adenine Thymine Thymine Thymine d: distance between 2 layers d: distance between 2 layers d: distance between 2 layers Adenine Adenine Adenine Thymine Thymine Thymine We applied DMC and DFT’s to this system…

  14. Computational details • DMC calculations using CASINO • Burkatzki-Fillippi-Dolg PPs with associated basis sets • VDZ and VTZ basis sets • Trial wavefunction = Slater-Jastrow type • LDA, PBE, B3LYP, HF trial nodes (generated by Gaussian 09) • Jastrow factor: one-, two-, and three-body terms (up to 8the poly.) • computational conditions of DMC • # of walkers = 1280, # of MC steps = 105 • DFT calculations using Gaussian 09 • recently developed various functionals examined • Hybrid, DFT+D, LC-DFT…

  15. DMC potential energy curves Adenine Adenine Thymine Thymine d: distance between 2 layers d: distance between 2 layers Adenine Adenine Thymine Thymine

  16. Potential energy curves Adenine Adenine Thymine Thymine d: distance between 2 layers d: distance between 2 layers Adenine Adenine Thymine Thymine

  17. DFT and QC potential energy curves Adenine Adenine Thymine Thymine d: distance between 2 layers d: distance between 2 layers Adenine Adenine Thymine Thymine

  18. DMC potential energy curves Adenine Adenine Thymine Thymine d: distance between 2 layers d: distance between 2 layers Adenine Adenine Thymine Thymine

  19. Basis set dependence Adenine Adenine Thymine Thymine d: distance between 2 layers d: distance between 2 layers Adenine Adenine Thymine Thymine

  20. Stacking energy at a reference geometry reference energy = CCSD(T)/CBS (Sponer et. al., Chem. Euro. J. 12, 2854, (2006)

  21. Summary • We applied DMC and DFT to the polymorphism of DIB. • The DMC result is consistent with the experiment. • The DFT results are not conclusive: The PBE, PW91, and B3LYP fail to give the correct sign of ΔE, while the LDA and B3LYP+D do, but strongly overestimate |ΔE|. • The base-pair-step stacking interaction in B-DNA (AT-AT stacking) were studied using DMC, HF, MP2, and DFT. • The DMC result is in good agreement with the reference result. • Commonly-used DFT methods fail to describe the stacking, but some new functionals show good performances.

  22. Thank you !

  23. a variational problem of minimizing the energy with respect to the many-electron wavefunction Quantum chemistry; Molecular orbital (MO) methods (HF, CI, CC, etc.) Quantum Monte Carlo (QMC)methods Electronic structure calculations The Schrodinger equation for many-electron systems Ab-initio calculations Density functional theory (DFT) LDA, PW91, BLYP, B3LYP, etc Solid state physics; a variational problem of minimizing the energy with respect to the electron density

  24. physical quantity : an expectation value of an operator Monte Carlo sampling according to 3N dimensional distribution Variational Monte Carlo (VMC); ; a choice of trial wavefunction Accuracy of results strongly depends on the trial wavefuntion adopted. Diffusion Monte Carlo (DMC); ; an imaginary time-dependent trial wavefunction that converges to the exact wavefunction after a long enough interval  Iterative operation of the imaginary-time evolution operator Purification of a trial wavefunction Quantum Monte Carlo methods

  25. Present study • DIB is a challenging benchmark system for ab initio methods. • Reliable reference calculations are necessary • But for periodic systems.... highly accurate quantum chemistry methods are not applicable due to their computational cost. • Diffusion Monte Carlo (DMC) is one of the most promising candidates for treating such systems from the viewpoint of accuracy and computational cost. • Purpose of our study • To investigate the relative stability of two polymorphs of DIB using DMC and several DFT methods. • LDA(SVWN), PW91, PBE, B3LYP • B3LYP+D: Semiempirical dispersion correction to DFT

  26. Stacking energy at a reference geometry Stacking energy (kcal/mol)

  27. The present method -1.17447(4) DMC Example: H2 molecule < Ground state Energy of H2 > (6-311++G** were used in all the calculations) -1.1330 HF Conventional methods -1.17203(1) VMC -1.1723 CISD (=Full CI) -1.17447 Exact J. Chem. Phys. 49, 404 (1965). Trial wfn. = (Ground state) (excited state) Signal(=target) Noise Trial wfn. = (Ground state) (excited state)

  28. Current topics in electronic structure calculations A variety of materials: e.g. molecular crystals, etc. Conventional exchange-correlation potentials (LDA, GGA, Hybrid…) sometimes give different predictions… universally reliable reference calculations Computational aspects: SOLIDS! Accurate quantum chemical calculations can NOT be applied to Although Full CI and QMC give results of almost the same accuracy, Full CI scales as O(N!), while QMC scales as O(N3 )-O(N4) QMC algorithm is intrinsically parallel Suitable for recent massive parallel computers

  29. Trial wfn. = (Ground state) (excited state) Signal(=target) Noise Total Energy (a.u.) Trial wfn. = (Ground state) (excited state) VMC energy “Purification” of the trial wavefuntion Example taken from NiO calculation. DMC energy Imaginary time steps The basic idea of DMC

  30. Computational details: QMC • DMC calc. using QMCPACK • Pseudo potential = Trail-Needs type • Taken from CASINO pseudo potential library • The total number of electrons = 168 • Trial wavefunction = Slater-Jastrow type • Slater determinant: LDA (cut off energy = 70 Ry) • Jastrow factor: one- and two-body terms • 6 parameters for the one-body and 4 parameters for the two-body • Optimized by minimizing the variance of local energy in VMC • computational conditions of DMC • # of walkers = 16384, # of MC steps = 3.2×107

  31. Computational details: DFT • All-electron DFT calculations using the Crystal 09 code • The total number of (core) electrons = 584 (416) • Basis set: • 6-31 G for the hydrogen and carbon atoms • 3-21 G for the iodine atom • K-point sampling = 1x3x4 (1-2 meV accuracy) • Exchange-correlation functionals • LDA(SVWN), PW91, PBE, B3LYP • B3LYP+D: Semiempirical dispersion correction to DFT Grimme Scheme: J. Comput. Chem. 27, 1787 (2006). • Vibration frequencies within the harmonic approximation J. Comput. Chem. 25, 888, 1873 (2004). Vibration frequencies calculated within the harmonic approximation using the Crystal 09 code Refs.: J. Comput. Chem. 25, 888 (2004); 25, 1873 (2004)

  32. DFT geometry optimizations:Lattice constants and volumes Length and volume are in units of Å and Å3, respectively. GGAs incorrectly predict that the α phase has a larger volume than the β one. LDA and B3LYP+D give the correct ordering and the most compact structures. B3LYP strongly overestimates the cell volume for both phases. The deviation from experiment is more than 10% (not satisfactory).

  33. Relative energies using the optimized geometries LDA B3LYP+D B3LYP PBE PW91 -12 (25) -16 (68) -128 (-123) -163 (-155) -165 (29) Numbers in the parenthesis are ΔE using the experimental geometries All five functionals predict the α phase to be more stable than the β phase at the optimized geometries.

  34. Importance of vibration analysis in polymorphism of molecular crystals The zero-point energy (ZPE) is of the same order of magnitude as the lattice energy difference between the polymorphs. In some cases, the ZPE can influence the relative stability. glycine case: S. A. Rvera, et al. Cryst. Growth Des. 8, 3905 (2008). ΔHexp= 0.27 Lattice energy difference = 1.28 kJ/mol Vibrational ΔH = 1.01 kJ/mol Figure: Diagram showing that a correct computation of the lattice energy difference for α- and γ-glycine must have γ more stable than α-glycine by ca. 1.28 kJ/mol. α-glycine γ-glycine

  35. Phonon calculations at the Γ point:Zero-point energy All energies are in units of meV. EL: Electronic energy; E0: Zero-point energy; ΔE := E (α) – E(β) Vibration frequencies calculated within the harmonic approximation using the Crystal 09 code Refs.: J. Comput. Chem. 25, 888 (2004); 25, 1873 (2004)

  36. Phonon calculations at the Γ point:Temperature dependence P = 0.101325 MPa (Experimental transition temperature: Tpt = 326 K) Tpt=471 K Tpt=78 K [meV] = G(α)-G(β) Tpt=1005 K [Kelvin] G = U + PV – TS; U = EL + E0 + ET G: Gibbs free energy; P: Pressure; V: Volume; T: Temperature; S: Entropy; U: Internal energy; EL: Electronic energy; E0: Zero-point energy; ET: Thermal contribution to the vibration energy

  37. Finite-size (FS) effects in QMC Phys. Rev. B 53, 1814 (1996). Phys. Rev. E 64, 016702 (2001). Phys. Rev. Lett. 100, 126404 (2008). All energies are in units of meV. The present DMC FS correction = 50 meV • QMC calculations for periodic systems: • Super cell calculations and their extrapolation to infinite size • Twist-averaged boundary conditions (for metallic systems) • Kwee, Zhang, and Krakauer scheme • QMC FS corrections are estimated from modified LDA calculations with finite-size functionals.

  38. Basis set dependence in DFT energiesin progress… • Commonly used basis sets in CRYSTAL09: Standard Gaussian (Pople) basis sets • STO-nG n=2-6 (H-Xe), 3-21G (H-Xe), 6-21G (H-Ar) • plus polarization and diffuse function extensions • Standard basis sets are optimized for molecular systems. • Are they appropriate for molecular crystal systems? • A number of ab initio studies on molecular crystals have employed standard basis sets. eg. α-quartz: 6-21G* and 6-31G* for Si and O, respectively [J. Comput. Chem. 25, 888 (2004)] • No investigation of the basis set dependence because of its computational cost. • The quality of basis sets significantly affects the final results in some cases. • It is crucial to investigate the basis set dependence.

  39. Basis sets to be considered: • Standard Gaussian basis sets • 6-31G & 3-21G (# of AOs = 664) • 6-311G (# of AOs = 816) • MDT: (# of AOs = 784) • from Dr. Mike D. Towler’s web site: http://www.tcm.phy.cam.ac.uk/~mdt26/crystal.html • H: (7s1p)/[3s1p], C: (9s3p1d)/[3s2p1d], I: (29s20p10d)/[8s7p3d] (# of AOs = 784) • Those exponents were optimized with respect to the HF total energy of a typical solid. • MDT+d: (# of AOs = 824) • add an extra diffuse function to the I atom • The exponent in the d function is optimized with respect to the DFT total energy. • opt-MDT: • optimize the exponent in the outermost function centered on each atom, i.e., a p function with the smallest exponent for H, d functions with the smallest exponents for C and I. • opt-MDT+d: • optimize exponents in the extra d functions and the outermost functions only.

  40. Energetics at the experimental geometry:The SVWN total energy of the α phase 243.92 8.33 0.45 0.31 0.31 0.28 The MDT group gives a better energy than standard Gaussian basis sets. Optimization of exponents with respect to a solid total energy is very important; Standard sets are optimal only for molecules, NOT SOLIDS.

  41. Energetics at the experimental geometry:The B3LYP total energy of the α phase 243.99 NA 6-311G 0.36 0.24 0.23 0.21 Results similar to SVWN were obtained except for 6-311G, i.e., the B3LYP energy with 6-311G did not converge. The exponent optimization also plays an important role in the SCF convergence.

  42. Energetics at the experimental geometries:relative stability ΔE = E(α) – E(β) All energies are in units of meV. • What is a criterion for the “optimal” basis set? • NOT CLEAR TO ME…. The energy difference depends on the basis set. The MDT group lowers the energy difference, compared to 6-311G. An extra d function on the I atom lowers the energy difference. For some cases in PW91 and B3LYP, the sign of the energy difference turns to be negative.

  43. Summary and future work • We investigated the relative stability of two polymorphs of DIB using the DMC and DFT methods. • The DMC result correctly predicts that the α phase is more stable than the β phase at zero temperature. • The DFT results were inconsistent and inconclusive. • A proper treatment of electron correlation is important. • More systematic study of basis set dependence • Pseudopotential DFT calculations • Phonon dispersion calculations

  44. Polymorphism in molecular crystals • Polymorphism: the existence of more than one form of a compound • Each polymorph has different physical and chemical properties • Important in, e.g., pharmaceutical industry and so on. • Computational studies of polymorphism in molecular crystals • mostly relies on molecular mechanics (MM) because of computational costs • ab-initio studies are also useful for verifying the reliability of MM Phonon calculations, etc. with empirical potentials e.g., S. L. Price, Adv. Drug Delivery Rev. 56, 301 (2004), For example: (JACS 127, 3038 (2005) ) - at zero temperature: ground state energy calculations by DFT - at finite temperature: Car-Parrinello MD simulation based on DFT

  45. Problems in ab-initio studies of Polymorphism in molecular crystals • Computational costs: • Unit cell contains a number of complicated molecules • Periodicity of molecular crystal structures • Accuracy of theoretical methods: • The energy difference between polymorphs is quite small (- order of 1 mH) • Noncovalent bond: hydrogen bond/ dispersion-dominated complexes Accurate quantum chemistry methods such as CI and CC are NOT applicable to periodic systems… Mostly, DFT Electron correlation effects play an important role in describing such chemical bondings! Many DFT exchange-correlation functionals fail… DMC is a good benchmark for molecular crystals

More Related