Major methods. The J matrix engine (Ahmadi & Almlöf, 1995; White & HeadGordon, 1996; Shao &, HeadGordon, 2000) Density fitting (RI) (Dunlap & Rösch, Ahlrichs) Numerical solution of the Poisson equation (Delley, Becke, Termath & Handy, Baker, Nemeth & Pulay, unpublished)
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.
Selected Chapters Budapest Fall 2011
r1r21 = (1/22)dk k2 exp[ik( r1  r2)]
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
ADVANTAGES vs. PURE PLANE WAVES
DISADVANTAGES
Selected Chapters Budapest Fall 2011
The application of plane wave expansion for the solution of the Coulomb problem in Gaussian basis was pioneered by the Parrinello group:
Selected Chapters Budapest Fall 2011
Let us assume first that the whole density can be represented by plane waves
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
a+D
a
molecule
ghost
L extended
L+D 2L
ghost...
Selected Chapters Budapest Fall 2011
(J0, J1: Bessel functions of the first kind; K0,K1: of the 2nd kind)
Selected Chapters Budapest Fall 2011
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 planewave expansion. E.g. contributions arising from (cd) [c=compact, d=diffuse charge distribution (AO product)] can be evaluated because the Coulomb operator is diagonal in kspace. However, the subsequent numerical quadrature requires that c> is subjected to a lowpass filter operation.
Currently: multipole approximation for (cc’) contributions
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
a In minutes on a 2.4 GHz Pentium 4
*Slightly decontracted
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
a In minutes on a 2.4 GHz Pentium 4
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
a In minutes on a 2.4 GHz Pentium 4
Selected Chapters Budapest Fall 2011
627 min
30.2 min
5.8 min
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
10
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
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
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
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.41E4 1.45E72.02E10
min(1S)2k asymptotically
Even in a 1D ring, if only nearest neighbors overlap, the overlap matrix is singular for even n if S0.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 BoysBernardi 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
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
Largest PW density
0.167 PW/a0 (0.3/Å
Largest PW density
0.133 PW/a0
Selected Chapters Budapest Fall 2011
Integrals for electronmolecule scattering cross sections
Selected Chapters Budapest Fall 2011
GaugeIncluding Atomic Orbitals (GIAOs)Polynomial(x,y,z) exp[a(r – A)2] exp[(i/2c) BAr]
Selected Chapters Budapest Fall 2011
The product of a Gaussian and a plane wave gives a Gaussian with complex center coordinates (G. G. Hall; also GaugeIncluding Atomic Orbitals needed for magnetic properties)
Integral types (g = Gaussian, w = plane wave):
(gggg) : use traditional Gaussian techniques (Rys polynomial, ObaraSaiko, McMurchieDavidson etc.) for complex coord
(gwg’g”), (gwg’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”),(gww’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 2D integrals (simple but numerous)
Change of variable: t2 = u2/(+u2), = pq/(p+q)
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
HeadGordon, 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 (pqrs)=Xqs
Y=CTXC (matrix mult.)
store Yij=(pirj) on disk
end do r
end do p
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
Hexapep=Nformyl pentalanine amide
For wellconditioned basis sets, the default threshold guarantees accuracy to a few Eh.
Selected Chapters Budapest Fall 2011
a In minutes on an 800 MHz Pentium III bTetramethoxycalix[4]arene
c In minutes on a 3 GHz Pentium 4 Xeon
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
Master: spawn slaves and bin writer daemons
Slave 1
Read halftransformed integrals from disk and sort them in bins by i,j. Send full bins to the appropriate nodered or blue
Slave 2
Read halftransformed 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
Calix[4]arene:
ccpVTZ basis,
1528 BF, C2h symmetry
Elapsed time on 24 1 GHz Pentium IIIs: 150.4 min for MP2
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
aCounterpoise
corrected
bDistance of sheets in graphite = 3.35 Å
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
Sandwich
P.D. 1/4
P.D. 3/4
Parallel Displaced 1/2
Selected Chapters Budapest Fall 2011
Vertical distance: 3.4 Å
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
SCSMP2 (CP corr.)
MP2 (CP corrected)
Selected Chapters Budapest Fall 2011
631G* 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
631G* 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
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
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
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
+

+

Leftright correlation
Updown correlation
Two sources of savings
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
SCF: N2.02
MP2:N1.33
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
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 = (iajb)/(a + b i  j)
Invoke the multipole expansion but not for approximating integrals: 1/r12 = R3[r1r2 3R2(r1R)(r2R) + …]
Substitute into the pair wavefunction ij = a,bTijabijab
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; ris 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 23 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, 1020 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
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
Selected Chapters Budapest Fall 2011
Based on the orbitalinvariant form of the MP2 energy
P. Pulay and S. Saebo, Orbitalinvariant formulation and gradient evaluation in MøllerPlesset Perturbation Theory, Theor. Chim. Acta 69, 357 (1986).
Also used to derived the equations for local MP2 gradient:
A. ElAzhary, G. Rauhut, P. Pulay and H.J. Werner, Analytical Energy Gradients for Local SecondOrder MollerPlesset Perturbation Theory, J. Chem. Phys., 108, 5185 (1998).
Selected Chapters Budapest Fall 2011
Normalized to 4 (ab) or to 2 (a=b), and pairwise nonorthogonal,
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 ClosedShell SelfConsistent Electron Pair Theory, J. Chem. Phys. 81, 1901 (1984).
Selected Chapters Budapest Fall 2011
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
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”=HartreeFock
Selected Chapters Budapest Fall 2011
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
It can be derived in several ways. For gradient evaluation the Hylleraas form is particularly useful:
Ec = 21HE00  1H0E01 = 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 (contravariant) amplitudes gives the following equation
Selected Chapters Budapest Fall 2011
*No symmetry
Selected Chapters Budapest Fall 2011
C32H32O10, 631G*, 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
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
Selected Chapters Budapest Fall 2011
Tinky
Selected Chapters Budapest Fall 2011