Computational Amperometry of Nanoscale Capacitors in Molecular Simulations.

In recent years, constant applied potential molecular dynamics has allowed researchers to study the structure and dynamics of the electrochemical double-layer of a large variety of nanoscale capacitors. Nevertheless, it has remained impossible to simulate polarized electrodes at fixed total charge. Here, we show that combining a constant potential electrode with a finite electric displacement fills this gap by allowing us to simulate open-circuit conditions. The method can be extended by applying an electric displacement ramp to perform computational amperometry experiments at different current intensities. As in experiments, the full capacitance of the system is obtained at low intensity, but this quantity decreases when the applied ramp becomes too fast with respect to the microscopic dynamics of the liquid.

Electric-field-based Poisson-Boltzmann theory: Treating mobile charge as polarization.

Mobile charge in an electrolytic solution can in principle be represented as the divergence of ionic polarization. After adding explicit solvent polarization a finite volume of an electrolyte can then be treated as a composite nonuniform dielectric body. Writing the electrostatic interactions as an integral over electric-field energy density we show that the Poisson-Boltzmann functional in this formulation is convex and can be used to derive the equilibrium equations for electric potential and ion concentration by a variational procedure developed by Ericksen for dielectric continua [J. L. Ericksen, Arch. Rational Mech. Anal. 183, 299 (2007)AVRMAW0003-952710.1007/s00205-006-0042-4]. The Maxwell field equations are enforced by extending the set of variational parameters by a vector potential representing the dielectric displacement which is fully transverse in a dielectric system without embedded external charge. The electric-field energy density in this representation is a function of the vector potential and the sum of ionic and solvent polarization making the mutual screening explicit. Transverse polarization is accounted for by construction, lifting the restriction to longitudinal polarization inherent in the electrostatic potential based formulation of Poisson-Boltzmann mean field theory.

Chemomechanical equilibrium at the interface between a simple elastic solid and its liquid phase.

Applying diffusion coupled deformation theory, we investigate how the elastic properties of a solid body are modified when forced to keep its chemical potential aligned with that of its melt. The theory is implemented at the classical level of continuum mechanics, treating materials as simple continua defined by uniform constitutive relations. A phase boundary is a sharp dividing surface separating two continua in mechanical and chemical equilibrium. We closely follow the continuum theory of the swelling of elastomers (gels) but now applied to a simple two phase one-component system. The liquid is modeled by a local free energy density defining a chemical potential and hydrostatic pressure as usual. The model is extended to a solid by adding a non-linear shear elastic energy term with an effective modulus depending on density. Imposing chemomechanical equilibrium with the liquid reservoir reduces the bulk modulus of the solid to zero. The shear modulus remains finite. The stability of the hyper-compressible solid is investigated in a thought experiment. A mechanical load is applied to a rectangular bar under the constraint of fixed lateral dimensions. The linear elastic modulus for axial loading is evaluated and found to be larger than zero, implying that the bar, despite the zero bulk modulus, can support a weight placed on its upper surface. The weight is stabilized by the induced shear stress. The density dependence of the shear modulus is found to be a second order effect reducing the density of the stressed solid (chemostriction).

Computing Surface Acidity Constants of Proton Hopping Groups from Density Functional Theory-Based Molecular Dynamics: Application to the SnO2(110)/H2O Interface.

Proton transfer at metal oxide/water interfaces plays an important role in electrochemistry, geochemistry, and environmental science. The key thermodynamic quantity to characterize this process is the surface acidity constant. An ab initio method that combines density functional theory-based molecular dynamics (DFTMD) and free energy perturbation theory has been established for computing surface acidity constants. However, it involves a reversible proton insertion procedure in which frequent proton hopping, e.g., for strong bases and some oxide surfaces (e.g., SnO2), can cause instability issues in electronic structure calculation. In the original implementation, harmonic restraining potentials are imposed on all O-H bonds (denoted by the VrH scheme) to prevent proton hopping and thus may not be applicable for systems involving spontaneous proton hopping. In this work, we introduce an improved restraining scheme with a repulsive potential Vrep to compute the surface acidities of systems in which proton hopping is spontaneous and fast. In this Vrep scheme, a Buckingham-type repulsive potential Vrep is applied between the deprotonation site and all other protons in DFTMD simulations. We first verify the Vrep scheme by calculating the pKa values of H2O and aqueous HS- solution (i.e., strong conjugate bases) and then apply it to the SnO2(110)/H2O interface. It is found that the Vrep scheme leads to a prediction of the point of zero charge (PZC) of 4.6, which agrees well with experiment. The intrinsic individual pKa values of the terminal five-coordinated Sn site (Sn5cOH2) and bridge oxygen site (Sn2ObrH+) are 4.4 and 4.7, respectively, both being almost the same as the PZC. The similarity of the two pKa values indicates that dissociation of terminal water has almost zero free energy at this proton hopping interface (i.e., partial water dissociation), as expected from the acid-base equilibrium on SnO2.

