Today, ab initio molecular dynamics (AIMD) relies on the locality of one-electron density matrices to achieve linear growth of computation time with the system size, crucial in large-scale simulations. While Kohn-Sham orbitals strictly localized within predefined radii can offer substantial computational advantages over density matrices, such compact orbitals are not used in AIMD because a compact representation of the electronic ground state is difficult to find. Here, a robust method for maintaining compact orbitals close to the ground state is coupled with a modified Langevin integrator to produce stable nuclear dynamics for molecular and ionic systems. This eliminates a density matrix optimization and enables first orbital-only linear-scaling AIMD. An application to liquid water demonstrates that low computational overhead of the new method makes it ideal for routine medium-scale simulations, while its linear-scaling complexity allows us to extend first-principle studies of molecular systems to completely new physical phenomena on previously inaccessible length scales.

Since the unification of molecular dynamics and density functional theory (DFT),1ab initio molecular dynamics (AIMD) has become an important tool to study processes in molecules and materials. Unfortunately, the computational cost of the conventional Kohn-Sham (KS) DFT grows cubically with the number of atoms, which severely limits the length scales accessible by AIMD. To address this issue, substantial efforts have been directed to the development of linear-scaling (LS) DFT.

In all LS DFT methods, the delocalized eigenstates of the effective KS Hamiltonian must be replaced with an alternative set of local electronic descriptors. Most LS methods2–5 explore the natural locality of the one-electron density matrix (DM). However, the DM DFT becomes advantageous only for impractically large systems when accurate multifunction basis sets are used.3,6–8 This issue is rectified in optimal-basis DM methods9–11 that contract large basis sets into a small number of new localized functions and then optimize the DM in the contracted basis. Despite becoming the most popular approach to LS DFT, the efficiency of these methods is hampered by the costly optimization of both the contracted orbitals and the DM.12 From this point of view, a direct variation of compact molecular orbitals—orbitals that are strictly localized within predefined regions—is preferable because LS can be achieved with significantly fewer variables. Advantages of the orbitals-only LS DFT are especially pronounced in accurate calculations that require many basis functions per atom. Unfortunately, the development of promising orbital-based LS methods has been all but abandoned13,14 because of the inherently difficult optimization of localized orbitals.2,13–17

Thus, despite impressive progress of the LS description of the electronic and atomic structure of large static systems,6,18 the high computational overhead of existing LS methods restricts their use in dynamical simulations to very short time scales, systems of low dimensions, and low-quality minimal basis sets.6,18–20 On typical length and time scales required in practical and accurate AIMD simulations, LS DFT still cannot compete with the straightforward low-cost cubically-scaling KS DFT.

In this work, we present an AIMD method that overcomes difficulties of DFT based on compact orbitals to achieve LS with extremely low computational overhead. To demonstrate advantages of the new method, we applied it here to systems of weakly interacting molecules. However, the same approach is readily applicable to systems of strongly interacting fragments that do not form strong covalent bonds such as ionic materials—salts, liquids, and semiconductors. A generalization of the method to all finite-gap systems, including covalently bonded atoms, will be reported later.

The new AIMD method utilizes a recently developed LS DFT8 based on absolutely localized molecular orbitals (ALMOs)—compact orbitals first described in Ref. 21. Unlike delocalized KS orbitals, each ALMO has its own localization center and a predefined localization radius Rc that typically includes nearby atoms or molecules.8,21 In the current implementation, a localization center is defined as a set of all Gaussian atomic orbitals of one molecule. However, the approach can use other local and nonlocal basis sets.22,23 The key feature of ALMO DFT is that its one-electron wavefunctions are constructed in a two-stage self-consistent-field (SCF) procedure8 to circumvent the problem of the sluggish variational optimization emphasized above. In the first stage, ALMOs are constrained to their localization centers24 whereas, in the second stage, ALMOs are relaxed to allow delocalization onto the neighbor molecules within their localization radius Rc. To achieve a robust optimization in the problematic second stage, it is important to keep the delocalization component of the trial wavefunction orthogonal to the fixed orbitals obtained in the first stage. For mathematical details, see the ALMO SCF method in Ref. 8. ALMO DFT algorithms are presented in the supplementary material.

