1 / 91

Major methods - PowerPoint PPT Presentation

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)

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.

Major methods

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

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

• 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 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

Selected Chapters Budapest Fall 2011

Calculating the Coulomb operator in a Gaussian basis by plane-wave expansion

• 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

• 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 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

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

Selected Chapters Budapest Fall 2011

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

• 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

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

• 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)

Selected Chapters Budapest Fall 2011

Examples

Selected Chapters Budapest Fall 2011

AspirinC9H8O4

Selected Chapters Budapest Fall 2011

Timings for aspirin (B-LYP)a

a In minutes on a 2.4 GHz Pentium 4

*Slightly decontracted

Selected Chapters Budapest Fall 2011

Sucrose, C12H22O11

Selected Chapters Budapest Fall 2011

Timings for sucrose (B-LYP)a

a In minutes on a 2.4 GHz Pentium 4

Selected Chapters Budapest Fall 2011

Taxol C47H51NO14

Selected Chapters Budapest Fall 2011

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)

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

Selected Chapters Budapest Fall 2011

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

Selected Chapters Budapest Fall 2011

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

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

• 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

• 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

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

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

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?

• 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

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

Largest PW density

0.167 PW/a0 (0.3/Å

Largest PW density

0.133 PW/a0

Selected Chapters Budapest Fall 2011

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

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

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

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

• 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

• 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

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

• 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

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]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

• 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

• 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

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

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

Selected Chapters Budapest Fall 2011

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)

• 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 (C24H12)2a,b

aCounterpoise

corrected

bDistance of sheets in graphite = 3.35 Å

Selected Chapters Budapest Fall 2011

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

Selected Chapters Budapest Fall 2011

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)

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

Vertical distance: 3.4 Å

Selected Chapters Budapest Fall 2011

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

SCS-MP2 (CP corr.)

MP2 (CP corrected)

Selected Chapters Budapest Fall 2011

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)

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

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

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)

Selected Chapters Budapest Fall 2011

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

+

-

+

-

Left-right correlation

Up-down correlation

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

• 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

SCF: N2.02

MP2:N1.33

Selected Chapters Budapest Fall 2011

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

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

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

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

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

Selected Chapters Budapest Fall 2011

Selected Chapters Budapest Fall 2011

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

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

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

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.

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

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

*No symmetry

Selected Chapters Budapest Fall 2011

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

* incl. storing the amplitudes

*Courtesy Dr. Re Suyong

Selected Chapters Budapest Fall 2011

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

Selected Chapters Budapest Fall 2011

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

Selected Chapters Budapest Fall 2011

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

Selected Chapters Budapest Fall 2011

The end

Tinky

Selected Chapters Budapest Fall 2011