Constant-μ method

The constant-μ method performs first-principles calculations at a fixed electrode potential, that is, at a fixed electron chemical potential \(\mu_\mathrm{e}\). It treats the number of electrons in the system as a dynamical variable and introduces a fictitious equation of motion through which the system exchanges electrons with an external potentiostat. The electronic structure at each step is still obtained by an ordinary self-consistent calculation at a fixed number of electrons, so the calculation remains stable. The counter charge of the electrode is carried by the metal electrode of the ESM method. Atomic units (\(e=1\)) are used below.

Constant charge versus constant potential

Ordinary first-principles calculations fix the number of electrons \(N\) in the system (constant charge). In electrochemical experiments, by contrast, a potentiostat controls the electrode potential (constant potential). The difference becomes significant for reactions that involve electron transfer.

Consider hydrogen adsorption on a Pt electrode, \(\mathrm{H^+}+e^-\to\mathrm{H_{ad}}\). In a constant-charge molecular dynamics simulation with an excess charge on the electrode, the Fermi level (electrode potential) changes substantially as electrons move from the electrode to the adsorbate. Under experimental conditions, the potentiostat replenishes the electrons and the electrode potential stays fixed. Constant-charge calculations thus fail to match experiment, because the electrode potential shifts with every reaction event.

Grand-potential formulation

Treating the number of electrons \(N\) as a variable, consider the grand potential at an externally imposed electron chemical potential \(\mu_\mathrm{ext}\):

\[\begin{equation} \Omega(\{\boldsymbol{r}\};\mu_\mathrm{ext})=\min_N\left[E(\{\boldsymbol{r}\};N)-\mu_\mathrm{ext}N\right] \label{eq:grand} \end{equation}\]

where \(E(\{\boldsymbol{r}\};N)\) is the DFT total energy for atomic configuration \(\{\boldsymbol{r}\}\) and \(N\) electrons. The stationarity condition

\[\begin{equation} \mu_\mathrm{e}\equiv\frac{\partial E}{\partial N}=\mu_\mathrm{ext} \label{eq:equilibrium} \end{equation}\]

states that the Fermi level of the system equals the chemical potential of the potentiostat. Since \(\partial^2E/\partial N^2=\partial\mu_\mathrm{e}/\partial N\) is the inverse of the differential capacitance \(C\) of the system and is positive, the stationary point is a minimum. The electrode potential is defined as \(\Phi=-\mu_\mathrm{e}/e\).

\(\Omega\) is the Legendre transform of \(E(N)\) with respect to \(N\). Obtaining \(\Omega\) from constant-charge calculations alone requires calculations at several values of \(N\), interpolation, and a manual Legendre transform.

Adjusting the electron number inside the SCF loop and its problems

Lozovoi et al. (2001) proposed adjusting the number of electrons within the self-consistent field (SCF) calculation. An ordinary SCF calculation fills the Kohn-Sham levels from the bottom until the electron count reaches \(N\) in each iteration. In Lozovoi's method, levels are filled until the Fermi level reaches \(\mu\), so the number of electrons changes from iteration to iteration.

This approach has the following problems.

• The variational principle does not hold for the SCF iterations, so convergence is not guaranteed. While the total energy of an ordinary SCF calculation converges steadily, this method often diverges no matter how small the density mixing is made.
• The initial charge density must be chosen carefully, which makes the method unsuitable for applications such as molecular dynamics that repeat many SCF calculations.

The SCF calculation itself should therefore be done at constant charge, and the number of electrons should be adjusted outside the SCF loop so that the Fermi level approaches its target. The extended-system approach of molecular dynamics serves this purpose.

Analogy with the constant-pressure method

The constant-μ method has the same structure as the constant-pressure molecular dynamics of Andersen (1980). In the constant-pressure method, the cell volume \(V\) is a dynamical variable, and atomic positions are written with scaled coordinates \(\tilde{\boldsymbol{r}}_i=V^{-1/3}\boldsymbol{r}_i\). The fictitious Lagrangian

