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.

The performance of density functional theory methods for the modeling of condensed aqueous systems is hard to predict and validation by ab initio molecular simulation of liquid water is absolutely necessary. In order to assess the reliability of these tests, the effect of temperature on the structure and dynamics of liquid water has been characterized with 16 simulations of 20 ps in the temperature range of 280–380 K. We find a pronounced influence of temperature on the pair correlation functions and on the diffusion constant including nonergodic behavior on the time scale of the simulation in the lower temperature range (which includes ambient temperature). These observations were taken into account in a consistent comparison of a series of density functionals (BLYP, PBE, TPSS, OLYP, HCTH120, HCTH407). All simulations were carried out using an ab initio molecular dynamics approach in which wave functions are represented using Gaussians and the density is expanded in an auxiliary basis of plane waves. Whereas the first three functionals show similar behavior, it is found that the latter three functionals yield more diffusive dynamics and less structure.

A series of first principles molecular dynamics and Monte Carlo simulations were carried out for liquid water to investigate the reproducibility of different sampling approaches. These simulations include Car−Parrinello molecular dynamics simulations using the program cpmd with different values of the fictitious electron mass in the microcanonical and canonical ensembles, Born−Oppenheimer molecular dynamics using the programs cpmd and cp2k in the microcanonical ensemble, and Metropolis Monte Carlo using cp2k in the canonical ensemble. With the exception of one simulation for 128 water molecules, all other simulations were carried out for systems consisting of 64 molecules. Although the simulations yield somewhat fortuitous agreement in structural properties, analysis of other properties demonstrate that one should exercise caution when assuming the reproducibility of Car−Parrinello and Born−Oppenheimer molecular dynamics simulations for small system sizes in the microcanonical ensemble. In contrast, the molecular dynamics and Monte Carlo simulations in the canonical ensemble appear to be more reliable. Furthermore, in the case of canonical Car−Parrinello molecular dynamics simulations the application of Nosé−Hoover chain thermostats allows the use of larger fictitious electron masses. For the Becke−Lee−Yang−Parr exchange and correlation energy functionals and norm-conserving Troullier−Martins or Goedecker−Teter−Hutter pseudopotentials, these simulations at a fixed density of 1.0 g/cm3 and a temperature close to 315 K point to an overstructured liquid with a height of the first peak in the oxygen−oxygen radial distribution function of about 3.0, an underestimated value of the classical constant-volume heat capacity of about 70 J/(mol K), and an underestimated self-diffusion constant of about 0.04 Å2/ps.

We present a new method which combines Car-Parrinello and Born-Oppenheimer molecular dynamics in order to accelerate density functional theory based ab initio simulations. Depending on the system a gain in efficiency of 1 to 2 orders of magnitude has been observed, which allows ab initio molecular dynamics of much larger time and length scales than previously thought feasible. It will be demonstrated that the dynamics is correctly reproduced and that high accuracy can be maintained throughout for systems ranging from insulators to semiconductors and even to metals in condensed phases. This development considerably extends the scope of ab initio simulations.

A popular strategy for simulating large systems where quantum chemical effects are important is the use of mixed quantum mechanical/molecular mechanics methods (QM/MM). While the cost of solving the Schrödinger equation in the QM part is the bottleneck of these calculations, evaluating the Coulomb interaction between the QM and the MM part is surprisingly expensive. In fact it can be just as time-consuming as solving the QM part. We present here a novel real space multigrid approach that handles Coulomb interactions very effectively and implement it in the CP2K code. This novel scheme cuts the cost of this part of the calculation by 2 orders of magnitude. The method does not need very fine-tuning or adjustable parameters, and it is quite accurate, leading to a dynamics with very good energy conservation. We exemplify the validity of our algorithms with simulations of water and of a zwitterionic dipeptide solvated in water.

Efficient Monte Carlo algorithms and a mixed-basis set electronic structure program were used to compute from first principles the vapor-liquid coexistence curve of water. A water representation based on the Becke-Lee-Yang-Parr exchange and correlation functionals yields a saturated liquid density of 900 kg/m3 at 323 K and normal boiling and critical temperatures of 350 and 550 K, respectively. An analysis of the structural and electronic properties of the saturated liquid phase shows an increase of the asymmetry of the local hydrogen-bonded structure despite the persistence of a 4-fold coordination and decreases of the molecular dipole moment and of the spread of the lowest unoccupied molecular orbital with increasing temperature.

A new linear-scaling method based on a multigrid approach to treat long-range electrostatic interactions in hybrid quantum mechanics/molecular mechanics (QM/MM) simulations is described. The scheme has been implemented in the context of a QM calculation based on density functional theory (DFT). The method is tested on an analytical model to validate the new algorithm. Two realistic problems in α-quartz crystals and a zwitterionic dipeptide (GLY-ALA) in water have been chosen as further tests. Results from QM/MM calculations with periodic boundary conditions (PBC) show that the use of PBC is essential when studying highly ordered crystal structures, unless a carefully designed MM crystal is used for the calculation. With a general shaped MM subsystem, the absence of PBC leads to an incorrect description of Kohn-Sham band gaps and charge density. The present method allows periodic boundary conditions to be used in molecular simulations of biological and material science systems.

A series of first-principles Monte Carlo simulations in the isobaric-isothermal ensemble were carried out for liquid water at ambient conditions (T=298 K and p=1 atm). The Becke-Lee-Yang-Parr (BLYP) exchange and correlation energy functionals and norm-conserving Goedecker-Teter-Hutter (GTH) pseudopotentials were employed with the CP2 K simulation package to examine systems consisting of 64 water molecules. The fluctuations in the system volume encountered in simulations in the isobaric-isothermal ensemble require a reconsideration of the suitability of the typical charge-density cutoff and the regular grid-generation method previously used for the computation of the electrostatic energy in first-principles simulations in the microcanonical or canonical ensembles. In particular, it is noted that a much higher cutoff is needed and that the most computationally efficient method of creating grids can result in poor simulations. Analysis of the simulation trajectories using a very large charge-density cutoff at 1200 Ry and four different grid-generation methods point to a significantly underestimated liquid density of about 0.8 g cm-3 resulting in a somewhat understructured liquid (with a value of about 2.7 for the height of the first peak in the oxygen-oxygen radial distribution function) for BLYP-GTH water at ambient conditions. In addition, a simulation using a charge-density cutoff at 280 Ry yields a higher density of 0.9 g cm-3, showing the sensitivity of the simulation outcome to this parameter.