Electromechanics of the liquid water vapour interface.

Two collective properties distinguishing the thin liquid water vapour interface from the bulk liquid are the anisotropy of the pressure tensor giving rise to surface tension and the orientational alignment of the molecules leading to a finite dipolar surface potential. Both properties can be regarded as capillary phenomena and are likely to be coupled. We have investigated this coupling by determining the response of the tangential component of the surface tension to the application of an electric field normal to the surface using finite field molecular dynamics simulations. We find an upside down parabola with a maximum shifted away from zero field. Comparing the molecular dynamics results to a phenomenological electromechanical model we relate the zero field derivative of the tangential part of the surface tension to the electrostatic potential generated by the spontaneous interface polarization. When interpreted with this model our simulations also indicate that Kelvin forces due to electric field gradients at a polarized interface play an important role in the effective dielectric response.

Band positions of anatase (001) and (101) surfaces in contact with water from density functional theory.

Titanium dioxide in the anatase configuration plays an increasingly important role in photo(electro)catalytic applications due to its superior electronic properties when compared to rutile. In aqueous environments, the surface chemistry and energetic band positions upon contact with water determine charge-transfer processes over solid-solid or solid-electrolyte interfaces. Here, we study the interaction of anatase (001) and (101) surfaces with water and the resulting energetic alignment by means of hybrid density functional theory. While the alignment of band positions favors charge-transfer processes between the two facets for the pristine surfaces, we find the magnitude of this underlying driving force to crucially depend on the water coverage and the degree of dissociation. It can be largely alleviated for intermediate water coverages. Surface states and their passivation by dissociatively adsorbed water play an important role here. Our results suggest that anatase band positions can be controlled over a range of almost 1 eV via its surface chemistry.

Simulating Electrochemical Systems by Combining the Finite Field Method with a Constant Potential Electrode.

A better understanding of interfacial mechanisms is needed to improve the performances of electrochemical devices. Yet, simulating an electrode surface at fixed electrolyte composition remains a challenge. Here, we apply a finite electric field to a single electrode held at constant potential and in contact with an aqueous ionic solution, using classical molecular dynamics. The polarization yields two electrochemical interfaces on opposite sides of the same metal slab. The net charge on one electrode surface is the opposite of the net charge on the other, maintaining overall charge neutrality of the metal. The electrode surface charge fluctuations are compensated by the adsorption of ions from the electrolyte, forming a pair of electric double layers with aligned dipoles. This opens the way towards the efficient simulation of electrochemical interfaces using any flavor of molecular dynamics, from classical to first-principles-based methods.

Coupling of Surface Chemistry and Electric Double Layer at TiO2 Electrochemical Interfaces.

Surfaces of metal oxides at working conditions are usually electrified because of the acid-base chemistry. The charged interface compensated with counterions forms the so-called electric double layer. The coupling of surface chemistry and the electric double layer is considered to be crucial but is poorly understood because of the lack of information at the atomistic scale. Here, we used the latest development in density functional theory-based finite-field molecular dynamics simulation to investigate the pH dependence of the Helmholtz capacitance at electrified rutile TiO2(110)-NaCl electrolyte interfaces. It is found that, because of competing forces from surface adsorption and from the electric double layer, water molecules have a stronger structural fluctuation at high pH, and this leads to a much larger capacitance. It is also seen that interfacial proton transfers at low pH increase significantly the capacitance value. These findings elucidate the microscopic origin of the same trend observed in titration experiments.

Finite electric displacement simulations of polar ionic solid-electrolyte interfaces: Application to NaCl(111)/aqueous NaCl solution.

