Nuclear Quantum Effects on Aqueous Electron Attachment and Redox

Letter pubs.acs.org/JPCL

Nuclear Quantum Effects on Aqueous Electron Attachment and Redox Properties Vladimir V. Rybkin* and Joost VandeVondele* Nanoscale Simulations, Department of Materials, ETH Zürich, Wolfgang-Pauli-Str. 27, CH-8093 Zürich, Switzerland S Supporting Information *

ABSTRACT: Nuclear quantum effects (NQEs) on the reduction and oxidation properties of small aqueous species (CO2, HO2, and O2) are quantified and rationalized by first-principles molecular dynamics and thermodynamic integration. Vertical electron attachment, or electron affinity, and detachment energies (VEA and VDE) are strongly affected by NQEs, decreasing in absolute value by 0.3 eV going from a classical to a quantum description of the nuclei. The effect is attributed to NQEs that lessen the solvent response upon oxidation/reduction. The reduction of solvent reorganization energy is expected to be general for small solutes in water. In the thermodynamic integral that yields the free energy of oxidation/reduction, these large changes enter with opposite sign, and only a small net effect (0.1 eV) remains. This is not obvious for CO2, where the integrand is strongly influenced by NQEs due to the onset of interaction of the reduced orbital with the conduction band of the liquid during thermodynamic integration. We conclude that NQEs might not have to be included in the computation of redox potentials, unless high accuracy is needed, but are important for VEA and VDE calculations.

N

uclear quantum effects (NQEs) have recently attracted significant attention in the field of molecular modeling due to considerable conceptual and technical progress. Indeed, as computational accuracy increases, ignoring NQEs can be the largest source of error in atomistic simulation. At the same time, efficient methods for the computation of the NQE have become available. The availability of NQEs in molecular dynamics (MD) simulations results in a paradigmatic shift requiring some reconsideration of our understanding of molecular simulations that have been performed for almost half a century in a classical fashion. Although strictly speaking NQEs are not observable because there are no classical nuclei in nature, they are observable via isotopic effects.1,2 Some manifestations of NQEs are easy to foresee, although not easy to quantify, for example, the presence of zero-point vibrational energy (ZPE)3 and proton tunneling.4 Others are less obvious: lowering of the conduction band (0.6 eV in the liquid water5 and also found for ice6) analogous to the established isotope band gap renormalization in semiconductors and insulators,7 competitive NQEs in different vibrational modes,8 destabilization of the first solvation shell on aqueous ions,9 and selective solvation effects on analogous ions (experiments with HOD showed that fluoride preferentially binds to OH and iodine to OD10). Here we present another nonstraightforward manifestation of NQEs, namely, on redox properties in solution, a field that has recently become amendable to first-principles simulation.11−16 In this work, the reduction/oxidation of three aqueous systems, CO2, HO2 radical, and O2, is studied by DFT MD in explicit solvent,16 using quantum and classical thermostats (see Computational Details) © XXXX American Chemical Society

CO2 ⇋ CO2•−

(1)

HO2• ⇋ HO2−

(2)

O2•• ⇋ O2•−

(3)

For all systems, vertical reduction/oxidation is studied, providing access to the experimentally observable vertical electron attachment and detachment energies (VEA and VDE, respectively). For CO2 and HO2 , full thermodynamic integration is performed to obtain free energies of reduction, giving access to redox potentials. For thermodynamic integration, the mixed Hamiltonian has been utilized Hλ = (1 − λ)Hred + λHox

(4)

where Hred is the Hamiltonian of the reduced state, Hox is the Hamiltonian of the oxidized system, and λ is the electron transfer parameter. The thermodynamic integral, which yields the free-energy difference (ΔA) between these states, is calculated according to the following equation ΔA =

∫0

1

∂Hλ ∂λ

dλ =

∫0

1

⟨ΔEλ⟩ dλ

(5)

Mixed Hamiltonians formally correspond to fictitious systems with noninteger oxidation states except for λ = 0 and λ = 1, where Hλ becomes equal to the reduced and the oxidized states, respectively. ⟨E0⟩ is the negative vertical electron detachment Received: February 16, 2017 Accepted: March 15, 2017 Published: March 15, 2017 1424

