Mass-Zero constrained dynamics and statistics for the shell model in magnetic field

In several domains of physics, including first principle simulations and classical models for polarizable systems, the minimization of an energy function with respect to a set of auxiliary variables must be performed to define the dynamics of physical degrees of freedom. In this paper, we discuss a recent algorithm proposed to efficiently and rigorously simulate this type of systems: the Mass-Zero (MaZe) Constrained Dynamics. In MaZe the minimum condition is imposed as a constraint on the auxiliary variables treated as degrees of freedom of zero inertia driven by the physical system. The method is formulated in the Lagrangian framework, enabling the properties of the approach to emerge naturally from a fully consistent dynamical and statistical viewpoint. We begin by presenting MaZe for typical minimization problems where the imposed constraints are holonomic and summarizing its key formal properties, notably the exact Born-Oppenheimer dynamics followed by the physical variables and the exact sampling of the corresponding physical probability density. We then generalize the approach to the case of conditions on the auxiliary variables that linearly involve their velocities. Such conditions occur, for example, when describing systems in external magnetic field and they require to adapt MaZe to integrate semiholonomic constraints. The new development is presented in the second part of this paper and illustrated via a proof-of-principle calculation of the charge transport properties of a simple classical polarizable model of NaCl.


Introduction
In this paper, we discuss, focusing on recent developments, the Mass-Zero (MaZe) constrained dynamics and further extend it to a simple but interesting model of classical polarizable systems in constant external magnetic field. MaZe is a general simulation approach to study the motion of a set of physical degrees of freedom (dofs) whose evolution depends on parameters subject to given conditions. The method considers an extended system in which the parameters appear as (auxiliary) dynamical variables together with the original dofs and the conditions are interpreted as constraints. The coupled evolution equations for the overall constrained system are then conveniently obtained in the Lagrangian formalism. From these, the original parametric dynamics for the physical dofs is rigorously recovered by taking the limit of zero mass for the auxiliary variables. In practical implementations, use of a e-mail: alessandro.coretti@polito.it (corresponding author) b e-mail: sara.bonella@epfl.ch the SHAKE algorithm enables symplectic and efficient numerical integration of the extended dynamical system.
The method of mass-zero constraints was originally introduced in the early 1980s to study the rotationaltranslational coupling in diatomic molecules [1]. Recently, it has undergone a new set of developments when adiabatic systems were identified as an important area where MaZe dynamics can provide an original formal approach and an effective integration algorithm [2][3][4]. In adiabatic systems, the substantial timescale separation of the motion of two sets of interacting degrees of freedom justifies adopting the Born-Oppenheimer approximation for the evolution. The timescale gap is typically due to the disparate masses of the two sets and full adiabatic separation is achieved in the limit of zero mass for the fast dofs. In the context of nuclear and electronic motion, for which the Born-Oppenheimer approximation was originally introduced, adiabaticity also requires the existence of a finite gap between the HOMO and LUMO electronic states. In the full adiabatic regime, evolution equations, typically of classical form, are given for the slow degrees of freedom. The forces on the slow variables, however, depend parametrically on the values of the fast variables. These values are obtained, for each configuration of the slow degrees of freedom along the trajectory, enforcing the condition that the interaction potential (a function of both sets) is at a minimum with respect to the fast dofs. In Molecular Dynamics (MD), a typical example of adiabatic dynamics is the evolution of ionic (slow) and electronic (fast) degrees of freedom in first principle calculations based on Kohn-Sham or orbital-free Density Functional Theory. Another important example is given by classical models of polarization in which electrons do not appear directly, but the dynamical system is extended to include sets of classical auxiliary variables of null mass mimicking different polarization effects.
Current methods adopted for the MD simulation of such systems combine standard propagation schemes for the slow variables-that we shall indicate as the ions for simplicity-with algorithms to find, or approximate, the minimum of the potential with respect to the fast dofs at each ionic configuration. Depending on the specific system, additional conditions, such as orthonormality or sum rules, may be imposed on the fast dofs, affecting the minimum search. Traditional schemes for finding the minimum include iterative methods, notably the conjugate-gradient approach [5,6] adopted in Born-Oppenheimer MD. In Car-Parrinello MD, on the other hand, an extended system in which the auxiliary dofs are treated as dynamical variables with a small mass is defined [7,8]. In this scheme, the minimum condition is approximately tracked, with a precision that improves with smaller mass for the auxiliary dofs [9][10][11][12], via the dynamics itself thus avoiding iterations. More recently, alternative ad hoc dynamics for the auxiliary variables like the so-called always stable predictorcorrector approach [13,14] have also been employed. All these methods, however, suffer from practical or conceptual limitations. Conjugate-gradient minimization is guaranteed to converge only in the case of a quadratic function to be minimized, and, for the general minimization problems typically associated with realistic condensed-phase models, can be unstable [15] or expensive [16] to fully converge. Incomplete convergence of the iterative minimization in Born-Oppenheimer MD, also known as the self-consistent-field optimization, has been shown to cause energy transfer between the slow and fast dofs [17,18], leading to energy drift in the ionic propagation and hindering access to long simulation timescales. Energy transfer, and the consequent violation of the adiabatic separation in the system, affects also Car-Parrinello propagation [5] due to the finite ratio of the masses associated with the fast and slow dofs. This pathology is often mitigated via thermostats that, however, may affect or, in the worst cases, compromise correct statistical sampling. Furthermore, the Car-Parrinello algorithm requires a very small timestep to integrate accurately the dynamics of the fast variables. The always stable predictor-corrector scheme is only approximately time-reversible [13,14] leading again to energy drifts that are usually quenched via a Berendsen thermostat (thus raising questions on the ensemble sampled by the dynamics), and it contains system-dependent parameters that can only be determined by trial and error. An alternative and quite successful scheme, that combines an extended Lagrangian approach with efficient self-consistent minimization, was recently proposed to address these shortcomings in the so-called extended Lagrangian Born-Oppenheimer approach [19,20]. The method is formally fully timereversible and has been used in several interesting applications. The practical integration of the evolution equations for the auxiliary variables, however, requires, in particular for longer simulations, to dampen propagation of numerical noise that causes a divergence of the auxiliary dofs from the exact ground-state density [20,21]. The modified Verlet algorithm used for this purpose breaks full-time reversibility of the propagation.
The MaZe approach avoids many of the difficulties described above. Adopting the framework of constrained MD [22], in combination with the SHAKE algorithm [23], MaZe enables to derive and numerically solve classical evolution equations for an extended system that rigorously enforces exact adiabatic propagation. Exploiting and adapting the formalism of Lagrange multipliers, the method can incorporate easily additional constraints that depend only on the auxiliary variables, such as electroneutrality in classical polarizable models or orthonormality in first principles calculations. The MaZe dynamical system is solved via a fully symplectic, time-reversible algorithm that guarantees stability of the evolution on the same timescale and with the same timestep size of standard MD for the physical dofs. The algorithm prevents, by construction, propagation of the error when imposing the constraints. Furthermore, the approach avoids standard self-consistent cycles for the minimization along the dynamics and uses iterations only to solve the equation of constraints, a process that usually has fast convergence, in particular in nonlinear problems, enabling to reach essentially the numerical precision limit at an affordable cost [2][3][4]. Importantly, rigorous sampling of the target probability for the physical dofs (which coincides with the one usually associated to Born-Oppenheimer dynamics) is also guaranteed [3] ensuring that not only the adiabatic evolution, but also the exact statistical mechanics of the system is obtained. In the context of adiabatic systems, MaZe has been used to simulate simple classical polarizable systems [2], and more recently for state-of-the-art classical modeling of electrode charges in electrochemical systems [4] (the method has also been implemented in MetalWalls [24], a high-performance community software in this area). The generalization to first principle MD based on orbital-free density functional theory was also derived and tested with very good results [3].
In the following, we first discuss the key aspects of MaZe summarizing recent work on the dynamical and statistical properties of the approach. For simplicity, the formalism will be presented using as reference application the adiabatic evolution of classical polarizable models. This choice is motivated also by the new MaZe development presented in the second part of the paper: the generalization of the approach for adiabatic propagation of classical polarizable models in external magnetic field. Although the model system considered in this work has some limitations (see detailed description in the following), this generalization is motivated by the problem of simulating ionic charge transport in systems subject to a magnetic field, with particular focus on the ionic Hall effect [25][26][27], an area that has attracted considerable interest to investigate the properties of superionic conductors [28,29] and, more recently, to enhance the capacitance of batteries [30]; including the magnetic field in classical polarizable models requires, even with the simplified model adopted here, some nontrivial adaptations of the MaZe approach. First, due to the Lorentz force, the condition of null force on the fast degrees of freedom involves both the coordinates and the velocities of these particles, thus leading to a set of nonholonomic constraints associated with the MaZe dynamics. These constraints are, however, linear in the velocities, enabling to adapt the approach via a relatively standard generalization of the Lagrangian equations of motion for the system. Second, again due to the dependence of the constraint on the velocities, the standard SHAKE algorithm-not the idea at the basis of the approach-cannot be directly applied. While some methods for the numerical integration of systems subject to nonholonomic constraints exist [31,32], they are not as consolidated as SHAKE and often rely on nonsymplectic algorithms. In this work, we then propose an appropriate integration algorithm and test its properties.
The paper is organized as follows. In Sect. 2, the derivation of the MaZe dynamical system and the proof of exact sampling of the Born-Oppenheimer probability density for the slow degrees of freedom are summarized. The presentation will be self-contained, also with the support of Appendix A, but we limit the proofs to their key steps, referring to previous work for more details. Our focus, in fact, will be to illustrate the most interesting formal and practical aspects of the approach. In Sect. 3, we then introduce the generalization of the MaZe framework to the case of classical polarizable simulations in constant magnetic field, together with the new algorithm for the solution of the nonholonomic constrained dynamics. We consider, for convenience and given the exploratory purposes of these first MaZe developments in the presence of magnetic field, the simplest model of classical polarization, i.e., the shell model [33], and neglect possible magnetization effects. Section 4 reports an illustrative calculation investigating the combined effect of polarization and of the magnetic field on ionic charge transport properties in liquid NaCl.