\[\begin{equation} L_P=\frac{1}{2}\sum_i m_iV^{2/3}\dot{\tilde{\boldsymbol{r}}}_i^2-E(\{V^{1/3}\tilde{\boldsymbol{r}}\})+\frac{1}{2}W\dot{V}^2-P_\mathrm{ext}V \label{eq:andersen-lagrangian} \end{equation}\]

yields the equation of motion for the volume

\[\begin{equation} W\ddot{V}=\frac{1}{3V}\left[\sum_i m_iV^{2/3}\dot{\tilde{\boldsymbol{r}}}_i^2-\sum_i\boldsymbol{r}_i\cdot\frac{\partial E}{\partial\boldsymbol{r}_i}\right]-P_\mathrm{ext}=P_\mathrm{in}-P_\mathrm{ext} \label{eq:andersen-eom} \end{equation}\]

where \(W\) is a fictitious mass of the cell and \(P_\mathrm{in}\) is the instantaneous internal pressure from the virial theorem. A feedback mechanism operates: the volume grows when the internal pressure exceeds the external one and shrinks otherwise, so it oscillates around its equilibrium value and \(\langle P_\mathrm{in}\rangle=P_\mathrm{ext}\).

The constant-μ method replaces the volume by the number of electrons and the pressure by the chemical potential.

	Constant pressure	Constant μ
Extensive variable	volume \(V\)	electron number \(N\)
Intensive variable	pressure \(P_\mathrm{in}\)	chemical potential \(\mu_\mathrm{in}\)
Imposed value	\(P_\mathrm{ext}=\langle P_\mathrm{in}\rangle\)	\(\mu_\mathrm{ext}=\langle\mu_\mathrm{in}\rangle\)
Fictitious mass	\(W\)	\(M\)
Restoring force near equilibrium	bulk modulus \(B=-V_0\,\partial P/\partial V\)	inverse capacitance \(1/C=\partial\mu/\partial N\)

This analogy had not been applied to the electron number because, under periodic boundary conditions, there was no physically meaningful way to change the charge of the cell (see below).

Equations of motion for the fictitious charge particle

Lagrangian and equations of motion

Connect the system to a potentiostat and make the electron number \(N\) a dynamical variable. Moving one electron into the system lowers the energy of the potentiostat by \(\mu_\mathrm{ext}\), so the potential energy of the extended system is \(\tilde{E}=E(\{\boldsymbol{r}\};N)-\mu_\mathrm{ext}N\). Take the fictitious Lagrangian

\[\begin{equation} L_\mu=\frac{1}{2}\sum_i m_i\dot{\boldsymbol{r}}_i^2+\frac{1}{2}M\dot{N}^2-\left[E(\{\boldsymbol{r}\};N)-\mu_\mathrm{ext}N\right] \label{eq:lagrangian} \end{equation}\]

where \(M\) is a fictitious mass for the charge. The Euler-Lagrange equations are

\[\begin{align} m_i\ddot{\boldsymbol{r}}_i&=-\frac{\partial E(\{\boldsymbol{r}\};N)}{\partial\boldsymbol{r}_i},\label{eq:eom-r}\\ M\ddot{N}&=-\left(\frac{\partial E(\{\boldsymbol{r}\};N)}{\partial N}-\mu_\mathrm{ext}\right)=-\left(\mu_\mathrm{in}-\mu_\mathrm{ext}\right)\label{eq:eom-n} \end{align}\]

where \(\mu_\mathrm{in}=\partial E/\partial N\) is the instantaneous chemical potential (Fermi level) of the electrons. Electrons leave the system when the Fermi level is above the target and enter when it is below. This fictitious degree of freedom representing the charge is called the fictitious charge particle (FCP).

Hamiltonian and conserved quantity

With the conjugate momenta \(\boldsymbol{p}_i=m_i\dot{\boldsymbol{r}}_i\) and \(p_N=M\dot{N}\), the Hamiltonian is

\[\begin{equation} H_\mu=\sum_i\frac{\boldsymbol{p}_i^2}{2m_i}+\frac{p_N^2}{2M}+E(\{\boldsymbol{r}\};N)-\mu_\mathrm{ext}N \label{eq:hamiltonian} \end{equation}\]

and Hamilton's equations are

