Adiabatic motion and statistical mechanics via mass zero constrained dynamics
Sara Bonella, Alessandro Coretti, Rodolphe Vuilleumier, Giovanni Ciccotti
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] J a n Adiabatic motion and statistical mechanics via mass zero constrained dynamics
Sara Bonella, ∗ a Alessandro Coretti, a,b
Rodolphe Vuilleumier, c and Giovanni Ciccotti d,e,f a Centre Europ´een de Calcul Atomique et Mol´eculaire (CECAM),´Ecole Polytechnique F´ed´erale de Lausanne, Batochime, Avenue Forel 2, 1015 Lausanne,Switzerland. Fax: 0041 21 693 19 85; Tel: 0041 21 693 19 79; E-mail: sara.bonella@epfl.ch b Department of Mathematical Sciences, Politecnico di Torino,Corso Duca degli Abruzzi 24, I-10129 Torino, Italy. c PASTEUR, D´epartement de chimie, ´Ecole normale sup´erieure,PSL University, Sorbonne Universit´e, CNRS, 75005 Paris, France. d Institute for Applied Computing “Mauro Picone” (IAC), CNR Via dei Taurini 19, 00185 Rome, Italy. e Universit`a di Roma La Sapienza, Ple. A. Moro 5, 00185 Roma, Italy. f School of Physics, University College of Dublin UCD-Belfield, Dublin 4, Ireland.
In recent work [Coretti et al., The Journal of Chemical Physics , 2018, , 191102], a newalgorithm to solve numerically the dynamics of the shell model for polarization was presented.The approach, broadly applicable to systems involving adiabatically separated dynamical variables,employs constrained molecular dynamics to strictly enforce the condition that the force on thefast degrees of freedom, modeled as having zero mass, is null at each time step. The algorithm issymplectic and fully time reversible, and results in stable and efficient propagation. In this paperwe complete the discussion of the mechanics of mass zero constrained dynamics by showing howto adapt it to problems where the fast degrees of freedom must satisfy additional conditions. Thisextension includes, in particular, the important case of first principles molecular dynamics. We thenconsider the statistical mechanics of the mass zero constrained dynamical system demonstrating thatthe marginal probability sampled by the dynamics in the physical phase space recovers the form ofthe Born-Oppenheimer probability density.
PACS numbers:
I. INTRODUCTION
In this paper, we reconsider the simulation of adiabatically separated systems using the recently introduced masszero constrained molecular dynamics (MD) approach [1] by first extending the dynamics to the case of kinemati-cally constrained mass zero evolution and then determining the statistical mechanics associated to the dynamicalsystem. The inclusion of additional constraints enables to apply the method to new problems including first princi-ples molecular dynamics, while the statistical analysis defines the adiabatic probability for the physical system underconsideration.In adiabatically separated systems, the evolution of a set of physical degrees of freedom (dofs) is coupled to that ofa separate, auxiliary, set which is usually introduced to account for some specific physics or chemistry. The mass ofthe auxiliary dofs is assumed much smaller than the mass(es) associated to the physical variables, and the adiabaticregime is fully achieved in the limit of null auxiliary mass. Due to the time-scale separation of the motion of thetwo sets of dofs, in this regime the physical variables evolve in a potential computed at the minimum with respect tothe auxiliary variables. Depending on the specific system, additional (kinematic) conditions, such as orthonormalityor sum rules, may have to be imposed on the auxiliary dofs, affecting the minimum search. In MD simulations, thechoice of the algorithm to find the (free or conditioned) minimum is delicate. Typical schemes adopted so far includeiterative methods, such as conjugate-gradient [2, 3], extended Lagrangian schemes `a la Car-Parrinello in which asmall but finite mass is assigned to the auxiliary variables [4, 5] or, more recently, alternative ad hoc dynamics for theauxiliary variables like the so-called always stable predictor corrector approach [6, 7]. These methods, however, sufferfrom practical or conceptual limitations. Conjugate-gradient minimization is guaranteed to converge only in the caseof a quadratic function to be minimized, and, for the general minimization problems typically associated with realisticcondensed phase models, can be unstable [8] or expensive [9] to fully converge. In Car-Parrinello propagation, it canbe difficult or impossible to maintain adiabatic separation between the physical and auxiliary variables for sufficientlylong simulation times [2] and the algorithm requires a very small timestep to integrate accurately the dynamics ofthe auxiliary variables. The always stable predictor corrector scheme is only approximately time reversible [6, 7]leading to energy drifts that are usually quenched via a Berendsen thermostat (thus raising questions on the ensemblesampled by the dynamics), and it contains system dependent parameters that can only be determined by trial anderror. In spite of several interesting and successful applications of the approaches mentioned above, these limitationsjustify efforts to develop new algorithms for adiabatic dynamics.In recent work [1], we proposed to enforce the minimum energy condition using constrained molecular dynamics [10].In this approach, evolution equations for the physical and auxiliary degrees of freedom are derived under the constraintthat the derivative of the Hamiltonian with respect to the auxiliary variables is zero along the trajectory of the physicaldegrees of freedom. Consistent with the adiabatic prescription, this is achieved by writing the method of constraintsin the special case of mass zero constrained dofs [11]. The resulting dynamical system is then integrated combiningVerlet propagation with the SHAKE algorithm [12] to satisfy the constraint, thus ensuring that the overall numericalintegration of the system is symplectic and time reversible [13]. These formal properties guarantee exact energyconservation and stability of the evolution on the time-scales and with the time-step of standard MD for the physicaldegrees of freedom without the use of thermostats. Furthermore, the SHAKE algorithm, which requires iterationsto solve the constraint equation, usually converges very fast in particular in non-linear problems, enabling to fulfillthe minimum condition at the level of the numerical precision limit at an affordable cost. Importantly, SHAKE alsoprevents propagation of the error in each minimization along the trajectory. Proof of principle calculations [1] on aclassical model for polarizable systems have indeed shown that the mass zero constrained dynamics is an interestingand cost-effective alternative to currently adopted schemes.In this paper, we extend the approach to models where the fast degrees of freedom must satisfy additional conditions,such as specific normalization or orthogonality requirements. This development enables, in particular, to adapt themethod to the very interesting case of first principles molecular dynamics, but is relevant also for classical polarizablemodels, for example in connection with the global electroneutrality constraint in simulation of capacitors. Thesimilarities between classical polarizable models and first principles molecular dynamics have been long recognized,and algorithms have often migrated from one area to the other, a notable example being the seminal contributions ofSprik et al. [4]. Here we show that these similarities fully persist within the mass zero constrained approach. Havinggeneralized the mechanics of the system, we consider the question of the statistical mechanics associated with masszero constrained dynamics. As mentioned above, this approach entails an extended phase space in which both thephysical and the auxiliary variables are treated as dynamical degrees of freedom. In the following, we derive theprobability density sampled in this extended phase space, focusing in particular, on the marginal in the phase space ofthe physical variables, thus identifying unambiguously the full adiabatic limit for the probability associated to thesedegrees of freedom.The paper is organized as follows. For the reader’s convenience, in section II we begin by summarizing the derivationof the mass zero constrained dynamics for the case of classical polarization of shell-like models where only theminimization of the potential is required. We then show, using orbital-free Density Functional Theory [14, 15] as asimple but significative illustrative example, how to extend the dynamical system when an additional condition, in thiscase the conservation of the number of electrons, must be enforced. The derivation we present can be generalized easilyto situations, such as Kohn-Sham orbitals [16, 17], in which more than one additional condition must be satisfied.In section III, we discuss the statistical mechanics of mass zero constrained dynamics showing that the marginalprobability density in the physical degrees of freedom recovers the form usually assumed for the Born-Oppenheimerprobability.
II. ADIABATIC DYNAMICS VIA THE MASS ZERO CONSTRAINED EVOLUTIONA. Unconditioned minimization
To provide a specific example for the mass zero constrained approach in the case of simple minimization of the po-tential with respect to the fast variables, we shall illustrate our (general) formalism via the case of classical polarizablemodels. These models are crucial in the simulation of ionic systems of theoretical and technological interest such as,for example, energy storage devices (batteries and supercapacitors) [18–20]. An accurate description of polarizationis essential [21], for example, to compute with sufficient precision phonon dispersion curves in crystals or transportand structural properties in fluids, and to account for some apparently anomalous ordered structures observed insolids and melts without the need to invoke charge transfer effects.
First principles
MD can, in principle, includethese effects directly via the quantum treatment of the electronic charge distribution. However, since the sizes andtime-scales necessary to compute the properties of ionic systems are substantial, semi-empirical potentials are stillthe most convenient modeling tool in this field. In these potentials, appropriate sets of auxiliary dynamical variablesmimic changes in the electronic charge density. These variables might appear in the terms of the potential repre-senting short-range repulsion, Coulomb and dispersion interactions, and their coupling to the ionic dofs is controlledby parameters obtained from comparisons with experiments or with accurate ab initio calculations. One of the firstexamples of this kind was the shell model [22–24], which accounts for dipole polarization, followed by improvementstaking into account interactions due to quadrupoles [25] and changes in the ions size and shape [26, 27]. More recently,models for capacitors have been proposed that include the mutual polarization of the elements combining a multipoledescription of the electrolyte with the so called fluctuating charge model [4] for the electrodes.Let us indicate with R i ( i = 1 , ..., N ) the coordinates associated with the N physical variables in the system, andwith s α ( α = 1 , ..., M ) the M auxiliary degrees of freedom introduced to represent specific polarization effects. The s α variables may represent positions (as in the shell model) or different types of degrees of freedom (e.g. dipoles orquadrupoles, or charges) and their physical dimensions and number vary accordingly. The interaction potential of thesystem has the generic form V ( R, s ), where we have used the notation R = { R , ..., R N } and s = { s , ..., s M } . Thedynamics of the system is defined as follows. The auxiliary variables, having zero mass, adapt instantaneously to theposition of the physical degrees of freedom so as to satisfy σ α ( R, s ) ≡ ∂V ( R, s ) ∂s α = 0 α = 1 , . . . , M (1)In the following, we indicate with ˆ s the values of the auxiliary variables satisfying the conditions above. Note that,due to the dependence of the potential on the physical degrees of freedom, ˆ s = ˆ s ( R ). The evolution of the physicalvariables, of mass m i , is then given by m i ¨ R i = − (cid:20) ∂V ( R, s ) ∂R i (cid:21) s =ˆ s (2)The key idea of the mass zero constrained dynamics is to interpret eq. 1 as a set of M holonomic constraints andderive the evolution equations for the overall constrained system. This is done most conveniently by assigning at firsta finite mass µ to the auxiliary variables and considering the extended Lagrangian: L ( R, ˙ R, s, ˙ s ) = 12 N X i =1 m i ˙ R i + 12 M X α =1 µ ˙ s α − V ( R, s ) (3)together with eq. 1 to obtain the constrained evolution equations m i ¨ R i = − ∂V ( R, s ) ∂R i − M X β =1 λ β ∂σ β ( R, s ) ∂R i µ ¨ s α = − ∂V ( R, s ) ∂s α − M X β =1 λ β ∂σ β ( R, s ) ∂s α (4)where the λ β are the Lagrange multipliers associated with the constraints. The system above can be simplified byobserving that the gradient of the potential in the second equation is null due to the minimum condition. The equationfor the auxiliary degrees of freedom can then be rearranged as¨ s α = − M X β =1 λ β µ ∂σ β ( R, s ) ∂s α (5)We now consider the limit µ →
0. The equation above shows that, in order for the auxiliary variables to havefinite acceleration, the ratio γ α = lim µ → λ α µ must remain finite. This implies that the Lagrange multipliers λ α areproportional to µ . In the limit then, the constraint forces on the physical degrees of freedom vanish and we have m i ¨ R i = − ∂V ( R, s ) ∂R i ¨ s α = − M X β =1 γ β ∂σ β ( R, s ) ∂s α (6)Eq. 6 is the main result of ref. [1] and defines the mass zero constrained dynamics. It shows that, consistent witheq. 2, the evolution of the physical variables does not depend directly on the constraints. Due to the constraints, onthe other hand, the (free) s variables satisfy the minimum condition. This implies that this condition is automaticallysatisfied also in the first equation which is then equivalent to eq. 2. The system above is an exact, classical, evolutionfor all degrees of freedom that rigorously enforces the zero mass limit, and therefore leads to full adiabatic separation,for the auxiliary variables. The numerical integration of the first equation can be performed with any standard MDalgorithms (e.g. Verlet) with a time-step determined only by the explicit physical force term. In addition, at eachtime-step, the Lagrange multipliers γ α , that appear as unknown, time-dependent parameters in the dynamical system,must be determined. This is done enforcing the constraint, σ α ( R ( t + dt ) , s ( t + dt )) = 0, at the position predicted by theMD algorithm as described in [10, 12]. In this way propagation of the error is prevented. In current implementationsof the approach, the constraint is satisfied via the SHAKE iterative algorithm, which was proven to be symplecticand time reversible [13, 28]. B. Conditioned minimization
So far, we have assumed that the only conditions to be satisfied by the auxiliary dofs were those expressed byeq. 1. While this is often the case, e.g. the shell model or classical treatment of multipoles, there are importantproblems in which additional requirements are placed on these degrees of freedom. The most interesting exampleis probably given by first principles
MD, that we shall therefore adopt hereafter to illustrate how the mass zeroconstrained dynamics can be extended to include generic additional conditions. In first principles
MD, classicalpropagation of the nuclei is combined with a quantum evaluation of the forces, obtained on-the-fly from accurateelectronic structure calculations. In this approach, adiabatic separation of the nuclear and electronic dofs is invokedto claim that, at each nuclear configuration, the electronic subsystem has relaxed in its ground state. The groundstate energy is found by minimizing the quantum expectation value of the electronic Hamiltonian with respect tothe (normalized) electronic state. This energy is a functional of the electronic dofs that depends parametricallyon the nuclear positions, thus providing the force on the nuclei via the appropriate derivatives. Currently, the bestcompromise between cost and accuracy in determining the ground state energy is offered by Density Functional Theory(DFT). DFT-
First principles
MD has indeed seen an impressive growth since the mid 1980s and is adopted in areasranging from materials modeling to drug design [17, 29, 30]. The success of available methods, however, still leavesroom for methodological and algorithmic improvements to secure stable, unbiased, and efficient long-time evolutionas witnessed by several recent efforts [31–40]. To demonstrate how mass zero constrained MD can be extended tothis problem, we specify our discussion to the so-called orbital-free DFT scheme [14, 15] in which the electronic stateis represented directly via the electronic density, n ( r ). We have chosen this representation because it permits toillustrate the key ingredients of our approach with minimal notational complexity since it requires to include only oneadditional condition, i.e. conservation of the number or electrons, to the minimization. The extension to the case ofmultiple conditions is straightforward and will be discussed at the end of the section. The electronic density is a fieldin Cartesian space, but in the following we shall restrict our considerations to an appropriate discretization of thisquantity. This discretization must, of course, always be enforced in numerical implementations and can be realizedvia a decomposition on a basis set or a representation on a grid. Below, we shall consider a grid decomposition butthe discussion can be similarly performed in a basis. Let us then indicate with r k the points on the grid and use thenotation n ( r k ) ≡ n k ( k = 1 , . . . , M ) for the discretized density. These n k variables play the role of the auxiliary dofs s α of the previous section. The Born-Oppenheimer evolution equations for the nuclear degrees of freedom are givenby m i ¨ R i = − (cid:20) ∂E ( R, n ) ∂R i (cid:21) n =ˆ n (7)where E ( R, n ) is the ground state expectation value of the electronic Hamiltonian, h ψ | H el | ψ i , expressed as afunction of the discretized density variables n = { n , . . . , n M } . In the equation above, we have indicated with ˆ n thevalues of these variables that minimize E ( R, n ) subject to the condition f ( n ) = M X k =1 n k δV − N el = 0 (8)which expresses, in discrete form, the fact that the integral of the electronic density must be equal to the total numberof electrons, N el , in the system ( δV is the elementary volume). Eq. 8 implies that not all variations of the n k areindependent, and this must be accounted for in the statement of the minimum condition. This can be done mostconveniently via the Lagrange multiplier method by introducing the auxiliary function E ( R, n, ν ) ≡ E ( R, n ) − νf ( n ) (9)where ν is the Lagrange multiplier associated to the condition above. The solution ˆ n for n satisfying eq. 8 is thengiven by the stationary point (ˆ n, ˆ ν ) of E ( R, n, ν ).[41, 42] This leads to the M + 1 conditions σ k ( R, n, ν ) = ∂ E ( R, n, ν ) ∂n k = 0 ( k = 1 , . . . , M ) σ M +1 ( R, n, ν ) = ∂ E ( R, n, ν ) ∂ν = 0 . (10)Note that, as usual in the Lagrange multiplier scheme, the last equation above is in fact the additional conditioneq. 8 but now obtained as a result of an optimization problem in the space that includes the Lagrange multiplier ν as a variable. Within this framework, we can proceed in analogy with the previous section by defining the followingextended Lagrangian, which includes also ν as an auxiliary variable, L ( R, ˙ R, n, ˙ n, ν, ˙ ν ) = 12 N X i =1 m i ˙ R i + 12 µ n M X k =1 ˙ n k + 12 µ ν ˙ ν − E ( R, n, ν ) , (11)and interpreting eq. 10 as M + 1 holonomic constraints on the system. In the equation above, we have introduced, fordimensional consistency, two finite masses µ n and µ ν related to the auxiliary variables associated with the discretizeddensity and with the Lagrange multiplier ν , respectively. Proceeding as in the previous section, and imposing thatthe masses of both sets of auxiliary variables tend to zero, the mass zero constrained dynamical system is given by m i ¨ R i = − ∂E ( R, n ) ∂R i , ¨ n i = − M +1 X α =1 γ α ∂σ α ( R, n, ν ) ∂n i , ¨ ν = − M +1 X α =1 γ α ∂σ α ( R, n, ν ) ∂ν (12)where, in the first equation, we have used the fact that ∂ E ( R,n,ν ) ∂R i = ∂E ( R,n ) ∂R i , while the first derivatives of E in thesecond and third line vanish due to the stationary point conditions. In this extended system, we now recover asituation similar to that of the polarizable models and the numerical integration of the equations above can again beperformed by combining standard algorithms such as Verlet with SHAKE to enforce the constraints. In this respect,note that the potential energy term E ( R, n, ν ) is linear in ν with the consequence that the second derivative of E iszero in this direction. Thus the matrix of second derivatives at the stationary point (ˆ n, ˆ ν ) is non-positive definite and(ˆ n, ˆ ν ) is a saddle point.[41, 42] However, this does not hinder the procedure since SHAKE only requires non-singularitybut not positive definiteness of the matrix of second-derivatives.The framework described above is quite general and can be easily extended to the case of G additional conditions, f j ( s ) (with j = 1 , . . . , G ), for the auxiliary degrees of freedom. This can be done by introducing, and subsequentlytreating as an additional auxiliary variable, one Lagrange multiplier, ν j , per additional condition and modifying theLagrangian by subtracting the sum P Gj =1 ν j f j ( s ) to the energy functional. This procedure enables, for example, toapply the mass zero dynamics procedure to Kohn-Sham orbital based first principle MD thus further illustrating theoverall interest of this approach.
III. STATISTICAL MECHANICS OF THE MASS ZERO CONSTRAINED EVOLUTION
We now investigate the probability density sampled by the dynamical system eq. 6 . This is an interesting questionbecause the use of constraints may induce a non-trivial metrics in the phase space of the system [43] and appropriatereweighting of averaging might be necessary when computing statistical properties [44] in the physical phase space.To set the stage for the discussion, and proceeding in analogy with [45], it is useful to introduce the change of variables(we avoid considering also the ν variable(s) for notational simplicity) R i ρ i = R i s α σ α = σ α ( R, s ) (13)Observables computed in the new variables will be denoted in calligraphic font. In the following, we shall also use thenotation υ = ( ρ, σ ) where ρ = { ρ , ..., ρ N } and σ = { σ , ..., σ M } . The Lagrangian in the new set of coordinates isobtained from eq. 3 as L ( υ, ˙ υ ) = L (cid:0) R ( υ ) , s ( υ ) , ˙ R ( υ, ˙ υ ) , ˙ s ( υ, ˙ υ ) (cid:1) . The Hamiltonian of the system in the new coordinatesis then derived via a standard Legendre transform and can be written as H ( υ, π υ ) = K ( υ, π υ ) + V ( υ ) ≡ π υ T M − ( υ ) π υ + V ( υ ) (14)where K ( υ, π υ ) is the kinetic energy in the new variables, and the momentum π υ is given by π υk = ∂ L ( υ, ˙ υ ) ∂ ˙ υ k = ( π ρi = ∂ L ( ρ, ˙ ρ,σ, ˙ σ ) ∂ ˙ ρ i if k = 1 , . . . , Nπ σα = ∂ L ( ρ, ˙ ρ,σ, ˙ σ ) ∂ ˙ σ α if k = 3 N + 1 , . . . , N + M (15)(as always, i = 1 , . . . , N and α = 1 , . . . , M ) and M − ( υ ) is the inverse of the metric matrix associated with the newvariables M kk ′ = N X i =1 m i ∂R i ∂υ k ∂R i ∂υ k ′ + M X α =1 µ ∂s α ∂υ k ∂s α ∂υ k ′ (16)The metric matrix and its inverse can be conveniently expressed in block form as M = (cid:20) A BB T Γ (cid:21) with inverse M − = (cid:20) ∆ EE T Z (cid:21) (17)where A and ∆ are 3 N × N matrices, B and E are 3 N × M matrices, and Γ and Z are M × M matrices, whoseexpressions are given in the Appendix. The Hamiltonian in eq. 14 can be evaluated explicitly on the constrainedhypersurface in the full phase space. To that end, it is important to recognize that, since the set of constraints ineq. 1 have to be fulfilled at all times, the additional conditions ˙ σ = 0 hold. The definition ˙ υ = M − π υ then impliesthat, when the constraints are imposed, the momenta π σ must satisfyˆ π σ = − e Z − e E T π ρ (18)where the tildes indicate that all matrices are evaluated at σ = 0. The constrained motion then occurs on thehypersurface σ = 0, π σ = ˆ π σ and we have H ( ρ, σ = 0 , π ρ , π σ = ˆ π σ )= 12 ( π ρ , ˆ π σ ) (cid:20) e ∆ e E e E T e Z (cid:21) (cid:18) π ρ ˆ π σ (cid:19) + V ( ρ, σ = 0) (19)The Hamiltonian can be written in a more suitable form by observing that the block expression of the product MM − = imposes e A e ∆ + e B e E T = and e A e E + e B e Z = . These two relationships, in turn, imply e ∆ − e E e Z − e E T = e A − sothat H ( ρ, σ = 0 , π ρ , π σ = ˆ π σ )= 12 [ π ρ T ( e ∆ − e E e Z − e E T ) π ρ ] + V ( ρ, σ = 0)= 12 π ρ T e A − π ρ + V ( ρ, σ = 0) (20)The microcanonical partition function in the full phase space ( ρ, σ, π ρ , π σ ) can now be written as Z = Z d N ρ d N π ρ d M σ d M π σ δ M ( σ ) δ M ( π σ − ˆ π σ ) ×× δ ( H ( ρ, σ, π ρ , π σ ) − E ) (21)while the expression of the microcanonical average of observable O is given by hOi = 1 Z Z d N ρ d N π ρ d M σ d M π σ δ M ( σ ) δ M ( π σ − ˆ π σ ) × O ( ρ, σ, π ρ , π σ ) δ ( H ( ρ, σ, π ρ , π σ ) − E ) (22)The expression above can be usefully simplified by first integrating over the π σ variables hOi = 1 Z ′ Z d N ρ d N π ρ d M σδ M ( σ ) × O ( ρ, σ, π ρ , ˆ π σ ) δ ( H ( ρ, σ, π ρ , ˆ π σ ) − E ) (23)and then performing the change of variables σ α s α , ρ i R i to obtain at first hOi = 1 Z ′ Z d N R d N π ρ d M s | J ( R ) | δ M ( σ ( R, s )) × O ( R, σ ( R, s ) , π ρ , ˆ π σ ) δ ( H ( R, σ ( R, s ) , π ρ , ˆ π σ ) − E ) (24)where | J ( R ) | is the Jacobian of the coordinate transformation, which reduces to det (cid:2) ∂σ∂s (cid:3) . Then, making the depen-dence on s of the delta explicit, we get hOi = 1 Z ′ Z d N R d N π ρ d M s | J ( R ) || J ( R ) | − δ M ( s − ˆ s ( R )) × O ( R, s, π ρ , ˆ π σ ) δ ( H ( R, s, π ρ , ˆ π σ ) − E ) (25)where we have O ( R, s, π ρ , ˆ π σ ) = O ( R, σ ( R, s ) , π ρ , ˆ π σ ). In this last equality we have used the properties of the deltaof a vector function of the integration variable to express the constraint condition directly as a function of the s , withˆ s ( R ) such that σ ( R, ˆ s ) = 0 (we assume, as commonly done in the Born-Oppenheimer framework that this expressionhas, for any R , a single root). Finally, performing the integral over the s variables, and noting that the product ofJacobians in the integrand simplifies, we obtain hOi = 1 Z ′ Z d N R d N π ρ O ( R, ˆ s, π ρ , ˆ π σ ) δ ( H ( R, ˆ s, π ρ , ˆ π σ ) − E ) (26)where (see eq. 20) H ( R, ˆ s, π ρ , ˆ π σ ) = 12 π ρ T e A − π ρ + V ( R, ˆ s ) (27)Let us now consider the limit µ → e A = D + µ e R = D h + µ D − e R i (28)where D jj ′ = m j δ jj ′ and e R jj ′ = P Mα =1 ∂s α ∂ρ j ∂s α ∂ρ j ′ so that e A − = h + µ D − e R i − D − (29)Furthermore, the relation π υ = M ˙ υ implies, π ρ = A ˙ ρ + B ˙ σ giving, on the constrained hypersurface π ρ = e A ˙ ρ = e A ˙ R (30)where in the last equality we used ˙ ρ = ˙ R , as implied by eq. 13. Finally, from MM − = , the identities e A e E + e B e Z = and e B T e E + e Γ e Z = follow for the involved submatrices. Using these identities, we obtain e Z − = e Γ − e B T e A − e B (31)so that ˆ π σ = − e Z − e E T π ρ = he Γ − e B T e A − e B i e E T π ρ (32)We can now take the limit µ → e Γ and e B are proportional to µ and therefore vanish in thelimit, while lim µ → e A = e D . In the zero auxiliary mass limit then, ˆ π σ = 0 and the Hamiltonian of the system becomes H ( R, ˆ s, π R , ˆ π σ = 0) = 12 π R T D − π R + V ( R, ˆ s ) (33)with D − = m j δ jj ′ and where we have used the fact that, see eq. 30, in the null auxiliary mass limit, π ρ = D ˙ R = π R .Substituting in the expression for the average, we obtain hOi = 1 Z ′ Z d N R d N π R O ( R, ˆ s, π R , ˆ π σ = 0) δ ( H ( R, ˆ s, π R , ˆ π σ = 0) − E ) (34)The result above implicitly defines the microcanonical marginal probability in the physical phase space in the fulladiabatic limit. Interestingly, this definition is in agreement with the form usually assumed for the Born-Oppenheimerprobability. The derivation presented here, however, indicates that the dynamical systems rigorously samples thisdensity only in the full µ → IV. CONCLUSIONS
In this paper, we have shown how the mass zero constrained dynamics can be adapted to models that imposeadditional conditions on the dynamics of the auxiliary variables. The interesting case of Born-Oppenheimer firstprinciples molecular dynamics for the so-called orbital free approach to DFT was used to illustrate in some detailthe extension of the approach for these systems when a single additional condition must be fulfilled. The procedurecan be easily extended to the case of multiple additional conditions, thus enabling to tackle also Kohn-Sham DFTpropagation. We have also analyzed the statistical mechanics associated to the mass zero constrained dynamics foradiabatically separated systems, deriving the fully adiabatic marginal probability in the physical phase space. Thisprobability coincides with the one typically assumed in Born-Oppenheimer sampling. The discussion presented herefocused on the microcanonical ensemble, but the generalization to other ensembles is straightforward.
Appendix
Here we provide the explicit expressions of the block matrices appearing in eq. 17. Remembering that ρ = R and σ = σ ( R, s ), we have A jj ′ = N X i =1 m i ∂R i ∂ρ j ∂R i ∂ρ j ′ + M X α =1 µ ∂s α ∂ρ j ∂s α ∂ρ j ′ = m j δ jj ′ + M X α =1 µ ∂s α ∂ρ j ∂s α ∂ρ j ′ for j, j ′ = 1 , . . . , N B jβ = N X i =1 m i ∂R i ∂ρ j ∂R i ∂σ β + M X α =1 µ ∂s α ∂ρ j ∂s α ∂σ β = M X α =1 µ ∂s α ∂ρ j ∂s α ∂σ β for j = 1 , . . . , N , β = 1 , . . . , M Γ ββ ′ = N X i =1 m i ∂R i ∂σ β ∂R i ∂σ β ′ + M X α =1 µ ∂s α ∂σ β ∂s α ∂σ β ′ = M X α =1 µ ∂s α ∂σ β ∂s α ∂σ β ′ for β, β ′ = 1 , . . . , M (35)and ∆ jj ′ = N X i =1 m i ∂ρ j ∂R i ∂ρ j ′ ∂R i for j, j ′ = 1 , . . . , N E jβ = − N X i =1 m i ∂ρ j ∂R i ∂σ β ∂R i for j = 1 , . . . , N , β = 1 , . . . , M Z ββ ′ = N X i =1 m i ∂σ β ∂R i ∂σ β ′ ∂R i + M X α =1 µ ∂σ β ∂s α ∂σ β ′ ∂s α for β, β ′ = 1 , . . . , M (36) Conflicts of interest
There are no conflicts to declare.
Acknowledgements
The authors are grateful to M. Salanne and B. Rothenberg for useful exchanges on the statistics of Born-Oppenheimer dynamics. [1] A. Coretti, S. Bonella and G. Ciccotti,
The Journal of Chemical Physics , 2018, , 191102.[2] A. Aguado, L. Bernasconi and P. A. Madden,
The Journal of Chemical Physics , 2003, , 5704–5717.[3] S. Jahn, P. A. Madden and M. Wilson,
Phys. Rev. B , 2004, , 020106.[4] M. Sprik and M. L. Klein, The Journal of Chemical Physics , 1988, , 7556–7560.[5] M. Wilson and P. A. Madden, Journal of Physics: Condensed Matter , 1993, , 2687.[6] J. Kolafa, Journal of Computational Chemistry , 2004, , 335–342.[7] J. Genzer and J. Kolafa, Journal of Molecular Liquids , 2004, , 63 – 72.[8] F. Pacaud, M. Salanne, T. Charpentier, L. Cormier and J.-M. Delaye,
Journal of Non-Crystalline Solids , 2018, , 371– 379.[9] M. Pounds, S. Tazi, M. Salanne and P. A. Madden,
Journal of Physics: Condensed Matter , 2009, , 424109.[10] G. Ciccotti and J.-P. Ryckaert, Computer Physics Reports , 1986, , 346 – 392.[11] J.-P. Ryckaert, A. Bellemans and G. Ciccotti, Molecular Physics , 1981, , 979–996.[12] J.-P. Ryckaert, G. Ciccotti and H. J. Berendsen, Journal of Computational Physics , 1977, , 327 – 341.[13] B. J. Leimkuhler and R. D. Skeel, Journal of Computational Physics , 1994, , 117 – 125.[14] M. Pearson, E. Smargiassi and P. A. Madden,
Journal of Physics: Condensed Matter , 1993, , 3221–3240.[15] V. L. Lign`eres and E. A. Carter, in An Introduction to Orbital-Free Density Functional Theory , ed. S. Yip, SpringerNetherlands, Dordrecht, 2005, pp. 137–148.[16] W. Kohn and L. J. Sham,
Physical Review , 1965, , A1133–A1138.[17] D. Marx and J. Hutter,
Ab initio molecular dynamics: Basic theory and advanced methods , Cambridge University Press,Cambridge, 2012.[18] P. Simon and Y. Gogotsi,
Nature Materials , 2008, , 845 EP –.[19] M. Armand and J. M. Tarascon, Nature , 2008, , 652 EP –.[20] F. Beguin, V. Presser, A. Balducci and E. Frackowiak,
Advanced Materials , 2014, , 2219–2251.[21] A. M. Stoneham and J. H. Harding, Annual Review of Physical Chemistry , 1986, , 53–80.[22] B. G. Dick and A. W. Overhauser, Physical Review , 1958, , 90–103.[23] Jacucci, G., Klein, M.L. and McDonald, I.R.,
J. Physique Lett. , 1975, , 97–100.[24] G. Jacucci, I. R. McDonald and A. Rahman, Physical Review A , 1976, , 1581–1592.[25] M. Wilson, P. A. Madden and B. J. Costa-Cabral, The Journal of Physical Chemistry , 1996, , 1227–1237.[26] M. Wilson, P. A. Madden, N. C. Pyper and J. H. Harding,
The Journal of Chemical Physics , 1996, , 8068–8081.[27] A. J. Rowley, P. ¨Jemmer, M. Wilson and P. A. Madden,
The Journal of Chemical Physics , 1998, , 10209–10219.[28] B. Leimkuhler and S. Reich,
Simulating Hamiltonian dynamics , Cambridge University Press, Cambridge, UK New York,2004.[29] A. Cavalli, P. Carloni and M. Recanatini,
Chemical Reviews , 2006, , 3497–3519.[30] P. R. Briddon and M. J. Rayson,
Physica Status Solidi (b) , 2011, , 1309–1318.[31] A. M. N. Niklasson, C. J. Tymczak and M. Challacombe,
Physical Review Letters , 2006, , 123001.[32] A. M. N. Niklasson, C. J. Tymczak and M. Challacombe, The Journal of Chemical Physics , 2007, , 144103.[33] T. D. K¨uhne, M. Krack, F. R. Mohamed and M. Parrinello,
Physical Review Letters , 2007, , 066401.[34] A. M. N. Niklasson, Physical Review Letters , 2008, , 123004.[35] A. M. N. Niklasson, P. Steneteg, A. Odell, N. Bock, M. Challacombe, C. J. Tymczak, E. Holmstr¨om, G. Zheng andV. Weber,
The Journal of Chemical Physics , 2009, , 214109.[36] A. M. N. Niklasson and M. J. Cawkwell,
Physical Review B , 2012, , 174308.[37] L. Lin, J. Lu and S. Shao, Entropy , 2014, , 110–137.[38] T. D. K¨uhne, Wiley Interdisciplinary Reviews: Computational Molecular Science , 2014, , 391–406.[39] A. M. N. Niklasson, The Journal of Chemical Physics , 2017, , 054103.[40] A. Albaugh, T. Head-Gordon and A. M. N. Niklasson,
Journal of Chemical Theory and Computation , 2018, , 499–511.[41] C. Lanczos, The Variational Principles of Mechanics , Dover books in physics, New York, USA, 1970.[42] G. Allaire and A. Craig,
Numerical Analysis and Optimization: An Introduction to Mathematical Modelling and NumericalSimulation , Oxford University Press, Incorporated, Oxford, United Kingdom, 2007.[43] M. E. Tuckerman, Y. Liu, G. Ciccotti and G. J. Martyna,
The Journal of Chemical Physics , 2001, , 1678–1702. [44] G. Ciccotti, R. Kapral and E. Vanden-Eijnden, ChemPhysChem , 2005, , 1809–1814.[45] G. Ciccotti and M. Ferrario, Computation , 2018,6