DOI: 10.1021/acs.jpclett.7b00386 J. Phys. Chem. Lett. 2017, 8, 1424−1428

Letter

The Journal of Physical Chemistry Letters energy (VDE) and ⟨E1⟩ is the vertical electron attachment energy (VEA), both observable quantities. The three different solutes studied have interesting complementary properties. The oxidized solute CO2 is a singlet linear molecule, and the reduced state CO2•− is a doublet radical with a bent geometry. On the contrary, HO2• is a doublet, whereas HO2− is a singlet, without qualitative changes in geometric structure upon oxidation/reduction. O2 is a triplet in the ground state, while O2•− is a doublet. The reduction of HO2• completes the electron pairing on the HOMO (SOMO). The reduction of O2 pairs only one of the uncompensated spins. Conversely, the reduction of CO2 implies an electron populating the LUMO (turning into SOMO). Another key difference lays in the fact that the reduced energy level of CO2 is situated just below the lower edge of the bulk water conduction band, the one of O2 is somewhat lower, and the corresponding level of HO2 lies in the middle of the band gap. Thus one would anticipate the manifestation of NQEs in the following ways: •Thermally averaged geometric and electronic properties of the free solutes. Naively, HO2 is expected to be more effected by NQEs because it contains a light hydrogen atom. •Changes in the solvent shell structure and associated solvent rearrangement energy. •Effects due to the band gap renormalization (lowering of the water conduction band edge5). If the reduced state energy level is sufficiently close to the band gap edge, this state can delocalize due to hybridization with the solute orbitals. On the basis of the available data16 this could be expected for CO2 and possibly O2 but is not supposed to effect the HO2 solution. The thermodynamic integrands are shown in Figures 1 and 2 and summarized in Tables 1 and 2. λ = 0 corresponds to the reduced (negatively charged) state and λ = 1 to the oxidized (neutral) state. The striking feature of the plots is that despite the apparent differences between the systems and shape of the plots, NQEs on the thermodynamic integral (eq 5) partially compensate one another for all systems. Indeed, NQEs reduce the absolute values of VEA and VDE by up to 0.3 eV, which partially cancel each other in the integration to ca. 0.1 eV in the thermodynamic integral. To validate these observed trends, VEA and VDE of CO2 and HO2 are also computed using more accurate hybrid DFT calculations, as listed in Table 3. Furthermore, the accuracy of the quantum MD setup based on quantum thermostats was checked by a path-integral MD simulation, PIMD (CO2, λ = 1, see Table 2). Both hybrid DFT calculations and PIMD simulations result in very similar numbers for NQEs, confirming the validity and reliability of the main data sets. Thus NQEs on VEA and VDE are distinguishable within the MD statistical errors; the procedure used to obtain the latter is given in the Supporting Information. The integrand for HO2 exhibits a behavior as expected for a system in a linear Marcus regime,17 but the one for CO2 is distinctly different in that there is a nearly discontinuous jump in the curve. Furthermore, the location of this feature shifts clearly upon the introduction of NQEs, which we explain in the following. The nonlinearity of the integrand goes hand in hand with a linear-to-bend geometry change of the CO2 molecule upon reduction. This geometry change modifies the electron structure of the solute and triggers a large solvent rearrangement,15,16 strongly influencing ΔE. In the case of quantum MD sampling this transition is less sharp and displaced to the “more reduced” state (lower λ values). As shown in the lower panel of Figure 1, this effect is not visible in the gas phase, and as such

Figure 1. Thermodynamic integration for CO2: top, aqueous; bottom, gas phase. λ = 0 corresponding to the reduced (negatively charged) state and λ = 1 to the oxidized (neutral) state. ⟨ΔE⟩λ is the average energy difference between the oxidized and the reduced states at the given value of λ, in electronvolts.The inset in the top panel shows the thermally averaged CO2 bending angle in degrees.

Figure 2. Thermodynamic integration for HO2: top, aqueous; bottom, gas phase. λ = 0 corresponding to the reduced (negatively charged) state and λ = 1 to the oxidized (neutral) state. ⟨ΔE⟩λ is the average energy difference between the oxidized and the reduced states at the given value of λ, in electronvolts.

1425