\[\begin{equation} \dot{\boldsymbol{r}}_i=\frac{\boldsymbol{p}_i}{m_i},\qquad \dot{\boldsymbol{p}}_i=-\frac{\partial E}{\partial\boldsymbol{r}_i},\qquad \dot{N}=\frac{p_N}{M},\qquad \dot{p}_N=-\left(\mu_\mathrm{in}-\mu_\mathrm{ext}\right) \label{eq:canonical} \end{equation}\]

\(H_\mu\) is conserved. In a test that propagated only the FCP for an Al monolayer (four atoms per cell, atoms fixed), the Fermi level oscillated harmonically around the target; without a thermostat, the amplitude reflected the initial difference between the Fermi level and the target. The DFT total energy \(E\) oscillated with the Fermi level, while \(H_\mu\) was conserved.

Formulation of Bonnet et al. and thermostats

Bonnet et al. (2012) write the same equations in terms of the electrode potential. With the electrode potential \(\Phi=-\partial E/\partial N=-\mu_\mathrm{in}\) and the potentiostat potential \(\Phi_\mathrm{pot}=-\mu_\mathrm{ext}\), the force on the FCP is \(F_N=-\partial\tilde{E}/\partial N=\Phi-\Phi_\mathrm{pot}\), and

\[\begin{equation} \dot{N}=\frac{P_N}{M},\qquad\dot{P}_N=\Phi-\Phi_\mathrm{pot} \label{eq:bonnet} \end{equation}\]

This force plays the same role as a real electric force: it drives electrons into the system when \(\Phi>\Phi_\mathrm{pot}\). After equilibration, \(\Phi\) fluctuates around \(\Phi_\mathrm{pot}\) and \(\langle\Phi\rangle=\Phi_\mathrm{pot}\).

With a Nosé-Hoover thermostat on the particles only (scheme I), the dynamics yield the grand-canonical distribution at temperature \(T\) and potential \(\Phi_\mathrm{pot}\), assuming ergodicity. The coupling between the particle and charge degrees of freedom is weak, however, and equilibration takes a long time. The charge is therefore also coupled directly to a thermostat: either the particle thermostat (scheme II) or a separate thermostat at temperature \(T_\mathrm{c}\) (scheme III). With a thermostat variable \(\xi_\mathrm{c}\), its momentum \(P_\mathrm{c}\), and its fictitious mass \(M_\mathrm{c}\), the equations of motion for scheme III are

\[\begin{equation} \dot{N}=\frac{P_N}{M},\qquad \dot{P}_N=\Phi-\Phi_\mathrm{pot}-\frac{P_\mathrm{c}}{M_\mathrm{c}}P_N,\qquad \dot{\xi}_\mathrm{c}=\frac{P_\mathrm{c}}{M_\mathrm{c}},\qquad \dot{P}_\mathrm{c}=\frac{P_N^2}{M}-k_\mathrm{B}T_\mathrm{c} \label{eq:thermostat} \end{equation}\]

For \(T_\mathrm{c}=T\), the grand-canonical distribution is obtained. Choosing \(T_\mathrm{c}<T\) departs from any standard ensemble, but because the coupling is weak, the particles and the charge reach a quasi-equilibrium at temperatures \(T\) and \(T_\mathrm{c}\), respectively, and the fluctuations of the electrode potential are reduced.

Statistical-mechanical properties

The extended system with thermostats realizes the grand-canonical ensemble for the electrons. The partition function is

\[\begin{equation} \Xi_\mu=\iiint d\boldsymbol{r}\,d\boldsymbol{p}\,dN\,\exp\left[-\frac{H-\mu_\mathrm{ext}N}{k_\mathrm{B}T}\right],\qquad H=\sum_i\frac{\boldsymbol{p}_i^2}{2m_i}+E(\{\boldsymbol{r}\};N) \label{eq:partition} \end{equation}\]

With the canonical partition function \(Z_N\) at electron number \(N\), \(\Xi_\mu=\int dN\,Z_N\,e^{\mu_\mathrm{ext}N/k_\mathrm{B}T}\). Differentiating the mean electron number with respect to \(\mu_\mathrm{ext}\) expresses the differential capacitance through the fluctuations of the electron number:

