Major methods
This presentation is the property of its rightful owner.
Sponsored Links
1 / 91

Major methods PowerPoint PPT Presentation


  • 46 Views
  • Uploaded on
  • Presentation posted in: General

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)

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.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Presentation Transcript


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)

  • 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


Plane wave basis widely used in solid state physics

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


Plane wave vs gaussians basis sets

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


The fourier transform coulomb method the best of both worlds

The Fourier Transform Coulomb method: The best of both worlds

Selected Chapters Budapest Fall 2011


Calculating the coulomb operator in a gaussian basis by plane wave expansion

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


The gapw method of parrinello et al

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


Calculation of the coulomb operator in plane wave basis

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


Difficulties i ghost images

Difficulties I: Ghost Images

Selected Chapters Budapest Fall 2011


Eliminating ghost images

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


Generalization to systems periodic in 1 2 or 3 dimensions

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


Treatment of compact orbitals

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


Characteristics of the ftc method

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


Accuracy comparison with gapw error relative to accurate gaussian calculations for small molecules

Accuracy comparison with GAPW (error relative to accurate Gaussian calculations for small molecules)

Selected Chapters Budapest Fall 2011


Examples

Examples

Selected Chapters Budapest Fall 2011


Aspirin c 9 h 8 o 4

AspirinC9H8O4

Selected Chapters Budapest Fall 2011


Timings for aspirin b lyp a

Timings for aspirin (B-LYP)a

a In minutes on a 2.4 GHz Pentium 4

*Slightly decontracted

Selected Chapters Budapest Fall 2011


Sucrose c 12 h 22 o 11

Sucrose, C12H22O11

Selected Chapters Budapest Fall 2011


Timings for sucrose b lyp a

Timings for sucrose (B-LYP)a

a In minutes on a 2.4 GHz Pentium 4

Selected Chapters Budapest Fall 2011


Taxol c 47 h 51 no 14

Taxol C47H51NO14

Selected Chapters Budapest Fall 2011


Timings for taxol b lyp a

Timings for taxol (B-LYP)a

a In minutes on a 2.4 GHz Pentium 4

Selected Chapters Budapest Fall 2011


Scaling with respect to molecular size alanine n n 1 15 6 311g 2df 2pd

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


Scaling with respect to molecular size alanine n n 2 15 log log plot

Scaling with respect to molecular size - (alanine)n, n=2-15 log-log plot

Selected Chapters Budapest Fall 2011


Scaling with respect to basis set size alanine 5

Scaling with respect to basis set size - (alanine)5,

10

Selected Chapters Budapest Fall 2011


Average computational costs with and without multipole expansion in ftc

Average computational costs with and without multipole expansion in FTC

Selected Chapters Budapest Fall 2011


Conclusions

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


A new project genuine mixed gaussian and plane wave basis sets

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


Problems with atomic basis sets 1

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


Problems with atomic basis sets 2

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


An example

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-72.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


Problems with atomic basis sets 3

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


Major methods

The 3a1 MO in water along a line passing through the symmetry axis. Red= Gaussian basis, Green=small PW

Selected Chapters Budapest Fall 2011


How many plane waves are needed

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


A diffuse gaussian and its pw simulation

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


Using pws as basis functions is much less demanding than reproducing a gaussian

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