ALMO constraints imposed by Rc prohibit electron density transfer between distant molecules but retain all other types of interactions such as long-range electrostatic, exchange, polarization, and—if the exchange-correlation (XC) functional includes them—dispersion interactions.25 Since the importance of electron transfer decays exponentially with distance in finite-gap materials,2 the ALMO approximation is expected to provide a natural and accurate representation of the electronic structure of molecular systems. Because of the greatly reduced number of electronic descriptors and the robust optimization, the computational complexity of ALMO DFT grows linearly with the number of molecules, while its computational overhead remains very low. These features make ALMO DFT a promising method for accurate AIMD simulations of large molecular systems.

The challenge of adopting ALMO DFT for dynamical simulations arises from the slightly nonvariational character of the localized orbitals. While ALMOs are variationally optimized in both SCF stages, the occupied subspace defined in the first stage must remain fixed during the second stage to ensure convergence. In addition, electron transfer effects can suddenly become inactive in the course of a dynamical simulation when a neighboring molecule crosses the localization threshold Rc. Futhermore, the variational optimization in any AIMD method is never complete in practice and interrupted once the maximum norm of the gradient of the energy with respect to the electronic descriptors drops below small but nevertheless finite convergence threshold ϵSCF. These errors do not affect the accuracy of static ALMO DFT calculations, geometry optimization, and Monte Carlo simulations. Unfortunately they tend to accumulate in molecular dynamics trajectories leading to non-physical sampling and eventual failure. Traditional strategies to cope with these problems are computationally expensive and include computing the nonvariational contribution to the forces via a variational coupled-perturbed procedure,4,26 increasing Rc, and decreasing ϵSCF.

In this work, we propose another approach that obviates the need in a coupled-perturbed solver, relaxes tight constraints on Rc and ϵSCF, and thus enables us to maintain stable dynamics and to keep the algorithmic complexity and cost of simulations low. In our approach, the forces on atoms are calculated approximately after the two-stage ALMO SCF using a straightforward procedure that computes only the Hellmann-Feynman and Pulay components and neglects the computationally intense nonvariational component of the forces. The difference between these approximate ALMO forces and the reference forces that could be obtained from perfectly converged fully-delocalized KS orbitals is δf(t),

fiαKS(t)=fiαALMO(t)+δfiα(t),
(1)

where α is a Cartesian component of the force acting on atom i at time t. δf(t) comprises all neglected terms that originate from a finite localization radius Rc and incomplete SCF optimization. δf(t) can be reduced to zero systematically by increasing Rc and decreasing ϵSCF.

Our approach to compensate for the missing δf(t) term is inspired by the methodology introduced into AIMD by Krajewski and Parrinello,27 formalized by Kühne et al.,28 and rationalized by Dai and Yuan29 before becoming informally known as the second generation Car-Parrinello molecular dynamics.30 Adopting the principle of Refs. 27 and 28, ALMO AIMD is chosen to be governed by the Langevin equation of motion that can be written in terms of the unknown reference forces

mir̈iα=fiαKS(t)γmiṙiα+Riαγ(t),
(2)

where mi is the mass of atom i, r is its position along dimension α, γ is the Langevin scaling factor, and Riαγ(t) is the stochastic force represented by a zero-mean white Gaussian noise

Riαγ(t)=0,
(3)
Riαγ(t)Rjβγ(t)=2kBTγmiδijδαβδ(tt).
(4)

The last relation means that, for any value of γ, the damping and stochastic terms are in perfect balance, and trajectories generated with Eq. (2) will sample the canonical ensemble at a specified temperature T.31 In the limit γ → 0, the Newton equation is recovered and the microcanonical ensemble is sampled.

The main assumption of ALMO AIMD is that the error in the ALMO forces is well approximated by Gaussian noise RiαΔ(t)

δfiα(t)=RiαΔ(t)
(5)

that obeys

RiαΔ(t)=0,
(6)
RiαΔ(t)RjβΔ(t)=2kBTΔmiδijδαβδ(tt).
(7)

This assumption, shown to be well justified, allows us to rewrite the Langevin equation using the ALMO forces

mir̈iα=fiαALMO(t)γmiṙiα+Riαγ+Δ(t),
(8)

where the two stochastic terms are combined into one Riαγ+Δ=Riαγ+RiαΔ. The only missing piece in the modified Langevin equation is the value of Δ, which describes the strength of the newly introduced stochastic term. This term compensates for imperfections in ALMO forces and must be adjusted to re-balance the damping and stochastic components in ALMO AIMD.

