Propagators for molecular dynamics in a magnetic field

Ab initio molecular dynamics in a magnetic field requires solving equations of motion with velocity-dependent forces -- namely, the Lorentz force arising from the nuclear charges moving in a magnetic field and the Berry force arising from the shielding of these charges from the magnetic field by the surrounding electrons. In this work, we revisit two existing propagators for these equations of motion, the auxiliary-coordinates-and-momenta (ACM) propagator and the Tajima propagator (TAJ), and compare them with a new exponential (EXP) propagator based on the Magnus expansion. Additionally, we explore limits (for example, the zero-shielding limit), the implementation of higher-order integration schemes, and series truncation to reduce computational cost by carrying out simulations of a HeH$^+$ model system for a wide range of field strengths. While being as efficient as the TAJ propagator, the EXP propagator is the only propagator that converges to both the schemes of Spreiter and Walter (derived for systems without shielding of the Lorentz force) and to the exact cyclotronic motion of a charged particle. Since it also performs best in our model simulations, we conclude that the EXP propagator is the recommended propagator for molecules in magnetic fields.


I. INTRODUCTION
Since the first simulations more than 50 years ago, 1,2 molecular dynamics has become a ubiquitous tool in computational chemistry, allowing for the calculation of reaction rates and energies, 3 the exploration of reaction networks, 4,5 and the prediction of vibrational spectra, [6][7][8][9] including infrared, Raman, and circular dichroism spectroscopies.The principal idea is to solve the equations of motion to determine the motion of a set of N nuclei with masses m 1 , . . .m N and positions x 1 (t), . . .x N (t).The potential a) Electronic mail: laurens.peters@kjemi.uio.noE(x) and the corresponding gradient are usually determined using force-field or (for smaller systems) ab initio methods, while the equations of motion are solved using propagators such as the velocity Verlet algorithm 10,11 or higher-order schemes. 12,13n a magnetic field B, a classical particle of charge q I experiences an additional Lorentz force: More than 20 years ago, Spreiter and Walter recognized that the appearance of this velocity-dependent force requires a modification of the velocity Verlet scheme. 14][19] Taking into account the electrons as quantum particles within the Born-Oppenheimer approximation, the energy becomes dependent on the magnetic field, [20][21][22][23][24][25][26] while the velocity-dependent force assumes a more gen-eral form: [27][28][29] m I ẍI = − dE(x, B) dx I + N J=1 A IJ (x, B) ẋJ .
Here, the last term contains not only the bare Lorentz force acting on the nuclei [see eq. ( 2)], but also a contribution from the Berry curvature, 30,31 reflecting the shielding of the nuclei by the electrons 32,33 and introducing a coupling between the motion of different nuclei. 34,35ile the form and the implications of eq. ( 3) were discussed more than three decades ago, 27,36,37 there are only a few examples where it has been actually solved for a molecular system in a magnetic field.Ceresoli and coworkers simulated an H 2 model system in 2007, integrating the equations of motion with a Runge-Kutta scheme. 28In 2021, Peters and coworkers were the first to conduct accurate dynamics of H 2 using the auxiliarycoordinates-and-momenta (ACM) propagator, 29 which is derived from a general scheme for propagating nonseparable Hamiltonians by Tao. 38A more extensive study was conducted by Monzel and coworkers 39 in 2022, studying H 2 and LiH with the Tajima (TAJ) propagator, originating from particle physics. 40 this work, we introduce a new propagator that we refer to as the exponential (EXP) propagator.It is inspired by the fact that equations of type appearing, for example, in time-dependent Kohn-Sham theory 41 and in surface-hopping algorithms 42 , can be solved exactly using the Magnus expansion 43 and approximately using matrix exponential(s) of A II . 44We compare the EXP propagator to the previously published ACM and TAJ propagators, regarding their theoretical foundation, their time-step requirements, their performance during actual simulations, and properties, such as their behaviour in the zero-field case [see eq. ( 1)], in the zero-shielding case [see eq. ( 2)], and in the cyclotron limit [only the Lorentz force in eq. ( 2)].We use simulations of a HeH + model system to corroborate our theoretical findings and test schemes to further reduce the computational cost.
Having established the equations of motion, their propagation, and the weak-field limit for the step size in Sections II.A-C, respectively, we derive the working equations of propagators in the absence of a field as well as the ACM, TAJ, and EXP propagators in Sections II.D-G and compare them in Section II.H. Computational details for the HeH + simulations are given in Section III, while the results of these simulations are presented and discussed in Section IV. Conclusions and an outlook are given in Section V.