Tasker type III polar terminations of ionic crystals carry a net surface charge as well as a dipole moment and are fundamentally unstable. In contact with electrolytes, such polar surfaces can be stabilized by adsorption of counterions from the solution to form electric double layers. In a previous work [T. Sayer et al., J. Chem. Phys 147, 104702 (2017)], we reported on a classical force field based molecular dynamics study of a prototype model system, namely, a NaCl(111) slab interfaced with an aqueous NaCl solution on both sides. A serious hurdle in the simulation is that the finite width of the slab admits an electric field in the solid perturbing the theoretical charge balance at the interface of semi-infinite systems [half the surface charge density for NaCl(111)]. It was demonstrated that the application of a finite macroscopic field E canceling the internal electric field can recover the correct charge compensation at the interface. In the present work, we expand this method by applying a conjugate electric displacement field D. The benefits of using D instead of E as the control variable are two fold: it does not only speed up the convergence of the polarization in the simulation but also leads to a succinct expression for the biasing displacement field involving only structural parameters which are known in advance. This makes it feasible to study the charge compensating phenomenon of this prototype system with density functional theory based molecular dynamics, as shown in this work.

Effects of third-order susceptibility in sum frequency generation spectra: a molecular dynamics study in liquid water.

When sum frequency generation (SFG) spectroscopy is applied to charged solid/liquid interfaces, the observed SFG signals include both the second-order and third-order polarizations. The latter is called the χ(3) effect, which mainly includes induced molecular orientation by electric fields at charged interfaces. We theoretically evaluate the χ(3) effect on the SFG spectroscopy of liquid water using molecular dynamics (MD) simulations. The MD simulations enable us to definitely calculate the χ(3) susceptibility as a bulk property, and thereby separating it from the usual χ(2) effect shown in the SFG spectra. The calculated results of χ(3) for liquid water are fairly consistent with the experimental estimates. The present finding is utilized to analyze the spectral change of SFG at the air/water interface under electric fields and at the charged silica/water interface. The present analysis of the spectral changes allows for distinguishing the intrinsic change of the interface structure and the χ(3) effect from bulk liquid.

Charge compensation at the interface between the polar NaCl(111) surface and a NaCl aqueous solution.

Periodic supercell models of electric double layers formed at the interface between a charged surface and an electrolyte are subject to serious finite size errors and require certain adjustments in the treatment of the long-range electrostatic interactions. In a previous publication Zhang and Sprik [Phys. Rev. B 94, 245309 (2016)], we have shown how this can be achieved using finite field methods. The test system was the familiar simple point charge model of a NaCl aqueous solution confined between two oppositely charged walls. Here this method is extended to the interface between the (111) polar surface of a NaCl crystal and a high concentration NaCl aqueous solution. The crystal is kept completely rigid and the compensating charge screening the polarization can only be provided by the electrolyte. We verify that the excess electrolyte ionic charge at the interface conforms to the Tasker 1/2 rule for compensating charge in the theory of polar rock salt (111) surfaces. The interface can be viewed as an electric double layer with a net charge. We define a generalized Helmholtz capacitance CH which can be computed by varying the applied electric field. We find CH=8.23 µF cm-2, which should be compared to the 4.23 µF cm-2 for the (100) non-polar surface of the same NaCl crystal. This is rationalized by the observation that compensating ions shed their first solvation shell adsorbing as contact ions pairs on the polar surface.

Computing the Kirkwood g-Factor by Combining Constant Maxwell Electric Field and Electric Displacement Simulations: Application to the Dielectric Constant of Liquid Water.

In his classic 1939 paper, Kirkwood linked the macroscopic dielectric constant of polar liquids to the local orientational order as measured by the g-factor (later named after him) and suggested that the corresponding dielectric constant at short-range is effectively equal to the macroscopic value just after "a distance of molecular magnitude" [ Kirkwood, J. Chem. Phys., 1939, 7, 911 ]. Here, we show a simple approach to extract the short-ranged Kirkwood g-factor from molecular dynamics (MD) simulation by superposing the outcomes of constant electric field E and constant electric displacement D simulations [ Zhang and Sprik, Phys. Rev. B: Condens. Matter Mater. Phys., 2016, 93, 144201 ]. Rather than from the notoriously slow fluctuations of the dipole moment of the full MD cell, the dielectric constant can now be estimated from dipole fluctuations at short-range, accelerating the convergence. Exploiting this feature, we computed the bulk dielectric constant of liquid water modeled in the generalized gradient approximation (PBE) to density functional theory and found it to be at least 40% larger than the experimental value.

Density Functional Theory Calculation of the Band Alignment of (101Ì 0) In(x)Ga(1-x)N/Water Interfaces.