In principle, Δ can be calculated using the integral of Eq. (7) averaged over atoms with different mi

Δ=(2kBTmi)113δfi(0)δfi(τ)dτ,
(9)

if one can afford computing the reference forces fiαKS(t) (i.e., Rc and ϵSCF → 0) for a short representative AIMD trajectory. In practice, we found (see results below) that this approach is not sufficiently accurate because the δijδαβ assumption in Eq. (7) does not strictly hold. Nevertheless, the integral of the autocorrelation function (ACF) in Eq. (9) can provide a reasonable starting value of Δ. The estimated Δ can be further fine-tuned in a series of short trial-and-error ALMO AIMD runs until the average kinetic energy corresponds to the requested temperature, 12miṙi2=32kBT.

The inherently stochastic approach presented here does not aim to produce fully time-reversible dynamics for atomic nuclei. Nevertheless, it is capable to reproduce correct dynamical properties of a system as long as γ is set to a small value and partially optimized ALMOs remain close to the ground state resulting in Δ ≪ γ. An alternative approach to performing microcanonical MO-based LS AIMD would be to combine ALMO DFT with the extended Lagrangian method described in Ref. 32. This method has been shown to conserve the total energy in LS DM-based AIMD simulations.33 

ALMO AIMD was implemented in CP2K, an open source materials modeling package.34 Accuracy and efficiency of ALMO AIMD was tested using liquid water as an example. This system is challenging because intermolecular electron delocalization is a critical component of hydrogen bonding and must be described correctly to reproduce static and dynamical properties of liquid water. A periodic cell containing 125 molecules was simulated for 30 ps at T = 298 K and a constant density of 1.01 g cm−3. The initial configurations of the water boxes were created by equilibrating the system for 20 ps using fully converged orbitals. The Ricci-Ciccotti algorithm35 was used to integrate the Langevin equation (see the supplementary material). We found that γ = 10−3 fs−1 is large enough to thermostat the system efficiently and small enough not to significantly affect dynamical properties of liquid water. In the dual Gaussian and plane-wave scheme implemented in CP2K,36 the TZV2P basis set was used to represent molecular orbitals, and a plane-wave cutoff of 320 Ry was used to represent electron density. The XC energy was approximated using the dispersion-corrected Perdew-Burke-Ernzerhof (PBE) functional.37,38 Separable norm-conserving pseudopotentials were used,39 and the Brillouin zone was sampled at the Γ-point. The predictor of the Kolafa scheme40 was adopted to localized orbitals28 to generate a highly accurate initial ALMOs in both SCF stages, which can be brought close to the ground state with just a few SCF steps of the robust two-stage optimization procedure.

The reference forces were calculated with fully delocalized electrons using the tightly converged, ϵSCF = 10−6 a.u., orbital transformation (OT) method.41 In ALMO AIMD, the element-specific cutoff radius for electron delocalization Rc was set to 1.6 in units of the elements’ van der Waals radii (vdWR). This localization radius includes approximately two coordination shells of an average water molecule and was shown to reproduce the reference radial distribution function (RDF) perfectly in Monte Carlo simulations.8 To check the ability of the RΔ(t) term to compensate for imperfections in ALMO forces, we varied ϵSCF between tight 10−6 a.u. and loose 10−2 a.u.

Even with ϵSCF = 10−2, the simulation is stable with the correct average temperature and perfect Maxwell-Boltzmann distribution [Fig. 1(b)]. Δ was initially estimated at 2 × 10−5 fs−1 using Eq. (9) and then refined heuristically to 6 × 10−5 fs−1. We found that it is easier to optimize Δ when γ is set to zero because of reduced noise in the trial runs. Analysis of δfi(t) shows that the error indeed resembles Gaussian white noise. The mean of the error is zero [black line in Fig. 2(a)]. Its ACF decays rapidly [Fig. 2(c)] so that the errors can be considered uncorrelated on time scale of 50 fs. Thus the main assumption behind our approach to ALMO AIMD is justified for liquid water. We established that the main source of error in forces for this system is the loose convergence criterion and not the finite Rc: fully converged ALMO SCF calculations remove the oscillating component of δf [Fig. 2(b)]. We also verified that the ALMO forces converge to the reference forces in the limit Rc (Figure S7 in supplementary material).