A. Equations of motion
Within the Born-Oppenheimer approximation, the equations of motion of a molecule in a magnetic field are given by [27][28][29] Here, we use indices I, J, ... for the N nuc nuclei and R I , Z I , M I , and ṘI for the position, charge, mass, and velocity of nucleus I.The collective nuclear coordinates are denoted by the 3N nuc -dimensional column vector R, while we represent the uniform magnetic field B of strength B by the 3 × 3 antisymmetric matrix B in the manner so that The first term in eq. ( 5) is the gradient of the Born-Oppenheimer electronic energy, obtained at some ab initio level of theory, where H el (r, R, B) and φ(r; R, B) are the electronic Hamiltonian and wave function, respectively, both of which depend on the electronic coordinates r [over which the integration is performed in eq. ( 9)].The second term in eq. ( 5) is the Lorentz force arising from the nuclear charge screened by the surrounding electrons.This shielding arises from the Berry curvature calculated from derivatives of the electronic wave function with respect to the nuclear coordinates.For more details on the calculation (at the Hartree-Fock level of theory) and interpretation of the Berry curvature, see refs.30-33.With M being the 3N nuc × 3N nuc -dimensional matrix with the nuclear masses M I on the diagonal, we can define the kinetic momenta of the nuclei as We may now rewrite the nuclear equations of motion in eq. ( 5) in the more convenient form where F(R, B) is the 3N nuc -dimensional vector of the Born-Oppenheimer gradient forces and W(R, B) is a 3N nuc ×3N nuc matrix, whose IJ blocks of dimension 3×3 contain the Berry curvature as well as the contribution from the bare (unscreened) Lorentz force: The form of the equations of motion given in eq. ( 12) highlights that the effect of the magnetic field is twofold: It changes the potential energy surface leading to different Born-Oppenheimer gradient forces and it introduces an additional, velocity-dependent term.

B. Propagators: general considerations
The main objective of this theory section is to discuss how these modified equations of motion can be integrated efficiently using different propagators.We denote by the position and momentum, respectively, of the system at time a relative to time t in units of the time step ∆t.
Introducing the notation we may now express the equations of motion in terms of differentials with respect to the factor a The choice of the step length ∆t, here appearing as a prefactor, will be discussed in the next subsection.
In this notation, a propagator is defined as a scheme that updates the coordinates (x 0 → x 1 ) and momenta (π 0 → π 1 ) for a given ∆t.Since eqs.(18) and (19)  depend on each other, they are usually solved alternately for a series of substeps.To ease their reading, we will only derive working equations for a single set of those substeps, where x 0 and π 0 are initial values and a and b are arbitrary step lengths in units of ∆t.Equations for all other substeps can then be derived by adjusting the initial values and step lengths accordingly.
To evaluate the integrals in eqs.( 20) and ( 21), we will draw inspiration from the mean-value theorems.They state that, for real-valued functions ϕ(α) and ψ(α) that are continuous in [0, a], there exists a point 0 and The corresponding statements are not guaranteed to hold for vector-valued functions; however, we will at one point rely on them as approximations.Specifically, we use the forms and to approximate integrals over two vectors of continuous functions ϕ(α) and ψ(α).Note that both the midpoint rule and the trapezoidal rule can be derived from these mean-value approximations by choosing b = a/2 and, in case of latter, ψ(α) = 1: C. Choice of step size The error introduced by a given propagator depends on the applied step size ∆t.Before considering the different propagators, we consider briefly in this subsection the restriction on the time step imposed by an external magnetic field.As in field-free simulations, we demand that the time step ∆t is significantly smaller than the fastest molecular vibration.In addition, we demand that the matrix ∆t w a is convergent for all a, meaning that Introducing the spectral radius ρ which we loosely interpret as a cyclotron frequency, we have in some matrix norm • .Consequently, eq. ( 28) holds when This weak-field limit, as defined by Spreiter and Walter, 14 means that the time step is small enough to resolve the cyclotronic motion.For a neutral molecule, the chargemass ratio is roughly on the order of 10 −4 a.u., so that Time steps of up to 100 a.u.(2.4 fs) therefore allow for simulations in field strengths up to about 10 B 0 , covering the entire range of field strengths from the (particularly interesting) intermediate regime (< 1 B 0 ) up to the Landau regime.