\[\begin{equation} C=\frac{\partial\langle N\rangle}{\partial\mu_\mathrm{ext}}=\frac{\langle N^2\rangle-\langle N\rangle^2}{k_\mathrm{B}T}=\frac{\langle\delta N^2\rangle}{k_\mathrm{B}T} \label{eq:capacitance} \end{equation}\]

This corresponds to obtaining the compressibility from volume fluctuations in the constant-pressure method; in electronics it is known as Johnson-Nyquist noise. The standard deviation of the electrode potential is \(\sigma_\Phi=\sigma_N/C\propto\sqrt{T/C}\). Because the capacitance scales with system size, small systems show large potential fluctuations. An electrode connected to an ideal potentiostat has a non-fluctuating potential, so these fluctuations are an artifact of not treating the potentiostat explicitly.

Linearizing Eq. \(\eqref{eq:eom-n}\) around equilibrium with \(\mu_\mathrm{in}-\mu_\mathrm{ext}\simeq\delta N/C\) gives

\[\begin{equation} M\,\delta\ddot{N}\simeq-\frac{1}{C}\,\delta N,\qquad \omega=\frac{1}{\sqrt{MC}} \label{eq:linear} \end{equation}\]

so the charge oscillates around its mean at angular frequency \(\omega\). The fictitious mass \(M\) is chosen so that this frequency couples appropriately to the atomic motion. For an Al electrode, Bonnet et al. (2012) chose \(M\) so that the charge oscillates at \(300\ \mathrm{cm^{-1}}\), near the top of the Al phonon spectrum.

Combination with ESM and the reference of the electrode potential

The constant-μ method itself is general and can be run under ordinary periodic boundary conditions. Changing the charge of a periodic cell, however, requires a uniform background charge to keep the electrostatic energy finite. The uniform background is a computational device that removes the meaning of the absolute Fermi level; calculations of charged defects in semiconductors, for example, need corrections to restore it. A constant-μ calculation under periodic boundary conditions is therefore hard to interpret physically.

In the ESM method, the Poisson equation is solved with a Green's function, so the whole system stays neutral and the compensating charge induced on the ESM electrode is distributed in a physically meaningful way. The absolute Fermi level retains its meaning, which gives constant-μ calculations a clear interpretation. The electrochemical cell is modeled with bc3 (vacuum/slab/metal), in which the ESM metal on one side of the slab (the working electrode) serves as the counter electrode.

The electrode potential is referenced as follows.

• In bc3 the electrostatic potential of the ESM metal is set to zero, so the Fermi level \(\mu_\mathrm{e}\) is measured from the potential of the counter electrode. The electrode potential \(\Phi=-\mu_\mathrm{e}/e\) is the potential of the working electrode with the counter electrode at 0 V.
• At zero excess charge (the potential of zero charge, pzc), there is no field between the slab and the counter electrode, and the potential in the vacuum gap equals that of the counter electrode (zero). Hence \(\mu_\mathrm{e}=-W\), where \(W\) is the work function of the surface facing the counter electrode, and \(\Phi_\mathrm{pzc}=W/e\).
• The relative electrode potential is defined as \(U=\Phi_\mathrm{pot}-\Phi_\mathrm{pzc}\). For the Al electrode of Bonnet et al. (2012) (36 atoms, counter electrode 10 bohr away), \(\Phi_\mathrm{pzc}=4.31\) V; when \(U\) was switched among 0 V, \(-0.5\) V, and \(+0.5\) V, the electrode potential quickly settled into fluctuations around each new target.

Because the electron number \(N\) is treated as a continuous variable, the energy function \(E(\{\boldsymbol{r}\};N)\) must be differentiable for non-integer electron numbers. With approximate exchange-correlation functionals, self-interaction makes orbital energies depend on their occupations, and the behavior at fractional electron numbers can be poor.

Structure optimization and reaction paths

A molecular dynamics run repeats the following steps at each time step.