FIG. 1.

Calculated properties of water using ALMO AIMD with ϵSCF = 10−2 a.u. and Rc = 1.6 vdWR (red line) and fully converged OT reference (black line). (a) RDF, (b) kinetic energy distribution (the gray curve shows the theoretical Maxwell-Boltzmann distribution), and (c) IR spectrum.

FIG. 1.

Calculated properties of water using ALMO AIMD with ϵSCF = 10−2 a.u. and Rc = 1.6 vdWR (red line) and fully converged OT reference (black line). (a) RDF, (b) kinetic energy distribution (the gray curve shows the theoretical Maxwell-Boltzmann distribution), and (c) IR spectrum.

Close modal
FIG. 2.

The red line is the Euclidean norm of the instantaneous error δfi(t)i, the black line is the magnitude of the time average of the instantaneous error vector, and the green line is the time average of the red line. (a) Rc = 1.6 vdWR and ϵSCF = 10−2 a.u., (b) Rc = 1.6 vdWR and fully converged ALMO SCF, and (c) Normalized ACF 13δfi(t)δfi(t+τ)it of the instantaneous error in panel (a).

FIG. 2.

The red line is the Euclidean norm of the instantaneous error δfi(t)i, the black line is the magnitude of the time average of the instantaneous error vector, and the green line is the time average of the red line. (a) Rc = 1.6 vdWR and ϵSCF = 10−2 a.u., (b) Rc = 1.6 vdWR and fully converged ALMO SCF, and (c) Normalized ACF 13δfi(t)δfi(t+τ)it of the instantaneous error in panel (a).

Close modal

To test the accuracy of ALMO AIMD, we used the trajectory analyzer TRAVIS42 to compute the infrared (IR) spectrum, RDF, and diffusion coefficient of liquid water from both the ALMO trajectory (ϵSCF = 10−2 a.u. and Rc = 1.6 vdWR) and from the reference trajectory. The diffusion coefficients DOT = 1.7(1) × 10−10 and DALMO = 1.8(4) × 10−10 and RDFs [Fig. 1(a)] are in good agreement. The quality of the ALMO IR spectrum [Fig. 1(c)] is good despite minor errors in the intensity of the OH stretching mode, which is sensitive to the precise positions of the centers of localized orbitals. These stringent tests show that despite noticeable errors in the ALMO forces [Fig. 2(a)], the compensating RΔ(t) term in the modified Langevin equation makes it possible to recover atomic dynamics properly. We would like to note that ALMO AIMD could not be stabilized with Δ = 0. Neither were we able to find any values of Δ that stabilize trajectories generated using perturbative versions of ALMO DFT.8 

To demonstrate the computational efficiency of ALMO AIMD, we compared the average wall-time per MD step for a variety of methods in Fig. 3. It is important to emphasize that the comparison is performed for a three-dimensional condensed phase system described with an accurate triple-ζ basis set with polarized functions—a particularly challenging case for DM-based LS methods. ALMO AIMD shows clear LS behavior for all values of ϵSCF, even for medium-size systems. While the second generation Car-Parrinello method decreases the computational overhead of the cubically-scaling AIMD for small systems,28 ALMO AIMD exploits the modified Langevin concept to substantially reduce the simulation cost for systems of all sizes. The crossover point between ALMO AIMD and cubically scaling methods lies in the region of 256 molecules—scale routinely accessible with AIMD today. For comparison, the crossover point between the cubically-scaling DFT and a representative LS method that computes a sparse DM as the matrix sign function of the effective Hamiltonian6 is expected to lie well above 10 000 water molecules [Fig. 5(a) and Fig. S5 of the supplementary material in Ref. 8].

FIG. 3.

Timing benchmarks for PBE/TZV2P simulations of liquid water on 256 compute cores. For ALMO methods, Rc = 1.6 vdWR. Cyan lines represent perfect cubic scaling, whereas gray lines represent perfect linear scaling.

FIG. 3.

Timing benchmarks for PBE/TZV2P simulations of liquid water on 256 compute cores. For ALMO methods, Rc = 1.6 vdWR. Cyan lines represent perfect cubic scaling, whereas gray lines represent perfect linear scaling.

Close modal

