Massively parallel linear-scaling Hartree-Fock exchange and hybrid exchange-correlation functionals with plane wave basis set accuracy.

We extend our linear-scaling approach for the calculation of Hartree-Fock exchange energy using localized in situ optimized orbitals [Dziedzic et al., J. Chem. Phys. 139, 214103 (2013)] to leverage massive parallelism. Our approach has been implemented in the onetep (Order-N Electronic Total Energy Package) density functional theory framework, which employs a basis of non-orthogonal generalized Wannier functions (NGWFs) to achieve linear scaling with system size while retaining controllable near-complete-basis-set accuracy. For the calculation of Hartree-Fock exchange, we use a resolution-of-identity approach, where an auxiliary basis set of truncated spherical waves is used to fit products of NGWFs. The fact that the electrostatic potential of spherical waves (SWs) is known analytically, combined with the use of a distance-based cutoff for exchange interactions, leads to a calculation cost that scales linearly with the system size. Our new implementation, which we describe in detail, combines distributed memory parallelism (using the message passing interface) with shared memory parallelism (OpenMP threads) to efficiently utilize numbers of central processing unit cores comparable to, or exceeding, the number of atoms in the system. We show how the use of multiple time-memory trade-offs substantially increases performance, enabling our approach to achieve superlinear strong parallel scaling in many cases and excellent, although sublinear, parallel scaling otherwise. We demonstrate that in scenarios with low available memory, which preclude or limit the use of time-memory trade-offs, the performance degradation of our algorithm is graceful. We show that, crucially, linear scaling with system size is maintained in all cases. We demonstrate the practicability of our approach by performing a set of fully converged production calculations with a hybrid functional on large imogolite nanotubes up to over 1400 atoms. We finish with a brief study of how the employed approximations (exchange cutoff and the quality of the SW basis) affect the calculation walltime and the accuracy of the obtained results.

The ONETEP linear-scaling density functional theory program.

We present an overview of the onetep program for linear-scaling density functional theory (DFT) calculations with large basis set (plane-wave) accuracy on parallel computers. The DFT energy is computed from the density matrix, which is constructed from spatially localized orbitals we call Non-orthogonal Generalized Wannier Functions (NGWFs), expressed in terms of periodic sinc (psinc) functions. During the calculation, both the density matrix and the NGWFs are optimized with localization constraints. By taking advantage of localization, onetep is able to perform calculations including thousands of atoms with computational effort, which scales linearly with the number or atoms. The code has a large and diverse range of capabilities, explored in this paper, including different boundary conditions, various exchange-correlation functionals (with and without exact exchange), finite electronic temperature methods for metallic systems, methods for strongly correlated systems, molecular dynamics, vibrational calculations, time-dependent DFT, electronic transport, core loss spectroscopy, implicit solvation, quantum mechanical (QM)/molecular mechanical and QM-in-QM embedding, density of states calculations, distributed multipole analysis, and methods for partitioning charges and interactions between fragments. Calculations with onetep provide unique insights into large and complex systems that require an accurate atomic-level description, ranging from biomolecular to chemical, to materials, and to physical problems, as we show with a small selection of illustrative examples. onetep has always aimed to be at the cutting edge of method and software developments, and it serves as a platform for developing new methods of electronic structure simulation. We therefore conclude by describing some of the challenges and directions for its future developments and applications.

DL_MG: A Parallel Multigrid Poisson and Poisson-Boltzmann Solver for Electronic Structure Calculations in Vacuum and Solution.

The solution of the Poisson equation is a crucial step in electronic structure calculations, yielding the electrostatic potential-a key component of the quantum mechanical Hamiltonian. In recent decades, theoretical advances and increases in computer performance have made it possible to simulate the electronic structure of extended systems in complex environments. This requires the solution of more complicated variants of the Poisson equation, featuring nonhomogeneous dielectric permittivities, ionic concentrations with nonlinear dependencies, and diverse boundary conditions. The analytic solutions generally used to solve the Poisson equation in vacuum (or with homogeneous permittivity) are not applicable in these circumstances, and numerical methods must be used. In this work, we present DL_MG, a flexible, scalable, and accurate solver library, developed specifically to tackle the challenges of solving the Poisson equation in modern large-scale electronic structure calculations on parallel computers. Our solver is based on the multigrid approach and uses an iterative high-order defect correction method to improve the accuracy of solutions. Using two chemically relevant model systems, we tested the accuracy and computational performance of DL_MG when solving the generalized Poisson and Poisson-Boltzmann equations, demonstrating excellent agreement with analytic solutions and efficient scaling to ∼109 unknowns and 100s of CPU cores. We also applied DL_MG in actual large-scale electronic structure calculations, using the ONETEP linear-scaling electronic structure package to study a 2615 atom protein-ligand complex with routinely available computational resources. In these calculations, the overall execution time with DL_MG was not significantly greater than the time required for calculations using a conventional FFT-based solver.

