1 / 91

Major methods

Major methods. The J matrix engine (Ahmadi & Almlöf, 1995; White & Head-Gordon, 1996; Shao &, Head-Gordon, 2000) Density fitting (RI) (Dunlap & Rösch, Ahlrichs) Numerical solution of the Poisson equation (Delley, Becke, Termath & Handy, Baker, Nemeth & Pulay, unpublished)

xanti
Download Presentation

Major methods

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. Major methods • The J matrix engine (Ahmadi & Almlöf, 1995; White & Head-Gordon, 1996; Shao &, Head-Gordon, 2000) • Density fitting (RI) (Dunlap & Rösch, Ahlrichs) • Numerical solution of the Poisson equation (Delley, Becke, Termath & Handy, Baker, Nemeth & Pulay, unpublished) • Gaussian fast multipoles and related methods (White, Johnson, Gill, Head-Gordon, 1994; Strain, Scuseria and Frisch, 1996; Challacombe & Schwegler, 1997) • Pseudospectral methods (Friesner) • Intermediate Plane Wave expansion (GAPW and the Fourier Transform Coulomb method) Selected Chapters Budapest Fall 2011

  2. Plane wave basis: widely used in solid-state physics • The Coulomb operator is diagonal in momentum space: |r1-r2|-1 = (1/22)dk k-2 exp[ik( r1 - r2)] • Include the molecule in a sufficiently large box (for simplicity, it is treated as a cube with sides L) • Plane wave basis exp(iagr); g integer |g|< gmax; a=2/L • Introduce a real-space grid with spacing h = L/2gmax • Functions which can be represented by the plane wave basis are exactly determined by the real-space grid values • Calculate each quantity in the appropriate (momentum or real space) representation, using fast Fourier transformation to switch between representations Selected Chapters Budapest Fall 2011

  3. Plane wave vs. Gaussians basis sets • Plane wave basis sets have important advantages: flexibility, efficient calculation of forces (Hellmann-Feynman), extensible to crystals, no BSSE, no near-singularities due to overcompleteness • They also have problems: periodic ghost images, divergences for charged systems, cannot describe compact orbitals, very large basis sets are neededpoor scaling, exact exchange (for hybrid DFT) is extremely costly Selected Chapters Budapest Fall 2011

  4. The Fourier Transform Coulomb method: The best of both worlds Selected Chapters Budapest Fall 2011

  5. Calculating the Coulomb operator in a Gaussian basis by plane-wave expansion ADVANTAGES vs. PURE PLANE WAVES • No new model chemistry: strictly comparable to existing Gaussian-based codes • Conventional (pseudo)diagonalization is applicable • Can be combined with the conventional calculation of the exchange operator - hybrid DFT is possible DISADVANTAGES • The extra flexibility of the plane wave basis is lost • Calculation of the forces is more expensive • BSSE reappears • For large diffuse basis sets, Gaussians are ill-conditioned Selected Chapters Budapest Fall 2011

  6. The GAPW method of Parrinello et al. The application of plane wave expansion for the solution of the Coulomb problem in Gaussian basis was pioneered by the Parrinello group: • G. Lippert, J. Hutter and M. Parinello, Theor. Chem. Acc. 103 (1999) 124. (Pseudopotentials; Blöchl’s PAW) • M. Krack and M. Parrinello, Phys. Chem. Chem. Phys. 2 (2000) 2105. (All electron: for the first timecomparison of absolute energies with Gaussian results) • Our motivations and goals are similar and both are based on plane waves. However, the methods are very different technically. In particular, we aim at much higher precision Selected Chapters Budapest Fall 2011

  7. Calculation of the Coulomb operator in plane wave basis Let us assume first that the whole density can be represented by plane waves • Calculate the density (r) on the real space rectangular grid (asymptotically O(N)) • Fast Fourier Transform to k space: (k) O(NlogN) • Calculate the potential V(k)=C(k)/k2 O(N) • Reverse FFT yields V(r): the Coulomb potential on the grid • Add VXC(r) (currently done separately) • Calculate the Coulomb matrix elements J by quadrature (exact if the basis functions can be represented on the PW basis), asymptotically O(N) Selected Chapters Budapest Fall 2011

  8. Difficulties I: Ghost Images Selected Chapters Budapest Fall 2011

  9. a+D a molecule ghost Eliminating ghost images • Periodic ghost images  use a truncated Coulomb operator • f(r) = 1/|r| if r<D, = 0 if r>D (D=diameter of the molecule) • analytical Fourier transform (4/k2)(1-cos(kD)) • contrary to accepted wisdom, its cost is negligibl • Can be generalized to 1, 2, 3-dimensional periodic systems L extended  L+D  2L ghost... Selected Chapters Budapest Fall 2011

  10. Generalization to systems periodic in 1, 2 or 3 dimensions • This concept can be generalized to polymers, layers or crystals by using an appropriately modified form of the Coulomb operator. • For 1-D: f(r)=1/r if x2+y2<(D/2)2, 0 otherwise • For 2-D: f(r)=1/r if |z|<D/2, 0 otherwise • For 3-D” f(r)=1/r • Exact Fourier transforms: • {1+D[kJ1(Dk )K0(Dk║) - k║J0(Dk )K1(Dk║)]}4k-2 • {1-sign(k )exp(-kD/2)} 4k-2 • 4k-2 (J0, J1: Bessel functions of the first kind; K0,K1: of the 2nd kind) Selected Chapters Budapest Fall 2011

  11. Treatment of compact orbitals They are necessary also for the inner valence shells (O,F…) At first we have retained the usual Gaussian Electron Repulsion Integral (ERI) evaluation for integrals which cannot be calculated in plane wave basis Some integrals involving compact charge distributions can be calculated by a plane-wave expansion. E.g. contributions arising from (c|d) [c=compact, d=diffuse charge distribution (AO product)] can be evaluated because the Coulomb operator is diagonal in k-space. However, the subsequent numerical quadrature requires that |c> is subjected to a low-pass filter operation. Currently: multipole approximation for (c|c’) contributions Selected Chapters Budapest Fall 2011

  12. Characteristics of the FTC method • Results are essentially exact with a sufficient PW basis • Asymptotic scaling with mol. size (at constant basis quality) is Nlog N • Scaling with the basis set size at constant mol. size is O(N2) • The onset point is favorable for moderate (6-31G*) and larger basis sets; already for a small molecule like aspirin and the 6-31G** basis, the calculation is faster than the conventional one • Without the multipoles no savings for the 3-21G basis in systems with ~100 atoms • Exact exchange is calculated by Electron Repulsion Integrals Selected Chapters Budapest Fall 2011

  13. Accuracy comparison with GAPW (error relative to accurate Gaussian calculations for small molecules) Selected Chapters Budapest Fall 2011

  14. Examples Selected Chapters Budapest Fall 2011

  15. AspirinC9H8O4 Selected Chapters Budapest Fall 2011

  16. Timings for aspirin (B-LYP)a a In minutes on a 2.4 GHz Pentium 4 *Slightly decontracted Selected Chapters Budapest Fall 2011

  17. Sucrose, C12H22O11 Selected Chapters Budapest Fall 2011

  18. Timings for sucrose (B-LYP)a a In minutes on a 2.4 GHz Pentium 4 Selected Chapters Budapest Fall 2011

  19. Taxol C47H51NO14 Selected Chapters Budapest Fall 2011

  20. Timings for taxol (B-LYP)a a In minutes on a 2.4 GHz Pentium 4 Selected Chapters Budapest Fall 2011

  21. Scaling with respect to molecular size(alanine)n, n=1-15, 6-311G(2df,2pd) 627 min 30.2 min 5.8 min Selected Chapters Budapest Fall 2011

  22. Scaling with respect to molecular size - (alanine)n, n=2-15 log-log plot Selected Chapters Budapest Fall 2011

  23. Scaling with respect to basis set size - (alanine)5, 10 Selected Chapters Budapest Fall 2011

  24. Average computational costs with and without multipole expansion in FTC Selected Chapters Budapest Fall 2011

  25. Conclusions • The FTC method, i.e. the calculation of the Coulomb operator by plane-wave expansion of the valence AOs has significant advantages, particularly for large basis sets: • Accurate • Low scaling with molecular size • Low scaling with basis set size • Early onset • Main projected applications: DFT calculations • on large clusters • of NMR chemical shifts • direct molecular dynamics Selected Chapters Budapest Fall 2011

  26. A new project: Genuine Mixed Gaussian and Plane Wave Basis Sets Both the Basis Set Overcompleteness (BSO) and the Basis Set Superposition Error (BSSE) are caused by diffuse functions Diffuse functions are essential for intermolecular (van der Waals) forces, for negative ions and negatively charged fragments, for excited states, polarizabilities, hyperpolarizabilities… Replace the diffuse basis functions (and only these) by a plane wave basis is a box The number of plane waves remains manageable. The truncated Coulomb operator trick eliminates periodic ghost images. However: new kinds of integrals are needed Selected Chapters Budapest Fall 2011

  27. Problems with atomic basis sets 1 • Small and medium-sized atomic basis sets are probably the best choice for routine molecular calculations. Large basis sets, particularly if augmented with diffuse functions, have serious problems in larger 2-D and 3-D systems • In principle, AO expansions around one center are complete if all states (including continuum states) are included. Therefore, a full AO expansion around all centers is overcomplete: it contains redundant functions • For finite basis sets, the basis becomes nearly linearly dependent. There is a linear combination of basis functions with normalized coefficients (i |cik|2 = 1) which has a very low norm in physical space (ijcik*Sijcjk= k << 1); Sij =<i |j > This leads to large MO coefficients (up to 102 – 103). Selected Chapters Budapest Fall 2011

  28. Problems with atomic basis sets 2 • Explicitly (in wavefunction-based correlation theories), or implicitly (in SCF), we always transform to MOs: this transformation becomes numerically unstable in the presence of large, cancelling MO coefficients. This is always the case for diffuse basis functions and larger systems which became accessible recently. Particularly bad for wavefunction-based correlation theories (virtuals!). • Possible remedies: • Eliminate linear combinations of basis functions which belong to low eigenvalues of the overlap (Gram) matrix. This may cause discontinuities on potential energy surfaces and does not help in really serious cases. Penalty function? • Use an orthogonal basis set Selected Chapters Budapest Fall 2011

  29. An example Consider a fragment of a crystal lattice with 1 Å lattice constant and diffuse Gaussians (=0.036). Overlap of neighbors S=0.93774 Number of atoms (k3) 8 27 64 Lowest eigenvalue of S 2.41E-4 1.45E-7 2.02E-10 min(1-S)2k asymptotically Even in a 1-D ring, if only nearest neighbors overlap, the overlap matrix is singular for even n if S0.5. Most atomic basis sets, and all diffuse basis sets become singular in the limit, and have very low eigenvalues for larger systems. Overlap of two 1s Slater functions (=1.24) is 0.66 at R=0.74 Å Selected Chapters Budapest Fall 2011

  30. Problems with atomic basis sets 3 Another ubiquitous problem with atomic basis sets is the basis set superposition error (BSSE). Particularly significant for weak intermolecular interactions but also present within the molecule and is nearly impossible to eliminate. BSSE: Artificial attraction between molecular fragments due to the improved description of the orbitals as two fragments approach each other BSSE is usually corrected for by the Boys-Bernardi Counterpoise correction (CP). There are strong theoretical arguments for the validity of the original CP methodology (J. van Lenthe, van Duijneveldt). However, in practice CP overcorrects. Selected Chapters Budapest Fall 2011

  31. The 3a1 MO in water along a line passing through the symmetry axis. Red= Gaussian basis, Green=small PW Selected Chapters Budapest Fall 2011

  32. How many plane waves are needed? • To replace diffuse AOs, a plane wave density of about 0.4-1 PW/Å is needed. Not economical for small systems because even a hydrogen atom needs a sizable box. However, the box size does not grow proportionately with the molecular size • We don’t have to be as accurate as in FTC where the plane waves are used to evaluate integrals over Gaussian orbitals. Approximating orbitals means changing the Hamiltonian  the error is first-order. Changing a (nearly optimum) basis function introduces only a second-order error in the energy. Selected Chapters Budapest Fall 2011

  33. A diffuse Gaussian and its PW simulation Gaussian with exponent 0.036 (recommended for H diffuse S orbital by Radom, P. Schleyer) Largest PW density 0.2 PW/a0 distance/a0 Selected Chapters Budapest Fall 2011

  34. Using PWs as basis functions is much less demanding than reproducing a Gaussian Largest PW density 0.167 PW/a0 (0.3/Å Largest PW density 0.133 PW/a0 Selected Chapters Budapest Fall 2011

  35. Previous work: Scattering Integrals for electron-molecule scattering cross sections • McWeeny, Acta Cryst. 1953; Miller and Krauss, JCP 1967, Stewart, JCP 1969; Liu, JCP, 1973; Ostlund CPL 1975; Pulay et al. JCP 1983; Kohl and Pulay, J. Mol. Struct. 1984 (elastic) • Polašek and Čarsky Int. J. Quantum Chem. 2007, J. Phys. B. 2010 (Accuracy is OK for scattering, not sufficient for electronic structure; Kuchitsu, Okuda, Tachikawa Int. J. Quantum Chem. 2009 (exchange included) • More general integrals: Čarsky and Polašek, J. Comput. Phys. 1998. The accuracy of the routine used is sufficient for scattering but unfortunately not sufficient for molecular calculations Selected Chapters Budapest Fall 2011

  36. PreviousWork - GIAOs Gauge-Including Atomic Orbitals (GIAOs)Polynomial(x,y,z) exp[-a(r – A)2] exp[(-i/2c) BAr] • Usually, only the derivatives of the GIAOs with respect to the magnetic induction B are needed. In this case, no complex arithmetic is necessary. Ditchfield, Mol. Phys.1974 • In our case, the full complex integrals are needed. GIAO Integrals: Ishida JCP 2003. Useful for very strong magnetic fields (cannot be created in the lab) • Also H. Fukui, 1998, 1999 • Genuine mixed basis proposal: Colle, Fortunelli, Simonucci, Nuovo Cimento 1987, 1988 (no continuation) • Problem: probably divergent terms in the expansion Selected Chapters Budapest Fall 2011

  37. Integral types The product of a Gaussian and a plane wave gives a Gaussian with complex center coordinates (G. G. Hall; also Gauge-Including Atomic Orbitals needed for magnetic properties) Integral types (g = Gaussian, w = plane wave): (gg|gg) : use traditional Gaussian techniques (Rys polynomial, Obara-Saiko, McMurchie-Davidson etc.) for complex coord (gw|g’g”), (gw|g’w’): traditional integral methods extended to complex origins/GIAOs (we use the Rys polynomial technique of H. King and M. Dupuis). Need complex Fm(t). (gg’|w’w”),(gw|w’w”): expand |g> in a plane wave basis. Compact functions are OK: Only components of |g> with |k|<|k’- k”| areneeded. Use the truncated Coulomb operator. Same for <ww’|w”w’”>: only diagonal terms Selected Chapters Budapest Fall 2011

  38. Integral evaluation: Rys quadrature H. F. King, M. Dupuis, J. Rys, 1976- The Gaussian transform of the inverse distance factorizes The product of 2 Gaussians or a PW and a Gaussian is a Gaussian (possibly with complex center coordinates) which factorizes to x,y,z factors: P, Q are Gaussians (exponents p,q) with polynomial factors in x,y,z; Ix, Iy, Iz are 2-D integrals (simple but numerous) Change of variable: t2 = u2/(+u2),  = pq/(p+q) Selected Chapters Budapest Fall 2011

  39. Motivation • Needed canonical to check approximate methods (local MP2, FMO=Fragment Molecular Orbital, Many-Body Expansion, Density Fitting) against exact results • Commonly used programs failed or were too expensive • Local electron correlation becomes less efficient for systems with aromatic delocalization (porphyrins, fullerenes, graphenes). • There are excellent local correlation implementations (the Saebo-Pulay technique is implemented in MOLPRO and JAGUAR; TRIM in QChem, Scuseria has a method) • An efficient integral transformation is necessary for higher configuration-based correlation methods Selected Chapters Budapest Fall 2011

  40. Four-index transformation: theory • Essentially all effort in canonical MP2 goes into the four-index integral transformation (i,j occupied, a,b virtual) (ai|bj)= pqrs (pq|rs)CpaCqiCrbCsj , usually broken up to four quarter transformations • Formal scaling of the first quarter transformation,(pi|rs) =  q (pq|rs)Cqi, ( i=1,…n)is O(nN4), where n is the number of correlated occupied orbitals, and N is the number of AOs • Subsequents steps scale as O(n2N3). As N>>n (for better basis sets N 6-10n), the first transformation is expected to dominate • Traditional: Saunders and Van Lenthe, Werner and Meyer Selected Chapters Budapest Fall 2011

  41. Head-Gordon, Pople, Frisch, Chem. Phys. Lett. 153 (1988) 503 Cubic memory demand. Calculate in batches (Gets expensive if there is insufficient memory) For a comparison of methods, and some new algorithms (two methods that require O(N) memory) see: M. Schütz, R. Lindh and H.-J. Werner, Mol. Phys. 96 (1999) 719 Saebo and Almlöf, Chem. Phys. Lett. 154(1989) 83 Quadratic memory; only one permutation used (integrals calculated 4x) do p do r calculate all (pq|rs)=Xqs Y=CTXC (matrix mult.) store Yij=(pi|rj) on disk end do r end do p Integral transformation for large calculations: memory limitations Selected Chapters Budapest Fall 2011

  42. Prescreening • Efficient prescreening in the integral evaluation phase based on local correlation ideas (Chem. Phys. Lett. 2001, 344, 543) • Basic idea: an integral which is negligible in local MP2 is also negligible in the canonical method • Pair correlation amplitudes between well-separated electron pairs decreases as R-3 • In large, well-localized molecules, only a small fraction of the AO integrals needs to be calculated and transformed • Dilemma: dense matrix multiplication is very fast but scales steeply; sparse matrix multiplication has good scaling but is slow. Solution: compact the matrices before transformation Selected Chapters Budapest Fall 2011

  43. Effect of the threshold on the accuracy of the calculated MP2 energy Hexapep=N-formyl pentalanine amide For well-conditioned basis sets, the default threshold guarantees accuracy to a few Eh. Selected Chapters Budapest Fall 2011

  44. Timings for calix[4]arenea,b, C60a and C74(C1)c a In minutes on an 800 MHz Pentium III bTetramethoxy-calix[4]arene c In minutes on a 3 GHz Pentium 4 Xeon Selected Chapters Budapest Fall 2011

  45. Scaling • The ultimate scaling is determined by the second half transformation which is performed just like in traditional MP2 • Routine calculations on a single processor are possible for molecules with ~1000 basis functions and ~300 correlated electrons in the absence of symmetry • Larger calculations are possible for symmetrical molecules Selected Chapters Budapest Fall 2011

  46. Parallelization • First half transformation: The Saebo-Almlöf algorithm naturally parallelizes over two fixed AO indices p and r • Yoshimine bin sort. Each node owns a subset of the pair indices (i,j). The parallel sort sends the half-transformed integrals to the appropriate node. Synchronization delays are avoided by spawning a set of independent bin receive/write daemon processes • The second half transformation naturally parallelizes over the orbital pair indices (i,j) Selected Chapters Budapest Fall 2011

  47. Parallel Yoshimine bin sort Master: spawn slaves and bin writer daemons Slave 1 Read half-transformed integrals from disk and sort them in bins by i,j. Send full bins to the appropriate nodered or blue Slave 2 Read half-transformed integrals from disk and sort them in bins by i,j. Send full bins to the appropriate nodered or blue Bin writer 1 (daemon) Receive sorted bins and store them on disk Bin writer 2 (daemon) Receive sorted bins and store them on disk spawn send Selected Chapters Budapest Fall 2011

  48. Parallel scaling: single CPUs Calix[4]arene: cc-pVTZ basis, 1528 BF, C2h symmetry Elapsed time on 24 1 GHz Pentium IIIs: 150.4 min for MP2 Selected Chapters Budapest Fall 2011

  49. Application to  stacking energies Selected Chapters Budapest Fall 2011

  50. General • Large planar aromatic systems (graphene sheets, porphyrins, phthalocyanines, DNA bases) are attracted by a considerable dispersion force. • Dispersion forces are also important in bucky onions and carbon nanotubes • The true dimerization energy is difficult to measure (heat of vaporization would give a clue) • In the limit of large parallel sheets, the dispersion force diminishes as 1/d4, not as 1/d6 Selected Chapters Budapest Fall 2011

More Related