Weak scaling benchmarks for very large systems show (Fig. 4) that localized orbitals are naturally suited for parallel execution: LS is retained for a wide range of systems and compute cores. We were able to successfully simulate systems as large as ∼105 atoms within reasonable wall-clock time using only moderate number of compute cores—an impressive feat for AIMD considering that accurate molecular orbitals and the idempotent DM are computed on each step. The horizontal line in Fig. 4 is shown as a rough guide to time and length scales accessible in a fixed wall-clock time given various computational resources. It indicates that ALMO AIMD can extend the range of routine simulations to ∼104 atoms on modern high-performance computing platforms.

FIG. 4.

Weak scalability benchmarks for PBE/TZV2P ALMO AIMD with Rc = 1.6 vdWR and ϵSCF = 10−2 a.u. Dashed gray lines connect systems simulated on the same number of cores to confirm LS behavior.

FIG. 4.

Weak scalability benchmarks for PBE/TZV2P ALMO AIMD with Rc = 1.6 vdWR and ϵSCF = 10−2 a.u. Dashed gray lines connect systems simulated on the same number of cores to confirm LS behavior.

Close modal

To summarize, we demonstrated—for the first time—that compact localized orbitals can be utilized to perform accurate and efficient LS AIMD without concomitant optimization of the DM. High efficiency of the presented method is achieved without sacrificing accuracy with a combination of two techniques: (1) on-the-fly calculation of approximate forces without lengthy self-consistent optimization of localized orbitals and (2) integration of a modified Langevin equation of motion that is fine-tuned to retain stable dynamics even with imperfect forces. By obviating the optimization of the DM, the method remains remarkably efficient even with large localized basis sets. Using liquid water as an example, we showed that the new approach enables simulations of molecular systems on previously inaccessible length scales. The developed method will have a significant impact on modeling of complex molecular systems (e.g., interfaces or nuclei) making completely new phenomena accessible to AIMD. A generalization of the methodology to systems of strongly interacting atoms (e.g., covalent crystals) is underway.

See supplementary material for the detailed description of ALMO AIMD algorithms and dependence of the atomic forces on the electron localization radius.

The authors are grateful to Thomas Kühne for insightful discussions. The research was funded by the Natural Sciences and Engineering Research Council of Canada through the Discovery Grant (No. RGPIN-2016-05059). The authors are grateful to Compute Canada and McGill HPC Centre for computer time.