DOI: 10.1021/acs.jpclett.7b00386 J. Phys. Chem. Lett. 2017, 8, 1424−1428

Letter

The Journal of Physical Chemistry Letters Table 1. Thermodynamic Data: Thermodynamic Integrals ΔA (TI) in eVa

Ereorg =

TI CO2 HO2 a

classical

quantum

ΔNQE

−1.466 ± 0.005 0.948 ± 0.035

−1.364 ± 0.007 1.057 ± 0.013

−0.102 −0.109

1 (VDE − VEA) 2

(6)

For all three systems it is ca. 1.8 eV in the classical case and 1.5 eV in the quantum case. The observed NQEs in VDE and VEA are thus attributed to a change of solvent reorganization energy, in particular, a reduced solvent reorganization energy in the presence of quantum nuclei. From noncontinuum statistical theory of Matyushov,24 MD simulations,25 and experimental measurements of the model systems in polar solvents,26 it is known that the solvent reorganization energy decreases with the temperature. We can connect to this insight by realizing that the key difference between classical and quantum thermostatting is that the latter preserves the vibrational ZPE, raising, by a few 100 K, the effective temperature computed from classical kinetic energy equipartitioning. This is also consistent with the electron transfer theories in the harmonic bath, which, for instance, predict for the Ru(II)/(III) self-exchange reaction the reduction of reorganization energy due to quantization of nuclei to be ∼0.2 eV,27 a similar order of magnitude as reported here. To conclude, for three aqueous redox pairs, CO2/CO2•−, HO2•/HO2−, and O2••/O2•−, we have found that NQEs have a significant effect (up to 0.33 eV) on the computed observables, VEA and VDE. These differences may likely be measured (e.g., by photoemission spectroscopy28,29) in experiments with normal and heavy water. However, the overall effect on the thermodynamic integrals (and consequently, on the redox potential differences) is very modest (on the order of 0.1 eV). This is explained by the origin of the NQEs, namely, the reduced solvent reorganization energy, which effectively cancels in the thermodynamic integral. Our data show that the effect is independent of the nature of the solute: composition, geometric structure, spin states, and energy level alignment. Therefore, we expect it to be a general trend for aqueous smallmolecule neutral/anion redox pairs.

ΔNQE is the difference between classical and quantum TI.

we rule out the “shrinkage effect”, well-known from gas electron diffraction,18,19 that the most likely geometry of a linear molecule is actually bend, due to thermal and zero-point vibrations. (For a discussion on gas-phase CS2, see ref 20.) Therefore, this effect is related to the solvent and, in particular, the interaction of the solute levels with the conduction band of the liquid,21 which leads to a delocalized anionic state upon reduction of neutral CO2. The quantum shift of the infliction point is than rationalized as the result of a conduction band edge lowering due to NQEs,5 which implies that the delocalized state, corresponding to the linear CO2 geometry, will persist to lower values of λ as compared with the classical sampling. Only for solutes with precisely tuned vertical levels, just above the conduction band edge in H2O but below in D2O, the NQEs of the solvent could make the difference between having a localized and delocalized reduced state, with a large effect on the redox potential. The quantum effect on VEA and VDE remains to be explained. First, we exclude NQEs on the electronic band structure to be the decisive factor. One indication for this is that the both the hybrid and the GGA calculations exhibit very similar shifts, despite the differences in band gap and electron delocalization tendencies of the two approaches. Noticeable delocalization is observed in the case of vertical reduction of CO2, as indicated by the Mulliken spin population on CO2 for classical and quantum sampling equal to 0.076 ± 0.002 and 0.071 ± 0.003, respectively. This is in agreement with recent spectroscopic observations on the hydrated CO2 anionradical.22 Such a difference is, however, too small to explain a noticeable change of the VEA. Moreover, a similar shift in VEA is observed for HO2 and O2 with localized anionic states for all values of λ. Second, NQEs of the gas-phase solute are very modest (bottom panels of Figures 1 and 2) and not sufficient to explain the large shifts observed in solution. We therefore suggest that there are NQEs on the structural solvent response, common to all solutes. On the basis of the linear response Marcus theory framework, we can approximate the effective solvent reorganization energy from the vertical quantities as23

■