D. Propagators in the absence of a field
For the field-free case, w is zero so that we can directly apply the mean-value approximation [eq.( 24)] to eqs. ( 20) and ( 21): As shown in Algorithm 1, both equations are solved alternately using the current forces and momenta as mean Algorithm 1: General algorithm for a propagator of order K, in the absence of a field.See eqs.( 14)-( 17) for definitions of the variables.a and b are predefined vectors of length K + 1 and K, respectively.In case of velocity Verlet, a = [0.5, 0.5] and b = [1.0].
values to propagate the momenta and positions, respectively.The remaining task is now to come up with a series of a's and b's (denoted by the vectors a and b) to approximate the mean values.The lengths of a and b (K + 1 and K) reflect the order K of the scheme.
In the standard velocity Verlet scheme (K = 1), 10,11 we set a = [0.5, 0.5] and b = [1.0],so that we obtain: x VV 1 = ∆t M −1 π VV 0.5 + x 0 (36) Note that this quadrature corresponds to solving the full integral (a = 1) in eq. ( 33) with the trapezoidal [see eq. ( 27)] and the full integral in eq. ( 34) (b = 1) with the midpoint rule [see eq. ( 26)].As a consequence of this the velocity Verlet algorithm is symplectic, timereversible, and of second-order accuracy.However, it has been shown that higher-order schemes, 12,13 can significantly increase the stability of the dynamics.
Unfortunately, there is no straightforward expansion of the upper scheme towards systems with velocitydependent forces since, taking the velocity Verlet algorithm as an example, the step of eq. ( 37) becomes The mismatch between the required (π 1 ) and the available (π 0.5 ) momenta for the propagation leads to a systematic error in the integration of the equations of motion. 29Clearly, there is a need to develop alternative propagators, which will be done in the following subsections.

E. Auxiliary-coordinates-and-momenta propagator
The auxiliary-coordinates-and-momenta (ACM) propagator method was proposed by Tao for general nonseparable Hamiltonians 38 and adapted by Peters and coworkers for molecular simulations in a magnetic field. 29The idea is to circumvent the mismatch in eq. ( 38), by introducing an additional pair of coordinates and momenta (x ′ , π ′ ) that are kept close to the original pair (x, π) during the dynamics (x ≈ x ′ and π ≈ π ′ ).If this coupling is sufficiently strong, we can write the propagation of π in terms of π ′ and x and the propagation of π ′ in terms of π and The only remaining step is the coupling of coordinates and momenta, which is done when all quantities x, π, x ′ , and π ′ are at the same time step by carrying out, for a given coupling constant s, the transformation 29 with the coupling matrix where Setting s = 0 reduces the coupling matrix to the unity matrix, while setting it to s = π/(c∆t) leads to a complete exchange of x and π with x ′ and π ′ and vice versa.
The pseudo code for a general ACM propagator of order K is given in Algorithm 2.

F. Tajima propagator
An alternative to the ACM propagator is the Tajima (TAJ) propagator.Initially used in particle physics, 40 it was rewritten for molecular applications by Monzel et al. 39 Here, we present a slightly different derivation of the working equations that starts with applying the meanvalue approximations [eqs.(24) and (25)] to eq. ( 20): w α π a dα + π 0 (46)   Algorithm 2: General algorithm for the auxiliary coordinates and momenta (ACM) method of order K. See eqs.( 14)- (17) for definitions of the variables.a and b are predefined vectors of length K + 1 and K, respectively.In the velocity Verlet variant, a = [0.5, 0.5] and b = [1.0].
Assuming that c = a/2 and that both integrals have the same mean value w b , we obtain: Introducing the inverted matrix we arrive at the working equation for the TAJ propagator, yielding the algorithm in Algorithm 3. In the weak-field limit [see eq. ( 28)], we can expand the matrix inversion in a Neumann series: In standard applications, we expect this series to converge fast.Using the estimate in eq. ( 32) with B = 0.1 and ∆t 10 in atomic units, the error is on the order of 10 −8 when using N = 2. Consequently, we can avoid the (comparably) high computational cost of matrix inversions during the simulations.
Algorithm 3: General algorithm for the Tajima (TAJ) method of order K. See eqs.( 14)-( 17) for definitions of the variables.a and b are predefined vectors of length K + 1 and K, respectively.In the velocity Verlet variant, a = [0.5, 0.5] and b = [1.0].