1. Perform a constant-charge SCF calculation for the current atomic configuration \(\{\boldsymbol{r}\}\) and electron number \(N\).
2. Obtain the atomic forces \(\boldsymbol{F}_i\) and the Fermi level \(\mu_\mathrm{in}\).
3. Update the atoms and the electron number with the same time step using Eq. \(\eqref{eq:canonical}\) (or Eq. \(\eqref{eq:thermostat}\)).

Because the charge dynamics and the electronic structure calculation are separate, the scheme can be added to an existing first-principles molecular dynamics code with minor changes.

The same idea gives structure optimization at constant potential. The potential energy of the extended system, \(\tilde{E}=E(\{\boldsymbol{r}\};N)-\mu_\mathrm{ext}N\), is minimized with the generalized force

\[\begin{equation} \tilde{\boldsymbol{F}}=\left\{\boldsymbol{F}_i,\,F_\mathrm{FCP}\right\},\qquad \boldsymbol{F}_i=-\frac{\partial\tilde{E}}{\partial\boldsymbol{r}_i}=-\frac{\partial E}{\partial\boldsymbol{r}_i},\qquad F_\mathrm{FCP}=-\frac{\partial\tilde{E}}{\partial N}=-\left(\mu_\mathrm{in}-\mu_\mathrm{ext}\right) \label{eq:generalized-force} \end{equation}\]

The atoms and the FCP (the charge) move together after each electronic structure calculation, so there is no extra loop, and a constant-potential structure costs about the same as an ordinary optimization.

The FCP update needs care. Steepest descent converges slowly near the end. Line minimization with a quadratic model converges well when the atoms are fixed, but when the atoms move at the same time, \(F_\mathrm{FCP}\) can be nearly equal in two consecutive steps, and the quadratic prediction of the next charge breaks down. A fixed upper limit on the charge change cannot adapt when the surface area of the cell changes. The limit is therefore set from a vacuum-capacitor model. With \(d\) the distance between the electrode surface and the ESM and \(S\) the area, the vacuum capacitance is \(C_0=S/4\pi d\) and the capacitance of the system is \(C=\varepsilon_\mathrm{sys}S/4\pi d\) (\(\varepsilon_\mathrm{sys}\ge1\)). The charge change needed to close the current potential gap is \(C|\mu_i-\mu_\mathrm{ext}|\), but \(C\) is not known in advance. Its lower bound

\[\begin{equation} \Delta N^\mathrm{upper}=C_0\left|\mu_i-\mu_\mathrm{ext}\right| \label{eq:dn-upper} \end{equation}\]

is used as the upper limit of the charge change per step. This scheme works well for molecules on metal surfaces and for electrode/solid-electrolyte interfaces.

For reaction paths, the nudged elastic band (NEB) method is applied in \(\tilde{E}\) space (constant-μ+NEB). With constant-charge calculations, the total energy at a given potential requires calculations at several electron numbers followed by a Legendre transform; the constant-μ method gives the value at the desired potential in a single calculation. If the reaction path did not depend on the potential, constant-charge NEB calculations with a manual Legendre transform would suffice. In real systems the path depends on the potential, so the constant-μ method is essential. For hydrogen diffusion on Al(001), the method yields the path along the grand potential \(\Delta\Omega=\Delta E_\mathrm{tot}-\mu_\mathrm{ext}\Delta N\): at constant charge the Fermi level varies along the path, whereas at constant μ the potential stays fixed and the charge varies along the path.

References

• N. Bonnet, T. Morishita, O. Sugino, and M. Otani, "First-Principles Molecular Dynamics at a Constant Electrode Potential," Phys. Rev. Lett. 109, 266101 (2012). doi
• M. Otani and O. Sugino, "First-principles calculations of charged surfaces and interfaces: A plane-wave nonrepeated slab approach," Phys. Rev. B 73, 115407 (2006). doi
• A. Y. Lozovoi, A. Alavi, J. Kohanoff, and R. M. Lynden-Bell, "Ab initio simulation of charged slabs at constant chemical potential," J. Chem. Phys. 115, 1661 (2001). doi
• H. C. Andersen, "Molecular dynamics simulations at constant pressure and/or temperature," J. Chem. Phys. 72, 2384 (1980). doi