COMPUTATIONAL DETAILS All MD simulations have been performed with the CP2K program30,31 at the GGA-DFT level of theory using the BLYP functional,32,33 including the second-generation, D2,34 dispersion corrections of Grimme. Goedecker−Teter−Hutter pseudopotentials35 together with triple-ζ quality basis sets TZV2P36 have been used. The Gaussian-and-Plane-Waves method37 employed an auxiliary plane-wave basis with a 400 Ry cutoff. The SCF guess was propagated using the always stable predictor-corrector method,38 with the SCF convergence criterion being set to 10−5 a.u.

Table 2. Thermodynamic Data from GGA DFT: Vertical Electron Affinity (VEA) and Vertical Detachment Energy (VDE) in eVa VEA classical CO2

−2.217 ± 0.018

HO2 O2

−1.185 ± 0.067 −1.492 ± 0.036

VDE

quantum −1.892 −1.952 −0.875 −1.344

± ± ± ±

0.011 0.007 0.076 0.015

ΔNQE

classical

quantum

ΔNQE

−0.325 −0.265 −0.310 −0.148

1.481 ± 0.054

1.305 ± 0.074

0.176

2.364 ± 0.011 2.104 ± 0.009

2.190 ± 0.022 1.810 ± 0.028

0.174 0.294

a ΔNQE is the difference between classical and quantum values. The upper number for CO2 quantum-computed VEA corresponds to the quantum CNGLE thermostat, while the lower one corresponds to the more accurate PIMD using PIGLET.

1426

DOI: 10.1021/acs.jpclett.7b00386 J. Phys. Chem. Lett. 2017, 8, 1424−1428

Letter

The Journal of Physical Chemistry Letters

Table 3. Thermodynamic Data from Hybrid DFT: Vertical Electron Affinity (VEA) and Vertical Detachment Energy (VDE) in eVa VEA CO2 HO2 a

VDE

classical

quantum

ΔNQE

classical

quantum

ΔNQE

−3.041 ± 0.022 −0.941 ± 0.103

−2.700 ± 0.019 −0.763 ± 0.117

−0.341 −0.178

2.512 ± 0.074 3.703 ± 0.029

2.216 ± 0.259 3.577 ± 0.043

0.296 0.126

ΔNQE is the difference between classical and quantum values.

Additional hybrid DFT calculations with PBE0 functional39 were performed for CO2 and HO2 using 400 to 500 frames from GGA-DFT-driven MD at λ = 0 and 1. The simulations were performed in cubic periodic cells containing 63 water molecules and the solute. The edge size of 12.22 Å corresponds to the bulk water density at the applied level of theory. Each MD trajectory has been integrated for at least 30 000 steps of 0.5 fs after equilibration, giving a total length of 15 ps. Function values at the integration points for the thermodynamic integrals have been computed in parallel. NQEs were simulated by the quantum colored-noise generalized Langevin equation thermostat (CNGLE),40 whereas the canonical CNGLE41 thermostat was used for the classical calculations. PIMD simulations also have been performed using the accelerated PIGLET42 algorithm with eight beads for the VEA of CO2 (λ = 1) to check the accuracy of the quantum thermostat: The two numbers are in a good agreement.

■