1.
R.
Car
and
M.
Parrinello
,
Phys. Rev. Lett.
55
,
2471
(
1985
).
2.
S.
Goedecker
,
Rev. Mod. Phys.
71
,
1085
(
1999
).
3.
D. R.
Bowler
and
T.
Miyazaki
,
Rep. Prog. Phys.
75
,
36503
(
2012
).
4.
J.
Kussmann
,
M.
Beer
, and
C.
Ochsenfeld
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
3
,
614
(
2013
).
5.
J.
Aarons
,
M.
Sarwar
,
D.
Thompsett
, and
C.-K.
Skylaris
,
J. Chem. Phys.
145
,
220901
(
2016
).
6.
J.
Vandevondele
,
U.
Borstnik
, and
J.
Hutter
,
J. Chem. Theory Comput.
8
,
3565
(
2012
).
7.
M.
Arita
,
D. R.
Bowler
, and
T.
Miyazaki
,
J. Chem. Theory Comput.
10
,
5419
(
2014
).
8.
R. Z.
Khaliullin
,
J.
VandeVondele
, and
J.
Hutter
,
J. Chem. Theory Comput.
9
,
4421
(
2013
).
9.
C.-K.
Skylaris
,
P. D.
Haynes
,
A. A.
Mostofi
, and
M. C.
Payne
,
J. Chem. Phys.
122
,
084119
(
2005
).
10.
A.
Nakata
,
D. R.
Bowler
, and
T.
Miyazaki
,
Phys. Chem. Chem. Phys.
17
,
31427
(
2015
).
11.
S.
Mohr
,
L. E.
Ratcliff
,
L.
Genovese
,
D.
Caliste
,
P.
Boulanger
,
S.
Goedecker
, and
T.
Deutsch
,
Phys. Chem. Chem. Phys.
17
,
31360
(
2015
).
12.
A. A.
Mostofi
,
P. D.
Haynes
,
C.-K.
Skylaris
, and
M. C.
Payne
,
J. Chem. Phys.
119
,
8842
(
2003
).
13.
H.
van Aggelen
,
Y.
Yang
, and
W.
Yang
,
Phys. Rev. A
88
,
030501
(
2013
).
14.
E.
Tsuchida
,
J. Phys. Soc. Jpn.
76
,
034708
(
2007
).
15.
F.
Mauri
,
G.
Galli
, and
R.
Car
,
Phys. Rev. B
47
,
9973
(
1993
).
16.
P.
Ordejon
,
D. A.
Drabold
,
R. M.
Martin
, and
M. P.
Grumbac
,
Phys. Rev. B
51
,
1456
(
1995
).
17.
J. L.
Fattebert
and
F.
Gygi
,
Comput. Phys. Commun.
162
,
24
(
2004
).
18.
D. R.
Bowler
and
T.
Miyazaki
,
J. Phys.: Condens. Matter
22
,
74207
(
2010
).
19.
T.
Otsuka
,
M.
Taiji
,
D. R.
Bowler
, and
T.
Miyazaki
,
Jpn. J. Appl. Phys., Part 1
55
,
1102B1
(
2016
).
20.
N. D. M.
Hine
,
M.
Robinson
,
P. D.
Haynes
,
C.-K.
Skylaris
,
M. C.
Payne
, and
A. A.
Mostofi
,
Phys. Rev. B
83
,
195102
(
2011
).
21.
H.
Stoll
,
G.
Wagenblast
, and
H.
Preuss
,
Theor. Chim. Acta
57
,
169
(
1980
).
22.
G.
Galli
and
M.
Parrinello
,
Phys. Rev. Lett.
69
,
3547
(
1992
).
23.
L.
Lin
,
J.
Lu
,
L.
Ying
, and
W.
E
,
J. Comput. Phys.
231
,
2140
(
2012
).
24.
R. Z.
Khaliullin
,
M.
Head-Gordon
, and
A. T.
Bell
,
J. Chem. Phys.
124
,
204105
(
2006
).
25.
R. Z.
Khaliullin
,
E. A.
Cobar
,
R. C.
Lochan
,
A. T.
Bell
, and
M.
Head-Gordon
,
J. Phys. Chem. A
111
,
8753
(
2007
).
26.
D. M.
Benoit
,
D.
Sebastiani
, and
M.
Parrinello
,
Phys. Rev. Lett.
87
,
226401
(
2001
).
27.
F. R.
Krajewski
and
M.
Parrinello
,
Phys. Rev. B
73
,
041105
(
2006
).
28.
T. D.
Kühne
,
M.
Krack
,
F. R.
Mohamed
, and
M.
Parrinello
,
Phys. Rev. Lett.
98
,
066401
(
2007
).
29.
J.
Dai
and
J.
Yuan
,
Europhys. Lett.
88
,
20001
(
2009
).
30.
T. D.
Kühne
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
4
,
391
(
2014
).
32.
A. M. N.
Niklasson
,
Phys. Rev. Lett.
100
,
123004
(
2008
).
33.
M. J.
Cawkwell
and
A. M.
Niklasson
,
J. Chem. Phys.
137
,
134105
(
2012
).
34.
See www.cp2k.org for CP2K Open Source Molecular Dynamics, accessed August 18, 2017.
35.
A.
Ricci
and
G.
Ciccotti
,
Mol. Phys.
101
,
1927
(
2003
).
36.
J.
Vandevondele
,
M.
Krack
,
F.
Mohamed
,
M.
Parrinello
,
T.
Chassaing
, and
J.
Hutter
,
Comput. Phys. Commun.
167
,
103
(
2005
).
37.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
38.
S.
Grimme
,
J.
Antony
,
S.
Ehrlich
, and
H.
Krieg
,
J. Chem. Phys.
132
,
154104
(
2010
).
39.
C.
Hartwigsen
,
S.
Goedecker
, and
J.
Hutter
,
Phys. Rev. B
58
,
3641
(
1998
).
40.
J.
Kolafa
,
J. Comput. Chem.
25
,
335
(
2004
).
41.
J.
Vandevondele
and
J.
Hutter
,
J. Chem. Phys.
118
,
4365
(
2003
).
42.
M.
Brehm
and
B.
Kirchner
,
J. Chem. Inf. Model.
51
,
2007
(
2011
).

Supplementary Material