Conduction band edge (CBE) and valence band edge (VBE) positions of InxGa1-xN photoelectrodes were computed using density functional theory methods. The band edges of fully solvated GaN and InN model systems were aligned with respect to the standard hydrogen electrode using a molecular dynamics hydrogen electrode scheme applied earlier to TiO2/water interfaces. Similar to the findings for TiO2, we found that the Purdew-Burke-Ernzerhof (PBE) functional gives a VBE potential which is too negative by 1 V. This cathodic bias is largely corrected by application of the Heyd-Scuseria-Ernzerhof (HSE06) hybrid functional containing a fraction of Hartree-Fock exchange. The effect of a change of composition was investigated using simplified model systems consisting of vacuum slabs covered on both sides by one monolayer of H2O. The CBE was found to vary linearly with In content. The VBE, in comparison, is much less sensitive to composition. The data show that the band edges straddle the hydrogen and oxygen evolution potentials for In fractions less than 47%. The band gap was found to exceed 2 eV for an In fraction less than 54%.

Electronic Energy Levels and Band Alignment for Aqueous Phenol and Phenolate from First Principles.

Electronic energy levels in phenol and phenolate solutions have been computed using density functional theory and many-body perturbation theory. The valence and conduction bands of the solvent and the ionization energies of the solutes have been aligned with respect to the vacuum level based on the concept of a computational standard hydrogen electrode. We have found significant quantitative differences between the generalized-gradient approximation, calculations with the HSE hybrid functional, and many-body perturbation theory in the G0W0 approximation. For phenol, two ionization energies below the photoionization threshold of bulk water have been assigned in the spectrum of Kohn-Sham eigenvalues of the solution. Deprotonation to phenolate was found to lift a third occupied energy level above the valence band maximum of the solvent which is characterized by an electronic lone pair at the hydroxyl group. The second and third ionization energies of phenolate were found to be very similar and explain the intensity pattern observed in recent experiments using liquid-microjet photoemission spectroscopy.

Molecular simulation study of hydrated Na-rectorite.

The swelling behavior of clay minerals is an important issue in industrial processes and environmental applications. Mixed-layer clay minerals containing a smectite fraction, such as rectorite, are neglected even though they could swell and exist in nature widely. The hydration of rectorite has not been well comprehended even though they are meaningful to mineralogy and industry. This study combines molecular dynamics (MD) and Monte Carlo (MC) simulations to disclose the swelling behavior of rectorite and compare with montmorillonite. From grand canonical Monte Carlo (GCMC) and MD simulations, we obtain swelling curves and swelling free-energy curves of rectorite with a relative humidity of 100%. With the comparisons of swelling free-energy minima, we find that the bilayer hydrate of Na-rectorite is more thermodynamically stable than the monolayer hydrate, which is similar to Na-montmorillonite. However, the interlayer sodium ions in rectorite show an asymmetrical distribution quite different from the symmetrical distribution in montmorillonite. Because of unequal layer charges between the smectite part and illite part of retorite, sodium ions prefer to distribute close to the illite part surface.

Aqueous transition-metal cations as impurities in a wide gap oxide: the Cu(2+)/Cu(+) and Ag(2+)/Ag(+) redox couples revisited.

The interactions of the d electrons of transition-metal aqua ions with the solvent are usually divided in short-range electronic interactions with ligand water molecules and long-range electrostatic interactions with molecules beyond the first coordination shell. This is the rationale behind the cluster continuum and QM/MM methods developed for the computation of the redox potentials. In the density functional theory based molecular dynamics (DFTMD) method, the electronic states of the complex are also allowed to mix with the extended band states of the solvent. Returning to the Cu(+) and Ag(+) oxidation reaction, which has been the subject of DFTMD simulation before, we show that coupling to the valence band states of water is greatly enhanced by the band gap error in the density functional approximation commonly used in DFTMD (the generalized gradient approximation). This effect is analyzed by viewing the solvent as a wide gap oxide and the redox active ions as electronic defects. The errors can be reduced significantly by application of hybrid functionals containing a fraction of Hartree-Fock exchange. These calculations make use of recent progress in DFTMD technology, enabling us to include sp core polarization and Hartree-Fock exchange in condensed-phase model systems.

Redox potentials and acidity constants from density functional theory based molecular dynamics.

