CECAM/Psi-k Workshop

Local orbitals and linear-scaling ab initio calculations

CECAM, Lyon, 3-7 September 2001

Abstracts for Tuesday, 4 September 2001

CRYSTAL [1] is a periodic ab initio program that expands the Crystalline Orbitals (CO) in Bloch functions built from local functions ("atomic orbitals", AOs). Each AO is a contraction of Gaussian Type Functions (GTFs), which in turn are the product of a Gaussian and a real solid spherical harmonic. It can treat systems periodic in one, two and three dimensions, and solve both the Hartree-Fock (HF) and Kohn-Sham (DFT) equations. All the most popular local and non-local exchange and correlation functionals can be used in any combination; also "hybrid" schemes, such as the so called B3-LYP, which combines the Hartree-Fock exchange term with the Becke [2] and Lee-Yang-Parr [3] functionals according to the formula proposed by Becke [4], are available.

The Coulomb infinite series are always evaluated analytically, by using for the long range interactions a scheme based on multipolar expansions, spherical harmonics and Hermite polynomial recursion relations and Ewald type summations.

In the DFT calculations two options are now available for the construction of the exchange-correlation matrix elements: either direct numerical integration or the use of an auxiliary basis set of GTFs in which the exchange-correlation potential is expanded. For the numerical integration, the atomic partition method proposed by Becke [5] has been adopted, combined with Gauss-Legendre (radial) and Lebedev (angular) quadrature.

