Simulation of many-electron systems that exchange matter with the environment

TL;DR

This paper reviews computational methods for simulating many-electron systems that exchange matter with their environment, highlighting current approaches, challenges, and future directions in materials science and chemical physics.

Contribution

It provides an overview of existing quantum and classical techniques for modeling open-boundary electronic systems, emphasizing their differences and potential for advancing the field.

Findings

Diverse computational approaches with varying efficiency and accuracy

Identification of conceptual and pragmatic strengths in current methods

Potential for combining approaches to push simulation frontiers

Abstract

The computational treatment of many-electron systems capable of exchanging {electrons and nuclei} with the environment represents one of the outermost frontiers in simulation methodology. The exchanging process occurs in a large variety of natural and artificially induced phenomena which are of major relevance to several leading fields of academic research and modern technology. In this progress report I will present an overview of problems in current materials science and chemical physics where the corresponding computational approaches require the concept of an electronic system with open boundaries. Quantum and Quantum/Classical computational techniques treat the exchange of electrons with the environment at different computational efficiency, conceptual rigorousness and numerical accuracy. The overall emerging picture shows a rich availability of interesting ideas, some with a higher weight on the pragmatic side, others with higher weight on the conceptual side; possible combinations, in perspective, may push the field much beyond its current frontiers.

Equations

$$I(V)=\frac{2e}{h}{\int }_{-\infty }^{\infty }T(E)[f_L(E)-f_R(E)]dE$$

$$T=tr{\mathbf{\Gamma }}_{\mathbf{L}}{\mathcal{G}}_{\mathbf{V}}{\mathbf{\Gamma }}_{\mathbf{R}}{\mathcal{G}}_{\mathbf{V}}^{†}.$$

$${\mathcal{G}}_{0V}^{-1}(\mathbf{x},{\mathbf{x}}^{{}^{\prime }},E)=\sum _n\frac{{\varphi }_n(\mathbf{x}){\varphi }_n^{\ast }({\mathbf{x}}^{{}^{\prime }})}{E-{ϵ}_n+i\eta /2}$$

$$\frac{d}{dt}\mathcal{P}\rho (t)+i\mathcal{P}L\mathcal{P}\rho (t)+{\int }_0^{t}d{t}^{{}^{\prime }}\mathcal{P}L(I-\mathcal{P})\mathcal{U}(t-t^{{}^{\prime }})(I-\mathcal{P})L\mathcal{P}\rho (t^{{}^{\prime }})=0$$

$${\rho }_R(0)\otimes \left(\frac{d}{dt}+i{L}_{eff}^S\right){\rho }_S(t)=-{\int }_0^{t}d{t}^{{}^{\prime }}\mathcal{K}(t-t^{{}^{\prime }}){\rho }_R(0)\otimes {\rho }_S(t^{{}^{\prime }})$$

$$\left(\frac{d}{dt}+i{L}_{eff}^S\right){\rho }_S(t)=-{\int }_0^{t}d{t}^{{}^{\prime }}{\mathcal{K}}^S(t-t^{{}^{\prime }}){\rho }_S(t^{{}^{\prime }}).$$

$$\mathcal{K}(t)=\mathcal{P}{\mathcal{U}}_S(t)L_I(t)(I-P){\mathcal{U}}^{{}^{\prime }}(t)L_I\mathcal{P}$$

$${\mathcal{U}}_S(t)=exp\{-i{L}_{S}t\},L_I(t)=exp[i(L_R+L_S)t]L_{I}exp[-i(L_R+L_S)t]$$