G. Exponential propagator
A general approach to solving first-order linear differential equations is the Magnus integrator method.Specifically, we may write where y α,a is the Magnus series, 43 y γ,a = ∆t From y a,a = 0 and it is straightforward to verify that eq. ( 51) solves eq. ( 18) and is thus an alternative to eq. (20).Since this holds for any γ, we will set it to zero from now on.Truncating the Magnus series after the first term and applying the meanvalue approximation [eq.( 24)], the matrix exponential becomes where we have introduced the matrix exponential Algorithm 4: General algorithm for the Exponential (EXP) method of order K. See eqs.( 14)-( 17) for definitions of the variables.a and b are predefined vectors of length K + 1 and K, respectively.In the velocity Verlet variant, a = [0.5, 0.5] and b = [1.0].
Equation (51) now becomes: Approximating the remaining integral via the midpoint rule [eq.( 26)], we obtain as the exponential propagator (EXP) for the nuclear momenta in magnetic fields; see Algorithm 4. The matrix exponential can be expanded as In contrast to the Neumann series, this series converges unconditionally, for all step sizes.

H. Comparison of the propagators
We close this section by discussing a few properties of the introduced propagators.The results are summarized in Table I.From Algorithms 2-4, we see that the ACM propagator is about three times more expensive than the other schemes, requiring 3K instead of K forces and Berry curvature calculations per step.In addition, it depends on the definition of a parameter (s), which has an influence on the stability of the dynamics.Moreover, unlike the ACM propagator, TAJ and EXP propagators converge to the correct zero-field solution when setting B and therefore w b to zero.
As the correct zero-shielding limit, we define the propagators derived by Spreiter and Walter, which were constructed for the special case where the Berry curvature is zero.While their working equations clearly differ from the ACM propagator, we can compare them to the TAJ and EXP propagators by assuming a single particle with charge Z 1 and mass M 1 , experiencing a time-dependent external force f and a magnetic field of strength B z in the z-direction.In this particular case, w and thus u a and v a are constants and the velocity Verlet variants of the EXP and TAJ propagator reduce to respectively.Expanding u a and v a up to second order and using the Taylor expansion of the forces in ydirection we obtain, for both EXP and TAJ, an expression that is identical to those obtained by Spreiter and Walter via inversion [see eq. ( 16) in ref. 14] and via Taylor expansion [see eq. ( 39) in ref. 14]: Their difference, however, becomes obvious when additionally setting the forces in eqs.( 60) and (61) to zero: In this cyclotron limit, the EXP propagator collapses to the exact result for the cyclotron motion of a single, charged particle No Yes Yes Correct zero-field limit?
No Yes Yes Correct zero-shielding, No Yes Yes velocity Verlet limit?Exact cyclotron limit?
No No Yes while TAJ introduces an error at the order of O([∆tω] 3 ).For the previously mentioned estimate [see eq. ( 32)] with B = 0.1 and ∆t = 10, this results in an error at the order of 10 −12 (all given in atomic units).Thus, in practice EXP and TAJ trajectories and their computational cost will be very similar.

