^{1,a)}and Anders M. N. Niklasson

^{1,b)}

### Abstract

Born-Oppenheimer molecular dynamics simulations with long-term conservation of the total energy and a computational cost that scales linearly with system size have been obtained simultaneously. Linear scaling with a low pre-factor is achieved using density matrix purification with sparse matrix algebra and a numerical threshold on matrix elements. The extended Lagrangian Born-Oppenheimer molecular dynamics formalism [A. M. N. Niklasson, Phys. Rev. Lett.100, 123004 (2008)10.1103/PhysRevLett.100.123004] yields microcanonical trajectories with the approximate forces obtained from the linear scaling method that exhibit no systematic drift over hundreds of picoseconds and which are indistinguishable from trajectories computed using exact forces.

This work was supported by the LDRD program at Los Alamos National Laboratory and by the U.S. Department of Energy Office of Basic Energy Sciences.

I. INTRODUCTION

II. METHODS

III. RESULTS

IV. CONCLUSIONS

### Key Topics

- Methane
- 13.0
- Molecular dynamics
- 12.0
- Interatomic forces
- 7.0
- Lagrangian mechanics
- 7.0
- Algebras
- 6.0

## Figures

Total energy versus time for an extended Lagrangian Born-Oppenheimer MD (XL BOMD) trajectory of liquid methane with 1 SCF cycle per 0.25 fs time step and regular Born-Oppenheimer MD trajectories with 1, 2, 4, and 10 SCF cycles per time step that employ the partial charges computed at each time step as a starting guess for the SCF procedure in the next.

Total energy versus time for an extended Lagrangian Born-Oppenheimer MD (XL BOMD) trajectory of liquid methane with 1 SCF cycle per 0.25 fs time step and regular Born-Oppenheimer MD trajectories with 1, 2, 4, and 10 SCF cycles per time step that employ the partial charges computed at each time step as a starting guess for the SCF procedure in the next.

Time per density matrix calculation as a function of the number of orbitals for (a) liquid methane, and (b) polyethylene. Diagonalization and the dense matrix SP2 computations were performed using the Intel MKL DSYEV and DGEMM, respectively on 12 cores of a dual X5660 Xeon node. The sparse matrix SP2 computations were performed on one core of the same node.

Time per density matrix calculation as a function of the number of orbitals for (a) liquid methane, and (b) polyethylene. Diagonalization and the dense matrix SP2 computations were performed using the Intel MKL DSYEV and DGEMM, respectively on 12 cores of a dual X5660 Xeon node. The sparse matrix SP2 computations were performed on one core of the same node.

Per-atom error in the potential energy, Δ*U*, for liquid methane and polyethylene as a function of the number of atoms and the size of the threshold on matrix elements, ɛ, in the scheme.

Per-atom error in the potential energy, Δ*U*, for liquid methane and polyethylene as a function of the number of atoms and the size of the threshold on matrix elements, ɛ, in the scheme.

Per-atom error in the interatomic forces, Δ*F*, for liquid methane (open symbols) and polyethylene (filled symbols). The diamond, circle, and square symbols correspond to thresholds on matrix elements in the scheme of ɛ = 10^{−5}, 10^{−6}, and 10^{−7}, respectively.

Per-atom error in the interatomic forces, Δ*F*, for liquid methane (open symbols) and polyethylene (filled symbols). The diamond, circle, and square symbols correspond to thresholds on matrix elements in the scheme of ɛ = 10^{−5}, 10^{−6}, and 10^{−7}, respectively.

Total energy versus time for extended Lagrangian Born-Oppenheimer MD trajectories with dissipation for liquid methane. (a) Exact forces, (b) sparse matrix SP2 with ɛ = 10^{−5}.

Total energy versus time for extended Lagrangian Born-Oppenheimer MD trajectories with dissipation for liquid methane. (a) Exact forces, (b) sparse matrix SP2 with ɛ = 10^{−5}.

Standard deviation of the total energy about its mean value as a function of the size of the MD time step for MD trajectories computed with exact forces and forces obtained from the scheme with thresholds ɛ = 10^{−5}, 10^{−6}, and 10^{−7}. The time step of 0.25 fs used throughout this work is marked with an arrow.

Standard deviation of the total energy about its mean value as a function of the size of the MD time step for MD trajectories computed with exact forces and forces obtained from the scheme with thresholds ɛ = 10^{−5}, 10^{−6}, and 10^{−7}. The time step of 0.25 fs used throughout this work is marked with an arrow.

Total energy as a function of simulation time for MD simulations of a C_{100}H_{202} polyethylene molecule in the gas phase with exact forces computed with the SP2 algorithm using dense matrix algebra and sparse matrix algebra with ɛ = 10^{−5}. The systematic drift in the calculation is less than 0.05 μeV atom^{−1} ps^{−1}.

Total energy as a function of simulation time for MD simulations of a C_{100}H_{202} polyethylene molecule in the gas phase with exact forces computed with the SP2 algorithm using dense matrix algebra and sparse matrix algebra with ɛ = 10^{−5}. The systematic drift in the calculation is less than 0.05 μeV atom^{−1} ps^{−1}.

Time average pair distribution functions for gas phase C_{100}H_{202} computed from MD trajectories with exact forces and those from calculations with ɛ = 10^{−5}. The plots have been offset for clarity.

Time average pair distribution functions for gas phase C_{100}H_{202} computed from MD trajectories with exact forces and those from calculations with ɛ = 10^{−5}. The plots have been offset for clarity.

Time average radial distribution functions for liquid methane computed from MD trajectories with exact forces and those from calculations with ɛ = 10^{−5}. The plots have been offset for clarity.

Time average radial distribution functions for liquid methane computed from MD trajectories with exact forces and those from calculations with ɛ = 10^{−5}. The plots have been offset for clarity.

## Tables

Levels of sparsity as measured using the fraction of non-zero elements in the self-consistent density matrices for liquid methane and C_{166}H_{332} computed using dense matrix algebra and sparse matrix algebra with ɛ = 10^{−5}, 10^{−6}, and 10^{−7}.

Levels of sparsity as measured using the fraction of non-zero elements in the self-consistent density matrices for liquid methane and C_{166}H_{332} computed using dense matrix algebra and sparse matrix algebra with ɛ = 10^{−5}, 10^{−6}, and 10^{−7}.

Average measures of the error in the self-consistent density matrices for polyethylene and liquid methane computed with the SP2 algorithm using sparse matrix algebra with a numerical threshold on matrix elements, ɛ, and with dense matrix algebra.

Average measures of the error in the self-consistent density matrices for polyethylene and liquid methane computed with the SP2 algorithm using sparse matrix algebra with a numerical threshold on matrix elements, ɛ, and with dense matrix algebra.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content