Computational methods for first-principles molecular dynamics with linear-scaling density functional theory

Abstract

Nowadays, Kohm-Sham density functional theory (KS-DFT) calculations are routinely employed in several research fields, due to the ability of KS-DFT to provide great accuracy for a wide class of molecular systems and materials. Unfortunately, conventional KS-DFT calculations, although very powerful, require a computational cost that goes with the cube of the system size N, also known as O (N³) scaling, undermining in practice the study of large systems. The advent of linear-scaling, or O (N), DFT (LS-DFT) methods, relying on the locality of the electronic matter, has enabled calculations on increasingly large systems, up to tens of thousands of atoms. The central tenet of linear-scaling methods is the exponential decay in real space of the single-particle reduced density matrix. This property allows to enforce localisation constraints on the electronic structure, significantly reducing the size of the matrices, such as the Hamiltonian matrix, and increasing their sparsity. The single-particle density matrix in the LS-DFT formalism is expanded in terms of a set of atom-centred, strictly localised functions. Employing periodic boundary conditions (PBCs), the energy is minimised with respect to all the degrees of freedom in the density matrix, which allows to attain chemical accuracy using a high-resolution minimal basis set. The combination of localisation constraints and sparse algebra form the substrate for O (N) calculations. In this thesis, we used ONETEP, a linear-scaling DFT program, to carry out our calculations. The aim of our research is to combine molecular dynamics simulations, within the Born-Oppenheimer approximation (BOMD), with linear-scaling DFT methods. In particular, our main goal is to advance current methodology by developing new algorithms to better exploit locality in BOMD and to reduce the computational load while maintaining DFT accuracy.

Dipole moment autocorrelation functions can directly be employed to obtain the IR spectrum of a given molecular system. DFT-MD simulations offer the perfect tool to generate accurate autocorrelation functions which automatically take into account the anharmonicity of the potential energy surface and temperature effects. Computational IR spectroscopy plays a pivotal role in the understanding of conformational changes of biomolecules, which tend to show several almost-degenerate conformers at room temperatures (floppy molecules). It is particularly valuable when interpreting their fingerprint in solution in combination with experimental spectra. We have implemented two algorithms for the computation of the local electronic dipole moments of molecules in solution. Both methods demand a strategy to partition the density. These methods enable the computation of IR spectra of large molecules, such as polypeptide, in explicit solvent. In the resulting IR spectrum the effect of the solvent on the target molecule is automatically captured, whereas its IR signature is removed. We expect these new functionalities to be very helpful in the understanding of how bio-molecules interact with the solvent at room temperature and the effect of these interactions on conformational changes.

Computationally, the most demanding step in molecular dynamics simulations is the evaluation of energies and forces. This has particular severe consequences on BOMD-based approaches. In fact, a self-consistent field (SCF) step is required at each MD step, which in turn requires multiple energy evaluations. As a consequence, the SCF loop has a major effect on the computational load and overall wall-time. MD Schemes that are capable to by-pass the SCF loop altogether, e.g. Car-Parrinello MD or fixed charge force fields in classical MD, are inherently faster in terms of wall-time per MD step, although they usually demand a much smaller time-step. Moreover, the quality of the converged density matrix is crucial for energy conservation and forces in the LS-DFT BOMD approaches. In theory, the self-consistent solution does not depend on the initial guess. In practice, the SCF optimisation is always incomplete, leading to memory effects and the breaking of time-reversal symmetry, which gives rise to systematic errors in energy gradients that manifest as a drift in microcanonical energy. To ameliorate this problem, we present two integration schemes based on an extended-Lagrangian (XL) approach which introduces extended or auxiliary electronic degrees of freedom to generate good quality time-reversible initial guesses in the SCF loop. Both schemes are improvements over the original XL formulation, which suffered from numerical noise accumulation. The first approach, known as dissipative XL (dXL), introduces a dissipative-like term in the Verlet integration (modified Verlet) of the auxiliary degrees of freedom; the second approach, known as inertial XL (iXL), introduces a thermostat, hence requiring a velocity Verlet integrator. We have implemented both schemes in ONETEP and studied their performance using liquid water as test case. In collaboration with the authors of the schemes, we have analysed the energy drift in both classical polarisable force field MD calculations and LS-DFT BOMD. We have found that both schemes are very efficient in reducing the number of SCF iterations while maintaining good energy drifts and have similar performance. We believe that our implementation and analysis will be very beneficial for future applications both with LS-DFT BOMD and classical polarisable force field MD since both schemes constitute important algorithmic improvements that markedly extend the timescales accessible to classical and LS-DFT MD simulations alike.

More information

Identifiers

ORCID for Chris-Kriton Skylaris: orcid.org/0000-0003-0258-3433

Catalogue record

Export record