June 04, 2026
We derive the first-principles Ehrenfest molecular dynamics describing non-adiabatic processes with the inclusion of the nuclear-velocity-dependent phases (also known as electron-translation factors) on the atomic-orbital basis. These phases, appearing when nuclei are treated dynamically, affect effective Hamiltonians constructed from localised orbitals. In this work, we focus on the effects in the first-principles pseudo-potential Hamiltonian, both for the norm-conserving and ultra-soft cases, derived within the Projector-Augmented-Wave (PAW) method framework. Peierls-like phases depending on the nuclear velocities appear in the non-local part of the potential, while additional nuclear velocity and acceleration-dependent corrections appear in the ultra-soft pseudo-potential case. The use of velocity-including atomic orbital basis enables a Galilean-invariant description of the non-adiabatic Ehrenfest molecular dynamics, removing spurious non-adiabatic couplings that arise from neglecting the nuclear velocity phases in the atomic orbitals.
Effective electronic Hamiltonians are widely used in condensed matter physics to describe the low-energy spectrum of the system and the related physical properties [1]. Within the ab-initio framework, the use of pseudo-potentials enables the achievement of large computational gains with respect to all-electron calculations [2]–[9]. In the Born-Oppenheimer approximation, effective electronic Hamiltonians are usually derived starting from the all-electron one, at fixed nuclei, and excluding non-adiabatic effects. Even when non-adiabatic effects are included in the nuclear dynamics, the effective electronic Hamiltonians are still usually constructed using atomic-like orbitals as if the nuclei were at rest, or at most considering the electronic wavefunction rigidly shifting with the nuclei. This approach introduces discrepancies between the dynamics of the effective models and of the all-electron Hamiltonian. For example, in the context of the non-adiabatic molecular dynamics, spurious couplings arise in the effective Hamiltonian governing the temporal evolution of the electronic systems [10]–[13]. They manifest in the paradox that a motion at constant velocity of an atom can induce electronic transitions [10], [11], clearly breaking the principle of Galilean invariance. Analogously, the breaking of Galilean invariance manifests in the violation of nonadiabatic sum rules relating Born effective charges and the interatomic force constant matrix to the frequency-dependent electromagnetic susceptibility, when calculated with pseudo-potential Hamiltonians obtained rigidly shifting the electronic orbitals with the nuclei, as opposed to all-electron ones, as recently pointed out in Ref. [14].
These issues arise from neglecting the nuclear-velocity dependent phases in the atomic-like orbitals [10], [10], [11], [11], [14]–[16], often called "electron-translation factors (ETFs)" in the literature. The use of nuclear-velocity-including atomic orbital basis has a long history in the study of atomic collisions, where it was first introduced [17] and saw several developments [15], [18]–[25], [25]–[27]. In the study of the properties of molecules and solids, velocity-including atomic orbitals were first introduced by Nafie to assess vibrational circular dichroism (VCD) in Ref. [28] using the complete adiabatic approach [29]. Nuclear velocity-dependent phases enable non-zero electronic currents and magnetic dipole moments, which, conversely, vanish for the standard real-valued Born-Oppenheimer electronic wavefunctions, resulting in a vanishing electronic contribution to VCD [30], [31]. In this framework, nuclear velocity perturbation theory (NVPT) has been developed and implemented, mainly for the purpose of determining VCD spectra [32]–[35]. It is worth noticing that the exact factorisation method [36] encodes these effects, as it can be appreciated in the classical limit for the nuclear wavefunction [30].
Besides being essential for correctly describing VCD spectra, the use of nuclear-velocity-including atomic orbital basis allows for the removal of spurious electronic transitions in the dynamics. This has been shown, at the linear order in the nuclear velocity, in the context of a time-dependent Hartree-Fock dynamics in Refs. [12], [13], [37], [38] and for localised atomic orbitals in Refs. [10], [11], [16], as well as in the context of TDDFT linear response theory in Refs. [39], [40]. Furthermore, they restore the frequency-dependent sum rules in the framework of pseudopotential approximation, correctly accounting for nonadiabatic couplings and frequency-dependent vibrational responses within linear-response theory [14]. Nonetheless, nuclear velocity effects on the wavefunctions are still neglected in current implementations of first-principles non-adiabatic molecular dynamics [41]–[48].
In this work, we determine the first-principles non-adiabatic semiclassical Ehrenfest dynamics - where nuclei are classical and electrons are quantum mechanical particles [49] - within the Projector Augmented Wave (PAW) method framework for both the ultra-soft and norm-conserving cases, using nuclear-velocity-including atomic orbitals. Our approach enables the removal of the spurious contribution in the non-adiabatic Ehrenfest dynamics, recovering the physical results of the all-electron dynamics with the pseudopotential one, where full Galilean invariance is respected. These nuclear velocity dependences manifest as Peierls-like phases appearing in the non-local part of the pseudopotential, capturing the response at any order in the nuclear velocity, as already discussed in Ref. [14]. We generalise the results of Ref. [14] to the case of ultrasoft pseudo-potentials and we consider the contributions from the temporal derivatives of the PAW transformation. They originate additional nuclear velocity and acceleration dependent terms, relevant in the ultrasoft case, always neglected in the previous literature.
The paper is organised as follows: in Sec. 2 we introduce the velocity-including atomic orbitals, studying the electronic problem for an isolated nucleus in motion, while in Sec. 3 we present the semiclassical Ehrenfest Lagrangian formalism used to derive the effective Hamiltonians and the equations of motion for velocity-including atomic orbitals. The standard PAW method is explained in Section 4, showing the adiabatic Hamiltonian and the related responses. These instruments allow us to review the currently implemented non-adiabatic PAW Ehrenfest dynamics [41]–[43] in 4.2. Conversely, Section 5 is devoted to explaining the PAW method with the velocity-including atomic orbitals, which enables the proper description of the non-adiabatic Ehrenfest dynamics 5.1. Finally, we draw our conclusions in Section 6. The Hamiltonians, the electronic temporal dynamics and the conserved energies are compared with the literature in Table [tab:hamiltonians95dynamics], that summarises the main findings of this work.
Localised atomic orbital bases are usually constructed from the orbitals of isolated atoms with the nucleus fixed in its position. In this paragraph, we establish a localised atomic orbital basis for the case of moving nuclei.
Consider an isolated atom \(s\) with the nucleus at rest, which, without loss of generality, is located at the origin. The electrons are described by a single particle all-electron Hamiltonian \[\hat{H}_{\mathbf{0}}^{\mathrm{AE}}(\hat{\mathbf{r}},\hat{\mathbf{p}})=\frac{\hat{\mathbf{p}}^2}{2m}+V_s(\hat{\mathbf{r}}). \label{eq:all-electron95hamiltonian95static}\tag{1}\] where the hat \(\hat{}\) denotes operators, \(V_s(\hat{\mathbf{r}})\) is the effective potential acting on the electron, the subscript \(s\) indicates the dependence of the potential on the atomic properties of \(s\). The Hamiltonian is diagonalised by the atomic orbitals \(\ket{\phi^{\rm \mathbf{0}}_{si}}\) with energy \(E_{si}\), where \(i\) indicates the electronic quantum number and the superscript \(\mathbf{0}\) indicates that the orbitals \(\ket{\phi^{\rm \mathbf{0}}_{si}}\) are solutions for the atom at rest at the origin, \[\hat{H}_{\mathbf{0}}^{\mathrm{AE}}(\hat{\mathbf{r}},\hat{\mathbf{p}})\ket{\phi^{\rm \mathbf{0}}_{si}}=E_{si}\ket{\phi^{\rm \mathbf{0}}_{si}}.\] Usually, localised atomic orbital bases are constructed using \(\ket{\phi^{\rm \mathbf{0}}_{si}}\) orbitals, translated to the positions of atoms of the system \[\ket{\phi^{\mathbf{R}_{s}}_{si}}=\hat{T}_{\mathbf{R}_{s}}\ket{\phi^{\mathbf{0}}_{si}},\] that form the basis set \(\{\ket{\phi^{\mathbf{R}_{s}}_{si}}\}_{si}\). We remark that the translation operator \(\hat{T}_{\mathbf{R}_{s}}\) acts on the position eigenstates as \(\hat{T}_{\mathbf{R}_s}\ket{\mathbf{r}}=\ket{\mathbf{r}+\mathbf{R}_s}\). The rigid translation of the orbitals assumes that nuclei are static. If the nuclei are moving, with their coordinate explicitly depending on time \(\mathbf{R}_{s}(t)\), the rigidly translated atomic orbitals, in spite of their common use, should be replaced by atomic orbitals accounting for the nuclear motion. In the following, we derive this basis set.
Suppose for \(t<0\) the \(s\) nucleus, identified by the coordinate \(\mathbf{R}_s(t)\), is at rest (\(\mathbf{R}_s(t<0)=\boldsymbol{0}\)) and that, at \(t=0\), the nucleus starts moving with velocity \(\dot{\mathbf{R}}_s(t)\) and acceleration \(\ddot{\mathbf{R}}_s(t)\). At \(t=0\) \(\ket{\phi_{si}(t)}=\ket{\phi^{\rm \mathbf{0}}_{si}}\), then the time evolution of the atomic orbital \(\ket{\phi_{si}(t)}\) for \(t>0\) is determined by the Schrödinger equation with the Hamiltonian \(\hat{H}^{\mathrm{AE}}(\hat{\mathbf{r}},\hat{\mathbf{p}};\mathbf{R}_s(t))\), that depends on time through the nuclear position \(\mathbf{R}_s(t)\), \[\begin{align} &i\hbar \frac{d\ket{\phi_{si}(t)}}{dt}=\hat{H}^{\mathrm{AE}}(\hat{\mathbf{r}},\hat{\mathbf{p}};\mathbf{R}_s(t)) \ket{\phi_{si}(t)},\tag{2}\\ & \hat{H}^{\mathrm{AE}}(\hat{\mathbf{r}},\hat{\mathbf{p}};\mathbf{R}_s(t))=\frac{\hat{\mathbf{p}}^2}{2m}+V_s(\hat{\mathbf{r}}-\mathbf{R}_s(t)). \tag{3} \end{align}\] For \(t>0\), we suppose that the solution of the above Schrödinger equation (Eq. 2 ) is expressed in the form of \[\begin{align} \ket{\phi_{si}(t)}&=e^{i\alpha_{s}(\hat{\mathbf{r}})}\hat{T}_{\mathbf{R}_{s}(t)}e^{i\Theta(t)}\ket{\phi'_{si} (t)},\label{eq:atomic95orbital95I}\\ \alpha_{s}(\hat{\mathbf{r}})&=\frac{m}{\hbar}\dot{\mathbf{R}}_s(t)\cdot(\hat{\mathbf{r}}-\mathbf{R}_s(t)), \\ \Theta(t)&=-\frac{1}{\hbar}\left(E_{si}t-\int_0^t dt' \frac{m|\dot{\mathbf{R}}_s(t')|^2}{2}\right). \nonumber \end{align}\tag{4}\] We remark that the velocity-dependent phase \(\alpha_{s}(\hat{\mathbf{r}})\) does not depend on the choice of the origin of the reference frame. Its time dependence is instantaneous. Conversely, the purely time-dependent phase \(\Theta(t)\) and \(\ket{\phi'_{si} (t)}\) depend on the system’s history. In particular, the phase \(\Theta(t)\) depends on the temporal integration of the electron’s kinetic energy moving at the velocity of the nucleus.
To obtain the differential equation for \(\ket{\phi'_{si} (t)}\), we plug Eq. 4 into Eq. 2 . The temporal derivative on the left-hand side is \[\begin{align} &\frac{d \ket{\phi_{si}(t)}}{dt} = e^{i\alpha_{s}(\hat{\mathbf{r}})}\hat{T}_{\mathbf{R}_{s}(t)}e^{i\Theta(t)}\Bigg[\frac{d \ket{\phi'_{si}(t)}}{dt}-\frac{i}{\hbar}\Bigg(E_{si}\\&+\dot{\mathbf{R}}_s(t)\cdot \hat{\mathbf{p}} -m \ddot{\mathbf{R}}_s(t)\cdot \hat{\mathbf{r}} +\frac{m|\dot{\mathbf{R}}_s(t)|^2}{2}\Bigg)\ket{\phi'_{si}(t)}\Bigg] \label{eq:derivative95velocity95including} \end{align}\tag{5}\] where we used that \(\hat{T}^{\dagger}_{\mathbf{R}_{s}(t)} \hat{\mathbf{r}}\hat{T}_{\mathbf{R}_{s}(t)}=\hat{\mathbf{r}}+\mathbf{R}_s(t)\). On the right-hand side of Eq. 2 , the application of the Hamiltonian to \(\ket{\phi_{si}(t)}\) corresponds to \[\begin{align} & \hat{H}^{\mathrm{AE}}(\hat{\mathbf{r}},\hat{\mathbf{p}};\mathbf{R}_s(t)) e^{i\alpha_{s}(\hat{\mathbf{r}})}\hat{T}_{\mathbf{R}_{s}(t)} \\ =&e^{i\alpha_{s}(\hat{\mathbf{r}})} \hat{T}_{\mathbf{R}_{s}(t)}\hat{H}^{\mathrm{AE}}(\hat{\mathbf{r}}+\mathbf{R}_s(t),\hat{\mathbf{p}}+m\dot{\mathbf{R}}_s(t);\mathbf{R}_s(t)) . \label{eq:Hamiltonian95isolated95atom95translated} \end{align}\tag{6}\] where we used \(e^{-i\alpha_{s}(\hat{\mathbf{r}})}\hat{\mathbf{p}}e^{i\alpha_{s}(\hat{\mathbf{r}})}=\hat{\mathbf{p}}+m\dot{\mathbf{R}}_s(t)\). Conversely, the purely time-dependent phase factor commutes with the Hamiltonian. By putting together Eq. 5 and Eq. 6 , we obtain that the Schrödinger equation for moving nuclei of Eq. 2 is equivalent to the following equations for \(\ket{\phi'_{si}(t)}\) with the boundary condition \(\ket{\phi'_{si}(t=0)}=\ket{\phi^{\rm \mathbf{0}}_{si}}\) \[\begin{align} &i\hbar \frac{d\ket{\phi'_{si}(t)}}{dt}=\left(\hat{H}'^{\mathrm{AE}}(\hat{\mathbf{r}}, \hat{\mathbf{p}};\ddot{\mathbf{R}}_s(t))-E_{si}\right) \ket{\phi'_{si}(t)},\tag{7}\\ &\hat{H}'^{\mathrm{AE}}(\hat{\mathbf{r}}, \hat{\mathbf{p}};\ddot{\mathbf{R}}_s(t))=\frac{\hat{\mathbf{p}}^2}{2m}+V_s(\hat{\mathbf{r}})+m\ddot{\mathbf{R}}_s(t)\cdot \hat{\mathbf{r}}\tag{8}. \end{align}\] \(\hat{H}'^{\mathrm{AE}}(\hat{\mathbf{r}}, \hat{\mathbf{p}};\ddot{\mathbf{R}}_s(t))\) differs from the \(\hat{H}_{\mathbf{0}}^{\mathrm{AE}}(\hat{\mathbf{r}}, \hat{\mathbf{p}})\) only by the term \(m\ddot{\mathbf{R}}_s(t)\cdot \hat{\mathbf{r}}\), that can be interpreted as the non-inertial force acting on the electrons in the frame where the nucleus is at rest. Its effect on the Hamiltonian is analogous to a time-dependent electric field, described by the Stark effect.
By neglecting the nuclear acceleration term \(m\ddot{\mathbf{R}}_s(t)\cdot \hat{\mathbf{r}}\) in the Hamiltonian \(\hat{H}'^{\mathrm{AE}}(\hat{\mathbf{r}}, \hat{\mathbf{p}};\ddot{\mathbf{R}}_s(t))\), \(\ket{\phi'_{si}(t)}=\ket{\phi^{\rm \mathbf{0}}_{si}}\), implying that \(\ket{\phi'_{si}(t)}\) does not depend on the history of the system and that the atomic orbitals for static nuclei \(\ket{\phi^{\rm \mathbf{0}}_{si}}\) satisfy Eq. 7 .
Therefore, for moving nuclei the atomic orbitals should be expressed in the form of Eq. 4 . The purely time-dependent phases, including those depending on the history of the system, do not affect the basis and can be removed. Consequently, we define the velocity-including localised atomic orbital as \[\begin{align} \ket{\phi^{\dot{\mathbf{R}}_s,\mathbf{R}_{s}}_{si}}=e^{i\frac{m}{\hbar}\dot{\mathbf{R}}_s(t)\cdot(\hat{\mathbf{r}}-\mathbf{R}_s(t))}\ket{\phi^{\mathbf{R}_{s}(t)}_{si}}.\label{eq:vi95atomic95orbital} \end{align}\tag{9}\] In the superscript \(\dot{\mathbf{R}}_s,\mathbf{R}_{s}\) we omit the temporal dependence for brevity since it is clear from the context. The velocity-including basis set \(\{{\ket{\phi^{\dot{\mathbf{R}}_s,\mathbf{R}_{s}}_{si}}}\}_{si}\) depend on time instantaneously through the nuclear position \(\mathbf{R}_s(t)\) and velocity \(\dot{\mathbf{R}}_s(t)\). The temporal derivative of the states is \[\begin{align} &\frac{d\ket{\phi^{\dot{\mathbf{R}}_s,\mathbf{R}_{s}}_{si}}}{dt}=\frac{i}{\hbar}e^{i\alpha_s(\hat{\mathbf{r}})}\Bigg(m\ddot{\mathbf{R}}_s(t)\cdot(\hat{\mathbf{r}}-\mathbf{R}_s(t))\\&-\dot{\mathbf{R}}_s(t)\cdot \hat{\mathbf{p}}-m|\dot{\mathbf{R}}_s(t)|^2\Bigg) \ket{\phi^{\mathbf{R}_{s}(t)}_{si}},\label{eq:derivative95vi95atomic95orbital} \end{align}\tag{10}\] and, by using that \(e^{i\alpha_{s}(\hat{\mathbf{r}})}\hat{\mathbf{p}}e^{-i\alpha_{s}(\hat{\mathbf{r}})}=\hat{\mathbf{p}}-m\dot{\mathbf{R}}_s(t)\), \[\begin{align} \frac{d\ket{\phi^{\dot{\mathbf{R}}_s,\mathbf{R}_{s}}_{si}}}{dt}=&\frac{i}{\hbar}\Bigg(m\ddot{\mathbf{R}}_s(t)\cdot(\hat{\mathbf{r}}-\mathbf{R}_s(t))\\&-\dot{\mathbf{R}}_s(t)\cdot \hat{\mathbf{p}}\Bigg) \ket{\phi^{\dot{\mathbf{R}}_s,\mathbf{R}_{s}}_{si}}.\label{eq:derivative95vi95atomic95orbital95v2} \end{align}\tag{11}\] In the construction of the localised atomic orbital basis set, neglecting the nuclear acceleration term, which induces a Stark effect, causes a small error. Indeed, we are excluding the mixing, due to the Stark effect, of the occupied states with orbitals that are not included in the basis set. These contributions are small because of the large energy difference between the occupied states and those excluded from the basis set, becoming increasingly smaller as the basis set is enlarged.
The notation used to indicate the atomic orbitals is summarised in Table 1. From now on, we omit the \(\hat{}\) unless it is important for context.
| Atomic orbital for level \(i\) of atom \(s\) | |
|---|---|
| \(\ket{\phi^{\mathbf{0}}_{si}}\) | for the atom centred at the origin |
| \(\ket{\phi^{\mathbf{R}_{s}}_{si}}\) | translated to the fixed atomic equilibrium position |
| \(\ket{\phi^{\mathbf{R}_{s}(t)}_{si}}\) | translated to the time-dependent atomic position |
| \(\ket{\phi^{\dot{\mathbf{R}}_s,\mathbf{R}_{s}}_{si}}\) | \(\ket{\phi^{\mathbf{R}_{s}(t)}_{si}}\) times nuclear velocity-dependent phase |
Consider a system of classical nuclei with mass \(M_s\), located at the positions \(\{\mathbf{R}_s(t)\}\), and quantum electrons. The all-electron single-particle mean-field Hamiltonian \(\hat{H}\) presents an effective self-consistent local potential \(V(\hat{\mathbf{r}})\), obtained with a density functional theory (DFT) local or semi-local approximation for the exchange-correlation functional. We do not explicitly treat the self-consistent potential since it is local, adding no complications compared with the non-self-consistent case for our purposes. Therefore, the single-particle Hamiltonian of the electron interacting with many nuclei is \[\hat{H}^{\mathrm{AE}}(\hat{\mathbf{r}},\hat{\mathbf{p}};\{\mathbf{R}_s(t)\})=\frac{\hat{\mathbf{p}}^2}{2m}+V(\hat{\mathbf{r}}). \label{eq:HAEsingle}\tag{12}\] At zero temperature, the system of electrons and nuclei can be described through the real-valued Ehrenfest Lagrangian [49], depending on the independent variables \(\mathbf{q}=\left(\{\mathbf{R}_s(t)\}_s, \{\ket{\psi_I(t)}\}_I, \{\bra{\psi_I(t)}\}_I\right)\), \[\begin{align} & \mathcal{L}(\mathbf{q})=\sum_s\frac{M_s|\dot{\mathbf{R}}_s(t)|^2}{2}-\sum_{I=1}^{N_{\rm el}}\Bigg(\braket{\psi_I(t)|\hat{H}^{\mathrm{AE}}|\psi_I(t)}\\&-\frac{i\hbar}{2}\left(\bra{\psi_I(t)} \frac{d\ket{\psi_I(t)}}{dt}- \frac{d\bra{\psi_I(t)}}{dt}\ket{\psi_I(t)}\right)\Bigg) \end{align} \label{eq:lagrangian95general}\tag{13}\] where \(\ket{\psi_I(t)}\) are the solutions of Eq. 12 , and \(N_{\rm el}\) is the number of electrons. In addition, if we apply a transformation that depends on time to the wavefunction \[\ket{\psi'(t)}=\mathcal{R}(t)\ket{\psi(t)},\] the Ehrenfest Lagrangian of Eq. 13 is transformed to
\[\begin{align} \mathcal{L}'\left(\mathbf{q}'\right)=\sum_s\frac{M_s|\dot{\mathbf{R}}_s(t)|^2}{2}-&\sum_{I=1}^{N_{\rm el}}\Bigg[\braket{\psi'_I(t)|\mathcal{R}^{\dagger}(t)\hat{H}^{\rm AE}\mathcal{R}(t)|\psi'_I(t)}+\frac{i\hbar}{2}\Bigg(\braket{\psi'_I(t)| \frac{d\mathcal{R}^{\dagger}(t)}{dt}\mathcal{R}(t)-\mathcal{R}^{\dagger}(t)\frac{d\mathcal{R}(t)}{dt}| \psi'_I(t)}\\&+\frac{d\bra{\psi'_I(t)}}{dt}\mathcal{R}^{\dagger}(t)\mathcal{R}(t)\ket{\psi'_I(t)}-\bra{\psi'_I(t)}\mathcal{R}^{\dagger}(t)\mathcal{R}(t)| \frac{d\ket{\psi'_I(t)}}{dt}\Bigg)\Bigg], \end{align} \label{eq:transformed95lagrangian95general}\tag{14}\]
where the independent variables are \(\mathbf{q}'=\left(\{\mathbf{R}_s(t)\}_s, \{\ket{\psi'_I(t)}\}_I, \{\bra{\psi'_I(t)}\}_I\right)\). We do not treat the finite temperature case directly in the Lagrangian, but it can be obtained with straightforward generalisations. In the following sections, we use Eq. 14 to determine the Ehrenfest Lagrangian for the Projector Augmented Wave (PAW) method with the inclusion of the nuclear velocity phases in the atomic orbitals.
The equations of motion for the nuclear and electronic variables are determined by imposing the action \(\mathcal{A}=1/T\int_0^T \mathcal{L}dt\) to be stationary for a variation of an independent variable of the Lagrangian, i.e. \(\frac{\delta \mathcal{A}}{\delta q_i}=0\). Usually, the Lagrangian depends on \(\mathbf{q} \textrm{ and } \dot{\mathbf{q}}\), implying that the equations of motion are the Euler-Lagrange equations \[\frac{\partial \mathcal{L}(\mathbf{q},\dot{\mathbf{q}})}{\partial q_i}-\frac{d}{dt}\frac{\partial \mathcal{L}(\mathbf{q},\dot{\mathbf{q}})}{\partial \dot{q}_i}=0, \label{eq:Euler-Lagrange}\tag{15}\] which are always used in this work for the electronic degrees of freedom (wavefunctions). The conserved energy is obtained as the Legendre transformation of the Lagrangian \[E=\sum_{i}\frac{\partial \mathcal{L}(\mathbf{q},\dot{\mathbf{q}})}{\partial \dot{q}_i}\dot{q}_i-\mathcal{L}(\mathbf{q},\dot{\mathbf{q}}). \label{eq:conserved95energy}\tag{16}\] Conversely, the nuclear-velocity dependence of the atomic orbitals may lead to a nuclear-acceleration dependence in the effective model’s Lagrangian. Consequently, the Euler-Lagrange equations need to be generalised to higher orders to describe nuclear dynamics. Under the assumption that the highest order temporal derivative is the second order one, the Euler-Lagrange equations are generalised as [50]–[52] \[\frac{\partial \mathcal{L}(\mathbf{q},\dot{\mathbf{q}},\ddot{\mathbf{q}})}{\partial q_i}-\frac{d}{dt}\frac{\partial \mathcal{L}(\mathbf{q},\dot{\mathbf{q}},\ddot{\mathbf{q}})}{\partial \dot{q}_i}+\frac{d^2}{dt^2}\frac{\partial \mathcal{L}(\mathbf{q},\dot{\mathbf{q}},\ddot{\mathbf{q}})}{\partial \ddot{q}_i}=0. \label{eq:Euler-Lagrange95general}\tag{17}\] Specifically, the equations for the nuclear system are obtained by setting \(q_i \to \mathbf{R}_s\). Explicitly, if the Lagrangian depends also on the nuclear acceleration, the forces governing the nuclear dynamics are \[M_s \ddot{\mathbf{R}}_s(t)=\frac{\partial \mathcal{L}_{\mathrm{el}}}{\partial \mathbf{R}_s}-\frac{d}{dt}\frac{\partial \mathcal{L}_{\mathrm{el}}}{\partial \dot{\mathbf{R}}_s}+\frac{d^2}{dt^2}\frac{\partial \mathcal{L}_{\mathrm{el}}}{\partial \ddot{\mathbf{R}}_s} \label{eq:forces95general}\tag{18}\] where \(\displaystyle \mathcal{L}_{\mathrm{el}}=\mathcal{L}-\sum_s\frac{M_s|\dot{\mathbf{R}}_s(t)|^2}{2}\). If Lagrangian depends linearly on the acceleration, the higher order derivative in the equations of motion is the second order one, avoiding issues of Ostrogradsky’s instabilities. In the presence of a dependence on second order derivatives of the variables \(\mathbf{q}\) in the Lagrangian, the energy is obtained through Ostrogradsky’s construction [50]–[52] as \[\begin{align} E=\sum_{i}\frac{\partial \mathcal{L}(\mathbf{q},\dot{\mathbf{q}},\ddot{\mathbf{q}})}{\partial \ddot{q}_i}\ddot{q}_i-\mathcal{L}(\mathbf{q},\dot{\mathbf{q}},\ddot{\mathbf{q}})\nonumber\\ +\sum_{i}\Bigg[\left(\frac{\partial \mathcal{L}(\mathbf{q},\dot{\mathbf{q}},\ddot{\mathbf{q}})}{\partial \dot{q}_i}-\frac{d}{dt}\frac{\partial \mathcal{L}(\mathbf{q},\dot{\mathbf{q}},\ddot{\mathbf{q}})}{\partial \ddot{q}_i}\right)\dot{q}_i\Bigg]. \label{eq:conserved95energy95Ostrograsky} \end{align}\tag{19}\] As for the equations of motion, the linearity of the Lagrangian on \(\ddot{\mathbf{q}}\) forbids any derivatives of order higher than the second in the conserved energy.
In the absence of an explicit time-dependence in the Lagrangian, the total energy is conserved because the time-derivative of \(E\) is zero along the trajectory of the system. All the Lagrangians presented in this work do not have an explicit time dependence. This implies that the energy, either obtained in the usual Hamiltonian formulation of Eq. 16 or with the Ostrogradsky generalisation of Eq. 19 , is always conserved.
The Projector Augmented Wave (PAW) method, developed in Ref. [6], determines the single-particle pseudo-wavefunction \(\ket{\tilde{\psi}}\) through a linear transformation \(\hat{\mathcal{T}}\) applied to the all-electron wavefunction \[\ket{\psi}=\hat{\mathcal{T}}\ket{\tilde{\psi}}.\] In the original formulation, using atomic orbitals for nuclei at rest centred at the positions of the nuclei, the linear transformation is defined by selecting a set of all-electron partial waves \(\ket{\phi^{\mathbf{R}_s}_{si}}\) obtained by applying the transformation \(\hat{\mathcal{T}}\) to the set of pseudo-partial waves \(\ket{\tilde{\phi}^{\mathbf{R}_s}_{si}}\) \[\begin{align} \hat{\mathcal{T}}=\hat{\mathbb{1}}+\sum_{s}\hat{t}^s,\\ \hat{t}^s=\sum_i\left(\ket{\phi^{\mathbf{R}_s}_{si}}- \ket{\tilde{\phi}^{\mathbf{R}_s}_{si}}\right)\bra{\tilde{p}^{\mathbf{R}_s}_{si}} \end{align} \label{eq:PAW95transformation}\tag{20}\] where the projectors \(\ket{\tilde{p}^{\mathbf{R}_s}_{si}}\) satisfy the orthogonality relation with the pseudo-partial waves \(\braket{\tilde{p}^{\mathbf{R}_s}_{si}|\tilde{\phi}^{\mathbf{R}_{s'}}_{s'j}}=\delta_{ss'}\delta_{ij}\). We further assume that there exists an augmentation region \(\Omega_s\) around each site \(\mathbf{R}_s\), where the all-electron partial waves \(\ket{\phi^{\mathbf{R}_s}_{si}}\) form a complete set for the valence all-electron wavefunction. Outside the region, the all-electron and pseudo-partial waves coincide and the pseudo-projectors vanish. Finally, the augmentation regions do not overlap. For a complete basis set, if \(\mathbf{r} \in \Omega_s\), the following relations hold \[\begin{align} \sum_{i}\braket{\mathbf{r}|\tilde{\phi}^{\mathbf{R}_{s}}_{si}}\braket{\tilde{p}^{\mathbf{R}_s}_{si}|\tilde{\psi}}=\braket{\mathbf{r}|\tilde{\psi}}, \\ \sum_{i}\braket{\tilde{\psi}|\tilde{p}^{\mathbf{R}_{s}}_{si}}\braket{\tilde{\phi}^{\mathbf{R}_s}_{si}|\mathbf{r}}=\braket{\tilde{\psi}|\mathbf{r}}. \end{align} \label{eq:identity95augmentation95region}\tag{21}\]
Any local or semi-local operator transforms as [6] \[\begin{align} \tilde{O}=\hat{\mathcal{T}}^{\dagger}O\hat{\mathcal{T}}= O+\sum_{s,ij} \ket{\tilde{p}^{\mathbf{R}_s}_{si}} \Delta O^s_{ij}\bra{\tilde{p}^{\mathbf{R}_s}_{sj}}, \tag{22}\\ \Delta O^s_{ij}=\braket{\phi^{\mathbf{R}_s}_{si}|O|\phi^{\mathbf{R}_s}_{sj}}-\braket{\tilde{\phi}^{\mathbf{R}_s}_{si}|O|\tilde{\phi}^{\mathbf{R}_s}_{sj}}, \tag{23} \end{align}\] where the second equality requires the use of the identity in the augmentation region as given in Eq. 21 , which holds for a complete basis set. The transformation rule for the operators of Eq. 22 implies that the identity operator in the all-electron wavefunction space transforms as \[\begin{align} S=\hat{\mathcal{T}}^{\dagger}\hat{\mathcal{T}}=&\mathbb{1}+\sum_{s,ij}\ket{\tilde{p}^{\mathbf{R}_s}_{si}}Q^s_{ij}\bra{\tilde{p}^{\mathbf{R}_s}_{sj}}\\ Q^s_{ij}=&\braket{\phi^{\mathbf{R}_s}_{si}|\phi^{\mathbf{R}_s}_{sj}}-\braket{\tilde{\phi}^{\mathbf{R}_s}_{si}|\tilde{\phi}^{\mathbf{R}_s}_{sj}}, \end{align} \label{eq:TdaggerT95identity}\tag{24}\] where \(\mathbb{1}\) indicates the identity in the pseudo-wavefunction space. As a consequence, the pseudo-wavefunction are normalised as \(\braket{\psi|S|\psi}=1\). In the norm-conserving case, where \(\braket{\phi^{\mathbf{R}_s}_{si}|\phi^{\mathbf{R}_s}_{sj}}=\braket{\tilde{\phi}^{\mathbf{R}_s}_{si}|\tilde{\phi}^{\mathbf{R}_s}_{sj}}\), the identity in the all-electron wave-function space is mapped into the identity in the pseudo-wavefunction space \(S=\mathbb{1}\). The PAW Hamiltonian is analogously obtained as \[\begin{align} \hat{H}^{\rm PAW}=&\hat{\mathcal{T}}^{\dagger}\hat{H}^{\rm AE}\hat{\mathcal{T}} = \frac{\mathbf{p}^2}{2m}+V^{\mathrm{loc}}+\sum_{s}v_{s}^{\mathrm{nl}}. \end{align} \label{eq:PAW95Hamiltoniana95standard}\tag{25}\] \(V^{\mathrm{loc}}\) is the local part of the pseudo-potential, plus eventually a DFT self-consistent potential; \(v_{s}^{\mathrm{nl}}\) is the non-local part of the potential obtained as in Ref. [6]
In the adiabatic limit, the electronic wavefunction remains in the ground state, satisfying the eigenvalue equation [5], [6], [53], [54], where \(H^{\mathrm{PAW}}\) and \(S\) are time-dependent through the nuclear motion as well as the instantaneous eigenvalues and eigenvectors, \[H^{\mathrm{PAW}}\ket{\tilde{\psi}_I}=E_IS\ket{\tilde{\psi}_I}\]
The forces on the ions are [53], [55]–[57] \[M_s\ddot{\mathbf{R}}_s(t)=\mathbf{F}_s^{\mathrm{HF}}+\mathbf{F}_s^{\mathrm{S}}\] \[\begin{align} &\mathbf{F}_s^{\mathrm{HF}}=-\sum_{I=1}^{N_{\rm el}}\braket{\tilde{\psi}_I|\frac{\partial H^{\mathrm{PAW}}}{\partial \mathbf{R}_s}|\tilde{\psi}_I}, \end{align} \label{eq:HF95force95adiabatic}\tag{26}\] \[\begin{align} &\mathbf{F}_s^{S}=\sum_{I=1}^{N_{\rm el}}E_I\bra{\tilde{\psi}_I}\frac{\partial S}{\partial \mathbf{R}_s}\ket{\tilde{\psi}_I}. \end{align} \label{eq:S95force95adiabatic}\tag{27}\] where the first contribution is the standard Hellmann-Feynman force \(\mathbf{F}_s^{\mathrm{HF}}\), i.e. the quantum mechanical average over the electronic ground state of the variation of the electronic Hamiltonian due to the nuclear displacement; the second one \(\mathbf{F}_s^{S}\) is the force originating from the transformation of the identity operator in the pseudo-wavefunction space. In the case of a norm-conserving pseudo-potential \(\mathbf{F}_s^{S}\) vanishes, leaving only the \(\mathbf{F}_s^{\mathrm{HF}}\) contribution to the forces.
In this section, we present the equations of the PAW non-adiabatic Ehrenfest dynamics obtained with the electronic wavefunction rigidly shifting with the nuclei. These equations are presented in Refs. [41], [42], [58] and currently implemented in the GPAW code [43].
The Ehrenfest Lagrangian is obtained using Eq. 14 with the transformation \(\hat{\mathcal{T}}\), defined in Eq. 20 , with the nuclear positions explicitly depending on time \(\mathbf{R}_s(t)\), \[\begin{align} &\mathcal{L}^{\mathrm{PAW}}_{\rm R}=\sum_s\frac{M_s|\dot{\mathbf{R}}_s(t)|^2}{2}-\sum_{I=1}^{N_{\rm el}}\Bigg[\braket{\tilde{\psi}_I(t)|\hat{H}^{\mathrm{PAW}}_{\rm R}|\tilde{\psi}_I(t)}\\&+\frac{i\hbar}{2}\left(\frac{d\bra{\tilde{\psi}_I(t)}}{dt}S\ket{\tilde{\psi}_I(t)}-\bra{\tilde{\psi}_I(t)}S| \frac{d\ket{\tilde{\psi}_I(t)}}{dt}\right)\Bigg], \end{align} \label{eq:PAW95lagrangian95NOVI}\tag{28}\] where \[\begin{align} \hat{H}^{\mathrm{PAW}}_{\rm R}&=\hat{H}^{\rm PAW}+\frac{i\hbar}{2} \left(\frac{d\hat{\mathcal{T}}^{\dagger}}{dt}\hat{\mathcal{T}}-\hat{\mathcal{T}}^{\dagger}\frac{d\hat{\mathcal{T}}}{dt}\right). \label{eq:paw95hamiltonian95R} \end{align}\tag{29}\] The \(S\) matrix is here time dependent, differently from the adiabatic case, due to the time dependence of the nuclear positions \(\mathbf{R}_s(t)\) to be used inside Eq. 24 . As shown in detail in Appendix 9, the temporal derivative of the transformation, entering in Eq. 29 , is [41], [42] \[\frac{d\hat{\mathcal{T}}}{dt}=- \frac{i}{\hbar} \sum_s \left[\dot{\mathbf{R}}_s(t)\cdot \mathbf{\hat{p}},\hat{t}^s\right], \label{eq:derivative95T95velocity95including}\tag{30}\] implying that \[\begin{align} & \frac{i\hbar}{2} \left(\frac{d\hat{\mathcal{T}}^{\dagger}}{dt}\hat{\mathcal{T}}-\hat{\mathcal{T}}^{\dagger}\frac{d\hat{\mathcal{T}}}{dt}\right)=-\sum_{s}\dot{\mathbf{R}}_s(t) \cdot \mathbf{C}_{s}, \\&\mathbf{C}_{s}=\sum_{ij}\Bigg(-\left\{\frac{\hat{\mathbf{p}}}{2},\ket{\tilde{p}^{\mathbf{R}_s(t)}_{si}}Q^{s}_{ij}\bra{\tilde{p}^{\mathbf{R}_s(t)}_{sj}} \right\}\\ &+\ket{\tilde{p}^{\mathbf{R}_s(t)}_{si}} \Delta \hat{\mathbf{p}}^s_{ij}\bra{\tilde{p}^{\mathbf{R}_s(t)}_{sj}}\Bigg), \label{eq:M95no95vel} \end{align}\tag{31}\] where \(\{,\}\) is the anti-commutator, and \(\Delta\) is used in the sense of Eq. 23 . Therefore, the electronic and nuclear dynamics are obtained from the standard Euler-Lagrange equations (Eq. 15 ), \[i\hbar S\frac{d\ket{\tilde{\psi}_I(t)}}{dt}=\left(H^{\mathrm{PAW}}_{\rm R}-\frac{1}{2}\frac{dS}{dt}\right) \ket{\tilde{\psi}_I(t)}, \label{eq:Schrodinger95PAW95no95vel}\tag{32}\] where the effective Hamiltonian governing the electronic dynamics is \[\begin{align} H^{\mathrm{PAW}}_{\rm R}-\frac{1}{2}\frac{dS(t)}{dt}=H^{\mathrm{PAW}}-i\hbar \hat{\mathcal{T}}^{\dagger}\frac{d\hat{\mathcal{T}}}{dt},\\ -i\hbar \hat{\mathcal{T}}^{\dagger}\frac{d\hat{\mathcal{T}}}{dt}= -\sum_s \dot{\mathbf{R}}_s(t)\cdot \Bigg(\ket{\tilde{p}^{\mathbf{R}_s(t)}_{si}} \Delta \hat{\mathbf{p}}^s_{ij}\bra{\tilde{p}^{\mathbf{R}_s(t)}_{sj}}\\-\ket{\tilde{p}^{\mathbf{R}_s(t)}_{si}}Q^{s}_{ij}\bra{\tilde{p}^{\mathbf{R}_s(t)}_{sj}}\frac{\hat{\mathbf{p}}}{2} \Bigg). \end{align}\] Compared to the adiabatic dynamics, there is the additional time-dependent term \(\displaystyle -i\hbar \hat{\mathcal{T}}^{\dagger}\frac{d\hat{\mathcal{T}}}{dt}\). This term introduces coupling between different electronic states, even when the entire system undergoes a global translation at constant speed, breaking Galilean invariance [11]. As shown in Appendix 7, the presence of time-dependent \(S\) makes it so the evolution of the pseudo-wavefunction is determined by a non-hermitian operator \(\displaystyle H^{\mathrm{PAW}}_{\rm R}-\frac{1}{2}\frac{dS}{dt}\). This allows for the norm conservation of \(\displaystyle \braket{\tilde{\psi}(t)|S|\tilde{\psi}(t)}\) along the temporal evolution.
According to Eq. 18 , the nuclear dynamics is governed by the following equation
\[M_s\ddot{\mathbf{R}}_s(t)=\mathbf{F}_s^{\mathrm{HF}}+\mathbf{F}_s^{\mathrm{HF-\dot{R}}}+\mathbf{F}_s^{\mathrm{S}} \label{eq:Ehrenfest95ions95no95vel}\tag{33}\] where \[\begin{align} \mathbf{F}_s^{\mathrm{HF}}=&-\sum_{I=1}^{N_{\rm el}}\braket{\tilde{\psi}_I(t)|\frac{\partial H^{\mathrm{PAW}}_{\rm R}}{\partial \mathbf{R}_s}|\tilde{\psi}_I(t)},\\ \mathbf{F}_s^{\mathrm{HF-\dot{R}}}=&\sum_{I=1}^{N_{\rm el}}\frac{d}{dt}\braket{\tilde{\psi}_I(t)|\frac{\partial H^{\mathrm{PAW}}_{\rm R}}{\partial \dot{\mathbf{R}}_s}|\tilde{\psi}_I(t)},\\ & \text{with }\\ &\frac{\partial H^{\mathrm{PAW}}_{\rm R}}{\partial \dot{\mathbf{R}}_s}=-\mathbf{C}_s,\\ \end{align} \label{eq:HF95force95pseudo95novel}\tag{34}\]
\[\begin{align} \mathbf{F}_s^{S}=&\frac{i\hbar}{2}\sum_{I=1}^{N_{\rm el}}\Bigg(\bra{\tilde{\psi}_I(t)}\frac{\partial S}{\partial \mathbf{R}_s}\frac{d \ket{\tilde{\psi}_I(t)}}{dt}\\ & -\frac{d\bra{\tilde{\psi}_I(t)}}{dt}\frac{\partial S}{\partial \mathbf{R}_s}\ket{\tilde{\psi}_I(t)}\Bigg). \end{align} \label{eq:HF95force95pseudo95novel95S}\tag{35}\] \(\mathbf{F}_s^{\mathrm{HF}}\) and \(\mathbf{F}_s^{S}\) reduce to Eqs. 26 and 27 in the adiabatic limit where \(\displaystyle H^{\mathrm{PAW}}_{\rm R}\) reduces to \(\displaystyle H^{\mathrm{PAW}}\). In addition, in the non-adiabatic dynamics, an additional Hellmann-Feynman contribution \(\mathbf{F}_s^{\mathrm{HF-\dot{R}}}\) appears due to the nuclear velocity dependence in the electronic part of the Lagrangian, expressed as the time derivative of a Hellmann-Feynman term entailing the nuclear velocity derivative [41].
Finally, according to Eq. 16 , the conserved energy is \[\begin{align} E_{\rm R}^{\mathrm{PAW}}=\sum_s\frac{M_s}{2}|\dot{\mathbf{R}}_s(t)|^2+\sum_{I=1}^{N_{\rm el}}\braket{\tilde{\psi}_I(t)|\hat{H}^{\rm PAW}|\tilde{\psi}_I(t)}, \end{align} \label{eq:conserved95energy95PAW}\tag{36}\] where all the terms that are linear in the derivatives of the independent variables of the Lagrangian cancel.
In the norm-conserving case, where \(S=\mathbb{1}\), \(Q^s_{ij}=0\) and \[\begin{align} \mathbf{C}^{\rm NC}_{s}=\sum_{ij}\ket{\tilde{p}^{\mathbf{R}_s(t)}_{si}} \Delta \hat{\mathbf{p}}^s_{ij}\bra{\tilde{p}^{\mathbf{R}_s(t)}_{sj}}, \label{eq:M95no95vel95NC} \end{align}\tag{37}\] the Hamiltonian reduces to \[H^{\mathrm{PAW-NC}}_{\mathrm{R}}=H^{\mathrm{PAW}}-\sum_{s}\dot{\mathbf{R}}_s(t) \cdot \mathbf{C}^{\rm NC}_{s}. \label{eq:paw95hamiltonian95R95NC}\tag{38}\] The Schrödinger equation for the pseudo-wavefunctions is \[i\hbar \frac{d\ket{\tilde{\psi}_I(t)}}{dt}=H^{\mathrm{PAW-NC}}_{\mathrm{R}}\ket{\tilde{\psi}_I(t)}. \label{eq:Schrodinger95PAW95no95vel95NC}\tag{39}\] In addition, the nuclear dynamics of Eq. 33 is governed by \(\mathbf{F}_s^{\mathrm{HF}}\) and \(\mathbf{F}_s^{\mathrm{HF-\dot{R}}}\), obtained from the derivative of the norm-conserving Lagrangian, since \(\mathbf{F}_s^{S}\) is zero.
With an analogous procedure to the one followed in Refs. [59]–[61] for the inclusion of magnetic fields in the PAW method, we define the nuclear velocity-including PAW transformation operator by substituting the atomic orbitals for fixed nuclei with the velocity-including ones, defined in Eq. 9 . We assume that the atomic orbitals centred on a nucleus feel only the motion of that nucleus; in other words, atomic orbitals are those obtained as if each nucleus were isolated, and then, we exploit the superposition principle to construct the transformation operator for the entire wavefunction. The pseudo-wavefunction \(\ket{\tilde{\psi}}\) is mapped to the velocity-including all-electron one by the transformation operator for velocity-including orbitals \(\hat{\mathcal{T}}_{\dot{\mathrm{R}}}\), \[\begin{align} &\ket{\psi}=\hat{\mathcal{T}}_{\dot{\mathrm{R}}}\ket{\tilde{\psi}} \tag{40}\\ &\hat{\mathcal{T}}_{\dot{\mathrm{R}}}=\hat{\mathbb{1}}+\sum_{s,i}\hat{t}^s_{\dot{\mathbf{R}}_s}, \tag{41}\\ &\hat{t}^s_{\dot{\mathbf{R}}_s}=e^{i\alpha_s(\hat{\mathbf{r}})} \left(\ket{\phi^{\mathbf{R}_s(t)}_{si}}- \ket{\tilde{\phi}^{\mathbf{R}_s(t)}_{si}}\right)\bra{\tilde{p}^{\mathbf{R}_s(t)}_{si}}e^{-i\alpha_s(\hat{\mathbf{r}})}. \nonumber \end{align}\] This expression coincides with the one given in the Supplementary of Ref. [14]. In the following, for brevity we indicate the nuclear velocity phase using \(\alpha_{s}(\hat{\mathbf{r}})=\frac{m}{\hbar}\dot{\mathbf{R}}_s(t)\cdot(\hat{\mathbf{r}}-\mathbf{R}_s(t))\), as in Eq. 4 .
The transformation operator depends on both the time-dependent nuclear positions and velocities. For a complete basis set, the relations for the standard PAW of Eq. 21 , assuming that \(\mathbf{r} \in \Omega_s\), are transformed as \[\begin{align} \sum_{i}\bra{\mathbf{r}}\left(e^{i\alpha_s(\hat{\mathbf{r}})}\ket{\tilde{\phi}^{\mathbf{R}_{s}}_{si}}\bra{\tilde{p}^{\mathbf{R}_s}_{si}}e^{-i\alpha_s(\hat{\mathbf{r}})}\right)e^{i\alpha_s(\hat{\mathbf{r}})}\ket{\tilde{\psi}}\\=e^{i\alpha_s(\mathbf{r})}\braket{\mathbf{r}|\tilde{\psi}}, \\ \sum_{i}\bra{\tilde{\psi}}\left(e^{i\alpha_s(\hat{\mathbf{r}})}\ket{\tilde{p}^{\mathbf{R}_{s}}_{si}}\bra{\tilde{\phi}^{\mathbf{R}_s}_{si}}e^{-i\alpha_s(\hat{\mathbf{r}})}\right)e^{i\alpha_s(\hat{\mathbf{r}})}\ket{\mathbf{r}}\\=\braket{\tilde{\psi}|\mathbf{r}} e^{-i\alpha_s(\mathbf{r})}. \end{align} \label{eq:identity95augmentation95region95vel}\tag{42}\] The transformation of the operators of Eq. 22 - also given in the Supplementary of Ref. [14] - is generalised to an expression analogous to the case of the presence of a magnetic field [59] \[\begin{align} &\tilde{O}=O+\sum_{s,ij} e^{i\alpha_s(\hat{\mathbf{r}})}\ket{\tilde{p}^{\mathbf{R}_s(t)}_{si}} \Delta O^{s, \mathrm{\dot{R}}}_{ij}\bra{\tilde{p}^{\mathbf{R}_s(t)}_{sj}}e^{-i\alpha_s(\hat{\mathbf{r}})},\\ &\Delta O^{s, \mathrm{\dot{R}}}_{ij}=\braket{\phi^{\mathbf{R}_s(t)}_{si}|e^{-i\alpha_s(\hat{\mathbf{r}})}Oe^{i\alpha_s(\hat{\mathbf{r}})}|\phi^{\mathbf{R}_s(t)}_{sj}}\\& -\braket{\tilde{\phi}^{\mathbf{R}_s(t)}_{si}|e^{-i\alpha_s(\hat{\mathbf{r}})}Oe^{i\alpha_s(\hat{\mathbf{r}})}|\tilde{\phi}^{\mathbf{R}_s(t)}_{sj}}. \label{eq:otildevipaw} \end{align}\tag{43}\] The identity operator for the velocity-including case is transformed as the generalisation of the standard expression given in Eq. 24 \[\begin{align} S_{\rm R, \dot{R}}=\hat{\mathcal{T}}_{\dot{\mathrm{R}}}^{\dagger}\hat{\mathcal{T}}_{\dot{\mathrm{R}}}=&\mathbb{1}+\sum_{s,ij}e^{i\alpha_s(\hat{\mathbf{r}})}\ket{\tilde{p}^{\mathbf{R}_s(t)}_{si}}Q^s_{ij}\bra{\tilde{p}^{\mathbf{R}_s(t)}_{sj}}e^{-i\alpha_s(\hat{\mathbf{r}})},\\ Q^s_{ij}=&\braket{\phi^{\mathbf{R}_s(t)}_{si}|\phi^{\mathbf{R}_s(t)}_{sj}}-\braket{\tilde{\phi}^{\mathbf{R}_s(t)}_{si}|\tilde{\phi}^{\mathbf{R}_s(t)}_{sj}}. \end{align} \label{eq:TdaggerT95identity95vel}\tag{44}\] \(S_{\rm R, \dot{R}}\) depends on both the nuclear positions and velocities.
In this section, we derive the non-adiabatic Ehrenfest dynamics by using the velocity-including PAW method. The Ehrenfest Lagrangian is obtained using Eq. 14 with the transformation \(\hat{\mathcal{T}}_{\dot{\mathrm{R}}}\), defined in Eq. 41 , as \[\begin{align} &\mathcal{L}_{\rm R, \dot{R}}^{\mathrm{PAW}}=\sum_s\frac{M_s|\dot{\mathbf{R}}_s(t)|^2}{2}-\sum_{I=1}^{N_{\rm el}}\Bigg[\braket{\tilde{\psi}_I(t)|\hat{H}^{\rm PAW}_{\mathrm{R,\dot{R}}}|\tilde{\psi}_I(t)}\\&+\frac{i\hbar}{2}\left(\frac{d\bra{\tilde{\psi}_I(t)}}{dt}S_{\rm R, \dot{R}}\ket{\tilde{\psi}_I(t)}-\bra{\tilde{\psi}_I(t)}S_{\rm R, \dot{R}}| \frac{d\ket{\tilde{\psi}_I(t)}}{dt}\right)\Bigg], \end{align} \label{eq:PAW95lagrangian95VI}\tag{45}\] where \[\begin{align} \hat{H}^{\rm PAW}_{\mathrm{R,\dot{R}}}= &\hat{\mathcal{T}}_{\dot{\mathrm{R}}}^{\dagger}\hat{H}^{\mathrm{AE}}(\hat{\mathbf{r}},\hat{\mathbf{p}};\{\mathbf{R}_s(t)\}) \hat{\mathcal{T}}_{\dot{\mathrm{R}}}+\frac{i\hbar}{2}\left( \frac{d\hat{\mathcal{T}}_{\dot{\mathrm{R}}}^{\dagger}}{dt}\hat{\mathcal{T}}_{\dot{\mathrm{R}}}-\hat{\mathcal{T}}_{\dot{\mathrm{R}}}^{\dagger}\frac{d\hat{\mathcal{T}}_{\dot{\mathrm{R}}}}{dt}\right). \label{eq:HPAWRdotR} \end{align}\tag{46}\] We are going to drop \(\{\mathbf{R}_s(t)\}\) in the following. The contribution from the temporal derivative is crucial to obtain the correct Hamiltonian. By using \(e^{-i\alpha_s(\hat{\mathbf{r}})}H^{\mathrm{AE}}(\hat{\mathbf{r}}, \hat{\mathbf{p}})e^{i\alpha_s(\hat{\mathbf{r}})}=H^{\mathrm{AE}}(\hat{\mathbf{r}},\hat{\mathbf{p}}+m\dot{\mathbf{R}}_s(t))\) and the transformation of the operators expressed as in Eq. 43 , it follows that (see App. 8) \[\begin{align} &\hat{\mathcal{T}}_{\dot{\mathrm{R}}}^{\dagger}\hat{H}^{\mathrm{AE}}(\hat{\mathbf{r}},\hat{\mathbf{p}}) \hat{\mathcal{T}}_{\dot{\mathrm{R}}}=\\ &\frac{\hat{\mathbf{p}}^2}{2m}+V^{\mathrm{loc}}(\hat{\mathbf{r}})+\sum_s e^{i\alpha_s(\hat{\mathbf{r}})} v_s^{\mathrm{nl}}e^{-i\alpha_s(\hat{\mathbf{r}})} \\ &+\sum_{s,ij} e^{i\alpha_s(\hat{\mathbf{r}})}\ket{\tilde{p}^{\mathbf{R}_s(t)}_{si}} \Delta \mathcal{H}^s_{ij}\bra{\tilde{p}^{\mathbf{R}_s(t)}_{sj}}e^{-i\alpha_s(\hat{\mathbf{r}})}, \end{align} \label{eq:PAW95Hamiltonian}\tag{47}\] where \[\begin{align} & \mathcal{H}= H^{\mathrm{AE}}(\hat{\mathbf{r}},\hat{\mathbf{p}}+m\dot{\mathbf{R}}_s(t))-H^{\mathrm{AE}}(\hat{\mathbf{r}},\hat{\mathbf{p}})\\ &=\dot{\mathbf{R}}_s(t)\cdot \hat{\mathbf{p}}+\frac{m}{2}|\dot{\mathbf{R}}_s(t)|^2, \nonumber \end{align}\]
and we stress that \(\Delta \mathcal{H}^s_{ij}\) is computed as in Eq. 23 . To determine the contribution to the Hamiltonian in Eq. 46 originating in the temporal derivatives of the transformation, as expressed in detail in Appendix 9, we compute \(\displaystyle \frac{d\hat{\mathcal{T}}_{\dot{\mathrm{R}}}}{dt}\) that, using Eq. 11 and defining \[\mathbf{D}_s(\hat{\mathbf{r}},\hat{\mathbf{p}})=\frac{i}{\hbar}\left(m\mathbf{\ddot{R}}_s(t)\cdot \left(\hat{\mathbf{r}}-\mathbf{R}_s(t)\right)-\dot{\mathbf{R}}_s(t)\cdot \mathbf{\hat{p}}\right), \label{eq:orbital95derivative}\tag{48}\] is expressed as \[\frac{d\hat{\mathcal{T}}_{\dot{\mathrm{R}}}}{dt}= \sum_s \left[\mathbf{D}_s(\hat{\mathbf{r}},\hat{\mathbf{p}}),\hat{t}^s_{\dot{\mathbf{R}}_s}\right]. \label{eq:derivative95T95velocity95including95dotR}\tag{49}\] \(\displaystyle \frac{d\hat{\mathcal{T}}_{\dot{\mathrm{R}}}}{dt}\) vanishes outside each augmentation region \(\Omega_s\). As explained in detail in Appendix 9, we obtain
\[\begin{align} &\frac{i\hbar}{2}\left( \frac{d\hat{\mathcal{T}}_{\dot{\mathrm{R}}}^{\dagger}}{dt}\hat{\mathcal{T}}_{\dot{\mathrm{R}}}-\hat{\mathcal{T}}_{\dot{\mathrm{R}}}^{\dagger}\frac{d\hat{\mathcal{T}}_{\dot{\mathrm{R}}}}{dt}\right)= -i\hbar\sum_{s,ij} \Bigg(e^{i\alpha_s(\hat{\mathbf{r}})} \ket{\tilde{p}^{\mathbf{R}_s(t)}_{sj}} \Delta \mathbf{D}^{s,\mathrm{\dot{R}}}_{ij}\bra{\tilde{p}^{\mathbf{R}_s(t)}_{si}} e^{-i\alpha_s(\hat{\mathbf{r}})}-\left\{\frac{\mathbf{D}_s}{2},e^{i\alpha_s(\hat{\mathbf{r}})} \ket{p^{\mathbf{R}_s(t)}_{si}}Q^s_{ij}\bra{p^{\mathbf{R}_s(t)}_{sj}}e^{-i\alpha_s(\hat{\mathbf{r}})}\right\}\Bigg). \end{align} \label{eq:acceleration95PAW95lagrangian}\tag{50}\] Explicitly: \[\begin{align} \Delta \mathbf{D}^{s,\mathrm{\dot{R}}}_{ij}=&\frac{i}{\hbar}\Bigg(\braket{\phi^{\mathbf{R}_s(t)}_{si}|\left[ m\mathbf{\ddot{R}}_s(t)\cdot \left[\hat{\mathbf{r}}-\mathbf{R}_s(t)\right]-\dot{\mathbf{R}}_s(t)\cdot \mathbf{\hat{p}}-m|\dot{\mathbf{R}}_s(t)|^2\right]|\phi^{\mathbf{R}_s(t)}_{sj}} \\&-\braket{\tilde{\phi}^{\mathbf{R}_s(t)}_{si}|\left[m\mathbf{\ddot{R}}_s(t)\cdot \left[\hat{\mathbf{r}}-\mathbf{R}_s(t)\right]-\dot{\mathbf{R}}_s(t)\cdot \mathbf{\hat{p}}-m|\dot{\mathbf{R}}_s(t)|^2\right]|\tilde{\phi}^{\mathbf{R}_s(t)}_{sj}}\Bigg). \end{align}\] Summing Eqs. 47 and 50 , we obtain the effective Hamiltonian from Eq. 46 \[\begin{align} &\hat{H}^{\mathrm{PAW}}_{\mathrm{R,\dot{R}}}= \frac{\hat{\mathbf{p}}^2}{2m}+V^{\mathrm{loc}}(\hat{\mathbf{r}})+\sum_{s}e^{\frac{i}{\hbar}m\dot{\mathbf{R}}_s(t)\cdot \hat{\mathbf{r}}}v_{s}^{\mathrm{nl}}e^{-\frac{i}{\hbar}m\dot{\mathbf{R}}_s(t)\cdot \hat{\mathbf{r}}} +\sum_{si}\Bigg( m\ddot{\mathbf{R}}_s(t) \cdot e^{i\alpha_s(\hat{\mathbf{r}})} \ket{p^{\mathbf{R}_s(t)}_{si}}\mathbf{d}^s_{ij}\bra{p^{\mathbf{R}_s(t)}_{sj}}e^{-i\alpha_s(\hat{\mathbf{r}})} \\& -\frac{1}{2}\left\{m \ddot{\mathbf{R}}_s(t)\cdot\left[\hat{\mathbf{r}}-\mathbf{R}_s(t)\right]-\dot{\mathbf{R}}_s(t)\cdot \hat{\mathbf{p}}+\frac{m|\dot{\mathbf{R}}_s(t)|^2}{2},e^{i\alpha_s(\hat{\mathbf{r}})} \ket{p^{\mathbf{R}_s(t)}_{si}}Q^s_{ij}\bra{p^{\mathbf{R}_s(t)}_{sj}}e^{-i\alpha_s(\hat{\mathbf{r}})} \right\}\Bigg). \end{align} \label{eq:PAW95Hamiltonian-Rdot}\tag{51}\] We have defined, following the notation used in Ref. [62], the quantity \[\begin{align} \mathbf{d}^{s}_{ij}=&m\mathbf{\ddot{R}}_s(t)\cdot\Bigg(\braket{\phi^{\mathbf{R}_s(t)}_{si}| \left[ \hat{\mathbf{r}}-\mathbf{R}_s(t) \right]|\phi^{\mathbf{R}_s(t)}_{sj}}-\braket{\tilde{\phi}^{\mathbf{R}_s(t)}_{si}|\left[\hat{\mathbf{r}}-\mathbf{R}_s(t)\right] |\tilde{\phi}^{\mathbf{R}_s(t)}_{sj}}\Bigg). \end{align}\]
The last line of Eq. 51 are additional terms compared to the nuclear velocity-including pseudo-potential Hamiltonian given in Ref. [14]. While vanishing in the norm-conserving case, they can not be neglected for ultra-soft pseudo-potentials, however they do not raise any particular computational issues. Indeed, \(\displaystyle \mathbf{d}^s_{ij}\) is easily computed in the ultra-soft pseudo-potential implementations [57], [62]. Also the term \(\displaystyle \left\{\left[\hat{\mathbf{r}}-\mathbf{R}_s(t)\right],e^{i\alpha_s(\hat{\mathbf{r}})}\ket{p^{\mathbf{R}_s(t)}_{si}}Q^s_{ij}\bra{p^{\mathbf{R}_s(t)}_{sj}}e^{-i\alpha_s(\hat{\mathbf{r}})}\right\}\) is well-defined since it is localised in each augmentation region.
The \(\Delta \hat{\mathbf{p}}^s_{ij}\) terms, appearing in Eq. 31 , are cancelled by the momentum translation of the Hamiltonian due to the nuclear velocity-dependent phases. This removes the spurious couplings between electronic states present in Eq. 32 , emerging even when the entire system undergoes a global translation at constant speed. The difference between the two formulations can be appreciated in the comparison of Table [tab:hamiltonians95dynamics].
The electronic dynamics are described by the standard Euler-Lagrange equation (Eq. 15 ) \[i\hbar S\frac{d\ket{\tilde{\psi}_I(t)}}{dt}=\left(\hat{H}^{\mathrm{PAW}}_{\mathrm{R,\dot{R}}} -\frac{1}{2}\frac{dS_{\mathrm{R,\dot{R}}}}{dt}\right)\ket{\tilde{\psi}_I(t)}, \label{eq:Schrodinger95PAW}\tag{52}\] where \[\begin{align} & \hat{H}^{\mathrm{PAW}}_{\mathrm{R,\dot{R}}} -\frac{1}{2}\frac{dS_{\mathrm{R,\dot{R}}}}{dt}=\\&\hat{\mathcal{T}}_{\dot{\mathrm{R}}}^{\dagger}\hat{H}^{\mathrm{AE}}(\hat{\mathbf{r}},\hat{\mathbf{p}}) \hat{\mathcal{T}}_{\dot{\mathrm{R}}}-i\hbar \hat{\mathcal{T}}_{\dot{\mathrm{R}}}^{\dagger}\frac{d\hat{\mathcal{T}}_{\dot{\mathrm{R}}}}{dt}. \end{align}\] In the presence of nuclear acceleration contributions in the Lagrangian, the nuclear dynamics is governed by Eq. 18 with the nuclear acceleration terms \[M_s\ddot{\mathbf{R}}_s(t)=\mathbf{F}_s^{\mathrm{HF}}+\mathbf{F}_s^{\mathrm{HF-\dot{R}}}+\mathbf{F}_s^{\mathrm{HF-\ddot{R}}}+\mathbf{F}_s^{\mathrm{S}}+\mathbf{F}_s^{\mathrm{S-\dot{R}}}. \label{eq:Ehrenfest95ions951}\tag{53}\] A nuclear acceleration Hellmann-Feynman force \(\mathbf{F}_s^{\mathrm{HF-\ddot{R}}}\) appear in the equations for the ions compared to Eq. 33 . The other Hellmann-Feynman terms have analogous expressions to the ones given in Eqs. 34 , \[\begin{align} &\mathbf{F}_s^{\mathrm{HF}}=-\sum_{I=1}^{N_{\rm el}}\braket{\tilde{\psi}_I(t)|\frac{\partial \hat{H}^{\mathrm{PAW}}_{\mathrm{R,\dot{R}}}}{\partial \mathbf{R}_s}|\tilde{\psi}_I(t)},\\ &\mathbf{F}_s^{\mathrm{HF-\dot{R}}}=\sum_{I=1}^{N_{\rm el}}\frac{d}{dt}\braket{\tilde{\psi}_I(t)|\frac{\partial \hat{H}^{\mathrm{PAW}}_{\mathrm{R,\dot{R}}}}{\partial \dot{\mathbf{R}}_s}|\tilde{\psi}_I(t)},\\ \end{align} \label{eq:HF95force95pseudo95vel}\tag{54}\] \[\begin{align} &\mathbf{F}_s^{\mathrm{HF-\ddot{R}}}=-\sum_{I=1}^{N_{\rm el}}\frac{d^2}{dt^2}\braket{\tilde{\psi}_I(t)|\frac{\partial \hat{H}^{\mathrm{PAW}}_{\mathrm{R,\dot{R}}}}{\partial \ddot{\mathbf{R}}_s}|\tilde{\psi}_I(t)}. \end{align} \label{eq:HF95force95pseudo95vel95bis}\tag{55}\]
Crucially, the nuclear velocity derivative entering in \(\mathbf{F}_s^{\mathrm{HF-\dot{R}}}\) differs from that of Eq. 34 , being \[\begin{align} & \frac{\partial \hat{H}^{\mathrm{PAW}}_{\mathrm{R,\dot{R}}}}{\partial \dot{\mathbf{R}}_{s}}=\frac{im}{\hbar}\left[\mathbf{r},v_s^{\mathrm{nl}}\right]+\frac{1}{2}\sum_{ij}\Bigg\{ \hat{\mathbf{p}}-\frac{m\dot{\mathbf{R}}_s(t)}{2},\nonumber\\ &e^{i\alpha_s(\hat{\mathbf{r}})}\ket{p^{\mathbf{R}_s(t)}_{si}}Q^s_{ij}\bra{p^{\mathbf{R}_s(t)}_{sj}}e^{-i\alpha_s(\hat{\mathbf{r}})}\Bigg\}\label{eq:PAW95derivative95dotR}. \end{align}\tag{56}\] The nuclear acceleration derivative is expressed as \[\begin{align} &\frac{\partial \hat{H}^{\mathrm{PAW}}_{\mathrm{R,\dot{R}}}}{\partial \ddot{\mathbf{R}}_{s}}=e^{i\alpha_s(\hat{\mathbf{r}})} \ket{p^{\mathbf{R}_s(t)}_{si}}\mathbf{d}^s_{ij}\bra{p^{\mathbf{R}_s(t)}_{sj}}e^{-i\alpha_s(\hat{\mathbf{r}})}\\ -&\frac{1}{2}\left\{\left[\hat{\mathbf{r}}-\mathbf{R}_s(t)\right],e^{i\alpha_s(\hat{\mathbf{r}})}\ket{p^{\mathbf{R}_s(t)}_{si}}Q^s_{ij}\bra{p^{\mathbf{R}_s(t)}_{sj}}e^{-i\alpha_s(\hat{\mathbf{r}})}\right\}\Bigg). \end{align} \label{eq:PAW95derivative95ddotR}\tag{57}\] Additional forces originating from \(S_{\mathrm{R,\dot{R}}}\) are \(\mathbf{F}_s^{S}\) (analogous to Eq. 35 and reducing to Eq. 27 in the adiabatic limit) and its nuclear velocity-generalisation \(\mathbf{F}_s^{\rm S-\mathrm{\dot{R}}}\), appearing only in the velocity-including PAW, \[\begin{align} \mathbf{F}_s^{\rm S}=&\frac{i\hbar}{2}\sum_{I=1}^{N_{\rm el}}\Bigg(\bra{\tilde{\psi}_I(t)}\frac{\partial S_{\mathrm{R,\dot{R}}}}{\partial \mathbf{R}_s}\frac{d \ket{\tilde{\psi}_I(t)}}{dt}\\ & -\frac{d\bra{\tilde{\psi}_I(t)}}{dt}\frac{\partial S_{\mathrm{R,\dot{R}}}}{\partial \mathbf{R}_s}\ket{\tilde{\psi}_I(t)}\Bigg),\\ \end{align}\] \[\begin{align} \mathbf{F}_s^{\mathrm{S-\dot{R}}}=&-\frac{i\hbar}{2}\frac{d}{dt}\sum_{I=1}^{N_{\rm el}}\Bigg(\bra{\tilde{\psi}_I(t)}\frac{\partial S_{\mathrm{R,\dot{R}}}}{\partial \dot{\mathbf{R}}_s}\frac{d \ket{\tilde{\psi}_I(t)}}{dt}\\ & -\frac{d\bra{\tilde{\psi}_I(t)}}{dt}\frac{\partial S_{\mathrm{R,\dot{R}}}}{\partial \dot{\mathbf{R}}_s}\ket{\tilde{\psi}_I(t)}\Bigg), \end{align}\] where, from Eq. 44 , the nuclear velocity derivative of \(S_{\mathrm{R,\dot{R}}}\) is \[\frac{\partial S_{\mathrm{R,\dot{R}}}}{\partial \dot{\mathbf{R}}_s}=\Big[ \hat{\mathbf{r}},e^{i\alpha_s(\hat{\mathbf{r}})} \ket{p^{\mathbf{R}_s(t)}_{si}}Q^s_{ij}\bra{p^{\mathbf{R}_s(t)}_{sj}}e^{-i\alpha_s(\hat{\mathbf{r}})} \Big]. \label{eq:S95derivative95dotR}\tag{58}\] Finally, the conserved energy is computed according to Eq. 19 that holds in the presence of nuclear acceleration terms in the Lagrangian. The linearity of the Lagrangian in the nuclear acceleration forbids the presence of derivatives of the position of order higher than the second in the conserved energy. The complete expression is \[\begin{align} &E^{\mathrm{PAW}}_{\rm R,\dot{R}}=\sum_s\frac{M_s|\dot{\mathbf{R}}_s(t)|^2}{2}+\sum_{I=1}^{N_{\rm el}}\Bigg\{\braket{\tilde{\psi}_I(t)|\hat{H}^{\rm PAW}_{\mathrm{R,\dot{R}}}|\tilde{\psi}_I(t)}\nonumber\\ &-\sum_s\Bigg[\braket{\tilde{\psi}_I(t)|\frac{\partial \hat{H}^{\rm PAW}_{\mathrm{R,\dot{R}}}}{\partial \dot{\mathbf{R}}_s}\cdot \dot{\mathbf{R}}_s(t)+\frac{\partial \hat{H}^{\rm PAW}_{\mathrm{R,\dot{R}}}}{\partial \ddot{\mathbf{R}}_s}\cdot \ddot{\mathbf{R}}_s(t)|\tilde{\psi}_I(t)} \nonumber\\ &-\frac{d}{dt}\left( \braket{\tilde{\psi}_I(t)|\frac{\partial \hat{H}^{\rm PAW}_{\mathrm{R,\dot{R}}}}{\partial \ddot{\mathbf{R}}_s}|\tilde{\psi}_I(t)}\right)\cdot \dot{\mathbf{R}}_s(t)+\label{eq:energy95VI}\\ &\frac{i\hbar}{2}\Bigg(\frac{d\bra{\tilde{\psi}_I(t)}}{dt}\frac{\partial S_{\rm R, \dot{R}}}{\partial \dot{\mathbf{R}}_s}\ket{\tilde{\psi}_I(t)}-\nonumber\\ &\bra{\tilde{\psi}_I(t)}\frac{\partial S_{\rm R, \dot{R}}}{\partial \dot{\mathbf{R}}_s}| \frac{d\ket{\tilde{\psi}_I(t)}}{dt}\Bigg)\cdot \dot{\mathbf{R}}_s(t)\Bigg]\Bigg\}. \nonumber \end{align}\tag{59}\] The nuclear velocity and acceleration derivatives are reported explicitly in Eqs. 56 , 57 and 58 . For a direct comparison between the main results, including the conserved energy, with nuclear velocity-including atomic orbitals and with the rigidly translated atomic orbitals see Table [tab:hamiltonians95dynamics].
In the case of a norm-conserving pseudo-potential, \(S=1\) and the augmentation charges are zero \(Q^s_{ij}=0\). Moreover, for a smooth operator within the augmentation region, such as the position, the matrix elements with the all-electron and pseudo-partial waves almost coincide \[\braket{\phi^{\mathbf{R}_s}_{si}|\hat{\mathbf{r}}|\phi^{\mathbf{R}_s}_{sj}}\approx \braket{\tilde{\phi}^{\mathbf{R}_s}_{si}|\hat{\mathbf{r}}|\tilde{\phi}^{\mathbf{R}_s}_{sj}}. \label{eq:position95pseudo95all95electron}\tag{60}\] As a consequence, the Hamiltonian reduces to the one given by Ref. [14], \[\begin{align} \hat{H}^{\mathrm{PAW-NC}}_{\mathrm{R},\mathrm{\dot{R}}}= &\frac{\hat{\mathbf{p}}^2}{2m}+V^{\mathrm{loc}}\\&+\sum_{s}e^{\frac{i}{\hbar}m\dot{\mathbf{R}}_s(t)\cdot \hat{\mathbf{r}}}v_{s}^{\mathrm{nl}}e^{-\frac{i}{\hbar}m\dot{\mathbf{R}}_s(t)\cdot \hat{\mathbf{r}}}. \end{align} \label{eq:PAW95Hamiltonian-NC}\tag{61}\] The evolution of the electronic system is determined by \[i\hbar \frac{d\ket{\tilde{\psi}_I(t)}}{dt}=\hat{H}^{\mathrm{PAW-NC}}_{\mathrm{R},\mathrm{\dot{R}}} \ket{\tilde{\psi}_I(t)}, \label{eq:Schrodinger95PAW95NC}\tag{62}\] which guarantees Galilean invariance since the motion at constant velocity of the entire system does not cause any spurious electronic transition. Indeed, there are no terms in the Hamiltonian mixing different electronic orbitals, overcoming the paradox of the non-adiabatic Ehrenfest dynamics without the nuclear velocity-dependent phases.
In the nuclear dynamics, the nuclear acceleration contribution to the force is approximately zero \(\mathbf{F}_s^{\mathrm{HF-\ddot{R}}}\approx 0\), whereas the other two reduce to \[\begin{align} &\mathbf{F}_s^{\mathrm{HF}}=-\sum_{I=1}^{N_{\rm el}}\braket{\tilde{\psi}_I(t)|\frac{\partial \hat{H}^{\mathrm{PAW-NC}}_{\mathrm{R,\dot{R}}}}{\partial \mathbf{R}_s}|\tilde{\psi}_I(t)},\tag{63}\\ &\mathbf{F}_s^{\mathrm{HF-\dot{R}}}=\sum_{I=1}^{N_{\rm el}}\frac{d}{dt}\braket{\tilde{\psi}_I(t)|\frac{\partial \hat{H}^{\rm PAW-NC}_{\mathrm{R},\dot{\mathrm{R}}}}{\partial \dot{\mathbf{R}}_s}|\tilde{\psi}_I(t)} \tag{64}. \end{align}\]
These expressions give the force-constant matrix and the Born effective charges presented in Ref. [14], which neglects the nuclear-acceleration contribution. The conserved energy is \[\begin{align} & E^{\mathrm{PAW-NC}}_{\mathrm{R},\dot{\mathrm{R}}}=\sum_s \frac{M_s|\dot{\mathbf{R}}_s(t)|^2}{2}+ \sum_{I=1}^{N_{\rm el}}\bra{\tilde{\psi}_I(t)}\\&\left(\hat{H}^{\rm PAW-\mathrm{NC}}_{\mathrm{R},\dot{\mathrm{R}}}-\sum_s\frac{\partial \hat{H}^{\rm PAW-NC}_{\mathrm{R},\dot{\mathrm{R}}}}{\partial \dot{\mathbf{R}}_s}\cdot \dot{\mathbf{R}}_s\right)\ket{\tilde{\psi}_I(t)}\Bigg], \end{align} \label{eq:energy95VI95NC}\tag{65}\] where, from Eq. 56 , \[\frac{\partial \hat{H}^{\mathrm{PAW-NC}}_{\mathrm{R},\mathrm{\dot{R}}}}{\partial \dot{\mathbf{R}}_{s}}=\frac{im}{\hbar}\left[\mathbf{r},v_s^{\mathrm{nl}}\right].\] that coincides with the results of Ref. [14].
In this work, we constructed the pseudo-potential Hamiltonian through a PAW transformation that employs velocity including atomic orbitals, in both the norm-conserving and ultrasoft cases, with the aim of obtaining the non-adiabatic Ehrenfest dynamics, as summarised in Table [tab:hamiltonians95dynamics]. While the local part of the Hamiltonian remains unaffected, Peierls-like phases that depend on the nuclear velocity appear in the non-local part of the potential. Additional nuclear velocity and acceleration contributions, depending on the augmentation charge, are also present in the ultrasoft pseudopotential case. Indeed, the nuclear acceleration contributions, previously overlooked in the literature, cannot be neglected in the ultrasoft case, while in the norm-conserving case they are small and can be disregarded. Thus, the use of nuclear-velocity including atomic orbitals in the PAW transformation enables the proper reproduction of the all -electron properties in the pseudo-potential Hamiltonian when nuclei are moving, at any order in the nuclear velocity and addressing specific issues that arise in standard approaches. First of all, since electronic wavefunctions acquire spherical harmonics components at any order when nuclei are moving, neglecting the nuclear-velocity dependent phases in the atomic orbitals implies that the static orbital basis set is incomplete, making it inadequate for the description of the all-electron problem. Furthermore, our equations fully restore Galilean invariance through nuclear-velocity-dependent phases that allow for the removal of spurious couplings between different electronic states, which would otherwise appear also when the entire system undergoes a global translation. This is highlighted in Table [tab:hamiltonians95dynamics] in the comparison between the electronic dynamics obtained with the rigid translation of electronic orbitals to the time-dependent nuclear positions and the use of the nuclear-velocity dependent orbitals. Finally, we generalise to ultrasoft pseudo-potentials the results of Ref. [14], which is limited to the norm-conserving case, showing that additional nuclear velocity and acceleration dependent contributions have to be accounted for in the assessment of dynamical vibrational responses as the force constant matrix and the Born effective charge tensor.
We acknowledge the MORE-TEM ERC-SYN project, Grant Agreement No. 951215. PF acknowledge also the funding from the project Ateneo 2025 by Sapienza - University of Rome (code: B83C25004300005). We thank Massimiliano Stengel, Raffaele Resta, Antimo Marrazzo and Giorgio Sangiovanni for useful discussions and suggestions.
Consider a system where the electronic part of the Lagrangian is \[\begin{align} &\mathcal{L}_{\rm el}=-\sum_{I=1}^{N_{\rm el}}\Bigg[\braket{\tilde{\psi}_I(t)|\hat{O}(t)|\tilde{\psi}_I(t)}+\frac{i\hbar}{2}\left(\frac{d\bra{\tilde{\psi}_I(t)}}{dt}S\ket{\tilde{\psi}_I(t)}-\bra{\tilde{\psi}_I(t)}S| \frac{d\ket{\tilde{\psi}_I(t)}}{dt}\right)\Bigg] \end{align}\] and the the wavefunction is normalized such that \[\braket{\psi(t)|S|\psi(t)}=1,\] then the temporal evolution operator \(\hat{O}(t)\) is non-hermitian operator \[i\hbar S {\ket{\dot{\psi}(t)}}=O(t)\ket{\psi(t)}.\]
In order to prove this, we compute the temporal derivative of the norm \[\begin{align} &\frac{d}{dt}\left(\braket{\psi(t)|S|\psi(t)}\right)=0,\\ &\braket{\dot{\psi}(t)|S|\psi(t)}+\braket{\psi(t)|S|\dot{\psi}(t)}+\braket{\psi(t)|\frac{dS(t)}{dt}|\psi(t)}=0.\\ \end{align}\] Substituting the Schrödinger equations for \(\ket{\psi}\) and \(\bra{\psi}\), we obtain \[\begin{align} \braket{\psi|-O(t)+O^{\dagger}(t)+\frac{dS(t)}{dt}|\psi}=0,\\ \end{align}\] following that \[\frac{dS(t)}{dt}=O(t)-O^{\dagger}(t).\] In conclusion, a time-dependent \(S\) implies that the evolution is not determined by a Hermitian operator. We remark that the evolution is unitary because of the definition of the norm. A non-Hermitian evolution operator is needed to conserve the norm.
According to Eq. 43 of the main text, here reported, the operators transforms with the velocity-including PAW transformation as \[\begin{align} &\tilde{O}=O+\sum_{s,ij} e^{i\alpha_s(\hat{\mathbf{r}})}\ket{\tilde{p}^{\mathbf{R}_s(t)}_{si}} \Delta O^{s, \mathrm{\dot{R}}}_{ij}\bra{\tilde{p}^{\mathbf{R}_s(t)}_{sj}}e^{-i\alpha_s(\hat{\mathbf{r}})},\\ &\Delta O^{s, \mathrm{\dot{R}}}_{ij}=\braket{\phi^{\mathbf{R}_s(t)}_{si}|e^{-i\alpha_s(\hat{\mathbf{r}})}Oe^{i\alpha_s(\hat{\mathbf{r}})}|\phi^{\mathbf{R}_s(t)}_{sj}} -\braket{\tilde{\phi}^{\mathbf{R}_s(t)}_{si}|e^{-i\alpha_s(\hat{\mathbf{r}})}Oe^{i\alpha_s(\hat{\mathbf{r}})}|\tilde{\phi}^{\mathbf{R}_s(t)}_{sj}}. \label{eq:otildevipaw95app} \end{align}\tag{66}\] We remark that \(\Delta O^{s, \mathrm{\dot{R}}}_{ij}\) and \(\Delta O^{s}_{ij}\) of Eq. 23 of the main text indicate the difference between all-electron and pseudo wavefunctions matrix elements with and without the nuclear velocity dependent phases in the orbitals, respectively. Adding and subtracting \(\displaystyle \sum_{s,ij} e^{i\alpha_s(\hat{\mathbf{r}})}\ket{\tilde{p}^{\mathbf{R}_s(t)}_{si}} \Delta O^{s}_{ij}\bra{\tilde{p}^{\mathbf{R}_s(t)}_{sj}}e^{-i\alpha_s(\hat{\mathbf{r}})}\), we can rewrite Eq. 66 as \[\begin{align} &\tilde{O}=O+\sum_{s,ij} \Bigg(e^{i\alpha_s(\hat{\mathbf{r}})}\ket{\tilde{p}^{\mathbf{R}_s(t)}_{si}} \Delta O^{s}_{ij}\bra{\tilde{p}^{\mathbf{R}_s(t)}_{sj}}e^{-i\alpha_s(\hat{\mathbf{r}})}+e^{i\alpha_s(\hat{\mathbf{r}})}\ket{\tilde{p}^{\mathbf{R}_s(t)}_{si}}\left( \Delta O^{s,\mathrm{\dot{R}}}_{ij}-\Delta O^{s}_{ij}\right)\bra{\tilde{p}^{\mathbf{R}_s(t)}_{sj}}e^{-i\alpha_s(\hat{\mathbf{r}})}\Bigg). \end{align} \label{eq:delta95O95dotR}\tag{67}\] where \(\Delta O^{s}_{ij}\) is the difference between the pseudo and all-electron matrix elements without the nuclear velocity-dependent phases, as defined in Eq. 23 . The last term contains the difference between the velocity-including and the standard cases. It is expressed as \[\begin{align} \Delta O^{s,\mathrm{\dot{R}}}_{ij}- \Delta O^s_{ij}=\braket{\phi^{\mathbf{R}_s(t)}_{si}|e^{-i\alpha_s(\hat{\mathbf{r}})}Oe^{i\alpha_s(\hat{\mathbf{r}})}-O|\phi^{\mathbf{R}_s(t)}_{sj}}-\braket{\tilde{\phi}^{\mathbf{R}_s(t)}_{si}|e^{-i\alpha_s(\hat{\mathbf{r}})}Oe^{i\alpha_s(\hat{\mathbf{r}})}-O|\tilde{\phi}^{\mathbf{R}_s(t)}_{sj}} \end{align}\] For instance, for an operator depending on the momentum \(O(\hat{\mathbf{p}})\), the nuclear velocity-dependent phase factor cause a shift in momentum, according to the commutation rules between position and momentum operators \(e^{-i\alpha_s(\hat{\mathbf{r}})}O(\hat{\mathbf{p}})e^{is\alpha_s(\hat{\mathbf{r}})}=O(\hat{\mathbf{p}}+m\dot{\mathbf{R}}_s(t))\), implying that \[\begin{align} \Delta O^{s,\mathrm{\dot{R}}}_{ij}- \Delta O^s_{ij}=\braket{\phi^{\mathbf{R}_s(t)}_{si}|O(\hat{\mathbf{p}}+m\dot{\mathbf{R}}_s(t))-O|\phi^{\mathbf{R}_s(t)}_{sj}}-\braket{\tilde{\phi}^{\mathbf{R}_s(t)}_{si}|O(\hat{\mathbf{p}}+m\dot{\mathbf{R}}_s(t))-O|\tilde{\phi}^{\mathbf{R}_s(t)}_{sj}}. \end{align}\]
The PAW transformation of Eq. 20 for moving nuclei, but without the nuclear velocity-dependent phase factors, can be expressed as \[\begin{align} \hat{\mathcal{T}}=\hat{\mathbb{1}}+\sum_{s}\hat{t}^s,\qquad \hat{t}^s=\sum_i\left(\ket{\phi^{\mathbf{R}_s(t)}_{si}}- \ket{\tilde{\phi}^{\mathbf{R}_s(t)}_{si}}\right)\bra{\tilde{p}^{\mathbf{R}_s(t)}_{si}}, \end{align}\] whereas with the nuclear velocity-dependent phases, it is expressed as in Eq. 41 \[\begin{align} &\hat{\mathcal{T}}_{\dot{\mathrm{R}}}=\hat{\mathbb{1}}+\sum_{s}\hat{t}^s_{\dot{\mathbf{R}}_s},\qquad \hat{t}^s_{\dot{\mathbf{R}}_s}=\sum_ie^{i\alpha_s(\hat{\mathbf{r}})} \left(\ket{\phi^{\mathbf{R}_s(t)}_{si}}- \ket{\tilde{\phi}^{\mathbf{R}_s(t)}_{si}}\right)\bra{\tilde{p}^{\mathbf{R}_s(t)}_{si}}e^{-i\alpha_s(\hat{\mathbf{r}})}. \label{eq:Trdot95app} \end{align}\tag{68}\] By using the equations above, we can also express the difference between the pseudo and the all-electron (local o semi-local) operators as \[\tilde{O}-O=\sum_s \left[(\hat{t}_{\dot{\mathbf{R}}_s}^s)^{\dagger} O+O\hat{t}_{\dot{\mathbf{R}}_s}^s+(\hat{t}_{\dot{\mathbf{R}}_s}^s)^{\dagger}O\hat{t}_{\dot{\mathbf{R}}_s}^s\right]. \label{eq:transformation95detail}\tag{69}\] By applying the identity for the ket and the bra of the pseudo-wavefunction given in Eq. 21 to the first and the second term of the right-hand side, respectively, the expression for the transformation of the PAW operators, given in Eq. 22 , is obtained.
In the calculation of the temporal derivatives of the transformation, we consider the nuclear velocity-including case since the other case is recovered by setting \(\alpha(\hat{\mathbf{r}})=0\). Therefore, using Eq. 11 and defining \[\mathbf{D}_s(\hat{\mathbf{r}},\hat{\mathbf{p}})=\frac{i}{\hbar}\left(m\mathbf{\ddot{R}}_s(t)\cdot (\hat{\mathbf{r}}-\mathbf{R}_s(t))-\dot{\mathbf{R}}_s(t)\cdot \mathbf{\hat{p}}\right), \label{eq:orbital95derivative95app}\tag{70}\] we obtain that \[\frac{d\hat{\mathcal{T}}_{\dot{\mathrm{R}}}}{dt}= \sum_s \left[\mathbf{D}_s(\hat{\mathbf{r}},\hat{\mathbf{p}}),\hat{t}^s_{\dot{\mathbf{R}}_s}\right]. \label{eq:derivative95T95velocity95including95dotR95app}\tag{71}\] In the absence of the nuclear velocity-dependent phase factors, we have the same expression with \(\mathbf{\ddot{R}}_s(t)=0\). From Eq. 71 , it follows that \[\begin{align} &\hat{\mathcal{T}}_{\dot{\mathrm{R}}}^{\dagger}\frac{d\hat{\mathcal{T}}_{\dot{\mathrm{R}}}}{dt}=\hat{\mathcal{T}}_{\dot{\mathrm{R}}}^{\dagger}\sum_s \left[\mathbf{D}_s(\hat{\mathbf{r}},\hat{\mathbf{p}}),\hat{t}^s_{\dot{\mathbf{R}}_s}\right] =\sum_s \left(\mathbf{D}_s\hat{t}^s_{\dot{\mathbf{R}}_s}-\hat{t}^s_{\dot{\mathbf{R}}_s}\mathbf{D}_s+ (\hat{t}^s_{\dot{\mathbf{R}}_s})^{\dagger}\mathbf{D}_s\hat{t}^s_{\dot{\mathbf{R}}_s}-(\hat{t}^s_{\dot{\mathbf{R}}_s})^{\dagger}\hat{t}^s_{\dot{\mathbf{R}}_s}\mathbf{D}_s\right). \end{align}\] By summing and subtracting \(\displaystyle \sum_s (\hat{t}^s_{\dot{\mathbf{R}}_s})^{\dagger}\mathbf{D}_s\)
\[\begin{align} &\hat{\mathcal{T}}_{\dot{\mathrm{R}}}^{\dagger}\frac{d\hat{\mathcal{T}}_{\dot{\mathrm{R}}}}{dt} =\sum_s \bigg[\left((\hat{t}^s_{\dot{\mathbf{R}}_s})^{\dagger}\mathbf{D}_s+\mathbf{D}_s\hat{t}^s_{\dot{\mathbf{R}}_s}+ (\hat{t}^s_{\dot{\mathbf{R}}_s})^{\dagger}\mathbf{D}_s\hat{t}^s_{\dot{\mathbf{R}}_s} \right)-\left((\hat{t}^s_{\dot{\mathbf{R}}_s})^{\dagger}+\hat{t}^s_{\dot{\mathbf{R}}_s}+ (\hat{t}^s_{\dot{\mathbf{R}}_s})^{\dagger}\hat{t}^s_{\dot{\mathbf{R}}_s}\right)\mathbf{D}_s\bigg]. \end{align}\] According to Eq. 69 , the two terms of the right-hand side correspond to the augmentation region resolved difference between the all-electron and pseudo- operators \(\mathbf{D}_s\) and \(\mathbb{1}\). Therefore, by properly applying the identities in the augmentation region given in Eq. 21 , as for the transformation of the operators, we obtain with the hypothesis that augmentations regions do not overlap and performing the same calculation for \(\displaystyle \frac{d\hat{\mathcal{T}}^{\dagger}_{\mathrm{\dot{R}}}}{dt}\hat{\mathcal{T}}_{\dot{\mathrm{R}}}\), \[\begin{align} \hat{\mathcal{T}}_{\dot{\mathrm{R}}}^{\dagger}\frac{d\hat{\mathcal{T}}_{\dot{\mathrm{R}}}}{dt}= &\sum_{sij} e^{i\alpha_s(\hat{\mathbf{r}})} \Big(\ket{p^{\mathbf{R}_s(t)}_{si}}\Delta \mathbf{D}^{s,\mathrm{\dot{R}}}_{ij}\bra{p^{\mathbf{R}_s(t)}_{sj}}e^{-i\alpha_s(\hat{\mathbf{r}})} -\ket{p^{\mathbf{R}_s(t)}_{si}}Q^s_{ij}\bra{p^{\mathbf{R}_s(t)}_{sj}}e^{-i\alpha_s(\hat{\mathbf{r}})} \mathbf{D}_s\Big)\\ \frac{d\hat{\mathcal{T}}^{\dagger}_{\mathrm{\dot{R}}}}{dt}\hat{\mathcal{T}}_{\dot{\mathrm{R}}}=& \sum_{sij} \Big(-e^{i\alpha_s(\hat{\mathbf{r}})} \ket{p^{\mathbf{R}_s(t)}_{si}}\Delta \mathbf{D}^{s,\mathrm{\dot{R}}}_{ij}\bra{p^{\mathbf{R}_s(t)}_{sj}}+\mathbf{D}_se^{i\alpha_s(\hat{\mathbf{r}})} \ket{p^{\mathbf{R}_s(t)}_{si}}Q^s_{ij}\bra{p^{\mathbf{R}_s(t)}_{sj}}\Big)e^{i\alpha_s(\hat{\mathbf{r}})}. \end{align} \label{eq:TdotT95app}\tag{72}\] where, by using that \(e^{-i\alpha_{s}(\hat{\mathbf{r}})}\hat{\mathbf{p}}e^{i\alpha_{s}(\hat{\mathbf{r}})}=\hat{\mathbf{p}}+m\dot{\mathbf{R}}_s(t)\) in combination with the definition of \(\mathbf{D}_s\) given in Eq. 70 , \[\begin{align} \Delta \mathbf{D}^{s,\mathrm{\dot{R}}}_{ij}=&\frac{i}{\hbar}\Bigg(\braket{\phi^{\mathbf{R}_s(t)}_{si}|m\mathbf{\ddot{R}}_s(t)\cdot (\hat{\mathbf{r}}-\mathbf{R}_s(t))-\dot{\mathbf{R}}_s(t)\cdot \mathbf{\hat{p}}-m|\dot{\mathbf{R}}_s(t)|^2|\phi^{\mathbf{R}_s(t)}_{sj}} \\&-\braket{\tilde{\phi}^{\mathbf{R}_s(t)}_{si}|m\mathbf{\ddot{R}}_s(t)\cdot (\hat{\mathbf{r}}-\mathbf{R}_s(t))-\dot{\mathbf{R}}_s(t)\cdot \mathbf{\hat{p}}-m|\dot{\mathbf{R}}_s(t)|^2|\tilde{\phi}^{\mathbf{R}_s(t)}_{sj}}\Bigg). \end{align}\]
Finally, the contribution to the Lagrangian originating from the temporal derivative of the transformation is \[\begin{align} &\frac{i\hbar}{2}\left( \frac{d\hat{\mathcal{T}}_{\dot{\mathrm{R}}}^{\dagger}}{dt}\hat{\mathcal{T}}_{\dot{\mathrm{R}}}-\hat{\mathcal{T}}_{\dot{\mathrm{R}}}^{\dagger}\frac{d\hat{\mathcal{T}}_{\dot{\mathrm{R}}}}{dt}\right)= -i\hbar\sum_{s,ij} \Bigg(e^{i\alpha_s(\hat{\mathbf{r}})} \ket{\tilde{p}^{\mathbf{R}_s(t)}_{sj}} \Delta \mathbf{D}^{s,\mathrm{\dot{R}}}_{ij}\bra{\tilde{p}^{\mathbf{R}_s(t)}_{si}} e^{-i\alpha_s(\hat{\mathbf{r}})}-\left\{\frac{\mathbf{D}_s}{2},e^{i\alpha_s(\hat{\mathbf{r}})} \ket{p^{\mathbf{R}_s(t)}_{si}}Q^s_{ij}\bra{p^{\mathbf{R}_s(t)}_{sj}}e^{-i\alpha_s(\hat{\mathbf{r}})}\right\}\Bigg). \end{align} \label{eq:acceleration95PAW95lagrangian95app}\tag{73}\] In the norm-conserving case, where \(Q^s_{ij}=0\) and \(\hat{\mathcal{T}}_{\dot{\mathrm{R}}}^{\dagger}\hat{\mathcal{T}}_{\dot{\mathrm{R}}}=\mathbb{1}\), the temporal derivative terms reduces to \[\begin{align} \hat{\mathcal{T}}_{\dot{\mathrm{R}}}^{\dagger}\frac{d\hat{\mathcal{T}}_{\dot{\mathrm{R}}}}{dt}=-\frac{d\hat{\mathcal{T}^{\dagger}}_{\{\dot{\mathbf{R}}_s\}}}{dt}\hat{\mathcal{T}}_{\dot{\mathrm{R}}}= &\sum_s e^{i\alpha_s(\hat{\mathbf{r}})} \ket{p^{\mathbf{R}_s(t)}_{si}}\Delta \mathbf{D}^{s,\mathrm{\dot{R}}}_{ij}\bra{p^{\mathbf{R}_s(t)}_{sj}}e^{-i\alpha_s(\hat{\mathbf{r}})} \end{align} \label{eq:TdotT95app95norm95conserving}\tag{74}\] In the absence of nuclear velocity-dependent phase factors, the proof is the same with \(\mathbf{\ddot{R}}_s(t)=0\) and \(\Delta \mathbf{D}^{s,\mathrm{\dot{R}}}_{ij}\) replaced by \(\Delta \mathbf{D}^{s}_{ij}\), which gives a result Eq. 31 of the main text.