(3) Tikhonov, D. S.; Otlyotov, A. A.; Rybkin, V. V. The effect of molecular dynamics sampling on the calculated observable gas-phase structures. Phys. Chem. Chem. Phys. 2016, 18, 18237−18245. (4) Marx, D.; Tuckerman, M. E.; Hutter, J.; Parrinello, M. The nature of the hydrated excess proton in water. Nature 1999, 397, 601−604. (5) Del Ben, M.; Hutter, J.; VandeVondele, J. Probing the structural and dynamical properties of liquid water with models including nonlocal electron correlation. J. Chem. Phys. 2015, 143, 054506. (6) Engel, E. A.; Monserrat, B.; Needs, R. J. Vibrational renormalisation of the electronic band gap in hexagonal and cubic ice. J. Chem. Phys. 2015, 143, 244708. (7) Zollner, S.; Cardona, M.; Gopalan, S. Isotope and temperature shifts of direct and indirect band gaps in diamond-type semiconductors. Phys. Rev. B: Condens. Matter Mater. Phys. 1992, 45, 3376−3385. (8) Kapil, V.; VandeVondele, J.; Ceriotti, M. Accurate molecular dynamics and nuclear quantum effects at low cost by multiple steps in real and imaginary time: Using density functional theory to accelerate wavefunction methods. J. Chem. Phys. 2016, 144, 054111. (9) Wilkins, D. M.; Manolopoulos, D. E.; Dang, L. X. Nuclear quantum effects in water exchange around lithium and fluoride ions. J. Chem. Phys. 2015, 142, 064509. (10) Diken, E. G.; Shin, J.-W.; Price, E. A.; Johnson, M. A. Isotopic fractionation and zero-point effects in anionic H-bonded complexes: a comparison of the I−• HDO and F−• HDO ion−molecule clusters. Chem. Phys. Lett. 2004, 387, 17−22. (11) Tazhigulov, R. N.; Bravaya, K. B. Free Energies of Redox HalfReactions from First-Principles Calculations. J. Phys. Chem. Lett. 2016, 7, 2490−2495. (12) Marenich, A. V.; Ho, J.; Coote, M. L.; Cramer, C. J.; Truhlar, D. G. Computational electrochemistry: prediction of liquid-phase reduction potentials. Phys. Chem. Chem. Phys. 2014, 16, 15068−15106. (13) Guerard, J. J.; Tentscher, P. R.; Seijo, M.; Samuel Arey, J. Explicit solvent simulations of the aqueous oxidation potential and reorganization energy for neutral molecules: gas phase, linear solvent response, and non-linear response contributions. Phys. Chem. Chem. Phys. 2015, 17, 14811−14826. (14) Ambrosio, F.; Miceli, G.; Pasquarello, A. Redox levels in aqueous solution: Effect of van der Waals interactions and hybrid functionals. J. Chem. Phys. 2015, 143, 244508. (15) Cheng, J.; VandeVondele, J. Calculation of Electrochemical Energy Levels in Water Using the Random Phase Approximation and a Double Hybrid Functional. Phys. Rev. Lett. 2016, 116, 086402. (16) Cheng, J.; Liu, X.; VandeVondele, J.; Sulpizi, M.; Sprik, M. Redox Potentials and Acidity Constants from Density Functional Theory Based Molecular Dynamics. Acc. Chem. Res. 2014, 47, 3522− 3529. (17) Marcus, R. A. On the Theory of Oxidation-Reduction Reactions Involving Electron Transfer. I. J. Chem. Phys. 1956, 24, 966−978. (18) Bastiansen, O.; Trætteberg, M. The influence of thermal motion on structure determination of linear molecules using the electrondiffraction method. Acta Crystallogr. 1960, 13, 1108. (19) Morino, Y.; Nakamura, J.; Moore, P. W. Shrinkage Effect by Thermal Vibrations in Linear-Skeleton Molecules. J. Chem. Phys. 1962, 36, 1050−1056. (20) Tikhonov, D. S.; Sharapa, D. I.; Schwabedissen, J.; Rybkin, V. V. Application of classical simulations for the computation of vibrational properties of free molecules. Phys. Chem. Chem. Phys. 2016, 18, 28325−28338.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpclett.7b00386. The methodology for MD error calculation with an example. (PDF)

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]; [email protected]. *E-mail: [email protected]. ORCID

Vladimir V. Rybkin: 0000-0001-5136-6035 Joost VandeVondele: 0000-0002-0902-5111 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS J.V. acknowledges financial support by the European Union FP7 in the form of an ERC Starting Grant under Contract No. 277910. This research was also supported by NCCR MARVEL funded by the Swiss National Science Foundation. Calculations were enabled by the Swiss National Supercomputer Centre (CSCS) under Project ID ch5 and mr15.

■

REFERENCES

(1) Ceriotti, M.; Fang, W.; Kusalik, P. G.; McKenzie, R. H.; Michaelides, A.; Morales, M. A.; Markland, T. E. Nuclear Quantum Effects in Water and Aqueous Systems: Experiment, Theory, and Current Challenges. Chem. Rev. 2016, 116, 7529−7550. (2) Gonzalez-Lafont, A.; Lluch, J. M. Kinetic isotope effects in chemical and biochemical reactions: physical basis and theoretical methods of calculation. WIREs: Comput. Mol. Sci. 2016, 6, 584−603. 1427