Previous work scattering

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


    Previouswork giaos

    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


    Integral types

    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


    Integral evaluation rys quadrature

    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


    Motivation

    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


    Four index transformation theory

    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


    Integral transformation for large calculations memory limitations

    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


    Prescreening

    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


    Effect of the threshold on the accuracy of the calculated mp2 energy

    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


    Timings for calix 4 arene a b c 60 a and c 74 c 1 c

    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


    Scaling

    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


    Parallelization

    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


    Parallel yoshimine bin sort

    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


    Parallel scaling single cpus

    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


    Application to stacking energies

    Application to  stacking energies

    Selected Chapters Budapest Fall 2011


    General

    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


    Coronene dimer parallel displaced

    Coronene dimer(parallel displaced)

    • Counterpoise correction is important

    • At the 6-311G* level (936 BF, C2h)

      • Energy minimum at 3.27 Å (graphite 3.35 Å room temp.)

      • Minimum energy = -0.0442 Eh = -27.7 kcal/mol

      • At 3.35 Å: SCF=+22.1, corr. -49.6, B3LYP=+15.7 kcal/mol

      • Elapsed time on 4 dual-processor computers: 2hrs 1 min

    • At the 6-311G(2d,f) level (1512 BF):

      • Energy at 3.35 Å (close to minimum) is -0.0556 Eh=-35 kcal

      • Elapsed time on 10 800 MHz P III processors: 7.1 hrs

    Selected Chapters Budapest Fall 2011


    Potential curve for c 24 h 12 2 a b

    Potential curve for (C24H12)2a,b

    aCounterpoise

    corrected

    bDistance of sheets in graphite = 3.35 Å

    Selected Chapters Budapest Fall 2011


    Porphine dimer

    Porphine dimer

    • Porphyrins show a tendency to dimerize in solution

    • Porphyrin dimer: C40H28N8, 228 correlated electrons

    • 6-311G* basis: 948 basis functions 73 min total on 8 proc.

    • Dimerization energy at 3.35 Å is comparable to coronene with the same basis, -0.043 Eh = -26.9 kcal, SCF=+20.5, correlation -47.4

    Selected Chapters Budapest Fall 2011


    Coronene dimer

    Coronene dimer

    Selected Chapters Budapest Fall 2011


    Circumcoronene dimer c 108 h 36

    Circumcoronene dimer (C108H36)

    • Two circumcoronene molecules separated by 3.35 Å

    • 6-31G* binding energy = 0.138 Eh = 86.3 kcal/mol; with soft d functions (P. Hobza) 0.157 Eh=98.5 kcal/mol

    • Per carbon atom = 0.80 kcal/mol (significantly more than in coronene, 0.57 kcal/mol (0.91 with soft d functions)

    Selected Chapters Budapest Fall 2011


    Circumcoronene dimer 6 31g 0 25d

    Circumcoronene dimer 6-31G(0.25d)

    Sandwich

    P.D. 1/4

    P.D. 3/4

    Parallel Displaced 1/2

    Selected Chapters Budapest Fall 2011


    Circumcoronene dimer the effect of horizontal displacement

    Circumcoronene dimer: the effect of horizontal displacement

    Vertical distance: 3.4 Å

    Selected Chapters Budapest Fall 2011


    Spin component scaled mp2

    Spin-component scaled MP2

    • MP2 overestimates -stacking energies, e.g. in benzene. The well depth is the result of a delicate balance between the Pauli repulsion and dispersion attraction. Errors in the latter are magnified.

    • S. Grimme (JCP, 2003). An empirical correction: increase slightly the opposite spin contribution, diminish strongly the parallel spin part. Note that these are NOT identical with singlet and triplet pairs.

    • Diminishes the absolute value of the dispersion interaction. Yields good values for the benzene dimer

    Selected Chapters Budapest Fall 2011


    Circumcoronene dimer mp2 versus scs mp2

    Circumcoronene dimer: MP2 versus SCS-MP2

    SCS-MP2 (CP corr.)

    MP2 (CP corrected)

    Selected Chapters Budapest Fall 2011


    Pyrazine in a cram carceplex courtesy dr suyong re

    Pyrazine in a Cram carceplex (courtesy Dr. Suyong Re)

    6-31G* basis, 1476 basis functions, D2 symm, 8 2.4 GHz Pentium Xeon processors.

    Elapsed times:

    SCF: 36.9 min

    First ½ tr. 51.0 min

    Sort and

    Second ½ tr.295.9 min

    Total job383.9 min

    C72H52N2O24

    Selected Chapters Budapest Fall 2011


    N methyl pyrrolidone in a cram carceplex courtesy dr suyong re

    N-methyl pyrrolidone in a Cram carceplex (courtesy Dr. Suyong Re)

    6-31G* basis, 1500 basis functions, no symm, 4 or 48 3 GHz processors.

    Elapsed times in MP2 (in min):

    # processors 4 48 ratio

    Total MP2 1923 159 12.1

    C73H57NO25

    Selected Chapters Budapest Fall 2011


    Can semi local dft describe dispersion

    Infinite potential barrier

    Electrons can’t penetrate

    Electromagnetic forces move

    freely through the barrier

    Can (semi)local DFT describe dispersion?

    A thought experiment shows that (semi)local DFT cannot describe genuine dispersion, even if there is an artificial attraction resembling dispersion

    e-

    Selected Chapters Budapest Fall 2011


    Can semi local dft describe dispersion1

    Infinite potential barrier

    Electrons can’t penetrate

    Electromagnetic forces move

    freely through the barrier

    Can (semi)local DFT describe dispersion?

    A thought experiment shows that (semi)local DFT cannot describe genuine dispersion, even if there is an artificial attraction resembling dispersion

    e-

    Selected Chapters Budapest Fall 2011


    Full accuracy local mp2 svein saebo

    Full-accuracy local MP2(Svein Saebo)

    Selected Chapters Budapest Fall 2011


    Scaling of configuration based correlation methods

    Scaling of configuration-based correlation methods

    • Traditional correlation theories scale steeply and this is not reduced automatically by sparsity

    • This steep scaling is an artifact of using delocalized canonical MOs, as dynamical correlation is very local

    • To get rid of it, correlation theories must be formulated in terms of localized quantities: either localized molecular orbitals or atomic orbitals (basis functions)

    Selected Chapters Budapest Fall 2011


    Advantages of localized orbitals

    +

    -

    +

    -

    Left-right correlation

    Up-down correlation

    Advantages of localized orbitals

    Two sources of savings

    • Distant localized orbitals interact weakly (~ R-6): neglect

    • Correlation orbitals must have nodes dissecting an occupied orbital. Only the basis functions localized in the neighborhood of the occupied orbital are effective: truncate the correlation space

    Selected Chapters Budapest Fall 2011


    New local correlation program

    New local correlation program

    • S. Saebo, P. Pulay, J. Chem. Phys. 2001, 115, 3975

    • Simultaneous transformation of 2 indices (P. R. Taylor, Int. J. Quantum Chem.1987, 31, 521). This is not the essence of the method but makes it very memory-efficient

    • Memory and disk space requirements are minimal (E.g. S. Sabo evaluated the MP2 energy of (glycine)50, 6-31G(d): 3118 BF on a single PC. Similarly, (glycine)30, 6-311G**: 2638 BF. However, distant pairs were approximated in these calculations

    • With distant pairs added, the results are essentially identical with canonical theory

    Selected Chapters Budapest Fall 2011


    Major methods

    SCF: N2.02

    MP2:N1.33

    Selected Chapters Budapest Fall 2011


    Problems with the current code

    Problems with the current code

    • Pair domains are large. If distant pairs are included, computational times approach the traditional timings for medium-sized molecules. The major computational task is MP2 iteration

    • Even weak correlation requires a considerable number of AOs for full description. The current program is wasteful because it uses the same local basis for all pairs: strong, weak, distant. A version which uses a smaller local basis for distant pairs would be much more efficient.

    • A more fundamental approach is to switch to a molecular orbital representation (W. Meyer)

    Selected Chapters Budapest Fall 2011


    Mo treatment of weak correlation

    MO treatment of weak correlation

    By changing to an MO basis, a much smaller basis is sufficient.

    Consider the pair correlation between two localized orbitals |i> and |j> (use canonical virtuals)

    Tijab = -(ia|jb)/(a + b -i - j)

    Invoke the multipole expansion but not for approximating integrals: 1/r12 = R-3[r1r2 -3R-2(r1R)(r2R) + …]

    Substitute into the pair wavefunction ij = a,bTijabijab

    If the energy denominator is assumed constant then, in the dipole approximation, three virtual orbitals on each center,

    Selected Chapters Budapest Fall 2011


    Major methods

    are sufficient to account for all correlation. ,=x,y,z; ris the position vector relative to the centroid of  i in the  direction, r=-Ri ; r is analogous for jOf course, the orbital energies are not constant. Numerical experiments by W. Meyer suggest that 2-3 denominator shifts are sufficient to generate essentially all the virtual space

    Selected Chapters Budapest Fall 2011


    Major methods

    Moreover, the dipole approximation holds only for very distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute.

    These considerations suggest that a small molecular orbital basis, 10-20 orbitals, determined individually for each localized internal orbital, is capable of describing essentially all dispersion interaction in the virtual space.

    The MOs belonging to different orbitals are not orthogonal. Note that this is essentially the pseudonatural orbital method of Meyer.

    It is best to introduce a pseudocanonical basis in the small subspace, separately for each localized occupied orbital, by diagonalizing the Fock matrix in the subspace

    Selected Chapters Budapest Fall 2011


    Natural orbitals of distant pairs

    Natural Orbitals of Distant Pairs

    A typical distribution of the natural orbital occupation numbers for a weak pair in hexane:

    Eigenvalues of S1/2T†STS1/2 for pair 13 4

    0.000000259 0.000000063 0.000000045 0.000000008 0.000000002

    0.000000001 0.000000001 0.000000001 0.000000001 0.000000000

    0.000000000 0.000000000 0.000000000 0.000000000 0.000000000

    0.000000000 0.000000000 0.000000000 0.000000000 0.000000000

    0.000000000 0.000000000 0.000000000 0.000000000 0.000000000

    Sum of all eigenvalues

    0.000000385

    Sum of the 3 largest eigenvalues

    0.000000368

    Sum of the 8 largest eigenvalues

    0.000000381

    Selected Chapters Budapest Fall 2011


    Major methods

    Selected Chapters Budapest Fall 2011


    Major methods

    Selected Chapters Budapest Fall 2011


    An efficient atomic orbital based mp2 gradient program

    An efficient atomic-orbital based MP2 gradient program

    Selected Chapters Budapest Fall 2011


    Efficient ao formulation of mp2 gradients j chem phys 120 11423 2004

    Efficient AO formulation of MP2 gradients J. Chem. Phys. 120, 11423 (2004).

    Based on the orbital-invariant form of the MP2 energy

    P. Pulay and S. Saebo, Orbital-invariant formulation and gradient evaluation in Møller-Plesset Perturbation Theory, Theor. Chim. Acta 69, 357 (1986).

    Also used to derived the equations for local MP2 gradient:

    A. El-Azhary, G. Rauhut, P. Pulay and H.-J. Werner, Analytical Energy Gradients for Local Second-Order Moller-Plesset Perturbation Theory, J. Chem. Phys., 108, 5185 (1998).

    Selected Chapters Budapest Fall 2011


    Generator state formulation

    Generator State Formulation

    Normalized to 4 (ab) or to 2 (a=b), and pairwise non-orthogonal,

    they allow a substantial simplification of the MPn, CI and Coupled Cluster equations, see

    P. Pulay, S. Saebo, and W. Meyer: An Efficient Reformulation of the Closed-Shell Self-Consistent Electron Pair Theory, J. Chem. Phys. 81, 1901 (1984).

    Selected Chapters Budapest Fall 2011


    The self consistent electron pair theory w meyer j chem phys 64 2901 1976

    The Self-Consistent Electron Pair TheoryW. Meyer, J. Chem. Phys. 64, 2901 (1976).

    An efficient matrix formulation of the singles and doubles CI problem and related methods. Both the amplitudes T

    and the integrals, e.g. the exchange and Coulomb integrals

    are collected in matrices. It is useful to introduce the contravariant amplitude matrices (Tji is the transpose of Tij)

    Selected Chapters Budapest Fall 2011


    Scep in ao basis

    SCEP in AO basis

    The canonical virtual orbitals are replaced by AOs:

    This complicates the formalism but is essential in local correlation, and is also advantageous in gradient evaluation. E.g. no problem with frozen cores. With a little extra work, “QM/MM” is possible: QM=MP2, “MM”=Hartree-Fock

    Selected Chapters Budapest Fall 2011


    Scep cont

    SCEP, cont.

    For canonical occupied orbitals this reduces to

    In a canonical MO basis for virtual orbitals (S=1; F=diag{}) it becomes the usual canonical formula

    Selected Chapters Budapest Fall 2011


    Orbital invariant form of mp2

    Orbital-invariant form of MP2

    It can be derived in several ways. For gradient evaluation the Hylleraas form is particularly useful:

    Ec = 21|H-E0|0 - 1|H0-E0|1 = minimum

    The MP2 energy can be expressed as Ec=ijeij

    (S is the overlap matrix, F is the Fock matrix).

    Differentiating the pair energy eij with respect to the (contra-variant) amplitudes gives the following equation

    Selected Chapters Budapest Fall 2011


    Timing benchmarks minutes on a 2 8 ghz pentium 4 xeon processor

    Timing benchmarks (minutes) on a 2.8 GHz Pentium 4 Xeon processor

    *No symmetry

    Selected Chapters Budapest Fall 2011


    Cavitand with 2 oxo 1 butanol

    Cavitand with 2-oxo-1-butanol*

    C32H32O10, 6-31G*, 652 BF

    A single 2.8 GHz P 4 Xeon

    processor (parallelization in progress).

    Elapsed times:

    SCF= 99.6 min

    MP2*= 309.0 min

    MP2 gradient= 883.5 min

    * incl. storing the amplitudes

    *Courtesy Dr. Re Suyong

    Selected Chapters Budapest Fall 2011


    Nucleic acid bases

    Nucleic acid bases

    • In connection with a recent plane-wave based DFT paper (M. Preuss et al., J. Comput. Chem. 2004, 25, 112) which predicts essentially planar geometries for these molecules, explicit correlation methods predict significant non-planarity for the amino (-NH2) group.

    • We have optimized the geometries of adenine, cytosine, guanine and thymine at the MP2 level using the cc-pVTZ, aug-cc-pvTZ and cc-pVQZ basis sets; the latter has the composition of 5s4p3d2f1g (755 basis functions for guanine). All three large basis sets give very similar geometries, showing that the results are converged. These are among the largest MP2 optimizations performed.

    Selected Chapters Budapest Fall 2011


    Guanine optimized

    Guanine, optimized

    Selected Chapters Budapest Fall 2011


    Summary

    Summary

    • Canonical MP2 is possible for large systems, by using natural localization

    • Parallel implementation can handle >2000 basis functions:  stacking

    • Full-accuracy local MP2 has been implemented

    • MP2 gradients up to ~atoms and ~1000 basis functions are feasible on a single PC

    • Array Files have been implemented; will be released soon.

    Selected Chapters Budapest Fall 2011


    The end

    The end

    Selected Chapters Budapest Fall 2011


    Saebo pulay local correlation

    Saebo-Pulay local correlation

    • Recent very efficient implementation by H.-J. Werner and his coworkers (M. Schütz, G. Rauhut) in MOLPRO

    • Generalization to Coupled Cluster methods (Werner group)

    • Other major local correlation methods (still in the formative stage):

      • Scuseria group (Laplace transform)

      • Head-Gordon group (diatomics in molecules, triatomics in molecules)

    • Still occasional difficulties with localization artifacts

    • We decided to develop a full accuracy MP2

    Selected Chapters Budapest Fall 2011


    Major methods

    Selected Chapters Budapest Fall 2011


    The end1

    The end

    Tinky

    Selected Chapters Budapest Fall 2011


  • Login