Electronically Excited States in Solution via a Smooth Dielectric Model Combined with Equation-of-Motion Coupled Cluster Theory.

We present a method for computing excitation energies for molecules in solvent, based on the combination of a minimal parameter implicit solvent model and the equation-of-motion coupled-cluster singles and doubles method (EOM-CCSD). In this method, the solvent medium is represented by a smoothly varying dielectric function, constructed directly from the quantum mechanical electronic density using only two tunable parameters. The solvent-solute electrostatic interactions are computed by numerical solution of the nonhomogeneous Poisson equation and incorporated at the Hartree-Fock stage of the EOM-CCSD calculation by modification of the electrostatic potential. We demonstrate the method by computing excited state transition energies and solvent shifts for several small molecules in water. Results are presented for solvated H2O, formaldehyde, acetone, and trans-acrolein, which have low-lying n → π* transitions and associated blue shifts in aqueous solution. Comparisons are made with experimental data and other theoretical approaches, including popular implicit solvation models and QM/MM methods. We find that our approach provides surprisingly good agreement with both experiment and the other models, despite its comparative simplicity. This approach only requires modification of the Fock operator and total energy expressions at the Hartree-Fock level-solvation effects enter into the EOM-CCSD calculation only through the Hartree-Fock orbitals. Our model provides a theoretically and computationally simple route for accurate simulations of excited state spectra of molecules in solution, paving the way for studies of larger and more complex molecules.

Use of the rVV10 Nonlocal Correlation Functional in the B97M-V Density Functional: Defining B97M-rV and Related Functionals.

The VV10 and rVV10 nonlocal correlation functionals are consistently implemented and assessed, with the goal of determining if the rVV10 nonlocal correlation functional can replace the VV10 nonlocal correlation functional in the recently developed B97M-V density functional, to give the B97M-rV density functional. Along the way, four density functionals are simultaneously tested: VV10, rVV10, B97M-V, and B97M-rV. An initial assessment is carried out across the S22 data set, and the short-range damping variable, b, is varied for all four density functionals in order to determine the sensitivity of the functionals to the empirical parameter. The results of this test indicate that a value of b = 6 (fortuitously the same as that in B97M-V) is suitable for B97M-rV. The functionals are then compared across an extensive database of interaction energies, and it is demonstrated that B97M-rV either matches or outperforms B97M-V for all of the tests considered. Finally, the optimization of b across the S22 data set is extended to two range-separated hybrid density functionals, ωB97X-V and ωB97M-V, and a value of b = 6.2 is recommended for both ωB97X-rV and ωB97M-rV.

Self-consistent implementation of meta-GGA functionals for the ONETEP linear-scaling electronic structure package.

Accurate and computationally efficient exchange-correlation functionals are critical to the successful application of linear-scaling density functional theory (DFT). Local and semi-local functionals of the density are naturally compatible with linear-scaling approaches, having a general form which assumes the locality of electronic interactions and which can be efficiently evaluated by numerical quadrature. Presently, the most sophisticated and flexible semi-local functionals are members of the meta-generalized-gradient approximation (meta-GGA) family, and depend upon the kinetic energy density, τ, in addition to the charge density and its gradient. In order to extend the theoretical and computational advantages of τ-dependent meta-GGA functionals to large-scale DFT calculations on thousands of atoms, we have implemented support for τ-dependent meta-GGA functionals in the ONETEP program. In this paper we lay out the theoretical innovations necessary to implement τ-dependent meta-GGA functionals within ONETEP's linear-scaling formalism. We present expressions for the gradient of the τ-dependent exchange-correlation energy, necessary for direct energy minimization. We also derive the forms of the τ-dependent exchange-correlation potential and kinetic energy density in terms of the strictly localized, self-consistently optimized orbitals used by ONETEP. To validate the numerical accuracy of our self-consistent meta-GGA implementation, we performed calculations using the B97M-V and PKZB meta-GGAs on a variety of small molecules. Using only a minimal basis set of self-consistently optimized local orbitals, we obtain energies in excellent agreement with large basis set calculations performed using other codes. Finally, to establish the linear-scaling computational cost and applicability of our approach to large-scale calculations, we present the outcome of self-consistent meta-GGA calculations on amyloid fibrils of increasing size, up to tens of thousands of atoms.

Density fitting for three-electron integrals in explicitly correlated electronic structure theory.

The principal challenge in using explicitly correlated wavefunctions for molecules is the evaluation of nonfactorizable integrals over the coordinates of three or more electrons. Immense progress was made in tackling this problem through the introduction of a single-particle resolution of the identity. Decompositions of sufficient accuracy can be achieved, but only with large auxiliary basis sets. Density fitting is an alternative integral approximation scheme, which has proven to be very reliable for two-electron integrals. Here, we extend density fitting to the treatment of all three-electron integrals that appear at the MP2-F12/3*A level of theory. We demonstrate that the convergence of energies with respect to auxiliary basis size is much more rapid with density fitting than with the traditional resolution-of-the-identity approach.