A basis set of symmetry adapted crystalline orbitals (SACO's) is used for obtaining a block-diagonal representation of the Hamiltonian matrix at each k-point in reciprocal space, with a consequent reduction in computational time when large unit cell and high symmetry systems are considered [6,7,8]. In these conditions, it is possible to numerically optimise the geometry of complex structure at a relatively low cost. Recent improvements with respect to the public version of the code [1] include the evaluation of the analytic gradients of the HF Energy with respect to the internal coordinates [9,10,11] and the construction of well localised Wannier functions [12].

Three different examples will be discussed during the workshop:

1. F-centre defects in LiF, by using supercells containing 8, 16, 32, 64, 128, 256 and 512 atoms; the geometry is optimised at the HF level; the convergence of the defect properties with respect to the defect-defect distance is documented.
2. Perfect MgO supercells containing 2-512 atoms; it is shown that all the intensive properties are stable as a function of the supercell size; the scaling properties of the code are discussed;
3. the optimisation of the structure of various zeolites with and without adsorbates.

The following tables and figure refer to the first example, treated at the HF level with a basis set containing 5 and 13 atomic orbitals/atom for Li and F, respectively.

Table 1: Effect of relaxation for an F-centre in LiF as a function of supercell size S_n. In the Li₁₀₀, F₁₁₀ and Li₁₁₁ columns, the variation of the distances of the first three neighbors from the defect is reported; a positive value indicates larger distance. and E_form (in hartree) are the energy gain with respect to the unrelaxed geometry and the defect formation energy. N_var and N_cyc are the number of variables to be optimised and of optimisation cycles.

Sn	Li100	F110	Li111		Eform	Nvar	Ncyc
8	-	-	-	-0.00000	0.25554	0	-
16	0.024	-	-	-0.00063	0.25481	1	5
32	0.039	0.013	0.000	-0.00103	0.25386	2	5
54	0.039	0.011	-0.003	-0.00105	0.25319	6	7
64	0.043	0.015	-0.003	-0.00115	0.25265	6	12
128	0.043	0.014	-0.004	-0.00122	0.25039	12	9
216	0.045	0.011	-0.005	-0.00124	0.24739	20	22
250	0.045	0.011	-0.004	-0.00128	0.24618	24	17
256	0.045	0.016	-0.004	-0.00149	0.24575	21	15

Figure 1: Electron charge (left) and spin (right) density of an F-centre in LiF, investigated at the HF level and using three different supercells: S₈ (top), S₆₄ (middle) and S₂₁₆ (bottom). The section is parallel to a (100) plane through the defect, labelled as XX. The dashed squares indicate the supercell sizes in the plane. The map side is 16.08 Å long. The separation between contiguous isodensity curves is 0.01 and 0.001 e/bohr³, for the electron charge and spin densities, respectively. The density range is 0.01/0.1 (left) and -0.01/0.01 (right) e/bohr³; continuous, dashed and dot-dashed lines denote positive, negative, and zero values, respectively.

1. V R Saunders, R Dovesi, C Roetti, M Causà, N M Harrison, R Orlando and C M Zicovich-Wilson, CRYSTAL98 User's Manual, Università di Torino (Torino, 1998).
2. A D Becke, Phys. Rev. A 38, 3098 (1988).
3. C Lee, W Yang and R G Parr, Phys. Rev. B 37, 785 (1988).
4. A D Becke, "Density-functional thermochemistry. III The role of exact exchange", J. Chem. Phys. 98, 5648 (1993).
5. A D Becke, J. Chem. Phys. 88, 1053 (1988).
6. C Zicovich-Wilson and R Dovesi, "On the use of Symmetry Adapted Crystalline Orbitals in SCF-LCAO periodic calculations. I. The construction of the Symmetrized Orbitals", Int. J. Quantum Chem. 67, 299-309 (1998).
7. C Zicovich-Wilson and R Dovesi, "On the use of Symmetry Adapted Crystalline Orbitals in SCF-LCAO periodic calculations. II. Implementation of the Self-Consistent-Field scheme and examples", Int. J. Quantum Chem. 67, 311-320 (1998).
8. C Zicovich-Wilson and R Dovesi, "Periodic ab-initio study of silico-faujasite", Chem. Phys. Lett. 277, 227-233 (1997).
9. K Doll, "Implementation of analytical Hartree-Fock gradients for periodic systems", Computer Physics Communications 137, 74-88 (2001).
10. K Doll, V R Saunders, N M Harrison, Int. J. Quantum Chem. 82, 1 (2001).
11. V R Saunders, R Orlando and R Dovesi, unpublished.
12. C M Zicovich-Wilson, R Dovesi and V R Saunders, "A general method to obtain well localized Wannier Functions in LCAO periodic calculations", J. Chem. Phys. (2001), submitted.

In recent years progress has been made towards density functional calculations with computational costs proportional to system size [1]. The electrostatic problem is solved by fast multipole or fast Fourier/wavelet transform methods. Direct search methods for the density matrix and projection methods have replaced diagonalization of the Kohn-Sham matrix.

However, accuracy, numerical stability, and performance of these methods have to be further improved. We will discuss these items in connection with our implementation of the Gaussian and Augmented-Plane Wave (GAPW) method [2,3].

In the first part we will discuss two methods applied in the calculation of elements of the Kohn-Sham matrix. Most atomic orbital based methods use fast multipole methods [4,5] to achieve linear scaling in calculating electrostatics or make use of auxiliary basis functions to expand the electron density [6,7]. The GAPW method uses a special combination of all these methods together with techniques transfered from solid state calculations. In a first step projector techniques [8,9] are used to separate the electron density in atomic and interstitial contributions

In the second part the problem of calculating the density matrix from a given Kohn-Sham matrix is discussed. Many different methods have been proposed to solve this task with linear scaling [1]. All of these methods suffer from the slow decay of off-diagonal matrix elements of the density operator. A problem greatly enhanced by ill-conditioned overlap matrices encountered in calculations with larger localized basis sets [11]. The polarized atomic orbital (PAO) method [12,13] provides a partial solution to these problems. However, the accuracy of the calculation is reduced. Test calculations and possible new fields of applications will be presented.

1. S Goedecker, Rev. Mod. Phys. 71, 1085 (1999).
2. G Lippert, J Hutter and M Parrinello, Molec. Phys. 92, 477 (1997).
3. G Lippert, J Hutter and M Parrinello, Theor. Chem. Acc. 103, 124 (1999).
4. C A White, B G Johnson, P M W Gill and M Head-Gordon, Chem. Phys. Lett. 230, 8 (1994).
5. M C Strain, G E Scuseria and M J Frisch, Science 271, 51 (1996).
6. B I Dunlap, J W D Connolly and J R Sabin, J. Chem. Phys. 71, 3396 (1979).
7. O Vahtras, J Almlöf and M W Feyereisen, Chem. Phys. Lett. 208, 359 (1993).
8. O K Andersen, Phys. Rev. B 12, 3060 (1975).
9. P Blöchl, Phys. Rev. B 50, 17953 (1994).
10. M Krack and M Parrinello, Phys. Chem. Chem. Phys. 2, 2105 (2000).
11. P E Maslen, C Ochsenfeld, C A White, M S Lee and M Head-Gordon, J. Phys. Chem. A 102, 2215 (1998).
12. M S Lee and M Head-Gordon, J. Chem. Phys. 107, 9085 (1997).
13. G Berghold, M Parrinello and J Hutter, submitted.

This presentation will address our recent efforts towards developing linear scaling electronic structure methods for the quantum mechanical modelling of large molecules and periodic systems, using gaussian-type orbitals. In particular, we will discuss our recent advances in density functional theory methods and seconder-order perturbation theory (MP2). Time permitting, applications to fluorinated carbon nanotubes will be presented.

In Epstein-Nesbet perturbation theory several aspects of localization and delocalization can be studied. In canonical orbitals (completely delocalized) Epstein-Nesbet perturbation theory presents some defaults, which lead to the fact that nowadays the method is rarely employed. With respect to Møller-Plesset perturbation theory, Epstein-Nesbet perturbation theory needs Coulomb and Exchange integrals between hole-hole, hole-particle, and particle-particle pairs, which, when concentrated on only a few terms, lead to the inclusion of essential interactions found in CEPA-0 or MP4, but with much less effort. In canonical orbitals these integals are dispersed on many individual terms, and, compared to Møller-Plesset perturbation theory, no further correlation correction can be calculated in the limit of large systems. In order to obtain a useful method, a mixture of Møller-Plesset perturbation in canonical orbitals (or the orbital-invariant formulation) may be coupled to the additional interactions from the Epstein-Nesbet perturbation series, the latter being expressed in a localized basis.

A corresponding intermediate between Epstein-Nesbet and Møller-Plesset H₀ is investigated, and applied to a few test cases.

Linear scaling methods based on Gaussian orbitals vary greatly in their maturity and their extent of relevance to chemists who perform routine applications of electronic structure methods. In this talk, I shall focus first on the construction of effective Hamiltonian matrices, which is without doubt the most mature area. The status of fast Coulomb and exchange methods will be reviewed for energies and analytical gradients, and compared with grid-based quadrature for exchange-correlation contributions. In the second part of the talk, some new developments in diagonalization-free methods for updating the one-particle density matrix will be discussed. These methods are at this stage significantly less relevant to most chemical applications. Finally, if time permits, I shall discuss the use of localized orbital methods in describing highly correlated wavefunctions, and an interesting theorem about the structure of the exact wavefunction.