III. COMPUTATIONAL DETAILS A. Diatomic model system
To conduct a study of the performance of the previously discussed propagators, we carry out simulations of a diatomic model system perpendicular to the magnetic field.It has been shown that, in this particular case, the Berry curvature depends solely on the Berry charges Q IJ (R, B) and Berry charge fluctuations P IJ (R, B) as well as the orientation of the molecule relative to the field vector: 30,33 where RIJ is the normalized interatomic distance vector: The resulting bond-length dependence of the energy and charges is shown in Fig. 1.We note that HeH + dissociates into He with effective charge Z He + Q HeHe = 0 and a proton H + with effective charge Z H + Q HH = 1.In the bonding region, both nuclear charges are partially screened by the electrons.Note that, in our calculations, we neglect the magneticfield dependence of the electronic energy itself E(d HeH ),  69).This approach allows for a more systematic study of the propagators, since only the velocity-dependent forces increase linearly with B, while the potential-energy surface and the electronic structure remain unaffected.Because of its overall charge, small mass, and significant geometry dependence of the screening [see Fig. 1(b)], the HeH + model system can be regarded as an extreme case, featuring (comparably) large cyclotron frequencies (ω).

B. Molecular simulations
All NVE simulations (constant number of particles, volume, and energy) were conducted in the xy-plane with the magnetic field aligned along the z-axis.Each simulation began at the equilibrium geometry, with random initial momenta.We investigated the ACM propagator (Algorithm 2), the TAJ propagator (Algorithm 3), and the exponential operator (Algorithm 4) at three different orders K: the velocity-Verlet (VV) scheme 10,11 with K = 1, the OM scheme of Omelyan et al. 12 with K = 4, and the RK4 scheme (S 6 /O4 in ref. 13) with K = 6.The corresponding factors a and b are given in Tab.II.
Each trajectory was simulated for 24 ps using an effective time step (∆t/K) of 0.1 fs, storing energies, geometries, and momenta every 2.4 fs.For the ACM propagator, we used an optimized coupling constant s = 0.013.Where N is not given, we used the standard matrixinversion and exponential algorithms of python3.6 to calculate v and u in the TAJ and EXP propagators, re-spectively, and the truncated series of eqs.( 50) and (58) otherwise.Estimating that ∆t ω = 4.5 × 10 −4 B, we apply a magnetic-field range of 10 −2 -10 3 in units of B 0 to ensure that the weak-field condition in eq. ( 31) holds for all simulations.
We use the standard deviation of the total energy (ε tot ) as a criterion for the stability of the dynamics.In addition, vibrational spectra 8,29 and the center-of-mass motion are used to compare trajectories of the different propagators.

IV. RESULTS AND DISCUSSION
In Fig. 2, we have plotted ε tot for HeH + as a function of field strength B for the ACM, TAJ, and EXP propagators with K = 1.The results for K = 4 and K = 6 are very similar and therefore not included in the figure .In general, all three propagators become less stable with increasing field strength-in particular the ACM propagator, for which the simulation at 10 3 B 0 becomes unstable.At this high field strength (characteristic of neutron stars rather than white dwarfs), a reoptimization of the coupling-strength parameter s may be required.In agreement with a previous comparison of the ACM and TAJ propagators, 39 we note that the ACM propagator performs less well than the TAJ and EXP propagators.Since, in addition, the TAJ and EXP propagators are parameter-free and since they require only one rather than three force and Berry-curvature calculations per step, we conclude that the the TAJ and EXP propagators are to be preferred over the ACM propagator.
While the EXP and TAJ propagators are indistinguishable at low field strengths, the EXP propagator becomes marginally more stable for B > 10 B 0 , as the cyclotronic center-of-mass motion of the charged HeH + system begins to dominate the trajectories.However, even at 10 B 0 , the center-of-mass motions and the vibrational spectra of the two propagators are still indistinguishable and in fact identical to those obtained with the ACM propagator; see Fig. 3.In line with previous studies on diatomics, 29,39 the peaks in the vibrational spectrum of HeH + (Fig. 3b) correspond to the cyclotronic center-ofmass motion (∼ 30 cm −1 , also visible in Fig. 3a), rotation (∼ 370 cm −1 ), and vibration (∼ 3350 cm −1 ), and feature splitting patterns that arise from the couplings between the different modes.
Since the EXP propagator is the most stable propagator, we restrict our attention to this propagator when investigating dependence on the propagator order K and on the truncation level N in the exponential series; see Fig. 4(a) and (b), respectively.We expect the TAJ propagator to show a similar behaviour, while an investigation of the influence of K on the ACM propagator can be found in ref. 29.
At low field strengths, the higher-order schemes (i.e., OM with K = 4 and RK4 with K = 6) improve the stability of the EXP propagator by up to two orders of magnitude; see Fig. 4.These higher-order schemes therefore allow for a significantly larger effective time step ∆t/K (constant in Fig. 4), reducing the computational cost of the dynamics.Perhaps surprisingly, the OM and RK4 schemes become less stable than the lower-order VV scheme as the velocity-dependent forces become dominant in the Landau regime (beyond 10 B 0 ).This might be due to the fact that the coefficients a and b (see Algorithm 4) were not optimized for such propagations.
Simulations with a truncated exponential series at N = 2 in eq. ( 58) reproduce the results obtained with the exact matrix exponential for field strengths up 10 B 0 ; see Fig. 4. Since calculating a matrix exponential is (comparably) expensive, replacing it with a series will reduce the computation time.For N = 2 at higher field strengths, the truncation error exceeds O [ω∆t] 2 = 10 −5 , giving unstable dynamics.This indicates that N must be chosen with care.