CONSPECTUS: All-atom methods treat solute and solvent at the same level of electronic structure theory and statistical mechanics. All-atom computation of acidity constants (pKa) and redox potentials is still a challenge. In this Account, we review such a method combining density functional theory based molecular dynamics (DFTMD) and free energy perturbation (FEP) methods. The key computational tool is a FEP based method for reversible insertion of a proton or electron in a periodic DFTMD model system. The free energy of insertion (work function) is computed by thermodynamic integration of vertical energy gaps obtained from total energy differences. The problem of the loss of a physical reference for ionization energies under periodic boundary conditions is solved by comparing with the proton work function computed for the same supercell. The scheme acts as a computational hydrogen electrode, and the DFTMD redox energies can be directly compared with experimental redox potentials. Consistent with the closed shell nature of acid dissociation, pKa estimates computed using the proton insertion/removal scheme are found to be significantly more accurate than the redox potential calculations. This enables us to separate the DFT error from other sources of uncertainty such as finite system size and sampling errors. Drawing an analogy with charged defects in solids, we trace the error in redox potentials back to underestimation of the energy gap of the extended states of the solvent. Accordingly the improvement in the redox potential as calculated by hybrid functionals is explained as a consequence of the opening up of the bandgap by the Hartree-Fock exchange component in hybrids. Test calculations for a number of small inorganic and organic molecules show that the hybrid functional implementation of our method can reproduce acidity constants with an uncertainty of 1-2 pKa units (0.1 eV). The error for redox potentials is in the order of 0.2 V.

Aligning electronic and protonic energy levels of proton-coupled electron transfer in water oxidation on aqueous TiO2.

The high overpotential in water oxidation on anodes is a limiting factor for the large-scale application of photoelectrochemical cells. To overcome this limitation, it is essential to understand the four proton-coupled electron transfer (PCET) steps in the reaction mechanism and their implications to the overpotential. Herein, a simple scheme to compute the energies of the PCET steps in water oxidation on the aqueous TiO2 surface using a hybrid density functional is described. An energy level diagram for fully decoupled electron- and proton-transfer reactions in which both electronic and protonic levels are placed on the same potential scale is also described. The level diagram helps to visualize the electronic and protonic components of the overpotential, and points out what are needed to improve. For TiO2, it is found that its catalytic activity is due to aligning the protonic energy levels in the PCET steps, while improving the activity requires also aligning the electronic levels.

The ionization potential of aqueous hydroxide computed using many-body perturbation theory.

The ionization potentials of electrolyte solutions provide important information about the electronic structure of liquids and solute-solvent interactions. We analyzed the positions of solute and solvent bands of aqueous hydroxide and the influence of the solvent environment on the ionization potential of hydroxide ions. We used the concept of a computational hydrogen electrode to define absolute band positions with respect to vacuum. We found that many-body perturbation theory in the G0 W0 approximation substantially improves the relative and absolute positions of the band edges of solute and solvent with respect to those obtained within Density Functional Theory, using semi-local functionals, yielding results in satisfactory agreement with recent experiments.

Solute-Solvent Charge-Transfer Excitations and Optical Absorption of Hydrated Hydroxide from Time-Dependent Density-Functional Theory.

The electronic structure of simple hydrated ions represents one of the most challenging problems in electronic-structure theory. Spectroscopic experiments identified the lowest excited state of the solvated hydroxide as a charge-transfer-to-solvent (CTTS) state. In the present work we report computations of the absorption spectrum of the solvated hydroxide ion, treating both solvent and solute strictly at the same level of theory. The average absorption spectrum up to 25 eV has been computed for samples taken from periodic ab initio molecular dynamics simulations. The experimentally observed CTTS state near the onset of the absorption threshold has been analyzed at the generalized-gradient approximation (GGA) and with a hybrid density-functional. Based on results for the lowest excitation energies computed with the HSE hybrid functional and a Davidson diagonalization scheme, the CTTS transition has been found 0.6 eV below the first absorption band of liquid water. The transfer of an electron to the solvent can be assigned to an excitation from the solute 2pπ orbitals, which are subject to a small energetic splitting due to the asymmetric solvent environment, to the significantly delocalized lowest unoccupied orbital of the solvent. The distribution of the centers of the excited state shows that CTTS along the OH(-) axis of the hydroxide ion is avoided. Furthermore, our simulations indicate that the systematic error arising in the calculated spectrum at the GGA originates from a poor description of the valence band energies in the solution.