MaZe dynamics and statistical mechanics
As mentioned above, to simplify the presentation and tackle an interesting physical case, we illustrate the MaZe approach via its application to classical polarizable models. These models are commonly employed to simulate systems of theoretical and technological interest such as devices for electrochemical energy storage [34][35][36] in which large sizes and long timescale prevent direct calculation of polarization effects via the quantum treatment of the electronic density in first principles MD. In this area then, polarization effects are described by constructing empirical potentials that include sets of auxiliary variables that mimic changes in the electronic density. Because they represent electronic properties, these auxiliary variables are assumed to adapt instantaneously to the ionic configuration in the spirit of the Born-Oppenheimer approximation and are assigned a null mass. An early example of this type of models is the shell model [33,37,38], which accounts for dipole polarization. Potentials that take into account interactions due to quadrupoles [39] and changes in the ions size and shape [40,41] have also been introduced. More recently, models for capacitors have been proposed that include the mutual polarization of the elements combining a multipole description of the electrolyte with the so-called fluctuating charge model [7,42] for the electrodes.
To be more specific, let us indicate with R ∈ R 3N the Cartesian coordinates of the N physical dofs in the system (ions), and with s ∈ R M the M adiabatically separated auxiliary variables. Depending on the specific polarization model, the s variables may represent positions (as in the shell model) or different types of degrees of freedom (e.g., dipoles or quadrupoles, or charges) and their physical dimensions and number vary accordingly. The two sets of variables interact via the potential V (R, s). The adiabatic dynamics of the system is obtained by imposing that the auxiliary variables satisfy, for all values of R along the trajectory, the condition ∂V (R, s) ∂s α = 0 α = 1, . . . , M.
Additional conditions, such as the charge neutrality constraint for classical models of electrodes, may be associated with the auxiliary dofs. These conditions are typically expressed as where C is the number of the additional constraints that we shall assume to be, as it is often the case, functions of the s variables alone. Equation (2) implies that not all variations of the s are independent, and this must be accounted for in the search of the minimum of the potential. For example, in Born-Oppenheimer dynamics, the conjugate-gradient minimization is conducted via a constrained search, while in Car-Parrinello schemes, the additional constraints are added in the evolution equations derived from the Lagrangian. In the following, we indicate withs the values of the auxiliary variables satisfying both the minimum condition on the potential and the additional constraints, if they are present. Due to the dependence of the potential on R,s =s(R). In the adiabatic regime, then, the evolution of the physical variables is given by where m i is the mass of ion i. As discussed in the Introduction, current approaches for the numerical solution of the equation above have limitations that justify the development of alternative schemes. Starting from the next subsection, the MaZe approach is presented.

