Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach

https://doi.org/10.1016/j.cpc.2004.12.014Get rights and content

Abstract

We present the Gaussian and plane waves (GPW) method and its implementation in Quickstep which is part of the freely available program package CP2K. The GPW method allows for accurate density functional calculations in gas and condensed phases and can be effectively used for molecular dynamics simulations. We show how derivatives of the GPW energy functional, namely ionic forces and the Kohn–Sham matrix, can be computed in a consistent way. The computational cost of computing the total energy and the Kohn–Sham matrix is scaling linearly with the system size, even for condensed phase systems of just a few tens of atoms. The efficiency of the method allows for the use of large Gaussian basis sets for systems up to 3000 atoms, and we illustrate the accuracy of the method for various basis sets in gas and condensed phases. Agreement with basis set free calculations for single molecules and plane wave based calculations in the condensed phase is excellent. Wave function optimisation with the orbital transformation technique leads to good parallel performance, and outperforms traditional diagonalisation methods. Energy conserving Born–Oppenheimer dynamics can be performed, and a highly efficient scheme is obtained using an extrapolation of the density matrix. We illustrate these findings with calculations using commodity PCs as well as supercomputers.

Introduction

Density functional theory [1], [2] (DFT) is a well established method to perform electronic structure calculations. The accuracy of the method is such that many properties of systems of interest to chemistry, physics, material science, and biology can be predicted in a parameter free way. The standard computational approach to DFT is already efficient and thus appropriate for fairly large systems, currently about 100 atoms. Nevertheless, the computation of the Hartree (Coulomb) energy and the orthogonalisation of the wave functions are not scaling linearly with system size, and these terms therefore dominate the computational cost for larger systems [3]. The hybrid Gaussian and plane waves (GPW) method [4] provides an efficient way to treat these terms accurately at a significantly reduced cost. We present here the implementation of this method in Quickstep, which is part of the freely available program package CP2K [5].

The method uses an atom-centred Gaussian-type basis to describe the wave functions, but uses an auxiliary plane wave basis to describe the density. With a density represented as plane waves or on a regular grid, the efficiency of Fast Fourier Transforms (FFT) can be exploited to solve the Poisson equation and to obtain the Hartree energy in a time that scales linearly with the system size. Fast Fourier Transforms and regular grids are well established in plane wave codes [6] and their efficiency has recently been exploited in a similar method [7], [8], [9], [10], [11]. The use of an auxiliary basis set to represent the density goes back to the seventies [12], [13] and has become increasingly popular as resolution of the identity (RI) method or density fitting method. Contrary to the GPW method, most RI methods expand the density in an auxiliary basis of the same nature as the primary basis, but optimised specifically for this purpose [14], [15], [16]. Furthermore, since a density expanded in plane waves can be represented on a real space grid, there is a direct connection to methods that use numerical calculation of matrix elements [17], [18] or grid discretisation and finite element methods (for a recent review see Ref. [19]). The GPW method is most similar to methods that employ auxiliary real space grids but differ by the choice of localised primary basis functions used to represent the wave functions [20], [21], [22], [23], [24], [25].

Periodic boundary conditions follow naturally from the FFT based treatment of the Poisson equation, and the GPW method scales linearly for three-dimensional systems with a small prefactor and an early onset. The GPW method seems therefore best suited for the simulation of large and dense systems, such as liquids and solids, and all recent applications of the method fall in this category [26], [27], [28], [29], [30]. For these systems, it is important to be able to efficiently perform stable molecular dynamics simulations, in order to address finite temperature effects. Plane wave codes and the basic GPW implementation presented here require that the nuclei are described using pseudopotential. This approximation is highly accurate if, e.g., Goedecker–Teter–Hutter (GTH) pseudopotentials are employed [31], [32]. An extension of the GPW method, the Gaussian and augmented-plane-wave (GAPW) method [26], [33] allows for all electron calculations.