V. CONCLUSION AND OUTLOOK
In this work, we have studied three propagators for molecular dynamics in a magnetic field-namely, the auxiliary-coordinates-and-momenta (ACM) propagator, the Tajima (TAJ) propagator, and the new exponential (EXP) propagator, testing their performance using simulations of a HeH + model system at different field strengths (10 −2 -10 3 B 0 ).
The EXP propagator, which was derived from a truncated Magnus expansion, correctly reduces to the standard velocity Verlet propagator in the zero-field limit, the propagators of Spreiter and Walter 14 , neglecting electron shielding, and the cyclotronic motion of a charged particle.Since it also performed best in our model simulations and showed the best stability, especially at higher field strengths, we recommend the EXP propagator for simulations of molecules in a magnetic field.However, the TAJ propagator delivers identical results at field strengths ≤ 1 B 0 and the more expensive ACM propagator might be useful for cases where the energy and/or the electron shielding depend on the nuclear momenta.
In the study of molecules in strong magnetic fields, we are usually interested in field strengths below 1 B 0 and  apply time steps of ∆t = 0.1 fs to resolve the molecular vibrations.In practice, these time steps are small enough to resolve the cyclotronic motion (weak field limit ), while, additionally, allowing for the use of higher-order schemes and truncation of the exponential series to reduce the computational cost of the EXP propagator.
Our model system has been designed to reproduce the properties of HeH + at field strength 0.1 B 0 , with energies E(d HeH ), Berry charges Q IJ (d HeH ), and Berry charge fluctuations P IJ (d HeH ) obtained by spline fitting of Hartree-Fock (HF)/lu-aug-cc-pVTZ[45][46][47] results for different bond lengths (d HeH ) at this field strength.The prefix "lu-" indicates the use of an uncontracted London orbital basis set, as suggested in ref.48.All calculations were performed with the London 49 program package.

FIG. 1 :
FIG. 1: Bond-length (d HeH ) dependent energies (a) and charges (b) derived from the Berry curvature of HeH + at the HF/lu-aug-cc-pVTZ level of theory perpendicular to a field of B = 0.1 B 0 .We consistently use atomic units.

FIG. 3 :
FIG. 3: Center-of-mass motion (a) and vibrational spectra of He (b) obtained from the simulations of the HeH + model system at B = 10 B 0 using different propagators with K = 1 (Velocity Verlet): Auxiliary coordinates and momenta (ACM), Tajima (TAJ), and exponential (EXP) propagator.

FIG. 4 :
FIG. 4: Influence of the order of the EXP propagator (K, a) and the truncation of the exponential series (N , b) on the standard deviation of the total energy (ε tot ) during simulations of the HeH + model system using different magnetic field strengths.The reference (exact exponential with K = 1) is shown as a dashed black line in both plots.

TABLE I :
Comparison of the auxiliary coordinates and momenta (ACM), Tajima (TAJ), and exponential (EXP) propagator.

TABLE II :
Coefficients a and b for the different integrators of order K used in this work.