MaZe dynamical system
The key idea of the mass-zero constrained dynamics is to construct a Lagrangian that includes the s as dynamical variables and to interpret Eq. (1) as a set of holonomic constraints. In this section, MaZe is presented for the general case in which additional conditions must be satisfied by the fast dofs. In this case, the system is further extended to include also the Lagrange multipliers associated with these conditions as auxiliary variables.
To set the stage, let us introduce the auxiliary function where ν = {ν 1 , ..., ν C } are Lagrange multipliers associated to the additional constraints. The solution for s and ν satisfying Eqs. (1) and (2) is then given by the stationary point (ŝ,ν) of W (R, s, ν) [43,44]. This leads to the M + C conditions (5) Note that the last set of equations above represents in fact the additional conditions of Eq. (2), now obtained as the result of an optimization problem in the space that includes the Lagrange multipliers ν as additional variables. In the absence of additional conditions, on the other hand, W (R, s, ν) ≡ V (R, s) and the M surviving conditions above reduce to Eq. (1). At this stage, a finite fictitious "mass" is assigned to both sets of auxiliary variables. 1 Indicating the fictitious mass for the s and ν variables as μ s and μ ν , respectively, the Lagrangian for the extended system is defined as From this, the constrained evolution equations are obtained as The equations above can be simplified by observing that, in the first line, and that the forces acting on the auxiliary variables are null, because they coincide with the constraints. Thus, dividing both sides of the equations for the auxiliary variables by their masses Let us now consider the limit of zero mass for both sets of auxiliary variables. For this limit to be well defined, it needs to be taken in a homogeneous way. We then set μ ν ≡ μ s /κ with κ a nonzero constant of dimensions given by the ratio of the masses, so that the last equation in the system above becomes We can now take the limit μ s → 0 in Eq. (8). In order for the auxiliary variables to have finite acceleration, the ratio γ β = lim μs→0 λ β μs must remain finite, implying that the Lagrange multipliers λ β are proportional to μ s . In the limit of zero mass for the auxiliary variables, then Equation (10) defines the mass-zero constrained dynamics and it enables to recognize most of the interesting properties of the approach mentioned in the Introduction. First, since the Lagrange multipliers λ β go to zero with the auxiliary masses μ s , the evolution of the physical variables does not depend directly on the constraints. Second, the dynamics of the s and ν, controlled only by the constraint forces, satisfy by construction all the conditions imposed on the system. This implies that these conditions are automatically fulfilled also in the first of Eqs. (10) which is then equivalent to Eq. (3). Thus, by rigorously enforcing the mass-zero limit for the auxiliary variables, the system above provides a classical evolution for all degrees of freedom that leads to the exact adiabatic dynamics for the physical ones. Third, the numerical integration of the first equation can be performed with any standard MD algorithms (e.g., Verlet) with a timestep determined only by the force acting on the physical dofs. In addition, at each timestep, 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), ν(t + dt)) = 0, at the position predicted by the MD algorithm as described in Refs. [22,23]. This approach prevents propagation of the error between values of the variables at different timesteps. In the current implementations of the approach, the constraints are satisfied via the SHAKE iterative algorithm, which was proven to be symplectic and time-reversible [45,46]. Note that the homogeneous mass-zero limit has introduced an unknown scaling factor κ in the equations. The choice of the numerical value of this parameter depends on the specific system and is discussed more in detail in Refs. [3,4], where it is also shown that MaZe results are very stable with respect to this choice.
To conclude this section, note that when no additional constraints are present, i.e., in the absence of the ν variables, the MaZe system reduces to This form of the evolution equations is appropriate, for example, to simulate the shell model and, in view of the specific application considered in the results section and to simplify the notation, we shall adopt it in what follows. In particular, in the next subsection, we shall prove Eq. (11) samples exactly the Born-Oppenheimer probability density for the physical variables, i.e., the last MaZe property mentioned in the Introduction. This is an interesting result, because the use of constraints may induce a nontrivial metric in the phase space of the system [47,48] and require appropriate reweighting of statistical properties in the physical phase space. Furthermore, importantly, we shall show that approximate adiabatic separation, i.e., performing a dynamics with μ = 0, can induce a bias in the statistical properties of the system.