DOI: 10.1021/acs.jpclett.7b00386 J. Phys. Chem. Lett. 2017, 8, 1424−1428

Letter

The Journal of Physical Chemistry Letters (21) Adriaanse, C.; Cheng, J.; Chau, V.; Sulpizi, M.; VandeVondele, J.; Sprik, M. Aqueous Redox Chemistry and the Electronic Band Structure of Liquid Water. J. Phys. Chem. Lett. 2012, 3, 3411−3415. (22) Janik, I.; Tripathi, G. N. R. The nature of the CO2- radical anion in water. J. Chem. Phys. 2016, 144, 154307. (23) VandeVondele, J.; Ayala, R.; Sulpizi, M.; Sprik, M. Redox free energies and one-electron energy levels in density functional theory based ab initio molecular dynamics. J. Electroanal. Chem. 2007, 607, 113−120. (24) Matyushov, D. V. Reorganization energy of electron transfer in polar liquids. Dependence on reactant size, temperature and pressure. Chem. Phys. 1993, 174, 199−218. (25) Ghorai, P. K.; Matyushov, D. V. Solvent Reorganization Entropy of Electron Transfer in Polar Solvents. J. Phys. Chem. A 2006, 110, 8857−8863. (26) Derr, D. L.; Elliott, C. M. Temperature Dependence of the Outer-Sphere Reorganization Energy. J. Phys. Chem. A 1999, 103, 7888−7893. (27) Blumberger, J.; Lamoureux, G. Reorganization free energies and quantum corrections for a model electron self-exchange reaction: comparison of polarizable and non-polarizable solvent models. Mol. Phys. 2008, 106, 1597−1611. (28) Tentscher, P. R.; Seidel, R.; Winter, B.; Guerard, J. J.; Arey, J. S. Exploring the Aqueous Vertical Ionization of Organic Molecules by Molecular Simulation and Liquid Microjet Photoelectron Spectroscopy. J. Phys. Chem. B 2015, 119, 238−256. (29) Moens, J.; Seidel, R.; Geerlings, P.; Faubel, M.; Winter, B.; Blumberger, J. Energy Levels and Redox Properties of Aqueous Mn2+/3+ from Photoemission Spectroscopy and Density Functional Molecular Dynamics Simulation. J. Phys. Chem. B 2010, 114, 9173− 9182. (30) The CP2K developers group, CP2K is freely available from: http://www.cp2k.org/, 2017. (31) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. Comput. Phys. Commun. 2005, 167, 103−128. (32) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−3100. (33) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (34) Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27, 1787−1799. (35) Goedecker, S.; Teter, M.; Hutter, J. Separable dual-space Gaussian pseudopotentials. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 1703−1710. (36) VandeVondele, J.; Hutter, J. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. J. Chem. Phys. 2007, 127, 114105. (37) Lippert, G.; Hutter, J.; Parrinello, M. A hybrid Gaussian and plane wave density functional scheme. Mol. Phys. 1997, 92, 477−488. (38) Kolafa, J. Time-reversible always stable predictor-corrector method for molecular dynamics of polarizable molecules. J. Comput. Chem. 2004, 25, 335−342. (39) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. J. Chem. Phys. 1999, 110, 6158−6170. (40) Ceriotti, M.; Bussi, G.; Parrinello, M. Nuclear Quantum Effects in Solids Using a Colored-Noise Thermostat. Phys. Rev. Lett. 2009, 103, 030603. (41) Ceriotti, M.; Bussi, G.; Parrinello, M. Langevin Equation with Colored Noise for Constant-Temperature Molecular Dynamics Simulations. Phys. Rev. Lett. 2009, 102, 020601. (42) Ceriotti, M.; Manolopoulos, D. E. Efficient First-Principles Calculation of the Quantum Kinetic Energy and Momentum Distribution of Nuclei. Phys. Rev. Lett. 2012, 109, 100604. 1428

DOI: 10.1021/acs.jpclett.7b00386 J. Phys. Chem. Lett. 2017, 8, 1424−1428