The extensive experience with Gaussian-type basis sets shows that basis set sequences that increase rapidly in accuracy can be constructed in a systematic way [34]. At the same time, a compact description of the wave functions is maintained, and this opens the way for efficient methods to solve for the self-consistent field (SCF) equations. Furthermore, as Gaussian functions are localised, the representations of the Kohn–Sham, overlap and density matrix in this basis become sparse with increasing system size [3]. This eventually allows for solving the Kohn–Sham (KS) equations using computational resources that scale linearly with system size. We have currently only implemented methods that are scaling cubically with system size, but these have been designed to reach high efficiency for Gaussian basis sets [35].

In this paper, we provide an up-to-date description of the method and its implementation. We review the GPW energy functional and illustrate its linear scaling nature in Section 2, whereas the derivatives of the functional are described in Section 3. In Section 4 details on the program structure and implementation are provided. The construction of a sequence of systematically improving basis sets is discussed in Section 5. In Section 6, two methods to perform wave function optimisation are presented. Molecular dynamics simulations and improved wave function extrapolation method are the subject of Section 7. The accuracy of the method is illustrated for gas phase and condensed phase systems in Section 8, and the efficiency of the code for serial and parallel calculations is shown in Section 9.

Section snippets

Energy functional

Central in the Gaussian and plane wave (GPW) method [4] is the use of two representations of the electron density. Such a dual representation allows for an efficient treatment of the electrostatic interactions, and leads to a scheme that has a linear scaling cost for the computation of the total energy and Kohn–Sham matrix with respect to the system size. The first representation of the electron density n(r) is based on an expansion in atom centred, contracted Gaussian functions n(r)=∑μνPμνφμ(r)

Deriving the Kohn–Sham matrix from the GPW energy

In this section we present how the exact derivative Hμν=∂E/∂Pμν of the total energy is computed taking into account all approximations that lead from the Gaussian based density n({Pμν}) to the density represented on the grid n˜({Pμν}). This includes the mapping from the Gaussian basis to the grid using finite radii, the use of multi-grids, and of grid based methods to compute ∇n˜(r). We use the notation n˜i to denote a value of n˜(r) on a particular grid point with coordinates ri using a single

Program structure and implementation

Quickstep is tightly integrated in a larger program named CP2K and the Fortran 95 sources are made freely available [5]. Currently, good portability and efficiency is achieved by relying on a few commonly available libraries such as BLAS [59], LAPACK [60], FFTW [61], MPI [62] and ScaLAPACK [63], and restricting ourselves to standard Fortran with OpenMP [64] extensions. In this section, we present our approach to structural aspects of the code, and describe how good efficiency can be obtained

Basis sets

The Kohn–Sham orbitals are expanded in Gaussian orbital functions in the framework of the GPW method as described in Section 2. Significant experience exists with Gaussian basis sets and they are available in a number of formats [67], [68]. Whereas polarisation and diffuse functions can normally be adopted from published basis sets, the valence part of the basis has to be generated for the usage with the GTH pseudopotentials. A systematically improving sequence of basis sets for use with the

Wavefunction optimisation

The calculation of the electronic ground state amounts to the minimisation of the electronic energy (Eq. (3)) with respect to the orthonormal one-particle orbitals or the one-particle density matrix. Even though techniques exist to perform this minimisation in a time scaling linearly with the system size [3], no such technique is currently implemented in Quickstep, and the algorithms show a cubic scaling with respect to the number of atoms in the unit cell. Nevertheless, a careful design of the

Ab initio molecular dynamics

Classical molecular dynamics simulations are well established as a powerful technique to study dynamic and thermodynamic properties of atomic and molecular systems at a finite temperature [42]. Similar studies can be performed with ab initio molecular dynamics in which explicit electronic structure calculations are employed to compute potential energies and forces [6]. A significant advantage of ab initio molecular dynamics is that no parametrisation of an empirical potential is needed and as

Small molecules

As a first accuracy test for Quickstep we have employed the new basis sets described in Section 5 for the geometry optimisation of small molecules using the local density approximation (LDA). The CP2K geometry optimiser works with analytic first derivatives whereas the second derivatives are obtained via an updated Hessian method. In that way, each molecule of the following test set of 39 small molecules:

H2, Li2, LiH, BH3, CH4, C2H2, C2H4, C2H6, N2, NH3, HCN, H2O, H2O2, CO, CO2, CH3OH, N2O, F2

Benchmarks