$${\mathcal{U}}^{{}^{\prime }}(t)=exp\{-{\int }_0^{t}d{t}^{{}^{\prime }}(I-\mathcal{P}{L}_I(t^{{}^{\prime }})(I-\mathcal{P})\}.$$

$$\dot{\rho (t)}=L(\rho )=i[H,\rho ]+\frac{1}{2}\sum _j([{\mathbb{L}}_j\rho ,{\mathbb{L}}_j^{+}]+[{\mathbb{L}}_j,{\mathbb{L}}_j^{+}\rho ])$$

$$A=\underset{\{n(\mathbf{r})\}}{\min }\left(A_{HKM}[n(\mathbf{r})]+\int V(\mathbf{r})n(\mathbf{r})d\mathbf{r}\right)$$

$$A_{HKM}[n]=A_{ni}[n]+E_H[n(\mathbf{r})]+E_{XC}[n(\mathbf{r})]$$

$$A_{ni}[n]=\underset{\{({\psi }_i(\mathbf{r}),f_i)\to n(\mathbf{r})\}}{\min }\sum _i\left(\frac{f_i}{2}\int |\nabla {\psi }_i(\mathbf{r}){|}^{2}d\mathbf{r}-TS(f_i)\right)$$

$$S(f)=-f\ln f-(1-f)\ln (1-f).$$

$$-\frac{{\nabla }^2}{2}{\psi }_i(\mathbf{r})+V_{KS}(\mathbf{r}){\psi }_i(\mathbf{r})={ϵ}_i{\psi }_i(\mathbf{r})$$

$$f_i=\frac{1}{1+e^{\frac{{ϵ}_i-\mu }{T}}}$$

$$\Phi =A-\mu N$$

$$\Phi [n(\mathbf{r})]=-\frac{2}{\beta }\ln det\left(1+e^{-\beta (\mathcal{H}-\mu )}\right)-\int d\mathbf{r}n(\mathbf{r})\left(\frac{\varphi (\mathbf{r})}{2}+\frac{\delta {\Phi }_{xc}}{\delta n(\mathbf{r})}\right)+{\Phi }_{xc}$$

$$A=\underset{\{n(\mathbf{r}),N_{\alpha }(\mathbf{r})\}}{\min }\left(A_{JDFT}[n(\mathbf{r}),N_{\alpha }(\mathbf{r})]+\int V(\mathbf{r})n(\mathbf{r})d\mathbf{r}+\sum _{\alpha }{V}_{\alpha }(\mathbf{r})N_{\alpha }(\mathbf{r})d\mathbf{r}\right)$$

$$A_{JDFT}[n,N_{\alpha }]=A_{HKM}[n]+A_{diel}[n,N_{\alpha }].$$

$$U(\mathbf{r})=\sum _i^M{f}_i(\mathbf{r})U_i(\mathbf{r}).$$

$$F_{\alpha \beta }=w(X_{\alpha })w(X_{\alpha })F_{\alpha \beta }^{AT}+[1-w(X_{\alpha })w(X_{\alpha })]F_{\alpha \beta }^{CG}$$

Full text

1 Introduction

Open systems that exchange matter with the environment represent a major challenge for theoreticians and simulators. In fact, the variation of the number of particles during the evolution of the system corresponds to a sudden change of the total amount of microscopic information that one must process and analyze within a consistent system-environment physical framework. In a previous work co-authored with Matej Praprotnik, [1] an overview of theoretical principles and (mostly classical) simulation techniques for systems with open boundaries available in the literature has been discussed. In [1] it was concluded that when modeling matter classically, there is a satisfactory understanding of how to treat the case of varying the number of particles on a numerical level as well as on a conceptual level, with the corresponding dimensional change of the phase space. On the other hand, departing from a classical representation and instead modeling this process in the quantum case becomes far more complex. In standard situations, e.g. in equilibrium, a classical particle arriving from the external environment into a system needs only to accommodate locally and does not modify abruptly the overall state of the system. The same concept does not hold in the case of quantum particles where the change in the number of particles completely redefines the quantum state of the system. Quantum particles are characterized by quantum correlations or better, by their entanglement, which implies that the most accurate knowledge of a system does not imply the most accurate knowledge of its parts [2, 3, 4, 5], thus the gain/loss of information for particles entering/leaving a subsystem (of specific interest) must be treated and interpreted with special care. In this work I will focus on the treatment of one specific, albeit relevant, class of quantum systems, that is many-electron systems that exchange particles, i.e. electrons, or atoms/molecules (electrons + classical nuclei), with a large environment. The description of the environment and its coupling to the (sub)system are the two main ingredients for the construction of a computational procedure that can simulate the process of exchange. The environment can be considered with its full electronic structure, that is, if rigorously treated, the system of interest is merely a subsystem of a fully resolved quantum large system. In reality, pragmatic approximations are used to fully resolve the environment at a reasonable computational price. In general, the conceptual advantage of this model is that the electronic coupling (system-environment) is explicitly taken into consideration, and as a consequence the electronic correlations between the subsystem and the environment are, with some degree of accuracy, explicitly considered into the electronic properties of the subsystem. The computational disadvantage is that one needs to treat either large systems, which would be in most of the cases prohibitive for current computational resources, or a small environment which is likely to suffer from artifacts due to the reduced size. Alternatively, the environment can be considered as an ideal statistical reservoir which is assumed to provide or adsorb particles according to the electronic chemical potential. It follows that the system is described as a quantum Grand Canonical ensemble without the need of explicitly calculating the electronic properties of the environment. The advantage in such a case is that one can focus only on the (sub-)system of interest, but the obvious disadvantage is that information regarding quantum correlations with the environment cannot be derived in any manner. From the computational point of view, the two categories outlined above can be differentiated in terms of calculations at a fixed number of electrons (energy minimization of the system+environment) and calculations at a fixed chemical potential and a varying number of electrons for the system of interest. The general conceptual framework outlined above leads us to the backbone of the paper. I have taken as guiding examples the theoretical and computational procedures employed in two relevant classes of applications: (i) in nanotechnology and (ii) in chemical and biological physics. More specifically, to the first category belong the subjects of nanoelectronics and electrochemistry, where the varying number of particles corresponds to the transport of electrons between a (molecular) system and the environment acting as an electrode. To the second category belongs the field of solvation chemistry, where a molecular subsystem exchanges molecules with the large thermodynamic bath in which the subsystem is embedded, for example biomolecules in water. Here the exchange of electrons occurs through the exchange of molecules (i.e. electrons and nuclei) and from a computational point of view, in order to make calculations more feasible, the environment shall be treated classically [8, 9, 10]. In this case, the combination of quantum mechanics of electrons and the thermodynamics/statistical mechanics of the classical scale becomes particularly delicate [1]. Regarding the methods of calculation, among the electronic structure techniques, Density Functional Theory (DFT) [6], due to its intrinsically statistical mechanical structure, facilitates the extension of its principles to the treatments of electrons in various ensembles and at the same time it allows for the use of functionals beyond the total energy. The previous statement can be verified by consulting the seminal book of Parr and Yang [7] where DFT is derived in terms of density matrix and corresponding statistical ensembles. As a consequence, DFT represents the most flexible and popular approach to the Grand Canonical treatment of electrons; a large amount of studies of chemical and physical systems where the DFT Grand Canonical approach has been used are present in the literature (see e.g. [11, 12, 13, 14, 15, 16] and references therein). In general, all of the applications reported above use DFT, in various forms and in combination with other techniques, as the electronic structure method of preference. Beyond DFT, among the most advanced (wavefunction-based) electronic structure techniques, one finds high level quantum chemical methods in Fock space. In fact, this is a natural approach to create and destroy electrons in a system [17]. In addition, Quantum Monte Carlo (QMC), given its stochastic/statistical nature [18], is flexible enough to treat the case of a varying number of particles and thus, treat the Grand Canonical ensemble [19, 20, 21, 22]. I will discuss the principles of the abovementioned QMC for one specific example, although it must be underlined that such an approach has not been explicitly used to treat the physical exchange of electrons, but has been employed mostly as a numerical trick that minimizes the size effects of the calculations, i.e. it improves the convergence of calculated quantities to their corresponding value in the thermodynamic limit. Nevertheless, as an electronic structure Grand Canonical technique, the QMC approach may be embedded in multiscale methodologies of the near future that will treat the physical process of varying the number of electrons and molecules. Finally, the paper is concluded with a discussion where the current ideas and related techniques are put in perspective. Possible combinations and modifications are suggested, which may be likely to optimize the methodology and hopefully push it beyond the current frontiers of applicability to physical systems.

2 Electron flow: nanoelectronics and electrochemistry

The passage of electrons from a molecule to its environment that acts as an electrode (and vice versa) is the prototypical situation occurring in nanoelectronics and electrochemistry. In nanoelectronics, the most popular case treated is that of a molecule that acts as a junction/bridge between two surfaces (metals or semiconductors) and allows the passage of a current of electrons from one surface to the other [23, 24, 25, 26]. In electrochemistry, a typical example is the flow of electrons between the reactants and a catalytic surface as a chemical reaction proceeds [27]. In the next two sections, I will analyze how these two situations are treated at methodological level.

2.1 Electron Transport in Molecular Junctions

Figure 1 depicts the typical structure of a simple molecular junction attached to a right and left electrode.

Differently from the other examples treated in this work (systems in equilibrium), the theoretical description of the transport of electrons from one electrode to the other is particularly challenging. In fact in such a case one deals with the more complex situation of non-equilibrium for a quantum many-body system. In the analysis of the methodological aspects relevant to our current discussion, I will follow the recent perspective paper of Thoss and Evers [25], which exhaustively traces the current state of the art in the field (see also references in the special issue about frontiers in molecular electronics [28]). They make a major distinction between weakly correlated systems, that is, systems where the energy transfer between charge carriers and the molecule can be neglected, and strongly correlated systems, where electron-electron or electron-phonon correlation have a significant impact. The weakly correlated regime, for its methodological aspects related to the system-environment treatment, is of particular interest to our specific analysis and will be discussed in the next section. For the strongly correlated regime I will not discuss the specific aspects of the strong correlations and its physical implications, instead I will restrict the attention to its treatment in the framework of a dynamic model of quantum open systems that exchange electrons with a reservoir.

2.1.1 Weakly correlated regime

The weakly correlated regime allows the description of a molecular junction within the approximation of quasi-static scattering potential and is valid in the regime of linear response, that is, the transport is treated as an elastic process. Within this picture, the electronic properties of interest can be calculated by Kohn-Sham Density Functional theory (KS DFT)[6]. The current-voltage characteristic, following the Landauer formula is [24, 29, 30]:

[TABLE]

where $V$ is the applied voltage and $E$ the energy levels, $f_{L,R}(E)=f(E-\mu_{L,R})$ is the Fermi distribution function of the lead at the corresponding chemical potentials $\mu_{L,R}$. In this context, the key quantity is the (molecular) transmission function, $T(E)$ which can be written as (see also [31]):

[TABLE]

The quantities involved in the formula are: ${\bf\mathcal{G}}^{-1}_{V}={\bf\mathcal{G}}^{-1}_{0V}-{\bf\Sigma}_{L}-{\bf\Sigma}_{R}$, where the resolvent matrix is:

[TABLE]

the sum is over the KS energies and orbitals, $\epsilon_{n}$ and $\phi_{n}$ for the uncoupled molecule calculated with a standard DFT approach. ${\bf\Sigma}_{L,R}$ are the self-energies, so that in the approximation of non-interacting particles, ${\bf\Gamma_{L,R}}=i({\bf\Sigma}_{L,R}-{\bf\Sigma}^{\dagger}_{L,R})$. The self-energy is the coupling term between molecule and leads, and can be formulated in terms of a hopping matrix: ${\bf\Sigma}_{L,R}={\bf t}_{L,R}{\bf g}_{L,R}{\bf t}^{\dagger}_{L,R}$. The elements of the hopping matrix are defined as: $t_{L,R}({\bf X}N,{\bf x}n)$, describing the hopping process of an electron in an orbital $N$ of an atom at position ${\bf X}$ of the molecule to an orbital $n$ of an atom at position x, in the lead. These elements are approximated by their bulk values and are calculated using independent DFT simulations of the lead (i.e. usually a large metal cluster without the molecule) and the molecule. Finally, the surface Green’s function, $<{\bf x}n|{\bf g}_{L,R}|{\bf x}^{{}^{\prime}}n^{{}^{\prime}}>$, that is the effect of broken translational symmetry on the density of states of the metal occurring at the interface plane perpendicular to the lead-molecule direction, is also calculated in the same DFT simulation of the metal cluster used to find the elements of the hopping matrix. This description reports the very essential aspects to give an idea of the methodological procedure used in the field; there are of course technical improvements (e.g. “the extended molecule” model, that is the molecule and some atoms/layers of the leads) and further conceptual developments of this basic picture. Advantages, limitations and open problems are well described and discussed in Reference [25]. For the current focus the relevant methodological aspects are: (i) the model of environment/lead, is fully resolved in its electronic structure, although in the pragmatic approximation of large clusters, (ii) its coupling to the system/molecule is done through the hopping process by which electrons are exchanged. This latter aspect requires the construction of a matrix obtained by (iii) modeling the system of interest/junction as an “isolated molecule” treated, as for the bulk of the metal, with standard DFT approaches at a constant number of electrons. The approach underlined in this section makes use of a static picture where the dynamics of electron exchange between molecule and lead is expressed through the probabilistic event of the hopping. An interesting methodological variation for the calculation of $T(E)$ where the leads are characterized only by their chemical potentials (given as an input), without the need of knowing their electronic structure, has been proposed by Arnold, Weigend and Evers [32]. They develop a general scheme, which despite requiring standard quantum chemistry calculations for the bridging molecule (or extended molecule), in practice gives as a result electronic properties calculated at fixed chemical potentials (of the left and right lead) and a variable number of electrons. This is an effective and computationally robust way to realize the Grand Canonical Ensemble in an implicit manner from a series of (properly looped) standard fixed-particle quantum calculations. Its potential application can go far beyond the case of electron transport and can certainly be used for systems in equilibrium that will be described later on. Beyond the static approaches discussed so far, there exists a class of methods that explicitly considers the electron dynamics and can treat also the strongly correlated regime. They are based on the Liouville equation for the density matrix of the subsystem of interest (reduced density matrix of the junction). In essence, such an equation describes the dynamical behaviour of electrons in a molecule considered as an open quantum system; such an approach is described in the next section.

2.1.2 Electron Dynamics: Liouville Equation for the Density Matrix and related methods

In an interesting work of Emch and Sewell [33], the Liouville equation for the time evolution of a quantum subsystem (S) embedded in an ideal bath (R) is rigorously derived in terms of the evolution of the density matrix of the subsystem. The interesting aspects for the focus of this work is that the authors suggest a coupling term between the S and R also for the case in which the bath is a reservoir of particles (in addition to energy). The essence of the Emch-Sewell model is the filtering of the microscopic degrees of freedom of R through the Zwanzig projector [34]. The process of filtering leads to effective actions of such microscopic degrees of freedom in terms of actions of statistically averaged macroscopic quantities. The standard Liouville-Neumann equation: $\frac{d}{dt}\rho=-iL\rho$, with $\rho$ the density matrix of the system and $L=[H,*]$ the Liouville operator with Hamiltonian, $H$, is transformed through the projector operator $\mathcal{P}$ into the equivalent master equation for $\mathcal{P}\rho$:

[TABLE]

$I$ is the identity operator and $\mathcal{U}(t)=exp[-i(I-\mathcal{P})L(I-\mathcal{P})t]$. S and R are initially uncorrelated/independent $\rho(0)=\rho_{R}(0)\otimes\rho_{S}(0)$. The initial state of R is given by the measurement of the set of macroscopic variables, thus $\mathcal{P}\rho(0)=\rho(0)$ and since the average/macroscopic quantities of the reservoir are time independent, then $\mathcal{P}\rho(t)=\rho_{R}(0)\otimes\rho_{S}(t)$ thus Eq.4 can now be written as a self-contained equation for S:

[TABLE]

considering the trace with respect to R on both sides one obtains the master equation for $\rho_{S}(t)$:

[TABLE]

$\mathcal{K}^{S}(t)\rho_{S}(t)=Tr_{R}\{\mathcal{K}(t)\rho_{R}(0)\otimes\rho_{S}(t)\}$. As anticipated above, the interesting point for our focus is that $L^{S}_{eff}$ and $\mathcal{K}(t)$ can be used to define the exchange of particles between R and S. The Hamiltonian of the system can be divided in three parts: $H=H_{R}+H_{S}+H_{I}$, the relevant part is the latter term, that is the Hamiltonian of interaction between R and S. Correspondingly, one has the related Liouville operators $L_{R},L_{S},L_{I}$. The authors define $H_{I}=\int_{\Omega_{R}}dx\int_{\Omega_{S}}dyV(x,y)J_{R}(x)\otimes J_{S}(y)$, with $x,y$ configuration coordinates of R and S respectively, $\Omega_{R},\Omega_{S}$ the volumes occupied and $J_{R}(x),J_{S}(y)$ operators acting on the Hilbert subspace of R and S respectively. Such operators represents intensive variables, for example particle number, which could be function of the creation and annihilation operators for particles in R and in S. The intensity of their action is regulated by $V(x,y)$, a potential of (direct) interaction between R and S (e.g Voltage in a junction-lead system). Furthermore, $L^{S}_{eff}=L_{S}+L^{S}_{I}$, with $L^{S}_{I}\rho=[V_{S},\rho]$ and $V_{S}=\int_{\Omega_{S}}dy\langle V(y)\rangle_{0}J_{S}(y)$. The kernel is defined as follows:

[TABLE]

with

[TABLE]

and

[TABLE]

The treatment is (in principle) exact, but the memory kernel , $\mathcal{K}_{S}(t-t^{{}^{\prime}})$, cannot be analytically determined thus it must be approximated. As suggested in [25], recent research has brought substantial advances in the field [35, 36]. A similar formulation in which the Liouville-Neumann equation considers a stochastic coupling to a bath is the so-called Kossakowski-Lindblad equation [37, 38]:

[TABLE]

with $H$ the system Hamiltonian, $\mathbb{L}_{j}$,$\mathbb{L}_{j}^{+}$ operators that carry the interaction of the system with a reservoir (Lindblad operators). The term $\sum_{j}([\mathbb{L}_{j}\rho,\mathbb{L}^{+}_{j}]+[\mathbb{L}_{j},\mathbb{L}^{+}_{j}\rho])$, gives Eq.10 the form of a rate equation (quantum jumps of the system under the action of the environment), $[\mathbb{L}_{j}\rho,\mathbb{L}^{+}_{j}]$ and $[\mathbb{L}_{j},\mathbb{L}^{+}_{j}\rho]$ can be interpreted as transition rates between two events (e.g. particles exchange with a reservoir and transition from $N$ to $N^{{}^{\prime}}$). The use of the creation and annihilation operators inevitably leads to the treatment of the problem in Fock space, within standard one-particle orbitals of quantum chemistry, or Bloch states and plane waves [39, 40], further routes to describe the electron transport originating from the master equation and references for specific determinations of the memory kernel are discussed in Ref.[25, 41]. It must be reported that methods based on the direct calculation of the time evolution of the wavefunction of the system, rather than its density matrix, are also present in the literature, however in such a case, the wavefunction must be calculated for the whole system (molecule+leads). In order for this calculation to be computationally feasible, the leads must be approximated by a finite system [42]. The most popular approach in such a context is the time-dependent Hartree method (MCTDH) [43, 44]. For electron transport, the exchange of electrons between system (molecule) and environment (leads) occurs through the change of orbital occupations (from lead to molecule and molecule to lead) [45]. In such case correlations are explicitly taken into account, however, the computational effort may be sizable [25]. Beyond electron transport, within the same wavefunction approach, attempts to restrict the treatment only to a subsystem and add the rest (environment) as e.g. a sink of electrons (i.e. for treating ionization processes), led to different technical and conceptual problems; the most relevant is the fact that the modified Hamiltonian, with a sink: $H+i\Gamma$ with $\Gamma$ a one body potential vanishing on the domain of the molecule, is a non-Hermitian operator [46]. Overcoming such a limitation implies the fulfillment of the requirement that the quantum dynamical semigroup for the time evolution must be trace preserving, Markovian and strictly positive. In practice, one follows the application of the Gorini-Lindblad theorem [47, 37], and thus arrives at the Kossakowski-Lindblad equation discussed before. The previous statement implies that the density matrix, when the number of electrons is varying in time, is the proper quantity to consider rather than the direct wavefunction of the system. The use of the Liouville equation discussed so far implies two distinct computational strategies for the partitioning and definition of system and environment. In the approach of the projector operator, the density matrix should, in principle, be defined for the whole space and the physics of the subsystem is obtained by the action of Zwanzig projector operator. It is clear that while the accuracy of the results would be very high, its computational costs are likely to be prohibitive for systems of reasonable size (see also comments about Fock space), and of course reasonable approximations are possible and needed (see e.g. [25] and references therein). The Liouville equation in the approach of Kossakowski-Lindblad, requires instead only the definition of a transition probability without the necessity of explicit calculations of the electronic properties of the reservoir. Though, of course, its accuracy is directly linked to the construction of a physically well funded stochastic term, for example from transition rates determined experimentally. In general, as underlined before, the treatment of a problem in Fock space is the natural approach to describe the passage of electrons from one system to another, but in practice the computation may be difficult. In fact, it requires the construction of the density matrix through the treatment of a large number of states to which electrons can be located according to the action of the creation and annihilation operators. This means that one needs a predefined reasonable “active space” [17]. If the explicit time dependence of the electron flow is not of primary relevance, an alternative to the use of Fock space is given by DFT treated in the Grand Canonical ensemble. Calculations for electrochemistry are often done following this approach, and this subject will be treated in the next section.

2.2 Grand Canonical Density Functional Theory and its application to Electrochemistry

Electrochemical systems can typically be schematized as an electrolyte, usually consisting of ions in a solvent, and a surface, acting as an electrode, which promotes an electrochemical reaction [27, 48]. Electrons are absorbed from or injected into the solution as the reaction proceeds (see also Figure 2 for a pictorial example).

From the point of view of simulation, modeling such systems are extremely challenging: the full quantum treatment of the system is computationally prohibitive, while a pragmatic partitioning of the total system as subsystem of interest (where the chemical reaction occurs), and the environment (electrolyte and electrode), allows multiscale simulation techniques to be applied with success [15]. Later on I will describe a specific model/example that treats the electrolyte, but for the moment, the interesting aspect is the treatment of the electrode as a Grand Canonical reservoir of electrons for the (sub)system of interest. In fact, such a view can be generalized to the problem of a charged surface or a surface in an external electric field which can be treated as a system at constant chemical potential, $\mu$ and varying number of electrons $N$ (for methodologically-oriented work see e.g. [49, 50, 14, 15, 16]). The use of DFT calculations with constant $\mu$ and varying $N$ leads to the introduction of the (exact) Helmholtz free energy, $A$, for interacting electrons at finite temperature in an external potential $V({\bf r})$. $A$ satisfies the so-called Hohenberg-Kohn-Mermin variational theorem [6, 51] (see also Reference [15], of which I will follow the formalism):

[TABLE]

with $A_{HKM}[n]$ the universal free energy functional of the electron density $n({\bf r})$, for any external potential $V({\bf r})$ (in atomic units). As for the universal energy functional of Hohenberg and Kohn, $A_{HKM}[n({\bf r})]$ is not known and can be only approximated:

[TABLE]

where $A_{ni}[n]$ is the non-interacting free energy, $E_{H}[n]$ is the Hartree term and $E_{XC}[n]$ is the exchange and correlation energy. The relevant quantity here (compared to standard DFT) is:

[TABLE]

that is, the minimum over all the single-particle orbitals $\psi_{i}({\bf r})$ and the corresponding occupation factors $f_{i}\in[0,1]$ which leads to the electron density: $n({\bf r})=\sum_{i}f_{i}|\psi_{i}({\bf r})|^{2}$. $T$ is the temperature and $S(f_{i})$ is the single-particle entropy:

[TABLE]

The minimization leads to the single-particle KS equations:

[TABLE]

for the stationarity w.r.t. $\psi_{i}({\bf r})$ and to the Fermi Occupation condition:

[TABLE]

for the stationarity w.r.t. $f_{i}$, with $V_{KS}({\bf r})=V({\bf r})+\frac{\delta}{\delta n({\bf r})}(E_{H}[n]+E_{XC}[n])$. For standard DFT, $\mu$ is a Lagrange multiplier and the fixed number of electrons is assured by the condition: $\sum_{i}f_{i}=N$. In a Grand Canonical ensemble instead, the free energy to minimize is not $A$ but the Grand Free Energy:

[TABLE]

In such a case, for the optimization procedure applied above, one modification is required, that is the Lagrange multiplier term $-\mu(\sum_{i}f_{i}-N)$ is replaced by $-\mu\sum_{i}f_{i}$. Such modification removes the constraint on fixed $N$ and implements the Legendre transformation from $A$ to $\Phi$; $\mu$ is now given as an input while $N$ is variable. The grand potential of the electrons can be expressed also in a different formalism, for example following Alavi et al. [49]:

[TABLE]

with $\Phi_{xc}$ the finite-temperature exchange-correlation grand potential, $\phi({\bf r})$ the Hartree potential, $\beta=\frac{1}{k_{B}T_{e}}$ the electronic temperature parameter, $\mathcal{H}=-\frac{1}{2}\nabla^{2}+V({\bf r})$ the one-electron Hamiltonian, and the effective density-dependent potential, $V({\bf r})=V_{ext}({\bf r})+\phi({\bf r})+\frac{\delta\Phi_{xc}}{\delta n({\bf r})}$, with $V_{ext}({\bf r})$ the external potential. The exponential form can then be efficiently evaluated through the Trotter approximation as a product of $P$ high temperature matrices: $e^{-\beta\mathcal{H}}=\left(e^{-\beta\mathcal{H}/P}\right)^{P}$, with $P$ a large integer, so that $\epsilon=\beta/P$ is small and one can write: $e^{-\epsilon(K+V)}=e^{-\epsilon V/2}e^{-\epsilon K}e^{-\epsilon V/2}+\mathcal{O}(\epsilon^{3})$; based on such an idea, it was possible to devise efficient linear system-size scaling schemes [52, 53, 50]. In general, compared to fixed $N$ calculations, those at constant $\mu$ carry several technical problems. For example, the large fluctuations of $N$ (and $n({\bf r})$) at the initial stage, the need of compensating charges because of the finiteness of the slabs representing the surface and periodic boundary conditions [11, 14, 15], or the fact that $E_{XC}[n]$ are defined for integer numbers of electrons [54, 55, 56]. However, technical solutions were made available and thus the method can be routinely used nowadays for applications. The approach outlined above makes possible the introduction of a reservoir of electrons that adds or removes such particles as, e.g., a chemical reaction on a surface/electrode proceeds, so that our system of interest consists of the reactants while the surface represents the (active) environment. However, as underlined before, the effect of the solvent often plays a key role and its quantum treatment would be computationally prohibitive. The most popular approximation consists of employing continuum solvation models [57, 58], however an interesting proposal based on the so-called joint density-functional theory (JDFT) [59, 60] has been elaborated in the context of Grand Canonical DFT reported above [15]. JDFT consists of a variational theorem similar to the Hohenberg-Kohn theorem which allows for the description of the free energy of a solvated system in terms of the electron density $n({\bf r})$ of the solute and of nuclear densities for the solvent, $N_{\alpha}({\bf r})$, where $\alpha$ indicates a nuclear species:

[TABLE]

similarly to the Hohenberg-Kohn-Mermin functional, $V({\bf r})$ is the external electron potential, $V_{\alpha}({\bf r})$ is the external potential for the nuclei of the solvent and $A_{JDFT}[n({\bf r}),N_{\alpha}({\bf r})]$ is a universal functional independent of the two potentials and can be separated as:

[TABLE]

The solvent is described in terms of average density and not in terms of individual atomistic configurations which would imply expensive sampling computations with e.g. molecular dynamics or Monte Carlo, while the system of interest is treated at full quantum (Grand Canonical) accuracy. However, often the solvent and the corresponding configuration space that it can access are of crucial importance for the chemical or physical event of interest and thus the continuum approximation is far too drastic and it requires an explicit molecular treatment of the solvent [61]. In general, the region of interest is usually localized in space, e.g. the solvation region of a solute, thus the quantum accuracy may be needed only in a small portion of the system while the rest of the system can be treated at a coarser (classical) level. This is a typical scenario for, e.g. biophysical systems such as proteins and membranes in water, whose relevance for the current academic research as well as for technology has stimulated the development of multiscale computational methodologies, such as Adaptive Quantum Mechanics/Molecular Mechanics (A-QM/MM) [9] and Adaptive Resolutions Simulation (AdResS) [1]. In such cases, the environment is a classical system and provides/takes molecules from the quantum region, thus the exchange of electrons is strictly related to the exchange of entire molecules. This kind of approach is the subject of the next section.

3 Molecules in Solution: Exchange of Molecules between system and reservoir

Many interesting systems in chemical physics and in particular in biochemistry are characterized by events that happen locally, e.g. in solution. In such a context, the interesting event takes place in a specific subregion of the system (where electrons play a key role) embedded in a larger thermodynamic bath of, e.g., solvating molecules. The proper (full quantum/electronic treatment) in space and time, of such a class of systems is in most of the cases prohibitive; however, a decisive step forward was taken following the idea of space partitioning of Warshel and Karplus [62] and of Birge, Sullivan and Kohler [63], that is, each region is treated at the most convenient resolution, i.e. quantum if electrons play a major role or coarser otherwise. Later on, the idea of partitioning was extended further by Warshel and Levitt [64], who described the enzymatic reaction of lysozyme; their seminal work is nowadays recognized as the first example of the so-called Quantum Mechanics/Molecular Mechanics (QM/MM) method. This method has represented a technical revolution in the field because it allows for the treatment of systems that, up to that moment, were thought to be intractable at quantum level. For this reason, the pioneers of the technique, Martin Karplus, Michael Levitt and Arieh Warschel were awarded the Nobel Prize in Chemistry in 2013. Today, QM/MM represents a very precious computational tool for tackling problems in various field of chemistry and physics with applications in a broad range of subjects and disciplines (see e.g. [65, 66]). However, in the standard QM/MM method, the quantum region and the classical region are fixed, that is, there is exchange of energy, but not of matter. The physical implication is that the QM region represents an artificial Canonical ensemble whose accuracy increases by increasing the size of the QM region and checking that the effect of the environment is negligible. The most recent development in the field attempts to technically remove the constraints of a fixed number of molecules in the QM region and thus, allow the exchange of molecules between the quantum and the classical region (see e.g. [9]). I have already mentioned such an approach, referring to it as: “Adaptive Quantum Mechanics/Molecular Mechanics (A-QM/MM)”; its advantages and limitations will be discussed in the next sections. In parallel, within the context of classical multiscale simulation, (molecular) adaptive resolution methods for classical systems have been developed [1]. Typically such methods concurrently couple regions at (classical) atomistic resolution with regions at coarse-grained resolution in such a way that when a molecule passes from one region to the other, it slowly changes resolution without creating sizable physical artifacts. The next frontier for classical adaptive resolution methods is the extension to molecules with electrons. In this context, a computational framework that I have recently proposed, electronic Quantum Adaptive Resolution Simulation (el-QM-AdResS) [10], where the classical adaptive resolution technique is combined with the electronic scale, is discussed in the following sections. Finally, a discussion of the similarities, differences and possible synergies of A-QM/MM and el-QM-AdResS concludes the chapter.

3.1 A-QM/MM: Partitioning the system on-the-fly

In the following discussion, I will take as a main reference the exhaustive review by Hai Lin and collaborators [9] where the state of the art of the field is critically analyzed (for an additional complementary view see also Ref.[8]). A-QM/MM, as anticipated above, is an extension of the standard QM/MM method. The essence of A-QM/MM consists of an on-the-fly reclassification of QM and MM atoms/molecules during the simulation so that the molecules included in the QM or in the MM region can be dynamically updated according to the evolution of the system in space. A straightforward approach consists of an “abrupt” scheme of reclassification (see Figure 3): classical molecules that at time $t_{0}$ lie at the border of the QM region may directly enter into the QM region at time $t_{1}=t_{0}+\Delta t$, ($\Delta t$ time step of the simulation) and be reclassified as QM molecules for the calculations of the next time step; equivalently, QM molecules lying at the border with the MM region and moving at time $t_{1}$ into the MM region are reclassified as classical for the calculations of the next time step.

It was found out that such a sudden change implies the hopping between different energy surfaces whose discontinuity causes numerical instabilities and artificial results [71]. For such a reason, the latest generation of A-QM/MM methods are centered around more involved computational schemes based on the introduction of a buffer region. At each time-step of the simulation, the buffer region is partitioned in different subsets and the standard QM/MM interactions are defined for all possible subsets (see Figure 4). Next, for each partition, the molecules of the corresponding subset of the buffer are included in the QM region.

Finally, for each of these “extended” QM regions a standard QM/MM calculation is done. The total potential is then defined as a weighted average of the individual potentials from each individual simulation:

[TABLE]

Here, each $U_{i}({\bf r})$ corresponds to one of the $M$ partitions of the system in a group of QM molecules and a group of MM molecules and $f_{i}({\bf r})$ is a switching function depending on the coordinates of the single molecules. The switching function is constructed with the intuitive but well justified idea that the quantum energy of molecules in the buffer at larger distances from the center of the region of interest (active site) contributes less than the energy of molecules which are closer to the active site. This is the essential technical and conceptual point of the method. In fact, the dynamical evolution produced by the energy surface, can now be updated and it creates a new configuration on which, in turn, the partitioning step and the corresponding weighted potential (or forces) are applied for the time evolution of the next time step. From the technical point of view, the empirical approach to the smoothing of the coupling at the interface (either via potential or via inter-molecular forces) turned out to be very powerful and, in my view, can be considered a very relevant step forward in the development of truly multiscale simulation technologies for condensed matter and chemical physics, as testified by a large number of successful applications [72, 73, 74, 75, 76, 77] (see also additional references in the topical reviews [8, 9]). Yet, Lin and coworkers point out a series of technical problems, but above all conclude that despite the large number of encouraging results, A-QM/MM cannot be considered a truly predictive tool since its results always require a case-by-case validation, i.e. reproduction of some reference results (from experiment or larger QM/MM calculations). In the next section, I will discuss conceptual aspects that are mandatory for the construction of a A-QM/MM approach with automatic a priori physical control criteria that can transform A-QM/MM into a predictive tool without the need for external validation.

3.2 A-QM/MM: Physical validation and the necessity of a Grand Canonical view

In Reference [9], it is explicitly stated that, in principle, the QM region can be made as small as possible. This statement is certainly true from a technical point of view, however it carries the first conceptual bug of the current QM/MM methods (in general). In fact, if one is interested in a realistic quantum description of a specific region of interest, then the coupling energy between the QM region and the MM region must be negligible compared to the energy of the QM region, otherwise the electronic spectrum is essentially determined by the classical part. It is obvious that in QM/MM the electronic wavefunction and corresponding energy spectrum of the QM region cannot be the same as if the region was in a full QM environment. In good approximation, this situation holds only in the case of localized systems/properties, that is, when there is a negligible coupling energy with the environment (i.e. a reasonably separable Hamiltonian). This would be the first internal criterion of control for the validity of a QM/MM calculation. If we now take a step forward and move to the A-QM/MM approach, there are further conceptual problems. A very relevant one concerns the physical validity of the approach used in the partitioning scheme. Specifically: as said before, at quantum mechanical level, in terms of wavefunctions and electronic spectra, an open QM region embedded in a MM reservoir, in general, cannot be equivalent to a QM region embedded in a full QM reservoir. The equivalence can only hold from the quantum statistical mechanics point of view, if (and only if) the QM region of the A-QM/MM, is interpreted as a Grand Canonical region at the same macroscopic (thermodynamic) and electronic chemical potential of a QM region embedded in a full quantum environment. In this respect, the concept of “Nearsightedness of electronic Matter” introduced by W.Kohn [78, 79], provides a clear physical principle for the definition of a subsystem embedded in a larger environment. In fact it states that, at fixed chemical potential, local electronic properties depends on the effective external potential (i.e. the environment) only at nearby points. The effects of the external potential beyond a certain distance are negligible for the local properties. Furthermore, in a recent work of Fias, Heidar-Zadeh, Geerlings, and Ayers, it has been shown that the response kernel to the environment is local at constant chemical potential, thus in any partitioning scheme one should use the concept of constant chemical potential [80]. To my knowledge, the current partitioning schemes of A-QM/MM do not use any of the criteria of control listed above. Moreover, in a recent work by Miranda-Quintana and Ayers about electronic systems with a varying number of electrons in DFT, it is discussed the possibility that interpolations of property-values between electron numbers is not consistent with a physical ensemble average [81].

Their finding may be extended to the current discussion and would imply that any average property obtained by averaging over the different $M$ partitions of the A-QM/MM system, where each calculation is done minimizing the energy $E_{i}$ at a given (fixed) number of electrons $N_{i}$, is inconsistent with statistical mechanics.

An extension of the A-QM/MM method that includes the internal criteria of control described in this section, has been proposed by myself and is based on the inclusion of molecules with electrons within the classical atomistic/coarse-grained Adaptive Resolution scheme (AdResS) in its Grand Canonical version (GC-AdResS) [10]. The essence of the idea is reported in the next section.

3.3 el-QM-AdResS: A Grand Canonical electronic system embedded in a classical reservoir

The Adaptive Resolution Simulation technique (AdResS) [67, 68] for classical systems allows molecules to transform their resolution according to the region in space where they are instantaneously located. The root model of the computational algorithm consists of a space dependent interpolation for the force between two molecules, $\alpha,\beta$ (see Fig.5):

[TABLE]

where $F_{\alpha\beta}^{AT}$ is the atomistic force and $F_{\alpha\beta}^{CG}$ is the coarse-grained force.

$w(x)$ ($x$ is the coordinate of the center of mass of a molecule) smoothly goes from 0 to 1 in a transition region $\Delta$. The system is embedded in a thermostat that assures a target thermodynamic equilibrium. The conceptual solidity of the original AdResS method is then enhanced by the addition of a force, ${\bf F}_{th}(x)$ (thermodynamic force), which acts on the center of mass of the molecule in $\Delta$. Such a force assures that, in situations of equilibrium, the effective chemical potential of the whole system corresponds to that of the (target) atomistic resolution [82, 69, 70]. ${\bf F}_{th}(x)$, has been derived within a rigorous statistical mechanics model for molecular systems with open boundaries [83, 84, 69, 85, 70, 86, 1] (thus AdResS became Grand Canonical (GC-) AdResS) and it was found out that a numerically convenient way to calculate ${\bf F}_{th}(x)$ during a simulation (equilibration run), consists of expressing it as the gradient of the number density of the molecules in an iterative form [83]: $F_{k+1}^{th}(x)=F_{k}^{th}(x)-\frac{M_{\alpha}}{[\rho_{ref}]^{2}\kappa}\nabla\rho_{k}(x)$, $M_{\alpha}$ is the mass of the molecule, $\kappa$ a (conveniently) tunable constant, $\rho_{k}(x)$ is the molecular density as a function of the position in $\Delta$ at the $k$-th iteration and $\rho_{ref}$ is the density of reference, decided a priori according to the thermodynamic state point at which we wish to do the simulation. The convergence criterion depends on the accuracy required for the simulation but, based upon experience, $|\rho_{final}-\rho_{ref}|$ should always be below $10\%$ in $\Delta$. The accuracy of this method has been proven over the last ten years over a large range of systems and problems, see e.g. [1] for an overview, and Refs.[87, 88] for the most recent applications. In this context, it is important to notice that the coupling scheme allows one to interface any classical molecular representations (e.g. two different atomistic models, atomistic and path integral representation of molecules [89, 90], to mention a few). There also exists a version based on a global Hamiltonian [91, 92], which is, in essence, technically equivalent to the method outlined above (although conceptually confusing, see discussion in [1]). Given the technical equivalence, it is not a surprise that calculations repeated with the Hamiltonian scheme give the same results obtained years before with the (GC-)AdResS method (compare [89, 90], with [93, 94]; [95] with [96]; and [85] with [97]).

In order to include electrons, the theoretical framework of el-QM-AdResS, proposed in Reference [10], is based on two simple concepts (see also pictorial representation in Figure 6):

• •

Use the classical GC-AdResS scheme as a computational interface for the quantum region (i.e. equivalent to the MM region of a QM/MM scheme).

GC-AdResS will take care of providing and removing the (classical) nuclei (i.e. the skeleton of the molecule) in the QM region according to the macroscopic (thermodynamic) chemical potential.

• •

The electrons of the QM region are treated as a system in contact with an ideal (Grand Canonical) reservoir of electrons whose exchange of particles is regulated by a properly chosen (see explanation later) electronic chemical potential.

The scheme follows three automatic criteria of control for assuring physical consistency of the calculation:

• •

The coupling between the QM region and the classical reservoir must be such that: $\langle H_{QM-Ad}\rangle<<\langle H_{QM}\rangle$ at any time, with $H_{QM-Ad}$ the coupling Hamiltonian of the QM region with the GC-AdResS interface (as an example of application of this criterion in GC-AdResS, see [70, 89]).

• •

The macroscopic chemical potential must be such that $\mu_{macro}^{QM}=\mu_{macro}^{GC-AdResS}$, i.e. the whole system is at the same macroscopic chemical potential at any time.

• •

The electronic chemical potential in the QM region must satisfy the condition: $\mu_{el}^{QM}=\mu_{el}^{ref}$, where $\mu_{el}^{ref}$ is the electronic chemical potential the QM region would have if the whole system was treated at quantum resolution. For example, for a molecule solvated in water, the QM region would be composed of the solute and the water molecules of the first solvation shells embedded in bulk water, thus the electronic chemical potential of reference would be that of bulk water.

The first condition is a simple check that assures the choice of a physically consistent size of the QM region. The second condition is fulfilled by the thermodynamic force, calculated with the same formula used in standard GC-AdResS for the region $\Delta$, but now calculated for (and applied in) the MM region at the interface with the QM region. Finally, the third condition can be implemented by performing electronic structure calculations for the QM region at constant $\mu_{el}$, with the number of electrons being the variable in the (Grand Canonical) energy functional minimization (in the same spirit of the Grand Canonical DFT calculations of the previous sections). The electronic chemical potential can be determined from separate calculations of prototype systems that represent a reasonable bulk/environment (see example of liquid water before). The main advantage of the el-QM-AdResS scheme consists in the assurance regarding the statistical mechanical consistency of the results; however, it carries with it open questions regarding the integration of its conceptual aspects and its technical implementation. In fact, from the technical point of view, one major challenge is the abrupt interface between the QM and the GC-AdResS region. It is possible that the thermodynamic force, acting in the MM region directly adjacent to the QM region which, in essence, equalizes the molecular number density in an iterative process (equilibration run), may not converge as expected. It is encouraging though that recent tests have proved that this may not be a problem; in fact, an abrupt interface can be turned into an advantage [98]. Another challenge is the fact that $N$, the number of electrons in the QM region, may not be integer, which implies that some fractional charge is delocalized in the classical region and it must be taken into account. In the meantime, some technical solutions have been proposed [10] and probably the optimal development of this idea is to combine it with the already well developed technical A-QM/MM schemes. A brief discussion about this possibility is considered in the next section.

3.4 el-QM-AdResS + A-QM/MM: A path to a technically powerful and conceptually self-validating innovative scheme

el-QM-AdResS is, as a matter of fact, technically similar to A-QM/MM, and some work in merging AdResS with A-QM/MM has been already made in the group of Rosa Bulo [99]. A proposed way to proceed further would be to extend the criteria of physical validity used in el-QM-AdResS to A-QM/MM and, from the other side, use a buffer region and the corresponding partitioning scheme of A-QM/MM in el-QM-AdResS. The resulting method would be a new A-QM/MM where each extended QM region, corresponding to each partition, is treated within the (same) $\mu_{el}$ constant scheme (and implicitly, through the thermodynamic force, at the same $\mu_{macro}$). As a consequence, the average potential (or forces) calculated over all partitioning schemes in the buffer region would correspond to an average over different, but equivalent (i.e. same $\mu_{el}$ and same $\mu_{macro}$), Grand Canonical systems, and thus, it would have a solid statistical mechanics justification. At the same time the technical problem of the abrupt interface in el-QM-AdResS could be solved by the smoothing technique of the buffer region, where the thermodynamic force can also be applied and smoothly calculated. Moreover, GC-AdResS allows a further coupling of the atomistic resolution to coarse-grained models and beyond (continuum models) [100, 101], thus, the fully-realized A-QM/MM approach would represent a truly multiscale method.

Most of the QM calculations are done within the DFT scheme, however, one can, in principle, go beyond DFT and aim for a higher level of accuracy. For example, Hofer and Hünenberger have recently explored the idea of going beyond DFT and using the resolution-of-identity second-order Møller-Plesset perturbation (RIMP2), which is a sort of first level beyond DFT for describing electron correlations [102]. In this perspective, it is important to also consider other electronic structure methods that can accurately describe electron correlations but, given the context of this work, also assure exchange of matter, and thus calculations involving a varying number of electrons. In the previous chapters, quantum chemistry approaches in Fock space have already been discussed, however, one may also extend the discussion to other methods such as Quantum Monte Carlo (QMC). As previously anticipated, QMC techniques in the Grand Canonical ensemble are developed and used mostly for technical reasons (alleviating size effects), but in any case, in this context, they may be included in a scheme such as el-QM-AdResS or A-QM/MM. A discussion of QMC methods where the number of electrons can vary is reported in the next section.

4 Quantum Monte Carlo and a varying number of electrons

Quantum Monte Carlo methods can be classified among the most accurate techniques for describing electron correlations. There are multiple QMC approaches, each with its own positive and negative aspects; however, in general they are all closely related. Here, I will not provide an overview of specific QMC techniques (for an updated discussion see Reference [18] and references therein) but I will instead focus on some generic features that can be applied to different QMC approaches and involve the possibility of carrying out QMC calculations in a Grand Canonical ensemble. The work of reference which I will follow is the paper by Lin, Zong and Ceperley who introduced the so-called “twist average Grand Canonical ensemble” (TA-GCE) [21].

TA-GCE has been developed in order to reproduce results of the thermodynamic limit with calculations involving a (relatively small) finite number of electrons, thus circumventing finite size effects. The starting point of this idea is the observation that in QMC the phase factor of the many-electron wavefunction is usually assumed to return to the same value if electrons cross the periodic boundary and return to the same position. The technical consequence of this assumption is that, for delocalized electrons, the convergence of properties to their thermodynamic limit is very slow. In Reference [21], Lin, Zong and Ceperley proposed an alternative that allows electrons to take up a phase angle when they cross the periodic boundary. This is the basic concept of twist averaged boundary conditions (TABC); the pre-existing literature background they refer to is a work by Gros [19] who has shown that TABC, applied to the Hubbard model, gives exact results in the Grand Canonical ensemble for non-interacting particles. Thus, for the focus of this paper, this is a procedure of interest for the treatment of a varying number of electrons. In the following, I will report the basic ingredients of the idea of Lin, Zong and Ceperley. In non-interacting homogeneous systems with periodic boundary conditions, plane waves describe single particle states. For simplicity of, exposition the spinless case is considered, so that the single state of the system in a, e.g., cubic box of linear size $L$, can be written as: $\phi({\bf r},{\bf k})\propto e^{i{\bf k}{\bf r}}$. In order to satisfy the TABC one must have: ${\bf k}_{\bf n}=(2\pi+\theta)/L$, with ${\bf n}$ an integer vector. In the Canonical ensemble, the $N$ lowest energy states form the ground state and a Slater determinant of such states gives the ground state wavefunction of the system. In the Grand Canonical ensemble one has that: $E_{\alpha,N}=E_{\alpha,N}(\theta)$, where $(\alpha,N)$ labels the quantum states. The probability of a given state is: $P(\alpha,N)\propto exp\{-\beta[E_{\alpha,N}(\theta)-N\mu]\}$, with $\mu$ the chemical potential of the system. In the ground state, that is for $\beta\to\infty$, the wavefunction must not optimize the energy $E_{\alpha,N}$, as in the Canonical ensemble, but the quantity: $E_{\alpha,N}(\theta)-\mu N$. Next, they show that the occupation number is exactly what one gets if the calculation were to be done in the thermodynamic limit, and thus TABC (now renamed TA-GCE) carries no size effects. The procedure implies that the number of electrons varies, since, for a given $\theta$ and Fermi wave vector, the number of occupied states will vary. The fluctuations in the particle number is then derived and it is shown to be $\propto N^{1/4}$ in the limit for $N\to\infty$. The treatment discussed so far concerns non-interacting systems; the extension to interacting systems is based on the argument that the theory of Fermi Liquids states that low lying excited states of a non-interacting system are in a one-to-one correspondence with non-interacting states, thus TA-GCE is likely to reduce finite size effects for interacting systems, as well. There are technical difficulties, though, that must be reported. One is that the wavefunction must be optimized at each value of $\theta$, which implies additional computational resources. Another problem is the treatment of net charges that are not balanced, but one can show that the average charge of the supercell (averaged over the BC) is exactly given by the density; so, in many cases, this is sufficient to get rid of the problem of a charged supercell. In general, reasonable solutions to such difficulties have been found, and the method has been applied to calculate several properties of interacting electrons [22, 103, 104, 105]. In the schematic picture of “system/environment”, with which I have analyzed methods and applications so far, the TA-GCE approach can be visualized as the construction of a generic “infinite” environment (reservoir of the Grand Canonical ensemble) via the introduction of the phase factor of the wavefunction that breaks the symmetry of the periodicity. This characteristic suggests that in principle, TA-GCE, in the not too distant future, may find use beyond its current technical utility of circumventing finite size effects. In this perspective, TA-GCE could be employed, for example, to describe the (small) QM part of A-QM/MM or el-QM-AdResS treated in the previous section, or for the QM description of the electrode in electrochemistry. Besides TA-GCE, other routes have been explored within the QMC scheme. An interesting recent example concerns the possibility of considering systems with a variable number of electrons for studying photoemission and inverse photoemission phenomena [109]. The method used is the Full Configuration Interaction QMC method (FCIQMC) developed in the group of Ali Alavi [106, 107, 108]. In essence, the method uses a Monte Carlo sampling of Slater determinants and calculates the real-time propagation of the wave function of the system as an electron is added or removed. Differently from TA-GCE, electrons are not introduced in a Grand Canonical fashion, but the method is specifically designed to compute electronic spectra upon sudden addition or removal of an electron (as in photoemission). In a qualitative explanation, this is equivalent to saying that the removed electron has been put in a stationary, non-interacting state somewhere at infinity. It must be added that nowadays it is even possible to run small QMC calculations on desktop computers; however, the computational costs of any QMC approach for large-scale systems routinely treated by, e.g., DFT, requires high-performance resources that may be beyond the means of research groups of a moderate size. Nevertheless, or, actually, because of such limitations, there is the need to significantly intensify the amount of theoretical work being done in this field.

5 Conclusions and Perspectives

The overall analysis of the various situations discussed in this paper leads me to the conclusion that what lies at the core of the problem that is simulating open, electronic systems, is not the lack of accurate theoretical tools, but rather the difficulty of merging the corresponding simulation techniques, specific to multiple aspects of the problem, in a physically consistent way. A proper, non-empirical merging is a necessary step for building simulation tools with predictive power, without the mandatory need for a case-by-case external validation. In the last decade we have witnessed the exponential increase of multiscale simulation techniques, mostly based on pragmatic empiricism for solving specific problems in the short term. The challenge of the next decade, in my view, should be the construction of multiscale techniques where the scale-coupling interfaces are based only on solid physical principles with corresponding well defined formulas for possible errors induced by approximations. Such a project requires methodological and theoretical work with, as much as possible, mathematical rigor, thus strengthening even further the concept that methodology (across fields and disciplines) is a self-standing, relevant field. On a broader scale, the project would require a combined effort of applied mathematicians, theoretical physicists and simulators (including physicists, chemists and materials scientists) with constructive exchanges is every direction. In this perspective, the aim of this paper is to offer to the molecular simulation and related communities a bird’s-eye view (of a theoretical physicist) on recent progress in the computational treatment of many-electron systems embedded in a fully-open larger environment. This paper shall not be intended as a comprehensive review but rather as a stimulating discussion for future developments, taking, as possible starting points, the methods, ideas and applications discussed in the various sections. These examples discussed have not been chosen because they are necessarily more relevant than others, but because together they sample the field in a rather uniform way. Their overall analysis leads to the conclusion that, as a QM technique, DFT clearly plays a key role, however quantum chemical methods (beyond Hartree-Fock) and QMC methods have been developed, at least at a conceptual level, up to the point of allowing calculations at a varying number of electrons that mimic the exchange of matter with the environment. Yet, their computational price is still prohibitive for treating systems such as molecules in solution or the metal surfaces associated with electrochemistry. The hope is that with the rapid evolution of the technology related to computational capabilities, the application of such methods to the above-mentioned systems may become possible. An intermediate step, that I can foresee, may be represented by the development of methods similar to the QM/MM approach with open boundaries, but with different levels of QM resolution (for a similar idea within the density matrix embedding approach see [110]). An example could be Grand Canonical QMC embedded in a quantum reservoir treated at the (Grand Canonical) DFT level, itself embedded in a reservoir of classical molecules. Analogous to the coupling criteria of el-QM-AdResS, one would require that the electronic and macroscopic chemical potentials are equal at the various interfaces. In summary, the future challenge lies in the construction of a solid theory of boundary conditions, at the interface of varying resolutions, in order to achieve a global physical consistency. A selection of starting points has been reported here and these are now ready to be developed even further.

Acknowledgments

This research has been funded by Deutsche Forschungsgemeinschaft (DFG) through the grant CRC 1114: “Scaling Cascades in Complex Systems”, project C01. I would like to thank Paul Ayers, Ferdinand Evers and Christian Krekeler for useful suggestions that have clarified relevant concepts discussed in this work. I am grateful to Ali Alavi, David Ceperley, Markus Holzmann and Ravishankar Sundararaman, for a critical reading of the sections of the manuscript related to their field of expertise. I wish to thank John Whittaker for a critical reading of the manuscript and for his precious suggestions.