Statistical mechanics of the Mass-Zero constrained evolution
The discussion in this subsection summarizes the proof presented in Ref. [4] and is reported here for completeness and for the reader's convenience. Let us start by reconsidering the extended system before the mass-zero limit is taken. In the absence of additional conditions, the Lagrangian is given by The statistical mechanics of the system is described more naturally using (at first) a convenient set of generalized coordinates and in the Hamiltonian formalism. Proceeding in analogy with Ref. [48], we then start by performing the change of variables In the following, we shall use the notation υ = (ρ, σ) where ρ = {ρ 1 , ..., ρ N } and σ = {σ 1 , ..., σ M } and observables expressed in the new variables will be denoted in calligraphic font. The Hamiltonian of the system can be obtained via standard Legendre transform of the Lagrangian L(υ,υ) = L R(υ), s(υ),Ṙ(υ,υ), s(υ,υ) , and is given by where the momentum π υ is (i = 1, . . . , N and α = 1, . . . , M), and we have also introduced M −1 (υ), i.e., the inverse of the metric matrix associated with the new variables. For future convenience, the metric matrix and its inverse can also be expressed in block form as where A and´are 3N × 3N matrices, B and E are 3N × M matrices, and`and Z are M × M matrices, whose expressions are given in Appendix A.1.
The average of an observable O(ρ, σ, π ρ , π σ ) in the constrained microcanonical ensemble is given by where Z is the partition function. The delta functions in the equation above express the constant energy condition (second line) and the constraints (first line). Note that, in addition to the delta function associated with the holonomic constraints, δ M (σ), the integrand contains a delta function involving the momenta π σ . This delta originates from the fact that, in order for the constraints to be satisfied at all times, the additional conditionσ = 0 must hold. Using the relationυ = M −1 π υ , this implies (see also Appendix A.2) that, when the constraints are imposed, the momenta must satisfy where the tildes indicate that all matrices are evaluated at σ = 0. We now move to determine the expression for the average after integration over the variables associated with the constraints. To that end, the integral over π σ is evaluated first, followed by the change of variables σ α → s α , ρ i → R i , and a last integration over the s variables. The last two steps are discussed more in detail in Ref. [3] and in Appendix A.2. Here, we report the result of these operations, which is given by In the equation above,s =s(R) is such that σ(R,s) ≡ 0 (we assume, as commonly done in the Born-Oppenheimer framework that this expression has, for any R, a single root), andπ σ is defined in Eq. (19). The Hamiltonian is also evaluated on the hypersurface σ = 0, π σ =π σ where the constrained motion takes place. Its explicit form is derived in Appendix A.3 and is equal to Equation (20) defines the average of an observable with respect to a marginal probability where the constraints (or equivalently the auxiliary variables) have been integrated over. This marginal probability, however, still depends on the value of the mass, μ s , associated with the auxiliary dofs. This is apparent in the definition of the mass matrix A (and its inverse) and therefore of the generalized momenta π ρ , and from the dependence of the observable onπ σ . Indeed, from the definition in Appendix A.1, we have where D jj = m j δ jj and R jj = M α=1 ∂sα ∂ρ j · ∂sα ∂ρ j , so that Furthermore, the relation π υ = Mυ implies, π ρ = A · ρ + Bσ giving, on the constrained hypersurface where in the last equality we usedρ =Ṙ, as implied by the change of variable R → ρ. Finally (see Appendix A.1), on the constrained hypersurface, we also havẽ that also carries a dependence on the mass μ s due to the definition of the matrices `and B. Let us now consider the limit μ s → 0. `and B are proportional to μ s (see Appendix A) and therefore vanish in the limit, while lim μs→0 A = D. In the zero auxiliary mass limit then,π σ = 0 and the Hamiltonian of the system becomes with D −1 = 1 mj δ jj and where we have used the fact that, from Eq. (24), in the null auxiliary mass limit, Substituting in the expression for the average, we obtain The result above implicitly defines the microcanonical marginal probability in the physical phase space in the full adiabatic limit. This definition is in agreement with the form usually assumed for the Born-Oppenheimer probability. The discussion above also indicates that the dynamical system rigorously samples this density only in the full μ s → 0 limit and that, for finite auxiliary masses, corrections to the mass matrix associated with the momenta would be required, as indicated by Eqs. (23) and (24).
The shell model is one of the first attempts to represent polarization effects via empirical potentials, with specific focus in taking into account dipole polarization. This is described by assuming that the ion's total charge is divided between a core (representing the nucleus) and a massless shell (representing the electronic charge density). Each core-shell pair is bound by a harmonic potential and feels electrostatic interactions with the other particles. The ions evolve according to full adiabatic dynamics, i.e., subject to a force computed with the auxiliary variables at the minimum of the potential. Although several refinements have been proposed, the shell model still represents a valid benchmark for classical polarizable potentials and was chosen in this work, focused on exploratory calculations for a new development of MaZe, due to its simplicity. Polarization effects are important to capture accurately features of ionic systems ranging from phonon dispersion curves to structural and transport properties [49]. In the following, we shall consider how they influence charge diffusion in the presence of an external magnetic field. The model we consider does not take into account magnetization effects and therefore identifies the magnetic induction field B and the applied external magnetic field, usually noted as H. The generalization of the empirical Hamiltonian and of the integration algorithm that we describe in Sect. 3 to account for magnetization is nontrivial and beyond the developments presented in this work. To the best of our knowledge, however, this is the first time in which the effect of a magnetic field on the transport properties of a classical charged system is explored with any classical model of polarization and we consider it a first test on the way to more realistic simulations. Furthermore, in Sect. 4, we shall investigate a simple molecular liquid-molten sodium chloride-for which magnetization is likely to play a minor role.
Assigning a charge Q i to core i and q α to the shell α, the MaZe dynamics for the system in magnetic field is conveniently obtained, in analogy with the discussion in Sect. 2.1, by first considering the Lagrangian where a finite mass μ has been (temporarily) associated with the shell variables. In the equation above, V (R, S) is the total interaction potential, whose form is detailed in Appendix B, and we have introduced the at position r. We shall consider a system in a constant magnetic field parallel to the z axis: B = (0, 0, B z ). In the Coulomb gauge (∇ r · A(r) = 0), a valid choice for the vector potential is then A(r) = B z /2(−y, x, 0). The dynamics of the system is defined as fully adiabatic: the shells are assumed to adapt instantaneously to the positions of the cores, so that the force on each shell variable is null In the equations above, α = 1, ..., N and we have written explicitly the conditions of zero force for the components on the xy plane, which include the Lorentz force, and the component along the z axis, i.e., parallel to the magnetic field. Equation (29) can still be interpreted a set of 3N constraints, but, for this system, the components of the force on the plane orthogonal to the field depend on the velocity and the corresponding constraints are therefore no longer holonomic (the σ z α , on the other hand, are holonomic). The linear dependence on the velocity of these constraints, 2 however, still enables to write the Lagrangian equations of motion. In fact, for systems with mixed holonomic and semiholonomic constraints, these equations can be expressed as [31,32,[50][51][52][53][54] where we have indicated all the dynamical variables with the notation R 3N +3N ξ = (R, S). Following the same steps described in Sect. 2.1, the MaZe dynamical system is derived by first obtaining the evolution equations for the system on the basis of Eq. (30), then rearranging the evolution equations for the shell variables exploiting the condition of null force as in Eq. (8). In the μ → 0 limit, the resulting dynamical system then is where, as in the previous case, γ h β = lim μ→0 λ h β /μ.

The MaZe algorithm for the shell model
The numerical integration of Eqs. (31) must take into account two nontrivial features: the velocity dependence of the Lorentz force, which prevents direct use of standard integration algorithms (e.g., velocity Verlet), and the presence of nonholonomic constraints. These difficulties are solved combining a symplectic algorithm recently introduced to integrate the dynamics of ions in constant external magnetic field [55], with an adapted SHAKE algorithm to update the shell's positions and velocities. To set the stage, we introduce the auxiliary dynamical systeṁ where V ≡ V (R, S) and where we have introduced the notation (33) Taking the time derivative of all positions, it is immediate to show that the system above is equivalent to Eq. (31). A convenient integration algorithm can now be obtained by exploiting the Liouvillian formalism and writing the single timestep evolution operator associated with Eq. (32) as where the Liouvillian at the exponent is defined as iL X =Ẋ · ∇ X with i the imaginary unit and X = {R, P , S, p}, and where, for example, P = {P 1 , ..., P N }. To proceed, U(δt) is approximated via the following Trotter splitting In the equation above, we have separated, as commonly done [56], the differential operators acting on the coordinates R = {R 1 , ..., R N } and momenta of the ions and the shells. Note that, in the exploratory calculations presented here, we have employed a mixed Trotter break up, which is symmetric for the cores' momenta and simple for the shell variables and the cores' positions. The overall error in the approximation of the propagator is then of order δt 2 . One more observation is necessary to proceed. Focusing on the physical dynamical variables, we have In the absence of an external magnetic field, the exponentials of the Liouvillians in Eq. (35) correspond to simple translation operators for the components of the momenta (e i δt 2 L P ) and coordinates (e iδtL R ). Furthermore, always in the absence of an external magnetic field, the Liouvillians for the different Cartesian components of these variables commute among themselves so, for example, e i δt 2 L P = e i δt 2 L P x e i δt 2 L P y e i δt 2 L P z and the corresponding translations can be applied sequentially. Operating from left to right on the phase-space variables with the full set of translation operators, the velocity Verlet algorithm is recovered. When a (constant) magnetic field is present, however, the Liouvillians corresponding to the x and y components of the dynamical variables no longer commute. For example, by applying the commutator of the Liouvillians to a generic function of X, it can be seen that A nonzero result is obtained also for [iL P x , iL P y ], while the remaining commutators are zero. For this reason, two more Trotter break-ups become necessary to write the single-step propagator as a sequence of translations for the different variables. In particular, following Ref. [55], we choose the splitting: The actions of these operators can be directly translated into a set of instructions taking the system from time t to time t + δt and this is detailed in the insets in Fig. 1. In particular, the action of the momentum translations induced by the operators in the first line of the equation above corresponds to the "Update Core Momenta (1)" in the figure, while "Update Core Positions" shows the effect of the coordinates translations induced by the operators in the second line, and "Update Core Momenta (2)" derives from application of the operators in the last line. As indicated in the figure, the evaluation of the forces required to implement the second update of the core's momenta (these forces depend on the shell variables at time t + δt) must be performed under the condition that the updated shell velocities and positions satisfy the constraints Eq. (29).
These updated shell variables are determined via a straightforward generalization of the standard SHAKE algorithm (see also Refs. [31,32]) that preserves the key conceptual steps of the method. After the update of the positions of the cores, and based on the propagators associated with S and p (which is equal toṠ) in Eq. (38), the shell variables are advanced as The value of Γ , at this stage, is yet undetermined. Following the SHAKE strategy, these Lagrange multipliers are computed a posteriori imposing that the advanced shell positions and velocities satisfy the constraints at time t + δt: with, as usual, α = 1, ..., N and l = x, y. The expressions above are a system of nonlinear equations in the unknowns Γ that is conveniently solved using the SHAKE algorithm. This is an adapted Newton-Raphson method [44] in which the Lagrange multipliers are determined iteratively according to where the superscript (n) indicates the iteration step. In the equation above, the vector of parameters η = (η x , η y , η z ) was introduced. In standard SHAKE calculations, this vector is not present. However, following a common practice in minimization algorithms [57], it has been shown [58] that using a scaling factor to modulate the magnitude of the SHAKE update can improve convergence. As discussed more in detail in Sect. 4, for the particular problem considered in this paper, it proved useful to use a different scaling factor for the components of the constraints perpendicular and parallel to the field. This is due to the different nature (holonomic and semiholonomic) of the constraints shown in Eq. (29). In Eq. (42), we have also adopted the notation to define the vector of the constraints, and J d for the (i.e., the diagonal approximation of the Jacobian matrix of the standard Newton-Raphson method), with a, b = 1, ..., 3N . This matrix and the vector of the constraints are updated at each iteration step with the updated shell's positions and velocities that are computed according to in [23]. The iteration process is stopped when the modulus of the largest constraint, i.e., max a |Σ a |, becomes smaller than a predefined tolerance, typically chosen as close to the numerical precision achievable on the computer.

Simulation setup and results
The algorithm detailed in the previous section is applied to compute static and transport properties in a shell model simulation of molten NaCl in external magnetic field. The simulated system contains 108 Na + and 108 Cl − , placed in a cubic box of side L = 19.87Å, corresponding to a density ρ = 1.3113gcm −3 . In the presence of our vector potential, standard periodic boundary conditions can be preserved by imposing that images enter the box with their velocities [55]. This procedure ensures continuity of our Hamiltonian and of the relevant observables. The temperature of the system is set to T ≈ 1350K. The specific form of the interatomic potential, V (R, S), is given in Appendix B, and it is similar to the one adopted in Ref. [38]. We present results for the polarizable system in the presence and in the absence of a constant magnetic field directed along the z axis. The intensity of the field is chosen, so that the magnitude of the Lorentz force on each particle is comparable to that of the forces originating from the other interparticle interactions. As in our previous work on the shell model [2], we do not use the method of Ewald sums in the simulations, but rather truncate all interactions at a cut-off radius r c = 0.5L. The truncation of the long-range forces, sometimes adopted in simulations of large systems or when accuracy on the energy is not critical [59][60][61], was enforced for convenience. Our calculations are intended as a proof-ofprinciple validation of the MaZe dynamics in magnetic field and, while qualitative trends in the observables will be described, we are not focused on a realistic description of the system. Note that incorporating Ewald sums in the algorithm does not pose a conceptual problem nor has a significant effect on the numerical cost of the approach, as shown in Ref. [4] where a state-of-the-art classical model of polarization was considered.
The simulations are initialized as follows. The Cl − and Na + cores are placed on the sites of a simple cubic lattice and then displaced by a small uniform random amount in all directions. Initial shell positions are then found via a conjugate-gradient minimization of the potential energy with respect to these degrees of freedom. Velocities for all the degrees of freedom are set to zero. After the calculation of the interatomic forces, the first half of the evolution algorithm is applied to the cores to compute R(δt) and P (δt/2). New values of the shell variables are determined by applying SHAKE starting from the shell positions found by the conjugate-gradient minimization at step zero. Finally, interatomic forces are computed for this new configuration and the second half of the evolution algorithm is applied to the cores to obtain P (δt). Note that, because the shell velocities are set to zero at initialization, this procedure ensures that constraints are satisfied at t = 0 both in the absence and in the presence of the magnetic field. MaZe integration as described in Sect. 2.1 (no magnetic field) and Sect. 3 (magnetic field present) is then started. The timestep for the standard MaZe simulations is set to δt = 1fs, while for the runs with B z = 0, δt = 0.25fs (see below for a discussion of the reasons for the smaller timestep). Equilibration to the target temperature is achieved in all runs by simulating the system for 5ps. During this equilibration, the velocities are rescaled if the temperature differs from the target more than 10%. NVE runs of total length of 10ps are then performed to compute the properties reported in the following. A strict convergence criterion for the constraints is enforced in all runs by imposing that the maximum magnitude of the constraints is less than 10 −10 units of force. The relaxation parameter η, see Eq. (42) and discussion in the previous section, is set to η = (1, 1, 1) for the calculations with B z = 0 and to η = (0.33, 0.33, 1) when the magnetic field is present. While the value of η can be set via an automatic search and adapted during the run [58], here it was chosen via manual search by optimizing the number of iterations necessary to converge SHAKE for a typical configuration of the cores. Previous experience, confirmed by the simulations reported here, has shown that this is sufficient to provide a stable number of iterations along the whole trajectory. In Fig. 2, we show the convergence paths of the constraints in the presence (left panel) and absence (right panel) of the magnetic field for a few randomly chosen configurations along the trajectory of the cores. The figures show the magnitude of the largest constraint as a function of the number of SHAKE iterations. In agreement with previous calculations, the number of iterations needed to converge in MaZe calculations for fully holonomic constraints (no magnetic field) is very small and similar for different configurations. The path shows monotonic convergence, with a single slope on the semilogarithmic scale employed in the figure. This fast convergence is facilitated by the already small value of the largest constraint at the start of the iterative process, indicating that the provisional values for the shell positions are quite close to the minimizers of the potential. The convergence of MaZe for the mixed set of holonomic and semiholonomic constraints, on the other hand, is about four times slower. This may be related to the fact that the magnetic force results in larger nondiagonal terms in the Jacobian matrix of the Newton-Raphson procedure, implying that Eq. (43) provides a less effective approximant of its inverse. Furthermore, the value of the maximum constraint at the beginning of the minimization is now larger than in the holonomic case, indicating that our provisional shell positions and velocities at the zeroth iteration are farther from the final solution. This behavior is sensitive to the choice of the timestep δt, with worse performance (and even-tually lack of convergence) with larger timesteps. This suggests that the basin of convergence of SHAKE in the presence of semiholonomic constraints may be smaller than the one for standard applications, an issue that will be further investigated in future studies. Finally, the paths to convergence now present a double slope pattern: a fast initial decay is followed by a slower decrease. Closer inspection of the decrease of individual constraints suggests that this is due to the different speed of convergence of the holonomic and semiholonomic constraints, with the latter evolving faster toward the threshold. This is most likely also related to the different values for the components in the scaling vector η. In spite of the differences in the convergence pattern, and of the need to further investigate the behavior of the new algorithm, this first implementation of SHAKE for mixed constraints performs well for the nontrivial interactions of the model.
Inspection of typical dynamical indicators and structural properties confirms the reliability of the MaZe approach for classical polarizable models in external magnetic field. In particular, the fluctuations of the total energy relative to the fluctuations of the potential energy along the trajectory of the cores are ΔE/ΔV ≈ 2 · 10 −4 , where, for example, ΔE = E 2 − E 2 . The stability of the new algorithm is visible also in the calculation of the instantaneous temperature of the cores. Figure 3 shows the fluctuations of this quantity in a 10ps simulation in the presence of the magnetic field, which are-again-perfectly compatible with typical results for classical simulations. As a test of the reliability of MaZe semiholonomic dynamics, we consider the radial distribution function of the ionic species as obtained in the simulation with and without external magnetic field. Results are shown in Fig. 4, where we report as symbols the output of the runs in the absence of the magnetic field and as solid lines that of the calculations with the magnetic field. The position and shape of the peaks for all g(R) are in good agreement with experimental results [64] and with previous calculations [2,55] in spite of the somewhat crude treatment of the electrostatic interactions, of a different temperature (T ≈ 1350 K in this work and T ≈ 1550 K in Refs. [2,55]) and of some differences in the parameters in the shell model detailed in Appendix B. Perhaps more importantly for our purposes, the curves and the symbols are superimposed. This provides strong validation for the MaZe algorithm presented in Sect. 3: it is in fact known (see, for example, Ref. [55]) that timeindependent averages are not affected by the presence of the magnetic field. As shown by the results on a nontrivial observable, this property is respected by the MaZe algorithm.
We now move to the calculation of time-dependent statistical properties of the system. In this case, the presence of a magnetic field is expected to affect the results in nontrivial ways providing further and more interesting testing ground for our approach. We consider, in particular, the velocity correlation functions of the ionic species. In Figs. 5 and 6, we show results In the insets, we also show the corresponding elements of the diffusion tensor. As expected, when B z = 0, the three diagonal components of the velocity correlation function are equal for each species and show the characteristic initial decay followed by one or more minima before going to zero at longer times. For this system, the integral of the autocorrelation function of the velocity yields diffusion coefficients equal to D Na = (1.78 ± 0.08) · 10 −4 cm 2 s −1 and D Cl = (1.51 ± 0.06) · 10 −4 cm 2 s −1 , both obtained averaging the three components on the diagonal of the diffusion tensor. The error on the values of the diffusion is estimated from the off-diagonal components of the tensor for B z = 0. These off-diagonal components must be zero based on time-reversal symmetry arguments [65,66], and the results of the simulation can be used to estimate the statistical noise in computing the integrals. All the results discussed above are compatible with previous studies [64] performed on the same system at zero magnetic field. The presence of the magnetic field breaks the isotropy of space and this implies that the components of the correlation on the plane orthogonal to the field are now different from that in the direction parallel to it, and show an oscillatory behavior that reflects the rotatory motion induced by the Lorentz force. Consistently, the diffusion coefficients are also affected by the presence of the magnetic field. The observed reduction of their values is in fact a known phenomenon, the so-called magnetoresistance, which is commonly observed for the electrons in semiconductors in the presence of magnetic field and was also reported for ions in previous simulations [55,62]. In particular, the diffusion coefficients are reduced to D ⊥ Na = (0.96 ± 0.10) · 10 −4 cm 2 s −1 , D zz Na = (1.36 ± 0.14) · 10 −4 cm 2 s −1 for the Sodium ions and to D ⊥ Cl = (0.91 ± 0.07) · 10 −4 cm 2 s −1 , D zz Cl = (1.24 ± 0.10) · 10 −4 cm 2 s −1 for the Chlorine ions, where D ⊥ = (D xx + D yy )/2. The effect of the magnetic field is even more striking when considering the offdiagonal components of the velocity correlation functions. In Fig. 6, we present results for the xy and yx cross-correlations. In the absence of the field, timereversal invariance leads to null values of these quantities. On the other hand, when B z = 0 a characteristic oscillatory pattern is observed. As detailed in Ref. [66], the antisymmetry of these two observables is dictated by their properties under generalized timereversal symmetries and well reproduced by our simulations (all other off-diagonal components remain zero, for symmetry reasons). The behavior of the correlation function is reflected in the values obtained for the xy and yx components of the diffusion tensor, which are now equal to D xy Na = (0.46 ± 0.14) · 10 −4 cm 2 s −1 and D yx Na = (−0.43 ± 0.14) · 10 −4 cm 2 s −1 for Sodium and to D xy Cl = (−0.32 ± 0.10) · 10 −4 cm 2 s −1 and D yx Cl = (0.32 ± 0.10) · 10 −4 cm 2 s −1 for Chloride.
Finally, it is interesting to explore the effects of polarization on the transport properties of this model of molten NaCl. To assess the relevance of these effects, we compare the elements of the diffusion tensor discussed above with those from a simulation of an unpolarized model of the system. The unpolarized (or rigid ion) model is defined by removing the shell variables from The presence of polarization globally enhances transport in the system, but-due to the different polarizability of the two ions-leads to different values, in particular, of the component of the diffusion tensor parallel to the magnetic field. Similarly, the cross components of the diffusion for Na + and Cl − in the plane orthogonal to the magnetic field (bottom panel of the figure) are essentially identical for rigid ions, but polarization separates them. In particular, for the rigid ion case, we obtain values for the xy and yx components of the diffusion coefficients given byD xy Na = (0.25 ± 0.14) · 10 −4 cm 2 s −1 andD yx Na = (−0.24 ± 0.14) · 10 −4 cm 2 s −1 for Sodium andD xy Cl = (−0.26 ± 0.10) · 10 −4 cm 2 s −1 andD yx Cl = (0.25 ± 0.10) · 10 −4 cm 2 s −1 for Chloride, where the sym-bolD is used to indicate that the diffusion coefficient is computed for the rigid ion model, at difference with the notation D which indicates diffusion coefficient computed for the shell model. This has an interesting implication for the detection of the ionic Hall effect in molted NaCl. In fact, the key indicator of this phenomenon, the Hall mobility, in the Nerst-Einstein approximation, is given by [62] When the off-diagonal components of the diffusion tensor of the two species are equal and opposite, as in previous more refined calculations on a rigid ion model for the system [55,62] and, within errors, in the results shown in Fig. 7, the mobility is obviously null. In particular, the values of the diffusion coefficients obtained from the rigid ion simulations performed in this work yield a value for the Hall mobility given byμ H ≈ −0.9 · 10 −5 cm 2 V −1 s −1 . On the other hand, the diffusions obtained with the shell model result in μ H ≈ 8.5 · 10 −5 cm 2 V −1 s −1 . While error bars (not reported) are still quite large with our level of statistics, the noticeably different values for these mobilities suggest that Hall effect is absent or hardly measurable for the rigid ion model but quite appreciable when polarization is accounted for. This observation needs to be confirmed via more accurate calculations, but it clearly underlines the relevance of polarization on observables affected by relatively subtle effects in transport processes.

Conclusions
In this paper, we described the MaZe dynamics for the simulation of systems where the evolution of a set of physical dofs depends on parameters subject to assigned conditions. Fully adiabatic dynamics, in first principle or classical polarizable models, is perhaps the most relevant example of such systems. The derivation of MaZe and of its key properties was presented, using the classical shell model for polarization as a reference case. MaZe exploits the Lagrangian formulation of classical mechanics to define an extended system in which the external parameters evolve as auxiliary dynamical variables of zero mass. These variables are subject to constraints that strictly enforce the conditions on the parameters for each configuration of the physical dofs in the evolution. The mass-zero value for the auxiliary variables results in rigorous fully adiabatic evolution for the physical dofs and, consequently, on exact statistical sampling of the associated probability density. From a numerical point of view, the integration of the constrained dynamics is efficiently performed using the SHAKE algorithm in its standard form for holonomic systems.
A new development extending this approach to the physically interesting case of the shell model in external magnetic field was also presented. Although the model we have chosen neglects magnetization effects, it is still nontrivial and represents an interesting extension of the MaZe framework. In this case, in fact, the presence of the Lorentz force requires to generalize the formalism and the associated algorithm to systems with constraints that depend linearly on the velocities. This generalization was described in the second part of the paper and used in illustrative calculations on a shell model of molten NaCl. These calculations demonstrate the effectiveness of the new algorithm and provide interesting qualitative information on the effect of polarization on ionic transport in magnetic field. In particular, we showed indications that-within the model adopted -polarization is critical to obtain a nonnull value of the Hall mobility for the system. Future work will consider extending MaZe to more general models for electric and magnetic interactions.
Acknowledgements The authors are grateful to Rodolphe Vuilleumier, Benjamin Rotenberg e Mathieu Salanne for enlightening discussions. Jean-Paul Ryckaert also deserves special thanks and credit for his role in the birth of the zero mass constrained scheme.
Funding Open access funding provided by Politecnico di Torino within the CRUI-CARE Agreement.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.   a more sophisticated polarizable potential proposed in [67] (in particular, we used that work to define the elastic constant k) and the Sodium core and shell charges were assigned by mirroring (with opposite signs) those used for Chlorine where the parameters Aij, Cij, Dij, and λ are the same as the shell model and the ionic charges Zi = Qi + qi are used in place of the ones summarized in Table 2.