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.

## PowerPoint Slideshow about ' Major methods' - xanti

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

• 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

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

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. plane-wave expansion

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 plane-wave expansion

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 plane-wave expansion

Selected Chapters Budapest Fall 2011

a+D plane-wave expansion

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 plane-wave expansion

• 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 plane-wave expansion

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 plane-wave expansion

• 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 Gaussian calculations for small molecules)

Selected Chapters Budapest Fall 2011

Aspirin Gaussian calculations for small molecules) C9H8O4

Selected Chapters Budapest Fall 2011

Timings for aspirin (B-LYP) Gaussian calculations for small molecules) a

a In minutes on a 2.4 GHz Pentium 4

*Slightly decontracted

Selected Chapters Budapest Fall 2011

Sucrose, C Gaussian calculations for small molecules) 12H22O11

Selected Chapters Budapest Fall 2011

Timings for sucrose (B-LYP) Gaussian calculations for small molecules) a

a In minutes on a 2.4 GHz Pentium 4

Selected Chapters Budapest Fall 2011

Taxol Gaussian calculations for small molecules) C47H51NO14

Selected Chapters Budapest Fall 2011

Timings for taxol (B-LYP) Gaussian calculations for small molecules) a

a In minutes on a 2.4 GHz Pentium 4

Selected Chapters Budapest Fall 2011

Scaling with respect to molecular size Gaussian calculations for small molecules) (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) Gaussian calculations for small molecules) n, n=2-15 log-log plot

Selected Chapters Budapest Fall 2011

Scaling with respect to basis set size - (alanine) Gaussian calculations for small molecules) 5,

10

Selected Chapters Budapest Fall 2011

Average computational costs with and without multipole expansion in FTC

Selected Chapters Budapest Fall 2011

Conclusions expansion in FTC

• 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

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

• 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

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

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

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 3a Sets1 MO in water along a line passing through the symmetry axis. Red= Gaussian basis, Green=small PW

Selected Chapters Budapest Fall 2011

• 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

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 reproducing a Gaussian

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 reproducing a Gaussian

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 reproducing a Gaussian

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 reproducing a Gaussian

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 reproducing a Gaussian

• 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 reproducing a Gaussian

• 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. reproducing a Gaussian 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 reproducing a Gaussian

• 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

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 MP2 energya,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 MP2 energy

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

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

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

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 MP2 energy stacking energies

Selected Chapters Budapest Fall 2011

General MP2 energy

• 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 MP2 energy(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 MP2 energy24H12)2a,b

aCounterpoise

corrected

bDistance of sheets in graphite = 3.35 Å

Selected Chapters Budapest Fall 2011

Porphine dimer MP2 energy

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

Selected Chapters Budapest Fall 2011

Circumcoronene dimer (C MP2 energy108H36)

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

Sandwich

P.D. 1/4

P.D. 3/4

Parallel Displaced 1/2

Selected Chapters Budapest Fall 2011

Circumcoronene dimer: MP2 energythe effect of horizontal displacement

Vertical distance: 3.4 Å

Selected Chapters Budapest Fall 2011

Spin-component scaled MP2 MP2 energy

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

SCS-MP2 (CP corr.)

MP2 (CP corrected)

Selected Chapters Budapest Fall 2011

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 job 383.9 min

C72H52N2O24

Selected Chapters Budapest Fall 2011

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 Suyong Re)

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 Suyong Re)

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 Suyong Re)(Svein Saebo)

Selected Chapters Budapest Fall 2011

• 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

+ Suyong Re)

-

+

-

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 Suyong Re)

• 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: N Suyong Re)2.02

MP2:N1.33

Selected Chapters Budapest Fall 2011

Problems with the current code Suyong Re)

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

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. Suyong Re),=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 distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute.

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 distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute.

Selected Chapters Budapest Fall 2011 distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute.

An efficient atomic-orbital based MP2 gradient program distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute.

Selected Chapters Budapest Fall 2011

Efficient AO formulation of MP2 gradients distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute. 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 distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute.

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 distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute.W. 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 distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute.

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. distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute.

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 distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute.

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) distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute. on a 2.8 GHz Pentium 4 Xeon processor

*No symmetry

Selected Chapters Budapest Fall 2011

Cavitand with 2-oxo-1-butanol distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute.*

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 distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute.

• 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 distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute.

Selected Chapters Budapest Fall 2011

Summary distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute.

• 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 distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute.

Selected Chapters Budapest Fall 2011

Saebo-Pulay local correlation distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute.

• 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 distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute.

The end distant pairs; at shorter distances quadrupole, possibly octupole terms may also contribute.

Tinky

Selected Chapters Budapest Fall 2011