In this section, we show that the accuracy of Quickstep can be achieved with high computational efficiency. We illustrate both the serial performance and the scalability on parallel computers using a high-end supercomputer as well as a modern cluster based on a PC-like architecture.

Summary and outlook

In this paper, we have described the current implementation of the GPW method in the Quickstep program. The dual representation of the charge density allows for an efficient treatment of the Hartree terms, while maintaining the compact representation of the wave functions in the atomic orbital basis. Furthermore, since the description of the wave function is completely in a localised basis, linear scaling Kohn–Sham matrix construction is obtained using screening techniques that lead to sparse

Acknowledgments

This work was partially supported by the Bundesministerium für Bildung und Forschung (BMBF) in the framework of the project High-Performance Computing in Chemistry (HPC-Chem) [80]. One of the authors (J.V.) obtained financial support from a Marie Curie Fellowship, and would like to thank M. Sprik for supporting development of the code. T. Chassaing and J. Hutter acknowledge support from the Swiss National Science Foundation (Project 200020-100417). Computer resources were in part provided by

References (80)

  • L. Füsti-Molnár et al.

    J. Mol. Struct. (Theochem)

    (2003)
  • O. Vahtras et al.

    Chem. Phys. Lett.

    (1993)
  • V. Termath et al.

    Chem. Phys. Lett.

    (1994)
  • C.M. Goringe et al.

    Comput. Phys. Comm.

    (1997)
  • A.A. Mostofi et al.

    Comput. Phys. Comm.

    (2002)
  • B. Miehlich et al.

    Chem. Phys. Lett.

    (1989)
  • P. Pulay

    Chem. Phys. Lett.

    (1980)
  • P. Pulay et al.

    Chem. Phys. Lett.

    (2004)
  • H.M. Berman et al.

    Biophys. J.

    (1992)
  • P. Hohenberg et al.

    Phys. Rev. B

    (1964)
  • W. Kohn et al.

    Phys. Rev.

    (1965)
  • S. Goedecker

    Rev. Mod. Phys.

    (1999)
  • G. Lippert et al.

    Mol. Phys.

    (1997)
  • D. Marx et al.
  • L. Füsti-Molnár et al.

    J. Chem. Phys.

    (2002)
  • L. Füsti-Molnár et al.

    J. Chem. Phys.

    (2002)
  • L. Füsti-Molnár

    J. Chem. Phys.

    (2003)
  • J. Baker et al.

    J. Phys. Chem. A

    (2004)
  • J.L. Whitten

    J. Chem. Phys.

    (1973)
  • B.I. Dunlap et al.

    J. Chem. Phys.

    (1979)
  • K. Eichorn et al.

    Chem. Phys. Lett.

    (1995)
  • K. Eichorn et al.

    Theor. Chem. Acc.

    (1997)
  • B. Delley

    J. Chem. Phys.

    (1989)
  • T.L. Beck

    Rev. Mod. Phys.

    (2000)
  • X. Chen et al.

    Phys. Rev. B

    (1995)
  • P. Ordejón et al.

    Phys. Rev. B

    (1996)
  • S.D. Kenny et al.

    Phys. Rev. B

    (2000)
  • Y. Liu et al.

    Phys. Rev. B

    (2003)
  • M. Krack et al.

    Phys. Chem. Chem. Phys.

    (2000)
  • G. Hura et al.

    Phys. Chem. Chem. Phys.

    (2003)
  • I.-F.W. Kuo et al.

    Science

    (2004)
  • I.-F.W. Kuo et al.

    J. Phys. Chem. B

    (2004)
  • J. VandeVondele, F. Mohamed, M. Krack, J. Hutter, M. Sprik, M. Parrinello, J. Chem. Phys. (2004)...
  • S. Goedecker et al.

    Phys. Rev. B

    (1996)
  • C. Hartwigsen et al.

    Phys. Rev. B

    (1998)
  • G. Lippert et al.

    Theor. Chem. Acc.

    (1999)
  • T.H. Dunning

    J. Chem. Phys.

    (1989)
  • J. VandeVondele et al.

    J. Chem. Phys.

    (2003)
  • A.D. Becke

    Phys. Rev. A

    (1988)

    • Cited by (0)

    View full text
    Copyright © 2005 Elsevier B.V. All rights reserved.